math.cc

很好用的库

  if (fabs(x.h())<1.0e-11) { sinx=x; cosx=1.0-0.5*sqr(x); return; }
  int a,b; doubledouble sins, coss, k1, k3, t2, s, s2, sinb, cosb;
  k1=x/doubledouble::TwoPi(); k3=k1-rint(k1);
  t2=4*k3;
  a=int(rint(t2));
  b=int(rint((8*(t2-a))));
  s=doubledouble::Pi()*(k3+k3-(8*a+b)/16.0);
  s2=sqr(s);
  sins=s*(sinsTab[0]+(sinsTab[1]+(sinsTab[2]+(sinsTab[3]+(sinsTab[4]+
         (sinsTab[5]+sinsTab[6]*s2)*s2)*s2)*s2)*s2)*s2);
  coss=sqrt(1.0-sqr(sins)); // ok, sins is small
  sinb= (b>=0) ? tab[b] : -tab[-b]; cosb=tab[8-abs(b)];
  // sin(x)=
  // sin(s)(cos(1/2 a Pi) cos(1/16 b Pi) - sin(1/2 a Pi) sin(1/16 b Pi))
  // cos(s)(sin(1/2 a Pi) cos(1/16 b Pi) + cos(1/2 a Pi) sin(1/16 b Pi))
  // cos(x)=
  // cos(s)(cos(1/2 a Pi) cos(1/16 b Pi) - sin(1/2 a Pi) sin(1/16 b Pi))
  //-sin(s)(sin(1/2 a Pi) cos(1/16 b Pi) + cos(1/2 a Pi) sin(1/16 b Pi))
  if (0==a) {
    sinx= sins*cosb+coss*sinb;
    cosx= coss*cosb-sins*sinb;
  } else if (1==a) {
    sinx=-sins*sinb+coss*cosb;
    cosx=-coss*sinb-sins*cosb;
  } else if (-1==a) {
    sinx= sins*sinb-coss*cosb;
    cosx= coss*sinb+sins*cosb;
  } else  { // |a|=2
    sinx=-sins*cosb-coss*sinb;
    cosx=-coss*cosb+sins*sinb;
  }
  return;
}

// cos
inline doubledouble cos(const doubledouble& x) { return sin(doubledouble::Pion2()-x); }

// hyperbolic
inline doubledouble sinh(const doubledouble& x) { 
  if (fabs(x.h())<1.0e-5) { // avoid cancellation in e^x-e^(-x), use Taylor series...
    doubledouble q=sqr(x);  return x*(1+q/6*(1+q/20*(1+q/42))); }
  doubledouble t=exp(x); 
  return 0.5*(t-recip(t)); 
}
inline doubledouble cosh(const doubledouble& x) { doubledouble t=exp(x); return 0.5*(t+recip(t)); }
inline doubledouble tanh(const doubledouble& z) {
  doubledouble e;
  if (z.h()>0.0) { e=exp(-2.0*z); return (1.0-e)/(1.0+e); }
            else { e=exp( 2.0*z); return (e-1.0)/(1.0+e); }
}

inline doubledouble atan(const doubledouble& x) {
  double xh=x.h();
  if (0.0==xh) return doubledouble(0.0);
  // Asymptotic formula for large |x|...
  if (fabs(xh)>1.0e6) 
    { doubledouble r=recip(x), r2=sqr(r); return doubledouble::Pion2()-r+r2*r*(1.0-"0.6"*r2)/3; }
  doubledouble s,c,a=atan(xh); // a=double approx to result
  sincos(a,s,c);
  return a+c*(c*x-s);  // Newton step
}

