diff options
Diffstat (limited to 'sal/rtl/source/math.cxx')
| -rw-r--r-- | sal/rtl/source/math.cxx | 250 |
1 files changed, 250 insertions, 0 deletions
diff --git a/sal/rtl/source/math.cxx b/sal/rtl/source/math.cxx index 3f41d5a9c105..a255ca21b13a 100644 --- a/sal/rtl/source/math.cxx +++ b/sal/rtl/source/math.cxx @@ -79,6 +79,130 @@ static double getN10Exp( int nExp ) return 1.0; } +/** Approximation algorithm for erf for 0 < x < 0.65. */ +void lcl_Erf0065( double x, double& fVal ) +{ + static const double pn[] = { + 1.12837916709551256, + 1.35894887627277916E-1, + 4.03259488531795274E-2, + 1.20339380863079457E-3, + 6.49254556481904354E-5 + }; + static const double qn[] = { + 1.00000000000000000, + 4.53767041780002545E-1, + 8.69936222615385890E-2, + 8.49717371168693357E-3, + 3.64915280629351082E-4 + }; + double fPSum = 0.0; + double fQSum = 0.0; + double fXPow = 1.0; + for ( unsigned int i = 0; i <= 4; ++i ) + { + fPSum += pn[i]*fXPow; + fQSum += qn[i]*fXPow; + fXPow *= x*x; + } + fVal = x * fPSum / fQSum; +} + +/** Approximation algorithm for erfc for 0.65 < x < 6.0. */ +void lcl_Erfc0600( double x, double& fVal ) +{ + double fPSum = 0.0; + double fQSum = 0.0; + double fXPow = 1.0; + const double *pn; + const double *qn; + + if ( x < 2.2 ) + { + static const double pn22[] = { + 9.99999992049799098E-1, + 1.33154163936765307, + 8.78115804155881782E-1, + 3.31899559578213215E-1, + 7.14193832506776067E-2, + 7.06940843763253131E-3 + }; + static const double qn22[] = { + 1.00000000000000000, + 2.45992070144245533, + 2.65383972869775752, + 1.61876655543871376, + 5.94651311286481502E-1, + 1.26579413030177940E-1, + 1.25304936549413393E-2 + }; + pn = pn22; + qn = qn22; + } + else /* if ( x < 6.0 ) this is true, but the compiler does not know */ + { + static const double pn60[] = { + 9.99921140009714409E-1, + 1.62356584489366647, + 1.26739901455873222, + 5.81528574177741135E-1, + 1.57289620742838702E-1, + 2.25716982919217555E-2 + }; + static const double qn60[] = { + 1.00000000000000000, + 2.75143870676376208, + 3.37367334657284535, + 2.38574194785344389, + 1.05074004614827206, + 2.78788439273628983E-1, + 4.00072964526861362E-2 + }; + pn = pn60; + qn = qn60; + } + + for ( unsigned int i = 0; i < 6; ++i ) + { + fPSum += pn[i]*fXPow; + fQSum += qn[i]*fXPow; + fXPow *= x; + } + fQSum += qn[6]*fXPow; + fVal = exp( -1.0*x*x )* fPSum / fQSum; +} + +/** Approximation algorithm for erfc for 6.0 < x < 26.54 (but used for all + x > 6.0). */ +void lcl_Erfc2654( double x, double& fVal ) +{ + static const double pn[] = { + 5.64189583547756078E-1, + 8.80253746105525775, + 3.84683103716117320E1, + 4.77209965874436377E1, + 8.08040729052301677 + }; + static const double qn[] = { + 1.00000000000000000, + 1.61020914205869003E1, + 7.54843505665954743E1, + 1.12123870801026015E2, + 3.73997570145040850E1 + }; + + double fPSum = 0.0; + double fQSum = 0.0; + double fXPow = 1.0; + + for ( unsigned int i = 0; i <= 4; ++i ) + { + fPSum += pn[i]*fXPow; + fQSum += qn[i]*fXPow; + fXPow /= x*x; + } + fVal = exp(-1.0*x*x)*fPSum / (x*fQSum); +} namespace { @@ -991,3 +1115,129 @@ double SAL_CALL rtl_math_atanh( double fValue ) SAL_THROW_EXTERN_C() { return 0.5 * rtl_math_log1p( 2.0 * fValue / (1.0-fValue) ); } + + +/** Parent error function (erf) that calls different algorithms based on the + value of x. It takes care of cases where x is negative as erf is an odd + function i.e. erf(-x) = -erf(x). + + Kramer, W., and Blomquist, F., 2000, Algorithms with Guaranteed Error Bounds + for the Error Function and the Complementary Error Function + + http://www.math.uni-wuppertal.de/wrswt/literatur_en.html + + @author Kohei Yoshida <kohei@openoffice.org> + + @see #i55735# + */ +double SAL_CALL rtl_math_erf( double x ) SAL_THROW_EXTERN_C() +{ + if( x == 0.0 ) + return 0.0; + + bool bNegative = false; + if ( x < 0.0 ) + { + x = fabs( x ); + bNegative = true; + } + + double fErf = 1.0; + if ( x < 1.0e-10 ) + fErf = (double) (x*1.1283791670955125738961589031215452L); + else if ( x < 0.65 ) + lcl_Erf0065( x, fErf ); + else + fErf = 1.0 - rtl_math_erfc( x ); + + if ( bNegative ) + fErf *= -1.0; + + return fErf; +} + + +/** Parent complementary error function (erfc) that calls different algorithms + based on the value of x. It takes care of cases where x is negative as erfc + satisfies relationship erfc(-x) = 2 - erfc(x). See the comment for Erf(x) + for the source publication. + + @author Kohei Yoshida <kohei@openoffice.org> + + @see #i55735#, moved from module scaddins (#i97091#) + + */ +double SAL_CALL rtl_math_erfc( double x ) SAL_THROW_EXTERN_C() +{ + if ( x == 0.0 ) + return 1.0; + + bool bNegative = false; + if ( x < 0.0 ) + { + x = fabs( x ); + bNegative = true; + } + + double fErfc = 0.0; + if ( x >= 0.65 ) + { + if ( x < 6.0 ) + lcl_Erfc0600( x, fErfc ); + else + lcl_Erfc2654( x, fErfc ); + } + else + fErfc = 1.0 - rtl_math_erf( x ); + + if ( bNegative ) + fErfc = 2.0 - fErfc; + + return fErfc; +} + +/** improved accuracy of asinh for |x| large and for x near zero + @see #i97605# + */ +double SAL_CALL rtl_math_asinh( double fX ) SAL_THROW_EXTERN_C() +{ + double fSign = 1.0; + if ( fX == 0.0 ) + return 0.0; + else + { + if ( fX < 0.0 ) + { + fX = - fX; + fSign = -1.0; + } + if ( fX < 0.125 ) + return fSign * rtl_math_log1p( fX + fX*fX / (1.0 + sqrt( 1.0 + fX*fX))); + else if ( fX < 1.25e7 ) + return fSign * log( fX + sqrt( 1.0 + fX*fX)); + else + return fSign * log( 2.0*fX); + } +} + +/** improved accuracy of acosh for x large and for x near 1 + @see #i97605# + */ +double SAL_CALL rtl_math_acosh( double fX ) SAL_THROW_EXTERN_C() +{ + volatile double fZ = fX - 1.0; + if ( fX < 1.0 ) + { + double fResult; + ::rtl::math::setNan( &fResult ); + return fResult; + } + else if ( fX == 1.0 ) + return 0.0; + else if ( fX < 1.1 ) + return rtl_math_log1p( fZ + sqrt( fZ*fZ + 2.0*fZ)); + else if ( fX < 1.25e7 ) + return log( fX + sqrt( fX*fX - 1.0)); + else + return log( 2.0*fX); +} |