Back to index

glibc  2.9
Classes | Defines | Functions
math_private.h File Reference
#include <endian.h>
#include <sys/types.h>
#include <math_ldbl.h>

Go to the source code of this file.

Classes

union  ieee_double_shape_type
union  ieee_double_shape_type
union  ieee_float_shape_type
struct  ieee_double_shape_type.parts
struct  ieee_double_shape_type.parts

Defines

#define EXTRACT_WORDS(ix0, ix1, d)
#define GET_HIGH_WORD(i, d)
#define GET_LOW_WORD(i, d)
#define INSERT_WORDS(d, ix0, ix1)
#define SET_HIGH_WORD(d, v)
#define SET_LOW_WORD(d, v)
#define GET_FLOAT_WORD(i, d)
#define SET_FLOAT_WORD(d, i)
#define math_opt_barrier(x)   ({ __typeof (x) __x = x; __asm ("" : "+m" (__x)); __x; })
#define math_force_eval(x)   __asm __volatile ("" : : "m" (x))

Functions

double __ieee754_sqrt (double)
double __ieee754_acos (double)
double __ieee754_acosh (double)
double __ieee754_log (double)
double __ieee754_atanh (double)
double __ieee754_asin (double)
double __ieee754_atan2 (double, double)
double __ieee754_exp (double)
double __ieee754_exp2 (double)
double __ieee754_exp10 (double)
double __ieee754_cosh (double)
double __ieee754_fmod (double, double)
double __ieee754_pow (double, double)
double __ieee754_lgamma_r (double, int *)
double __ieee754_gamma_r (double, int *)
double __ieee754_lgamma (double)
double __ieee754_gamma (double)
double __ieee754_log10 (double)
double __ieee754_log2 (double)
double __ieee754_sinh (double)
double __ieee754_hypot (double, double)
double __ieee754_j0 (double)
double __ieee754_j1 (double)
double __ieee754_y0 (double)
double __ieee754_y1 (double)
double __ieee754_jn (int, double)
double __ieee754_yn (int, double)
double __ieee754_remainder (double, double)
int32_t __ieee754_rem_pio2 (double, double *)
double __ieee754_scalb (double, double)
double __kernel_standard (double, double, int)
double __kernel_sin (double, double, int)
double __kernel_cos (double, double)
double __kernel_tan (double, double, int)
int __kernel_rem_pio2 (double *, double *, int, int, int, const int32_t *)
double __copysign (double x, double __y)
float __ieee754_sqrtf (float)
float __ieee754_acosf (float)
float __ieee754_acoshf (float)
float __ieee754_logf (float)
float __ieee754_atanhf (float)
float __ieee754_asinf (float)
float __ieee754_atan2f (float, float)
float __ieee754_expf (float)
float __ieee754_exp2f (float)
float __ieee754_exp10f (float)
float __ieee754_coshf (float)
float __ieee754_fmodf (float, float)
float __ieee754_powf (float, float)
float __ieee754_lgammaf_r (float, int *)
float __ieee754_gammaf_r (float, int *)
float __ieee754_lgammaf (float)
float __ieee754_gammaf (float)
float __ieee754_log10f (float)
float __ieee754_log2f (float)
float __ieee754_sinhf (float)
float __ieee754_hypotf (float, float)
float __ieee754_j0f (float)
float __ieee754_j1f (float)
float __ieee754_y0f (float)
float __ieee754_y1f (float)
float __ieee754_jnf (int, float)
float __ieee754_ynf (int, float)
float __ieee754_remainderf (float, float)
int32_t __ieee754_rem_pio2f (float, float *)
float __ieee754_scalbf (float, float)
float __kernel_sinf (float, float, int)
float __kernel_cosf (float, float)
float __kernel_tanf (float, float, int)
int __kernel_rem_pio2f (float *, float *, int, int, int, const int32_t *)
float __copysignf (float x, float __y)
long double __ieee754_sqrtl (long double)
long double __ieee754_acosl (long double)
long double __ieee754_acoshl (long double)
long double __ieee754_logl (long double)
long double __ieee754_atanhl (long double)
long double __ieee754_asinl (long double)
long double __ieee754_atan2l (long double, long double)
long double __ieee754_expl (long double)
long double __ieee754_exp2l (long double)
long double __ieee754_exp10l (long double)
long double __ieee754_coshl (long double)
long double __ieee754_fmodl (long double, long double)
long double __ieee754_powl (long double, long double)
long double __ieee754_lgammal_r (long double, int *)
long double __ieee754_gammal_r (long double, int *)
long double __ieee754_lgammal (long double)
long double __ieee754_gammal (long double)
long double __ieee754_log10l (long double)
long double __ieee754_log2l (long double)
long double __ieee754_sinhl (long double)
long double __ieee754_hypotl (long double, long double)
long double __ieee754_j0l (long double)
long double __ieee754_j1l (long double)
long double __ieee754_y0l (long double)
long double __ieee754_y1l (long double)
long double __ieee754_jnl (int, long double)
long double __ieee754_ynl (int, long double)
long double __ieee754_remainderl (long double, long double)
int __ieee754_rem_pio2l (long double, long double *)
long double __ieee754_scalbl (long double, long double)
long double __kernel_sinl (long double, long double, int)
long double __kernel_cosl (long double, long double)
long double __kernel_tanl (long double, long double, int)
void __kernel_sincosl (long double, long double, long double *, long double *, int)
int __kernel_rem_pio2l (long double *, long double *, int, int, int, const int *)
int __finitel (long double)
int __ilogbl (long double)
int __isinfl (long double)
int __isnanl (long double)
long double __atanl (long double)
long double __copysignl (long double, long double)
long double __expm1l (long double)
long double __floorl (long double)
long double __frexpl (long double, int *)
long double __ldexpl (long double, int)
long double __log1pl (long double)
long double __nanl (const char *)
long double __rintl (long double)
long double __scalbnl (long double, int)
long double __sqrtl (long double x)
long double fabsl (long double x)
void __sincosl (long double, long double *, long double *)
long double __logbl (long double x)
long double __significandl (long double x)
double __exp1 (double __x, double __xx, double __error)
double __sin (double __x)
double __cos (double __x)
int __branred (double __x, double *__a, double *__aa)
void __doasin (double __x, double __dx, double __v[])
void __dubsin (double __x, double __dx, double __v[])
void __dubcos (double __x, double __dx, double __v[])
double __halfulp (double __x, double __y)
double __sin32 (double __x, double __res, double __res1)
double __cos32 (double __x, double __res, double __res1)
double __mpsin (double __x, double __dx)
double __mpcos (double __x, double __dx)
double __mpsin1 (double __x)
double __mpcos1 (double __x)
double __slowexp (double __x)
double __slowpow (double __x, double __y, double __z)
void __docos (double __x, double __dx, double __v[])

Class Documentation

union ieee_double_shape_type

Definition at line 38 of file math_private.h.

Class Members
struct ieee_double_shape_type parts
struct ieee_double_shape_type parts
double value
union ieee_double_shape_type

Definition at line 38 of file math_private.h.

Class Members
struct ieee_double_shape_type parts
struct ieee_double_shape_type parts
double value
union ieee_float_shape_type

Definition at line 125 of file math_private.h.

Class Members
float value
u_int32_t word
struct ieee_double_shape_type.parts

Definition at line 41 of file math_private.h.

Class Members
u_int32_t lsw
u_int32_t msw
struct ieee_double_shape_type.parts

Definition at line 55 of file math_private.h.

Class Members
u_int32_t lsw
u_int32_t msw

Define Documentation

#define EXTRACT_WORDS (   ix0,
  ix1,
  d 
)
Value:
do {                                                    \
  ieee_double_shape_type ew_u;                                 \
  ew_u.value = (d);                                     \
  (ix0) = ew_u.parts.msw;                               \
  (ix1) = ew_u.parts.lsw;                               \
} while (0)

Definition at line 66 of file math_private.h.

#define GET_FLOAT_WORD (   i,
  d 
)
Value:
do {                                                    \
  ieee_float_shape_type gf_u;                                  \
  gf_u.value = (d);                                     \
  (i) = gf_u.word;                                      \
} while (0)

Definition at line 133 of file math_private.h.

#define GET_HIGH_WORD (   i,
  d 
)
Value:
do {                                                    \
  ieee_double_shape_type gh_u;                                 \
  gh_u.value = (d);                                     \
  (i) = gh_u.parts.msw;                                        \
} while (0)

Definition at line 76 of file math_private.h.

#define GET_LOW_WORD (   i,
  d 
)
Value:
do {                                                    \
  ieee_double_shape_type gl_u;                                 \
  gl_u.value = (d);                                     \
  (i) = gl_u.parts.lsw;                                        \
} while (0)

Definition at line 85 of file math_private.h.

#define INSERT_WORDS (   d,
  ix0,
  ix1 
)
Value:
do {                                                    \
  ieee_double_shape_type iw_u;                                 \
  iw_u.parts.msw = (ix0);                               \
  iw_u.parts.lsw = (ix1);                               \
  (d) = iw_u.value;                                     \
} while (0)

Definition at line 94 of file math_private.h.

#define math_force_eval (   x)    __asm __volatile ("" : : "m" (x))

Definition at line 338 of file math_private.h.

#define math_opt_barrier (   x)    ({ __typeof (x) __x = x; __asm ("" : "+m" (__x)); __x; })

Definition at line 336 of file math_private.h.

#define SET_FLOAT_WORD (   d,
  i 
)
Value:
do {                                                    \
  ieee_float_shape_type sf_u;                                  \
  sf_u.word = (i);                                      \
  (d) = sf_u.value;                                     \
} while (0)

Definition at line 142 of file math_private.h.

#define SET_HIGH_WORD (   d,
 
)
Value:
do {                                                    \
  ieee_double_shape_type sh_u;                                 \
  sh_u.value = (d);                                     \
  sh_u.parts.msw = (v);                                        \
  (d) = sh_u.value;                                     \
} while (0)

Definition at line 104 of file math_private.h.

#define SET_LOW_WORD (   d,
 
)
Value:
do {                                                    \
  ieee_double_shape_type sl_u;                                 \
  sl_u.value = (d);                                     \
  sl_u.parts.lsw = (v);                                        \
  (d) = sl_u.value;                                     \
} while (0)

Definition at line 114 of file math_private.h.


Function Documentation

long double __atanl ( long  double)

Definition at line 6 of file s_atanl.c.

{
  fputs ("__atanl not implemented\n", stderr);
  __set_errno (ENOSYS);
  return 0.0;
}

Here is the call graph for this function:

Here is the caller graph for this function:

int __branred ( double  __x,
double *  __a,
double *  __aa 
)

Definition at line 48 of file branred.c.

{
  int i,k;
#if 0
  int n;
#endif
  mynumber  u,gor;
#if 0
  mynumber v;
#endif
  double r[6],s,t,sum,b,bb,sum1,sum2,b1,bb1,b2,bb2,x1,x2,t1,t2;

  x*=tm600.x;
  t=x*split;   /* split x to two numbers */
  x1=t-(t-x);
  x2=x-x1;
  sum=0;
  u.x = x1;
  k = (u.i[HIGH_HALF]>>20)&2047;
  k = (k-450)/24;
  if (k<0)
    k=0;
  gor.x = t576.x;
  gor.i[HIGH_HALF] -= ((k*24)<<20);
  for (i=0;i<6;i++)
    { r[i] = x1*toverp[k+i]*gor.x; gor.x *= tm24.x; }
  for (i=0;i<3;i++) {
    s=(r[i]+big.x)-big.x;
    sum+=s;
    r[i]-=s;
  }
  t=0;
  for (i=0;i<6;i++)
    t+=r[5-i];
  bb=(((((r[0]-t)+r[1])+r[2])+r[3])+r[4])+r[5];
  s=(t+big.x)-big.x;
  sum+=s;
  t-=s;
  b=t+bb;
  bb=(t-b)+bb;
  s=(sum+big1.x)-big1.x;
  sum-=s;
  b1=b;
  bb1=bb;
  sum1=sum;
  sum=0;

  u.x = x2;
  k = (u.i[HIGH_HALF]>>20)&2047;
  k = (k-450)/24;
  if (k<0)
    k=0;
  gor.x = t576.x;
  gor.i[HIGH_HALF] -= ((k*24)<<20);
  for (i=0;i<6;i++)
    { r[i] = x2*toverp[k+i]*gor.x; gor.x *= tm24.x; }
  for (i=0;i<3;i++) {
    s=(r[i]+big.x)-big.x;
    sum+=s;
    r[i]-=s;
  }
  t=0;
  for (i=0;i<6;i++)
    t+=r[5-i];
  bb=(((((r[0]-t)+r[1])+r[2])+r[3])+r[4])+r[5];
  s=(t+big.x)-big.x;
 sum+=s;
 t-=s;
 b=t+bb;
 bb=(t-b)+bb;
 s=(sum+big1.x)-big1.x;
 sum-=s;

 b2=b;
 bb2=bb;
 sum2=sum;

 sum=sum1+sum2;
 b=b1+b2;
 bb = (ABS(b1)>ABS(b2))? (b1-b)+b2 : (b2-b)+b1;
 if (b > 0.5)
   {b-=1.0; sum+=1.0;}
 else if (b < -0.5)
   {b+=1.0; sum-=1.0;}
 s=b+(bb+bb1+bb2);
 t=((b-s)+bb)+(bb1+bb2);
 b=s*split;
 t1=b-(b-s);
 t2=s-t1;
 b=s*hp0.x;
 bb=(((t1*mp1.x-b)+t1*mp2.x)+t2*mp1.x)+(t2*mp2.x+s*hp1.x+t*hp0.x);
 s=b+bb;
 t=(b-s)+bb;
 *a=s;
 *aa=t;
 return ((int) sum)&3; /* return quater of unit circle */
}

Here is the caller graph for this function:

double __copysign ( double  x,
double  __y 
)

Definition at line 24 of file s_copysign.c.

{
  __asm ("cpys %1, %2, %0" : "=f" (x) : "f" (y), "f" (x));
  return x;
}

Here is the caller graph for this function:

float __copysignf ( float  x,
float  __y 
)

Definition at line 23 of file s_copysignf.c.

{
  __asm ("cpys %1, %2, %0" : "=f" (x) : "f" (y), "f" (x));
  return x;
}

Here is the caller graph for this function:

long double __copysignl ( long  double,
long  double 
)

Definition at line 32 of file s_copysignl.c.

{
       u_int64_t hx,hy;
       GET_LDOUBLE_MSW64(hx,x);
       GET_LDOUBLE_MSW64(hy,y);
       SET_LDOUBLE_MSW64(x,(hx&0x7fffffffffffffffULL)
                         |(hy&0x8000000000000000ULL));
        return x;
}

Here is the caller graph for this function:

double __cos ( double  __x)

Definition at line 342 of file s_sin.c.

{
  double y,xx,res,t,cor,s,c,sn,ssn,cs,ccs,xn,a,da,db,eps,xn1,xn2;
  mynumber u,v;
  int4 k,m,n;

  u.x = x;
  m = u.i[HIGH_HALF];
  k = 0x7fffffff&m;

  if (k < 0x3e400000 ) return 1.0; /* |x|<2^-27 => cos(x)=1 */

  else if (k < 0x3feb6000 ) {/* 2^-27 < |x| < 0.855469 */
    y=ABS(x);
    u.x = big.x+y;
    y = y-(u.x-big.x);
    xx=y*y;
    s = y + y*xx*(sn3 +xx*sn5);
    c = xx*(cs2 +xx*(cs4 + xx*cs6));
    k=u.i[LOW_HALF]<<2;
    sn=sincos.x[k];
    ssn=sincos.x[k+1];
    cs=sincos.x[k+2];
    ccs=sincos.x[k+3];
    cor=(ccs-s*ssn-cs*c)-sn*s;
    res=cs+cor;
    cor=(cs-res)+cor;
    return (res==res+1.020*cor)? res : cslow2(x);

}    /*   else  if (k < 0x3feb6000)    */

  else if (k <  0x400368fd ) {/* 0.855469  <|x|<2.426265  */;
    y=hp0.x-ABS(x);
    a=y+hp1.x;
    da=(y-a)+hp1.x;
    xx=a*a;
    if (xx < 0.01588) {
      t = (((((s5.x*xx + s4.x)*xx + s3.x)*xx + s2.x)*xx + s1.x)*a - 0.5*da)*xx+da;
      res = a+t;
      cor = (a-res)+t;
      cor = (cor>0)? 1.02*cor+1.0e-31 : 1.02*cor -1.0e-31;
      return (res == res + cor)? res : csloww(a,da,x);
    }
    else  {
      if (a>0) {m=1;t=a;db=da;}
      else {m=0;t=-a;db=-da;}
      u.x=big.x+t;
      y=t-(u.x-big.x);
      xx=y*y;
      s = y + (db+y*xx*(sn3 +xx*sn5));
      c = y*db+xx*(cs2 +xx*(cs4 + xx*cs6));
      k=u.i[LOW_HALF]<<2;
      sn=sincos.x[k];
      ssn=sincos.x[k+1];
      cs=sincos.x[k+2];
      ccs=sincos.x[k+3];
      cor=(ssn+s*ccs-sn*c)+cs*s;
      res=sn+cor;
      cor=(sn-res)+cor;
      cor = (cor>0)? 1.035*cor+1.0e-31 : 1.035*cor-1.0e-31;
      return (res==res+cor)? ((m)?res:-res) : csloww1(a,da,x);
}

}    /*   else  if (k < 0x400368fd)    */


  else if (k < 0x419921FB ) {/* 2.426265<|x|< 105414350 */
    t = (x*hpinv.x + toint.x);
    xn = t - toint.x;
    v.x = t;
    y = (x - xn*mp1.x) - xn*mp2.x;
    n =v.i[LOW_HALF]&3;
    da = xn*mp3.x;
    a=y-da;
    da = (y-a)-da;
    eps = ABS(x)*1.2e-30;

    switch (n) {
    case 1:
    case 3:
      xx = a*a;
      if (n == 1) {a=-a;da=-da;}
      if (xx < 0.01588) {
       t = (((((s5.x*xx + s4.x)*xx + s3.x)*xx + s2.x)*xx + s1.x)*a - 0.5*da)*xx+da;
       res = a+t;
       cor = (a-res)+t;
       cor = (cor>0)? 1.02*cor+eps : 1.02*cor -eps;
       return (res == res + cor)? res : csloww(a,da,x);
      }
      else  {
       if (a>0) {m=1;t=a;db=da;}
       else {m=0;t=-a;db=-da;}
       u.x=big.x+t;
       y=t-(u.x-big.x);
       xx=y*y;
       s = y + (db+y*xx*(sn3 +xx*sn5));
       c = y*db+xx*(cs2 +xx*(cs4 + xx*cs6));
       k=u.i[LOW_HALF]<<2;
       sn=sincos.x[k];
       ssn=sincos.x[k+1];
       cs=sincos.x[k+2];
       ccs=sincos.x[k+3];
       cor=(ssn+s*ccs-sn*c)+cs*s;
       res=sn+cor;
       cor=(sn-res)+cor;
       cor = (cor>0)? 1.035*cor+eps : 1.035*cor-eps;
       return (res==res+cor)? ((m)?res:-res) : csloww1(a,da,x);
      }
      break;

  case 0:
    case 2:
      if (a<0) {a=-a;da=-da;}
      u.x=big.x+a;
      y=a-(u.x-big.x)+da;
      xx=y*y;
      k=u.i[LOW_HALF]<<2;
      sn=sincos.x[k];
      ssn=sincos.x[k+1];
      cs=sincos.x[k+2];
      ccs=sincos.x[k+3];
      s = y + y*xx*(sn3 +xx*sn5);
      c = xx*(cs2 +xx*(cs4 + xx*cs6));
      cor=(ccs-s*ssn-cs*c)-sn*s;
      res=cs+cor;
      cor=(cs-res)+cor;
      cor = (cor>0)? 1.025*cor+eps : 1.025*cor-eps;
      return (res==res+cor)? ((n)?-res:res) : csloww2(a,da,x,n);

           break;

    }

  }    /*   else  if (k <  0x419921FB )    */


  else if (k < 0x42F00000 ) {
    t = (x*hpinv.x + toint.x);
    xn = t - toint.x;
    v.x = t;
    xn1 = (xn+8.0e22)-8.0e22;
    xn2 = xn - xn1;
    y = ((((x - xn1*mp1.x) - xn1*mp2.x)-xn2*mp1.x)-xn2*mp2.x);
    n =v.i[LOW_HALF]&3;
    da = xn1*pp3.x;
    t=y-da;
    da = (y-t)-da;
    da = (da - xn2*pp3.x) -xn*pp4.x;
    a = t+da;
    da = (t-a)+da;
    eps = 1.0e-24;

    switch (n) {
    case 1:
    case 3:
      xx = a*a;
      if (n==1) {a=-a;da=-da;}
      if (xx < 0.01588) {
       t = (((((s5.x*xx + s4.x)*xx + s3.x)*xx + s2.x)*xx + s1.x)*a - 0.5*da)*xx+da;
       res = a+t;
       cor = (a-res)+t;
       cor = (cor>0)? 1.02*cor+eps : 1.02*cor -eps;
       return (res == res + cor)? res : bsloww(a,da,x,n);
      }
      else  {
       if (a>0) {m=1;t=a;db=da;}
       else {m=0;t=-a;db=-da;}
       u.x=big.x+t;
       y=t-(u.x-big.x);
       xx=y*y;
       s = y + (db+y*xx*(sn3 +xx*sn5));
       c = y*db+xx*(cs2 +xx*(cs4 + xx*cs6));
       k=u.i[LOW_HALF]<<2;
       sn=sincos.x[k];
       ssn=sincos.x[k+1];
       cs=sincos.x[k+2];
       ccs=sincos.x[k+3];
       cor=(ssn+s*ccs-sn*c)+cs*s;
       res=sn+cor;
       cor=(sn-res)+cor;
       cor = (cor>0)? 1.035*cor+eps : 1.035*cor-eps;
       return (res==res+cor)? ((m)?res:-res) : bsloww1(a,da,x,n);
      }
      break;

    case 0:
    case 2:
      if (a<0) {a=-a;da=-da;}
      u.x=big.x+a;
      y=a-(u.x-big.x)+da;
      xx=y*y;
      k=u.i[LOW_HALF]<<2;
      sn=sincos.x[k];
      ssn=sincos.x[k+1];
      cs=sincos.x[k+2];
      ccs=sincos.x[k+3];
      s = y + y*xx*(sn3 +xx*sn5);
      c = xx*(cs2 +xx*(cs4 + xx*cs6));
      cor=(ccs-s*ssn-cs*c)-sn*s;
      res=cs+cor;
      cor=(cs-res)+cor;
      cor = (cor>0)? 1.025*cor+eps : 1.025*cor-eps;
      return (res==res+cor)? ((n)?-res:res) : bsloww2(a,da,x,n);
      break;

    }

  }    /*   else  if (k <  0x42F00000 )    */

  else if (k < 0x7ff00000) {/* 281474976710656 <|x| <2^1024 */

    n = __branred(x,&a,&da);
    switch (n) {
    case 1:
      if (a*a < 0.01588) return bsloww(-a,-da,x,n);
      else return bsloww1(-a,-da,x,n);
      break;
              case 3:
                if (a*a < 0.01588) return bsloww(a,da,x,n);
                else return bsloww1(a,da,x,n);
                break;

    case 0:
    case 2:
      return  bsloww2(a,da,x,n);
      break;
    }

  }    /*   else  if (k <  0x7ff00000 )    */




  else return x / x; /* |x| > 2^1024 */
  return 0;

}

Here is the call graph for this function:

Here is the caller graph for this function:

double __cos32 ( double  __x,
double  __res,
double  __res1 
)

Definition at line 161 of file sincos32.c.

                                                  {
  int p;
  mp_no a,b,c;
  p=32;
  __dbl_mp(res,&a,p);
  __dbl_mp(0.5*(res1-res),&b,p);
  __add(&a,&b,&c,p);
  if (x>2.4)
  { __sub(&pi,&c,&a,p);
    __c32(&a,&b,&c,p);
    b.d[0]=-b.d[0];
  }
  else if (x>0.8)
       { __sub(&hp,&c,&a,p);
         __c32(&a,&c,&b,p);
       }
  else __c32(&c,&b,&a,p);     /* b=cos(0.5*(res+res1))  */
  __dbl_mp(x,&c,p);    /* c = x                  */
  __sub(&b,&c,&a,p);
             /* if a>0 return max(res,res1), otherwise return min(res,res1) */
  if (a.d[0]>0)  return (res>res1)?res:res1;
  else  return (res<res1)?res:res1;
}

Here is the caller graph for this function:

void __doasin ( double  __x,
double  __dx,
double  __v[] 
)

