e_j0l.c

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

  3.539652320122661637429658698954748337223E-1L,
  9.571721066119617436343740541777014319695E-1L,
  1.196258777828426399432550698612171955305E0L,
  6.069388659458926158392384709893753793967E-1L,
  9.026746127269713176512359976978248763621E-2L,
  5.317668723070450235320878117210807236375E-4L,
};
#define NQ2_2r3D 8
static const long double Q2_2r3D[NQ2_2r3D + 1] = {
  3.846924354014260866793741072933159380158E-5L,
  3.017562820057704325510067178327449946763E-3L,
  8.356305620686867949798885808540444210935E-2L,
  1.068314930499906838814019619594424586273E0L,
  6.900279623894821067017966573640732685233E0L,
  2.307667390886377924509090271780839563141E1L,
  3.921043465412723970791036825401273528513E1L,
  3.167569478939719383241775717095729233436E1L,
  1.051023841699200920276198346301543665909E1L,
 /* 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 zero.  */

long double
__ieee754_j0l (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 1.0L;

  xx = fabsl (x);
  if (xx <= 2.0L)
    {
      /* 0 <= x <= 2 */
      z = xx * xx;
      p = z * z * neval (z, J0_2N, NJ0_2N) / deval (z, J0_2D, NJ0_2D);
      p -= 0.25L * z;
      p += 1.0L;
      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 * xinv * q;
  q = q - 0.125L * xinv;
  /* X = x - pi/4
     cos(X) = cos(x) cos(pi/4) + sin(x) sin(pi/4)
     = 1/sqrt(2) * (cos(x) + sin(x))
     sin(X) = sin(x) cos(pi/4) - cos(x) sin(pi/4)
     = 1/sqrt(2) * (sin(x) - cos(x))
     sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
     cf. Fdlibm.  */
  c = cosl (xx);
  s = sinl (xx);
  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) / sqrtl (xx);
  return z;
}


/* Y0(x) = 2/pi * log(x) * J0(x) + R(x^2)
   Peak absolute error 1.7e-36 (relative where Y0 > 1)
   0 <= x <= 2   */
#define NY0_2N 7
static long double Y0_2N[NY0_2N + 1] = {
 -1.062023609591350692692296993537002558155E19L,
  2.542000883190248639104127452714966858866E19L,
 -1.984190771278515324281415820316054696545E18L,
  4.982586044371592942465373274440222033891E16L,
 -5.529326354780295177243773419090123407550E14L,
  3.013431465522152289279088265336861140391E12L,
 -7.959436160727126750732203098982718347785E9L,
  8.230845651379566339707130644134372793322E6L,
};
#define NY0_2D 7
static long double Y0_2D[NY0_2D + 1] = {
  1.438972634353286978700329883122253752192E20L,
  1.856409101981569254247700169486907405500E18L,
  1.219693352678218589553725579802986255614E16L,
  5.389428943282838648918475915779958097958E13L,
  1.774125762108874864433872173544743051653E11L,
  4.522104832545149534808218252434693007036E8L,
  8.872187401232943927082914504125234454930E5L,
  1.251945613186787532055610876304669413955E3L,
 /* 1.000000000000000000000000000000000000000E0 */
};


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

long double
 __ieee754_y0l(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);
      return 1.0L / zero;
    }
  xx = fabsl (x);
  if (xx <= 2.0L)
    {
      /* 0 <= x <= 2 */
      z = xx * xx;
      p = neval (z, Y0_2N, NY0_2N) / deval (z, Y0_2D, NY0_2D);
      p = TWOOPI * logl(x) * __ieee754_j0l(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 * xinv * q;
  q = q - 0.125L * xinv;
  /* X = x - pi/4
     cos(X) = cos(x) cos(pi/4) + sin(x) sin(pi/4)
     = 1/sqrt(2) * (cos(x) + sin(x))
     sin(X) = sin(x) cos(pi/4) - cos(x) sin(pi/4)
     = 1/sqrt(2) * (sin(x) - cos(x))
     sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
     cf. Fdlibm.  */
  c = cosl (x);
  s = sinl (x);
  ss = s - c;
  cc = s + c;
  z = -cosl (x + x);
  if ((s * c) < 0)
    cc = z / ss;
  else
    ss = z / cc;
  z = ONEOSQPI * (p * ss + q * cc) / sqrtl (x);
  return z;
}