inline doubledouble atan2(const doubledouble& qy, const doubledouble& qx) { // Based on GNU libc atan2.c
  static const double one=1.0, zero=0.0; double x, y;
  x=qx.h(); y=qy.h();
  double signx, signy;
  if (x!=x) return qx; // x=NaN
  if (y!=y) return qy;
  signy=copysign(one,y);
  signx=copysign(one,x);
  if (y==zero) return signx==one ? qy : Qcopysign(doubledouble::Pi(),signy);
  if (x==zero) return Qcopysign(doubledouble::Pion2(),signy);
  if (__isinf(x)) {
    if (__isinf(y)) return Qcopysign(signx==one ? doubledouble::Pion4() : 3.0*doubledouble::Pion4(),signy);
    else            return Qcopysign(signx==one ? doubledouble(0.0) : doubledouble::Pi(),signy);
  }
  if (__isinf(y)) return Qcopysign(doubledouble::Pion2(),signy);
  doubledouble aqy=fabs(qy);
  if (x<0.0) // X is negative.
    return Qcopysign(doubledouble::Pi()-atan(aqy/(-qx)),signy);
  return Qcopysign(atan(aqy/qx),signy);
}

// arcsin
inline doubledouble asin(const doubledouble& x) {
  if (fabs(x)>doubledouble(1.0)) 
    { cerr<<"\ndoubledouble |Argument|>1 in asin!\n"; DOMAIN_ERROR; }
  return atan2(x,sqrt(1.0-sqr(x)));
}

// KMB 96 Oct 28 version 0.3 97 Dec 22 OK now.
// erfc written 97 Dec 22, NOT CHECKED!!!!!!!
// Based on series and continued fraction for incomplete gamma function.

inline doubledouble erf(const doubledouble x) {
  if (0.0==x.h()) return 0.0;
  int i; doubledouble y,r;
  // const below is evalf(1/sqrt(Pi),40)
  static const doubledouble oneonrootpi="0.564189583547756286948079451560772585844";
  const double cut=1.5;  // switch to continued fraction here
  y=fabs(x);
  if (y.h()>26.0) { // erf is +or- 1 to 300 dp.
    r=1;
  } else if (y.h()<cut) { // use power series
    doubledouble ap=0.5,s,t,x2;
    s=t=2.0; x2=sqr(x);
    for (i=0; i<200; i++) {
      ap+=1;
      t*=x2/ap;
      s+=t;
      if (fabs(t.h())<1e-35*fabs(s.h())) {
        r=x*oneonrootpi*s/exp(x2);
	break;
      }
    }
    if (i>199) 
      { cerr<<"\ndoubledouble: no convergence in erf power series, x="<<x<<endl; exit(1); }
  } else { // |x|>=cut, use continued fraction, Lentz method
    double an,small=1e-300;
    doubledouble b,c,d,h,del,x2;
    x2=sqr(x);
    b=x2+0.5;
    c=1.0e300;
    d=recip(b);
    h=d;
    for (i=1; i<300; i++) {
      an=i*(0.5-i);
      b+=2.0;
      d=an*d+b;
      if (fabs(d.h())<small) d=small;
      c=b+an/c;
      if (fabs(c.h())<small) c=small;
      d=recip(d);
      del=d*c;
      h*=del;
      if (del.h()==1.0 && del.l()<1.0e-30) break;
      if (299==i)
        { cerr<<"\ndoubledouble: no convergence in erf continued fraction, x="<<x<<endl; exit(1); }
    }
    r=1.0-oneonrootpi*sqrt(x2)/exp(x2)*h;
  } 
  if (x.h()>0.0) return r; else return -r;
}

inline doubledouble erfc(const doubledouble x) {
  if (0.0==x.h()) return 1.0;
  if (x.h()<0.0) return 1.0-erf(x);
  int i; doubledouble y,r;
  // const below is evalf(1/sqrt(Pi),40)
  static const doubledouble oneonrootpi="0.564189583547756286948079451560772585844";
  const double cut=1.5;  // switch to continued fraction here
  y=fabs(x);
  if (y.h()<cut) { // use power series
    doubledouble ap=0.5,s,t,x2;
    s=t=2.0; x2=sqr(x);
    for (i=0; i<200; i++) {
      ap+=1;
      t*=x2/ap;
      s+=t;
      if (fabs(t.h())<1e-35*fabs(s.h())) {
        r=1.0-x*oneonrootpi*s/exp(x2);
	break;
      }
    }
    if (i>199) { cerr<<"\ndoubledouble: no convergence in erfc power series, x="<<x<<endl; exit(1); }
  } else { // |x|>=cut, use continued fraction
    double an,small=1e-300;
    doubledouble b,c,d,h,del,x2;
    x2=sqr(x);
    b=x2+0.5;
    c=1e300;
    h=d=recip(b);
    for (i=1; i<300; i++) {
      an=i*(0.5-i);
      b+=2.0;
      d=an*d+b;
      if (fabs(d.h())<small) d=small;
      c=b+an/c;
      if (fabs(c.h())<small) c=small;
      d=recip(d);
      del=d*c;
      h*=del;
      if (del.h()==1.0 && del.l()<1.0e-30) break;
      if (299==i) { cerr<<"\ndoubledouble: no convergence in erfc continued fraction, x="<<x<<endl; exit(1);}
    }
    r=oneonrootpi*sqrt(x2)/exp(x2)*h;
  } 
  return r;
}

