s_sin.c

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

    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;

}

/************************************************************************/
/*  Routine compute sin(x) for  2^-26 < |x|< 0.25 by  Taylor with more   */
/* precision  and if still doesn't accurate enough by mpsin   or dubsin */
/************************************************************************/

static double slow(double x) {
static const double th2_36 = 206158430208.0;   /*    1.5*2**37   */
 double y,x1,x2,xx,r,t,res,cor,w[2];
 x1=(x+th2_36)-th2_36;
 y = aa.x*x1*x1*x1;
 r=x+y;
 x2=x-x1;
 xx=x*x;
 t = (((((s5.x*xx + s4.x)*xx + s3.x)*xx + s2.x)*xx + bb.x)*xx + 3.0*aa.x*x1*x2)*x +aa.x*x2*x2*x2;
 t=((x-r)+y)+t;
 res=r+t;
 cor = (r-res)+t;
 if (res == res + 1.0007*cor) return res;
 else {
   __dubsin(ABS(x),0,w);
   if (w[0] == w[0]+1.000000001*w[1]) return (x>0)?w[0]:-w[0];
   else return (x>0)?__mpsin(x,0):-__mpsin(-x,0);
 }
}
/*******************************************************************************/
/* Routine compute sin(x) for   0.25<|x|< 0.855469 by  sincos.tbl   and Taylor */
/* and if result still doesn't accurate enough by mpsin   or dubsin            */
/*******************************************************************************/

static double slow1(double x) {
  mynumber u;
  double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,c1,c2,xx,cor,res;
  static const double t22 = 6291456.0;
  int4 k;
  y=ABS(x);
  u.x=big.x+y;
  y=y-(u.x-big.x);
  xx=y*y;
  s = y*xx*(sn3 +xx*sn5);
  c = xx*(cs2 +xx*(cs4 + xx*cs6));
  k=u.i[LOW_HALF]<<2;
  sn=sincos.x[k];          /* Data          */
  ssn=sincos.x[k+1];       /*  from         */
  cs=sincos.x[k+2];        /*   tables      */
  ccs=sincos.x[k+3];       /*    sincos.tbl */
  y1 = (y+t22)-t22;
  y2 = y - y1;
  c1 = (cs+t22)-t22;
  c2=(cs-c1)+ccs;
  cor=(ssn+s*ccs+cs*s+c2*y+c1*y2)-sn*c;
  y=sn+c1*y1;
  cor = cor+((sn-y)+c1*y1);
  res=y+cor;
  cor=(y-res)+cor;
  if (res == res+1.0005*cor) return (x>0)?res:-res;
  else {
    __dubsin(ABS(x),0,w);
    if (w[0] == w[0]+1.000000005*w[1]) return (x>0)?w[0]:-w[0];
    else return (x>0)?__mpsin(x,0):-__mpsin(-x,0);
  }
}
/**************************************************************************/
/*  Routine compute sin(x) for   0.855469  <|x|<2.426265  by  sincos.tbl  */
/* and if result still doesn't accurate enough by mpsin   or dubsin       */
/**************************************************************************/
static double slow2(double x) {
  mynumber u;
  double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,e1,e2,xx,cor,res,del;
  static const double t22 = 6291456.0;
  int4 k;
  y=ABS(x);
  y = hp0.x-y;
  if (y>=0) {
    u.x = big.x+y;
    y = y-(u.x-big.x);
    del = hp1.x;
  }
  else {
    u.x = big.x-y;
    y = -(y+(u.x-big.x));
    del = -hp1.x;
  }
  xx=y*y;
  s = y*xx*(sn3 +xx*sn5);
  c = y*del+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];
  y1 = (y+t22)-t22;
  y2 = (y - y1)+del;
  e1 = (sn+t22)-t22;
  e2=(sn-e1)+ssn;
  cor=(ccs-cs*c-e1*y2-e2*y)-sn*s;
  y=cs-e1*y1;
  cor = cor+((cs-y)-e1*y1);
  res=y+cor;
  cor=(y-res)+cor;
  if (res == res+1.0005*cor) return (x>0)?res:-res;
  else {
    y=ABS(x)-hp0.x;
    y1=y-hp1.x;
    y2=(y-y1)-hp1.x;
    __docos(y1,y2,w);
    if (w[0] == w[0]+1.000000005*w[1]) return (x>0)?w[0]:-w[0];
    else return (x>0)?__mpsin(x,0):-__mpsin(-x,0);
  }
}
/***************************************************************************/
/*  Routine compute sin(x+dx) (Double-Length number) where x is small enough*/
/* to use Taylor series around zero and   (x+dx)                            */
/* in first or third quarter of unit circle.Routine receive also            */
/* (right argument) the  original   value of x for computing error of      */
/* result.And if result not accurate enough routine calls mpsin1 or dubsin */
/***************************************************************************/

static double sloww(double x,double dx, double orig) {
  static const double th2_36 = 206158430208.0;   /*    1.5*2**37   */
  double y,x1,x2,xx,r,t,res,cor,w[2],a,da,xn;
  union {int4 i[2]; double x;} v;
  int4 n;
  x1=(x+th2_36)-th2_36;
  y = aa.x*x1*x1*x1;
  r=x+y;
  x2=(x-x1)+dx;
  xx=x*x;
  t = (((((s5.x*xx + s4.x)*xx + s3.x)*xx + s2.x)*xx + bb.x)*xx + 3.0*aa.x*x1*x2)*x +aa.x*x2*x2*x2+dx;
  t=((x-r)+y)+t;
  res=r+t;
  cor = (r-res)+t;
  cor = (cor>0)? 1.0005*cor+ABS(orig)*3.1e-30 : 1.0005*cor-ABS(orig)*3.1e-30;
  if (res == res + cor) return res;
  else {
    (x>0)? __dubsin(x,dx,w) : __dubsin(-x,-dx,w);
    cor = (w[1]>0)? 1.000000001*w[1] + ABS(orig)*1.1e-30 : 1.000000001*w[1] - ABS(orig)*1.1e-30;
    if (w[0] == w[0]+cor) return (x>0)?w[0]:-w[0];
    else {
      t = (orig*hpinv.x + toint.x);
      xn = t - toint.x;
      v.x = t;
      y = (orig - xn*mp1.x) - xn*mp2.x;
      n =v.i[LOW_HALF]&3;
      da = xn*pp3.x;
      t=y-da;
      da = (y-t)-da;
      y = xn*pp4.x;
      a = t - y;
      da = ((t-a)-y)+da;
      if (n&2) {a=-a; da=-da;}
      (a>0)? __dubsin(a,da,w) : __dubsin(-a,-da,w);
      cor = (w[1]>0)? 1.000000001*w[1] + ABS(orig)*1.1e-40 : 1.000000001*w[1] - ABS(orig)*1.1e-40;
      if (w[0] == w[0]+cor) return (a>0)?w[0]:-w[0];
      else return __mpsin1(orig);
    }
  }
}
/***************************************************************************/
/*  Routine compute sin(x+dx)   (Double-Length number) where x in first or  */
/*  third quarter of unit circle.Routine receive also (right argument) the  */
/*  original   value of x for computing error of result.And if result not  */
/* accurate enough routine calls  mpsin1   or dubsin                       */
/***************************************************************************/

static double sloww1(double x, double dx, double orig) {
  mynumber u;
  double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,c1,c2,xx,cor,res;
  static const double t22 = 6291456.0;
  int4 k;
  y=ABS(x);
  u.x=big.x+y;
  y=y-(u.x-big.x);
  dx=(x>0)?dx:-dx;
  xx=y*y;
  s = y*xx*(sn3 +xx*sn5);