e_j1l.c

glibc 2.9,最新版的C语言库函数

  -1.059928728869218962607068840646564457980E-2L,
  -1.212070036005832342565792241385459023801E-1L,
  -6.688350110633603958684302153362735625156E-1L,
  -1.793587878197360221340277951304429821582E0L,
  -2.225407682237197485644647380483725045326E0L,
  -1.123402135458940189438898496348239744403E0L,
  -1.679187241566347077204805190763597299805E-1L,
  -1.458550613639093752909985189067233504148E-3L,
};
#define NQ2_2r3D 8
static const long double Q2_2r3D[NQ2_2r3D + 1] = {
  5.415024336507980465169023996403597916115E-5L,
  4.179246497380453022046357404266022870788E-3L,
  1.136306384261959483095442402929502368598E-1L,
  1.422640343719842213484515445393284072830E0L,
  8.968786703393158374728850922289204805764E0L,
  2.914542473339246127533384118781216495934E1L,
  4.781605421020380669870197378210457054685E1L,
  3.693865837171883152382820584714795072937E1L,
  1.153220502744204904763115556224395893076E1L,
  /* 1.000000000000000000000000000000000000000E0 */
};


/* Evaluate P[n] x^n  +  P[n-1] x^(n-1)  +  ...  +  P[0] */

static long double
neval (long double x, const long double *p, int n)
{
  long double y;

  p += n;
  y = *p--;
  do
    {
      y = y * x + *p--;
    }
  while (--n > 0);
  return y;
}


/* Evaluate x^n+1  +  P[n] x^(n)  +  P[n-1] x^(n-1)  +  ...  +  P[0] */

static long double
deval (long double x, const long double *p, int n)
{
  long double y;

  p += n;
  y = x + *p--;
  do
    {
      y = y * x + *p--;
    }
  while (--n > 0);
  return y;
}


/* Bessel function of the first kind, order one.  */