Definition at line 41 of file doasin.c.

                                               {

#include "doasin.h"

  static const double
    d5 =  0.22372159090911789889975459505194491E-01,
    d6 =  0.17352764422456822913014975683014622E-01,
    d7 =  0.13964843843786693521653681033981614E-01,
    d8 =  0.11551791438485242609036067259086589E-01,
    d9 =  0.97622386568166960207425666787248914E-02,
    d10 = 0.83638737193775788576092749009744976E-02,
    d11 = 0.79470250400727425881446981833568758E-02;

  double xx,p,pp,u,uu,r,s;
  double hx,tx,hy,ty,tp,tq,tc,tcc;


/* Taylor series for arcsin for Double-Length numbers         */
  xx = x*x+2.0*x*dx;
  p = ((((((d11*xx+d10)*xx+d9)*xx+d8)*xx+d7)*xx+d6)*xx+d5)*xx;
  pp = 0;

  MUL2(x,dx,x,dx,u,uu,tp,hx,tx,hy,ty,tq,tc,tcc);
  ADD2(p,pp,c4.x,cc4.x,p,pp,r,s);
  MUL2(p,pp,u,uu,p,pp,tp,hx,tx,hy,ty,tq,tc,tcc);
  ADD2(p,pp,c3.x,cc3.x,p,pp,r,s);
  MUL2(p,pp,u,uu,p,pp,tp,hx,tx,hy,ty,tq,tc,tcc);
  ADD2(p,pp,c2.x,cc2.x,p,pp,r,s);
  MUL2(p,pp,u,uu,p,pp,tp,hx,tx,hy,ty,tq,tc,tcc);
  ADD2(p,pp,c1.x,cc1.x,p,pp,r,s);
  MUL2(p,pp,u,uu,p,pp,tp,hx,tx,hy,ty,tq,tc,tcc);
  MUL2(p,pp,x,dx,p,pp,tp,hx,tx,hy,ty,tq,tc,tcc);
  ADD2(p,pp,x,dx,p,pp,r,s);
  v[0]=p;
  v[1]=pp; /* arcsin(x+dx)=v[0]+v[1] */
}

Here is the caller graph for this function:

void __docos ( double  __x,
double  __dx,
double  __v[] 
)

Definition at line 164 of file dosincos.c.

                                              {
  double y,yy,p,w[2];
  if (x>0) {y=x; yy=dx;}
     else {y=-x; yy=-dx;}
  if (y<0.5*hp0.x)                                 /*  y< PI/4    */
           {__dubcos(y,yy,w); v[0]=w[0]; v[1]=w[1];}
     else if (y<1.5*hp0.x) {                       /* y< 3/4 * PI */
       p=hp0.x-y;  /* p = PI/2 - y */
       yy=hp1.x-yy;
       y=p+yy;
       yy=(p-y)+yy;
       if (y>0) {__dubsin(y,yy,w); v[0]=w[0]; v[1]=w[1];}
                                       /* cos(x) = sin ( 90 -  x ) */
         else {__dubsin(-y,-yy,w); v[0]=-w[0]; v[1]=-w[1];
        }
     }
  else { /* y>= 3/4 * PI */
    p=2.0*hp0.x-y;    /* p = PI- y */
    yy=2.0*hp1.x-yy;
    y=p+yy;
    yy=(p-y)+yy;
    __dubcos(y,yy,w);
    v[0]=-w[0];
    v[1]=-w[1];
  }
}

Here is the caller graph for this function:

void __dubcos ( double  __x,
double  __dx,
double  __v[] 
)

Definition at line 101 of file dosincos.c.

                                               {
  double r,s,p,hx,tx,hy,ty,q,c,cc,d,dd,d2,dd2,e,ee,
    sn,ssn,cs,ccs,ds,dss,dc,dcc;
#if 0
  double xx,y,yy,z,zz;
#endif
  mynumber u;
  int4 k;
  u.x=x+big.x;
  k = u.i[LOW_HALF]<<2;
  x=x-(u.x-big.x);
  d=x+dx;
  dd=(x-d)+dx;  /* cos(x+dx)=cos(Xi+t)=cos(Xi)cos(t) - sin(Xi)sin(t) */
  MUL2(d,dd,d,dd,d2,dd2,p,hx,tx,hy,ty,q,c,cc);
  sn=sincos.x[k];     /*                                  */
  ssn=sincos.x[k+1];  /*      sin(Xi) and cos(Xi)         */
  cs=sincos.x[k+2];   /*                                  */
  ccs=sincos.x[k+3];  /*                                  */
  MUL2(d2,dd2,s7.x,ss7.x,ds,dss,p,hx,tx,hy,ty,q,c,cc);
  ADD2(ds,dss,s5.x,ss5.x,ds,dss,r,s);
  MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc);
  ADD2(ds,dss,s3.x,ss3.x,ds,dss,r,s);
  MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc);
  MUL2(d,dd,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc);
  ADD2(ds,dss,d,dd,ds,dss,r,s);

  MUL2(d2,dd2,c8.x,cc8.x,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
  ADD2(dc,dcc,c6.x,cc6.x,dc,dcc,r,s);
  MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
  ADD2(dc,dcc,c4.x,cc4.x,dc,dcc,r,s);
  MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
  ADD2(dc,dcc,c2.x,cc2.x,dc,dcc,r,s);
  MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);

  MUL2(cs,ccs,ds,dss,e,ee,p,hx,tx,hy,ty,q,c,cc);
  MUL2(dc,dcc,sn,ssn,dc,dcc,p,hx,tx,hy,ty,q,c,cc);

  MUL2(d2,dd2,s7.x,ss7.x,ds,dss,p,hx,tx,hy,ty,q,c,cc);
  ADD2(ds,dss,s5.x,ss5.x,ds,dss,r,s);
  MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc);
  ADD2(ds,dss,s3.x,ss3.x,ds,dss,r,s);
  MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc);
  MUL2(d,dd,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc);
  ADD2(ds,dss,d,dd,ds,dss,r,s);
  MUL2(d2,dd2,c8.x,cc8.x,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
  ADD2(dc,dcc,c6.x,cc6.x,dc,dcc,r,s);
  MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
  ADD2(dc,dcc,c4.x,cc4.x,dc,dcc,r,s);
  MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
  ADD2(dc,dcc,c2.x,cc2.x,dc,dcc,r,s);
  MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
  MUL2(sn,ssn,ds,dss,e,ee,p,hx,tx,hy,ty,q,c,cc);
  MUL2(dc,dcc,cs,ccs,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
  ADD2(e,ee,dc,dcc,e,ee,r,s);
  SUB2(cs,ccs,e,ee,e,ee,r,s);

  v[0]=e;
  v[1]=ee;
}

Here is the caller graph for this function:

void __dubsin ( double  __x,
double  __dx,
double  __v[] 
)

Definition at line 50 of file dosincos.c.

                                               {
  double r,s,p,hx,tx,hy,ty,q,c,cc,d,dd,d2,dd2,e,ee,
    sn,ssn,cs,ccs,ds,dss,dc,dcc;
#if 0
  double xx,y,yy,z,zz;
#endif
  mynumber u;
  int4 k;

  u.x=x+big.x;
  k = u.i[LOW_HALF]<<2;
  x=x-(u.x-big.x);
  d=x+dx;
  dd=(x-d)+dx;
         /* sin(x+dx)=sin(Xi+t)=sin(Xi)*cos(t) + cos(Xi)sin(t) where t ->0 */
  MUL2(d,dd,d,dd,d2,dd2,p,hx,tx,hy,ty,q,c,cc);
  sn=sincos.x[k];     /*                                  */
  ssn=sincos.x[k+1];  /*      sin(Xi) and cos(Xi)         */
  cs=sincos.x[k+2];   /*                                  */
  ccs=sincos.x[k+3];  /*                                  */
  MUL2(d2,dd2,s7.x,ss7.x,ds,dss,p,hx,tx,hy,ty,q,c,cc);  /* Taylor    */
  ADD2(ds,dss,s5.x,ss5.x,ds,dss,r,s);
  MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc);      /* series    */
  ADD2(ds,dss,s3.x,ss3.x,ds,dss,r,s);
  MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc);      /* for sin   */
  MUL2(d,dd,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc);
  ADD2(ds,dss,d,dd,ds,dss,r,s);                         /* ds=sin(t) */

  MUL2(d2,dd2,c8.x,cc8.x,dc,dcc,p,hx,tx,hy,ty,q,c,cc); ;/* Taylor    */
  ADD2(dc,dcc,c6.x,cc6.x,dc,dcc,r,s);
  MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);      /* series    */
  ADD2(dc,dcc,c4.x,cc4.x,dc,dcc,r,s);
  MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);      /* for cos   */
  ADD2(dc,dcc,c2.x,cc2.x,dc,dcc,r,s);
  MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);      /* dc=cos(t) */

  MUL2(cs,ccs,ds,dss,e,ee,p,hx,tx,hy,ty,q,c,cc);
  MUL2(dc,dcc,sn,ssn,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
  SUB2(e,ee,dc,dcc,e,ee,r,s);
  ADD2(e,ee,sn,ssn,e,ee,r,s);                    /* e+ee=sin(x+dx) */

  v[0]=e;
  v[1]=ee;
}

Here is the caller graph for this function:

double __exp1 ( double  __x,
double  __xx,
double  __error 
)

Definition at line 156 of file e_exp.c.

                                                 {
  double bexp, t, eps, del, base, y, al, bet, res, rem, cor;
  mynumber junk1, junk2, binexp  = {{0,0}};
#if 0
  int4 k;
#endif
  int4 i,j,m,n,ex;

  junk1.x = x;
  m = junk1.i[HIGH_HALF];
  n = m&hugeint;                 /* no sign */

  if (n > smallint && n < bigint) {
    y = x*log2e.x + three51.x;
    bexp = y - three51.x;      /*  multiply the result by 2**bexp        */

    junk1.x = y;

    eps = bexp*ln_two2.x;      /* x = bexp*ln(2) + t - eps               */
    t = x - bexp*ln_two1.x;

    y = t + three33.x;
    base = y - three33.x;      /* t rounded to a multiple of 2**-18      */
    junk2.x = y;
    del = (t - base) + (xx-eps);    /*  x = bexp*ln(2) + base + del      */
    eps = del + del*del*(p3.x*del + p2.x);

    binexp.i[HIGH_HALF] =(junk1.i[LOW_HALF]+1023)<<20;

    i = ((junk2.i[LOW_HALF]>>8)&0xfffffffe)+356;
    j = (junk2.i[LOW_HALF]&511)<<1;

    al = coar.x[i]*fine.x[j];
    bet =(coar.x[i]*fine.x[j+1] + coar.x[i+1]*fine.x[j]) + coar.x[i+1]*fine.x[j+1];

    rem=(bet + bet*eps)+al*eps;
    res = al + rem;
    cor = (al - res) + rem;
    if  (res == (res+cor*(1.0+error+err_1))) return res*binexp.x;
    else return -10.0;
  }

  if (n <= smallint) return 1.0; /*  if x->0 e^x=1 */

  if (n >= badint) {
    if (n > infint) return(zero/zero);    /* x is NaN,  return invalid */
    if (n < infint) return ( (x>0) ? (hhuge*hhuge) : (tiny*tiny) );
    /* x is finite,  cause either overflow or underflow  */
    if (junk1.i[LOW_HALF] != 0)  return (zero/zero);        /*  x is NaN  */
    return ((x>0)?inf.x:zero );   /* |x| = inf;  return either inf or 0 */
  }

  y = x*log2e.x + three51.x;
  bexp = y - three51.x;
  junk1.x = y;
  eps = bexp*ln_two2.x;
  t = x - bexp*ln_two1.x;
  y = t + three33.x;
  base = y - three33.x;
  junk2.x = y;
  del = (t - base) + (xx-eps);
  eps = del + del*del*(p3.x*del + p2.x);
  i = ((junk2.i[LOW_HALF]>>8)&0xfffffffe)+356;
  j = (junk2.i[LOW_HALF]&511)<<1;
  al = coar.x[i]*fine.x[j];
  bet =(coar.x[i]*fine.x[j+1] + coar.x[i+1]*fine.x[j]) + coar.x[i+1]*fine.x[j+1];
  rem=(bet + bet*eps)+al*eps;
  res = al + rem;
  cor = (al - res) + rem;
  if (m>>31) {
    ex=junk1.i[LOW_HALF];
    if (res < 1.0) {res+=res; cor+=cor; ex-=1;}
    if (ex >=-1022) {
      binexp.i[HIGH_HALF] = (1023+ex)<<20;
      if  (res == (res+cor*(1.0+error+err_1))) return res*binexp.x;
      else return -10.0;
    }
    ex = -(1022+ex);
    binexp.i[HIGH_HALF] = (1023-ex)<<20;
    res*=binexp.x;
    cor*=binexp.x;
    eps=1.00000000001+(error+err_1)*binexp.x;
    t=1.0+res;
    y = ((1.0-t)+res)+cor;
    res=t+y;
    cor = (t-res)+y;
    if (res == (res + eps*cor))
      {binexp.i[HIGH_HALF] = 0x00100000; return (res-1.0)*binexp.x;}
    else return -10.0;
  }
  else {
    binexp.i[HIGH_HALF] =(junk1.i[LOW_HALF]+767)<<20;
    if  (res == (res+cor*(1.0+error+err_1)))
      return res*binexp.x*t256.x;
    else return -10.0;
  }
}

Here is the caller graph for this function:

long double __expm1l ( long  double)

Definition at line 6 of file s_expm1l.c.

{
  fputs ("__expm1l not implemented\n", stderr);
  __set_errno (ENOSYS);
  return 0.0;
}

Here is the call graph for this function:

Here is the caller graph for this function:

int __finitel ( long  double)

Definition at line 31 of file s_finitel.c.

{
       int64_t hx;
       GET_LDOUBLE_MSW64(hx,x);
       return (int)((u_int64_t)((hx&0x7fffffffffffffffLL)
                             -0x7fff000000000000LL)>>63);
}
long double __floorl ( long  double)

Definition at line 41 of file s_floorl.c.

{
       int64_t i0,i1,j0;
       u_int64_t i,j;
       GET_LDOUBLE_WORDS64(i0,i1,x);
       j0 = ((i0>>48)&0x7fff)-0x3fff;
       if(j0<48) {
           if(j0<0) {       /* raise inexact if x != 0 */
              if(huge+x>0.0) {/* return 0*sign(x) if |x|<1 */
                  if(i0>=0) {i0=i1=0;}
                  else if(((i0&0x7fffffffffffffffLL)|i1)!=0)
                     { i0=0xbfff000000000000ULL;i1=0;}
              }
           } else {
              i = (0x0000ffffffffffffULL)>>j0;
              if(((i0&i)|i1)==0) return x; /* x is integral */
              if(huge+x>0.0) {     /* raise inexact flag */
                  if(i0<0) i0 += (0x0001000000000000LL)>>j0;
                  i0 &= (~i); i1=0;
              }
           }
       } else if (j0>111) {
           if(j0==0x4000) return x+x;     /* inf or NaN */
           else return x;          /* x is integral */
       } else {
           i = -1ULL>>(j0-48);
           if((i1&i)==0) return x; /* x is integral */
           if(huge+x>0.0) {               /* raise inexact flag */
              if(i0<0) {
                  if(j0==48) i0+=1;
                  else {
                     j = i1+(1LL<<(112-j0));
                     if(j<i1) i0 +=1 ;    /* got a carry */
                     i1=j;
                  }
              }
              i1 &= (~i);
           }
       }
       SET_LDOUBLE_WORDS64(x,i0,i1);
       return x;
}

Here is the call graph for this function:

Here is the caller graph for this function:

long double __frexpl ( long  double,
int  
)

Definition at line 22 of file s_frexpl.c.

{
  long double mantissa, exponent;
  int iexponent;
  unsigned long fpsr;

  __asm ("ftst%.x %1\n"
        "fmove%.l %/fpsr, %0"
        : "=dm" (fpsr) : "f" (value));
  if (fpsr & (7 << 24))
    {
      /* Not finite or zero.  */
      *expptr = 0;
      return value;
    }
  __asm ("fgetexp%.x %1, %0" : "=f" (exponent) : "f" (value));
  iexponent = (int) exponent + 1;
  *expptr = iexponent;
  /* Unnormalized numbers must be handled specially, otherwise fscale
     results in overflow.  */
  if (iexponent <= -16384)
    {
      value *= 0x1p16383L;
      iexponent += 16383;
    }
  else if (iexponent >= 16384)
    {
      value *= 0x1p-16383L;
      iexponent -= 16383;
    }

  __asm ("fscale%.l %2, %0"
        : "=f" (mantissa)
        : "0" (value), "dmi" (-iexponent));
  return mantissa;
}

Here is the caller graph for this function:

double __halfulp ( double  __x,
double  __y 
)

Definition at line 52 of file halfulp.c.

{
  mynumber v;
  double z,u,uu,j1,j2,j3,j4,j5;
  int4 k,l,m,n;
  if (y <= 0) {               /*if power is negative or zero */
    v.x = y;
    if (v.i[LOW_HALF] != 0) return -10.0;
    v.x = x;
    if (v.i[LOW_HALF] != 0) return -10.0;
    if ((v.i[HIGH_HALF]&0x000fffff) != 0) return -10;   /* if x =2 ^ n */
    k = ((v.i[HIGH_HALF]&0x7fffffff)>>20)-1023;         /* find this n */
    z = (double) k;
    return (z*y == -1075.0)?0: -10.0;
  }
                              /* if y > 0  */
  v.x = y;
    if (v.i[LOW_HALF] != 0) return -10.0;

  v.x=x;
                              /*  case where x = 2**n for some integer n */
  if (((v.i[HIGH_HALF]&0x000fffff)|v.i[LOW_HALF]) == 0) {
    k=(v.i[HIGH_HALF]>>20)-1023;
    return (((double) k)*y == -1075.0)?0:-10.0;
  }

  v.x = y;
  k = v.i[HIGH_HALF];
  m = k<<12;
  l = 0;
  while (m)
    {m = m<<1; l++; }
  n = (k&0x000fffff)|0x00100000;
  n = n>>(20-l);                       /*   n is the odd integer of y    */
  k = ((k>>20) -1023)-l;               /*   y = n*2**k                   */
  if (k>5) return -10.0;
  if (k>0) for (;k>0;k--) n *= 2;
  if (n > 34) return -10.0;
  k = -k;
  if (k>5) return -10.0;

                            /*   now treat x        */
  while (k>0) {
    z = __ieee754_sqrt(x);
    EMULV(z,z,u,uu,j1,j2,j3,j4,j5);
    if (((u-x)+uu) != 0) break;
    x = z;
    k--;
 }
  if (k) return -10.0;

  /* it is impossible that n == 2,  so the mantissa of x must be short  */

  v.x = x;
  if (v.i[LOW_HALF]) return -10.0;
  k = v.i[HIGH_HALF];
  m = k<<12;
  l = 0;
  while (m) {m = m<<1; l++; }
  m = (k&0x000fffff)|0x00100000;
  m = m>>(20-l);                       /*   m is the odd integer of x    */

            /*   now check whether the length of m**n is at most 54 bits */

  if  (m > tab54[n-3]) return -10.0;

             /* yes, it is - now compute x**n by simple multiplications  */

  u = x;
  for (k=1;k<n;k++) u = u*x;
  return u;
}

Here is the caller graph for this function:

double __ieee754_acos ( double  )

Definition at line 334 of file e_asin.c.

{
  double x1,x2,xx,s1,s2,res1,p,t,res,r,cor,cc,y,c,z,w[2],eps;
#if 0
  double fc;
#endif
  mynumber u,v;
  int4 k,m,n;
#if 0
  int4 nn;
#endif
  u.x = x;
  m = u.i[HIGH_HALF];
  k = 0x7fffffff&m;
  /*-------------------  |x|<2.77556*10^-17 ----------------------*/
  if (k < 0x3c880000) return hp0.x;

  /*-----------------  2.77556*10^-17 <= |x| < 2^-3 --------------*/
  else
  if (k < 0x3fc00000) {
    x2 = x*x;
    t = (((((f6*x2 + f5)*x2 + f4)*x2 + f3)*x2 + f2)*x2 + f1)*(x2*x);
    r=hp0.x-x;
    cor=(((hp0.x-r)-x)+hp1.x)-t;
    res = r+cor;
    cor = (r-res)+cor;
    if (res == res+1.004*cor) return res;
    else {
      x1 = x+big;
      xx = x*x;
      x1 -= big;
      x2 = x - x1;
      p = x1*x1*x1;
      s1 = a1.x*p;
      s2 = ((((((c7*xx + c6)*xx + c5)*xx + c4)*xx + c3)*xx + c2)*xx*xx*x +
           ((a1.x+a2.x)*x2*x2+ 0.5*x1*x)*x2) + a2.x*p;
      res1 = x+s1;
      s2 = ((x-res1)+s1)+s2;
      r=hp0.x-res1;
      cor=(((hp0.x-r)-res1)+hp1.x)-s2;
      res = r+cor;
      cor = (r-res)+cor;
      if (res == res+1.00004*cor) return res;
      else {
       __doasin(x,0,w);
       r=hp0.x-w[0];
       cor=((hp0.x-r)-w[0])+(hp1.x-w[1]);
       res=r+cor;
       cor=(r-res)+cor;
       if (res ==(res +1.00000001*cor)) return res;
       else {
         res1=res+1.1*cor;
         return __cos32(x,res,res1);
       }
      }
    }
  }    /*   else  if (k < 0x3fc00000)    */
  /*----------------------  0.125 <= |x| < 0.5 --------------------*/
  else
  if (k < 0x3fe00000) {
    if (k<0x3fd00000) n = 11*((k&0x000fffff)>>15);
    else n = 11*((k&0x000fffff)>>14)+352;
    if (m>0) xx = x - asncs.x[n];
    else xx = -x - asncs.x[n];
    t = asncs.x[n+1]*xx;
    p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
                   xx*(asncs.x[n+5]+xx*asncs.x[n+6]))))+asncs.x[n+7];
    t+=p;
    y = (m>0)?(hp0.x-asncs.x[n+8]):(hp0.x+asncs.x[n+8]);
    t = (m>0)?(hp1.x-t):(hp1.x+t);
    res = y+t;
    if (res == res+1.02*((y-res)+t)) return res;
    else {
      r=asncs.x[n+8]+xx*asncs.x[n+9];
      t=((asncs.x[n+8]-r)+xx*asncs.x[n+9])+(p+xx*asncs.x[n+10]);
      if (m>0)
       {p = hp0.x-r; t = (((hp0.x-p)-r)-t)+hp1.x; }
      else
       {p = hp0.x+r; t = ((hp0.x-p)+r)+(hp1.x+t); }
      res = p+t;
      cor = (p-res)+t;
      if (res == (res+1.0002*cor)) return res;
      else {
       res1=res+1.1*cor;
       z=0.5*(res1-res);
       __docos(res,z,w);
       z=(w[0]-x)+w[1];
       if (z>1.0e-27) return max(res,res1);
       else if (z<-1.0e-27) return min(res,res1);
       else return __cos32(x,res,res1);
      }
    }
  }    /*   else  if (k < 0x3fe00000)    */

  /*--------------------------- 0.5 <= |x| < 0.75 ---------------------*/
  else
  if (k < 0x3fe80000) {
    n = 1056+((k&0x000fe000)>>11)*3;
    if (m>0) {xx = x - asncs.x[n]; eps=1.04; }
    else {xx = -x - asncs.x[n]; eps=1.02; }
    t = asncs.x[n+1]*xx;
    p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
                   xx*(asncs.x[n+5]+xx*(asncs.x[n+6]+
                   xx*asncs.x[n+7])))))+asncs.x[n+8];
    t+=p;
   y = (m>0)?(hp0.x-asncs.x[n+9]):(hp0.x+asncs.x[n+9]);
   t = (m>0)?(hp1.x-t):(hp1.x+t);
   res = y+t;
   if (res == res+eps*((y-res)+t)) return res;
   else {
     r=asncs.x[n+9]+xx*asncs.x[n+10];
     t=((asncs.x[n+9]-r)+xx*asncs.x[n+10])+(p+xx*asncs.x[n+11]);
     if (m>0) {p = hp0.x-r; t = (((hp0.x-p)-r)-t)+hp1.x; eps=1.0004; }
     else   {p = hp0.x+r; t = ((hp0.x-p)+r)+(hp1.x+t); eps=1.0002; }
     res = p+t;
     cor = (p-res)+t;
     if (res == (res+eps*cor)) return res;
     else {
       res1=res+1.1*cor;
       z=0.5*(res1-res);
       __docos(res,z,w);
       z=(w[0]-x)+w[1];
       if (z>1.0e-27) return max(res,res1);
       else if (z<-1.0e-27) return min(res,res1);
       else return __cos32(x,res,res1);
     }
   }
  }    /*   else  if (k < 0x3fe80000)    */

