math.cc

很好用的库

// math.cc
// KMB 98 Feb 03
/*
Copyright (C) 1997 Keith Martin Briggs

Except where otherwise indicated,
this program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 2
of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
*/

//#include "config.h"
//#include "doubledouble.h"
#include <math.h>
#if (HAVE_COPYSIGN == 0)
extern "C" {
  double copysign(double x, double y);
}
#endif
#include <stdlib.h>
//JGS AIX doesn't have this
//   #include <nan.h>  // defines NAN, at least in gcc
#include <assert.h>

namespace std {

inline doubledouble exp(const doubledouble& x) { 
/* Uses some ideas by Alan Miller
   Method:
    x    x.log2(e)    nint[x.log2(e)] + frac[x.log2(e)]
   e  = 2          = 2
 
                      iy    fy
                   = 2   . 2
   Then
    fy    y.loge(2)
   2   = e
 
   Now y.loge(2) will be less than 0.3466 in absolute value.
   This is halved and a Pade' approximant is used to approximate e^x over
   the region (-0.1733, +0.1733).   This approximation is then squared.
*/
  if (x.h()<-744.4400719213812) return 0.0; // exp(x)<1e-300
  int iy;
  doubledouble y,temp,ysq,sum1,sum2;
  y=x/doubledouble::Log2();
  temp=iy=rint(y);
  y=(y-temp)*doubledouble::Log2();
  y=ldexp(y,-1);
  ysq=sqr(y);
  sum1=y*((((ysq+3960.)*ysq+2162160.)*ysq+302702400.)*ysq+8821612800.);
  sum2=(((90.*ysq+110880.)*ysq+30270240.)*ysq+2075673600.)*ysq+17643225600.;
/*
                      sum2 + sum1         2.sum1
  Now approximation = ----------- = 1 + ----------- = 1 + 2.temp
                      sum2 - sum1       sum2 - sum1
 
  Then (1 + 2.temp)^2 = 4.temp.(1 + temp) + 1
*/
  temp=sum1/(sum2-sum1);
  y=temp*(temp+1);
  y=ldexp(y,2);
  return ldexp(y+1,iy);
}

// See Higham "Accuracy and Stability of Numerical Algorithms", p 511
inline doubledouble hypot(const doubledouble a, const doubledouble b) { 
  doubledouble p,q,r,s,aa,ab,four=4.0;
  aa=fabs(a); ab=fabs(b);
  if (aa>ab) { p=aa; q=ab; } else { q=aa; p=ab; } // now p>=q
  if (0.0==p.h()) return q;
  while (1) {
     r=sqr(q/p);
     if (four==(aa=r+four)) return p;
     s=r/aa; p+=2*s*p; q*=s;
  } // Never get here
}

// square root
inline doubledouble sqrt(const doubledouble& y) {
  x86_FIX
  double c,p,q,hx,tx,u,uu,cc,hi=y.h();
  if (hi<0.0) {
    cerr << "\ndoubledouble: attempt to take sqrt of " << hi << endl;
    DOMAIN_ERROR;
  }
  if (0.0==hi) return y;
  c=sqrt(hi);
  p=doubledouble::Split()*c; hx=c-p; hx+=p; tx=c-hx;
  p=hx*hx;
  q=2.0*hx*tx;
  u=p+q;
  uu=(p-u)+q+tx*tx;
  cc=(((y.h()-u)-uu)+y.l())/(c+c);
  u=c+cc;
  doubledouble r=doubledouble(u,cc+(c-u));
  END_x86_FIX
  return r;
}

// Natural logarithm
inline doubledouble log(const doubledouble& x) { // Newton method
  if (x.h()<0.0) {
    cerr<<"\ndoubledouble: attempt to take log of negative argument!\n";
    DOMAIN_ERROR;
  }
  if (x.h()==0.0) {
    cerr<<"\ndoubledouble: attempt to take log(0)!\n";
    DOMAIN_ERROR;
  }
  doubledouble s=log(x.h()), e=exp(s); // s=double approximation to result
  return s+(x-e)/e;  // Newton correction, good enough
  //doubledouble d=(x-e)/e; return s+d-0.5*sqr(d);  // or 2nd order correction
}

// logarithm base 10
inline doubledouble log10(const doubledouble& t) { 
  static const doubledouble one_on_log10="0.4342944819032518276511289189166050822944";
  return one_on_log10*log(t);
}

// doubledouble^doubledouble
inline doubledouble pow(const doubledouble& a, const doubledouble& b) {
  return exp(b*log(a));
}

// doubledouble^int
inline doubledouble powint(const doubledouble& u, const int c) {
  if (c<0) return recip(powint(u,-c));
  switch (c) {
    case 0: return u.h()==0.0?doubledouble(pow(0.0,0.0)):doubledouble(1.0); // let math.h handle 0^0.
    case 1: return u;  
    case 2: return sqr(u);
    case 3: return sqr(u)*u;
    default: { // binary method
      int n=c, m; doubledouble y=1.0, z=u; if (n<0) n=-n;
      while (1) {
        m=n; n/=2;
        if (n+n!=m) { // if m odd
           y*=z; if (0==n) return y;
        }
        z=sqr(z); 
      }
    }
  }
}

// Like Miller's modr
// a=n*b+rem, |rem|<=b/2, exact result.
inline doubledouble modr(const doubledouble a, const doubledouble b, int& n, doubledouble& rem) {
  if (0.0==b.h()) {
    cerr<<"\ndoubledouble: divisor is zero in modr!\n";
    DOMAIN_ERROR;
  }
  doubledouble temp;
  temp=a/b;
  n=int(rint(temp.h()));
  temp=n*doubledouble(b.h());
  rem=doubledouble(a.h());
  temp=rem-temp;
  rem=doubledouble(a.l());
  temp=rem+temp;
  rem=n*doubledouble(b.l());
  rem=temp-rem;
  return rem;
}

inline doubledouble sin(const doubledouble& x) {
  static const doubledouble tab[9]={ // tab[b] := sin(b*Pi/16)...
    0.0,
    "0.1950903220161282678482848684770222409277",
    "0.3826834323650897717284599840303988667613",
    "0.5555702330196022247428308139485328743749",
    "0.7071067811865475244008443621048490392850",
    "0.8314696123025452370787883776179057567386",
    "0.9238795325112867561281831893967882868225",
    "0.9807852804032304491261822361342390369739",
    1.0
  };
  static const doubledouble sinsTab[7] = { // Chebyshev coefficients
    "0.9999999999999999999999999999993767021096",
    "-0.1666666666666666666666666602899977158461",
    "8333333333333333333322459353395394180616.0e-42",
    "-1984126984126984056685882073709830240680.0e-43",
    "2755731922396443936999523827282063607870.0e-45",
    "-2505210805220830174499424295197047025509.0e-47",
    "1605649194713006247696761143618673476113.0e-49"
  };
  if (fabs(x.h())<1.0e-7) return x*(1.0-sqr(x)/3);
  int a,b; doubledouble sins, coss, k1, k3, t2, s, s2, sinb, cosb;
  // reduce x: -Pi < x <= Pi
  k1=x/doubledouble::TwoPi(); k3=k1-rint(k1);
  // determine integers a and b to minimize |s|, where  s=x-a*Pi/2-b*Pi/16
  t2=4*k3;
  a=int(rint(t2));
  b=int(rint((8*(t2-a))));   // must have |a|<=2 and |b|<=7 now
  s=doubledouble::Pi()*(k3+k3-(8*a+b)/16.0); // s is now reduced argument. -Pi/32 < s < Pi/32
  s2=sqr(s); 
  // Chebyshev series on -Pi/32..Pi/32, max abs error 2^-98=3.16e-30, whereas
  // power series has error 6e-28 at Pi/32 with terms up to x^13 included.
  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 as -Pi/32 < s < Pi/32
  // here sinb=sin(b*Pi/16) etc.
  sinb= (b>=0) ? tab[b] : -tab[-b]; cosb=tab[8-abs(b)];
  if (0==a) {
    return  sins*cosb+coss*sinb;
  } else if (1==a) {
    return -sins*sinb+coss*cosb;
  } else if (-1==a) {
    return  sins*sinb-coss*cosb;
  } else  { // |a|=2
    return -sins*cosb-coss*sinb;
  }
  // i.e. return sins*(cosa*cosb-sina*sinb)+coss*(sina*cosb+cosa*sinb);
}

// sin and cos.   Faster than separate calls of sin and cos.
inline void sincos(const doubledouble x, doubledouble& sinx, doubledouble& cosx) { 
  static const doubledouble tab[9]={ // tab[b] := sin(b*Pi/16)...
    0.0,
    "0.1950903220161282678482848684770222409277",
    "0.3826834323650897717284599840303988667613",
    "0.5555702330196022247428308139485328743749",
    "0.7071067811865475244008443621048490392850",
    "0.8314696123025452370787883776179057567386",
    "0.9238795325112867561281831893967882868225",
    "0.9807852804032304491261822361342390369739",
    1.0
  };
  static const doubledouble sinsTab[7] = { // Chebyshev coefficients
    "0.9999999999999999999999999999993767021096",
    "-0.1666666666666666666666666602899977158461",
    "8333333333333333333322459353395394180616.0e-42",
    "-1984126984126984056685882073709830240680.0e-43",
    "2755731922396443936999523827282063607870.0e-45",
    "-2505210805220830174499424295197047025509.0e-47",
    "1605649194713006247696761143618673476113.0e-49"
  };