xref: /trunk/main/sal/rtl/source/math.cxx (revision 63b47d7f)
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]
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. */
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. */
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). */
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 
213     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 
219     static inline void createBuffer(rtl_String ** pBuffer,
220                                     sal_Int32 * pCapacity)
221     {
222         rtl_string_new_WithLength(pBuffer, *pCapacity);
223     }
224 
225     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 
232     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 
240     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 
255     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 
261     static inline void createBuffer(rtl_uString ** pBuffer,
262                                     sal_Int32 * pCapacity)
263     {
264         rtl_uString_new_WithLength(pBuffer, *pCapacity);
265     }
266 
267     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 
274     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 
282     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 >
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 
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 
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
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 >
698 inline bool isDigit( CharT c )
699 {
700     return 0x30 <= c && c <= 0x39;
701 }
702 
703 template< typename CharT >
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 
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 
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 
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 
1068 double SAL_CALL rtl_math_pow10Exp(double fValue, int nExp) SAL_THROW_EXTERN_C()
1069 {
1070     return fValue * getN10Exp( nExp );
1071 }
1072 
1073 
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 
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 
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 
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 
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  */
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  */
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  */
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  */
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