/*------------------------- 0.75 <= |x| < 0.921875 -------------*/
  else
  if (k < 0x3fed8000) {
    n = 992+((k&0x000fe000)>>13)*13;
    if (m>0) {xx = x - asncs.x[n]; eps = 1.04; }
    else {xx = -x - asncs.x[n]; eps = 1.01; }
    t = asncs.x[n+1]*xx;
    p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
                      xx*(asncs.x[n+5]+xx*(asncs.x[n+6]+xx*(asncs.x[n+7]+
                      xx*asncs.x[n+8]))))))+asncs.x[n+9];
    t+=p;
    y = (m>0)?(hp0.x-asncs.x[n+10]):(hp0.x+asncs.x[n+10]);
    t = (m>0)?(hp1.x-t):(hp1.x+t);
    res = y+t;
    if (res == res+eps*((y-res)+t)) return res;
    else {
      r=asncs.x[n+10]+xx*asncs.x[n+11];
      t=((asncs.x[n+10]-r)+xx*asncs.x[n+11])+(p+xx*asncs.x[n+12]);
      if (m>0) {p = hp0.x-r; t = (((hp0.x-p)-r)-t)+hp1.x; eps=1.0032; }
      else   {p = hp0.x+r; t = ((hp0.x-p)+r)+(hp1.x+t); eps=1.0008; }
      res = p+t;
      cor = (p-res)+t;
      if (res == (res+eps*cor)) return res;
      else {
       res1=res+1.1*cor;
       z=0.5*(res1-res);
       __docos(res,z,w);
       z=(w[0]-x)+w[1];
       if (z>1.0e-27) return max(res,res1);
       else if (z<-1.0e-27) return min(res,res1);
       else return __cos32(x,res,res1);
      }
    }
  }    /*   else  if (k < 0x3fed8000)    */

/*-------------------0.921875 <= |x| < 0.953125 ------------------*/
  else
  if (k < 0x3fee8000) {
    n = 884+((k&0x000fe000)>>13)*14;
    if (m>0) {xx = x - asncs.x[n]; eps=1.04; }
    else {xx = -x - asncs.x[n]; eps =1.005; }
    t = asncs.x[n+1]*xx;
    p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
                   xx*(asncs.x[n+5]+xx*(asncs.x[n+6]
                 +xx*(asncs.x[n+7]+xx*(asncs.x[n+8]+
                   xx*asncs.x[n+9])))))))+asncs.x[n+10];
    t+=p;
    y = (m>0)?(hp0.x-asncs.x[n+11]):(hp0.x+asncs.x[n+11]);
    t = (m>0)?(hp1.x-t):(hp1.x+t);
    res = y+t;
    if (res == res+eps*((y-res)+t)) return res;
    else {
      r=asncs.x[n+11]+xx*asncs.x[n+12];
      t=((asncs.x[n+11]-r)+xx*asncs.x[n+12])+(p+xx*asncs.x[n+13]);
      if (m>0) {p = hp0.x-r; t = (((hp0.x-p)-r)-t)+hp1.x; eps=1.0030; }
      else   {p = hp0.x+r; t = ((hp0.x-p)+r)+(hp1.x+t); eps=1.0005; }
      res = p+t;
      cor = (p-res)+t;
      if (res == (res+eps*cor)) return res;
      else {
       res1=res+1.1*cor;
       z=0.5*(res1-res);
       __docos(res,z,w);
       z=(w[0]-x)+w[1];
       if (z>1.0e-27) return max(res,res1);
       else if (z<-1.0e-27) return min(res,res1);
       else return __cos32(x,res,res1);
      }
    }
  }    /*   else  if (k < 0x3fee8000)    */

  /*--------------------0.953125 <= |x| < 0.96875 ----------------*/
  else
  if (k < 0x3fef0000) {
    n = 768+((k&0x000fe000)>>13)*15;
    if (m>0) {xx = x - asncs.x[n]; eps=1.04; }
    else {xx = -x - asncs.x[n]; eps=1.005;}
    t = asncs.x[n+1]*xx;
    p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
            xx*(asncs.x[n+5]+xx*(asncs.x[n+6]
           +xx*(asncs.x[n+7]+xx*(asncs.x[n+8]+xx*(asncs.x[n+9]+
            xx*asncs.x[n+10]))))))))+asncs.x[n+11];
    t+=p;
    y = (m>0)?(hp0.x-asncs.x[n+12]):(hp0.x+asncs.x[n+12]);
   t = (m>0)?(hp1.x-t):(hp1.x+t);
   res = y+t;
   if (res == res+eps*((y-res)+t)) return res;
   else {
     r=asncs.x[n+12]+xx*asncs.x[n+13];
     t=((asncs.x[n+12]-r)+xx*asncs.x[n+13])+(p+xx*asncs.x[n+14]);
     if (m>0) {p = hp0.x-r; t = (((hp0.x-p)-r)-t)+hp1.x; eps=1.0030; }
     else   {p = hp0.x+r; t = ((hp0.x-p)+r)+(hp1.x+t); eps=1.0005; }
     res = p+t;
     cor = (p-res)+t;
     if (res == (res+eps*cor)) return res;
     else {
       res1=res+1.1*cor;
       z=0.5*(res1-res);
       __docos(res,z,w);
       z=(w[0]-x)+w[1];
       if (z>1.0e-27) return max(res,res1);
       else if (z<-1.0e-27) return min(res,res1);
       else return __cos32(x,res,res1);
     }
   }
  }    /*   else  if (k < 0x3fef0000)    */
  /*-----------------0.96875 <= |x| < 1 ---------------------------*/

  else
  if (k<0x3ff00000)  {
    z = 0.5*((m>0)?(1.0-x):(1.0+x));
    v.x=z;
    k=v.i[HIGH_HALF];
    t=inroot[(k&0x001fffff)>>14]*powtwo[511-(k>>21)];
    r=1.0-t*t*z;
    t = t*(rt0+r*(rt1+r*(rt2+r*rt3)));
    c=t*z;
    t=c*(1.5-0.5*t*c);
    y = (t27*c+c)-t27*c;
    cc = (z-y*y)/(t+y);
    p=(((((f6*z+f5)*z+f4)*z+f3)*z+f2)*z+f1)*z;
    if (m<0) {
      cor = (hp1.x - cc)-(y+cc)*p;
      res1 = hp0.x - y;
      res =res1 + cor;
      if (res == res+1.002*((res1-res)+cor)) return (res+res);
      else {
       c=y+cc;
       cc=(y-c)+cc;
       __doasin(c,cc,w);
       res1=hp0.x-w[0];
       cor=((hp0.x-res1)-w[0])+(hp1.x-w[1]);
       res = res1+cor;
       cor = (res1-res)+cor;
       if (res==(res+1.000001*cor)) return (res+res);
       else {
         res=res+res;
         res1=res+1.2*cor;
         return __cos32(x,res,res1);
       }
      }
    }
    else {
      cor = cc+p*(y+cc);
      res = y + cor;
      if (res == res+1.03*((y-res)+cor)) return (res+res);
      else {
       c=y+cc;
       cc=(y-c)+cc;
       __doasin(c,cc,w);
       res = w[0];
       cor=w[1];
       if (res==(res+1.000001*cor)) return (res+res);
       else {
         res=res+res;
         res1=res+1.2*cor;
         return __cos32(x,res,res1);
       }
      }
    }
  }    /*   else  if (k < 0x3ff00000)    */

  /*---------------------------- |x|>=1 -----------------------*/
  else
  if (k==0x3ff00000 && u.i[LOW_HALF]==0) return (m>0)?0:2.0*hp0.x;
  else
  if (k>0x7ff00000 || (k == 0x7ff00000 && u.i[LOW_HALF] != 0)) return x;
  else {
    u.i[HIGH_HALF]=0x7ff00000;
    v.i[HIGH_HALF]=0x7ff00000;
    u.i[LOW_HALF]=0;
    v.i[LOW_HALF]=0;
    return u.x/v.x;
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

float __ieee754_acosf ( float  )

Definition at line 46 of file e_acosf.c.

{
       float z,p,q,r,w,s,c,df;
       int32_t hx,ix;
       GET_FLOAT_WORD(hx,x);
       ix = hx&0x7fffffff;
       if(ix==0x3f800000) {        /* |x|==1 */
           if(hx>0) return 0.0;    /* acos(1) = 0  */
           else return pi+(float)2.0*pio2_lo;    /* acos(-1)= pi */
       } else if(ix>0x3f800000) {  /* |x| >= 1 */
           return (x-x)/(x-x);            /* acos(|x|>1) is NaN */
       }
       if(ix<0x3f000000) {  /* |x| < 0.5 */
           if(ix<=0x23000000) return pio2_hi+pio2_lo;/*if|x|<2**-57*/
           z = x*x;
           p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
           q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
           r = p/q;
           return pio2_hi - (x - (pio2_lo-x*r));
       } else  if (hx<0) {         /* x < -0.5 */
           z = (one+x)*(float)0.5;
           p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
           q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
           s = __ieee754_sqrtf(z);
           r = p/q;
           w = r*s-pio2_lo;
           return pi - (float)2.0*(s+w);
       } else {                    /* x > 0.5 */
           int32_t idf;
           z = (one-x)*(float)0.5;
           s = __ieee754_sqrtf(z);
           df = s;
           GET_FLOAT_WORD(idf,df);
           SET_FLOAT_WORD(df,idf&0xfffff000);
           c  = (z-df*df)/(s+df);
           p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
           q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
           r = p/q;
           w = r*s+c;
           return (float)2.0*(df+w);
       }
}

Here is the call graph for this function:

Here is the caller graph for this function:

double __ieee754_acosh ( double  )

Definition at line 45 of file e_acosh.c.

{      
       double t;
       int32_t hx;
       u_int32_t lx;
       EXTRACT_WORDS(hx,lx,x);
       if(hx<0x3ff00000) {         /* x < 1 */
           return (x-x)/(x-x);
       } else if(hx >=0x41b00000) {       /* x > 2**28 */
           if(hx >=0x7ff00000) {   /* x is inf of NaN */
               return x+x;
           } else 
              return __ieee754_log(x)+ln2;       /* acosh(huge)=log(2x) */
       } else if(((hx-0x3ff00000)|lx)==0) {
           return 0.0;                    /* acosh(1) = 0 */
       } else if (hx > 0x40000000) {      /* 2**28 > x > 2 */
           t=x*x;
           return __ieee754_log(2.0*x-one/(x+__ieee754_sqrt(t-one)));
       } else {                    /* 1<x<2 */
           t = x-one;
           return __log1p(t+__sqrt(2.0*t+t*t));
       }
}

Here is the call graph for this function:

Here is the caller graph for this function:

float __ieee754_acoshf ( float  )

Definition at line 34 of file e_acoshf.c.

{      
       float t;
       int32_t hx;
       GET_FLOAT_WORD(hx,x);
       if(hx<0x3f800000) {         /* x < 1 */
           return (x-x)/(x-x);
       } else if(hx >=0x4d800000) {       /* x > 2**28 */
           if(hx >=0x7f800000) {   /* x is inf of NaN */
               return x+x;
           } else 
              return __ieee754_logf(x)+ln2;      /* acosh(huge)=log(2x) */
       } else if (hx==0x3f800000) {
           return 0.0;                    /* acosh(1) = 0 */
       } else if (hx > 0x40000000) {      /* 2**28 > x > 2 */
           t=x*x;
           return __ieee754_logf((float)2.0*x-one/(x+__ieee754_sqrtf(t-one)));
       } else {                    /* 1<x<2 */
           t = x-one;
           return __log1pf(t+__sqrtf((float)2.0*t+t*t));
       }
}

Here is the call graph for this function:

Here is the caller graph for this function:

long double __ieee754_acoshl ( long  double)

Definition at line 6 of file e_acoshl.c.

{
  fputs ("__ieee754_acoshl not implemented\n", stderr);
  __set_errno (ENOSYS);
  return 0.0;
}

Here is the call graph for this function:

Here is the caller graph for this function:

long double __ieee754_acosl ( long  double)

Definition at line 6 of file e_acosl.c.

{
  fputs ("__ieee754_acosl not implemented\n", stderr);
  __set_errno (ENOSYS);
  return 0.0;
}

Here is the call graph for this function:

Here is the caller graph for this function:

double __ieee754_asin ( double  )

Definition at line 56 of file e_asin.c.

                               {
  double x1,x2,xx,s1,s2,res1,p,t,res,r,cor,cc,y,c,z,w[2];
  mynumber u,v;
  int4 k,m,n;
#if 0
  int4 nn;
#endif

  u.x = x;
  m = u.i[HIGH_HALF];
  k = 0x7fffffff&m;              /* no sign */

  if (k < 0x3e500000) return x;  /* for x->0 => sin(x)=x */
  /*----------------------2^-26 <= |x| < 2^ -3    -----------------*/
  else
  if (k < 0x3fc00000) {
    x2 = x*x;
    t = (((((f6*x2 + f5)*x2 + f4)*x2 + f3)*x2 + f2)*x2 + f1)*(x2*x);
    res = x+t;         /*  res=arcsin(x) according to Taylor series  */
    cor = (x-res)+t;
    if (res == res+1.025*cor) return res;
    else {
      x1 = x+big;
      xx = x*x;
      x1 -= big;
      x2 = x - x1;
      p = x1*x1*x1;
      s1 = a1.x*p;
      s2 = ((((((c7*xx + c6)*xx + c5)*xx + c4)*xx + c3)*xx + c2)*xx*xx*x +
            ((a1.x+a2.x)*x2*x2+ 0.5*x1*x)*x2) + a2.x*p;
      res1 = x+s1;
      s2 = ((x-res1)+s1)+s2;
      res = res1+s2;
      cor = (res1-res)+s2;
      if (res == res+1.00014*cor) return res;
      else {
       __doasin(x,0,w);
       if (w[0]==(w[0]+1.00000001*w[1])) return w[0];
       else {
         y=ABS(x);
         res=ABS(w[0]);
         res1=ABS(w[0]+1.1*w[1]);
         return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1);
       }
      }
    }
  }
  /*---------------------0.125 <= |x| < 0.5 -----------------------------*/
  else if (k < 0x3fe00000) {
    if (k<0x3fd00000) n = 11*((k&0x000fffff)>>15);
    else n = 11*((k&0x000fffff)>>14)+352;
    if (m>0) xx = x - asncs.x[n];
    else xx = -x - asncs.x[n];
    t = asncs.x[n+1]*xx;
    p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+xx*(asncs.x[n+5]
     +xx*asncs.x[n+6]))))+asncs.x[n+7];
    t+=p;
    res =asncs.x[n+8] +t;
    cor = (asncs.x[n+8]-res)+t;
    if (res == res+1.05*cor) return (m>0)?res:-res;
    else {
      r=asncs.x[n+8]+xx*asncs.x[n+9];
      t=((asncs.x[n+8]-r)+xx*asncs.x[n+9])+(p+xx*asncs.x[n+10]);
      res = r+t;
      cor = (r-res)+t;
      if (res == res+1.0005*cor) return (m>0)?res:-res;
      else {
       res1=res+1.1*cor;
       z=0.5*(res1-res);
       __dubsin(res,z,w);
       z=(w[0]-ABS(x))+w[1];
       if (z>1.0e-27) return (m>0)?min(res,res1):-min(res,res1);
       else if (z<-1.0e-27) return (m>0)?max(res,res1):-max(res,res1);
       else {
         y=ABS(x);
         return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1);
       }
      }
    }
  }    /*   else  if (k < 0x3fe00000)    */
  /*-------------------- 0.5 <= |x| < 0.75 -----------------------------*/
  else
  if (k < 0x3fe80000) {
    n = 1056+((k&0x000fe000)>>11)*3;
    if (m>0) xx = x - asncs.x[n];
    else xx = -x - asncs.x[n];
    t = asncs.x[n+1]*xx;
    p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+xx*(asncs.x[n+5]
          +xx*(asncs.x[n+6]+xx*asncs.x[n+7])))))+asncs.x[n+8];
    t+=p;
    res =asncs.x[n+9] +t;
    cor = (asncs.x[n+9]-res)+t;
    if (res == res+1.01*cor) return (m>0)?res:-res;
    else {
      r=asncs.x[n+9]+xx*asncs.x[n+10];
      t=((asncs.x[n+9]-r)+xx*asncs.x[n+10])+(p+xx*asncs.x[n+11]);
      res = r+t;
      cor = (r-res)+t;
      if (res == res+1.0005*cor) return (m>0)?res:-res;
      else {
       res1=res+1.1*cor;
       z=0.5*(res1-res);
       __dubsin(res,z,w);
       z=(w[0]-ABS(x))+w[1];
       if (z>1.0e-27) return (m>0)?min(res,res1):-min(res,res1);
       else if (z<-1.0e-27) return (m>0)?max(res,res1):-max(res,res1);
       else {
         y=ABS(x);
         return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1);
       }
      }
    }
  }    /*   else  if (k < 0x3fe80000)    */
  /*--------------------- 0.75 <= |x|< 0.921875 ----------------------*/
  else
  if (k < 0x3fed8000) {
    n = 992+((k&0x000fe000)>>13)*13;
    if (m>0) xx = x - asncs.x[n];
    else xx = -x - asncs.x[n];
    t = asncs.x[n+1]*xx;
    p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+xx*(asncs.x[n+5]
     +xx*(asncs.x[n+6]+xx*(asncs.x[n+7]+xx*asncs.x[n+8]))))))+asncs.x[n+9];
    t+=p;
    res =asncs.x[n+10] +t;
    cor = (asncs.x[n+10]-res)+t;
    if (res == res+1.01*cor) return (m>0)?res:-res;
    else {
      r=asncs.x[n+10]+xx*asncs.x[n+11];
      t=((asncs.x[n+10]-r)+xx*asncs.x[n+11])+(p+xx*asncs.x[n+12]);
      res = r+t;
      cor = (r-res)+t;
      if (res == res+1.0008*cor) return (m>0)?res:-res;
      else {
       res1=res+1.1*cor;
       z=0.5*(res1-res);
       y=hp0.x-res;
       z=((hp0.x-y)-res)+(hp1.x-z);
       __dubcos(y,z,w);
       z=(w[0]-ABS(x))+w[1];
       if (z>1.0e-27) return (m>0)?min(res,res1):-min(res,res1);
       else if (z<-1.0e-27) return (m>0)?max(res,res1):-max(res,res1);
       else {
         y=ABS(x);
         return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1);
       }
      }
    }
  }    /*   else  if (k < 0x3fed8000)    */
  /*-------------------0.921875 <= |x| < 0.953125 ------------------------*/
  else
  if (k < 0x3fee8000) {
    n = 884+((k&0x000fe000)>>13)*14;
    if (m>0) xx = x - asncs.x[n];
    else xx = -x - asncs.x[n];
    t = asncs.x[n+1]*xx;
    p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
                      xx*(asncs.x[n+5]+xx*(asncs.x[n+6]
                    +xx*(asncs.x[n+7]+xx*(asncs.x[n+8]+
                      xx*asncs.x[n+9])))))))+asncs.x[n+10];
    t+=p;
    res =asncs.x[n+11] +t;
    cor = (asncs.x[n+11]-res)+t;
    if (res == res+1.01*cor) return (m>0)?res:-res;
    else {
      r=asncs.x[n+11]+xx*asncs.x[n+12];
      t=((asncs.x[n+11]-r)+xx*asncs.x[n+12])+(p+xx*asncs.x[n+13]);
      res = r+t;
      cor = (r-res)+t;
      if (res == res+1.0007*cor) return (m>0)?res:-res;
      else {
       res1=res+1.1*cor;
       z=0.5*(res1-res);
       y=(hp0.x-res)-z;
       z=y+hp1.x;
       y=(y-z)+hp1.x;
       __dubcos(z,y,w);
       z=(w[0]-ABS(x))+w[1];
       if (z>1.0e-27) return (m>0)?min(res,res1):-min(res,res1);
       else if (z<-1.0e-27) return (m>0)?max(res,res1):-max(res,res1);
       else {
         y=ABS(x);
         return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1);
       }
      }
    }
  }    /*   else  if (k < 0x3fee8000)    */

  /*--------------------0.953125 <= |x| < 0.96875 ------------------------*/
  else
  if (k < 0x3fef0000) {
    n = 768+((k&0x000fe000)>>13)*15;
    if (m>0) xx = x - asncs.x[n];
    else xx = -x - asncs.x[n];
    t = asncs.x[n+1]*xx;
    p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
                         xx*(asncs.x[n+5]+xx*(asncs.x[n+6]
                      +xx*(asncs.x[n+7]+xx*(asncs.x[n+8]+
                    xx*(asncs.x[n+9]+xx*asncs.x[n+10]))))))))+asncs.x[n+11];
    t+=p;
    res =asncs.x[n+12] +t;
    cor = (asncs.x[n+12]-res)+t;
    if (res == res+1.01*cor) return (m>0)?res:-res;
    else {
      r=asncs.x[n+12]+xx*asncs.x[n+13];
      t=((asncs.x[n+12]-r)+xx*asncs.x[n+13])+(p+xx*asncs.x[n+14]);
      res = r+t;
      cor = (r-res)+t;
      if (res == res+1.0007*cor) return (m>0)?res:-res;
      else {
       res1=res+1.1*cor;
       z=0.5*(res1-res);
       y=(hp0.x-res)-z;
       z=y+hp1.x;
       y=(y-z)+hp1.x;
       __dubcos(z,y,w);
       z=(w[0]-ABS(x))+w[1];
       if (z>1.0e-27) return (m>0)?min(res,res1):-min(res,res1);
       else if (z<-1.0e-27) return (m>0)?max(res,res1):-max(res,res1);
       else {
         y=ABS(x);
         return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1);
       }
      }
    }
  }    /*   else  if (k < 0x3fef0000)    */
  /*--------------------0.96875 <= |x| < 1 --------------------------------*/
  else
  if (k<0x3ff00000)  {
    z = 0.5*((m>0)?(1.0-x):(1.0+x));
    v.x=z;
    k=v.i[HIGH_HALF];
    t=inroot[(k&0x001fffff)>>14]*powtwo[511-(k>>21)];
    r=1.0-t*t*z;
    t = t*(rt0+r*(rt1+r*(rt2+r*rt3)));
    c=t*z;
    t=c*(1.5-0.5*t*c);
    y=(c+t24)-t24;
    cc = (z-y*y)/(t+y);
    p=(((((f6*z+f5)*z+f4)*z+f3)*z+f2)*z+f1)*z;
    cor = (hp1.x - 2.0*cc)-2.0*(y+cc)*p;
    res1 = hp0.x - 2.0*y;
    res =res1 + cor;
    if (res == res+1.003*((res1-res)+cor)) return (m>0)?res:-res;
    else {
      c=y+cc;
      cc=(y-c)+cc;
      __doasin(c,cc,w);
      res1=hp0.x-2.0*w[0];
      cor=((hp0.x-res1)-2.0*w[0])+(hp1.x-2.0*w[1]);
      res = res1+cor;
      cor = (res1-res)+cor;
      if (res==(res+1.0000001*cor)) return (m>0)?res:-res;
      else {
       y=ABS(x);
       res1=res+1.1*cor;
       return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1);
      }
    }
  }    /*   else  if (k < 0x3ff00000)    */
  /*---------------------------- |x|>=1 -------------------------------*/
  else if (k==0x3ff00000 && u.i[LOW_HALF]==0) return (m>0)?hp0.x:-hp0.x;
  else
  if (k>0x7ff00000 || (k == 0x7ff00000 && u.i[LOW_HALF] != 0)) return x;
  else {
    u.i[HIGH_HALF]=0x7ff00000;
    v.i[HIGH_HALF]=0x7ff00000;
    u.i[LOW_HALF]=0;
    v.i[LOW_HALF]=0;
    return u.x/v.x;  /* NaN */
 }
}

Here is the call graph for this function:

Here is the caller graph for this function:

float __ieee754_asinf ( float  )

Definition at line 69 of file e_asinf.c.

{
       float t,w,p,q,c,r,s;
       int32_t hx,ix;
       GET_FLOAT_WORD(hx,x);
       ix = hx&0x7fffffff;
       if(ix==0x3f800000) {
              /* asin(1)=+-pi/2 with inexact */
           return x*pio2_hi+x*pio2_lo;
       } else if(ix> 0x3f800000) { /* |x|>= 1 */
           return (x-x)/(x-x);            /* asin(|x|>1) is NaN */
       } else if (ix<0x3f000000) { /* |x|<0.5 */
           if(ix<0x32000000) {            /* if |x| < 2**-27 */
              if(huge+x>one) return x;/* return x with inexact if x!=0*/
           } else {
              t = x*x;
              w = t * (p0 + t * (p1 + t * (p2 + t * (p3 + t * p4))));
              return x+x*w;
           }
       }
       /* 1> |x|>= 0.5 */
       w = one-fabsf(x);
       t = w*0.5f;
       p = t * (p0 + t * (p1 + t * (p2 + t * (p3 + t * p4))));
       s = __ieee754_sqrtf(t);
       if(ix>=0x3F79999A) {        /* if |x| > 0.975 */
           t = pio2_hi-(2.0f*(s+s*p)-pio2_lo);
       } else {
           int32_t iw;
           w  = s;
           GET_FLOAT_WORD(iw,w);
           SET_FLOAT_WORD(w,iw&0xfffff000);
           c  = (t-w*w)/(s+w);
           r  = p;
           p  = 2.0f*s*r-(pio2_lo-2.0f*c);
           q  = pio4_hi-2.0f*w;
           t  = pio4_hi-(p-q);
       }
       if(hx>0) return t; else return -t;
}

