s_sin.cpp

这是整套横扫千军3D版游戏的源码

  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;
  struct {
inline Double& d() {return DOUBLE_FROM_INT_PTR(&i[0]);}
inline Double& x() {return DOUBLE_FROM_INT_PTR(&i[0]);}
inline Double& d(int idx) {return DOUBLE_FROM_INT_PTR(&i[idx*(sizeof(double)/sizeof(i))]);}
inline Double& x(int idx) {return DOUBLE_FROM_INT_PTR(&i[idx*(sizeof(double)/sizeof(i))]);}
inline const Double& d() const {return CONST_DOUBLE_FROM_INT_PTR(&i[0]);}
inline const Double& x() const {return CONST_DOUBLE_FROM_INT_PTR(&i[0]);}
inline const Double& d(int idx) const {return CONST_DOUBLE_FROM_INT_PTR(&i[idx*(sizeof(double)/sizeof(i))]);}
inline const Double& x(int idx) const {return CONST_DOUBLE_FROM_INT_PTR(&i[idx*(sizeof(double)/sizeof(i))]);}
int4 i[2];} 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;
  struct {
inline Double& d() {return DOUBLE_FROM_INT_PTR(&i[0]);}
inline Double& x() {return DOUBLE_FROM_INT_PTR(&i[0]);}
inline Double& d(int idx) {return DOUBLE_FROM_INT_PTR(&i[idx*(sizeof(double)/sizeof(i))]);}
inline Double& x(int idx) {return DOUBLE_FROM_INT_PTR(&i[idx*(sizeof(double)/sizeof(i))]);}
inline const Double& d() const {return CONST_DOUBLE_FROM_INT_PTR(&i[0]);}
inline const Double& x() const {return CONST_DOUBLE_FROM_INT_PTR(&i[0]);}
inline const Double& d(int idx) const {return CONST_DOUBLE_FROM_INT_PTR(&i[idx*(sizeof(double)/sizeof(i))]);}
inline const Double& x(int idx) const {return CONST_DOUBLE_FROM_INT_PTR(&i[idx*(sizeof(double)/sizeof(i))]);}
int4 i[2];} 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
}