e_j1l.cpp

来自「这是整套横扫千军3D版游戏的源码」· C++ 代码 · 共 655 行 · 第 1/2 页

  1.642719430496086618401091544113220340094E-3l,
  2.228688005300803935928733750456396149104E-2l,
  1.142773760804150921573259605730018327162E-1l,
  1.755576530055079253910829652698703791957E-1l,
  3.218803858282095929559165965353784980613E-2l,
};
#ifdef __STDC__
static const Extended ps5[6] = {
#else
static Extended ps5[6] = {
#endif
  3.685108812227721334719884358034713967557E-6l,
  4.069102509511177498808856515005792027639E-4l,
  1.449728676496155025507893322405597039816E-2l,
  2.058869213229520086582695850441194363103E-1l,
  1.164890985918737148968424972072751066553E0l,
  2.274776933457009446573027260373361586841E0l,
  /*  1.000000000000000000000000000000000000000E0l,*/
};

/* J1(x) cosX + Y1(x) sinX  =  sqrt( 2/(pi x)) P1(x)
   P1(x) = 1 + z^2 R(z^2), z=1/x
   2.85711669921875l <= x <= 4.54541015625l
   Peak relative error 6.5e-21l  */
#ifdef __STDC__
static const Extended pr3[7] = {
#else
static Extended pr3[7] = {
#endif
  1.265251153957366716825382654273326407972E-5l,
  8.031057269201324914127680782288352574567E-4l,
  1.581648121115028333661412169396282881035E-2l,
  1.179534658087796321928362981518645033967E-1l,
  3.227936912780465219246440724502790727866E-1l,
  2.559223765418386621748404398017602935764E-1l,
  2.277136933287817911091370397134882441046E-2l,
};
#ifdef __STDC__
static const Extended ps3[6] = {
#else
static Extended ps3[6] = {
#endif
  1.079681071833391818661952793568345057548E-4l,
  6.986017817100477138417481463810841529026E-3l,
  1.429403701146942509913198539100230540503E-1l,
  1.148392024337075609460312658938700765074E0l,
  3.643663015091248720208251490291968840882E0l,
  3.990702269032018282145100741746633960737E0l,
  /* 1.000000000000000000000000000000000000000E0l, */
};

/* J1(x) cosX + Y1(x) sinX  =  sqrt( 2/(pi x)) P1(x)
   P1(x) = 1 + z^2 R(z^2), z=1/x
   2 <= x <= 2.85711669921875l
   Peak relative error 3.5e-21l  */
#ifdef __STDC__
static const Extended pr2[7] = {
#else
static Extended pr2[7] = {
#endif
  2.795623248568412225239401141338714516445E-4l,
  1.092578168441856711925254839815430061135E-2l,
  1.278024620468953761154963591853679640560E-1l,
  5.469680473691500673112904286228351988583E-1l,
  8.313769490922351300461498619045639016059E-1l,
  3.544176317308370086415403567097130611468E-1l,
  1.604142674802373041247957048801599740644E-2l,
};
#ifdef __STDC__
static const Extended ps2[6] = {
#else
static Extended ps2[6] = {
#endif
  2.385605161555183386205027000675875235980E-3l,
  9.616778294482695283928617708206967248579E-2l,
  1.195215570959693572089824415393951258510E0l,
  5.718412857897054829999458736064922974662E0l,
  1.065626298505499086386584642761602177568E1l,
  6.809140730053382188468983548092322151791E0l,
 /* 1.000000000000000000000000000000000000000E0l, */
};


#ifdef __STDC__
static Extended
pone (Extended x)
#else
static Extended
pone (x)
     Extended x;
#endif
{
#ifdef __STDC__
  const Extended *p, *q;
#else
  Extended *p, *q;
#endif
  Extended z, r, s;
  int32_t ix;
  u_int32_t se, i0, i1;

  GET_LDOUBLE_WORDS (se, i0, i1, x);
  ix = se & 0x7fff;
  if (ix >= 0x4002) /* x >= 8 */
    {
      p = pr8;
      q = ps8;
    }
  else
    {
      i1 = (ix << 16) | (i0 >> 16);
      if (i1 >= 0x40019174)	/* x >= 4.54541015625l */
	{
	  p = pr5;
	  q = ps5;
	}
      else if (i1 >= 0x4000b6db)	/* x >= 2.85711669921875l */
	{
	  p = pr3;
	  q = ps3;
	}
      else if (ix >= 0x4000)	/* x better be >= 2 */
	{
	  p = pr2;
	  q = ps2;
	}
    }
  z = one / (x * x);
 r = p[0] + z * (p[1] +
		 z * (p[2] + z * (p[3] + z * (p[4] + z * (p[5] + z * p[6])))));
 s = q[0] + z * (q[1] + z * (q[2] + z * (q[3] + z * (q[4] + z * (q[5] + z)))));
  return one + z * r / s;
}


/* For x >= 8, the asymptotic expansions of qone is
 *	3/8 s - 105/1024 s^3 - ..., where s = 1/x.
 * We approximate pone by
 * 	qone(x) = s*(0.375l + (R/S))
 */

/* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
   Q1(x) = z(.375 + z^2 R(z^2)), z=1/x
   8 <= x <= inf
   Peak relative error 8.3e-22l */

#ifdef __STDC__
static const Extended qr8[7] = {
#else
static Extended qr8[7] = {
#endif
  -5.691925079044209246015366919809404457380E-10l,
  -1.632587664706999307871963065396218379137E-7l,
  -1.577424682764651970003637263552027114600E-5l,
  -6.377627959241053914770158336842725291713E-4l,
  -1.087408516779972735197277149494929568768E-2l,
  -6.854943629378084419631926076882330494217E-2l,
  -1.055448290469180032312893377152490183203E-1l,
};
#ifdef __STDC__
static const Extended qs8[7] = {
#else
static Extended qs8[7] = {
#endif
  5.550982172325019811119223916998393907513E-9l,
  1.607188366646736068460131091130644192244E-6l,
  1.580792530091386496626494138334505893599E-4l,
  6.617859900815747303032860443855006056595E-3l,
  1.212840547336984859952597488863037659161E-1l,
  9.017885953937234900458186716154005541075E-1l,
  2.201114489712243262000939120146436167178E0l,
  /* 1.000000000000000000000000000000000000000E0l, */
};

/* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
   Q1(x) = z(.375 + z^2 R(z^2)), z=1/x
   4.54541015625l <= x <= 8
   Peak relative error 4.1e-22l */
#ifdef __STDC__
static const Extended qr5[7] = {
#else
static Extended qr5[7] = {
#endif
  -6.719134139179190546324213696633564965983E-8l,
  -9.467871458774950479909851595678622044140E-6l,
  -4.429341875348286176950914275723051452838E-4l,
  -8.539898021757342531563866270278505014487E-3l,
  -6.818691805848737010422337101409276287170E-2l,
  -1.964432669771684034858848142418228214855E-1l,
  -1.333896496989238600119596538299938520726E-1l,
};
#ifdef __STDC__
static const Extended qs5[7] = {
#else
static Extended qs5[7] = {
#endif
  6.552755584474634766937589285426911075101E-7l,
  9.410814032118155978663509073200494000589E-5l,
  4.561677087286518359461609153655021253238E-3l,
  9.397742096177905170800336715661091535805E-2l,
  8.518538116671013902180962914473967738771E-1l,
  3.177729183645800174212539541058292579009E0l,
  4.006745668510308096259753538973038902990E0l,
  /* 1.000000000000000000000000000000000000000E0l, */
};

/* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
   Q1(x) = z(.375 + z^2 R(z^2)), z=1/x
   2.85711669921875l <= x <= 4.54541015625l
   Peak relative error 2.2e-21l */
#ifdef __STDC__
static const Extended qr3[7] = {
#else
static Extended qr3[7] = {
#endif
  -3.618746299358445926506719188614570588404E-6l,
  -2.951146018465419674063882650970344502798E-4l,
  -7.728518171262562194043409753656506795258E-3l,
  -8.058010968753999435006488158237984014883E-2l,
  -3.356232856677966691703904770937143483472E-1l,
  -4.858192581793118040782557808823460276452E-1l,
  -1.592399251246473643510898335746432479373E-1l,
};
#ifdef __STDC__
static const Extended qs3[7] = {
#else
static Extended qs3[7] = {
#endif
  3.529139957987837084554591421329876744262E-5l,
  2.973602667215766676998703687065066180115E-3l,
  8.273534546240864308494062287908662592100E-2l,
  9.613359842126507198241321110649974032726E-1l,
  4.853923697093974370118387947065402707519E0l,
  1.002671608961669247462020977417828796933E1l,
  7.028927383922483728931327850683151410267E0l,
  /* 1.000000000000000000000000000000000000000E0l, */
};

/* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
   Q1(x) = z(.375 + z^2 R(z^2)), z=1/x
   2 <= x <= 2.85711669921875l
   Peak relative error 6.9e-22l */
#ifdef __STDC__
static const Extended qr2[7] = {
#else
static Extended qr2[7] = {
#endif
  -1.372751603025230017220666013816502528318E-4l,
  -6.879190253347766576229143006767218972834E-3l,
  -1.061253572090925414598304855316280077828E-1l,
  -6.262164224345471241219408329354943337214E-1l,
  -1.423149636514768476376254324731437473915E0l,
  -1.087955310491078933531734062917489870754E0l,
  -1.826821119773182847861406108689273719137E-1l,
};
#ifdef __STDC__
static const Extended qs2[7] = {
#else
static Extended qs2[7] = {
#endif
  1.338768933634451601814048220627185324007E-3l,
  7.071099998918497559736318523932241901810E-2l,
  1.200511429784048632105295629933382142221E0l,
  8.327301713640367079030141077172031825276E0l,
  2.468479301872299311658145549931764426840E1l,
  2.961179686096262083509383820557051621644E1l,
  1.201402313144305153005639494661767354977E1l,
 /* 1.000000000000000000000000000000000000000E0l, */
};


#ifdef __STDC__
static Extended
qone (Extended x)
#else
static Extended
qone (x)
     Extended x;
#endif
{
#ifdef __STDC__
  const Extended *p, *q;
#else
  Extended *p, *q;
#endif
  static Extended s, r, z;
  int32_t ix;
  u_int32_t se, i0, i1;

  GET_LDOUBLE_WORDS (se, i0, i1, x);
  ix = se & 0x7fff;
  if (ix >= 0x4002)		/* x >= 8 */
    {
      p = qr8;
      q = qs8;
    }
  else
    {
      i1 = (ix << 16) | (i0 >> 16);
      if (i1 >= 0x40019174)	/* x >= 4.54541015625l */
	{
	  p = qr5;
	  q = qs5;
	}
      else if (i1 >= 0x4000b6db)	/* x >= 2.85711669921875l */
	{
	  p = qr3;
	  q = qs3;
	}
      else if (ix >= 0x4000)	/* x better be >= 2 */
	{
	  p = qr2;
	  q = qs2;
	}
    }
  z = one / (x * x);
  r =
    p[0] + z * (p[1] +
		z * (p[2] + z * (p[3] + z * (p[4] + z * (p[5] + z * p[6])))));
  s =
    q[0] + z * (q[1] +
		z * (q[2] +
		     z * (q[3] + z * (q[4] + z * (q[5] + z * (q[6] + z))))));
  return (.375 + z * r / s) / x;
}
}