Here is the call graph for this function:

Here is the caller graph for this function:

long double __ieee754_asinl ( long  double)

Definition at line 6 of file e_asinl.c.

{
  fputs ("__ieee754_asinl not implemented\n", stderr);
  __set_errno (ENOSYS);
  return 0.0;
}

Here is the call graph for this function:

Here is the caller graph for this function:

double __ieee754_atan2 ( double  ,
double   
)

Definition at line 58 of file e_atan2.c.

                                          {

  int i,de,ux,dx,uy,dy;
#if 0
  int p;
#endif
  static const int pr[MM]={6,8,10,20,32};
  double ax,ay,u,du,u9,ua,v,vv,dv,t1,t2,t3,t4,t5,t6,t7,t8,
         z,zz,cor,s1,ss1,s2,ss2;
#if 0
  double z1,z2;
#endif
  number num;
#if 0
  mp_no mperr,mpt1,mpx,mpy,mpz,mpz1,mpz2;
#endif

  static const int ep= 59768832,   /*  57*16**5   */
                   em=-59768832;   /* -57*16**5   */

  /* x=NaN or y=NaN */
  num.d = x;  ux = num.i[HIGH_HALF];  dx = num.i[LOW_HALF];
  if   ((ux&0x7ff00000)    ==0x7ff00000) {
    if (((ux&0x000fffff)|dx)!=0x00000000) return x+x; }
  num.d = y;  uy = num.i[HIGH_HALF];  dy = num.i[LOW_HALF];
  if   ((uy&0x7ff00000)    ==0x7ff00000) {
    if (((uy&0x000fffff)|dy)!=0x00000000) return y+y; }

  /* y=+-0 */
  if      (uy==0x00000000) {
    if    (dy==0x00000000) {
      if  ((ux&0x80000000)==0x00000000)  return ZERO;
      else                               return opi.d; } }
  else if (uy==0x80000000) {
    if    (dy==0x00000000) {
      if  ((ux&0x80000000)==0x00000000)  return MZERO;
      else                               return mopi.d;} }

  /* x=+-0 */
  if (x==ZERO) {
    if ((uy&0x80000000)==0x00000000)     return hpi.d;
    else                                 return mhpi.d; }

  /* x=+-INF */
  if          (ux==0x7ff00000) {
    if        (dx==0x00000000) {
      if      (uy==0x7ff00000) {
        if    (dy==0x00000000)  return qpi.d; }
      else if (uy==0xfff00000) {
        if    (dy==0x00000000)  return mqpi.d; }
      else {
        if    ((uy&0x80000000)==0x00000000)  return ZERO;
        else                                 return MZERO; }
    }
  }
  else if     (ux==0xfff00000) {
    if        (dx==0x00000000) {
      if      (uy==0x7ff00000) {
        if    (dy==0x00000000)  return tqpi.d; }
      else if (uy==0xfff00000) {
        if    (dy==0x00000000)  return mtqpi.d; }
      else                     {
        if    ((uy&0x80000000)==0x00000000)  return opi.d;
        else                                 return mopi.d; }
    }
  }

  /* y=+-INF */
  if      (uy==0x7ff00000) {
    if    (dy==0x00000000)  return hpi.d; }
  else if (uy==0xfff00000) {
    if    (dy==0x00000000)  return mhpi.d; }

  /* either x/y or y/x is very close to zero */
  ax = (x<ZERO) ? -x : x;    ay = (y<ZERO) ? -y : y;
  de = (uy & 0x7ff00000) - (ux & 0x7ff00000);
  if      (de>=ep)  { return ((y>ZERO) ? hpi.d : mhpi.d); }
  else if (de<=em)  {
    if    (x>ZERO)  {
      if  ((z=ay/ax)<TWOM1022)  return normalized(ax,ay,y,z);
      else                      return signArctan2(y,z); }
    else            { return ((y>ZERO) ? opi.d : mopi.d); } }

  /* if either x or y is extremely close to zero, scale abs(x), abs(y). */
  if (ax<twom500.d || ay<twom500.d) { ax*=two500.d;  ay*=two500.d; }

  /* x,y which are neither special nor extreme */
  if (ay<ax) {
    u=ay/ax;
    EMULV(ax,u,v,vv,t1,t2,t3,t4,t5)
    du=((ay-v)-vv)/ax; }
  else {
    u=ax/ay;
    EMULV(ay,u,v,vv,t1,t2,t3,t4,t5)
    du=((ax-v)-vv)/ay; }

  if (x>ZERO) {

    /* (i)   x>0, abs(y)< abs(x):  atan(ay/ax) */
    if (ay<ax) {
      if (u<inv16.d) {
        v=u*u;  zz=du+u*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d)))));
        if ((z=u+(zz-u1.d*u)) == u+(zz+u1.d*u))  return signArctan2(y,z);

        MUL2(u,du,u,du,v,vv,t1,t2,t3,t4,t5,t6,t7,t8)
        s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d))));
        ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2)
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
        ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2)
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
        ADD2(f5.d,ff5.d,s1,ss1,s2,ss2,t1,t2)
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
        ADD2(f3.d,ff3.d,s1,ss1,s2,ss2,t1,t2)
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
        MUL2(u,du,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
        ADD2(u,du,s2,ss2,s1,ss1,t1,t2)
        if ((z=s1+(ss1-u5.d*s1)) == s1+(ss1+u5.d*s1))  return signArctan2(y,z);
        return atan2Mp(x,y,pr);
      }
      else {
        i=(TWO52+TWO8*u)-TWO52;  i-=16;
        t3=u-cij[i][0].d;
        EADD(t3,du,v,dv)
        t1=cij[i][1].d;  t2=cij[i][2].d;
        zz=v*t2+(dv*t2+v*v*(cij[i][3].d+v*(cij[i][4].d+
                         v*(cij[i][5].d+v* cij[i][6].d))));
        if (i<112) {
          if (i<48)  u9=u91.d;    /* u < 1/4        */
          else       u9=u92.d; }  /* 1/4 <= u < 1/2 */
        else {
          if (i<176) u9=u93.d;    /* 1/2 <= u < 3/4 */
          else       u9=u94.d; }  /* 3/4 <= u <= 1  */
        if ((z=t1+(zz-u9*t1)) == t1+(zz+u9*t1))  return signArctan2(y,z);

        t1=u-hij[i][0].d;
        EADD(t1,du,v,vv)
        s1=v*(hij[i][11].d+v*(hij[i][12].d+v*(hij[i][13].d+
           v*(hij[i][14].d+v* hij[i][15].d))));
        ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2)
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
        ADD2(hij[i][7].d,hij[i][8].d,s1,ss1,s2,ss2,t1,t2)
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
        ADD2(hij[i][5].d,hij[i][6].d,s1,ss1,s2,ss2,t1,t2)
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
        ADD2(hij[i][3].d,hij[i][4].d,s1,ss1,s2,ss2,t1,t2)
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
        ADD2(hij[i][1].d,hij[i][2].d,s1,ss1,s2,ss2,t1,t2)
        if ((z=s2+(ss2-ub.d*s2)) == s2+(ss2+ub.d*s2))  return signArctan2(y,z);
        return atan2Mp(x,y,pr);
      }
    }

    /* (ii)  x>0, abs(x)<=abs(y):  pi/2-atan(ax/ay) */
    else {
      if (u<inv16.d) {
        v=u*u;
        zz=u*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d)))));
        ESUB(hpi.d,u,t2,cor)
        t3=((hpi1.d+cor)-du)-zz;
        if ((z=t2+(t3-u2.d)) == t2+(t3+u2.d))  return signArctan2(y,z);

        MUL2(u,du,u,du,v,vv,t1,t2,t3,t4,t5,t6,t7,t8)
        s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d))));
        ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2)
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
        ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2)
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
        ADD2(f5.d,ff5.d,s1,ss1,s2,ss2,t1,t2)
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
        ADD2(f3.d,ff3.d,s1,ss1,s2,ss2,t1,t2)
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
        MUL2(u,du,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
        ADD2(u,du,s2,ss2,s1,ss1,t1,t2)
        SUB2(hpi.d,hpi1.d,s1,ss1,s2,ss2,t1,t2)
        if ((z=s2+(ss2-u6.d)) == s2+(ss2+u6.d))  return signArctan2(y,z);
        return atan2Mp(x,y,pr);
      }
      else {
        i=(TWO52+TWO8*u)-TWO52;  i-=16;
        v=(u-cij[i][0].d)+du;
        zz=hpi1.d-v*(cij[i][2].d+v*(cij[i][3].d+v*(cij[i][4].d+
                                 v*(cij[i][5].d+v* cij[i][6].d))));
        t1=hpi.d-cij[i][1].d;
        if (i<112)  ua=ua1.d;  /* w <  1/2 */
        else        ua=ua2.d;  /* w >= 1/2 */
        if ((z=t1+(zz-ua)) == t1+(zz+ua))  return signArctan2(y,z);

        t1=u-hij[i][0].d;
        EADD(t1,du,v,vv)
        s1=v*(hij[i][11].d+v*(hij[i][12].d+v*(hij[i][13].d+
           v*(hij[i][14].d+v* hij[i][15].d))));
        ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2)
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
        ADD2(hij[i][7].d,hij[i][8].d,s1,ss1,s2,ss2,t1,t2)
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
        ADD2(hij[i][5].d,hij[i][6].d,s1,ss1,s2,ss2,t1,t2)
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
        ADD2(hij[i][3].d,hij[i][4].d,s1,ss1,s2,ss2,t1,t2)
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
        ADD2(hij[i][1].d,hij[i][2].d,s1,ss1,s2,ss2,t1,t2)
        SUB2(hpi.d,hpi1.d,s2,ss2,s1,ss1,t1,t2)
        if ((z=s1+(ss1-uc.d)) == s1+(ss1+uc.d))  return signArctan2(y,z);
        return atan2Mp(x,y,pr);
      }
    }
  }
  else {

    /* (iii) x<0, abs(x)< abs(y):  pi/2+atan(ax/ay) */
    if (ax<ay) {
      if (u<inv16.d) {
        v=u*u;
        zz=u*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d)))));
        EADD(hpi.d,u,t2,cor)
        t3=((hpi1.d+cor)+du)+zz;
        if ((z=t2+(t3-u3.d)) == t2+(t3+u3.d))  return signArctan2(y,z);

        MUL2(u,du,u,du,v,vv,t1,t2,t3,t4,t5,t6,t7,t8)
        s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d))));
        ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2)
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
        ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2)
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
        ADD2(f5.d,ff5.d,s1,ss1,s2,ss2,t1,t2)
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
        ADD2(f3.d,ff3.d,s1,ss1,s2,ss2,t1,t2)
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
        MUL2(u,du,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
        ADD2(u,du,s2,ss2,s1,ss1,t1,t2)
        ADD2(hpi.d,hpi1.d,s1,ss1,s2,ss2,t1,t2)
        if ((z=s2+(ss2-u7.d)) == s2+(ss2+u7.d))  return signArctan2(y,z);
        return atan2Mp(x,y,pr);
      }
      else {
        i=(TWO52+TWO8*u)-TWO52;  i-=16;
        v=(u-cij[i][0].d)+du;
        zz=hpi1.d+v*(cij[i][2].d+v*(cij[i][3].d+v*(cij[i][4].d+
                                 v*(cij[i][5].d+v* cij[i][6].d))));
        t1=hpi.d+cij[i][1].d;
        if (i<112)  ua=ua1.d;  /* w <  1/2 */
        else        ua=ua2.d;  /* w >= 1/2 */
        if ((z=t1+(zz-ua)) == t1+(zz+ua))  return signArctan2(y,z);

        t1=u-hij[i][0].d;
        EADD(t1,du,v,vv)
        s1=v*(hij[i][11].d+v*(hij[i][12].d+v*(hij[i][13].d+
           v*(hij[i][14].d+v* hij[i][15].d))));
        ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2)
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
        ADD2(hij[i][7].d,hij[i][8].d,s1,ss1,s2,ss2,t1,t2)
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
        ADD2(hij[i][5].d,hij[i][6].d,s1,ss1,s2,ss2,t1,t2)
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
        ADD2(hij[i][3].d,hij[i][4].d,s1,ss1,s2,ss2,t1,t2)
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
        ADD2(hij[i][1].d,hij[i][2].d,s1,ss1,s2,ss2,t1,t2)
        ADD2(hpi.d,hpi1.d,s2,ss2,s1,ss1,t1,t2)
        if ((z=s1+(ss1-uc.d)) == s1+(ss1+uc.d))  return signArctan2(y,z);
        return atan2Mp(x,y,pr);
      }
    }

    /* (iv)  x<0, abs(y)<=abs(x):  pi-atan(ax/ay) */
    else {
      if (u<inv16.d) {
        v=u*u;
        zz=u*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d)))));
        ESUB(opi.d,u,t2,cor)
        t3=((opi1.d+cor)-du)-zz;
        if ((z=t2+(t3-u4.d)) == t2+(t3+u4.d))  return signArctan2(y,z);

        MUL2(u,du,u,du,v,vv,t1,t2,t3,t4,t5,t6,t7,t8)
        s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d))));
        ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2)
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
        ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2)
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
        ADD2(f5.d,ff5.d,s1,ss1,s2,ss2,t1,t2)
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
        ADD2(f3.d,ff3.d,s1,ss1,s2,ss2,t1,t2)
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
        MUL2(u,du,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
        ADD2(u,du,s2,ss2,s1,ss1,t1,t2)
        SUB2(opi.d,opi1.d,s1,ss1,s2,ss2,t1,t2)
        if ((z=s2+(ss2-u8.d)) == s2+(ss2+u8.d))  return signArctan2(y,z);
        return atan2Mp(x,y,pr);
      }
      else {
        i=(TWO52+TWO8*u)-TWO52;  i-=16;
        v=(u-cij[i][0].d)+du;
        zz=opi1.d-v*(cij[i][2].d+v*(cij[i][3].d+v*(cij[i][4].d+
                                 v*(cij[i][5].d+v* cij[i][6].d))));
        t1=opi.d-cij[i][1].d;
        if (i<112)  ua=ua1.d;  /* w <  1/2 */
        else        ua=ua2.d;  /* w >= 1/2 */
        if ((z=t1+(zz-ua)) == t1+(zz+ua))  return signArctan2(y,z);

        t1=u-hij[i][0].d;
        EADD(t1,du,v,vv)
        s1=v*(hij[i][11].d+v*(hij[i][12].d+v*(hij[i][13].d+
           v*(hij[i][14].d+v* hij[i][15].d))));
        ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2)
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
        ADD2(hij[i][7].d,hij[i][8].d,s1,ss1,s2,ss2,t1,t2)
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
        ADD2(hij[i][5].d,hij[i][6].d,s1,ss1,s2,ss2,t1,t2)
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
        ADD2(hij[i][3].d,hij[i][4].d,s1,ss1,s2,ss2,t1,t2)
        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
        ADD2(hij[i][1].d,hij[i][2].d,s1,ss1,s2,ss2,t1,t2)
        SUB2(opi.d,opi1.d,s2,ss2,s1,ss1,t1,t2)
        if ((z=s1+(ss1-uc.d)) == s1+(ss1+uc.d))  return signArctan2(y,z);
        return atan2Mp(x,y,pr);
      }
    }
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

float __ieee754_atan2f ( float  ,
float   
)

Definition at line 38 of file e_atan2f.c.

{
       float z;
       int32_t k,m,hx,hy,ix,iy;

       GET_FLOAT_WORD(hx,x);
       ix = hx&0x7fffffff;
       GET_FLOAT_WORD(hy,y);
       iy = hy&0x7fffffff;
       if((ix>0x7f800000)||
          (iy>0x7f800000))  /* x or y is NaN */
          return x+y;
       if(hx==0x3f800000) return __atanf(y);   /* x=1.0 */
       m = ((hy>>31)&1)|((hx>>30)&2);     /* 2*sign(x)+sign(y) */

    /* when y = 0 */
       if(iy==0) {
           switch(m) {
              case 0:
              case 1: return y;    /* atan(+-0,+anything)=+-0 */
              case 2: return  pi+tiny;/* atan(+0,-anything) = pi */
              case 3: return -pi-tiny;/* atan(-0,-anything) =-pi */
           }
       }
    /* when x = 0 */
       if(ix==0) return (hy<0)?  -pi_o_2-tiny: pi_o_2+tiny;

    /* when x is INF */
       if(ix==0x7f800000) {
           if(iy==0x7f800000) {
              switch(m) {
                  case 0: return  pi_o_4+tiny;/* atan(+INF,+INF) */
                  case 1: return -pi_o_4-tiny;/* atan(-INF,+INF) */
                  case 2: return  (float)3.0*pi_o_4+tiny;/*atan(+INF,-INF)*/
                  case 3: return (float)-3.0*pi_o_4-tiny;/*atan(-INF,-INF)*/
              }
           } else {
              switch(m) {
                  case 0: return  zero  ; /* atan(+...,+INF) */
                  case 1: return -zero  ; /* atan(-...,+INF) */
                  case 2: return  pi+tiny  ;     /* atan(+...,-INF) */
                  case 3: return -pi-tiny  ;     /* atan(-...,-INF) */
              }
           }
       }
    /* when y is INF */
       if(iy==0x7f800000) return (hy<0)? -pi_o_2-tiny: pi_o_2+tiny;

    /* compute y/x */
       k = (iy-ix)>>23;
       if(k > 60) z=pi_o_2+(float)0.5*pi_lo;     /* |y/x| >  2**60 */
       else if(hx<0&&k<-60) z=0.0;        /* |y|/x < -2**60 */
       else z=__atanf(fabsf(y/x)); /* safe to do y/x */
       switch (m) {
           case 0: return       z  ;      /* atan(+,+) */
           case 1: {
                    u_int32_t zh;
                    GET_FLOAT_WORD(zh,z);
                    SET_FLOAT_WORD(z,zh ^ 0x80000000);
                  }
                  return       z  ;       /* atan(-,+) */
           case 2: return  pi-(z-pi_lo);/* atan(+,-) */
           default: /* case 3 */
                  return  (z-pi_lo)-pi;/* atan(-,-) */
       }
}

Here is the call graph for this function:

Here is the caller graph for this function:

long double __ieee754_atan2l ( long  double,
long  double 
)

Definition at line 6 of file e_atan2l.c.

{
  fputs ("__ieee754_atan2l not implemented\n", stderr);
  __set_errno (ENOSYS);
  return 0.0;
}

Here is the call graph for this function:

Here is the caller graph for this function:

double __ieee754_atanh ( double  )

Definition at line 53 of file e_atanh.c.

{
       double t;
       int32_t hx,ix;
       u_int32_t lx;
       EXTRACT_WORDS(hx,lx,x);
       ix = hx&0x7fffffff;
       if ((ix|((lx|(-lx))>>31))>0x3ff00000) /* |x|>1 */
           return (x-x)/(x-x);
       if(ix==0x3ff00000) 
           return x/zero;
       if(ix<0x3e300000&&(huge+x)>zero) return x;       /* x<2**-28 */
       SET_HIGH_WORD(x,ix);
       if(ix<0x3fe00000) {         /* x < 0.5 */
           t = x+x;
           t = 0.5*__log1p(t+t*x/(one-x));
       } else 
           t = 0.5*__log1p((x+x)/(one-x));
       if(hx>=0) return t; else return -t;
}

Here is the call graph for this function:

Here is the caller graph for this function:

float __ieee754_atanhf ( float  )

Definition at line 38 of file e_atanhf.c.

{
       float t;
       int32_t hx,ix;
       GET_FLOAT_WORD(hx,x);
       ix = hx&0x7fffffff;
       if (ix>0x3f800000)          /* |x|>1 */
           return (x-x)/(x-x);
       if(ix==0x3f800000) 
           return x/zero;
       if(ix<0x31800000&&(huge+x)>zero) return x;       /* x<2**-28 */
       SET_FLOAT_WORD(x,ix);
       if(ix<0x3f000000) {         /* x < 0.5 */
           t = x+x;
           t = (float)0.5*__log1pf(t+t*x/(one-x));
       } else 
           t = (float)0.5*__log1pf((x+x)/(one-x));
       if(hx>=0) return t; else return -t;
}

Here is the call graph for this function:

Here is the caller graph for this function:

long double __ieee754_atanhl ( long  double)

Definition at line 6 of file e_atanhl.c.

{
  fputs ("__ieee754_atanhl not implemented\n", stderr);
  __set_errno (ENOSYS);
  return 0.0;
}

Here is the call graph for this function:

Here is the caller graph for this function:

double __ieee754_cosh ( double  )

Definition at line 50 of file e_cosh.c.

{      
       double t,w;
       int32_t ix;
       u_int32_t lx;

    /* High word of |x|. */
       GET_HIGH_WORD(ix,x);
       ix &= 0x7fffffff;

    /* x is INF or NaN */
       if(ix>=0x7ff00000) return x*x;     

    /* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */
       if(ix<0x3fd62e43) {
           t = __expm1(fabs(x));
           w = one+t;
           if (ix<0x3c800000) return w;   /* cosh(tiny) = 1 */
           return one+(t*t)/(w+w);
       }

    /* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */
       if (ix < 0x40360000) {
              t = __ieee754_exp(fabs(x));
              return half*t+half/t;
       }

    /* |x| in [22, log(maxdouble)] return half*exp(|x|) */
       if (ix < 0x40862e42)  return half*__ieee754_exp(fabs(x));

    /* |x| in [log(maxdouble), overflowthresold] */
       GET_LOW_WORD(lx,x);
       if (ix<0x408633ce || ((ix==0x408633ce)&&(lx<=(u_int32_t)0x8fb9f87d))) {
           w = __ieee754_exp(half*fabs(x));
           t = half*w;
           return t*w;
       }

    /* |x| > overflowthresold, cosh(x) overflow */
       return huge*huge;
}

Here is the call graph for this function:

Here is the caller graph for this function:

float __ieee754_coshf ( float  )

Definition at line 33 of file e_coshf.c.

{
       float t,w;
       int32_t ix;

       GET_FLOAT_WORD(ix,x);
       ix &= 0x7fffffff;

    /* x is INF or NaN */
       if(ix>=0x7f800000) return x*x;

    /* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */
       if(ix<0x3eb17218) {
           t = __expm1f(fabsf(x));
           w = one+t;
           if (ix<0x24000000) return w;   /* cosh(tiny) = 1 */
           return one+(t*t)/(w+w);
       }

    /* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */
       if (ix < 0x41b00000) {
              t = __ieee754_expf(fabsf(x));
              return half*t+half/t;
       }

    /* |x| in [22, log(maxdouble)] return half*exp(|x|) */
       if (ix < 0x42b17180)  return half*__ieee754_expf(fabsf(x));

    /* |x| in [log(maxdouble), overflowthresold] */
       if (ix<=0x42b2d4fc) {
           w = __ieee754_expf(half*fabsf(x));
           t = half*w;
           return t*w;
       }

    /* |x| > overflowthresold, cosh(x) overflow */
       return huge*huge;
}

Here is the call graph for this function:

Here is the caller graph for this function:

long double __ieee754_coshl ( long  double)

Definition at line 6 of file e_coshl.c.

{
  fputs ("__ieee754_coshl not implemented\n", stderr);
  __set_errno (ENOSYS);
  return 0.0;
}

Here is the call graph for this function:

Here is the caller graph for this function:

double __ieee754_exp ( double  )

