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