s_erfl.c

Glibc 2.3.2源代码(解压后有100多M)

  3.132269062552392974833215844236160958502E1L,
  4.430131663290563523933419966185230513168E0L
 /* 1.0E0 */
};
/* erfc(1.125) = C20a + C20b to extra precision.  */
static const long double C20a = 0.111602783203125L;
static const long double C20b = 8.9850951672359304215530728365232161564636E-6L;

/* erfc(1/x) = 1/x exp (-1/x^2 - 0.5625 + R(1/x^2))
   7/8 <= 1/x < 1
   Peak relative error 1.4e-35  */
#define NRNr8 9
static const long double RNr8[NRNr8 + 1] =
{
  3.587451489255356250759834295199296936784E1L,
  5.406249749087340431871378009874875889602E2L,
  2.931301290625250886238822286506381194157E3L,
  7.359254185241795584113047248898753470923E3L,
  9.201031849810636104112101947312492532314E3L,
  5.749697096193191467751650366613289284777E3L,
  1.710415234419860825710780802678697889231E3L,
  2.150753982543378580859546706243022719599E2L,
  8.740953582272147335100537849981160931197E0L,
  4.876422978828717219629814794707963640913E-2L
};
#define NRDr8 8
static const long double RDr8[NRDr8 + 1] =
{
  6.358593134096908350929496535931630140282E1L,
  9.900253816552450073757174323424051765523E2L,
  5.642928777856801020545245437089490805186E3L,
  1.524195375199570868195152698617273739609E4L,
  2.113829644500006749947332935305800887345E4L,
  1.526438562626465706267943737310282977138E4L,
  5.561370922149241457131421914140039411782E3L,
  9.394035530179705051609070428036834496942E2L,
  6.147019596150394577984175188032707343615E1L
  /* 1.0E0 */
};

/* erfc(1/x) = 1/x exp (-1/x^2 - 0.5625 + R(1/x^2))
   0.75 <= 1/x <= 0.875
   Peak relative error 2.0e-36  */
#define NRNr7 9
static const long double RNr7[NRNr7 + 1] =
{
 1.686222193385987690785945787708644476545E1L,
 1.178224543567604215602418571310612066594E3L,
 1.764550584290149466653899886088166091093E4L,
 1.073758321890334822002849369898232811561E5L,
 3.132840749205943137619839114451290324371E5L,
 4.607864939974100224615527007793867585915E5L,
 3.389781820105852303125270837910972384510E5L,
 1.174042187110565202875011358512564753399E5L,
 1.660013606011167144046604892622504338313E4L,
 6.700393957480661937695573729183733234400E2L
};
#define NRDr7 9
static const long double RDr7[NRDr7 + 1] =
{
-1.709305024718358874701575813642933561169E3L,
-3.280033887481333199580464617020514788369E4L,
-2.345284228022521885093072363418750835214E5L,
-8.086758123097763971926711729242327554917E5L,
-1.456900414510108718402423999575992450138E6L,
-1.391654264881255068392389037292702041855E6L,
-6.842360801869939983674527468509852583855E5L,
-1.597430214446573566179675395199807533371E5L,
-1.488876130609876681421645314851760773480E4L,
-3.511762950935060301403599443436465645703E2L
 /* 1.0E0 */
};

/* erfc(1/x) = 1/x exp(-1/x^2 - 0.5625 + R(1/x^2))
   5/8 <= 1/x < 3/4
   Peak relative error 1.9e-35  */
#define NRNr6 9
static const long double RNr6[NRNr6 + 1] =
{
 1.642076876176834390623842732352935761108E0L,
 1.207150003611117689000664385596211076662E2L,
 2.119260779316389904742873816462800103939E3L,
 1.562942227734663441801452930916044224174E4L,
 5.656779189549710079988084081145693580479E4L,
 1.052166241021481691922831746350942786299E5L,
 9.949798524786000595621602790068349165758E4L,
 4.491790734080265043407035220188849562856E4L,
 8.377074098301530326270432059434791287601E3L,
 4.506934806567986810091824791963991057083E2L
};
#define NRDr6 9
static const long double RDr6[NRDr6 + 1] =
{
-1.664557643928263091879301304019826629067E2L,
-3.800035902507656624590531122291160668452E3L,
-3.277028191591734928360050685359277076056E4L,
-1.381359471502885446400589109566587443987E5L,
-3.082204287382581873532528989283748656546E5L,
-3.691071488256738343008271448234631037095E5L,
-2.300482443038349815750714219117566715043E5L,
-6.873955300927636236692803579555752171530E4L,
-8.262158817978334142081581542749986845399E3L,
-2.517122254384430859629423488157361983661E2L
 /* 1.00 */
};