Definition at line 49 of file e_exp.c.

                               {
  double bexp, t, eps, del, base, y, al, bet, res, rem, cor;
  mynumber junk1, junk2, binexp  = {{0,0}};
#if 0
  int4 k;
#endif
  int4 i,j,m,n,ex;

  junk1.x = x;
  m = junk1.i[HIGH_HALF];
  n = m&hugeint;

  if (n > smallint && n < bigint) {

    y = x*log2e.x + three51.x;
    bexp = y - three51.x;      /*  multiply the result by 2**bexp        */

    junk1.x = y;

    eps = bexp*ln_two2.x;      /* x = bexp*ln(2) + t - eps               */
    t = x - bexp*ln_two1.x;

    y = t + three33.x;
    base = y - three33.x;      /* t rounded to a multiple of 2**-18      */
    junk2.x = y;
    del = (t - base) - eps;    /*  x = bexp*ln(2) + base + del           */
    eps = del + del*del*(p3.x*del + p2.x);

    binexp.i[HIGH_HALF] =(junk1.i[LOW_HALF]+1023)<<20;

    i = ((junk2.i[LOW_HALF]>>8)&0xfffffffe)+356;
    j = (junk2.i[LOW_HALF]&511)<<1;

    al = coar.x[i]*fine.x[j];
    bet =(coar.x[i]*fine.x[j+1] + coar.x[i+1]*fine.x[j]) + coar.x[i+1]*fine.x[j+1];

    rem=(bet + bet*eps)+al*eps;
    res = al + rem;
    cor = (al - res) + rem;
    if  (res == (res+cor*err_0)) return res*binexp.x;
    else return __slowexp(x); /*if error is over bound */
  }

  if (n <= smallint) return 1.0;

  if (n >= badint) {
    if (n > infint) return(x+x);               /* x is NaN */
    if (n < infint) return ( (x>0) ? (hhuge*hhuge) : (tiny*tiny) );
    /* x is finite,  cause either overflow or underflow  */
    if (junk1.i[LOW_HALF] != 0)  return (x+x);                /*  x is NaN  */
    return ((x>0)?inf.x:zero );             /* |x| = inf;  return either inf or 0 */
  }

  y = x*log2e.x + three51.x;
  bexp = y - three51.x;
  junk1.x = y;
  eps = bexp*ln_two2.x;
  t = x - bexp*ln_two1.x;
  y = t + three33.x;
  base = y - three33.x;
  junk2.x = y;
  del = (t - base) - eps;
  eps = del + del*del*(p3.x*del + p2.x);
  i = ((junk2.i[LOW_HALF]>>8)&0xfffffffe)+356;
  j = (junk2.i[LOW_HALF]&511)<<1;
  al = coar.x[i]*fine.x[j];
  bet =(coar.x[i]*fine.x[j+1] + coar.x[i+1]*fine.x[j]) + coar.x[i+1]*fine.x[j+1];
  rem=(bet + bet*eps)+al*eps;
  res = al + rem;
  cor = (al - res) + rem;
  if (m>>31) {
    ex=junk1.i[LOW_HALF];
    if (res < 1.0) {res+=res; cor+=cor; ex-=1;}
    if (ex >=-1022) {
      binexp.i[HIGH_HALF] = (1023+ex)<<20;
      if  (res == (res+cor*err_0)) return res*binexp.x;
      else return __slowexp(x); /*if error is over bound */
    }
    ex = -(1022+ex);
    binexp.i[HIGH_HALF] = (1023-ex)<<20;
    res*=binexp.x;
    cor*=binexp.x;
    eps=1.0000000001+err_0*binexp.x;
    t=1.0+res;
    y = ((1.0-t)+res)+cor;
    res=t+y;
    cor = (t-res)+y;
    if (res == (res + eps*cor))
    { binexp.i[HIGH_HALF] = 0x00100000;
      return (res-1.0)*binexp.x;
    }
    else return __slowexp(x); /*   if error is over bound    */
  }
  else {
    binexp.i[HIGH_HALF] =(junk1.i[LOW_HALF]+767)<<20;
    if  (res == (res+cor*err_0)) return res*binexp.x*t256.x;
    else return __slowexp(x);
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

double __ieee754_exp10 ( double  )

Definition at line 25 of file e_exp10.c.

{
  /* This is a very stupid and inprecise implementation.  It'll get
     replaced sometime (soon?).  */
  return __ieee754_exp (M_LN10 * arg);
}

Here is the call graph for this function:

Here is the caller graph for this function:

float __ieee754_exp10f ( float  )

Definition at line 25 of file e_exp10f.c.

{
  /* This is a very stupid and inprecise implementation.  It'll get
     replaced sometime (soon?).  */
  return __ieee754_expf (M_LN10 * arg);
}

Here is the call graph for this function:

Here is the caller graph for this function:

long double __ieee754_exp10l ( long  double)

Definition at line 25 of file e_exp10l.c.

{
  /* This is a very stupid and inprecise implementation.  It'll get
     replaced sometime (soon?).  */
  return __ieee754_expl (M_LN10l * arg);
}

Here is the call graph for this function:

Here is the caller graph for this function:

double __ieee754_exp2 ( double  )

Definition at line 49 of file e_exp2.c.

{
  static const double himark = (double) DBL_MAX_EXP;
  static const double lomark = (double) (DBL_MIN_EXP - DBL_MANT_DIG - 1);

  /* Check for usual case.  */
  if (isless (x, himark) && isgreaterequal (x, lomark))
    {
      static const double THREEp42 = 13194139533312.0;
      int tval, unsafe;
      double rx, x22, result;
      union ieee754_double ex2_u, scale_u;
      fenv_t oldenv;

      feholdexcept (&oldenv);
#ifdef FE_TONEAREST
      /* If we don't have this, it's too bad.  */
      fesetround (FE_TONEAREST);
#endif

      /* 1. Argument reduction.
        Choose integers ex, -256 <= t < 256, and some real
        -1/1024 <= x1 <= 1024 so that
        x = ex + t/512 + x1.

        First, calculate rx = ex + t/512.  */
      rx = x + THREEp42;
      rx -= THREEp42;
      x -= rx;  /* Compute x=x1. */
      /* Compute tval = (ex*512 + t)+256.
        Now, t = (tval mod 512)-256 and ex=tval/512  [that's mod, NOT %; and
        /-round-to-nearest not the usual c integer /].  */
      tval = (int) (rx * 512.0 + 256.0);

      /* 2. Adjust for accurate table entry.
        Find e so that
        x = ex + t/512 + e + x2
        where -1e6 < e < 1e6, and
        (double)(2^(t/512+e))
        is accurate to one part in 2^-64.  */

      /* 'tval & 511' is the same as 'tval%512' except that it's always
        positive.
        Compute x = x2.  */
      x -= exp2_deltatable[tval & 511];

      /* 3. Compute ex2 = 2^(t/512+e+ex).  */
      ex2_u.d = exp2_accuratetable[tval & 511];
      tval >>= 9;
      unsafe = abs(tval) >= -DBL_MIN_EXP - 1;
      ex2_u.ieee.exponent += tval >> unsafe;
      scale_u.d = 1.0;
      scale_u.ieee.exponent += tval - (tval >> unsafe);

      /* 4. Approximate 2^x2 - 1, using a fourth-degree polynomial,
        with maximum error in [-2^-10-2^-30,2^-10+2^-30]
        less than 10^-19.  */

      x22 = (((.0096181293647031180
              * x + .055504110254308625)
             * x + .240226506959100583)
            * x + .69314718055994495) * ex2_u.d;

      /* 5. Return (2^x2-1) * 2^(t/512+e+ex) + 2^(t/512+e+ex).  */
      fesetenv (&oldenv);

      result = x22 * x + ex2_u.d;

      if (!unsafe)
       return result;
      else
       return result * scale_u.d;
    }
  /* Exceptional cases:  */
  else if (isless (x, himark))
    {
      if (__isinf (x))
       /* e^-inf == 0, with no error.  */
       return 0;
      else
       /* Underflow */
       return TWOM1000 * TWOM1000;
    }
  else
    /* Return x, if x is a NaN or Inf; or overflow, otherwise.  */
    return TWO1023*x;
}

Here is the call graph for this function:

Here is the caller graph for this function:

float __ieee754_exp2f ( float  )

Definition at line 45 of file e_exp2f.c.

{
  static const float himark = (float) FLT_MAX_EXP;
  static const float lomark = (float) (FLT_MIN_EXP - FLT_MANT_DIG - 1);

  /* Check for usual case.  */
  if (isless (x, himark) && isgreaterequal (x, lomark))
    {
      static const float THREEp14 = 49152.0;
      int tval, unsafe;
      float rx, x22, result;
      union ieee754_float ex2_u, scale_u;
      fenv_t oldenv;

      feholdexcept (&oldenv);
#ifdef FE_TONEAREST
      /* If we don't have this, it's too bad.  */
      fesetround (FE_TONEAREST);
#endif

      /* 1. Argument reduction.
        Choose integers ex, -128 <= t < 128, and some real
        -1/512 <= x1 <= 1/512 so that
        x = ex + t/512 + x1.

        First, calculate rx = ex + t/256.  */
      rx = x + THREEp14;
      rx -= THREEp14;
      x -= rx;  /* Compute x=x1. */
      /* Compute tval = (ex*256 + t)+128.
        Now, t = (tval mod 256)-128 and ex=tval/256  [that's mod, NOT %; and
        /-round-to-nearest not the usual c integer /].  */
      tval = (int) (rx * 256.0f + 128.0f);

      /* 2. Adjust for accurate table entry.
        Find e so that
        x = ex + t/256 + e + x2
        where -7e-4 < e < 7e-4, and
        (float)(2^(t/256+e))
        is accurate to one part in 2^-64.  */

      /* 'tval & 255' is the same as 'tval%256' except that it's always
        positive.
        Compute x = x2.  */
      x -= __exp2f_deltatable[tval & 255];

      /* 3. Compute ex2 = 2^(t/255+e+ex).  */
      ex2_u.f = __exp2f_atable[tval & 255];
      tval >>= 8;
      unsafe = abs(tval) >= -FLT_MIN_EXP - 1;
      ex2_u.ieee.exponent += tval >> unsafe;
      scale_u.f = 1.0;
      scale_u.ieee.exponent += tval - (tval >> unsafe);

      /* 4. Approximate 2^x2 - 1, using a second-degree polynomial,
        with maximum error in [-2^-9 - 2^-14, 2^-9 + 2^-14]
        less than 1.3e-10.  */

      x22 = (.24022656679f * x + .69314736128f) * ex2_u.f;

      /* 5. Return (2^x2-1) * 2^(t/512+e+ex) + 2^(t/512+e+ex).  */
      fesetenv (&oldenv);

      result = x22 * x + ex2_u.f;

      if (!unsafe)
       return result;
      else
       return result * scale_u.f;
    }
  /* Exceptional cases:  */
  else if (isless (x, himark))
    {
      if (__isinff (x))
       /* e^-inf == 0, with no error.  */
       return 0;
      else
       /* Underflow */
       return TWOM100 * TWOM100;
    }
  else
    /* Return x, if x is a NaN or Inf; or overflow, otherwise.  */
    return TWO127*x;
}

Here is the call graph for this function:

Here is the caller graph for this function:

long double __ieee754_exp2l ( long  double)

Definition at line 5 of file e_exp2l.c.

{
  /* This is a very stupid and inprecise implementation.  It'll get
     replaced sometime (soon?).  */
  return __ieee754_expl (M_LN2l * x);
}

Here is the call graph for this function:

Here is the caller graph for this function:

float __ieee754_expf ( float  )

Definition at line 67 of file e_expf.c.

{
  static const float himark = 88.72283935546875;
  static const float lomark = -103.972084045410;
  /* Check for usual case.  */
  if (isless (x, himark) && isgreater (x, lomark))
    {
      static const float THREEp42 = 13194139533312.0;
      static const float THREEp22 = 12582912.0;
      /* 1/ln(2).  */
#undef M_1_LN2
      static const float M_1_LN2 = 1.44269502163f;
      /* ln(2) */
#undef M_LN2
      static const double M_LN2 = .6931471805599452862;

      int tval;
      double x22, t, result, dx;
      float n, delta;
      union ieee754_double ex2_u;
      fenv_t oldenv;

      feholdexcept (&oldenv);
#ifdef FE_TONEAREST
      fesetround (FE_TONEAREST);
#endif

      /* Calculate n.  */
      n = x * M_1_LN2 + THREEp22;
      n -= THREEp22;
      dx = x - n*M_LN2;

      /* Calculate t/512.  */
      t = dx + THREEp42;
      t -= THREEp42;
      dx -= t;

      /* Compute tval = t.  */
      tval = (int) (t * 512.0);

      if (t >= 0)
       delta = - __exp_deltatable[tval];
      else
       delta = __exp_deltatable[-tval];

      /* Compute ex2 = 2^n e^(t/512+delta[t]).  */
      ex2_u.d = __exp_atable[tval+177];
      ex2_u.ieee.exponent += (int) n;

      /* Approximate e^(dx+delta) - 1, using a second-degree polynomial,
        with maximum error in [-2^-10-2^-28,2^-10+2^-28]
        less than 5e-11.  */
      x22 = (0.5000000496709180453 * dx + 1.0000001192102037084) * dx + delta;

      /* Return result.  */
      fesetenv (&oldenv);

      result = x22 * ex2_u.d + ex2_u.d;
      return (float) result;
    }
  /* Exceptional cases:  */
  else if (isless (x, himark))
    {
      if (__isinff (x))
       /* e^-inf == 0, with no error.  */
       return 0;
      else
       /* Underflow */
       return TWOM100 * TWOM100;
    }
  else
    /* Return x, if x is a NaN or Inf; or overflow, otherwise.  */
    return TWO127*x;
}

Here is the call graph for this function:

Here is the caller graph for this function:

long double __ieee754_expl ( long  double)

Definition at line 6 of file e_expl.c.

{
  fputs ("__ieee754_expl not implemented\n", stderr);
  __set_errno (ENOSYS);
  return 0.0;
}

Here is the call graph for this function:

Here is the caller graph for this function:

double __ieee754_fmod ( double  ,
double   
)

Definition at line 35 of file e_fmod.c.

{
       int32_t n,hx,hy,hz,ix,iy,sx,i;
       u_int32_t lx,ly,lz;

       EXTRACT_WORDS(hx,lx,x);
       EXTRACT_WORDS(hy,ly,y);
       sx = hx&0x80000000;         /* sign of x */
       hx ^=sx;             /* |x| */
       hy &= 0x7fffffff;    /* |y| */

    /* purge off exception values */
       if((hy|ly)==0||(hx>=0x7ff00000)||  /* y=0,or x not finite */
         ((hy|((ly|-ly)>>31))>0x7ff00000))       /* or y is NaN */
           return (x*y)/(x*y);
       if(hx<=hy) {
           if((hx<hy)||(lx<ly)) return x; /* |x|<|y| return x */
           if(lx==ly) 
              return Zero[(u_int32_t)sx>>31];    /* |x|=|y| return x*0*/
       }

    /* determine ix = ilogb(x) */
       if(hx<0x00100000) {  /* subnormal x */
           if(hx==0) {
              for (ix = -1043, i=lx; i>0; i<<=1) ix -=1;
           } else {
              for (ix = -1022,i=(hx<<11); i>0; i<<=1) ix -=1;
           }
       } else ix = (hx>>20)-1023;

    /* determine iy = ilogb(y) */
       if(hy<0x00100000) {  /* subnormal y */
           if(hy==0) {
              for (iy = -1043, i=ly; i>0; i<<=1) iy -=1;
           } else {
              for (iy = -1022,i=(hy<<11); i>0; i<<=1) iy -=1;
           }
       } else iy = (hy>>20)-1023;

    /* set up {hx,lx}, {hy,ly} and align y to x */
       if(ix >= -1022) 
           hx = 0x00100000|(0x000fffff&hx);
       else {        /* subnormal x, shift x to normal */
           n = -1022-ix;
           if(n<=31) {
               hx = (hx<<n)|(lx>>(32-n));
               lx <<= n;
           } else {
              hx = lx<<(n-32);
              lx = 0;
           }
       }
       if(iy >= -1022) 
           hy = 0x00100000|(0x000fffff&hy);
       else {        /* subnormal y, shift y to normal */
           n = -1022-iy;
           if(n<=31) {
               hy = (hy<<n)|(ly>>(32-n));
               ly <<= n;
           } else {
              hy = ly<<(n-32);
              ly = 0;
           }
       }

    /* fix point fmod */
       n = ix - iy;
       while(n--) {
           hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
           if(hz<0){hx = hx+hx+(lx>>31); lx = lx+lx;}
           else {
              if((hz|lz)==0)              /* return sign(x)*0 */
                  return Zero[(u_int32_t)sx>>31];
              hx = hz+hz+(lz>>31); lx = lz+lz;
           }
       }
       hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
       if(hz>=0) {hx=hz;lx=lz;}

    /* convert back to floating value and restore the sign */
       if((hx|lx)==0)                     /* return sign(x)*0 */
           return Zero[(u_int32_t)sx>>31];       
       while(hx<0x00100000) {             /* normalize x */
           hx = hx+hx+(lx>>31); lx = lx+lx;
           iy -= 1;
       }
       if(iy>= -1022) {     /* normalize output */
           hx = ((hx-0x00100000)|((iy+1023)<<20));
           INSERT_WORDS(x,hx|sx,lx);
       } else {             /* subnormal output */
           n = -1022 - iy;
           if(n<=20) {
              lx = (lx>>n)|((u_int32_t)hx<<(32-n));
              hx >>= n;
           } else if (n<=31) {
              lx = (hx<<(32-n))|(lx>>n); hx = sx;
           } else {
              lx = hx>>(n-32); hx = sx;
           }
           INSERT_WORDS(x,hx|sx,lx);
           x *= one;        /* create necessary signal */
       }
       return x;            /* exact output */
}

Here is the caller graph for this function:

float __ieee754_fmodf ( float  ,
float   
)

Definition at line 38 of file e_fmodf.c.

{
       int32_t n,hx,hy,hz,ix,iy,sx,i;

       GET_FLOAT_WORD(hx,x);
       GET_FLOAT_WORD(hy,y);
       sx = hx&0x80000000;         /* sign of x */
       hx ^=sx;             /* |x| */
       hy &= 0x7fffffff;    /* |y| */

    /* purge off exception values */
       if(hy==0||(hx>=0x7f800000)||              /* y=0,or x not finite */
          (hy>0x7f800000))                /* or y is NaN */
           return (x*y)/(x*y);
       if(hx<hy) return x;                /* |x|<|y| return x */
       if(hx==hy)
           return Zero[(u_int32_t)sx>>31];       /* |x|=|y| return x*0*/

    /* determine ix = ilogb(x) */
       if(hx<0x00800000) {  /* subnormal x */
           for (ix = -126,i=(hx<<8); i>0; i<<=1) ix -=1;
       } else ix = (hx>>23)-127;

    /* determine iy = ilogb(y) */
       if(hy<0x00800000) {  /* subnormal y */
           for (iy = -126,i=(hy<<8); i>=0; i<<=1) iy -=1;
       } else iy = (hy>>23)-127;

    /* set up {hx,lx}, {hy,ly} and align y to x */
       if(ix >= -126) 
           hx = 0x00800000|(0x007fffff&hx);
       else {        /* subnormal x, shift x to normal */
           n = -126-ix;
           hx = hx<<n;
       }
       if(iy >= -126) 
           hy = 0x00800000|(0x007fffff&hy);
       else {        /* subnormal y, shift y to normal */
           n = -126-iy;
           hy = hy<<n;
       }

    /* fix point fmod */
       n = ix - iy;
       while(n--) {
           hz=hx-hy;
           if(hz<0){hx = hx+hx;}
           else {
              if(hz==0)            /* return sign(x)*0 */
                  return Zero[(u_int32_t)sx>>31];
              hx = hz+hz;
           }
       }
       hz=hx-hy;
       if(hz>=0) {hx=hz;}

    /* convert back to floating value and restore the sign */
       if(hx==0)                   /* return sign(x)*0 */
           return Zero[(u_int32_t)sx>>31];       
       while(hx<0x00800000) {             /* normalize x */
           hx = hx+hx;
           iy -= 1;
       }
       if(iy>= -126) {             /* normalize output */
           hx = ((hx-0x00800000)|((iy+127)<<23));
           SET_FLOAT_WORD(x,hx|sx);
       } else {             /* subnormal output */
           n = -126 - iy;
           hx >>= n;
           SET_FLOAT_WORD(x,hx|sx);
           x *= one;        /* create necessary signal */
       }
       return x;            /* exact output */
}

Here is the caller graph for this function:

long double __ieee754_fmodl ( long  double,
long  double 
)

Definition at line 6 of file e_fmodl.c.

{
  fputs ("__ieee754_fmodl not implemented\n", stderr);
  __set_errno (ENOSYS);
  return 0.0;
}

Here is the call graph for this function:

Here is the caller graph for this function:

double __ieee754_gamma ( double  )
double __ieee754_gamma_r ( double  ,
int  
)

Definition at line 26 of file e_gamma_r.c.

{
  /* We don't have a real gamma implementation now.  We'll use lgamma
     and the exp function.  But due to the required boundary
     conditions we must check some values separately.  */
  int32_t hx;
  u_int32_t lx;

  EXTRACT_WORDS (hx, lx, x);

  if (((hx & 0x7fffffff) | lx) == 0)
    {
      /* Return value for x == 0 is Inf with divide by zero exception.  */
      *signgamp = 0;
      return 1.0 / x;
    }
  if (hx < 0 && (u_int32_t) hx < 0xfff00000 && __rint (x) == x)
    {
      /* Return value for integer x < 0 is NaN with invalid exception.  */
      *signgamp = 0;
      return (x - x) / (x - x);
    }
  if ((unsigned int) hx == 0xfff00000 && lx==0)
    {
      /* x == -Inf.  According to ISO this is NaN.  */
      *signgamp = 0;
      return x - x;
    }

  /* XXX FIXME.  */
  return __ieee754_exp (__ieee754_lgamma_r (x, signgamp));
}

Here is the call graph for this function:

Here is the caller graph for this function:

float __ieee754_gammaf ( float  )
float __ieee754_gammaf_r ( float  ,
int  
)

Definition at line 26 of file e_gammaf_r.c.

{
  /* We don't have a real gamma implementation now.  We'll use lgamma
     and the exp function.  But due to the required boundary
     conditions we must check some values separately.  */
  int32_t hx;

  GET_FLOAT_WORD (hx, x);

  if ((hx & 0x7fffffff) == 0)
    {
      /* Return value for x == 0 is Inf with divide by zero exception.  */
      *signgamp = 0;
      return 1.0 / x;
    }
  if (hx < 0 && (u_int32_t) hx < 0xff800000 && __rintf (x) == x)
    {
      /* Return value for integer x < 0 is NaN with invalid exception.  */
      *signgamp = 0;
      return (x - x) / (x - x);
    }
  if (hx == 0xff800000)
    {
      /* x == -Inf.  According to ISO this is NaN.  */
      *signgamp = 0;
      return x - x;
    }

  /* XXX FIXME.  */
  return __ieee754_expf (__ieee754_lgammaf_r (x, signgamp));
}

Here is the call graph for this function:

Here is the caller graph for this function:

long double __ieee754_gammal ( long  double)
long double __ieee754_gammal_r ( long  double,
int  
)

Definition at line 6 of file e_gammal_r.c.

{
  *signgamp = 0;
  fputs ("__ieee754_gammal_r not implemented\n", stderr);
  __set_errno (ENOSYS);
  return 0.0;
}

Here is the call graph for this function:

Here is the caller graph for this function:

double __ieee754_hypot ( double  ,
double   
)

Definition at line 55 of file e_hypot.c.

{
       double a,b,t1,t2,y1,y2,w;
       int32_t j,k,ha,hb;

       GET_HIGH_WORD(ha,x);
       ha &= 0x7fffffff;
       GET_HIGH_WORD(hb,y);
       hb &= 0x7fffffff;
       if(hb > ha) {a=y;b=x;j=ha; ha=hb;hb=j;} else {a=x;b=y;}
       SET_HIGH_WORD(a,ha); /* a <- |a| */
       SET_HIGH_WORD(b,hb); /* b <- |b| */
       if((ha-hb)>0x3c00000) {return a+b;} /* x/y > 2**60 */
       k=0;
       if(ha > 0x5f300000) {       /* a>2**500 */
          if(ha >= 0x7ff00000) {   /* Inf or NaN */
              u_int32_t low;
              w = a+b;                    /* for sNaN */
              GET_LOW_WORD(low,a);
              if(((ha&0xfffff)|low)==0) w = a;
              GET_LOW_WORD(low,b);
              if(((hb^0x7ff00000)|low)==0) w = b;
              return w;
          }
          /* scale a and b by 2**-600 */
          ha -= 0x25800000; hb -= 0x25800000;    k += 600;
          SET_HIGH_WORD(a,ha);
          SET_HIGH_WORD(b,hb);
       }
       if(hb < 0x20b00000) {       /* b < 2**-500 */
           if(hb <= 0x000fffff) {  /* subnormal b or 0 */
               u_int32_t low;
              GET_LOW_WORD(low,b);
              if((hb|low)==0) return a;
              t1=0;
              SET_HIGH_WORD(t1,0x7fd00000);      /* t1=2^1022 */
              b *= t1;
              a *= t1;
              k -= 1022;
           } else {         /* scale a and b by 2^600 */
               ha += 0x25800000;   /* a *= 2^600 */
              hb += 0x25800000;    /* b *= 2^600 */
              k -= 600;
              SET_HIGH_WORD(a,ha);
              SET_HIGH_WORD(b,hb);
           }
       }
    /* medium size a and b */
       w = a-b;
       if (w>b) {
           t1 = 0;
           SET_HIGH_WORD(t1,ha);
           t2 = a-t1;
           w  = __ieee754_sqrt(t1*t1-(b*(-b)-t2*(a+t1)));
       } else {
           a  = a+a;
           y1 = 0;
           SET_HIGH_WORD(y1,hb);
           y2 = b - y1;
           t1 = 0;
           SET_HIGH_WORD(t1,ha+0x00100000);
           t2 = a - t1;
           w  = __ieee754_sqrt(t1*y1-(w*(-w)-(t1*y2+t2*b)));
       }
       if(k!=0) {
           u_int32_t high;
           t1 = 1.0;
           GET_HIGH_WORD(high,t1);
           SET_HIGH_WORD(t1,high+(k<<20));
           return t1*w;
       } else return w;
}

Here is the call graph for this function:

Here is the caller graph for this function:

float __ieee754_hypotf ( float  ,
float   
)

Definition at line 26 of file e_hypotf.c.

{
       float a,b,t1,t2,y1,y2,w;
       int32_t j,k,ha,hb;

       GET_FLOAT_WORD(ha,x);
       ha &= 0x7fffffff;
       GET_FLOAT_WORD(hb,y);
       hb &= 0x7fffffff;
       if(hb > ha) {a=y;b=x;j=ha; ha=hb;hb=j;} else {a=x;b=y;}
       SET_FLOAT_WORD(a,ha);       /* a <- |a| */
       SET_FLOAT_WORD(b,hb);       /* b <- |b| */
       if((ha-hb)>0xf000000) {return a+b;} /* x/y > 2**30 */
       k=0;
       if(ha > 0x58800000) {       /* a>2**50 */
          if(ha >= 0x7f800000) {   /* Inf or NaN */
              w = a+b;                    /* for sNaN */
              if(ha == 0x7f800000) w = a;
              if(hb == 0x7f800000) w = b;
              return w;
          }
          /* scale a and b by 2**-60 */
          ha -= 0x1e000000; hb -= 0x1e000000;    k += 60;
          SET_FLOAT_WORD(a,ha);
          SET_FLOAT_WORD(b,hb);
       }
       if(hb < 0x26800000) {       /* b < 2**-50 */
           if(hb <= 0x007fffff) {  /* subnormal b or 0 */
               if(hb==0) return a;
              SET_FLOAT_WORD(t1,0x7e800000);     /* t1=2^126 */
              b *= t1;
              a *= t1;
              k -= 126;
           } else {         /* scale a and b by 2^60 */
               ha += 0x1e000000;   /* a *= 2^60 */
              hb += 0x1e000000;    /* b *= 2^60 */
              k -= 60;
              SET_FLOAT_WORD(a,ha);
              SET_FLOAT_WORD(b,hb);
           }
       }
    /* medium size a and b */
       w = a-b;
       if (w>b) {
           SET_FLOAT_WORD(t1,ha&0xfffff000);
           t2 = a-t1;
           w  = __ieee754_sqrtf(t1*t1-(b*(-b)-t2*(a+t1)));
       } else {
           a  = a+a;
           SET_FLOAT_WORD(y1,hb&0xfffff000);
           y2 = b - y1;
           SET_FLOAT_WORD(t1,ha+0x00800000);
           t2 = a - t1;
           w  = __ieee754_sqrtf(t1*y1-(w*(-w)-(t1*y2+t2*b)));
       }
       if(k!=0) {
           SET_FLOAT_WORD(t1,0x3f800000+(k<<23));
           return t1*w;
       } else return w;
}

Here is the call graph for this function:

Here is the caller graph for this function:

long double __ieee754_hypotl ( long  double,
long  double 
)

Definition at line 6 of file e_hypotl.c.

{
  fputs ("__ieee754_hypotl not implemented\n", stderr);
  __set_errno (ENOSYS);
  return 0.0;
}

Here is the call graph for this function:

Here is the caller graph for this function:

double __ieee754_j0 ( double  )

Definition at line 102 of file e_j0.c.

{
       double z, s,c,ss,cc,r,u,v,r1,r2,s1,s2,z2,z4;
       int32_t hx,ix;

       GET_HIGH_WORD(hx,x);
       ix = hx&0x7fffffff;
       if(ix>=0x7ff00000) return one/(x*x);
       x = fabs(x);
       if(ix >= 0x40000000) {      /* |x| >= 2.0 */
              __sincos (x, &s, &c);
              ss = s-c;
              cc = s+c;
              if(ix<0x7fe00000) {  /* make sure x+x not overflow */
                  z = -__cos(x+x);
                  if ((s*c)<zero) cc = z/ss;
                  else          ss = z/cc;
              }
       /*
        * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
        * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
        */
              if(ix>0x48000000) z = (invsqrtpi*cc)/__ieee754_sqrt(x);
              else {
                  u = pzero(x); v = qzero(x);
                  z = invsqrtpi*(u*cc-v*ss)/__ieee754_sqrt(x);
              }
              return z;
       }
       if(ix<0x3f200000) {  /* |x| < 2**-13 */
           if(huge+x>one) { /* raise inexact if x != 0 */
               if(ix<0x3e400000) return one;     /* |x|<2**-27 */
               else        return one - 0.25*x*x;
           }
       }
       z = x*x;
#ifdef DO_NOT_USE_THIS
       r =  z*(R02+z*(R03+z*(R04+z*R05)));
       s =  one+z*(S01+z*(S02+z*(S03+z*S04)));
#else
       r1 = z*R[2]; z2=z*z;
       r2 = R[3]+z*R[4]; z4=z2*z2;
       r  = r1 + z2*r2 + z4*R[5];
       s1 = one+z*S[1];
       s2 = S[2]+z*S[3];
       s = s1 + z2*s2 + z4*S[4];
#endif
       if(ix < 0x3FF00000) {       /* |x| < 1.00 */
           return one + z*(-0.25+(r/s));
       } else {
           u = 0.5*x;
           return((one+u)*(one-u)+z*(r/s));
       }
}

Here is the call graph for this function:

Here is the caller graph for this function:

float __ieee754_j0f ( float  )

Definition at line 57 of file e_j0f.c.

{
       float z, s,c,ss,cc,r,u,v;
       int32_t hx,ix;

       GET_FLOAT_WORD(hx,x);
       ix = hx&0x7fffffff;
       if(ix>=0x7f800000) return one/(x*x);
       x = fabsf(x);
       if(ix >= 0x40000000) {      /* |x| >= 2.0 */
              __sincosf (x, &s, &c);
              ss = s-c;
              cc = s+c;
              if(ix<0x7f000000) {  /* make sure x+x not overflow */
                  z = -__cosf(x+x);
                  if ((s*c)<zero) cc = z/ss;
                  else          ss = z/cc;
              }
       /*
        * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
        * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
        */
              if(ix>0x48000000) z = (invsqrtpi*cc)/__ieee754_sqrtf(x);
              else {
                  u = pzerof(x); v = qzerof(x);
                  z = invsqrtpi*(u*cc-v*ss)/__ieee754_sqrtf(x);
              }
              return z;
       }
       if(ix<0x39000000) {  /* |x| < 2**-13 */
           if(huge+x>one) { /* raise inexact if x != 0 */
               if(ix<0x32000000) return one;     /* |x|<2**-27 */
               else        return one - (float)0.25*x*x;
           }
       }
       z = x*x;
       r =  z*(R02+z*(R03+z*(R04+z*R05)));
       s =  one+z*(S01+z*(S02+z*(S03+z*S04)));
       if(ix < 0x3F800000) {       /* |x| < 1.00 */
           return one + z*((float)-0.25+(r/s));
       } else {
           u = (float)0.5*x;
           return((one+u)*(one-u)+z*(r/s));
       }
}

Here is the call graph for this function:

Here is the caller graph for this function:

long double __ieee754_j0l ( long  double)

Definition at line 7 of file e_j0l.c.

{
  fputs ("__ieee754_j0l not implemented\n", stderr);
  __set_errno (ENOSYS);
  return 0.0;
}

Here is the call graph for this function:

Here is the caller graph for this function:

double __ieee754_j1 ( double  )

Definition at line 103 of file e_j1.c.

{
       double z, s,c,ss,cc,r,u,v,y,r1,r2,s1,s2,s3,z2,z4;
       int32_t hx,ix;

       GET_HIGH_WORD(hx,x);
       ix = hx&0x7fffffff;
       if(ix>=0x7ff00000) return one/x;
       y = fabs(x);
       if(ix >= 0x40000000) {      /* |x| >= 2.0 */
              __sincos (y, &s, &c);
              ss = -s-c;
              cc = s-c;
              if(ix<0x7fe00000) {  /* make sure y+y not overflow */
                  z = __cos(y+y);
                  if ((s*c)>zero) cc = z/ss;
                  else          ss = z/cc;
              }
       /*
        * j1(x) = 1/sqrt(pi) * (P(1,x)*cc - Q(1,x)*ss) / sqrt(x)
        * y1(x) = 1/sqrt(pi) * (P(1,x)*ss + Q(1,x)*cc) / sqrt(x)
        */
              if(ix>0x48000000) z = (invsqrtpi*cc)/__ieee754_sqrt(y);
              else {
                  u = pone(y); v = qone(y);
                  z = invsqrtpi*(u*cc-v*ss)/__ieee754_sqrt(y);
              }
              if(hx<0) return -z;
              else    return  z;
       }
       if(ix<0x3e400000) {  /* |x|<2**-27 */
           if(huge+x>one) return 0.5*x;/* inexact if x!=0 necessary */
       }
       z = x*x;
#ifdef DO_NOT_USE_THIS
       r =  z*(r00+z*(r01+z*(r02+z*r03)));
       s =  one+z*(s01+z*(s02+z*(s03+z*(s04+z*s05))));
       r *= x;
#else
       r1 = z*R[0]; z2=z*z;
       r2 = R[1]+z*R[2]; z4=z2*z2;
       r = r1 + z2*r2 + z4*R[3];
       r *= x;
       s1 = one+z*S[1];
       s2 = S[2]+z*S[3];
       s3 = S[4]+z*S[5];
       s = s1 + z2*s2 + z4*s3;
#endif
       return(x*0.5+r/s);
}

Here is the call graph for this function:

Here is the caller graph for this function:

float __ieee754_j1f ( float  )

Definition at line 58 of file e_j1f.c.

{
       float z, s,c,ss,cc,r,u,v,y;
       int32_t hx,ix;

       GET_FLOAT_WORD(hx,x);
       ix = hx&0x7fffffff;
       if(ix>=0x7f800000) return one/x;
       y = fabsf(x);
       if(ix >= 0x40000000) {      /* |x| >= 2.0 */
              __sincosf (y, &s, &c);
              ss = -s-c;
              cc = s-c;
              if(ix<0x7f000000) {  /* make sure y+y not overflow */
                  z = __cosf(y+y);
                  if ((s*c)>zero) cc = z/ss;
                  else          ss = z/cc;
              }
       /*
        * j1(x) = 1/sqrt(pi) * (P(1,x)*cc - Q(1,x)*ss) / sqrt(x)
        * y1(x) = 1/sqrt(pi) * (P(1,x)*ss + Q(1,x)*cc) / sqrt(x)
        */
              if(ix>0x48000000) z = (invsqrtpi*cc)/__ieee754_sqrtf(y);
              else {
                  u = ponef(y); v = qonef(y);
                  z = invsqrtpi*(u*cc-v*ss)/__ieee754_sqrtf(y);
              }
              if(hx<0) return -z;
              else    return  z;
       }
       if(ix<0x32000000) {  /* |x|<2**-27 */
           if(huge+x>one) return (float)0.5*x;/* inexact if x!=0 necessary */
       }
       z = x*x;
       r =  z*(r00+z*(r01+z*(r02+z*r03)));
       s =  one+z*(s01+z*(s02+z*(s03+z*(s04+z*s05))));
       r *= x;
       return(x*(float)0.5+r/s);
}

Here is the call graph for this function:

Here is the caller graph for this function:

long double __ieee754_j1l ( long  double)

Definition at line 7 of file e_j1l.c.

{
  fputs ("__ieee754_j1l not implemented\n", stderr);
  __set_errno (ENOSYS);
  return 0.0;
}

Here is the call graph for this function:

Here is the caller graph for this function:

double __ieee754_jn ( int  ,
double   
)

Definition at line 64 of file e_jn.c.

{
       int32_t i,hx,ix,lx, sgn;
       double a, b, temp, di;
       double z, w;

    /* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
     * Thus, J(-n,x) = J(n,-x)
     */
       EXTRACT_WORDS(hx,lx,x);
       ix = 0x7fffffff&hx;
    /* if J(n,NaN) is NaN */
       if((ix|((u_int32_t)(lx|-lx))>>31)>0x7ff00000) return x+x;
       if(n<0){
              n = -n;
              x = -x;
              hx ^= 0x80000000;
       }
       if(n==0) return(__ieee754_j0(x));
       if(n==1) return(__ieee754_j1(x));
       sgn = (n&1)&(hx>>31);       /* even n -- 0, odd n -- sign(x) */
       x = fabs(x);
       if((ix|lx)==0||ix>=0x7ff00000)     /* if x is 0 or inf */
           b = zero;
       else if((double)n<=x) {
              /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
           if(ix>=0x52D00000) { /* x > 2**302 */
    /* (x >> n**2)
     *     Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
     *     Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
     *     Let s=sin(x), c=cos(x),
     *        xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
     *
     *           n   sin(xn)*sqt2  cos(xn)*sqt2
     *        ----------------------------------
     *           0    s-c           c+s
     *           1   -s-c          -c+s
     *           2   -s+c          -c-s
     *           3    s+c           c-s
     */
              double s;
              double c;
              __sincos (x, &s, &c);
              switch(n&3) {
                  case 0: temp =  c + s; break;
                  case 1: temp = -c + s; break;
                  case 2: temp = -c - s; break;
                  case 3: temp =  c - s; break;
              }
              b = invsqrtpi*temp/__ieee754_sqrt(x);
           } else {
               a = __ieee754_j0(x);
               b = __ieee754_j1(x);
               for(i=1;i<n;i++){
                  temp = b;
                  b = b*((double)(i+i)/x) - a; /* avoid underflow */
                  a = temp;
               }
           }
       } else {
           if(ix<0x3e100000) {     /* x < 2**-29 */
    /* x is tiny, return the first Taylor expansion of J(n,x)
     * J(n,x) = 1/n!*(x/2)^n  - ...
     */
              if(n>33)      /* underflow */
                  b = zero;
              else {
                  temp = x*0.5; b = temp;
                  for (a=one,i=2;i<=n;i++) {
                     a *= (double)i;             /* a = n! */
                     b *= temp;           /* b = (x/2)^n */
                  }
                  b = b/a;
              }
           } else {
              /* use backward recurrence */
              /*                   x      x^2      x^2
               *  J(n,x)/J(n-1,x) =  ----   ------   ------   .....
               *                   2n  - 2(n+1) - 2(n+2)
               *
               *                   1      1        1
               *  (for large x)   =  ----  ------   ------   .....
               *                   2n   2(n+1)   2(n+2)
               *                   -- - ------ - ------ -
               *                    x     x         x
               *
               * Let w = 2n/x and h=2/x, then the above quotient
               * is equal to the continued fraction:
               *                1
               *     = -----------------------
               *                   1
               *        w - -----------------
               *                     1
               *             w+h - ---------
               *                   w+2h - ...
               *
               * To determine how many terms needed, let
               * Q(0) = w, Q(1) = w(w+h) - 1,
               * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
               * When Q(k) > 1e4   good for single
               * When Q(k) > 1e9   good for double
               * When Q(k) > 1e17  good for quadruple
               */
           /* determine k */
              double t,v;
              double q0,q1,h,tmp; int32_t k,m;
              w  = (n+n)/(double)x; h = 2.0/(double)x;
              q0 = w;  z = w+h; q1 = w*z - 1.0; k=1;
              while(q1<1.0e9) {
                     k += 1; z += h;
                     tmp = z*q1 - q0;
                     q0 = q1;
                     q1 = tmp;
              }
              m = n+n;
              for(t=zero, i = 2*(n+k); i>=m; i -= 2) t = one/(i/x-t);
              a = t;
              b = one;
              /*  estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
               *  Hence, if n*(log(2n/x)) > ...
               *  single 8.8722839355e+01
               *  double 7.09782712893383973096e+02
               *  long double 1.1356523406294143949491931077970765006170e+04
               *  then recurrent value may overflow and the result is
               *  likely underflow to zero
               */
              tmp = n;
              v = two/x;
              tmp = tmp*__ieee754_log(fabs(v*tmp));
              if(tmp<7.09782712893383973096e+02) {
                  for(i=n-1,di=(double)(i+i);i>0;i--){
                      temp = b;
                     b *= di;
                     b  = b/x - a;
                      a = temp;
                     di -= two;
                  }
              } else {
                  for(i=n-1,di=(double)(i+i);i>0;i--){
                      temp = b;
                     b *= di;
                     b  = b/x - a;
                      a = temp;
                     di -= two;
                  /* scale b to avoid spurious overflow */
                     if(b>1e100) {
                         a /= b;
                         t /= b;
                         b  = one;
                     }
                  }
              }
              b = (t*__ieee754_j0(x)/b);
           }
       }
       if(sgn==1) return -b; else return b;
}

Here is the call graph for this function:

Here is the caller graph for this function:

float __ieee754_jnf ( int  ,
float   
)

Definition at line 40 of file e_jnf.c.

{
       int32_t i,hx,ix, sgn;
       float a, b, temp, di;
       float z, w;

    /* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
     * Thus, J(-n,x) = J(n,-x)
     */
       GET_FLOAT_WORD(hx,x);
       ix = 0x7fffffff&hx;
    /* if J(n,NaN) is NaN */
       if(ix>0x7f800000) return x+x;
       if(n<0){
              n = -n;
              x = -x;
              hx ^= 0x80000000;
       }
       if(n==0) return(__ieee754_j0f(x));
       if(n==1) return(__ieee754_j1f(x));
       sgn = (n&1)&(hx>>31);       /* even n -- 0, odd n -- sign(x) */
       x = fabsf(x);
       if(ix==0||ix>=0x7f800000)   /* if x is 0 or inf */
           b = zero;
       else if((float)n<=x) {
              /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
           a = __ieee754_j0f(x);
           b = __ieee754_j1f(x);
           for(i=1;i<n;i++){
              temp = b;
              b = b*((float)(i+i)/x) - a; /* avoid underflow */
              a = temp;
           }
       } else {
           if(ix<0x30800000) {     /* x < 2**-29 */
    /* x is tiny, return the first Taylor expansion of J(n,x)
     * J(n,x) = 1/n!*(x/2)^n  - ...
     */
              if(n>33)      /* underflow */
                  b = zero;
              else {
                  temp = x*(float)0.5; b = temp;
                  for (a=one,i=2;i<=n;i++) {
                     a *= (float)i;              /* a = n! */
                     b *= temp;           /* b = (x/2)^n */
                  }
                  b = b/a;
              }
           } else {
              /* use backward recurrence */
              /*                   x      x^2      x^2
               *  J(n,x)/J(n-1,x) =  ----   ------   ------   .....
               *                   2n  - 2(n+1) - 2(n+2)
               *
               *                   1      1        1
               *  (for large x)   =  ----  ------   ------   .....
               *                   2n   2(n+1)   2(n+2)
               *                   -- - ------ - ------ -
               *                    x     x         x
               *
               * Let w = 2n/x and h=2/x, then the above quotient
               * is equal to the continued fraction:
               *                1
               *     = -----------------------
               *                   1
               *        w - -----------------
               *                     1
               *             w+h - ---------
               *                   w+2h - ...
               *
               * To determine how many terms needed, let
               * Q(0) = w, Q(1) = w(w+h) - 1,
               * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
               * When Q(k) > 1e4   good for single
               * When Q(k) > 1e9   good for double
               * When Q(k) > 1e17  good for quadruple
               */
           /* determine k */
              float t,v;
              float q0,q1,h,tmp; int32_t k,m;
              w  = (n+n)/(float)x; h = (float)2.0/(float)x;
              q0 = w;  z = w+h; q1 = w*z - (float)1.0; k=1;
              while(q1<(float)1.0e9) {
                     k += 1; z += h;
                     tmp = z*q1 - q0;
                     q0 = q1;
                     q1 = tmp;
              }
              m = n+n;
              for(t=zero, i = 2*(n+k); i>=m; i -= 2) t = one/(i/x-t);
              a = t;
              b = one;
              /*  estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
               *  Hence, if n*(log(2n/x)) > ...
               *  single 8.8722839355e+01
               *  double 7.09782712893383973096e+02
               *  long double 1.1356523406294143949491931077970765006170e+04
               *  then recurrent value may overflow and the result is
               *  likely underflow to zero
               */
              tmp = n;
              v = two/x;
              tmp = tmp*__ieee754_logf(fabsf(v*tmp));
              if(tmp<(float)8.8721679688e+01) {
                  for(i=n-1,di=(float)(i+i);i>0;i--){
                      temp = b;
                     b *= di;
                     b  = b/x - a;
                      a = temp;
                     di -= two;
                  }
              } else {
                  for(i=n-1,di=(float)(i+i);i>0;i--){
                      temp = b;
                     b *= di;
                     b  = b/x - a;
                      a = temp;
                     di -= two;
                  /* scale b to avoid spurious overflow */
                     if(b>(float)1e10) {
                         a /= b;
                         t /= b;
                         b  = one;
                     }
                  }
              }
              b = (t*__ieee754_j0f(x)/b);
           }
       }
       if(sgn==1) return -b; else return b;
}

Here is the call graph for this function:

Here is the caller graph for this function:

long double __ieee754_jnl ( int  ,
long  double 
)

Definition at line 7 of file e_jnl.c.

{
  fputs ("__ieee754_jnl not implemented\n", stderr);
  __set_errno (ENOSYS);
  return 0.0;
}

Here is the call graph for this function:

Here is the caller graph for this function:

double __ieee754_lgamma ( double  )

Definition at line 60 of file w_lgamma.c.

{
#ifdef __POSIX__
    extern int    signgam;
#else
    int    signgam;
#endif
    return __libm_lgamma(x, &signgam, sizeof(signgam));
}

Here is the call graph for this function:

double __ieee754_lgamma_r ( double  ,
int  
)

Definition at line 59 of file e_lgamma_r.c.

{
    return __libm_lgamma(x, signgam, sizeof(*signgam));
}

Here is the call graph for this function:

Here is the caller graph for this function:

float __ieee754_lgammaf ( float  )

Definition at line 60 of file w_lgammaf.c.

{
#ifdef __POSIX__
    extern int    signgam;
#else
    int    signgam;
#endif
    return __libm_lgammaf(x, &signgam, sizeof(signgam));
}

Here is the call graph for this function:

float __ieee754_lgammaf_r ( float  ,
int  
)

Definition at line 59 of file e_lgammaf_r.c.

{
    return __libm_lgammaf(x, signgam, sizeof(*signgam));
}

Here is the call graph for this function:

Here is the caller graph for this function:

long double __ieee754_lgammal ( long  double)

Definition at line 59 of file w_lgammal.c.

{
#ifdef __POSIX__
    extern int    signgam;
#else
    int    signgam;
#endif
    return __libm_lgammal(x, &signgam, sizeof(signgam));
}

Here is the call graph for this function:

long double __ieee754_lgammal_r ( long  double,
int  
)

Definition at line 7 of file e_lgammal_r.c.

{
  *signgamp = 0;
  fputs ("__ieee754_lgammal_r not implemented\n", stderr);
  __set_errno (ENOSYS);
  return 0.0;
}

Here is the call graph for this function:

Here is the caller graph for this function:

double __ieee754_log ( double  )

Definition at line 50 of file e_log.c.

                               {
#define M 4
  static const int pr[M]={8,10,18,32};
  int i,j,n,ux,dx,p;
#if 0
  int k;
#endif
  double dbl_n,u,p0,q,r0,w,nln2a,luai,lubi,lvaj,lvbj,
         sij,ssij,ttij,A,B,B0,y,y1,y2,polI,polII,sa,sb,
         t1,t2,t3,t4,t5,t6,t7,t8,t,ra,rb,ww,
         a0,aa0,s1,s2,ss2,s3,ss3,a1,aa1,a,aa,b,bb,c;
  number num;
  mp_no mpx,mpy,mpy1,mpy2,mperr;

#include "ulog.tbl"
#include "ulog.h"

  /* Treating special values of x ( x<=0, x=INF, x=NaN etc.). */

  num.d = x;  ux = num.i[HIGH_HALF];  dx = num.i[LOW_HALF];
  n=0;
  if (ux < 0x00100000) {
    if (((ux & 0x7fffffff) | dx) == 0)  return MHALF/ZERO; /* return -INF */
    if (ux < 0) return (x-x)/ZERO;                         /* return NaN  */
    n -= 54;    x *= two54.d;                              /* scale x     */
    num.d = x;
  }
  if (ux >= 0x7ff00000) return x+x;                        /* INF or NaN  */

  /* Regular values of x */

  w = x-ONE;
  if (ABS(w) > U03) { goto case_03; }


  /*--- Stage I, the case abs(x-1) < 0.03 */

  t8 = MHALF*w;
  EMULV(t8,w,a,aa,t1,t2,t3,t4,t5)
  EADD(w,a,b,bb)

  /* Evaluate polynomial II */
  polII = (b0.d+w*(b1.d+w*(b2.d+w*(b3.d+w*(b4.d+
          w*(b5.d+w*(b6.d+w*(b7.d+w*b8.d))))))))*w*w*w;
  c = (aa+bb)+polII;

  /* End stage I, case abs(x-1) < 0.03 */
  if ((y=b+(c+b*E2)) == b+(c-b*E2))  return y;

  /*--- Stage II, the case abs(x-1) < 0.03 */

  a = d11.d+w*(d12.d+w*(d13.d+w*(d14.d+w*(d15.d+w*(d16.d+
            w*(d17.d+w*(d18.d+w*(d19.d+w*d20.d))))))));
  EMULV(w,a,s2,ss2,t1,t2,t3,t4,t5)
  ADD2(d10.d,dd10.d,s2,ss2,s3,ss3,t1,t2)
  MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
  ADD2(d9.d,dd9.d,s2,ss2,s3,ss3,t1,t2)
  MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
  ADD2(d8.d,dd8.d,s2,ss2,s3,ss3,t1,t2)
  MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
  ADD2(d7.d,dd7.d,s2,ss2,s3,ss3,t1,t2)
  MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
  ADD2(d6.d,dd6.d,s2,ss2,s3,ss3,t1,t2)
  MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
  ADD2(d5.d,dd5.d,s2,ss2,s3,ss3,t1,t2)
  MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
  ADD2(d4.d,dd4.d,s2,ss2,s3,ss3,t1,t2)
  MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
  ADD2(d3.d,dd3.d,s2,ss2,s3,ss3,t1,t2)
  MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
  ADD2(d2.d,dd2.d,s2,ss2,s3,ss3,t1,t2)
  MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
  MUL2(w,ZERO,s2,ss2,s3,ss3,t1,t2,t3,t4,t5,t6,t7,t8)
  ADD2(w,ZERO,    s3,ss3, b, bb,t1,t2)

  /* End stage II, case abs(x-1) < 0.03 */
  if ((y=b+(bb+b*E4)) == b+(bb-b*E4))  return y;
  goto stage_n;

  /*--- Stage I, the case abs(x-1) > 0.03 */
  case_03:

  /* Find n,u such that x = u*2**n,   1/sqrt(2) < u < sqrt(2)  */
  n += (num.i[HIGH_HALF] >> 20) - 1023;
  num.i[HIGH_HALF] = (num.i[HIGH_HALF] & 0x000fffff) | 0x3ff00000;
  if (num.d > SQRT_2) { num.d *= HALF;  n++; }
  u = num.d;  dbl_n = (double) n;

  /* Find i such that ui=1+(i-75)/2**8 is closest to u (i= 0,1,2,...,181) */
  num.d += h1.d;
  i = (num.i[HIGH_HALF] & 0x000fffff) >> 12;

  /* Find j such that vj=1+(j-180)/2**16 is closest to v=u/ui (j= 0,...,361) */
  num.d = u*Iu[i].d + h2.d;
  j = (num.i[HIGH_HALF] & 0x000fffff) >> 4;

  /* Compute w=(u-ui*vj)/(ui*vj) */
  p0=(ONE+(i-75)*DEL_U)*(ONE+(j-180)*DEL_V);
  q=u-p0;   r0=Iu[i].d*Iv[j].d;   w=q*r0;

  /* Evaluate polynomial I */
  polI = w+(a2.d+a3.d*w)*w*w;

  /* Add up everything */
  nln2a = dbl_n*LN2A;
  luai  = Lu[i][0].d;   lubi  = Lu[i][1].d;
  lvaj  = Lv[j][0].d;   lvbj  = Lv[j][1].d;
  EADD(luai,lvaj,sij,ssij)
  EADD(nln2a,sij,A  ,ttij)
  B0 = (((lubi+lvbj)+ssij)+ttij)+dbl_n*LN2B;
  B  = polI+B0;

  /* End stage I, case abs(x-1) >= 0.03 */
  if ((y=A+(B+E1)) == A+(B-E1))  return y;


  /*--- Stage II, the case abs(x-1) > 0.03 */

  /* Improve the accuracy of r0 */
  EMULV(p0,r0,sa,sb,t1,t2,t3,t4,t5)
  t=r0*((ONE-sa)-sb);
  EADD(r0,t,ra,rb)

  /* Compute w */
  MUL2(q,ZERO,ra,rb,w,ww,t1,t2,t3,t4,t5,t6,t7,t8)

  EADD(A,B0,a0,aa0)

  /* Evaluate polynomial III */
  s1 = (c3.d+(c4.d+c5.d*w)*w)*w;
  EADD(c2.d,s1,s2,ss2)
  MUL2(s2,ss2,w,ww,s3,ss3,t1,t2,t3,t4,t5,t6,t7,t8)
  MUL2(s3,ss3,w,ww,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
  ADD2(s2,ss2,w,ww,s3,ss3,t1,t2)
  ADD2(s3,ss3,a0,aa0,a1,aa1,t1,t2)

  /* End stage II, case abs(x-1) >= 0.03 */
  if ((y=a1+(aa1+E3)) == a1+(aa1-E3)) return y;


  /* Final stages. Use multi-precision arithmetic. */
  stage_n:

  for (i=0; i<M; i++) {
    p = pr[i];
    __dbl_mp(x,&mpx,p);  __dbl_mp(y,&mpy,p);
    __mplog(&mpx,&mpy,p);
    __dbl_mp(e[i].d,&mperr,p);
    __add(&mpy,&mperr,&mpy1,p);  __sub(&mpy,&mperr,&mpy2,p);
    __mp_dbl(&mpy1,&y1,p);       __mp_dbl(&mpy2,&y2,p);
    if (y1==y2)   return y1;
  }
  return y1;
}

Here is the call graph for this function:

Here is the caller graph for this function:

double __ieee754_log10 ( double  )

Definition at line 72 of file e_log10.c.

{
       double y,z;
       int32_t i,k,hx;
       u_int32_t lx;

       EXTRACT_WORDS(hx,lx,x);

        k=0;
        if (hx < 0x00100000) {                   /* x < 2**-1022  */
            if (((hx&0x7fffffff)|lx)==0)
                return -two54/(x-x);             /* log(+-0)=-inf */
            if (hx<0) return (x-x)/(x-x); /* log(-#) = NaN */
            k -= 54; x *= two54; /* subnormal number, scale up x */
           GET_HIGH_WORD(hx,x);
        }
       if (hx >= 0x7ff00000) return x+x;
       k += (hx>>20)-1023;
       i  = ((u_int32_t)k&0x80000000)>>31;
        hx = (hx&0x000fffff)|((0x3ff-i)<<20);
        y  = (double)(k+i);
       SET_HIGH_WORD(x,hx);
       z  = y*log10_2lo + ivln10*__ieee754_log(x);
       return  z+y*log10_2hi;
}

Here is the call graph for this function:

Here is the caller graph for this function:

float __ieee754_log10f ( float  )

Definition at line 42 of file e_log10f.c.

{
       float y,z;
       int32_t i,k,hx;

       GET_FLOAT_WORD(hx,x);

        k=0;
        if (hx < 0x00800000) {                   /* x < 2**-126  */
            if ((hx&0x7fffffff)==0)
                return -two25/(x-x);             /* log(+-0)=-inf */
            if (hx<0) return (x-x)/(x-x); /* log(-#) = NaN */
            k -= 25; x *= two25; /* subnormal number, scale up x */
           GET_FLOAT_WORD(hx,x);
        }
       if (hx >= 0x7f800000) return x+x;
       k += (hx>>23)-127;
       i  = ((u_int32_t)k&0x80000000)>>31;
        hx = (hx&0x007fffff)|((0x7f-i)<<23);
        y  = (float)(k+i);
       SET_FLOAT_WORD(x,hx);
       z  = y*log10_2lo + ivln10*__ieee754_logf(x);
       return  z+y*log10_2hi;
}

Here is the call graph for this function:

Here is the caller graph for this function:

long double __ieee754_log10l ( long  double)

Definition at line 6 of file e_log10l.c.

{
  fputs ("__ieee754_log10l not implemented\n", stderr);
  __set_errno (ENOSYS);
  return 0.0;
}

Here is the call graph for this function:

Here is the caller graph for this function:

double __ieee754_log2 ( double  )

Definition at line 84 of file e_log2.c.

{
       double hfsq,f,s,z,R,w,t1,t2,dk;
       int32_t k,hx,i,j;
       u_int32_t lx;

       EXTRACT_WORDS(hx,lx,x);

       k=0;
       if (hx < 0x00100000) {                    /* x < 2**-1022  */
           if (((hx&0x7fffffff)|lx)==0)
              return -two54/(x-x);        /* log(+-0)=-inf */
           if (hx<0) return (x-x)/(x-x);  /* log(-#) = NaN */
           k -= 54; x *= two54; /* subnormal number, scale up x */
           GET_HIGH_WORD(hx,x);
       }
       if (hx >= 0x7ff00000) return x+x;
       k += (hx>>20)-1023;
       hx &= 0x000fffff;
       i = (hx+0x95f64)&0x100000;
       SET_HIGH_WORD(x,hx|(i^0x3ff00000));       /* normalize x or x/2 */
       k += (i>>20);
       dk = (double) k;
       f = x-1.0;
       if((0x000fffff&(2+hx))<3) { /* |f| < 2**-20 */
           if(f==zero) return dk;
           R = f*f*(0.5-0.33333333333333333*f);
           return dk-(R-f)/ln2;
       }
       s = f/(2.0+f);
       z = s*s;
       i = hx-0x6147a;
       w = z*z;
       j = 0x6b851-hx;
       t1= w*(Lg2+w*(Lg4+w*Lg6));
       t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7)));
       i |= j;
       R = t2+t1;
       if(i>0) {
           hfsq=0.5*f*f;
           return dk-((hfsq-(s*(hfsq+R)))-f)/ln2;
       } else {
           return dk-((s*(f-R))-f)/ln2;
       }
}

