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