Skip to content

C
core
Kaydet (Commit) a62bc6a6 authored tarafından Eike Rathke's avatar Eike Rathke
Dosyalara gözat
Seçenekler

replace implementation of rtl_math_erf() and rtl_math_erfc() 

... with ::std::erf() and ::std::erfc() of C++11

Change-Id: I8ccc86ec4d6d71a92409770fc119f72e7084073a
Reviewed-on: https://gerrit.libreoffice.org/19583Tested-by: 's avatarJenkins <ci@libreoffice.org>
Reviewed-by: 's avatarEike Rathke <erack@redhat.com>
Tested-by: 's avatarEike Rathke <erack@redhat.com>
üst 14cd70ba
Hide whitespace changes
Inline Side-by-side
Showing with 4 additions and 225 deletions
sal/rtl/math.cxx
Dosyayı görüntüle @ a62bc6a6
......@@ -68,131 +68,6 @@ static double getN10Exp( int nExp )
return 1.0;
}
/** Approximation algorithm for erf for 0 < x < 0.65. */
static 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. */
static 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). */
static 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 {
double const nKorrVal[] = {
......@@ -1166,112 +1041,16 @@ 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#
*/
/** Parent error function (erf) */
double SAL_CALL rtl_math_erf( double x ) SAL_THROW_EXTERN_C()
{
if( x == 0.0 )
return 0.0;
// Otherwise we may end up in endless recursion through rtl_math_erfc().
if (!::rtl::math::isFinite(x))
{
// See http://en.cppreference.com/w/cpp/numeric/math/erf
if (::rtl::math::isInf(x))
{
if (::rtl::math::isSignBitSet(x))
return -1.0;
else
return 1.0;
}
// It is a NaN.
return x;
}
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;
return ::std::erf(x);
}
/** 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#)
*/
/** Parent complementary error function (erfc) */
double SAL_CALL rtl_math_erfc( double x ) SAL_THROW_EXTERN_C()
{
if ( x == 0.0 )
return 1.0;
// Otherwise we may end up in endless recursion through rtl_math_erf().
if (!::rtl::math::isFinite(x))
{
// See http://en.cppreference.com/w/cpp/numeric/math/erfc
if (::rtl::math::isInf(x))
{
if (::rtl::math::isSignBitSet(x))
return 2.0;
else
return 0.0;
}
// It is a NaN.
return x;
}
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;
return ::std::erfc(x);
}
/** improved accuracy of asinh for |x| large and for x near zero
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment