LCOV - code coverage report
Current view: top level - sal/rtl - math.cxx (source / functions) Hit Total Coverage
Test: commit e02a6cb2c3e2b23b203b422e4e0680877f232636 Lines: 62 531 11.7 %
Date: 2014-04-14 Functions: 4 31 12.9 %
Legend: Lines: hit not hit

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

Generated by: LCOV version 1.10