Here is the caller graph for this function:

float __ieee754_log2f ( float  )

Definition at line 45 of file e_log2f.c.

{
       float hfsq,f,s,z,R,w,t1,t2,dk;
       int32_t k,ix,i,j;

       GET_FLOAT_WORD(ix,x);

       k=0;
       if (ix < 0x00800000) {                    /* x < 2**-126  */
           if ((ix&0x7fffffff)==0)
              return -two25/(x-x);        /* log(+-0)=-inf */
           if (ix<0) return (x-x)/(x-x);  /* log(-#) = NaN */
           k -= 25; x *= two25; /* subnormal number, scale up x */
           GET_FLOAT_WORD(ix,x);
       }
       if (ix >= 0x7f800000) return x+x;
       k += (ix>>23)-127;
       ix &= 0x007fffff;
       i = (ix+(0x95f64<<3))&0x800000;
       SET_FLOAT_WORD(x,ix|(i^0x3f800000));      /* normalize x or x/2 */
       k += (i>>23);
       dk = (float)k;
       f = x-(float)1.0;
       if((0x007fffff&(15+ix))<16) {      /* |f| < 2**-20 */
           if(f==zero) return dk;
           R = f*f*((float)0.5-(float)0.33333333333333333*f);
           return dk-(R-f)/ln2;
       }
       s = f/((float)2.0+f);
       z = s*s;
       i = ix-(0x6147a<<3);
       w = z*z;
       j = (0x6b851<<3)-ix;
       t1= w*(Lg2+w*(Lg4+w*Lg6));
       t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7)));
       i |= j;
       R = t2+t1;
       if(i>0) {
           hfsq=(float)0.5*f*f;
           return dk-((hfsq-(s*(hfsq+R)))-f)/ln2;
       } else {
           return dk-((s*(f-R))-f)/ln2;
       }
}

