Branch data Line data Source code
1 : : /* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2 : : /*************************************************************************
3 : : *
4 : : * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
5 : : *
6 : : * Copyright 2000, 2010 Oracle and/or its affiliates.
7 : : *
8 : : * OpenOffice.org - a multi-platform office productivity suite
9 : : *
10 : : * This file is part of OpenOffice.org.
11 : : *
12 : : * OpenOffice.org is free software: you can redistribute it and/or modify
13 : : * it under the terms of the GNU Lesser General Public License version 3
14 : : * only, as published by the Free Software Foundation.
15 : : *
16 : : * OpenOffice.org is distributed in the hope that it will be useful,
17 : : * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 : : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 : : * GNU Lesser General Public License version 3 for more details
20 : : * (a copy is included in the LICENSE file that accompanied this code).
21 : : *
22 : : * You should have received a copy of the GNU Lesser General Public License
23 : : * version 3 along with OpenOffice.org. If not, see
24 : : * <http://www.openoffice.org/license.html>
25 : : * for a copy of the LGPLv3 License.
26 : : *
27 : : ************************************************************************/
28 : :
29 : :
30 : : #include "rtl/math.h"
31 : :
32 : : #include "osl/diagnose.h"
33 : : #include "rtl/alloc.h"
34 : : #include "rtl/math.hxx"
35 : : #include "rtl/strbuf.h"
36 : : #include "rtl/string.h"
37 : : #include "rtl/ustrbuf.h"
38 : : #include "rtl/ustring.h"
39 : : #include "sal/mathconf.h"
40 : : #include "sal/types.h"
41 : :
42 : : #include <algorithm>
43 : : #include <float.h>
44 : : #include <limits.h>
45 : : #include <math.h>
46 : : #include <stdlib.h>
47 : :
48 : :
49 : : static int const n10Count = 16;
50 : : static double const n10s[2][n10Count] = {
51 : : { 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8,
52 : : 1e9, 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16 },
53 : : { 1e-1, 1e-2, 1e-3, 1e-4, 1e-5, 1e-6, 1e-7, 1e-8,
54 : : 1e-9, 1e-10, 1e-11, 1e-12, 1e-13, 1e-14, 1e-15, 1e-16 }
55 : : };
56 : :
57 : : // return pow(10.0,nExp) optimized for exponents in the interval [-16,16]
58 : 272455 : static double getN10Exp( int nExp )
59 : : {
60 [ + + ]: 272455 : if ( nExp < 0 )
61 : : {
62 : : // && -nExp > 0 necessary for std::numeric_limits<int>::min()
63 : : // because -nExp = nExp
64 [ + + ][ + - ]: 28513 : if ( -nExp <= n10Count && -nExp > 0 )
65 : 28483 : return n10s[1][-nExp-1];
66 : : else
67 : 30 : return pow( 10.0, static_cast<double>( nExp ) );
68 : : }
69 [ + + ]: 243942 : else if ( nExp > 0 )
70 : : {
71 [ + + ]: 229908 : if ( nExp <= n10Count )
72 : 229900 : return n10s[0][nExp-1];
73 : : else
74 : 8 : return pow( 10.0, static_cast<double>( nExp ) );
75 : : }
76 : : else // ( nExp == 0 )
77 : 272455 : return 1.0;
78 : : }
79 : :
80 : : /** Approximation algorithm for erf for 0 < x < 0.65. */
81 : 0 : void lcl_Erf0065( double x, double& fVal )
82 : : {
83 : : static const double pn[] = {
84 : : 1.12837916709551256,
85 : : 1.35894887627277916E-1,
86 : : 4.03259488531795274E-2,
87 : : 1.20339380863079457E-3,
88 : : 6.49254556481904354E-5
89 : : };
90 : : static const double qn[] = {
91 : : 1.00000000000000000,
92 : : 4.53767041780002545E-1,
93 : : 8.69936222615385890E-2,
94 : : 8.49717371168693357E-3,
95 : : 3.64915280629351082E-4
96 : : };
97 : 0 : double fPSum = 0.0;
98 : 0 : double fQSum = 0.0;
99 : 0 : double fXPow = 1.0;
100 [ # # ]: 0 : for ( unsigned int i = 0; i <= 4; ++i )
101 : : {
102 : 0 : fPSum += pn[i]*fXPow;
103 : 0 : fQSum += qn[i]*fXPow;
104 : 0 : fXPow *= x*x;
105 : : }
106 : 0 : fVal = x * fPSum / fQSum;
107 : 0 : }
108 : :
109 : : /** Approximation algorithm for erfc for 0.65 < x < 6.0. */
110 : 0 : void lcl_Erfc0600( double x, double& fVal )
111 : : {
112 : 0 : double fPSum = 0.0;
113 : 0 : double fQSum = 0.0;
114 : 0 : double fXPow = 1.0;
115 : : const double *pn;
116 : : const double *qn;
117 : :
118 [ # # ]: 0 : if ( x < 2.2 )
119 : : {
120 : : static const double pn22[] = {
121 : : 9.99999992049799098E-1,
122 : : 1.33154163936765307,
123 : : 8.78115804155881782E-1,
124 : : 3.31899559578213215E-1,
125 : : 7.14193832506776067E-2,
126 : : 7.06940843763253131E-3
127 : : };
128 : : static const double qn22[] = {
129 : : 1.00000000000000000,
130 : : 2.45992070144245533,
131 : : 2.65383972869775752,
132 : : 1.61876655543871376,
133 : : 5.94651311286481502E-1,
134 : : 1.26579413030177940E-1,
135 : : 1.25304936549413393E-2
136 : : };
137 : 0 : pn = pn22;
138 : 0 : qn = qn22;
139 : : }
140 : : else /* if ( x < 6.0 ) this is true, but the compiler does not know */
141 : : {
142 : : static const double pn60[] = {
143 : : 9.99921140009714409E-1,
144 : : 1.62356584489366647,
145 : : 1.26739901455873222,
146 : : 5.81528574177741135E-1,
147 : : 1.57289620742838702E-1,
148 : : 2.25716982919217555E-2
149 : : };
150 : : static const double qn60[] = {
151 : : 1.00000000000000000,
152 : : 2.75143870676376208,
153 : : 3.37367334657284535,
154 : : 2.38574194785344389,
155 : : 1.05074004614827206,
156 : : 2.78788439273628983E-1,
157 : : 4.00072964526861362E-2
158 : : };
159 : 0 : pn = pn60;
160 : 0 : qn = qn60;
161 : : }
162 : :
163 [ # # ]: 0 : for ( unsigned int i = 0; i < 6; ++i )
164 : : {
165 : 0 : fPSum += pn[i]*fXPow;
166 : 0 : fQSum += qn[i]*fXPow;
167 : 0 : fXPow *= x;
168 : : }
169 : 0 : fQSum += qn[6]*fXPow;
170 : 0 : fVal = exp( -1.0*x*x )* fPSum / fQSum;
171 : 0 : }
172 : :
173 : : /** Approximation algorithm for erfc for 6.0 < x < 26.54 (but used for all
174 : : x > 6.0). */
175 : 0 : void lcl_Erfc2654( double x, double& fVal )
176 : : {
177 : : static const double pn[] = {
178 : : 5.64189583547756078E-1,
179 : : 8.80253746105525775,
180 : : 3.84683103716117320E1,
181 : : 4.77209965874436377E1,
182 : : 8.08040729052301677
183 : : };
184 : : static const double qn[] = {
185 : : 1.00000000000000000,
186 : : 1.61020914205869003E1,
187 : : 7.54843505665954743E1,
188 : : 1.12123870801026015E2,
189 : : 3.73997570145040850E1
190 : : };
191 : :
192 : 0 : double fPSum = 0.0;
193 : 0 : double fQSum = 0.0;
194 : 0 : double fXPow = 1.0;
195 : :
196 [ # # ]: 0 : for ( unsigned int i = 0; i <= 4; ++i )
197 : : {
198 : 0 : fPSum += pn[i]*fXPow;
199 : 0 : fQSum += qn[i]*fXPow;
200 : 0 : fXPow /= x*x;
201 : : }
202 : 0 : fVal = exp(-1.0*x*x)*fPSum / (x*fQSum);
203 : 0 : }
204 : :
205 : : namespace {
206 : :
207 : : double const nKorrVal[] = {
208 : : 0, 9e-1, 9e-2, 9e-3, 9e-4, 9e-5, 9e-6, 9e-7, 9e-8,
209 : : 9e-9, 9e-10, 9e-11, 9e-12, 9e-13, 9e-14, 9e-15
210 : : };
211 : :
212 : : struct StringTraits
213 : : {
214 : : typedef sal_Char Char;
215 : :
216 : : typedef rtl_String String;
217 : :
218 : 4810 : static inline void createString(rtl_String ** pString,
219 : : sal_Char const * pChars, sal_Int32 nLen)
220 : : {
221 : 4810 : rtl_string_newFromStr_WithLength(pString, pChars, nLen);
222 : 4810 : }
223 : :
224 : 0 : static inline void createBuffer(rtl_String ** pBuffer,
225 : : sal_Int32 * pCapacity)
226 : : {
227 : 0 : rtl_string_new_WithLength(pBuffer, *pCapacity);
228 : 0 : }
229 : :
230 : : static inline void appendChar(rtl_String ** pBuffer, sal_Int32 * pCapacity,
231 : : sal_Int32 * pOffset, sal_Char cChar)
232 : : {
233 : : rtl_stringbuffer_insert(pBuffer, pCapacity, *pOffset, &cChar, 1);
234 : : ++*pOffset;
235 : : }
236 : :
237 : 0 : static inline void appendChars(rtl_String ** pBuffer, sal_Int32 * pCapacity,
238 : : sal_Int32 * pOffset, sal_Char const * pChars,
239 : : sal_Int32 nLen)
240 : : {
241 : 0 : rtl_stringbuffer_insert(pBuffer, pCapacity, *pOffset, pChars, nLen);
242 : 0 : *pOffset += nLen;
243 : 0 : }
244 : :
245 : 0 : static inline void appendAscii(rtl_String ** pBuffer, sal_Int32 * pCapacity,
246 : : sal_Int32 * pOffset, sal_Char const * pStr,
247 : : sal_Int32 nLen)
248 : : {
249 : 0 : rtl_stringbuffer_insert(pBuffer, pCapacity, *pOffset, pStr, nLen);
250 : 0 : *pOffset += nLen;
251 : 0 : }
252 : : };
253 : :
254 : : struct UStringTraits
255 : : {
256 : : typedef sal_Unicode Char;
257 : :
258 : : typedef rtl_uString String;
259 : :
260 : 92977 : static inline void createString(rtl_uString ** pString,
261 : : sal_Unicode const * pChars, sal_Int32 nLen)
262 : : {
263 : 92977 : rtl_uString_newFromStr_WithLength(pString, pChars, nLen);
264 : 92977 : }
265 : :
266 : 20 : static inline void createBuffer(rtl_uString ** pBuffer,
267 : : sal_Int32 * pCapacity)
268 : : {
269 : 20 : rtl_uString_new_WithLength(pBuffer, *pCapacity);
270 : 20 : }
271 : :
272 : : static inline void appendChar(rtl_uString ** pBuffer, sal_Int32 * pCapacity,
273 : : sal_Int32 * pOffset, sal_Unicode cChar)
274 : : {
275 : : rtl_uStringbuffer_insert(pBuffer, pCapacity, *pOffset, &cChar, 1);
276 : : ++*pOffset;
277 : : }
278 : :
279 : 1473 : static inline void appendChars(rtl_uString ** pBuffer,
280 : : sal_Int32 * pCapacity, sal_Int32 * pOffset,
281 : : sal_Unicode const * pChars, sal_Int32 nLen)
282 : : {
283 : 1473 : rtl_uStringbuffer_insert(pBuffer, pCapacity, *pOffset, pChars, nLen);
284 : 1473 : *pOffset += nLen;
285 : 1473 : }
286 : :
287 : 20 : static inline void appendAscii(rtl_uString ** pBuffer,
288 : : sal_Int32 * pCapacity, sal_Int32 * pOffset,
289 : : sal_Char const * pStr, sal_Int32 nLen)
290 : : {
291 : : rtl_uStringbuffer_insert_ascii(pBuffer, pCapacity, *pOffset, pStr,
292 : 20 : nLen);
293 : 20 : *pOffset += nLen;
294 : 20 : }
295 : : };
296 : :
297 : :
298 : : // Solaris C++ 5.2 compiler has problems when "StringT ** pResult" is
299 : : // "typename T::String ** pResult" instead:
300 : : template< typename T, typename StringT >
301 : 99280 : inline void doubleToString(StringT ** pResult,
302 : : sal_Int32 * pResultCapacity, sal_Int32 nResultOffset,
303 : : double fValue, rtl_math_StringFormat eFormat,
304 : : sal_Int32 nDecPlaces, typename T::Char cDecSeparator,
305 : : sal_Int32 const * pGroups,
306 : : typename T::Char cGroupSeparator,
307 : : bool bEraseTrailingDecZeros)
308 : : {
309 : : static double const nRoundVal[] = {
310 : : 5.0e+0, 0.5e+0, 0.5e-1, 0.5e-2, 0.5e-3, 0.5e-4, 0.5e-5, 0.5e-6,
311 : : 0.5e-7, 0.5e-8, 0.5e-9, 0.5e-10,0.5e-11,0.5e-12,0.5e-13,0.5e-14
312 : : };
313 : :
314 : : // sign adjustment, instead of testing for fValue<0.0 this will also fetch
315 : : // -0.0
316 : 99280 : bool bSign = rtl::math::isSignBitSet( fValue );
317 [ + + + + ]: 99280 : if( bSign )
318 : 6654 : fValue = -fValue;
319 : :
320 [ + + ][ - + ]: 99280 : if ( rtl::math::isNan( fValue ) )
321 : : {
322 : : // #i112652# XMLSchema-2
323 : 20 : sal_Int32 nCapacity = RTL_CONSTASCII_LENGTH("NaN");
324 [ + - ][ # # ]: 20 : if (pResultCapacity == 0)
325 : : {
326 : 20 : pResultCapacity = &nCapacity;
327 : 20 : T::createBuffer(pResult, pResultCapacity);
328 : 20 : nResultOffset = 0;
329 : : }
330 [ + - ][ # # ]: 20 : T::appendAscii(pResult, pResultCapacity, &nResultOffset,
331 : : RTL_CONSTASCII_STRINGPARAM("NaN"));
332 : :
333 : : return;
334 : : }
335 : :
336 : 99260 : bool bHuge = fValue == HUGE_VAL; // g++ 3.0.1 requires it this way...
337 [ + - ][ - + ]: 99260 : if ( bHuge || rtl::math::isInf( fValue ) )
[ - + ][ + - ]
[ - + ][ - + ]
338 : : {
339 : : // #i112652# XMLSchema-2
340 : 0 : sal_Int32 nCapacity = RTL_CONSTASCII_LENGTH("-INF");
341 [ # # ][ # # ]: 0 : if (pResultCapacity == 0)
342 : : {
343 : 0 : pResultCapacity = &nCapacity;
344 : 0 : T::createBuffer(pResult, pResultCapacity);
345 : 0 : nResultOffset = 0;
346 : : }
347 [ # # ][ # # ]: 0 : if ( bSign )
348 [ # # ][ # # ]: 0 : T::appendAscii(pResult, pResultCapacity, &nResultOffset,
349 : : RTL_CONSTASCII_STRINGPARAM("-"));
350 [ # # ][ # # ]: 0 : T::appendAscii(pResult, pResultCapacity, &nResultOffset,
351 : : RTL_CONSTASCII_STRINGPARAM("INF"));
352 : :
353 : : return;
354 : : }
355 : :
356 : : // find the exponent
357 : 99260 : int nExp = 0;
358 [ + + ][ + + ]: 99260 : if ( fValue > 0.0 )
359 : : {
360 : 87094 : nExp = static_cast< int >( floor( log10( fValue ) ) );
361 : 87094 : fValue /= getN10Exp( nExp );
362 : : }
363 : :
364 [ + + + ]: 99260 : switch ( eFormat )
[ - + - ]
365 : : {
366 : : case rtl_math_StringFormat_Automatic :
367 : : { // E or F depending on exponent magnitude
368 : : int nPrec;
369 [ + - ][ - + ]: 1694 : if ( nExp <= -15 || nExp >= 15 ) // #58531# was <-16, >16
[ # # ][ # # ]
370 : : {
371 : 0 : nPrec = 14;
372 : 0 : eFormat = rtl_math_StringFormat_E;
373 : : }
374 : : else
375 : : {
376 [ + - ][ # # ]: 1694 : if ( nExp < 14 )
377 : : {
378 : 1694 : nPrec = 15 - nExp - 1;
379 : 1694 : eFormat = rtl_math_StringFormat_F;
380 : : }
381 : : else
382 : : {
383 : 0 : nPrec = 15;
384 : 0 : eFormat = rtl_math_StringFormat_F;
385 : : }
386 : : }
387 [ + - ][ # # ]: 1694 : if ( nDecPlaces == rtl_math_DecimalPlaces_Max )
388 : 1694 : nDecPlaces = nPrec;
389 : : }
390 : 1694 : break;
391 : : case rtl_math_StringFormat_G :
392 : : { // G-Point, similar to sprintf %G
393 [ - + ][ - + ]: 24820 : if ( nDecPlaces == rtl_math_DecimalPlaces_DefaultSignificance )
394 : 0 : nDecPlaces = 6;
395 [ + - ][ - + ]: 24820 : if ( nExp < -4 || nExp >= nDecPlaces )
[ + - ][ - + ]
396 : : {
397 [ # # ][ # # ]: 0 : nDecPlaces = std::max< sal_Int32 >( 1, nDecPlaces - 1 );
398 : 0 : eFormat = rtl_math_StringFormat_E;
399 : : }
400 : : else
401 : : {
402 [ + - ][ + - ]: 24820 : nDecPlaces = std::max< sal_Int32 >( 0, nDecPlaces - nExp - 1 );
403 : 24820 : eFormat = rtl_math_StringFormat_F;
404 : : }
405 : : }
406 : 24820 : break;
407 : : default:
408 : 72746 : break;
409 : : }
410 : :
411 : 99260 : sal_Int32 nDigits = nDecPlaces + 1;
412 : :
413 [ + + ][ + - ]: 99260 : if( eFormat == rtl_math_StringFormat_F )
414 : 99216 : nDigits += nExp;
415 : :
416 : : // Round the number
417 [ + + ][ + - ]: 99260 : if( nDigits >= 0 )
418 : : {
419 [ - + ][ - + ]: 99250 : if( ( fValue += nRoundVal[ nDigits > 15 ? 15 : nDigits ] ) >= 10 )
420 : : {
421 : 0 : fValue = 1.0;
422 : 0 : nExp++;
423 [ # # ][ # # ]: 0 : if( eFormat == rtl_math_StringFormat_F )
424 : 0 : nDigits++;
425 : : }
426 : : }
427 : :
428 : : static sal_Int32 const nBufMax = 256;
429 : : typename T::Char aBuf[nBufMax];
430 : : typename T::Char * pBuf;
431 : : sal_Int32 nBuf = static_cast< sal_Int32 >
432 : : ( nDigits <= 0 ? std::max< sal_Int32 >( nDecPlaces, abs(nExp) )
433 [ + + ][ + - ]: 99260 : : nDigits + nDecPlaces ) + 10 + (pGroups ? abs(nDigits) * 2 : 0);
[ - + ][ + + ]
[ # # ][ - + ]
[ # # ][ - + ]
[ - + ][ # # ]
434 [ + + ][ - + ]: 99260 : if ( nBuf > nBufMax )
435 : : {
436 : 4 : pBuf = reinterpret_cast< typename T::Char * >(
437 : : rtl_allocateMemory(nBuf * sizeof (typename T::Char)));
438 : : OSL_ENSURE(pBuf != 0, "Out of memory");
439 : : }
440 : : else
441 : 99256 : pBuf = aBuf;
442 : 99260 : typename T::Char * p = pBuf;
443 [ + + ][ + + ]: 99260 : if ( bSign )
444 : 6654 : *p++ = static_cast< typename T::Char >('-');
445 : :
446 : 99260 : bool bHasDec = false;
447 : :
448 : : int nDecPos;
449 : : // Check for F format and number < 1
450 [ + + ][ + - ]: 99260 : if( eFormat == rtl_math_StringFormat_F )
451 : : {
452 [ + + ][ + + ]: 99216 : if( nExp < 0 )
453 : : {
454 : 3047 : *p++ = static_cast< typename T::Char >('0');
455 [ + - ][ + - ]: 3047 : if ( nDecPlaces > 0 )
456 : : {
457 : 3047 : *p++ = cDecSeparator;
458 : 3047 : bHasDec = true;
459 : : }
460 [ + + ][ + - ]: 3047 : sal_Int32 i = ( nDigits <= 0 ? nDecPlaces : -nExp - 1 );
461 [ + + ][ + + ]: 3936 : while( (i--) > 0 )
462 : 889 : *p++ = static_cast< typename T::Char >('0');
463 : 3047 : nDecPos = 0;
464 : : }
465 : : else
466 : 96169 : nDecPos = nExp + 1;
467 : : }
468 : : else
469 : 44 : nDecPos = 1;
470 : :
471 : 99260 : int nGrouping = 0, nGroupSelector = 0, nGroupExceed = 0;
472 [ + + ][ - + ]: 99260 : if ( nDecPos > 1 && pGroups && pGroups[0] && cGroupSeparator )
[ # # ][ # # ]
[ + + ][ - + ]
[ # # ][ # # ]
473 : : {
474 [ # # ][ # # ]: 0 : while ( nGrouping + pGroups[nGroupSelector] < nDecPos )
475 : : {
476 : 0 : nGrouping += pGroups[ nGroupSelector ];
477 [ # # ][ # # ]: 0 : if ( pGroups[nGroupSelector+1] )
478 : : {
479 [ # # ][ # # ]: 0 : if ( nGrouping + pGroups[nGroupSelector+1] >= nDecPos )
480 : 0 : break; // while
481 : 0 : ++nGroupSelector;
482 : : }
483 [ # # ][ # # ]: 0 : else if ( !nGroupExceed )
484 : 0 : nGroupExceed = nGrouping;
485 : : }
486 : : }
487 : :
488 : : // print the number
489 [ + + ][ + - ]: 99260 : if( nDigits > 0 )
490 : : {
491 : 1376972 : for ( int i = 0; ; i++ )
492 : : {
493 [ + + ][ + + ]: 1277722 : if( i < 15 )
494 : : {
495 : : int nDigit;
496 [ + + ][ + + ]: 1220966 : if (nDigits-1 == 0 && i > 0 && i < 14)
[ + + ][ + + ]
[ + - ][ + - ]
497 : 68791 : nDigit = static_cast< int >( floor( fValue
498 : : + nKorrVal[15-i] ) );
499 : : else
500 : 1152175 : nDigit = static_cast< int >( fValue + 1E-15 );
501 [ - + ][ - + ]: 1220966 : if (nDigit >= 10)
502 : : { // after-treatment of up-rounding to the next decade
503 : 0 : sal_Int32 sLen = static_cast< long >(p-pBuf)-1;
504 [ # # ][ # # ]: 0 : if (sLen == -1)
505 : : {
506 : 0 : p = pBuf;
507 [ # # ][ # # ]: 0 : if ( eFormat == rtl_math_StringFormat_F )
508 : : {
509 : 0 : *p++ = static_cast< typename T::Char >('1');
510 : 0 : *p++ = static_cast< typename T::Char >('0');
511 : : }
512 : : else
513 : : {
514 : 0 : *p++ = static_cast< typename T::Char >('1');
515 : 0 : *p++ = cDecSeparator;
516 : 0 : *p++ = static_cast< typename T::Char >('0');
517 : 0 : nExp++;
518 : 0 : bHasDec = true;
519 : : }
520 : : }
521 : : else
522 : : {
523 [ # # ][ # # ]: 0 : for (sal_Int32 j = sLen; j >= 0; j--)
524 : : {
525 : 0 : typename T::Char cS = pBuf[j];
526 [ # # ][ # # ]: 0 : if (cS != cDecSeparator)
527 : : {
528 [ # # ][ # # ]: 0 : if ( cS != static_cast< typename T::Char >('9'))
529 : : {
530 : 0 : pBuf[j] = ++cS;
531 : 0 : j = -1; // break loop
532 : : }
533 : : else
534 : : {
535 : 0 : pBuf[j]
536 : : = static_cast< typename T::Char >('0');
537 [ # # ][ # # ]: 0 : if (j == 0)
538 : : {
539 [ # # ][ # # ]: 0 : if ( eFormat == rtl_math_StringFormat_F)
540 : : { // insert '1'
541 : 0 : typename T::Char * px = p++;
542 [ # # ][ # # ]: 0 : while ( pBuf < px )
543 : : {
544 : 0 : *px = *(px-1);
545 : 0 : px--;
546 : : }
547 : 0 : pBuf[0] = static_cast<
548 : : typename T::Char >('1');
549 : : }
550 : : else
551 : : {
552 : 0 : pBuf[j] = static_cast<
553 : : typename T::Char >('1');
554 : 0 : nExp++;
555 : : }
556 : : }
557 : : }
558 : : }
559 : : }
560 : 0 : *p++ = static_cast< typename T::Char >('0');
561 : : }
562 : 0 : fValue = 0.0;
563 : : }
564 : : else
565 : : {
566 : 1220966 : *p++ = static_cast< typename T::Char >(
567 : : nDigit + static_cast< typename T::Char >('0') );
568 : 1220966 : fValue = ( fValue - nDigit ) * 10.0;
569 : : }
570 : : }
571 : : else
572 : 56756 : *p++ = static_cast< typename T::Char >('0');
573 [ + + ][ + + ]: 1277722 : if( !--nDigits )
574 : 99250 : break; // for
575 [ + + ][ + + ]: 1178472 : if( nDecPos )
576 : : {
577 [ + + ][ + + ]: 231976 : if( !--nDecPos )
578 : : {
579 : 94864 : *p++ = cDecSeparator;
580 : 94864 : bHasDec = true;
581 : : }
582 [ - + ][ - + ]: 137112 : else if ( nDecPos == nGrouping )
583 : : {
584 : 0 : *p++ = cGroupSeparator;
585 : 0 : nGrouping -= pGroups[ nGroupSelector ];
586 [ # # ][ # # ]: 0 : if ( nGroupSelector && nGrouping < nGroupExceed )
[ # # ][ # # ]
587 : 0 : --nGroupSelector;
588 : : }
589 : : }
590 : : }
591 : : }
592 : :
593 [ + + ][ + - ]: 99260 : if ( !bHasDec && eFormat == rtl_math_StringFormat_F )
[ - + ][ # # ]
594 : : { // nDecPlaces < 0 did round the value
595 [ - + ][ # # ]: 1349 : while ( --nDecPos > 0 )
596 : : { // fill before decimal point
597 [ # # ][ # # ]: 0 : if ( nDecPos == nGrouping )
598 : : {
599 : 0 : *p++ = cGroupSeparator;
600 : 0 : nGrouping -= pGroups[ nGroupSelector ];
601 [ # # ][ # # ]: 0 : if ( nGroupSelector && nGrouping < nGroupExceed )
[ # # ][ # # ]
602 : 0 : --nGroupSelector;
603 : : }
604 : 0 : *p++ = static_cast< typename T::Char >('0');
605 : : }
606 : : }
607 : :
608 [ + + ][ + + ]: 99260 : if ( bEraseTrailingDecZeros && bHasDec && p > pBuf )
[ + - ][ + - ]
[ + - ][ + - ]
609 : : {
610 [ + + ][ + + ]: 1053117 : while ( *(p-1) == static_cast< typename T::Char >('0') )
611 : 956083 : p--;
612 [ + + ][ + + ]: 97034 : if ( *(p-1) == cDecSeparator )
613 : 84790 : p--;
614 : : }
615 : :
616 : : // Print the exponent ('E', followed by '+' or '-', followed by exactly
617 : : // three digits). The code in rtl_[u]str_valueOf{Float|Double} relies on
618 : : // this format.
619 [ + + ][ - + ]: 99260 : if( eFormat == rtl_math_StringFormat_E )
620 : : {
621 [ - + ][ # # ]: 44 : if ( p == pBuf )
622 : 0 : *p++ = static_cast< typename T::Char >('1');
623 : : // maybe no nDigits if nDecPlaces < 0
624 : 44 : *p++ = static_cast< typename T::Char >('E');
625 [ - + ][ # # ]: 44 : if( nExp < 0 )
626 : : {
627 : 0 : nExp = -nExp;
628 : 0 : *p++ = static_cast< typename T::Char >('-');
629 : : }
630 : : else
631 : 44 : *p++ = static_cast< typename T::Char >('+');
632 : : // if (nExp >= 100 )
633 : 44 : *p++ = static_cast< typename T::Char >(
634 : : nExp / 100 + static_cast< typename T::Char >('0') );
635 : 44 : nExp %= 100;
636 : 44 : *p++ = static_cast< typename T::Char >(
637 : : nExp / 10 + static_cast< typename T::Char >('0') );
638 : 44 : *p++ = static_cast< typename T::Char >(
639 : : nExp % 10 + static_cast< typename T::Char >('0') );
640 : : }
641 : :
642 [ + + ][ + - ]: 99260 : if (pResultCapacity == 0)
643 : 97787 : T::createString(pResult, pBuf, p - pBuf);
644 : : else
645 [ + - ][ # # ]: 1473 : T::appendChars(pResult, pResultCapacity, &nResultOffset, pBuf,
646 : : p - pBuf);
647 : :
648 [ + + ][ - + ]: 99260 : if ( pBuf != &aBuf[0] )
649 : 99280 : rtl_freeMemory(pBuf);
650 : : }
651 : :
652 : : }
653 : :
654 : 4810 : void SAL_CALL rtl_math_doubleToString(rtl_String ** pResult,
655 : : sal_Int32 * pResultCapacity,
656 : : sal_Int32 nResultOffset, double fValue,
657 : : rtl_math_StringFormat eFormat,
658 : : sal_Int32 nDecPlaces,
659 : : sal_Char cDecSeparator,
660 : : sal_Int32 const * pGroups,
661 : : sal_Char cGroupSeparator,
662 : : sal_Bool bEraseTrailingDecZeros)
663 : : SAL_THROW_EXTERN_C()
664 : : {
665 : : doubleToString< StringTraits, StringTraits::String >(
666 : : pResult, pResultCapacity, nResultOffset, fValue, eFormat, nDecPlaces,
667 : 4810 : cDecSeparator, pGroups, cGroupSeparator, bEraseTrailingDecZeros);
668 : 4810 : }
669 : :
670 : 94470 : void SAL_CALL rtl_math_doubleToUString(rtl_uString ** pResult,
671 : : sal_Int32 * pResultCapacity,
672 : : sal_Int32 nResultOffset, double fValue,
673 : : rtl_math_StringFormat eFormat,
674 : : sal_Int32 nDecPlaces,
675 : : sal_Unicode cDecSeparator,
676 : : sal_Int32 const * pGroups,
677 : : sal_Unicode cGroupSeparator,
678 : : sal_Bool bEraseTrailingDecZeros)
679 : : SAL_THROW_EXTERN_C()
680 : : {
681 : : doubleToString< UStringTraits, UStringTraits::String >(
682 : : pResult, pResultCapacity, nResultOffset, fValue, eFormat, nDecPlaces,
683 : 94470 : cDecSeparator, pGroups, cGroupSeparator, bEraseTrailingDecZeros);
684 : 94470 : }
685 : :
686 : :
687 : : namespace {
688 : :
689 : : // if nExp * 10 + nAdd would result in overflow
690 : 65 : inline bool long10Overflow( long& nExp, int nAdd )
691 : : {
692 [ + - ][ - + ]: 65 : if ( nExp > (LONG_MAX/10)
[ # # ]
693 : : || (nExp == (LONG_MAX/10) && nAdd > (LONG_MAX%10)) )
694 : : {
695 : 0 : nExp = LONG_MAX;
696 : 0 : return true;
697 : : }
698 : 65 : return false;
699 : : }
700 : :
701 : : // We are only concerned about ASCII arabic numerical digits here
702 : : template< typename CharT >
703 : 316656 : inline bool isDigit( CharT c )
704 : : {
705 [ + + ][ + + ]: 316656 : return 0x30 <= c && c <= 0x39;
[ + + ][ + - ]
706 : : }
707 : :
708 : : template< typename CharT >
709 : 110451 : inline double stringToDouble(CharT const * pBegin, CharT const * pEnd,
710 : : CharT cDecSeparator, CharT cGroupSeparator,
711 : : rtl_math_ConversionStatus * pStatus,
712 : : CharT const ** pParsedEnd)
713 : : {
714 : 110451 : double fVal = 0.0;
715 : 110451 : rtl_math_ConversionStatus eStatus = rtl_math_ConversionStatus_Ok;
716 : :
717 : 110451 : CharT const * p0 = pBegin;
718 [ + + ][ + + ]: 110471 : while (p0 != pEnd && (*p0 == CharT(' ') || *p0 == CharT('\t')))
[ - + ][ + + ]
[ + - ][ + - ]
[ - + ][ - + ]
719 : 20 : ++p0;
720 : : bool bSign;
721 [ + + ][ + + ]: 110451 : if (p0 != pEnd && *p0 == CharT('-'))
[ + - ][ + + ]
722 : : {
723 : 15238 : bSign = true;
724 : 15238 : ++p0;
725 : : }
726 : : else
727 : : {
728 : 95213 : bSign = false;
729 [ + + ][ + + ]: 95213 : if (p0 != pEnd && *p0 == CharT('+'))
[ + - ][ - + ]
730 : 10 : ++p0;
731 : : }
732 : 110451 : CharT const * p = p0;
733 : 110451 : bool bDone = false;
734 : :
735 : : // #i112652# XMLSchema-2
736 [ + + ][ + + ]: 110451 : if (3 >= (pEnd - p))
737 : : {
738 [ - + ][ # # ]: 63470 : if ((CharT('N') == p[0]) && (CharT('a') == p[1])
[ # # ][ - + ]
[ # # ][ # # ]
739 : : && (CharT('N') == p[2]))
740 : : {
741 : 0 : p += 3;
742 : 0 : rtl::math::setNan( &fVal );
743 : 0 : bDone = true;
744 : : }
745 [ - + ][ # # ]: 63470 : else if ((CharT('I') == p[0]) && (CharT('N') == p[1])
[ # # ][ - + ]
[ # # ][ # # ]
746 : : && (CharT('F') == p[2]))
747 : : {
748 : 0 : p += 3;
749 : 0 : fVal = HUGE_VAL;
750 : 0 : eStatus = rtl_math_ConversionStatus_OutOfRange;
751 : 63470 : bDone = true;
752 : : }
753 : : }
754 : :
755 [ + - ][ + - ]: 110451 : if (!bDone) // do not recognize e.g. NaN1.23
756 : : {
757 : : // leading zeros and group separators may be safely ignored
758 [ + + ][ + + ]: 144013 : while (p != pEnd && (*p == CharT('0') || *p == cGroupSeparator))
[ + + ][ + + ]
[ + + ][ + + ]
[ - + ][ + + ]
759 : 33562 : ++p;
760 : :
761 : 110451 : long nValExp = 0; // carry along exponent of mantissa
762 : :
763 : : // integer part of mantissa
764 [ + + ][ + + ]: 298302 : for (; p != pEnd; ++p)
765 : : {
766 : 229881 : CharT c = *p;
767 [ + + ][ + + ]: 229881 : if (isDigit(c))
768 : : {
769 : 187821 : fVal = fVal * 10.0 + static_cast< double >( c - CharT('0') );
770 : 187821 : ++nValExp;
771 : : }
772 [ + + ][ + - ]: 42060 : else if (c != cGroupSeparator)
773 : 42030 : break;
774 : : }
775 : :
776 : : // fraction part of mantissa
777 [ + + ][ + + ]: 110451 : if (p != pEnd && *p == cDecSeparator)
[ + + ][ + - ]
778 : : {
779 : 32123 : ++p;
780 : 32123 : double fFrac = 0.0;
781 : 32123 : long nFracExp = 0;
782 [ + + ][ + + ]: 84229 : while (p != pEnd && *p == CharT('0'))
[ + + ][ + + ]
[ + + ][ + + ]
783 : : {
784 : 52106 : --nFracExp;
785 : 52106 : ++p;
786 : : }
787 [ + + ][ + + ]: 32123 : if ( nValExp == 0 )
788 : 8687 : nValExp = nFracExp - 1; // no integer part => fraction exponent
789 : : // one decimal digit needs ld(10) ~= 3.32 bits
790 : : static const int nSigs = (DBL_MANT_DIG / 3) + 1;
791 : 32123 : int nDigs = 0;
792 [ + + ][ + + ]: 118309 : for (; p != pEnd; ++p)
793 : : {
794 : 86705 : CharT c = *p;
795 [ + + ][ - + ]: 86705 : if (!isDigit(c))
796 : 519 : break;
797 [ + - ][ + - ]: 86186 : if ( nDigs < nSigs )
798 : : { // further digits (more than nSigs) don't have any
799 : : // significance
800 : 86186 : fFrac = fFrac * 10.0 + static_cast<double>(c - CharT('0'));
801 : 86186 : --nFracExp;
802 : 86186 : ++nDigs;
803 : : }
804 : : }
805 [ + + ][ + + ]: 32123 : if ( fFrac != 0.0 )
806 : 23185 : fVal += rtl::math::pow10Exp( fFrac, nFracExp );
807 [ + + ][ + + ]: 8938 : else if ( nValExp < 0 )
808 : 4329 : nValExp = 0; // no digit other than 0 after decimal point
809 : : }
810 : :
811 [ + + ][ + + ]: 110451 : if ( nValExp > 0 )
812 : 80114 : --nValExp; // started with offset +1 at the first mantissa digit
813 : :
814 : : // Exponent
815 [ + + ][ + + ]: 110451 : if (p != p0 && p != pEnd && (*p == CharT('E') || *p == CharT('e')))
[ + + ][ - + ]
[ + - ][ - + ]
[ # # ][ # # ]
816 : : {
817 : 45 : ++p;
818 : : bool bExpSign;
819 [ + + ][ + + ]: 45 : if (p != pEnd && *p == CharT('-'))
[ # # ][ # # ]
820 : : {
821 : 20 : bExpSign = true;
822 : 20 : ++p;
823 : : }
824 : : else
825 : : {
826 : 25 : bExpSign = false;
827 [ + + ][ - + ]: 25 : if (p != pEnd && *p == CharT('+'))
[ # # ][ # # ]
828 : 0 : ++p;
829 : : }
830 [ - + ][ # # ]: 45 : if ( fVal == 0.0 )
831 : : { // no matter what follows, zero stays zero, but carry on the
832 : : // offset
833 [ # # ][ # # ]: 0 : while (p != pEnd && isDigit(*p))
[ # # ][ # # ]
[ # # ][ # # ]
834 : 0 : ++p;
835 : : }
836 : : else
837 : : {
838 : 45 : bool bOverFlow = false;
839 : 45 : long nExp = 0;
840 [ + + ][ # # ]: 110 : for (; p != pEnd; ++p)
841 : : {
842 : 70 : CharT c = *p;
843 [ + + ][ # # ]: 70 : if (!isDigit(c))
844 : 5 : break;
845 : 65 : int i = c - CharT('0');
846 [ - + ][ # # ]: 65 : if ( long10Overflow( nExp, i ) )
847 : 0 : bOverFlow = true;
848 : : else
849 : 65 : nExp = nExp * 10 + i;
850 : : }
851 [ + + ][ # # ]: 45 : if ( nExp )
852 : : {
853 [ + + ][ # # ]: 40 : if ( bExpSign )
854 : 20 : nExp = -nExp;
855 [ - + ][ # # ]: 40 : long nAllExp = ( bOverFlow ? 0 : nExp + nValExp );
856 [ + - ][ - + ]: 40 : if ( nAllExp > DBL_MAX_10_EXP || (bOverFlow && !bExpSign) )
[ # # ][ # # ]
[ # # ][ # # ]
857 : : { // overflow
858 : 0 : fVal = HUGE_VAL;
859 : 0 : eStatus = rtl_math_ConversionStatus_OutOfRange;
860 : : }
861 [ + - ][ - + ]: 40 : else if ((nAllExp < DBL_MIN_10_EXP) ||
[ # # ][ # # ]
[ # # ][ # # ]
862 : : (bOverFlow && bExpSign) )
863 : : { // underflow
864 : 0 : fVal = 0.0;
865 : 0 : eStatus = rtl_math_ConversionStatus_OutOfRange;
866 : : }
867 [ + - ][ - + ]: 40 : else if ( nExp > DBL_MAX_10_EXP || nExp < DBL_MIN_10_EXP )
[ # # ][ # # ]
868 : : { // compensate exponents
869 : 0 : fVal = rtl::math::pow10Exp( fVal, -nValExp );
870 : 0 : fVal = rtl::math::pow10Exp( fVal, nAllExp );
871 : : }
872 : : else
873 : 45 : fVal = rtl::math::pow10Exp( fVal, nExp ); // normal
874 : : }
875 : : }
876 : : }
877 [ + + ][ + + ]: 110406 : else if (p - p0 == 2 && p != pEnd && p[0] == CharT('#')
[ - + ][ # # ]
[ # # ][ + + ]
[ - + ][ # # ]
[ # # ][ # # ]
878 : : && p[-1] == cDecSeparator && p[-2] == CharT('1'))
879 : : {
880 [ # # ][ # # ]: 0 : if (pEnd - p >= 4 && p[1] == CharT('I') && p[2] == CharT('N')
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
881 : : && p[3] == CharT('F'))
882 : : {
883 : : // "1.#INF", "+1.#INF", "-1.#INF"
884 : 0 : p += 4;
885 : 0 : fVal = HUGE_VAL;
886 : 0 : eStatus = rtl_math_ConversionStatus_OutOfRange;
887 : : // Eat any further digits:
888 [ # # ][ # # ]: 0 : while (p != pEnd && isDigit(*p))
[ # # ][ # # ]
[ # # ][ # # ]
889 : 0 : ++p;
890 : : }
891 [ # # ][ # # ]: 0 : else if (pEnd - p >= 4 && p[1] == CharT('N') && p[2] == CharT('A')
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
892 : : && p[3] == CharT('N'))
893 : : {
894 : : // "1.#NAN", "+1.#NAN", "-1.#NAN"
895 : 0 : p += 4;
896 : 0 : rtl::math::setNan( &fVal );
897 [ # # # # ]: 0 : if (bSign)
898 : : {
899 : : union {
900 : : double sd;
901 : : sal_math_Double md;
902 : : } m;
903 : 0 : m.sd = fVal;
904 : 0 : m.md.w32_parts.msw |= 0x80000000; // create negative NaN
905 : 0 : fVal = m.sd;
906 : 0 : bSign = false; // don't negate again
907 : : }
908 : : // Eat any further digits:
909 [ # # ][ # # ]: 110451 : while (p != pEnd && isDigit(*p))
[ # # ][ # # ]
[ # # ][ # # ]
910 : 0 : ++p;
911 : : }
912 : : }
913 : : }
914 : :
915 : : // overflow also if more than DBL_MAX_10_EXP digits without decimal
916 : : // separator, or 0. and more than DBL_MIN_10_EXP digits, ...
917 : 110451 : bool bHuge = fVal == HUGE_VAL; // g++ 3.0.1 requires it this way...
918 [ - + ][ - + ]: 110451 : if ( bHuge )
919 : 0 : eStatus = rtl_math_ConversionStatus_OutOfRange;
920 : :
921 [ + + ][ + + ]: 110451 : if ( bSign )
922 : 15238 : fVal = -fVal;
923 : :
924 [ + + ][ - + ]: 110451 : if (pStatus != 0)
925 : 45765 : *pStatus = eStatus;
926 [ + + ][ + + ]: 110451 : if (pParsedEnd != 0)
927 [ + + ][ - + ]: 51345 : *pParsedEnd = p == p0 ? pBegin : p;
928 : :
929 : 110451 : return fVal;
930 : : }
931 : :
932 : : }
933 : :
934 : 38821 : double SAL_CALL rtl_math_stringToDouble(sal_Char const * pBegin,
935 : : sal_Char const * pEnd,
936 : : sal_Char cDecSeparator,
937 : : sal_Char cGroupSeparator,
938 : : rtl_math_ConversionStatus * pStatus,
939 : : sal_Char const ** pParsedEnd)
940 : : SAL_THROW_EXTERN_C()
941 : : {
942 : : return stringToDouble(pBegin, pEnd, cDecSeparator, cGroupSeparator, pStatus,
943 : 38821 : pParsedEnd);
944 : : }
945 : :
946 : 71630 : double SAL_CALL rtl_math_uStringToDouble(sal_Unicode const * pBegin,
947 : : sal_Unicode const * pEnd,
948 : : sal_Unicode cDecSeparator,
949 : : sal_Unicode cGroupSeparator,
950 : : rtl_math_ConversionStatus * pStatus,
951 : : sal_Unicode const ** pParsedEnd)
952 : : SAL_THROW_EXTERN_C()
953 : : {
954 : : return stringToDouble(pBegin, pEnd, cDecSeparator, cGroupSeparator, pStatus,
955 : 71630 : pParsedEnd);
956 : : }
957 : :
958 : 49529421 : double SAL_CALL rtl_math_round(double fValue, int nDecPlaces,
959 : : enum rtl_math_RoundingMode eMode)
960 : : SAL_THROW_EXTERN_C()
961 : : {
962 : : OSL_ASSERT(nDecPlaces >= -20 && nDecPlaces <= 20);
963 : :
964 [ + + ]: 49529421 : if ( fValue == 0.0 )
965 : 49365308 : return fValue;
966 : :
967 : : // sign adjustment
968 : 164113 : bool bSign = rtl::math::isSignBitSet( fValue );
969 [ + + ]: 164113 : if ( bSign )
970 : 1882 : fValue = -fValue;
971 : :
972 : 164113 : double fFac = 0;
973 [ + + ]: 164113 : if ( nDecPlaces != 0 )
974 : : {
975 : : // max 20 decimals, we don't have unlimited precision
976 : : // #38810# and no overflow on fValue*=fFac
977 [ + - ][ + - ]: 121820 : if ( nDecPlaces < -20 || 20 < nDecPlaces || fValue > (DBL_MAX / 1e20) )
[ - + ]
978 [ # # ]: 0 : return bSign ? -fValue : fValue;
979 : :
980 : 121820 : fFac = getN10Exp( nDecPlaces );
981 : 121820 : fValue *= fFac;
982 : : }
983 : : //else //! uninitialized fFac, not needed
984 : :
985 [ + - - + : 164113 : switch ( eMode )
- - - -
- ]
986 : : {
987 : : case rtl_math_RoundingMode_Corrected :
988 : : {
989 : : int nExp; // exponent for correction
990 [ + - ]: 162623 : if ( fValue > 0.0 )
991 : 162623 : nExp = static_cast<int>( floor( log10( fValue ) ) );
992 : : else
993 : 0 : nExp = 0;
994 : 162623 : int nIndex = 15 - nExp;
995 [ + + ]: 162623 : if ( nIndex > 15 )
996 : 11 : nIndex = 15;
997 [ + + ]: 162612 : else if ( nIndex <= 1 )
998 : 37030 : nIndex = 0;
999 : 162623 : fValue = floor( fValue + 0.5 + nKorrVal[nIndex] );
1000 : : }
1001 : 162623 : break;
1002 : : case rtl_math_RoundingMode_Down :
1003 : 0 : fValue = rtl::math::approxFloor( fValue );
1004 : 0 : break;
1005 : : case rtl_math_RoundingMode_Up :
1006 : 0 : fValue = rtl::math::approxCeil( fValue );
1007 : 0 : break;
1008 : : case rtl_math_RoundingMode_Floor :
1009 : : fValue = bSign ? rtl::math::approxCeil( fValue )
1010 [ - + ]: 1490 : : rtl::math::approxFloor( fValue );
1011 : 1490 : break;
1012 : : case rtl_math_RoundingMode_Ceiling :
1013 : : fValue = bSign ? rtl::math::approxFloor( fValue )
1014 [ # # ]: 0 : : rtl::math::approxCeil( fValue );
1015 : 0 : break;
1016 : : case rtl_math_RoundingMode_HalfDown :
1017 : : {
1018 : 0 : double f = floor( fValue );
1019 [ # # ]: 0 : fValue = ((fValue - f) <= 0.5) ? f : ceil( fValue );
1020 : : }
1021 : 0 : break;
1022 : : case rtl_math_RoundingMode_HalfUp :
1023 : : {
1024 : 0 : double f = floor( fValue );
1025 [ # # ]: 0 : fValue = ((fValue - f) < 0.5) ? f : ceil( fValue );
1026 : : }
1027 : 0 : break;
1028 : : case rtl_math_RoundingMode_HalfEven :
1029 : : #if defined FLT_ROUNDS
1030 : : /*
1031 : : Use fast version. FLT_ROUNDS may be defined to a function by some compilers!
1032 : :
1033 : : DBL_EPSILON is the smallest fractional number which can be represented,
1034 : : its reciprocal is therefore the smallest number that cannot have a
1035 : : fractional part. Once you add this reciprocal to `x', its fractional part
1036 : : is stripped off. Simply subtracting the reciprocal back out returns `x'
1037 : : without its fractional component.
1038 : : Simple, clever, and elegant - thanks to Ross Cottrell, the original author,
1039 : : who placed it into public domain.
1040 : :
1041 : : volatile: prevent compiler from being too smart
1042 : : */
1043 : : if ( FLT_ROUNDS == 1 )
1044 : : {
1045 : 0 : volatile double x = fValue + 1.0 / DBL_EPSILON;
1046 : 0 : fValue = x - 1.0 / DBL_EPSILON;
1047 : : }
1048 : : else
1049 : : #endif // FLT_ROUNDS
1050 : : {
1051 : : double f = floor( fValue );
1052 : : if ( (fValue - f) != 0.5 )
1053 : : fValue = floor( fValue + 0.5 );
1054 : : else
1055 : : {
1056 : : double g = f / 2.0;
1057 : : fValue = (g == floor( g )) ? f : (f + 1.0);
1058 : : }
1059 : : }
1060 : 0 : break;
1061 : : default:
1062 : : OSL_ASSERT(false);
1063 : 0 : break;
1064 : : }
1065 : :
1066 [ + + ]: 164113 : if ( nDecPlaces != 0 )
1067 : 121820 : fValue /= fFac;
1068 : :
1069 [ + + ]: 49529421 : return bSign ? -fValue : fValue;
1070 : : }
1071 : :
1072 : :
1073 : 26511 : double SAL_CALL rtl_math_pow10Exp(double fValue, int nExp) SAL_THROW_EXTERN_C()
1074 : : {
1075 : 26511 : return fValue * getN10Exp( nExp );
1076 : : }
1077 : :
1078 : :
1079 : 38891 : double SAL_CALL rtl_math_approxValue( double fValue ) SAL_THROW_EXTERN_C()
1080 : : {
1081 [ + + ][ + - ]: 38891 : if (fValue == 0.0 || fValue == HUGE_VAL || !::rtl::math::isFinite( fValue))
[ - + ][ + + ]
1082 : : // We don't handle these conditions. Bail out.
1083 : 1861 : return fValue;
1084 : :
1085 : 37030 : double fOrigValue = fValue;
1086 : :
1087 : 37030 : bool bSign = ::rtl::math::isSignBitSet( fValue);
1088 [ + + ]: 37030 : if (bSign)
1089 : 2024 : fValue = -fValue;
1090 : :
1091 : 37030 : int nExp = static_cast<int>( floor( log10( fValue)));
1092 : 37030 : nExp = 14 - nExp;
1093 : 37030 : double fExpValue = getN10Exp( nExp);
1094 : :
1095 : 37030 : fValue *= fExpValue;
1096 : : // If the original value was near DBL_MIN we got an overflow. Restore and
1097 : : // bail out.
1098 [ - + ]: 37030 : if (!rtl::math::isFinite( fValue))
1099 : 0 : return fOrigValue;
1100 : 37030 : fValue = rtl_math_round( fValue, 0, rtl_math_RoundingMode_Corrected);
1101 : 37030 : fValue /= fExpValue;
1102 : : // If the original value was near DBL_MAX we got an overflow. Restore and
1103 : : // bail out.
1104 [ - + ]: 37030 : if (!rtl::math::isFinite( fValue))
1105 : 0 : return fOrigValue;
1106 : :
1107 [ + + ]: 38891 : return bSign ? -fValue : fValue;
1108 : : }
1109 : :
1110 : :
1111 : 0 : double SAL_CALL rtl_math_expm1( double fValue ) SAL_THROW_EXTERN_C()
1112 : : {
1113 : 0 : double fe = exp( fValue );
1114 [ # # ]: 0 : if (fe == 1.0)
1115 : 0 : return fValue;
1116 [ # # ]: 0 : if (fe-1.0 == -1.0)
1117 : 0 : return -1.0;
1118 : 0 : return (fe-1.0) * fValue / log(fe);
1119 : : }
1120 : :
1121 : :
1122 : 0 : double SAL_CALL rtl_math_log1p( double fValue ) SAL_THROW_EXTERN_C()
1123 : : {
1124 : : // Use volatile because a compiler may be too smart "optimizing" the
1125 : : // condition such that in certain cases the else path was called even if
1126 : : // (fp==1.0) was true, where the term (fp-1.0) then resulted in 0.0 and
1127 : : // hence the entire expression resulted in NaN.
1128 : : // Happened with g++ 3.4.1 and an input value of 9.87E-18
1129 : 0 : volatile double fp = 1.0 + fValue;
1130 [ # # ]: 0 : if (fp == 1.0)
1131 : 0 : return fValue;
1132 : : else
1133 : 0 : return log(fp) * fValue / (fp-1.0);
1134 : : }
1135 : :
1136 : :
1137 : 0 : double SAL_CALL rtl_math_atanh( double fValue ) SAL_THROW_EXTERN_C()
1138 : : {
1139 : 0 : return 0.5 * rtl_math_log1p( 2.0 * fValue / (1.0-fValue) );
1140 : : }
1141 : :
1142 : :
1143 : : /** Parent error function (erf) that calls different algorithms based on the
1144 : : value of x. It takes care of cases where x is negative as erf is an odd
1145 : : function i.e. erf(-x) = -erf(x).
1146 : :
1147 : : Kramer, W., and Blomquist, F., 2000, Algorithms with Guaranteed Error Bounds
1148 : : for the Error Function and the Complementary Error Function
1149 : :
1150 : : http://www.math.uni-wuppertal.de/wrswt/literatur_en.html
1151 : :
1152 : : @author Kohei Yoshida <kohei@openoffice.org>
1153 : :
1154 : : @see #i55735#
1155 : : */
1156 : 0 : double SAL_CALL rtl_math_erf( double x ) SAL_THROW_EXTERN_C()
1157 : : {
1158 [ # # ]: 0 : if( x == 0.0 )
1159 : 0 : return 0.0;
1160 : :
1161 : 0 : bool bNegative = false;
1162 [ # # ]: 0 : if ( x < 0.0 )
1163 : : {
1164 : 0 : x = fabs( x );
1165 : 0 : bNegative = true;
1166 : : }
1167 : :
1168 : 0 : double fErf = 1.0;
1169 [ # # ]: 0 : if ( x < 1.0e-10 )
1170 : 0 : fErf = (double) (x*1.1283791670955125738961589031215452L);
1171 [ # # ]: 0 : else if ( x < 0.65 )
1172 : 0 : lcl_Erf0065( x, fErf );
1173 : : else
1174 : 0 : fErf = 1.0 - rtl_math_erfc( x );
1175 : :
1176 [ # # ]: 0 : if ( bNegative )
1177 : 0 : fErf *= -1.0;
1178 : :
1179 : 0 : return fErf;
1180 : : }
1181 : :
1182 : :
1183 : : /** Parent complementary error function (erfc) that calls different algorithms
1184 : : based on the value of x. It takes care of cases where x is negative as erfc
1185 : : satisfies relationship erfc(-x) = 2 - erfc(x). See the comment for Erf(x)
1186 : : for the source publication.
1187 : :
1188 : : @author Kohei Yoshida <kohei@openoffice.org>
1189 : :
1190 : : @see #i55735#, moved from module scaddins (#i97091#)
1191 : :
1192 : : */
1193 : 0 : double SAL_CALL rtl_math_erfc( double x ) SAL_THROW_EXTERN_C()
1194 : : {
1195 [ # # ]: 0 : if ( x == 0.0 )
1196 : 0 : return 1.0;
1197 : :
1198 : 0 : bool bNegative = false;
1199 [ # # ]: 0 : if ( x < 0.0 )
1200 : : {
1201 : 0 : x = fabs( x );
1202 : 0 : bNegative = true;
1203 : : }
1204 : :
1205 : 0 : double fErfc = 0.0;
1206 [ # # ]: 0 : if ( x >= 0.65 )
1207 : : {
1208 [ # # ]: 0 : if ( x < 6.0 )
1209 : 0 : lcl_Erfc0600( x, fErfc );
1210 : : else
1211 : 0 : lcl_Erfc2654( x, fErfc );
1212 : : }
1213 : : else
1214 : 0 : fErfc = 1.0 - rtl_math_erf( x );
1215 : :
1216 [ # # ]: 0 : if ( bNegative )
1217 : 0 : fErfc = 2.0 - fErfc;
1218 : :
1219 : 0 : return fErfc;
1220 : : }
1221 : :
1222 : : /** improved accuracy of asinh for |x| large and for x near zero
1223 : : @see #i97605#
1224 : : */
1225 : 0 : double SAL_CALL rtl_math_asinh( double fX ) SAL_THROW_EXTERN_C()
1226 : : {
1227 [ # # ]: 0 : if ( fX == 0.0 )
1228 : 0 : return 0.0;
1229 : : else
1230 : : {
1231 : 0 : double fSign = 1.0;
1232 [ # # ]: 0 : if ( fX < 0.0 )
1233 : : {
1234 : 0 : fX = - fX;
1235 : 0 : fSign = -1.0;
1236 : : }
1237 [ # # ]: 0 : if ( fX < 0.125 )
1238 : 0 : return fSign * rtl_math_log1p( fX + fX*fX / (1.0 + sqrt( 1.0 + fX*fX)));
1239 [ # # ]: 0 : else if ( fX < 1.25e7 )
1240 : 0 : return fSign * log( fX + sqrt( 1.0 + fX*fX));
1241 : : else
1242 : 0 : return fSign * log( 2.0*fX);
1243 : : }
1244 : : }
1245 : :
1246 : : /** improved accuracy of acosh for x large and for x near 1
1247 : : @see #i97605#
1248 : : */
1249 : 5 : double SAL_CALL rtl_math_acosh( double fX ) SAL_THROW_EXTERN_C()
1250 : : {
1251 : 5 : volatile double fZ = fX - 1.0;
1252 [ - + ]: 5 : if ( fX < 1.0 )
1253 : : {
1254 : : double fResult;
1255 : 0 : ::rtl::math::setNan( &fResult );
1256 : 0 : return fResult;
1257 : : }
1258 [ + - ]: 5 : else if ( fX == 1.0 )
1259 : 5 : return 0.0;
1260 [ # # ]: 0 : else if ( fX < 1.1 )
1261 : 0 : return rtl_math_log1p( fZ + sqrt( fZ*fZ + 2.0*fZ));
1262 [ # # ]: 0 : else if ( fX < 1.25e7 )
1263 : 0 : return log( fX + sqrt( fX*fX - 1.0));
1264 : : else
1265 : 5 : return log( 2.0*fX);
1266 : : }
1267 : :
1268 : : /* vim:set shiftwidth=4 softtabstop=4 expandtab: */
|
