s_sin.c

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

  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];
  y1 = (y+t22)-t22;
  y2 = (y - y1)+dx;
  c1 = (cs+t22)-t22;
  c2=(cs-c1)+ccs;
  cor=(ssn+s*ccs+cs*s+c2*y+c1*y2-sn*y*dx)-sn*c;
  y=sn+c1*y1;
  cor = cor+((sn-y)+c1*y1);
  res=y+cor;
  cor=(y-res)+cor;
  cor = (cor>0)? 1.0005*cor+3.1e-30*ABS(orig) : 1.0005*cor-3.1e-30*ABS(orig);
  if (res == res + cor) return (x>0)?res:-res;
  else {
    __dubsin(ABS(x),dx,w);
    cor = (w[1]>0)? 1.000000005*w[1]+1.1e-30*ABS(orig) : 1.000000005*w[1]-1.1e-30*ABS(orig);
    if (w[0] == w[0]+cor) return (x>0)?w[0]:-w[0];
  else  return __mpsin1(orig);
  }
}
/***************************************************************************/
/*  Routine compute sin(x+dx)   (Double-Length number) where x in second or */
/*  fourth quarter of unit circle.Routine receive also  the  original value */
/* and quarter(n= 1or 3)of x for computing error of result.And if result not*/
/* accurate enough routine calls  mpsin1   or dubsin                       */
/***************************************************************************/

static double sloww2(double x, double dx, double orig, int n) {
  mynumber u;
  double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,e1,e2,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);
  c = y*dx+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)+dx;
  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;
  cor = (cor>0)? 1.0005*cor+3.1e-30*ABS(orig) : 1.0005*cor-3.1e-30*ABS(orig);
  if (res == res + cor) return (n&2)?-res:res;
  else {
   __docos(ABS(x),dx,w);
   cor = (w[1]>0)? 1.000000005*w[1]+1.1e-30*ABS(orig) : 1.000000005*w[1]-1.1e-30*ABS(orig);
   if (w[0] == w[0]+cor) return (n&2)?-w[0]:w[0];
   else  return __mpsin1(orig);
  }
}
/***************************************************************************/
/*  Routine compute sin(x+dx) or cos(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 other routines    */
/***************************************************************************/

static double bsloww(double x,double dx, double orig,int n) {
  static const double th2_36 = 206158430208.0;   /*    1.5*2**37   */
  double y,x1,x2,xx,r,t,res,cor,w[2];
#if 0
  double a,da,xn;
  union {int4 i[2]; double x;} v;
#endif
  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+1.1e-24 : 1.0005*cor-1.1e-24;
  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] + 1.1e-24 : 1.000000001*w[1] - 1.1e-24;
    if (w[0] == w[0]+cor) return (x>0)?w[0]:-w[0];
    else return (n&1)?__mpcos1(orig):__mpsin1(orig);
  }
}

/***************************************************************************/
/*  Routine compute sin(x+dx)  or cos(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  other routines         */
/***************************************************************************/

static double bsloww1(double x, double dx, double orig,int n) {
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);
 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];
 y1 = (y+t22)-t22;
 y2 = (y - y1)+dx;
 c1 = (cs+t22)-t22;
 c2=(cs-c1)+ccs;
 cor=(ssn+s*ccs+cs*s+c2*y+c1*y2-sn*y*dx)-sn*c;
 y=sn+c1*y1;
 cor = cor+((sn-y)+c1*y1);
 res=y+cor;
 cor=(y-res)+cor;
 cor = (cor>0)? 1.0005*cor+1.1e-24 : 1.0005*cor-1.1e-24;
 if (res == res + cor) return (x>0)?res:-res;
 else {
   __dubsin(ABS(x),dx,w);
   cor = (w[1]>0)? 1.000000005*w[1]+1.1e-24: 1.000000005*w[1]-1.1e-24;
   if (w[0] == w[0]+cor) return (x>0)?w[0]:-w[0];
   else  return (n&1)?__mpcos1(orig):__mpsin1(orig);
 }
}

/***************************************************************************/
/*  Routine compute sin(x+dx)  or cos(x+dx) (Double-Length number) where x  */
/* in second or fourth quarter of unit circle.Routine receive also  the     */
/* original value and quarter(n= 1or 3)of x for computing error of result.  */
/* And if result not accurate enough routine calls  other routines          */
/***************************************************************************/

