Print this page
11210 libm should be cstyle(1ONBLD) clean

*** 20,29 **** --- 20,30 ---- */ /* * Copyright 2011 Nexenta Systems, Inc. All rights reserved. */ + /* * Copyright 2006 Sun Microsystems, Inc. All rights reserved. * Use is subject to license terms. */
*** 97,113 **** #pragma weak __erfcl = erfcl #include "libm.h" #include "longdouble.h" ! static const long double ! tiny = 1e-40L, nearunfl = 1e-4000L, half = 0.5L, one = 1.0L, onehalf = 1.5L, ! L16_3 = 16.0L/3.0L; /* * Coefficients for even polynomial P for erf(x)=x+x*P(x^2) on [0,0.84375] */ static const long double P[] = { /* 21 coeffs */ 1.283791670955125738961589031215451715556e-0001L, --- 98,114 ---- #pragma weak __erfcl = erfcl #include "libm.h" #include "longdouble.h" ! static const long double tiny = 1e-40L, nearunfl = 1e-4000L, half = 0.5L, one = 1.0L, onehalf = 1.5L, ! L16_3 = 16.0L / 3.0L; ! /* * Coefficients for even polynomial P for erf(x)=x+x*P(x^2) on [0,0.84375] */ static const long double P[] = { /* 21 coeffs */ 1.283791670955125738961589031215451715556e-0001L,
*** 149,158 **** --- 150,160 ---- 5.390833481581033423020320734201065475098e-0004L, -1.978853912815115495053119023517805528300e-0004L, 6.184234513953600118335017885706420552487e-0005L, -5.331802711697810861017518515816271808286e-0006L, }; + static const long double Q1[] = { /* 12 bottom coeffs with leading 1.0 hidden */ 9.081506296064882195280178373107623196655e-0001L, 6.821049531968204097604392183650687642520e-0001L, 4.067869178233539502315055970743271822838e-0001L, 1.702332233546316765818144723063881095577e-0001L,
*** 163,172 **** --- 165,175 ---- 3.185620255011299476196039491205159718620e-0004L, 1.273405072153008775426376193374105840517e-0005L, 4.753866999959432971956781228148402971454e-0006L, -1.002287602111660026053981728549540200683e-0006L, }; + /* * Rational erf(x) = ((float)0.95478588343) + P2(x-1.5)/Q2(x-1.5) * on [1.25,1.75] */ static const long double C2 = (long double)((float)0.95478588343);
*** 182,191 **** --- 185,195 ---- -4.289851942513144714600285769022420962418e-0005L, 8.304719841341952705874781636002085119978e-0005L, -1.040460226177309338781902252282849903189e-0005L, 2.122913331584921470381327583672044434087e-0006L, }; + static const long double Q2[] = { /* 13 bottom coeffs with leading 1.0 hidden */ 7.448815737306992749168727691042003832150e-0001L, 7.161813850236008294484744312430122188043e-0001L, 3.603134756584225766144922727405641236121e-0001L, 1.955811609133766478080550795194535852653e-0001L,
*** 197,206 **** --- 201,211 ---- 8.664587895570043348530991997272212150316e-0005L, 1.109201582511752087060167429397033701988e-0005L, 1.357834375781831062713347000030984364311e-0006L, 4.957746280594384997273090385060680016451e-0008L, }; + /* * erfc(x) = exp(-x*x)/x * R1(1/x)/S1(1/x) on [1.75, 16/3] */ static const long double R1[] = { /* 14 top coeffs */ 4.630195122654315016370705767621550602948e+0006L,
*** 216,225 **** --- 221,231 ---- 2.839793161868140305907004392890348777338e+0003L, 2.786687241658423601778258694498655680778e+0002L, 1.779177837102695602425897452623985786464e+0001L, 5.641895835477470769043614623819144434731e-0001L, }; + static const long double S1[] = { /* 15 bottom coeffs with leading 1.0 hidden */ 4.630195122654331529595606896287596843110e+0006L, 1.780411093345512024324781084220509055058e+0007L, 3.250113097051800703707108623715776848283e+0007L, 3.737857099176755050912193712123489115755e+0007L,
*** 237,247 **** }; /* * erfc(x) = exp(-x*x)/x * R2(1/x)/S2(1/x) on [16/3, 107] */ ! static const long double R2[] = { /* 15 top coeffs in reverse order!!*/ 2.447288012254302966796326587537136931669e+0005L, 8.768592567189861896653369912716538739016e+0005L, 1.552293152581780065761497908005779524953e+0006L, 1.792075924835942935864231657504259926729e+0006L, 1.504001463155897344947500222052694835875e+0006L, --- 243,253 ---- }; /* * erfc(x) = exp(-x*x)/x * R2(1/x)/S2(1/x) on [16/3, 107] */ ! static const long double R2[] = { /* 15 top coeffs in reverse order!! */ 2.447288012254302966796326587537136931669e+0005L, 8.768592567189861896653369912716538739016e+0005L, 1.552293152581780065761497908005779524953e+0006L, 1.792075924835942935864231657504259926729e+0006L, 1.504001463155897344947500222052694835875e+0006L,
*** 254,263 **** --- 260,270 ---- 7.362346487427048068212968889642741734621e+0002L, 9.980359714211411423007641056580813116207e+0001L, 9.426910895135379181107191962193485174159e+0000L, 5.641895835477562869480794515623601280429e-0001L, }; + static const long double S2[] = { /* 16 coefficients */ 2.447282203601902971246004716790604686880e+0005L, 1.153009852759385309367759460934808489833e+0006L, 2.608580649612639131548966265078663384849e+0006L, 3.766673917346623308850202792390569025740e+0006L,
*** 273,347 **** 1.773972700887629157006326333696896516769e+0002L, 1.670876451822586800422009013880457094162e+0001L, 1.000L, }; ! long double erfl(x) ! long double x; { ! long double s,y,t; if (!finitel(x)) { ! if (x != x) return x+x; /* NaN */ ! return copysignl(one,x); /* return +-1.0 is x=Inf */ } y = fabsl(x); if (y <= 0.84375L) { ! if (y<=tiny) return x+P[0]*x; ! s = y*y; ! t = __poly_libmq(s,21,P); ! return x+x*t; } ! if (y<=1.25L) { ! s = y-one; ! t = C1+__poly_libmq(s,12,P1)/(one+s*__poly_libmq(s,12,Q1)); ! return (signbitl(x))? -t: t; ! } else if (y<=1.75L) { ! s = y-onehalf; ! t = C2+__poly_libmq(s,12,P2)/(one+s*__poly_libmq(s,13,Q2)); ! return (signbitl(x))? -t: t; } ! if (y<=9.0L) t = erfcl(y); else t = tiny; ! return (signbitl(x))? t-one: one-t; } ! long double erfcl(x) ! long double x; { ! long double s,y,t; if (!finitel(x)) { ! if (x != x) return x+x; /* NaN */ /* return 2.0 if x= -inf; 0.0 if x= +inf */ ! if (x < 0.0L) return 2.0L; else return 0.0L; } if (x <= 0.84375L) { ! if (x<=0.25) return one-erfl(x); ! s = x*x; ! t = half-x; ! t = t - x*__poly_libmq(s,21,P); ! return half+t; } ! if (x<=1.25L) { ! s = x-one; ! t = one-C1; ! return t - __poly_libmq(s,12,P1)/(one+s*__poly_libmq(s,12,Q1)); ! } else if (x<=1.75L) { ! s = x-onehalf; ! t = one-C2; ! return t - __poly_libmq(s,12,P2)/(one+s*__poly_libmq(s,13,Q2)); } ! if (x>=107.0L) return nearunfl*nearunfl; /* underflow */ ! else if (x >= L16_3) { ! y = __poly_libmq(x,15,R2); ! t = y/__poly_libmq(x,16,S2); } else { ! y = __poly_libmq(x,14,R1); ! t = y/__poly_libmq(x,15,S1); } /* * Note that exp(-x*x+d) = exp(-x*x)*exp(d), so to compute * exp(-x*x) with a small relative error, we need to compute * -x*x with a small absolute error. To this end, we set y * equal to the leading part of x but with enough trailing --- 280,380 ---- 1.773972700887629157006326333696896516769e+0002L, 1.670876451822586800422009013880457094162e+0001L, 1.000L, }; ! long double ! erfl(long double x) { ! long double s, y, t; if (!finitel(x)) { ! if (x != x) ! return (x + x); /* NaN */ ! ! return (copysignl(one, x)); /* return +-1.0 is x=Inf */ } y = fabsl(x); + if (y <= 0.84375L) { ! if (y <= tiny) ! return (x + P[0] * x); ! ! s = y * y; ! t = __poly_libmq(s, 21, P); ! return (x + x * t); } ! ! if (y <= 1.25L) { ! s = y - one; ! t = C1 + __poly_libmq(s, 12, P1) / (one + s * __poly_libmq(s, ! 12, Q1)); ! return ((signbitl(x)) ? -t : t); ! } else if (y <= 1.75L) { ! s = y - onehalf; ! t = C2 + __poly_libmq(s, 12, P2) / (one + s * __poly_libmq(s, ! 13, Q2)); ! return ((signbitl(x)) ? -t : t); } ! ! if (y <= 9.0L) ! t = erfcl(y); ! else ! t = tiny; ! ! return ((signbitl(x)) ? t - one : one - t); } ! long double ! erfcl(long double x) { ! long double s, y, t; if (!finitel(x)) { ! if (x != x) ! return (x + x); /* NaN */ ! /* return 2.0 if x= -inf; 0.0 if x= +inf */ ! if (x < 0.0L) ! return (2.0L); ! else ! return (0.0L); } if (x <= 0.84375L) { ! if (x <= 0.25) ! return (one - erfl(x)); ! ! s = x * x; ! t = half - x; ! t = t - x * __poly_libmq(s, 21, P); ! return (half + t); } ! ! if (x <= 1.25L) { ! s = x - one; ! t = one - C1; ! return (t - __poly_libmq(s, 12, P1) / (one + s * __poly_libmq(s, ! 12, Q1))); ! } else if (x <= 1.75L) { ! s = x - onehalf; ! t = one - C2; ! return (t - __poly_libmq(s, 12, P2) / (one + s * __poly_libmq(s, ! 13, Q2))); } ! ! if (x >= 107.0L) { ! return (nearunfl * nearunfl); /* underflow */ ! } else if (x >= L16_3) { ! y = __poly_libmq(x, 15, R2); ! t = y / __poly_libmq(x, 16, S2); } else { ! y = __poly_libmq(x, 14, R1); ! t = y / __poly_libmq(x, 15, S1); } + /* * Note that exp(-x*x+d) = exp(-x*x)*exp(d), so to compute * exp(-x*x) with a small relative error, we need to compute * -x*x with a small absolute error. To this end, we set y * equal to the leading part of x but with enough trailing
*** 358,366 **** * small but not so large that the conversion to int overflows. * When long double arithmetic is slow, however, the following * non-portable code is preferable. */ y = x; ! *(2+(int*)&y) = *(3+(int*)&y) = 0; ! t *= expl(-y*y)*expl(-(x-y)*(x+y)); ! return t; } --- 391,399 ---- * small but not so large that the conversion to int overflows. * When long double arithmetic is slow, however, the following * non-portable code is preferable. */ y = x; ! *(2 + (int *)&y) = *(3 + (int *)&y) = 0; ! t *= expl(-y * y) * expl(-(x - y) * (x + y)); ! return (t); }