/* erfc(1/x) = 1/x exp(-1/x^2 - 0.5625 + R(1/x^2))
   1/2 <= 1/x < 5/8
   Peak relative error 4.6e-36  */
#define NRNr5 10
static const long double RNr5[NRNr5 + 1] =
{
-3.332258927455285458355550878136506961608E-3L,
-2.697100758900280402659586595884478660721E-1L,
-6.083328551139621521416618424949137195536E0L,
-6.119863528983308012970821226810162441263E1L,
-3.176535282475593173248810678636522589861E2L,
-8.933395175080560925809992467187963260693E2L,
-1.360019508488475978060917477620199499560E3L,
-1.075075579828188621541398761300910213280E3L,
-4.017346561586014822824459436695197089916E2L,
-5.857581368145266249509589726077645791341E1L,
-2.077715925587834606379119585995758954399E0L
};
#define NRDr5 9
static const long double RDr5[NRDr5 + 1] =
{
 3.377879570417399341550710467744693125385E-1L,
 1.021963322742390735430008860602594456187E1L,
 1.200847646592942095192766255154827011939E2L,
 7.118915528142927104078182863387116942836E2L,
 2.318159380062066469386544552429625026238E3L,
 4.238729853534009221025582008928765281620E3L,
 4.279114907284825886266493994833515580782E3L,
 2.257277186663261531053293222591851737504E3L,
 5.570475501285054293371908382916063822957E2L,
 5.142189243856288981145786492585432443560E1L
 /* 1.0E0 */
};

/* erfc(1/x) = 1/x exp(-1/x^2 - 0.5625 + R(1/x^2))
   3/8 <= 1/x < 1/2
   Peak relative error 2.0e-36  */
#define NRNr4 10
static const long double RNr4[NRNr4 + 1] =
{
 3.258530712024527835089319075288494524465E-3L,
 2.987056016877277929720231688689431056567E-1L,
 8.738729089340199750734409156830371528862E0L,
 1.207211160148647782396337792426311125923E2L,
 8.997558632489032902250523945248208224445E2L,
 3.798025197699757225978410230530640879762E3L,
 9.113203668683080975637043118209210146846E3L,
 1.203285891339933238608683715194034900149E4L,
 8.100647057919140328536743641735339740855E3L,
 2.383888249907144945837976899822927411769E3L,
 2.127493573166454249221983582495245662319E2L
};
#define NRDr4 10
static const long double RDr4[NRDr4 + 1] =
{
-3.303141981514540274165450687270180479586E-1L,
-1.353768629363605300707949368917687066724E1L,
-2.206127630303621521950193783894598987033E2L,
-1.861800338758066696514480386180875607204E3L,
-8.889048775872605708249140016201753255599E3L,
-2.465888106627948210478692168261494857089E4L,
-3.934642211710774494879042116768390014289E4L,
-3.455077258242252974937480623730228841003E4L,
-1.524083977439690284820586063729912653196E4L,
-2.810541887397984804237552337349093953857E3L,
-1.343929553541159933824901621702567066156E2L
 /* 1.0E0 */
};

/* erfc(1/x) = 1/x exp(-1/x^2 - 0.5625 + R(1/x^2))
   1/4 <= 1/x < 3/8
   Peak relative error 8.4e-37  */
#define NRNr3 11
static const long double RNr3[NRNr3 + 1] =
{
-1.952401126551202208698629992497306292987E-6L,
-2.130881743066372952515162564941682716125E-4L,
-8.376493958090190943737529486107282224387E-3L,
-1.650592646560987700661598877522831234791E-1L,
-1.839290818933317338111364667708678163199E0L,
-1.216278715570882422410442318517814388470E1L,
-4.818759344462360427612133632533779091386E1L,
-1.120994661297476876804405329172164436784E2L,
-1.452850765662319264191141091859300126931E2L,
-9.485207851128957108648038238656777241333E1L,
-2.563663855025796641216191848818620020073E1L,
-1.787995944187565676837847610706317833247E0L
};
#define NRDr3 10
static const long double RDr3[NRDr3 + 1] =
{
 1.979130686770349481460559711878399476903E-4L,
 1.156941716128488266238105813374635099057E-2L,
 2.752657634309886336431266395637285974292E-1L,
 3.482245457248318787349778336603569327521E0L,
 2.569347069372696358578399521203959253162E1L,
 1.142279000180457419740314694631879921561E2L,
 3.056503977190564294341422623108332700840E2L,
 4.780844020923794821656358157128719184422E2L,
 4.105972727212554277496256802312730410518E2L,
 1.724072188063746970865027817017067646246E2L,
 2.815939183464818198705278118326590370435E1L
 /* 1.0E0 */
};

