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