static double bsloww2(double x, double dx, double orig, int n) {
mynumber u;
 double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,e1,e2,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);
 c = y*dx+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)+dx;
 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;
 cor = (cor>0)? 1.0005*cor+1.1e-24 : 1.0005*cor-1.1e-24;
 if (res == res + cor) return (n&2)?-res:res;
 else {
   __docos(ABS(x),dx,w);
   cor = (w[1]>0)? 1.000000005*w[1]+1.1e-24 : 1.000000005*w[1]-1.1e-24;
   if (w[0] == w[0]+cor) return (n&2)?-w[0]:w[0];
   else  return (n&1)?__mpsin1(orig):__mpcos1(orig);
 }
}

/************************************************************************/
/*  Routine compute cos(x) for  2^-27 < |x|< 0.25 by  Taylor with more   */
/* precision  and if still doesn't accurate enough by mpcos   or docos  */
/************************************************************************/

static double cslow2(double x) {
  mynumber u;
  double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,e1,e2,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];
  ssn=sincos.x[k+1];
  cs=sincos.x[k+2];
  ccs=sincos.x[k+3];
  y1 = (y+t22)-t22;
  y2 = y - y1;
  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 res;
  else {
    y=ABS(x);
    __docos(y,0,w);
    if (w[0] == w[0]+1.000000005*w[1]) return w[0];
    else return __mpcos(x,0);
  }
}

/***************************************************************************/
/*  Routine compute cos(x+dx) (Double-Length number) where x is small enough*/
/* to use Taylor series around zero and   (x+dx) .Routine receive also      */
/* (right argument) the  original   value of x for computing error of      */
/* result.And if result not accurate enough routine calls other routines    */
/***************************************************************************/


static double csloww(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;
    /* Taylor series */
  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==1) {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 __mpcos1(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  other routines                            */
/***************************************************************************/

static double csloww1(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);
  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];
  y1 = (y+t22)-t22;
  y2 = (y - y1)+dx;
  c1 = (cs+t22)-t22;
  c2=(cs-c1)+ccs;
  cor=(ssn+s*ccs+cs*s+c2*y+c1*y2-sn*y*dx)-sn*c;
  y=sn+c1*y1;
  cor = cor+((sn-y)+c1*y1);
  res=y+cor;
  cor=(y-res)+cor;
  cor = (cor>0)? 1.0005*cor+3.1e-30*ABS(orig) : 1.0005*cor-3.1e-30*ABS(orig);
  if (res == res + cor) return (x>0)?res:-res;
  else {
    __dubsin(ABS(x),dx,w);
    cor = (w[1]>0)? 1.000000005*w[1]+1.1e-30*ABS(orig) : 1.000000005*w[1]-1.1e-30*ABS(orig);
    if (w[0] == w[0]+cor) return (x>0)?w[0]:-w[0];
    else  return __mpcos1(orig);
  }
}


/***************************************************************************/
/*  Routine compute sin(x+dx)   (Double-Length number) where x in second or */
/*  fourth quarter of unit circle.Routine receive also  the  original value */
/* and quarter(n= 1or 3)of x for computing error of result.And if result not*/
/* accurate enough routine calls  other routines                            */
/***************************************************************************/

static double csloww2(double x, double dx, double orig, int n) {
  mynumber u;
  double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,e1,e2,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);
  c = y*dx+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)+dx;
  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;
  cor = (cor>0)? 1.0005*cor+3.1e-30*ABS(orig) : 1.0005*cor-3.1e-30*ABS(orig);
  if (res == res + cor) return (n)?-res:res;
  else {
    __docos(ABS(x),dx,w);
    cor = (w[1]>0)? 1.000000005*w[1]+1.1e-30*ABS(orig) : 1.000000005*w[1]-1.1e-30*ABS(orig);
    if (w[0] == w[0]+cor) return (n)?-w[0]:w[0];
    else  return __mpcos1(orig);
  }
}

weak_alias (__cos, cos)
weak_alias (__sin, sin)

#ifdef NO_LONG_DOUBLE
strong_alias (__sin, __sinl)
weak_alias (__sin, sinl)
strong_alias (__cos, __cosl)
weak_alias (__cos, cosl)
#endif