inline doubledouble gamma(const doubledouble x) {
  const int n=43; // don't really need so many!
  static const doubledouble c[n]={ // Taylor coefficients for 1/gamma(1+x)-x...
    "+0.5772156649015328606065120900824024310421593359",
    "-0.6558780715202538810770195151453904812797663805",
    "-0.0420026350340952355290039348754298187113945004",
    "+0.1665386113822914895017007951021052357177815022",
    "-0.0421977345555443367482083012891873913016526841",
    "-0.0096219715278769735621149216723481989753629422",
    "+0.0072189432466630995423950103404465727099048009",
    "-0.0011651675918590651121139710840183886668093337",
    "-0.0002152416741149509728157299630536478064782419",
    "+0.0001280502823881161861531986263281643233948920",
    "-0.0000201348547807882386556893914210218183822948",
    "-0.0000012504934821426706573453594738330922423226",
    "+0.0000011330272319816958823741296203307449433240",
    "-0.0000002056338416977607103450154130020572836512",
    "+0.0000000061160951044814158178624986828553428672",
    "+0.0000000050020076444692229300556650480599913030",
    "-0.0000000011812745704870201445881265654365055777",
    "+1.0434267116911005104915403323122501914007098231E-10",
    "+7.7822634399050712540499373113607772260680861813E-12",
    "-3.6968056186422057081878158780857662365709634513E-12",
    "+5.1003702874544759790154813228632318027268860697E-13",
    "-2.0583260535665067832224295448552374197460910808E-14",
    "-5.3481225394230179823700173187279399489897154781E-15",
    "+1.2267786282382607901588938466224224281654557504E-15",
    "-1.1812593016974587695137645868422978312115572918E-16",
    "+1.1866922547516003325797772429286740710884940796E-18",
    "+1.4123806553180317815558039475667090370863507503E-18",
    "-2.2987456844353702065924785806336992602845059314E-19",
    "+1.7144063219273374333839633702672570668126560625E-20",
    "+1.3373517304936931148647813951222680228750594717E-22",
    "-2.0542335517666727893250253513557337966820379352E-22",
    "+2.7360300486079998448315099043309820148653116958E-23",
    "-1.7323564459105166390574284515647797990697491087E-24",
    "-2.3606190244992872873434507354275310079264135521E-26",
    "+1.8649829417172944307184131618786668989458684290E-26",
    "+2.2180956242071972043997169136268603797317795006E-27",
    "+1.2977819749479936688244144863305941656194998646E-28",
    "+1.1806974749665284062227454155099715185596846378E-30",
    "-1.1245843492770880902936546742614395121194117955E-30",
    "+1.2770851751408662039902066777511246477487720656E-31",
    "-7.3914511696151408234612893301085528237105689924E-33",
    "+1.1347502575542157609541652594693063930086121959E-35",
    "+4.6391346410587220299448049079522284630579686797E-35"
  };
  doubledouble ss=x,f=1.0,sum=0.0,one=1.0;
  while (ss>one) { ss-=1; f*=ss; }
  while (ss<one) { f/=ss; ss+=1; }
  if (ss==one) return f;
  ss-=1.0;
  for (int i=n-1; i>=0; i--) sum=c[i]+ss*sum;
  return f/(ss*sum+1);
}

// end of math.cc
}