Here is the caller graph for this function:

long double __ieee754_log2l ( long  double)

Definition at line 6 of file e_log2l.c.

{
  fputs ("__ieee754_log2l not implemented\n", stderr);
  __set_errno (ENOSYS);
  return 0.0;
}

Here is the call graph for this function:

Here is the caller graph for this function:

float __ieee754_logf ( float  )

Definition at line 48 of file e_logf.c.

{
       float hfsq,f,s,z,R,w,t1,t2,dk;
       int32_t k,ix,i,j;

       GET_FLOAT_WORD(ix,x);

       k=0;
       if (ix < 0x00800000) {                    /* x < 2**-126  */
           if ((ix&0x7fffffff)==0)
              return -two25/(x-x);        /* log(+-0)=-inf */
           if (ix<0) return (x-x)/(x-x);  /* log(-#) = NaN */
           k -= 25; x *= two25; /* subnormal number, scale up x */
           GET_FLOAT_WORD(ix,x);
       }
       if (ix >= 0x7f800000) return x+x;
       k += (ix>>23)-127;
       ix &= 0x007fffff;
       i = (ix+(0x95f64<<3))&0x800000;
       SET_FLOAT_WORD(x,ix|(i^0x3f800000));      /* normalize x or x/2 */
       k += (i>>23);
       f = x-(float)1.0;
       if((0x007fffff&(15+ix))<16) {      /* |f| < 2**-20 */
           if(f==zero) {
             if(k==0) return zero;  else {dk=(float)k;
                                      return dk*ln2_hi+dk*ln2_lo;}
           }
           R = f*f*((float)0.5-(float)0.33333333333333333*f);
           if(k==0) return f-R; else {dk=(float)k;
                   return dk*ln2_hi-((R-dk*ln2_lo)-f);}
       }
       s = f/((float)2.0+f);
       dk = (float)k;
       z = s*s;
       i = ix-(0x6147a<<3);
       w = z*z;
       j = (0x6b851<<3)-ix;
       t1= w*(Lg2+w*(Lg4+w*Lg6));
       t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7)));
       i |= j;
       R = t2+t1;
       if(i>0) {
           hfsq=(float)0.5*f*f;
           if(k==0) return f-(hfsq-s*(hfsq+R)); else
                   return dk*ln2_hi-((hfsq-(s*(hfsq+R)+dk*ln2_lo))-f);
       } else {
           if(k==0) return f-s*(f-R); else
                   return dk*ln2_hi-((s*(f-R)-dk*ln2_lo)-f);
       }
}

Here is the caller graph for this function:

long double __ieee754_logl ( long  double)

Definition at line 6 of file e_logl.c.

{
  fputs ("__ieee754_logl not implemented\n", stderr);
  __set_errno (ENOSYS);
  return 0.0;
}

Here is the call graph for this function:

Here is the caller graph for this function:

double __ieee754_pow ( double  ,
double   
)

Definition at line 58 of file e_pow.c.

                                         {
  double z,a,aa,error, t,a1,a2,y1,y2;
#if 0
  double gor=1.0;
#endif
  mynumber u,v;
  int k;
  int4 qx,qy;
  v.x=y;
  u.x=x;
  if (v.i[LOW_HALF] == 0) { /* of y */
    qx = u.i[HIGH_HALF]&0x7fffffff;
    /* Checking  if x is not too small to compute */
    if (((qx==0x7ff00000)&&(u.i[LOW_HALF]!=0))||(qx>0x7ff00000)) return NaNQ.x;
    if (y == 1.0) return x;
    if (y == 2.0) return x*x;
    if (y == -1.0) return 1.0/x;
    if (y == 0) return 1.0;
  }
  /* else */
  if(((u.i[HIGH_HALF]>0 && u.i[HIGH_HALF]<0x7ff00000)||        /* x>0 and not x->0 */
       (u.i[HIGH_HALF]==0 && u.i[LOW_HALF]!=0))  &&
                                      /*   2^-1023< x<= 2^-1023 * 0x1.0000ffffffff */
      (v.i[HIGH_HALF]&0x7fffffff) < 0x4ff00000) {              /* if y<-1 or y>1   */
    z = log1(x,&aa,&error);                                 /* x^y  =e^(y log (X)) */
    t = y*134217729.0;
    y1 = t - (t-y);
    y2 = y - y1;
    t = z*134217729.0;
    a1 = t - (t-z);
    a2 = (z - a1)+aa;
    a = y1*a1;
    aa = y2*a1 + y*a2;
    a1 = a+aa;
    a2 = (a-a1)+aa;
    error = error*ABS(y);
    t = __exp1(a1,a2,1.9e16*error);     /* return -10 or 0 if wasn't computed exactly */
    return (t>0)?t:power1(x,y);
  }

  if (x == 0) {
    if (((v.i[HIGH_HALF] & 0x7fffffff) == 0x7ff00000 && v.i[LOW_HALF] != 0)
       || (v.i[HIGH_HALF] & 0x7fffffff) > 0x7ff00000)
      return y;
    if (ABS(y) > 1.0e20) return (y>0)?0:INF.x;
    k = checkint(y);
    if (k == -1)
      return y < 0 ? 1.0/x : x;
    else
      return y < 0 ? 1.0/ABS(x) : 0.0;                               /* return 0 */
  }

  qx = u.i[HIGH_HALF]&0x7fffffff;  /*   no sign   */
  qy = v.i[HIGH_HALF]&0x7fffffff;  /*   no sign   */

  if (qx >= 0x7ff00000 && (qx > 0x7ff00000 || u.i[LOW_HALF] != 0)) return NaNQ.x;
  if (qy >= 0x7ff00000 && (qy > 0x7ff00000 || v.i[LOW_HALF] != 0))
    return x == 1.0 ? 1.0 : NaNQ.x;

  /* if x<0 */
  if (u.i[HIGH_HALF] < 0) {
    k = checkint(y);
    if (k==0) {
      if (qy == 0x7ff00000) {
       if (x == -1.0) return 1.0;
       else if (x > -1.0) return v.i[HIGH_HALF] < 0 ? INF.x : 0.0;
       else return v.i[HIGH_HALF] < 0 ? 0.0 : INF.x;
      }
      else if (qx == 0x7ff00000)
       return y < 0 ? 0.0 : INF.x;
      return NaNQ.x;                              /* y not integer and x<0 */
    }
    else if (qx == 0x7ff00000)
      {
       if (k < 0)
         return y < 0 ? nZERO.x : nINF.x;
       else
         return y < 0 ? 0.0 : INF.x;
      }
    return (k==1)?__ieee754_pow(-x,y):-__ieee754_pow(-x,y); /* if y even or odd */
  }
  /* x>0 */

  if (qx == 0x7ff00000)                              /* x= 2^-0x3ff */
    {if (y == 0) return NaNQ.x;
    return (y>0)?x:0; }

  if (qy > 0x45f00000 && qy < 0x7ff00000) {
    if (x == 1.0) return 1.0;
    if (y>0) return (x>1.0)?INF.x:0;
    if (y<0) return (x<1.0)?INF.x:0;
  }

  if (x == 1.0) return 1.0;
  if (y>0) return (x>1.0)?INF.x:0;
  if (y<0) return (x<1.0)?INF.x:0;
  return 0;     /* unreachable, to make the compiler happy */
}

Here is the call graph for this function:

Here is the caller graph for this function:

float __ieee754_powf ( float  ,
float   
)

Definition at line 63 of file e_powf.c.

{
       float z,ax,z_h,z_l,p_h,p_l;
       float y1,t1,t2,r,s,t,u,v,w;
       int32_t i,j,k,yisint,n;
       int32_t hx,hy,ix,iy,is;

       GET_FLOAT_WORD(hx,x);
       GET_FLOAT_WORD(hy,y);
       ix = hx&0x7fffffff;  iy = hy&0x7fffffff;

    /* y==zero: x**0 = 1 */
       if(iy==0) return one;

    /* x==+-1 */
       if(x == 1.0) return one;
       if(x == -1.0 && isinf(y)) return one;

    /* +-NaN return x+y */
       if(ix > 0x7f800000 ||
          iy > 0x7f800000)
              return x+y;

    /* determine if y is an odd int when x < 0
     * yisint = 0    ... y is not an integer
     * yisint = 1    ... y is an odd int
     * yisint = 2    ... y is an even int
     */
       yisint  = 0;
       if(hx<0) {
           if(iy>=0x4b800000) yisint = 2; /* even integer y */
           else if(iy>=0x3f800000) {
              k = (iy>>23)-0x7f;      /* exponent */
              j = iy>>(23-k);
              if((j<<(23-k))==iy) yisint = 2-(j&1);
           }
       }

    /* special value of y */
       if (iy==0x7f800000) {       /* y is +-inf */
           if (ix==0x3f800000)
               return  y - y;      /* inf**+-1 is NaN */
           else if (ix > 0x3f800000)/* (|x|>1)**+-inf = inf,0 */
               return (hy>=0)? y: zero;
           else                    /* (|x|<1)**-,+inf = inf,0 */
               return (hy<0)?-y: zero;
       }
       if(iy==0x3f800000) { /* y is  +-1 */
           if(hy<0) return one/x; else return x;
       }
       if(hy==0x40000000) return x*x; /* y is  2 */
       if(hy==0x3f000000) { /* y is  0.5 */
           if(hx>=0) /* x >= +0 */
           return __ieee754_sqrtf(x);
       }

       ax   = fabsf(x);
    /* special value of x */
       if(ix==0x7f800000||ix==0||ix==0x3f800000){
           z = ax;                 /*x is +-0,+-inf,+-1*/
           if(hy<0) z = one/z;     /* z = (1/|x|) */
           if(hx<0) {
              if(((ix-0x3f800000)|yisint)==0) {
                  z = (z-z)/(z-z); /* (-1)**non-int is NaN */
              } else if(yisint==1)
                  z = -z;          /* (x<0)**odd = -(|x|**odd) */
           }
           return z;
       }

    /* (x<0)**(non-int) is NaN */
       if(((((u_int32_t)hx>>31)-1)|yisint)==0) return (x-x)/(x-x);

    /* |y| is huge */
       if(iy>0x4d000000) { /* if |y| > 2**27 */
       /* over/underflow if x is not close to one */
           if(ix<0x3f7ffff8) return (hy<0)? huge*huge:tiny*tiny;
           if(ix>0x3f800007) return (hy>0)? huge*huge:tiny*tiny;
       /* now |1-x| is tiny <= 2**-20, suffice to compute
          log(x) by x-x^2/2+x^3/3-x^4/4 */
           t = x-1;         /* t has 20 trailing zeros */
           w = (t*t)*((float)0.5-t*((float)0.333333333333-t*(float)0.25));
           u = ivln2_h*t;   /* ivln2_h has 16 sig. bits */
           v = t*ivln2_l-w*ivln2;
           t1 = u+v;
           GET_FLOAT_WORD(is,t1);
           SET_FLOAT_WORD(t1,is&0xfffff000);
           t2 = v-(t1-u);
       } else {
           float s2,s_h,s_l,t_h,t_l;
           n = 0;
       /* take care subnormal number */
           if(ix<0x00800000)
              {ax *= two24; n -= 24; GET_FLOAT_WORD(ix,ax); }
           n  += ((ix)>>23)-0x7f;
           j  = ix&0x007fffff;
       /* determine interval */
           ix = j|0x3f800000;             /* normalize ix */
           if(j<=0x1cc471) k=0;    /* |x|<sqrt(3/2) */
           else if(j<0x5db3d7) k=1;       /* |x|<sqrt(3)   */
           else {k=0;n+=1;ix -= 0x00800000;}
           SET_FLOAT_WORD(ax,ix);

       /* compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
           u = ax-bp[k];           /* bp[0]=1.0, bp[1]=1.5 */
           v = one/(ax+bp[k]);
           s = u*v;
           s_h = s;
           GET_FLOAT_WORD(is,s_h);
           SET_FLOAT_WORD(s_h,is&0xfffff000);
       /* t_h=ax+bp[k] High */
           SET_FLOAT_WORD(t_h,((ix>>1)|0x20000000)+0x0040000+(k<<21));
           t_l = ax - (t_h-bp[k]);
           s_l = v*((u-s_h*t_h)-s_h*t_l);
       /* compute log(ax) */
           s2 = s*s;
           r = s2*s2*(L1+s2*(L2+s2*(L3+s2*(L4+s2*(L5+s2*L6)))));
           r += s_l*(s_h+s);
           s2  = s_h*s_h;
           t_h = (float)3.0+s2+r;
           GET_FLOAT_WORD(is,t_h);
           SET_FLOAT_WORD(t_h,is&0xfffff000);
           t_l = r-((t_h-(float)3.0)-s2);
       /* u+v = s*(1+...) */
           u = s_h*t_h;
           v = s_l*t_h+t_l*s;
       /* 2/(3log2)*(s+...) */
           p_h = u+v;
           GET_FLOAT_WORD(is,p_h);
           SET_FLOAT_WORD(p_h,is&0xfffff000);
           p_l = v-(p_h-u);
           z_h = cp_h*p_h;         /* cp_h+cp_l = 2/(3*log2) */
           z_l = cp_l*p_h+p_l*cp+dp_l[k];
       /* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */
           t = (float)n;
           t1 = (((z_h+z_l)+dp_h[k])+t);
           GET_FLOAT_WORD(is,t1);
           SET_FLOAT_WORD(t1,is&0xfffff000);
           t2 = z_l-(((t1-t)-dp_h[k])-z_h);
       }

       s = one; /* s (sign of result -ve**odd) = -1 else = 1 */
       if(((((u_int32_t)hx>>31)-1)|(yisint-1))==0)
           s = -one; /* (-ve)**(odd int) */

    /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */
       GET_FLOAT_WORD(is,y);
       SET_FLOAT_WORD(y1,is&0xfffff000);
       p_l = (y-y1)*t1+y*t2;
       p_h = y1*t1;
       z = p_l+p_h;
       GET_FLOAT_WORD(j,z);
       if (j>0x43000000)                         /* if z > 128 */
           return s*huge*huge;                          /* overflow */
       else if (j==0x43000000) {                 /* if z == 128 */
           if(p_l+ovt>z-p_h) return s*huge*huge; /* overflow */
       }
       else if ((j&0x7fffffff)>0x43160000)              /* z <= -150 */
           return s*tiny*tiny;                          /* underflow */
       else if ((u_int32_t) j==0xc3160000){             /* z == -150 */
           if(p_l<=z-p_h) return s*tiny*tiny;           /* underflow */
       }
    /*
     * compute 2**(p_h+p_l)
     */
       i = j&0x7fffffff;
       k = (i>>23)-0x7f;
       n = 0;
       if(i>0x3f000000) {          /* if |z| > 0.5, set n = [z+0.5] */
           n = j+(0x00800000>>(k+1));
           k = ((n&0x7fffffff)>>23)-0x7f; /* new k for n */
           SET_FLOAT_WORD(t,n&~(0x007fffff>>k));
           n = ((n&0x007fffff)|0x00800000)>>(23-k);
           if(j<0) n = -n;
           p_h -= t;
       }
       t = p_l+p_h;
       GET_FLOAT_WORD(is,t);
       SET_FLOAT_WORD(t,is&0xfffff000);
       u = t*lg2_h;
       v = (p_l-(t-p_h))*lg2+t*lg2_l;
       z = u+v;
       w = v-(z-u);
       t  = z*z;
       t1  = z - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
       r  = (z*t1)/(t1-two)-(w+z*w);
       z  = one-(r-z);
       GET_FLOAT_WORD(j,z);
       j += (n<<23);
       if((j>>23)<=0) z = __scalbnf(z,n); /* subnormal output */
       else SET_FLOAT_WORD(z,j);
       return s*z;
}