/* erfc(1/x) = 1/x exp(-1/x^2 - 0.5625 + R(1/x^2))
   1/8 <= 1/x < 1/4
   Peak relative error 1.5e-36  */
#define NRNr2 11
static const long double RNr2[NRNr2 + 1] =
{
-2.638914383420287212401687401284326363787E-8L,
-3.479198370260633977258201271399116766619E-6L,
-1.783985295335697686382487087502222519983E-4L,
-4.777876933122576014266349277217559356276E-3L,
-7.450634738987325004070761301045014986520E-2L,
-7.068318854874733315971973707247467326619E-1L,
-4.113919921935944795764071670806867038732E0L,
-1.440447573226906222417767283691888875082E1L,
-2.883484031530718428417168042141288943905E1L,
-2.990886974328476387277797361464279931446E1L,
-1.325283914915104866248279787536128997331E1L,
-1.572436106228070195510230310658206154374E0L
};
#define NRDr2 10
static const long double RDr2[NRDr2 + 1] =
{
 2.675042728136731923554119302571867799673E-6L,
 2.170997868451812708585443282998329996268E-4L,
 7.249969752687540289422684951196241427445E-3L,
 1.302040375859768674620410563307838448508E-1L,
 1.380202483082910888897654537144485285549E0L,
 8.926594113174165352623847870299170069350E0L,
 3.521089584782616472372909095331572607185E1L,
 8.233547427533181375185259050330809105570E1L,
 1.072971579885803033079469639073292840135E2L,
 6.943803113337964469736022094105143158033E1L,
 1.775695341031607738233608307835017282662E1L
 /* 1.0E0 */
};

/* erfc(1/x) = 1/x exp(-1/x^2 - 0.5625 + R(1/x^2))
   1/128 <= 1/x < 1/8
   Peak relative error 2.2e-36  */
#define NRNr1 9
static const long double RNr1[NRNr1 + 1] =
{
-4.250780883202361946697751475473042685782E-8L,
-5.375777053288612282487696975623206383019E-6L,
-2.573645949220896816208565944117382460452E-4L,
-6.199032928113542080263152610799113086319E-3L,
-8.262721198693404060380104048479916247786E-2L,
-6.242615227257324746371284637695778043982E-1L,
-2.609874739199595400225113299437099626386E0L,
-5.581967563336676737146358534602770006970E0L,
-5.124398923356022609707490956634280573882E0L,
-1.290865243944292370661544030414667556649E0L
};
#define NRDr1 8
static const long double RDr1[NRDr1 + 1] =
{
 4.308976661749509034845251315983612976224E-6L,
 3.265390126432780184125233455960049294580E-4L,
 9.811328839187040701901866531796570418691E-3L,
 1.511222515036021033410078631914783519649E-1L,
 1.289264341917429958858379585970225092274E0L,
 6.147640356182230769548007536914983522270E0L,
 1.573966871337739784518246317003956180750E1L,
 1.955534123435095067199574045529218238263E1L,
 9.472613121363135472247929109615785855865E0L
  /* 1.0E0 */
};


#ifdef __STDC__
long double
__erfl (long double x)
#else
double
__erfl (x)
     long double x;
#endif
{
  long double a, y, z;
  int32_t i, ix, sign;
  ieee854_long_double_shape_type u;

  u.value = x;
  sign = u.parts32.w0;
  ix = sign & 0x7fffffff;

  if (ix >= 0x7fff0000)
    {				/* erf(nan)=nan */
      i = ((sign & 0xffff0000) >> 31) << 1;
      return (long double) (1 - i) + one / x;	/* erf(+-inf)=+-1 */
    }

  if (ix >= 0x3fff0000) /* |x| >= 1.0 */
    {
      y = __erfcl (x);
      return (one - y);
      /*    return (one - __erfcl (x)); */
    }
  u.parts32.w0 = ix;
  a = u.value;
  z = x * x;
  if (ix < 0x3ffec000)  /* a < 0.875 */
    {
      if (ix < 0x3fc60000) /* |x|<2**-57 */
	{
	  if (ix < 0x00080000)
	    return 0.125 * (8.0 * x + efx8 * x);	/*avoid underflow */
	  return x + efx * x;
	}
      y = a + a * neval (z, TN1, NTN1) / deval (z, TD1, NTD1);
    }
  else
    {
      a = a - one;
      y = erf_const + neval (a, TN2, NTN2) / deval (a, TD2, NTD2);
    }

  if (sign & 0x80000000) /* x < 0 */
    y = -y;
  return( y );
}