long double
__ieee754_j1l (long double x)
{
  long double xx, xinv, z, p, q, c, s, cc, ss;

  if (! __finitel (x))
    {
      if (x != x)
	return x;
      else
	return 0.0L;
    }
  if (x == 0.0L)
    return x;
  xx = fabsl (x);
  if (xx <= 2.0L)
    {
      /* 0 <= x <= 2 */
      z = xx * xx;
      p = xx * z * neval (z, J0_2N, NJ0_2N) / deval (z, J0_2D, NJ0_2D);
      p += 0.5L * xx;
      if (x < 0)
	p = -p;
      return p;
    }

  xinv = 1.0L / xx;
  z = xinv * xinv;
  if (xinv <= 0.25)
    {
      if (xinv <= 0.125)
	{
	  if (xinv <= 0.0625)
	    {
	      p = neval (z, P16_IN, NP16_IN) / deval (z, P16_ID, NP16_ID);
	      q = neval (z, Q16_IN, NQ16_IN) / deval (z, Q16_ID, NQ16_ID);
	    }
	  else
	    {
	      p = neval (z, P8_16N, NP8_16N) / deval (z, P8_16D, NP8_16D);
	      q = neval (z, Q8_16N, NQ8_16N) / deval (z, Q8_16D, NQ8_16D);
	    }
	}
      else if (xinv <= 0.1875)
	{
	  p = neval (z, P5_8N, NP5_8N) / deval (z, P5_8D, NP5_8D);
	  q = neval (z, Q5_8N, NQ5_8N) / deval (z, Q5_8D, NQ5_8D);
	}
      else
	{
	  p = neval (z, P4_5N, NP4_5N) / deval (z, P4_5D, NP4_5D);
	  q = neval (z, Q4_5N, NQ4_5N) / deval (z, Q4_5D, NQ4_5D);
	}
    }				/* .25 */
  else /* if (xinv <= 0.5) */
    {
      if (xinv <= 0.375)
	{
	  if (xinv <= 0.3125)
	    {
	      p = neval (z, P3r2_4N, NP3r2_4N) / deval (z, P3r2_4D, NP3r2_4D);
	      q = neval (z, Q3r2_4N, NQ3r2_4N) / deval (z, Q3r2_4D, NQ3r2_4D);
	    }
	  else
	    {
	      p = neval (z, P2r7_3r2N, NP2r7_3r2N)
		  / deval (z, P2r7_3r2D, NP2r7_3r2D);
	      q = neval (z, Q2r7_3r2N, NQ2r7_3r2N)
		  / deval (z, Q2r7_3r2D, NQ2r7_3r2D);
	    }
	}
      else if (xinv <= 0.4375)
	{
	  p = neval (z, P2r3_2r7N, NP2r3_2r7N)
	      / deval (z, P2r3_2r7D, NP2r3_2r7D);
	  q = neval (z, Q2r3_2r7N, NQ2r3_2r7N)
	      / deval (z, Q2r3_2r7D, NQ2r3_2r7D);
	}
      else
	{
	  p = neval (z, P2_2r3N, NP2_2r3N) / deval (z, P2_2r3D, NP2_2r3D);
	  q = neval (z, Q2_2r3N, NQ2_2r3N) / deval (z, Q2_2r3D, NQ2_2r3D);
	}
    }
  p = 1.0L + z * p;
  q = z * q;
  q = q * xinv + 0.375L * xinv;
  /* X = x - 3 pi/4
     cos(X) = cos(x) cos(3 pi/4) + sin(x) sin(3 pi/4)
     = 1/sqrt(2) * (-cos(x) + sin(x))
     sin(X) = sin(x) cos(3 pi/4) - cos(x) sin(3 pi/4)
     = -1/sqrt(2) * (sin(x) + cos(x))
     cf. Fdlibm.  */
  __sincosl (xx, &s, &c);
  ss = -s - c;
  cc = s - c;
  z = __cosl (xx + xx);
  if ((s * c) > 0)
    cc = z / ss;
  else
    ss = z / cc;
  z = ONEOSQPI * (p * cc - q * ss) / __ieee754_sqrtl (xx);
  if (x < 0)
    z = -z;
  return z;
}


/* Y1(x) = 2/pi * (log(x) * J1(x) - 1/x) + x R(x^2)
   Peak relative error 6.2e-38
   0 <= x <= 2   */
#define NY0_2N 7
static long double Y0_2N[NY0_2N + 1] = {
  -6.804415404830253804408698161694720833249E19L,
  1.805450517967019908027153056150465849237E19L,
  -8.065747497063694098810419456383006737312E17L,
  1.401336667383028259295830955439028236299E16L,
  -1.171654432898137585000399489686629680230E14L,
  5.061267920943853732895341125243428129150E11L,
  -1.096677850566094204586208610960870217970E9L,
  9.541172044989995856117187515882879304461E5L,
};
#define NY0_2D 7
static long double Y0_2D[NY0_2D + 1] = {
  3.470629591820267059538637461549677594549E20L,
  4.120796439009916326855848107545425217219E18L,
  2.477653371652018249749350657387030814542E16L,
  9.954678543353888958177169349272167762797E13L,
  2.957927997613630118216218290262851197754E11L,
  6.748421382188864486018861197614025972118E8L,
  1.173453425218010888004562071020305709319E6L,
  1.450335662961034949894009554536003377187E3L,
  /* 1.000000000000000000000000000000000000000E0 */
};


/* Bessel function of the second kind, order one.  */