Here is the call graph for this function:

Here is the caller graph for this function:

long double __ieee754_powl ( long  double,
long  double 
)

Definition at line 6 of file e_powl.c.

{
  fputs ("__ieee754_powl not implemented\n", stderr);
  __set_errno (ENOSYS);
  return 0.0;
}

Here is the call graph for this function:

Here is the caller graph for this function:

int32_t __ieee754_rem_pio2 ( double  ,
double *   
)
int32_t __ieee754_rem_pio2f ( float  ,
float *   
)

Here is the caller graph for this function:

int __ieee754_rem_pio2l ( long  double,
long double *   
)

Definition at line 7 of file e_rem_pio2l.c.

{
  fputs ("__ieee754_rem_pio2l not implemented\n", stderr);
  __set_errno (ENOSYS);
  return 0;
}

Here is the call graph for this function:

Here is the caller graph for this function:

double __ieee754_remainder ( double  ,
double   
)

Definition at line 43 of file e_remainder.c.

{
  double z,d,xx;
#if 0
  double yy;
#endif
  int4 kx,ky,n,nn,n1,m1,l;
#if 0
  int4 m;
#endif
  mynumber u,t,w={{0,0}},v={{0,0}},ww={{0,0}},r;
  u.x=x;
  t.x=y;
  kx=u.i[HIGH_HALF]&0x7fffffff; /* no sign  for x*/
  t.i[HIGH_HALF]&=0x7fffffff;   /*no sign for y */
  ky=t.i[HIGH_HALF];
  /*------ |x| < 2^1023  and   2^-970 < |y| < 2^1024 ------------------*/
  if (kx<0x7fe00000 && ky<0x7ff00000 && ky>=0x03500000) {
    if (kx+0x00100000<ky) return x;
    if ((kx-0x01500000)<ky) {
      z=x/t.x;
      v.i[HIGH_HALF]=t.i[HIGH_HALF];
      d=(z+big.x)-big.x;
      xx=(x-d*v.x)-d*(t.x-v.x);
      if (d-z!=0.5&&d-z!=-0.5) return (xx!=0)?xx:((x>0)?ZERO.x:nZERO.x);
      else {
       if (ABS(xx)>0.5*t.x) return (z>d)?xx-t.x:xx+t.x;
       else return xx;
      }
    }   /*    (kx<(ky+0x01500000))         */
    else  {
      r.x=1.0/t.x;
      n=t.i[HIGH_HALF];
      nn=(n&0x7ff00000)+0x01400000;
      w.i[HIGH_HALF]=n;
      ww.x=t.x-w.x;
      l=(kx-nn)&0xfff00000;
      n1=ww.i[HIGH_HALF];
      m1=r.i[HIGH_HALF];
      while (l>0) {
       r.i[HIGH_HALF]=m1-l;
       z=u.x*r.x;
       w.i[HIGH_HALF]=n+l;
       ww.i[HIGH_HALF]=(n1)?n1+l:n1;
       d=(z+big.x)-big.x;
       u.x=(u.x-d*w.x)-d*ww.x;
       l=(u.i[HIGH_HALF]&0x7ff00000)-nn;
      }
      r.i[HIGH_HALF]=m1;
      w.i[HIGH_HALF]=n;
      ww.i[HIGH_HALF]=n1;
      z=u.x*r.x;
      d=(z+big.x)-big.x;
      u.x=(u.x-d*w.x)-d*ww.x;
      if (ABS(u.x)<0.5*t.x) return (u.x!=0)?u.x:((x>0)?ZERO.x:nZERO.x);
      else
        if (ABS(u.x)>0.5*t.x) return (d>z)?u.x+t.x:u.x-t.x;
        else
        {z=u.x/t.x; d=(z+big.x)-big.x; return ((u.x-d*w.x)-d*ww.x);}
    }

  }   /*   (kx<0x7fe00000&&ky<0x7ff00000&&ky>=0x03500000)     */
  else {
    if (kx<0x7fe00000&&ky<0x7ff00000&&(ky>0||t.i[LOW_HALF]!=0)) {
      y=ABS(y)*t128.x;
      z=__ieee754_remainder(x,y)*t128.x;
      z=__ieee754_remainder(z,y)*tm128.x;
      return z;
    }
  else {
    if ((kx&0x7ff00000)==0x7fe00000&&ky<0x7ff00000&&(ky>0||t.i[LOW_HALF]!=0)) {
      y=ABS(y);
      z=2.0*__ieee754_remainder(0.5*x,y);
      d = ABS(z);
      if (d <= ABS(d-y)) return z;
      else return (z>0)?z-y:z+y;
    }
    else { /* if x is too big */
      if (kx == 0x7ff00000 && u.i[LOW_HALF] == 0 && y == 1.0)
       return x / x;
      if (kx>=0x7ff00000||(ky==0&&t.i[LOW_HALF]==0)||ky>0x7ff00000||
         (ky==0x7ff00000&&t.i[LOW_HALF]!=0))
       return (u.i[HIGH_HALF]&0x80000000)?nNAN.x:NAN.x;
      else return x;
    }
   }
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

float __ieee754_remainderf ( float  ,
float   
)

Definition at line 33 of file e_remainderf.c.

{
       int32_t hx,hp;
       u_int32_t sx;
       float p_half;

       GET_FLOAT_WORD(hx,x);
       GET_FLOAT_WORD(hp,p);
       sx = hx&0x80000000;
       hp &= 0x7fffffff;
       hx &= 0x7fffffff;

    /* purge off exception values */
       if(hp==0) return (x*p)/(x*p);             /* p = 0 */
       if((hx>=0x7f800000)||                     /* x not finite */
         ((hp>0x7f800000)))               /* p is NaN */
           return (x*p)/(x*p);


       if (hp<=0x7effffff) x = __ieee754_fmodf(x,p+p);  /* now x < 2p */
       if ((hx-hp)==0) return zero*x;
       x  = fabsf(x);
       p  = fabsf(p);
       if (hp<0x01000000) {
           if(x+x>p) {
              x-=p;
              if(x+x>=p) x -= p;
           }
       } else {
           p_half = (float)0.5*p;
           if(x>p_half) {
              x-=p;
              if(x>=p_half) x -= p;
           }
       }
       GET_FLOAT_WORD(hx,x);
       SET_FLOAT_WORD(x,hx^sx);
       return x;
}

Here is the call graph for this function:

Here is the caller graph for this function:

long double __ieee754_remainderl ( long  double,
long  double 
)

Definition at line 37 of file e_remainderl.c.

{
       int64_t hx,hp;
       u_int64_t sx,lx,lp;
       long double p_half;

       GET_LDOUBLE_WORDS64(hx,lx,x);
       GET_LDOUBLE_WORDS64(hp,lp,p);
       sx = hx&0x8000000000000000ULL;
       hp &= 0x7fffffffffffffffLL;
       hx &= 0x7fffffffffffffffLL;

    /* purge off exception values */
       if((hp|lp)==0) return (x*p)/(x*p);        /* p = 0 */
       if((hx>=0x7fff000000000000LL)||                  /* x not finite */
         ((hp>=0x7fff000000000000LL)&&                  /* p is NaN */
         (((hp-0x7fff000000000000LL)|lp)!=0)))
           return (x*p)/(x*p);


       if (hp<=0x7ffdffffffffffffLL) x = __ieee754_fmodl(x,p+p);      /* now x < 2p */
       if (((hx-hp)|(lx-lp))==0) return zero*x;
       x  = fabsl(x);
       p  = fabsl(p);
       if (hp<0x0002000000000000LL) {
           if(x+x>p) {
              x-=p;
              if(x+x>=p) x -= p;
           }
       } else {
           p_half = 0.5L*p;
           if(x>p_half) {
              x-=p;
              if(x>=p_half) x -= p;
           }
       }
       GET_LDOUBLE_MSW64(hx,x);
       SET_LDOUBLE_MSW64(x,hx^sx);
       return x;
}

Here is the call graph for this function:

Here is the caller graph for this function:

double __ieee754_scalb ( double  ,
double   
)

Definition at line 38 of file e_scalb.c.

{
#ifdef _SCALB_INT
       return __scalbn(x,fn);
#else
       if (__isnan(x)||__isnan(fn)) return x*fn;
       if (!__finite(fn)) {
           if(fn>0.0) return x*fn;
           else if (x == 0)
             return x;
           else if (!__finite (x))
             {
# ifdef FE_INVALID
              feraiseexcept (FE_INVALID);
# endif
              return __nan ("");
             }
           else       return x/(-fn);
       }
       if (__rint(fn)!=fn)
         {
# ifdef FE_INVALID
           feraiseexcept (FE_INVALID);
# endif
           return __nan ("");
         }
       if ( fn > 65000.0) return __scalbn(x, 65000);
       if (-fn > 65000.0) return __scalbn(x,-65000);
       return __scalbn(x,(int)fn);
#endif
}

Here is the call graph for this function:

Here is the caller graph for this function:

float __ieee754_scalbf ( float  ,
float   
)

Definition at line 35 of file e_scalbf.c.

{
#ifdef _SCALB_INT
       return __scalbnf(x,fn);
#else
       if (__isnanf(x)||__isnanf(fn)) return x*fn;
       if (!__finitef(fn)) {
           if(fn>(float)0.0) return x*fn;
           else if (x == 0)
             return x;
           else if (!__finitef (x))
             {
# ifdef FE_INVALID
              feraiseexcept (FE_INVALID);
# endif
              return __nanf ("");
             }
           else       return x/(-fn);
       }
       if (__rintf(fn)!=fn)
         {
# ifdef FE_INVALID
           feraiseexcept (FE_INVALID);
# endif
           return __nanf ("");
         }
       if ( fn > (float)65000.0) return __scalbnf(x, 65000);
       if (-fn > (float)65000.0) return __scalbnf(x,-65000);
       return __scalbnf(x,(int)fn);
#endif
}

Here is the call graph for this function:

Here is the caller graph for this function:

long double __ieee754_scalbl ( long  double,
long  double 
)

Definition at line 42 of file e_scalbl.c.

{
#ifdef _SCALB_INT
       return __scalbnl(x,fn);
#else
       if (__isnanl(x)||__isnanl(fn)) return x*fn;
       if (!__finitel(fn)) {
           if(fn>0.0) return x*fn;
           else if (x == 0)
             return x;
           else if (!__finitel (x))
             {
# ifdef FE_INVALID
              feraiseexcept (FE_INVALID);
# endif
              return __nanl ("");
             }
           else       return x/(-fn);
       }
       if (__rintl(fn)!=fn)
         {
# ifdef FE_INVALID
           feraiseexcept (FE_INVALID);
# endif
           return __nanl ("");
         }
       if ( fn > 65000.0) return __scalbnl(x, 65000);
       if (-fn > 65000.0) return __scalbnl(x,-65000);
       return __scalbnl(x,(int)fn);
#endif
}

Here is the call graph for this function:

Here is the caller graph for this function:

double __ieee754_sinh ( double  )

Definition at line 47 of file e_sinh.c.

{      
       double t,w,h;
       int32_t ix,jx;
       u_int32_t lx;

    /* High word of |x|. */
       GET_HIGH_WORD(jx,x);
       ix = jx&0x7fffffff;

    /* x is INF or NaN */
       if(ix>=0x7ff00000) return x+x;     

       h = 0.5;
       if (jx<0) h = -h;
    /* |x| in [0,22], return sign(x)*0.5*(E+E/(E+1))) */
       if (ix < 0x40360000) {             /* |x|<22 */
           if (ix<0x3e300000)             /* |x|<2**-28 */
              if(shuge+x>one) return x;/* sinh(tiny) = tiny with inexact */
           t = __expm1(fabs(x));
           if(ix<0x3ff00000) return h*(2.0*t-t*t/(t+one));
           return h*(t+t/(t+one));
       }

    /* |x| in [22, log(maxdouble)] return 0.5*exp(|x|) */
       if (ix < 0x40862e42)  return h*__ieee754_exp(fabs(x));

    /* |x| in [log(maxdouble), overflowthresold] */
       GET_LOW_WORD(lx,x);
       if (ix<0x408633ce || ((ix==0x408633ce)&&(lx<=(u_int32_t)0x8fb9f87d))) {
           w = __ieee754_exp(0.5*fabs(x));
           t = h*w;
           return t*w;
       }

    /* |x| > overflowthresold, sinh(x) overflow */
       return x*shuge;
}

Here is the call graph for this function:

Here is the caller graph for this function:

float __ieee754_sinhf ( float  )

Definition at line 32 of file e_sinhf.c.

{      
       float t,w,h;
       int32_t ix,jx;

       GET_FLOAT_WORD(jx,x);
       ix = jx&0x7fffffff;

    /* x is INF or NaN */
       if(ix>=0x7f800000) return x+x;     

       h = 0.5;
       if (jx<0) h = -h;
    /* |x| in [0,22], return sign(x)*0.5*(E+E/(E+1))) */
       if (ix < 0x41b00000) {             /* |x|<22 */
           if (ix<0x31800000)             /* |x|<2**-28 */
              if(shuge+x>one) return x;/* sinh(tiny) = tiny with inexact */
           t = __expm1f(fabsf(x));
           if(ix<0x3f800000) return h*((float)2.0*t-t*t/(t+one));
           return h*(t+t/(t+one));
       }

    /* |x| in [22, log(maxdouble)] return 0.5*exp(|x|) */
       if (ix < 0x42b17180)  return h*__ieee754_expf(fabsf(x));

    /* |x| in [log(maxdouble), overflowthresold] */
       if (ix<=0x42b2d4fc) {
           w = __ieee754_expf((float)0.5*fabsf(x));
           t = h*w;
           return t*w;
       }

    /* |x| > overflowthresold, sinh(x) overflow */
       return x*shuge;
}

Here is the call graph for this function:

Here is the caller graph for this function:

long double __ieee754_sinhl ( long  double)

Definition at line 6 of file e_sinhl.c.

{
  fputs ("__ieee754_sinhl not implemented\n", stderr);
  __set_errno (ENOSYS);
  return 0.0;
}

Here is the call graph for this function:

Here is the caller graph for this function:

double __ieee754_sqrt ( double  )

Definition at line 20 of file e_sqrt.c.

{
  asm ("fdsqrt.d %1,%0" : "=f" (x) : "fm" (x));
  return x;
}
float __ieee754_sqrtf ( float  )

Definition at line 20 of file e_sqrtf.c.

{
  double result;
  asm ("fssqrt.s %1,%0" : "=f" (result) : "dm" (x));
  return result;
}

Here is the call graph for this function:

Here is the caller graph for this function:

long double __ieee754_sqrtl ( long  double)

Definition at line 6 of file e_sqrtl.c.

{
  fputs ("__ieee754_sqrtl not implemented\n", stderr);
  __set_errno (ENOSYS);
  return 0.0;
}

Here is the call graph for this function:

Here is the caller graph for this function:

double __ieee754_y0 ( double  )

Definition at line 179 of file e_j0.c.

{
       double z, s,c,ss,cc,u,v,z2,z4,z6,u1,u2,u3,v1,v2;
       int32_t hx,ix,lx;

       EXTRACT_WORDS(hx,lx,x);
        ix = 0x7fffffff&hx;
    /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0, y0(0) is -inf.  */
       if(ix>=0x7ff00000) return  one/(x+x*x);
        if((ix|lx)==0) return -HUGE_VAL+x; /* -inf and overflow exception.  */
        if(hx<0) return zero/(zero*x);
        if(ix >= 0x40000000) {  /* |x| >= 2.0 */
        /* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0))
         * where x0 = x-pi/4
         *      Better formula:
         *              cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
         *                      =  1/sqrt(2) * (sin(x) + cos(x))
         *              sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
         *                      =  1/sqrt(2) * (sin(x) - cos(x))
         * To avoid cancellation, use
         *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
         * to compute the worse one.
         */
              __sincos (x, &s, &c);
                ss = s-c;
                cc = s+c;
       /*
        * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
        * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
        */
                if(ix<0x7fe00000) {  /* make sure x+x not overflow */
                    z = -__cos(x+x);
                    if ((s*c)<zero) cc = z/ss;
                    else            ss = z/cc;
                }
                if(ix>0x48000000) z = (invsqrtpi*ss)/__ieee754_sqrt(x);
                else {
                    u = pzero(x); v = qzero(x);
                    z = invsqrtpi*(u*ss+v*cc)/__ieee754_sqrt(x);
                }
                return z;
       }
       if(ix<=0x3e400000) { /* x < 2**-27 */
           return(U[0] + tpi*__ieee754_log(x));
       }
       z = x*x;
#ifdef DO_NOT_USE_THIS
       u = u00+z*(u01+z*(u02+z*(u03+z*(u04+z*(u05+z*u06)))));
       v = one+z*(v01+z*(v02+z*(v03+z*v04)));
#else
       u1 = U[0]+z*U[1]; z2=z*z;
       u2 = U[2]+z*U[3]; z4=z2*z2;
       u3 = U[4]+z*U[5]; z6=z4*z2;
       u = u1 + z2*u2 + z4*u3 + z6*U[6];
       v1 = one+z*V[0];
       v2 = V[1]+z*V[2];
       v = v1 + z2*v2 + z4*V[3];
#endif
       return(u/v + tpi*(__ieee754_j0(x)*__ieee754_log(x)));
}

Here is the call graph for this function:

Here is the caller graph for this function:

float __ieee754_y0f ( float  )

Definition at line 125 of file e_j0f.c.

{
       float z, s,c,ss,cc,u,v;
       int32_t hx,ix;

       GET_FLOAT_WORD(hx,x);
        ix = 0x7fffffff&hx;
    /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0, y0(0) is -inf.  */
       if(ix>=0x7f800000) return  one/(x+x*x);
        if(ix==0) return -HUGE_VALF+x;  /* -inf and overflow exception.  */
        if(hx<0) return zero/(zero*x);
        if(ix >= 0x40000000) {  /* |x| >= 2.0 */
        /* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0))
         * where x0 = x-pi/4
         *      Better formula:
         *              cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
         *                      =  1/sqrt(2) * (sin(x) + cos(x))
         *              sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
         *                      =  1/sqrt(2) * (sin(x) - cos(x))
         * To avoid cancellation, use
         *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
         * to compute the worse one.
         */
              __sincosf (x, &s, &c);
                ss = s-c;
                cc = s+c;
       /*
        * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
        * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
        */
                if(ix<0x7f000000) {  /* make sure x+x not overflow */
                    z = -__cosf(x+x);
                    if ((s*c)<zero) cc = z/ss;
                    else            ss = z/cc;
                }
                if(ix>0x48000000) z = (invsqrtpi*ss)/__ieee754_sqrtf(x);
                else {
                    u = pzerof(x); v = qzerof(x);
                    z = invsqrtpi*(u*ss+v*cc)/__ieee754_sqrtf(x);
                }
                return z;
       }
       if(ix<=0x32000000) { /* x < 2**-27 */
           return(u00 + tpi*__ieee754_logf(x));
       }
       z = x*x;
       u = u00+z*(u01+z*(u02+z*(u03+z*(u04+z*(u05+z*u06)))));
       v = one+z*(v01+z*(v02+z*(v03+z*v04)));
       return(u/v + tpi*(__ieee754_j0f(x)*__ieee754_logf(x)));
}

Here is the call graph for this function:

Here is the caller graph for this function:

long double __ieee754_y0l ( long  double)

Definition at line 815 of file e_j0l.c.

{
  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 = neval (z, Y0_2N, NY0_2N) / deval (z, Y0_2D, NY0_2D);
      p = TWOOPI * __ieee754_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.  */
  __sincosl (x, &s, &c);
  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) / __ieee754_sqrtl (x);
  return z;
}

Here is the call graph for this function:

Here is the caller graph for this function:

double __ieee754_y1 ( double  )

Definition at line 182 of file e_j1.c.

{
       double z, s,c,ss,cc,u,v,u1,u2,v1,v2,v3,z2,z4;
       int32_t hx,ix,lx;

       EXTRACT_WORDS(hx,lx,x);
        ix = 0x7fffffff&hx;
    /* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */
       if(ix>=0x7ff00000) return  one/(x+x*x);
        if((ix|lx)==0) return -HUGE_VAL+x; /* -inf and overflow exception.  */;
        if(hx<0) return zero/(zero*x);
        if(ix >= 0x40000000) {  /* |x| >= 2.0 */
              __sincos (x, &s, &c);
                ss = -s-c;
                cc = s-c;
                if(ix<0x7fe00000) {  /* make sure x+x not overflow */
                    z = __cos(x+x);
                    if ((s*c)>zero) cc = z/ss;
                    else            ss = z/cc;
                }
        /* y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x0)+q1(x)*cos(x0))
         * where x0 = x-3pi/4
         *      Better formula:
         *              cos(x0) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
         *                      =  1/sqrt(2) * (sin(x) - cos(x))
         *              sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
         *                      = -1/sqrt(2) * (cos(x) + sin(x))
         * To avoid cancellation, use
         *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
         * to compute the worse one.
         */
                if(ix>0x48000000) z = (invsqrtpi*ss)/__ieee754_sqrt(x);
                else {
                    u = pone(x); v = qone(x);
                    z = invsqrtpi*(u*ss+v*cc)/__ieee754_sqrt(x);
                }
                return z;
        }
        if(ix<=0x3c900000) {    /* x < 2**-54 */
            return(-tpi/x);
        }
        z = x*x;
#ifdef DO_NOT_USE_THIS
        u = U0[0]+z*(U0[1]+z*(U0[2]+z*(U0[3]+z*U0[4])));
        v = one+z*(V0[0]+z*(V0[1]+z*(V0[2]+z*(V0[3]+z*V0[4]))));
#else
       u1 = U0[0]+z*U0[1];z2=z*z;
       u2 = U0[2]+z*U0[3];z4=z2*z2;
       u  = u1 + z2*u2 + z4*U0[4];
       v1 = one+z*V0[0];
       v2 = V0[1]+z*V0[2];
       v3 = V0[3]+z*V0[4];
       v = v1 + z2*v2 + z4*v3;
#endif
        return(x*(u/v) + tpi*(__ieee754_j1(x)*__ieee754_log(x)-one/x));
}

Here is the call graph for this function:

Here is the caller graph for this function:

float __ieee754_y1f ( float  )

Definition at line 126 of file e_j1f.c.

{
       float z, s,c,ss,cc,u,v;
       int32_t hx,ix;

       GET_FLOAT_WORD(hx,x);
        ix = 0x7fffffff&hx;
    /* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */
       if(ix>=0x7f800000) return  one/(x+x*x);
        if(ix==0) return -HUGE_VALF+x;  /* -inf and overflow exception.  */
        if(hx<0) return zero/(zero*x);
        if(ix >= 0x40000000) {  /* |x| >= 2.0 */
              __sincosf (x, &s, &c);
                ss = -s-c;
                cc = s-c;
                if(ix<0x7f000000) {  /* make sure x+x not overflow */
                    z = __cosf(x+x);
                    if ((s*c)>zero) cc