weak_alias (__erfl, erfl)
#ifdef __STDC__
     long double
     __erfcl (long double x)
#else
     long double
     __erfcl (x)
     double
       x;
#endif
{
  long double y, z, p, r;
  int32_t i, ix, sign;
  ieee854_long_double_shape_type u;

  u.value = x;
  sign = u.parts32.w0;
  ix = sign & 0x7fffffff;
  u.parts32.w0 = ix;

  if (ix >= 0x7fff0000)
    {				/* erfc(nan)=nan */
      /* erfc(+-inf)=0,2 */
      return (long double) (((u_int32_t) sign >> 31) << 1) + one / x;
    }

  if (ix < 0x3ffd0000) /* |x| <1/4 */
    {
      if (ix < 0x3f8d0000) /* |x|<2**-114 */
	return one - x;
      return one - __erfl (x);
    }
  if (ix < 0x3fff4000) /* 1.25 */
    {
      x = u.value;
      i = 8.0 * x;
      switch (i)
	{
	case 2:
	  z = x - 0.25L;
	  y = C13b + z * neval (z, RNr13, NRNr13) / deval (z, RDr13, NRDr13);
	  y += C13a;
	  break;
	case 3:
	  z = x - 0.375L;
	  y = C14b + z * neval (z, RNr14, NRNr14) / deval (z, RDr14, NRDr14);
	  y += C14a;
	  break;
	case 4:
	  z = x - 0.5L;
	  y = C15b + z * neval (z, RNr15, NRNr15) / deval (z, RDr15, NRDr15);
	  y += C15a;
	  break;
	case 5:
	  z = x - 0.625L;
	  y = C16b + z * neval (z, RNr16, NRNr16) / deval (z, RDr16, NRDr16);
	  y += C16a;
	  break;
	case 6:
	  z = x - 0.75L;
	  y = C17b + z * neval (z, RNr17, NRNr17) / deval (z, RDr17, NRDr17);
	  y += C17a;
	  break;
	case 7:
	  z = x - 0.875L;
	  y = C18b + z * neval (z, RNr18, NRNr18) / deval (z, RDr18, NRDr18);
	  y += C18a;
	  break;
	case 8:
	  z = x - 1.0L;
	  y = C19b + z * neval (z, RNr19, NRNr19) / deval (z, RDr19, NRDr19);
	  y += C19a;
	  break;
	case 9:
	  z = x - 1.125L;
	  y = C20b + z * neval (z, RNr20, NRNr20) / deval (z, RDr20, NRDr20);
	  y += C20a;
	  break;
	}
      if (sign & 0x80000000)
	y = 2.0L - y;
      return y;
    }
  /* 1.25 < |x| < 107 */
  if (ix < 0x4005ac00)
    {
      /* x < -9 */
      if ((ix >= 0x40022000) && (sign & 0x80000000))
	return two - tiny;

      x = fabsl (x);
      z = one / (x * x);
      i = 8.0 / x;
      switch (i)
	{
	default:
	case 0:
	  p = neval (z, RNr1, NRNr1) / deval (z, RDr1, NRDr1);
	  break;
	case 1:
	  p = neval (z, RNr2, NRNr2) / deval (z, RDr2, NRDr2);
	  break;
	case 2:
	  p = neval (z, RNr3, NRNr3) / deval (z, RDr3, NRDr3);
	  break;
	case 3:
	  p = neval (z, RNr4, NRNr4) / deval (z, RDr4, NRDr4);
	  break;
	case 4:
	  p = neval (z, RNr5, NRNr5) / deval (z, RDr5, NRDr5);
	  break;
	case 5:
	  p = neval (z, RNr6, NRNr6) / deval (z, RDr6, NRDr6);
	  break;
	case 6:
	  p = neval (z, RNr7, NRNr7) / deval (z, RDr7, NRDr7);
	  break;
	case 7:
	  p = neval (z, RNr8, NRNr8) / deval (z, RDr8, NRDr8);
	  break;
	}
      u.value = x;
      u.parts32.w3 = 0;
      u.parts32.w2 &= 0xfe000000;
      z = u.value;
      r = __ieee754_expl (-z * z - 0.5625) *
	__ieee754_expl ((z - x) * (z + x) + p);
      if ((sign & 0x80000000) == 0)
	return r / x;
      else
	return two - r / x;
    }
  else
    {
      if ((sign & 0x80000000) == 0)
	return tiny * tiny;
      else
	return two - tiny;
    }
}

weak_alias (__erfcl, erfcl)