Back to index

glibc  2.9
s_erfl.c
Go to the documentation of this file.
00001 /*
00002  * ====================================================
00003  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
00004  *
00005  * Developed at SunPro, a Sun Microsystems, Inc. business.
00006  * Permission to use, copy, modify, and distribute this
00007  * software is freely granted, provided that this notice
00008  * is preserved.
00009  * ====================================================
00010  */
00011 
00012 /* Modifications and expansions for 128-bit long double are
00013    Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov>
00014    and are incorporated herein by permission of the author.  The author
00015    reserves the right to distribute this material elsewhere under different
00016    copying permissions.  These modifications are distributed here under
00017    the following terms:
00018 
00019     This library is free software; you can redistribute it and/or
00020     modify it under the terms of the GNU Lesser General Public
00021     License as published by the Free Software Foundation; either
00022     version 2.1 of the License, or (at your option) any later version.
00023 
00024     This library is distributed in the hope that it will be useful,
00025     but WITHOUT ANY WARRANTY; without even the implied warranty of
00026     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00027     Lesser General Public License for more details.
00028 
00029     You should have received a copy of the GNU Lesser General Public
00030     License along with this library; if not, write to the Free Software
00031     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA */
00032 
00033 /* double erf(double x)
00034  * double erfc(double x)
00035  *                        x
00036  *                  2      |\
00037  *     erf(x)  =  ---------  | exp(-t*t)dt
00038  *               sqrt(pi) \|
00039  *                        0
00040  *
00041  *     erfc(x) =  1-erf(x)
00042  *  Note that
00043  *            erf(-x) = -erf(x)
00044  *            erfc(-x) = 2 - erfc(x)
00045  *
00046  * Method:
00047  *     1.  erf(x)  = x + x*R(x^2) for |x| in [0, 7/8]
00048  *        Remark. The formula is derived by noting
00049  *          erf(x) = (2/sqrt(pi))*(x - x^3/3 + x^5/10 - x^7/42 + ....)
00050  *        and that
00051  *          2/sqrt(pi) = 1.128379167095512573896158903121545171688
00052  *        is close to one.
00053  *
00054  *      1a. erf(x)  = 1 - erfc(x), for |x| > 1.0
00055  *          erfc(x) = 1 - erf(x)  if |x| < 1/4
00056  *
00057  *      2. For |x| in [7/8, 1], let s = |x| - 1, and
00058  *         c = 0.84506291151 rounded to single (24 bits)
00059  *     erf(s + c)  = sign(x) * (c  + P1(s)/Q1(s))
00060  *        Remark: here we use the taylor series expansion at x=1.
00061  *            erf(1+s) = erf(1) + s*Poly(s)
00062  *                    = 0.845.. + P1(s)/Q1(s)
00063  *        Note that |P1/Q1|< 0.078 for x in [0.84375,1.25]
00064  *
00065  *      3. For x in [1/4, 5/4],
00066  *     erfc(s + const)  = erfc(const)  + s P1(s)/Q1(s)
00067  *              for const = 1/4, 3/8, ..., 9/8
00068  *              and 0 <= s <= 1/8 .
00069  *
00070  *      4. For x in [5/4, 107],
00071  *     erfc(x) = (1/x)*exp(-x*x-0.5625 + R(z))
00072  *              z=1/x^2
00073  *         The interval is partitioned into several segments
00074  *         of width 1/8 in 1/x.
00075  *
00076  *      Note1:
00077  *        To compute exp(-x*x-0.5625+R/S), let s be a single
00078  *        precision number and s := x; then
00079  *            -x*x = -s*s + (s-x)*(s+x)
00080  *             exp(-x*x-0.5626+R/S) =
00081  *                   exp(-s*s-0.5625)*exp((s-x)*(s+x)+R/S);
00082  *      Note2:
00083  *        Here 4 and 5 make use of the asymptotic series
00084  *                     exp(-x*x)
00085  *            erfc(x) ~ ---------- * ( 1 + Poly(1/x^2) )
00086  *                     x*sqrt(pi)
00087  *
00088  *      5. For inf > x >= 107
00089  *     erf(x)  = sign(x) *(1 - tiny)  (raise inexact)
00090  *     erfc(x) = tiny*tiny (raise underflow) if x > 0
00091  *                   = 2 - tiny if x<0
00092  *
00093  *      7. Special case:
00094  *     erf(0)  = 0, erf(inf)  = 1, erf(-inf) = -1,
00095  *     erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2,
00096  *            erfc/erf(NaN) is NaN
00097  */
00098 
00099 #include "math.h"
00100 #include "math_private.h"
00101 #include <math_ldbl_opt.h>
00102 
00103 /* Evaluate P[n] x^n  +  P[n-1] x^(n-1)  +  ...  +  P[0] */
00104 
00105 static long double
00106 neval (long double x, const long double *p, int n)
00107 {
00108   long double y;
00109 
00110   p += n;
00111   y = *p--;
00112   do
00113     {
00114       y = y * x + *p--;
00115     }
00116   while (--n > 0);
00117   return y;
00118 }
00119 
00120 
00121 /* Evaluate x^n+1  +  P[n] x^(n)  +  P[n-1] x^(n-1)  +  ...  +  P[0] */
00122 
00123 static long double
00124 deval (long double x, const long double *p, int n)
00125 {
00126   long double y;
00127 
00128   p += n;
00129   y = x + *p--;
00130   do
00131     {
00132       y = y * x + *p--;
00133     }
00134   while (--n > 0);
00135   return y;
00136 }
00137 
00138 
00139 
00140 #ifdef __STDC__
00141 static const long double
00142 #else
00143 static long double
00144 #endif
00145 tiny = 1e-300L,
00146   half = 0.5L,
00147   one = 1.0L,
00148   two = 2.0L,
00149   /* 2/sqrt(pi) - 1 */
00150   efx = 1.2837916709551257389615890312154517168810E-1L,
00151   /* 8 * (2/sqrt(pi) - 1) */
00152   efx8 = 1.0270333367641005911692712249723613735048E0L;
00153 
00154 
00155 /* erf(x)  = x  + x R(x^2)
00156    0 <= x <= 7/8
00157    Peak relative error 1.8e-35  */
00158 #define NTN1 8
00159 static const long double TN1[NTN1 + 1] =
00160 {
00161  -3.858252324254637124543172907442106422373E10L,
00162   9.580319248590464682316366876952214879858E10L,
00163   1.302170519734879977595901236693040544854E10L,
00164   2.922956950426397417800321486727032845006E9L,
00165   1.764317520783319397868923218385468729799E8L,
00166   1.573436014601118630105796794840834145120E7L,
00167   4.028077380105721388745632295157816229289E5L,
00168   1.644056806467289066852135096352853491530E4L,
00169   3.390868480059991640235675479463287886081E1L
00170 };
00171 #define NTD1 8
00172 static const long double TD1[NTD1 + 1] =
00173 {
00174   -3.005357030696532927149885530689529032152E11L,
00175   -1.342602283126282827411658673839982164042E11L,
00176   -2.777153893355340961288511024443668743399E10L,
00177   -3.483826391033531996955620074072768276974E9L,
00178   -2.906321047071299585682722511260895227921E8L,
00179   -1.653347985722154162439387878512427542691E7L,
00180   -6.245520581562848778466500301865173123136E5L,
00181   -1.402124304177498828590239373389110545142E4L,
00182   -1.209368072473510674493129989468348633579E2L
00183 /* 1.0E0 */
00184 };
00185 
00186 
00187 /* erf(z+1)  = erf_const + P(z)/Q(z)
00188    -.125 <= z <= 0
00189    Peak relative error 7.3e-36  */
00190 static const long double erf_const = 0.845062911510467529296875L;
00191 #define NTN2 8
00192 static const long double TN2[NTN2 + 1] =
00193 {
00194  -4.088889697077485301010486931817357000235E1L,
00195   7.157046430681808553842307502826960051036E3L,
00196  -2.191561912574409865550015485451373731780E3L,
00197   2.180174916555316874988981177654057337219E3L,
00198   2.848578658049670668231333682379720943455E2L,
00199   1.630362490952512836762810462174798925274E2L,
00200   6.317712353961866974143739396865293596895E0L,
00201   2.450441034183492434655586496522857578066E1L,
00202   5.127662277706787664956025545897050896203E-1L
00203 };
00204 #define NTD2 8
00205 static const long double TD2[NTD2 + 1] =
00206 {
00207   1.731026445926834008273768924015161048885E4L,
00208   1.209682239007990370796112604286048173750E4L,
00209   1.160950290217993641320602282462976163857E4L,
00210   5.394294645127126577825507169061355698157E3L,
00211   2.791239340533632669442158497532521776093E3L,
00212   8.989365571337319032943005387378993827684E2L,
00213   2.974016493766349409725385710897298069677E2L,
00214   6.148192754590376378740261072533527271947E1L,
00215   1.178502892490738445655468927408440847480E1L
00216  /* 1.0E0 */
00217 };
00218 
00219 
00220 /* erfc(x + 0.25) = erfc(0.25) + x R(x)
00221    0 <= x < 0.125
00222    Peak relative error 1.4e-35  */
00223 #define NRNr13 8
00224 static const long double RNr13[NRNr13 + 1] =
00225 {
00226  -2.353707097641280550282633036456457014829E3L,
00227   3.871159656228743599994116143079870279866E2L,
00228  -3.888105134258266192210485617504098426679E2L,
00229  -2.129998539120061668038806696199343094971E1L,
00230  -8.125462263594034672468446317145384108734E1L,
00231   8.151549093983505810118308635926270319660E0L,
00232  -5.033362032729207310462422357772568553670E0L,
00233  -4.253956621135136090295893547735851168471E-2L,
00234  -8.098602878463854789780108161581050357814E-2L
00235 };
00236 #define NRDr13 7
00237 static const long double RDr13[NRDr13 + 1] =
00238 {
00239   2.220448796306693503549505450626652881752E3L,
00240   1.899133258779578688791041599040951431383E2L,
00241   1.061906712284961110196427571557149268454E3L,
00242   7.497086072306967965180978101974566760042E1L,
00243   2.146796115662672795876463568170441327274E2L,
00244   1.120156008362573736664338015952284925592E1L,
00245   2.211014952075052616409845051695042741074E1L,
00246   6.469655675326150785692908453094054988938E-1L
00247  /* 1.0E0 */
00248 };
00249 /* erfc(0.25) = C13a + C13b to extra precision.  */
00250 static const long double C13a = 0.723663330078125L;
00251 static const long double C13b = 1.0279753638067014931732235184287934646022E-5L;
00252 
00253 
00254 /* erfc(x + 0.375) = erfc(0.375) + x R(x)
00255    0 <= x < 0.125
00256    Peak relative error 1.2e-35  */
00257 #define NRNr14 8
00258 static const long double RNr14[NRNr14 + 1] =
00259 {
00260  -2.446164016404426277577283038988918202456E3L,
00261   6.718753324496563913392217011618096698140E2L,
00262  -4.581631138049836157425391886957389240794E2L,
00263  -2.382844088987092233033215402335026078208E1L,
00264  -7.119237852400600507927038680970936336458E1L,
00265   1.313609646108420136332418282286454287146E1L,
00266  -6.188608702082264389155862490056401365834E0L,
00267  -2.787116601106678287277373011101132659279E-2L,
00268  -2.230395570574153963203348263549700967918E-2L
00269 };
00270 #define NRDr14 7
00271 static const long double RDr14[NRDr14 + 1] =
00272 {
00273   2.495187439241869732696223349840963702875E3L,
00274   2.503549449872925580011284635695738412162E2L,
00275   1.159033560988895481698051531263861842461E3L,
00276   9.493751466542304491261487998684383688622E1L,
00277   2.276214929562354328261422263078480321204E2L,
00278   1.367697521219069280358984081407807931847E1L,
00279   2.276988395995528495055594829206582732682E1L,
00280   7.647745753648996559837591812375456641163E-1L
00281  /* 1.0E0 */
00282 };
00283 /* erfc(0.375) = C14a + C14b to extra precision.  */
00284 static const long double C14a = 0.5958709716796875L;
00285 static const long double C14b = 1.2118885490201676174914080878232469565953E-5L;
00286 
00287 /* erfc(x + 0.5) = erfc(0.5) + x R(x)
00288    0 <= x < 0.125
00289    Peak relative error 4.7e-36  */
00290 #define NRNr15 8
00291 static const long double RNr15[NRNr15 + 1] =
00292 {
00293  -2.624212418011181487924855581955853461925E3L,
00294   8.473828904647825181073831556439301342756E2L,
00295  -5.286207458628380765099405359607331669027E2L,
00296  -3.895781234155315729088407259045269652318E1L,
00297  -6.200857908065163618041240848728398496256E1L,
00298   1.469324610346924001393137895116129204737E1L,
00299  -6.961356525370658572800674953305625578903E0L,
00300   5.145724386641163809595512876629030548495E-3L,
00301   1.990253655948179713415957791776180406812E-2L
00302 };
00303 #define NRDr15 7
00304 static const long double RDr15[NRDr15 + 1] =
00305 {
00306   2.986190760847974943034021764693341524962E3L,
00307   5.288262758961073066335410218650047725985E2L,
00308   1.363649178071006978355113026427856008978E3L,
00309   1.921707975649915894241864988942255320833E2L,
00310   2.588651100651029023069013885900085533226E2L,
00311   2.628752920321455606558942309396855629459E1L,
00312   2.455649035885114308978333741080991380610E1L,
00313   1.378826653595128464383127836412100939126E0L
00314   /* 1.0E0 */
00315 };
00316 /* erfc(0.5) = C15a + C15b to extra precision.  */
00317 static const long double C15a = 0.4794921875L;
00318 static const long double C15b = 7.9346869534623172533461080354712635484242E-6L;
00319 
00320 /* erfc(x + 0.625) = erfc(0.625) + x R(x)
00321    0 <= x < 0.125
00322    Peak relative error 5.1e-36  */
00323 #define NRNr16 8
00324 static const long double RNr16[NRNr16 + 1] =
00325 {
00326  -2.347887943200680563784690094002722906820E3L,
00327   8.008590660692105004780722726421020136482E2L,
00328  -5.257363310384119728760181252132311447963E2L,
00329  -4.471737717857801230450290232600243795637E1L,
00330  -4.849540386452573306708795324759300320304E1L,
00331   1.140885264677134679275986782978655952843E1L,
00332  -6.731591085460269447926746876983786152300E0L,
00333   1.370831653033047440345050025876085121231E-1L,
00334   2.022958279982138755020825717073966576670E-2L,
00335 };
00336 #define NRDr16 7
00337 static const long double RDr16[NRDr16 + 1] =
00338 {
00339   3.075166170024837215399323264868308087281E3L,
00340   8.730468942160798031608053127270430036627E2L,
00341   1.458472799166340479742581949088453244767E3L,
00342   3.230423687568019709453130785873540386217E2L,
00343   2.804009872719893612081109617983169474655E2L,
00344   4.465334221323222943418085830026979293091E1L,
00345   2.612723259683205928103787842214809134746E1L,
00346   2.341526751185244109722204018543276124997E0L,
00347   /* 1.0E0 */
00348 };
00349 /* erfc(0.625) = C16a + C16b to extra precision.  */
00350 static const long double C16a = 0.3767547607421875L;
00351 static const long double C16b = 4.3570693945275513594941232097252997287766E-6L;
00352 
00353 /* erfc(x + 0.75) = erfc(0.75) + x R(x)
00354    0 <= x < 0.125
00355    Peak relative error 1.7e-35  */
00356 #define NRNr17 8
00357 static const long double RNr17[NRNr17 + 1] =
00358 {
00359   -1.767068734220277728233364375724380366826E3L,
00360   6.693746645665242832426891888805363898707E2L,
00361   -4.746224241837275958126060307406616817753E2L,
00362   -2.274160637728782675145666064841883803196E1L,
00363   -3.541232266140939050094370552538987982637E1L,
00364   6.988950514747052676394491563585179503865E0L,
00365   -5.807687216836540830881352383529281215100E0L,
00366   3.631915988567346438830283503729569443642E-1L,
00367   -1.488945487149634820537348176770282391202E-2L
00368 };
00369 #define NRDr17 7
00370 static const long double RDr17[NRDr17 + 1] =
00371 {
00372   2.748457523498150741964464942246913394647E3L,
00373   1.020213390713477686776037331757871252652E3L,
00374   1.388857635935432621972601695296561952738E3L,
00375   3.903363681143817750895999579637315491087E2L,
00376   2.784568344378139499217928969529219886578E2L,
00377   5.555800830216764702779238020065345401144E1L,
00378   2.646215470959050279430447295801291168941E1L,
00379   2.984905282103517497081766758550112011265E0L,
00380   /* 1.0E0 */
00381 };
00382 /* erfc(0.75) = C17a + C17b to extra precision.  */
00383 static const long double C17a = 0.2888336181640625L;
00384 static const long double C17b = 1.0748182422368401062165408589222625794046E-5L;
00385 
00386 
00387 /* erfc(x + 0.875) = erfc(0.875) + x R(x)
00388    0 <= x < 0.125
00389    Peak relative error 2.2e-35  */
00390 #define NRNr18 8
00391 static const long double RNr18[NRNr18 + 1] =
00392 {
00393  -1.342044899087593397419622771847219619588E3L,
00394   6.127221294229172997509252330961641850598E2L,
00395  -4.519821356522291185621206350470820610727E2L,
00396   1.223275177825128732497510264197915160235E1L,
00397  -2.730789571382971355625020710543532867692E1L,
00398   4.045181204921538886880171727755445395862E0L,
00399  -4.925146477876592723401384464691452700539E0L,
00400   5.933878036611279244654299924101068088582E-1L,
00401  -5.557645435858916025452563379795159124753E-2L
00402 };
00403 #define NRDr18 7
00404 static const long double RDr18[NRDr18 + 1] =
00405 {
00406   2.557518000661700588758505116291983092951E3L,
00407   1.070171433382888994954602511991940418588E3L,
00408   1.344842834423493081054489613250688918709E3L,
00409   4.161144478449381901208660598266288188426E2L,
00410   2.763670252219855198052378138756906980422E2L,
00411   5.998153487868943708236273854747564557632E1L,
00412   2.657695108438628847733050476209037025318E1L,
00413   3.252140524394421868923289114410336976512E0L,
00414   /* 1.0E0 */
00415 };
00416 /* erfc(0.875) = C18a + C18b to extra precision.  */
00417 static const long double C18a = 0.215911865234375L;
00418 static const long double C18b = 1.3073705765341685464282101150637224028267E-5L;
00419 
00420 /* erfc(x + 1.0) = erfc(1.0) + x R(x)
00421    0 <= x < 0.125
00422    Peak relative error 1.6e-35  */
00423 #define NRNr19 8
00424 static const long double RNr19[NRNr19 + 1] =
00425 {
00426  -1.139180936454157193495882956565663294826E3L,
00427   6.134903129086899737514712477207945973616E2L,
00428  -4.628909024715329562325555164720732868263E2L,
00429   4.165702387210732352564932347500364010833E1L,
00430  -2.286979913515229747204101330405771801610E1L,
00431   1.870695256449872743066783202326943667722E0L,
00432  -4.177486601273105752879868187237000032364E0L,
00433   7.533980372789646140112424811291782526263E-1L,
00434  -8.629945436917752003058064731308767664446E-2L
00435 };
00436 #define NRDr19 7
00437 static const long double RDr19[NRDr19 + 1] =
00438 {
00439   2.744303447981132701432716278363418643778E3L,
00440   1.266396359526187065222528050591302171471E3L,
00441   1.466739461422073351497972255511919814273E3L,
00442   4.868710570759693955597496520298058147162E2L,
00443   2.993694301559756046478189634131722579643E2L,
00444   6.868976819510254139741559102693828237440E1L,
00445   2.801505816247677193480190483913753613630E1L,
00446   3.604439909194350263552750347742663954481E0L,
00447   /* 1.0E0 */
00448 };
00449 /* erfc(1.0) = C19a + C19b to extra precision.  */
00450 static const long double C19a = 0.15728759765625L;
00451 static const long double C19b = 1.1609394035130658779364917390740703933002E-5L;
00452 
00453 /* erfc(x + 1.125) = erfc(1.125) + x R(x)
00454    0 <= x < 0.125
00455    Peak relative error 3.6e-36  */
00456 #define NRNr20 8
00457 static const long double RNr20[NRNr20 + 1] =
00458 {
00459  -9.652706916457973956366721379612508047640E2L,
00460   5.577066396050932776683469951773643880634E2L,
00461  -4.406335508848496713572223098693575485978E2L,
00462   5.202893466490242733570232680736966655434E1L,
00463  -1.931311847665757913322495948705563937159E1L,
00464  -9.364318268748287664267341457164918090611E-2L,
00465  -3.306390351286352764891355375882586201069E0L,
00466   7.573806045289044647727613003096916516475E-1L,
00467  -9.611744011489092894027478899545635991213E-2L
00468 };
00469 #define NRDr20 7
00470 static const long double RDr20[NRDr20 + 1] =
00471 {
00472   3.032829629520142564106649167182428189014E3L,
00473   1.659648470721967719961167083684972196891E3L,
00474   1.703545128657284619402511356932569292535E3L,
00475   6.393465677731598872500200253155257708763E2L,
00476   3.489131397281030947405287112726059221934E2L,
00477   8.848641738570783406484348434387611713070E1L,
00478   3.132269062552392974833215844236160958502E1L,
00479   4.430131663290563523933419966185230513168E0L
00480  /* 1.0E0 */
00481 };
00482 /* erfc(1.125) = C20a + C20b to extra precision.  */
00483 static const long double C20a = 0.111602783203125L;
00484 static const long double C20b = 8.9850951672359304215530728365232161564636E-6L;
00485 
00486 /* erfc(1/x) = 1/x exp (-1/x^2 - 0.5625 + R(1/x^2))
00487    7/8 <= 1/x < 1
00488    Peak relative error 1.4e-35  */
00489 #define NRNr8 9
00490 static const long double RNr8[NRNr8 + 1] =
00491 {
00492   3.587451489255356250759834295199296936784E1L,
00493   5.406249749087340431871378009874875889602E2L,
00494   2.931301290625250886238822286506381194157E3L,
00495   7.359254185241795584113047248898753470923E3L,
00496   9.201031849810636104112101947312492532314E3L,
00497   5.749697096193191467751650366613289284777E3L,
00498   1.710415234419860825710780802678697889231E3L,
00499   2.150753982543378580859546706243022719599E2L,
00500   8.740953582272147335100537849981160931197E0L,
00501   4.876422978828717219629814794707963640913E-2L
00502 };
00503 #define NRDr8 8
00504 static const long double RDr8[NRDr8 + 1] =
00505 {
00506   6.358593134096908350929496535931630140282E1L,
00507   9.900253816552450073757174323424051765523E2L,
00508   5.642928777856801020545245437089490805186E3L,
00509   1.524195375199570868195152698617273739609E4L,
00510   2.113829644500006749947332935305800887345E4L,
00511   1.526438562626465706267943737310282977138E4L,
00512   5.561370922149241457131421914140039411782E3L,
00513   9.394035530179705051609070428036834496942E2L,
00514   6.147019596150394577984175188032707343615E1L
00515   /* 1.0E0 */
00516 };
00517 
00518 /* erfc(1/x) = 1/x exp (-1/x^2 - 0.5625 + R(1/x^2))
00519    0.75 <= 1/x <= 0.875
00520    Peak relative error 2.0e-36  */
00521 #define NRNr7 9
00522 static const long double RNr7[NRNr7 + 1] =
00523 {
00524  1.686222193385987690785945787708644476545E1L,
00525  1.178224543567604215602418571310612066594E3L,
00526  1.764550584290149466653899886088166091093E4L,
00527  1.073758321890334822002849369898232811561E5L,
00528  3.132840749205943137619839114451290324371E5L,
00529  4.607864939974100224615527007793867585915E5L,
00530  3.389781820105852303125270837910972384510E5L,
00531  1.174042187110565202875011358512564753399E5L,
00532  1.660013606011167144046604892622504338313E4L,
00533  6.700393957480661937695573729183733234400E2L
00534 };
00535 #define NRDr7 9
00536 static const long double RDr7[NRDr7 + 1] =
00537 {
00538 -1.709305024718358874701575813642933561169E3L,
00539 -3.280033887481333199580464617020514788369E4L,
00540 -2.345284228022521885093072363418750835214E5L,
00541 -8.086758123097763971926711729242327554917E5L,
00542 -1.456900414510108718402423999575992450138E6L,
00543 -1.391654264881255068392389037292702041855E6L,
00544 -6.842360801869939983674527468509852583855E5L,
00545 -1.597430214446573566179675395199807533371E5L,
00546 -1.488876130609876681421645314851760773480E4L,
00547 -3.511762950935060301403599443436465645703E2L
00548  /* 1.0E0 */
00549 };
00550 
00551 /* erfc(1/x) = 1/x exp(-1/x^2 - 0.5625 + R(1/x^2))
00552    5/8 <= 1/x < 3/4
00553    Peak relative error 1.9e-35  */
00554 #define NRNr6 9
00555 static const long double RNr6[NRNr6 + 1] =
00556 {
00557  1.642076876176834390623842732352935761108E0L,
00558  1.207150003611117689000664385596211076662E2L,
00559  2.119260779316389904742873816462800103939E3L,
00560  1.562942227734663441801452930916044224174E4L,
00561  5.656779189549710079988084081145693580479E4L,
00562  1.052166241021481691922831746350942786299E5L,
00563  9.949798524786000595621602790068349165758E4L,
00564  4.491790734080265043407035220188849562856E4L,
00565  8.377074098301530326270432059434791287601E3L,
00566  4.506934806567986810091824791963991057083E2L
00567 };
00568 #define NRDr6 9
00569 static const long double RDr6[NRDr6 + 1] =
00570 {
00571 -1.664557643928263091879301304019826629067E2L,
00572 -3.800035902507656624590531122291160668452E3L,
00573 -3.277028191591734928360050685359277076056E4L,
00574 -1.381359471502885446400589109566587443987E5L,
00575 -3.082204287382581873532528989283748656546E5L,
00576 -3.691071488256738343008271448234631037095E5L,
00577 -2.300482443038349815750714219117566715043E5L,
00578 -6.873955300927636236692803579555752171530E4L,
00579 -8.262158817978334142081581542749986845399E3L,
00580 -2.517122254384430859629423488157361983661E2L
00581  /* 1.00 */
00582 };
00583 
00584 /* erfc(1/x) = 1/x exp(-1/x^2 - 0.5625 + R(1/x^2))
00585    1/2 <= 1/x < 5/8
00586    Peak relative error 4.6e-36  */
00587 #define NRNr5 10
00588 static const long double RNr5[NRNr5 + 1] =
00589 {
00590 -3.332258927455285458355550878136506961608E-3L,
00591 -2.697100758900280402659586595884478660721E-1L,
00592 -6.083328551139621521416618424949137195536E0L,
00593 -6.119863528983308012970821226810162441263E1L,
00594 -3.176535282475593173248810678636522589861E2L,
00595 -8.933395175080560925809992467187963260693E2L,
00596 -1.360019508488475978060917477620199499560E3L,
00597 -1.075075579828188621541398761300910213280E3L,
00598 -4.017346561586014822824459436695197089916E2L,
00599 -5.857581368145266249509589726077645791341E1L,
00600 -2.077715925587834606379119585995758954399E0L
00601 };
00602 #define NRDr5 9
00603 static const long double RDr5[NRDr5 + 1] =
00604 {
00605  3.377879570417399341550710467744693125385E-1L,
00606  1.021963322742390735430008860602594456187E1L,
00607  1.200847646592942095192766255154827011939E2L,
00608  7.118915528142927104078182863387116942836E2L,
00609  2.318159380062066469386544552429625026238E3L,
00610  4.238729853534009221025582008928765281620E3L,
00611  4.279114907284825886266493994833515580782E3L,
00612  2.257277186663261531053293222591851737504E3L,
00613  5.570475501285054293371908382916063822957E2L,
00614  5.142189243856288981145786492585432443560E1L
00615  /* 1.0E0 */
00616 };
00617 
00618 /* erfc(1/x) = 1/x exp(-1/x^2 - 0.5625 + R(1/x^2))
00619    3/8 <= 1/x < 1/2
00620    Peak relative error 2.0e-36  */
00621 #define NRNr4 10
00622 static const long double RNr4[NRNr4 + 1] =
00623 {
00624  3.258530712024527835089319075288494524465E-3L,
00625  2.987056016877277929720231688689431056567E-1L,
00626  8.738729089340199750734409156830371528862E0L,
00627  1.207211160148647782396337792426311125923E2L,
00628  8.997558632489032902250523945248208224445E2L,
00629  3.798025197699757225978410230530640879762E3L,
00630  9.113203668683080975637043118209210146846E3L,
00631  1.203285891339933238608683715194034900149E4L,
00632  8.100647057919140328536743641735339740855E3L,
00633  2.383888249907144945837976899822927411769E3L,
00634  2.127493573166454249221983582495245662319E2L
00635 };
00636 #define NRDr4 10
00637 static const long double RDr4[NRDr4 + 1] =
00638 {
00639 -3.303141981514540274165450687270180479586E-1L,
00640 -1.353768629363605300707949368917687066724E1L,
00641 -2.206127630303621521950193783894598987033E2L,
00642 -1.861800338758066696514480386180875607204E3L,
00643 -8.889048775872605708249140016201753255599E3L,
00644 -2.465888106627948210478692168261494857089E4L,
00645 -3.934642211710774494879042116768390014289E4L,
00646 -3.455077258242252974937480623730228841003E4L,
00647 -1.524083977439690284820586063729912653196E4L,
00648 -2.810541887397984804237552337349093953857E3L,
00649 -1.343929553541159933824901621702567066156E2L
00650  /* 1.0E0 */
00651 };
00652 
00653 /* erfc(1/x) = 1/x exp(-1/x^2 - 0.5625 + R(1/x^2))
00654    1/4 <= 1/x < 3/8
00655    Peak relative error 8.4e-37  */
00656 #define NRNr3 11
00657 static const long double RNr3[NRNr3 + 1] =
00658 {
00659 -1.952401126551202208698629992497306292987E-6L,
00660 -2.130881743066372952515162564941682716125E-4L,
00661 -8.376493958090190943737529486107282224387E-3L,
00662 -1.650592646560987700661598877522831234791E-1L,
00663 -1.839290818933317338111364667708678163199E0L,
00664 -1.216278715570882422410442318517814388470E1L,
00665 -4.818759344462360427612133632533779091386E1L,
00666 -1.120994661297476876804405329172164436784E2L,
00667 -1.452850765662319264191141091859300126931E2L,
00668 -9.485207851128957108648038238656777241333E1L,
00669 -2.563663855025796641216191848818620020073E1L,
00670 -1.787995944187565676837847610706317833247E0L
00671 };
00672 #define NRDr3 10
00673 static const long double RDr3[NRDr3 + 1] =
00674 {
00675  1.979130686770349481460559711878399476903E-4L,
00676  1.156941716128488266238105813374635099057E-2L,
00677  2.752657634309886336431266395637285974292E-1L,
00678  3.482245457248318787349778336603569327521E0L,
00679  2.569347069372696358578399521203959253162E1L,
00680  1.142279000180457419740314694631879921561E2L,
00681  3.056503977190564294341422623108332700840E2L,
00682  4.780844020923794821656358157128719184422E2L,
00683  4.105972727212554277496256802312730410518E2L,
00684  1.724072188063746970865027817017067646246E2L,
00685  2.815939183464818198705278118326590370435E1L
00686  /* 1.0E0 */
00687 };
00688 
00689 /* erfc(1/x) = 1/x exp(-1/x^2 - 0.5625 + R(1/x^2))
00690    1/8 <= 1/x < 1/4
00691    Peak relative error 1.5e-36  */
00692 #define NRNr2 11
00693 static const long double RNr2[NRNr2 + 1] =
00694 {
00695 -2.638914383420287212401687401284326363787E-8L,
00696 -3.479198370260633977258201271399116766619E-6L,
00697 -1.783985295335697686382487087502222519983E-4L,
00698 -4.777876933122576014266349277217559356276E-3L,
00699 -7.450634738987325004070761301045014986520E-2L,
00700 -7.068318854874733315971973707247467326619E-1L,
00701 -4.113919921935944795764071670806867038732E0L,
00702 -1.440447573226906222417767283691888875082E1L,
00703 -2.883484031530718428417168042141288943905E1L,
00704 -2.990886974328476387277797361464279931446E1L,
00705 -1.325283914915104866248279787536128997331E1L,
00706 -1.572436106228070195510230310658206154374E0L
00707 };
00708 #define NRDr2 10
00709 static const long double RDr2[NRDr2 + 1] =
00710 {
00711  2.675042728136731923554119302571867799673E-6L,
00712  2.170997868451812708585443282998329996268E-4L,
00713  7.249969752687540289422684951196241427445E-3L,
00714  1.302040375859768674620410563307838448508E-1L,
00715  1.380202483082910888897654537144485285549E0L,
00716  8.926594113174165352623847870299170069350E0L,
00717  3.521089584782616472372909095331572607185E1L,
00718  8.233547427533181375185259050330809105570E1L,
00719  1.072971579885803033079469639073292840135E2L,
00720  6.943803113337964469736022094105143158033E1L,
00721  1.775695341031607738233608307835017282662E1L
00722  /* 1.0E0 */
00723 };
00724 
00725 /* erfc(1/x) = 1/x exp(-1/x^2 - 0.5625 + R(1/x^2))
00726    1/128 <= 1/x < 1/8
00727    Peak relative error 2.2e-36  */
00728 #define NRNr1 9
00729 static const long double RNr1[NRNr1 + 1] =
00730 {
00731 -4.250780883202361946697751475473042685782E-8L,
00732 -5.375777053288612282487696975623206383019E-6L,
00733 -2.573645949220896816208565944117382460452E-4L,
00734 -6.199032928113542080263152610799113086319E-3L,
00735 -8.262721198693404060380104048479916247786E-2L,
00736 -6.242615227257324746371284637695778043982E-1L,
00737 -2.609874739199595400225113299437099626386E0L,
00738 -5.581967563336676737146358534602770006970E0L,
00739 -5.124398923356022609707490956634280573882E0L,
00740 -1.290865243944292370661544030414667556649E0L
00741 };
00742 #define NRDr1 8
00743 static const long double RDr1[NRDr1 + 1] =
00744 {
00745  4.308976661749509034845251315983612976224E-6L,
00746  3.265390126432780184125233455960049294580E-4L,
00747  9.811328839187040701901866531796570418691E-3L,
00748  1.511222515036021033410078631914783519649E-1L,
00749  1.289264341917429958858379585970225092274E0L,
00750  6.147640356182230769548007536914983522270E0L,
00751  1.573966871337739784518246317003956180750E1L,
00752  1.955534123435095067199574045529218238263E1L,
00753  9.472613121363135472247929109615785855865E0L
00754   /* 1.0E0 */
00755 };
00756 
00757 
00758 #ifdef __STDC__
00759 long double
00760 __erfl (long double x)
00761 #else
00762 double
00763 __erfl (x)
00764      long double x;
00765 #endif
00766 {
00767   long double a, y, z;
00768   int32_t i, ix, sign;
00769   ieee854_long_double_shape_type u;
00770 
00771   u.value = x;
00772   sign = u.parts32.w0;
00773   ix = sign & 0x7fffffff;
00774 
00775   if (ix >= 0x7ff00000)
00776     {                       /* erf(nan)=nan */
00777       i = ((sign & 0xfff00000) >> 31) << 1;
00778       return (long double) (1 - i) + one / x;    /* erf(+-inf)=+-1 */
00779     }
00780 
00781   if (ix >= 0x3ff00000) /* |x| >= 1.0 */
00782     {
00783       y = __erfcl (x);
00784       return (one - y);
00785       /*    return (one - __erfcl (x)); */
00786     }
00787   u.parts32.w0 = ix;
00788   a = u.value;
00789   z = x * x;
00790   if (ix < 0x3fec0000)  /* a < 0.875 */
00791     {
00792       if (ix < 0x3c600000) /* |x|<2**-57 */
00793        {
00794          if (ix < 0x00800000)
00795            {
00796              /* erf (-0) = -0.  Unfortunately, for IBM extended double
00797                0.125 * (8.0 * x + efx8 * x) for x = -0 evaluates to 0.  */
00798              if (x == 0)
00799               return x;
00800              return 0.125 * (8.0 * x + efx8 * x);       /*avoid underflow */
00801            }
00802          return x + efx * x;
00803        }
00804       y = a + a * neval (z, TN1, NTN1) / deval (z, TD1, NTD1);
00805     }
00806   else
00807     {
00808       a = a - one;
00809       y = erf_const + neval (a, TN2, NTN2) / deval (a, TD2, NTD2);
00810     }
00811 
00812   if (sign & 0x80000000) /* x < 0 */
00813     y = -y;
00814   return( y );
00815 }
00816 
00817 long_double_symbol (libm, __erfl, erfl);
00818 #ifdef __STDC__
00819      long double
00820      __erfcl (long double x)
00821 #else
00822      long double
00823      __erfcl (x)
00824      double
00825        x;
00826 #endif
00827 {
00828   long double y, z, p, r;
00829   int32_t i, ix, sign;
00830   ieee854_long_double_shape_type u;
00831 
00832   u.value = x;
00833   sign = u.parts32.w0;
00834   ix = sign & 0x7fffffff;
00835   u.parts32.w0 = ix;
00836 
00837   if (ix >= 0x7ff00000)
00838     {                       /* erfc(nan)=nan */
00839       /* erfc(+-inf)=0,2 */
00840       return (long double) (((u_int32_t) sign >> 31) << 1) + one / x;
00841     }
00842 
00843   if (ix < 0x3fd00000) /* |x| <1/4 */
00844     {
00845       if (ix < 0x38d00000) /* |x|<2**-114 */
00846        return one - x;
00847       return one - __erfl (x);
00848     }
00849   if (ix < 0x3ff40000) /* 1.25 */
00850     {
00851       x = u.value;
00852       i = 8.0 * x;
00853       switch (i)
00854        {
00855        case 2:
00856          z = x - 0.25L;
00857          y = C13b + z * neval (z, RNr13, NRNr13) / deval (z, RDr13, NRDr13);
00858          y += C13a;
00859          break;
00860        case 3:
00861          z = x - 0.375L;
00862          y = C14b + z * neval (z, RNr14, NRNr14) / deval (z, RDr14, NRDr14);
00863          y += C14a;
00864          break;
00865        case 4:
00866          z = x - 0.5L;
00867          y = C15b + z * neval (z, RNr15, NRNr15) / deval (z, RDr15, NRDr15);
00868          y += C15a;
00869          break;
00870        case 5:
00871          z = x - 0.625L;
00872          y = C16b + z * neval (z, RNr16, NRNr16) / deval (z, RDr16, NRDr16);
00873          y += C16a;
00874          break;
00875        case 6:
00876          z = x - 0.75L;
00877          y = C17b + z * neval (z, RNr17, NRNr17) / deval (z, RDr17, NRDr17);
00878          y += C17a;
00879          break;
00880        case 7:
00881          z = x - 0.875L;
00882          y = C18b + z * neval (z, RNr18, NRNr18) / deval (z, RDr18, NRDr18);
00883          y += C18a;
00884          break;
00885        case 8:
00886          z = x - 1.0L;
00887          y = C19b + z * neval (z, RNr19, NRNr19) / deval (z, RDr19, NRDr19);
00888          y += C19a;
00889          break;
00890        case 9:
00891          z = x - 1.125L;
00892          y = C20b + z * neval (z, RNr20, NRNr20) / deval (z, RDr20, NRDr20);
00893          y += C20a;
00894          break;
00895        }
00896       if (sign & 0x80000000)
00897        y = 2.0L - y;
00898       return y;
00899     }
00900   /* 1.25 < |x| < 107 */
00901   if (ix < 0x405ac000)
00902     {
00903       /* x < -9 */
00904       if ((ix >= 0x40220000) && (sign & 0x80000000))
00905        return two - tiny;
00906 
00907       x = fabsl (x);
00908       z = one / (x * x);
00909       i = 8.0 / x;
00910       switch (i)
00911        {
00912        default:
00913        case 0:
00914          p = neval (z, RNr1, NRNr1) / deval (z, RDr1, NRDr1);
00915          break;
00916        case 1:
00917          p = neval (z, RNr2, NRNr2) / deval (z, RDr2, NRDr2);
00918          break;
00919        case 2:
00920          p = neval (z, RNr3, NRNr3) / deval (z, RDr3, NRDr3);
00921          break;
00922        case 3:
00923          p = neval (z, RNr4, NRNr4) / deval (z, RDr4, NRDr4);
00924          break;
00925        case 4:
00926          p = neval (z, RNr5, NRNr5) / deval (z, RDr5, NRDr5);
00927          break;
00928        case 5:
00929          p = neval (z, RNr6, NRNr6) / deval (z, RDr6, NRDr6);
00930          break;
00931        case 6:
00932          p = neval (z, RNr7, NRNr7) / deval (z, RDr7, NRDr7);
00933          break;
00934        case 7:
00935          p = neval (z, RNr8, NRNr8) / deval (z, RDr8, NRDr8);
00936          break;
00937        }
00938       u.value = x;
00939       u.parts32.w3 = 0;
00940       u.parts32.w2 &= 0xffffe000;
00941       z = u.value;
00942       r = __ieee754_expl (-z * z - 0.5625) *
00943        __ieee754_expl ((z - x) * (z + x) + p);
00944       if ((sign & 0x80000000) == 0)
00945        return r / x;
00946       else
00947        return two - r / x;
00948     }
00949   else
00950     {
00951       if ((sign & 0x80000000) == 0)
00952        return tiny * tiny;
00953       else
00954        return two - tiny;
00955     }
00956 }
00957 
00958 long_double_symbol (libm, __erfcl, erfcl);