long double
__ieee754_y1l (long double x)
{
  long double xx, xinv, z, p, q, c, s, cc, ss;

  if (! __finitel (x))
    {
      if (x != x)
	return x;
      else
	return 0.0L;
    }
  if (x <= 0.0L)
    {
      if (x < 0.0L)
	return (zero / (zero * x));
      return -HUGE_VALL + x;
    }
  xx = fabsl (x);
  if (xx <= 2.0L)
    {
      /* 0 <= x <= 2 */
      z = xx * xx;
      p = xx * neval (z, Y0_2N, NY0_2N) / deval (z, Y0_2D, NY0_2D);
      p = -TWOOPI / xx + p;
      p = TWOOPI * __ieee754_logl (x) * __ieee754_j1l (x) + p;
      return p;
    }

  xinv = 1.0L / xx;
  z = xinv * xinv;
  if (xinv <= 0.25)
    {
      if (xinv <= 0.125)
	{
	  if (xinv <= 0.0625)
	    {
	      p = neval (z, P16_IN, NP16_IN) / deval (z, P16_ID, NP16_ID);
	      q = neval (z, Q16_IN, NQ16_IN) / deval (z, Q16_ID, NQ16_ID);
	    }
	  else
	    {
	      p = neval (z, P8_16N, NP8_16N) / deval (z, P8_16D, NP8_16D);
	      q = neval (z, Q8_16N, NQ8_16N) / deval (z, Q8_16D, NQ8_16D);
	    }
	}
      else if (xinv <= 0.1875)
	{
	  p = neval (z, P5_8N, NP5_8N) / deval (z, P5_8D, NP5_8D);
	  q = neval (z, Q5_8N, NQ5_8N) / deval (z, Q5_8D, NQ5_8D);
	}
      else
	{
	  p = neval (z, P4_5N, NP4_5N) / deval (z, P4_5D, NP4_5D);
	  q = neval (z, Q4_5N, NQ4_5N) / deval (z, Q4_5D, NQ4_5D);
	}
    }				/* .25 */
  else /* if (xinv <= 0.5) */
    {
      if (xinv <= 0.375)
	{
	  if (xinv <= 0.3125)
	    {
	      p = neval (z, P3r2_4N, NP3r2_4N) / deval (z, P3r2_4D, NP3r2_4D);
	      q = neval (z, Q3r2_4N, NQ3r2_4N) / deval (z, Q3r2_4D, NQ3r2_4D);
	    }
	  else
	    {
	      p = neval (z, P2r7_3r2N, NP2r7_3r2N)
		  / deval (z, P2r7_3r2D, NP2r7_3r2D);
	      q = neval (z, Q2r7_3r2N, NQ2r7_3r2N)
		  / deval (z, Q2r7_3r2D, NQ2r7_3r2D);
	    }
	}
      else if (xinv <= 0.4375)
	{
	  p = neval (z, P2r3_2r7N, NP2r3_2r7N)
	      / deval (z, P2r3_2r7D, NP2r3_2r7D);
	  q = neval (z, Q2r3_2r7N, NQ2r3_2r7N)
	      / deval (z, Q2r3_2r7D, NQ2r3_2r7D);
	}
      else
	{
	  p = neval (z, P2_2r3N, NP2_2r3N) / deval (z, P2_2r3D, NP2_2r3D);
	  q = neval (z, Q2_2r3N, NQ2_2r3N) / deval (z, Q2_2r3D, NQ2_2r3D);
	}
    }
  p = 1.0L + z * p;
  q = z * q;
  q = q * xinv + 0.375L * xinv;
  /* X = x - 3 pi/4
     cos(X) = cos(x) cos(3 pi/4) + sin(x) sin(3 pi/4)
     = 1/sqrt(2) * (-cos(x) + sin(x))
     sin(X) = sin(x) cos(3 pi/4) - cos(x) sin(3 pi/4)
     = -1/sqrt(2) * (sin(x) + cos(x))
     cf. Fdlibm.  */
  __sincosl (xx, &s, &c);
  ss = -s - c;
  cc = s - c;
  z = __cosl (xx + xx);
  if ((s * c) > 0)
    cc = z / ss;
  else
    ss = z / cc;
  z = ONEOSQPI * (p * ss + q * cc) / __ieee754_sqrtl (xx);
  return z;
}