// quad.cc 
// Version 1.9 for g++ 2.7.1
// Keith Briggs.   Last revised 96 Mar 7
// TODO: 
//    copysign not Qcopysign
//    Check ieee 754 for more things to emulate
// C++ functions for quad (double+double) precision.
// Use with quad.h
// Change log:
//  95dec21 added base_and_prec.
//  95nov24 fixed sqrt(0)
//  95nov24 fixed operators <= and >=
//  95nov24 Fixed floor
//  96mar01 changes to sin and added sincos
//  96mar04 further changes to sin and added sincos
//  96mar06 added tanh.  Not tested yet.
//  96mar06 new version of exp using Chebyshev series.
//  96mar06 new version of atan using Newton method.
//  96mar07 debugged ostream& operator>>.

#include "quad.h"
#include <stdlib.h>
#include <nan.h>  // defines NAN, at least in gcc

// DOMAIN_ERROR=action to take on function domain errors 
#ifdef NAN
#define DOMAIN_ERROR return(Quad(NAN,0.0))
#else
#define DOMAIN_ERROR exit(1)
#endif

void base_and_prec(void) { // Linnainmaa ACM TOMS 7, 272 Thm 3
  double U,R,u,r,beta; int p;
  U=4.0/3.0;  U-=1.0; U*=3.0; U-=1.0; U=fabs(U);
  R=U/2.0+1.0; R-=1.0;
  if (R!=0.0) U=R;
  u=2.0/3.0;  u-=0.5; u*=3.0; u-=0.5; u=fabs(u);
  r=u/2.0+0.5; r-=0.5;
  if (r!=0.0) u=r;
  beta=U/u; p=int(-log(u)/log(beta)+0.5);
  cout<<"Type double: base is "<<beta<<",  precision is "<<p<<endl;
}

// number of decimal digits to which x and y agree
int digits(const Quad& x, const Quad& y){ 
 Quad diff=fabs(x-y);
 if (diff.h()==0.0) return 32;
 int d=-intq(floor((0.4*log((diff/fabs(x))))).h());
 return d<32?d:32;
}

inline double copysign(const double, const double); // ieee
inline double copysign(const double x, const double y){
  if (y>=0.0) return fabs(x); else return -fabs(x); }

inline Quad Qcopysign(const Quad& x, const double y){
  if (y>=0.0) return fabs(x); else return -fabs(x); }

Quad atoq(const char *s) {
  Quad result = 0.0; int n, sign, ex = 0;
  /* eat whitespace */ while (*s == ' ' || *s == '\t' || *s == '\n') s++;
  switch (*s) { // get sign of mantissa
    case '-': { sign = -1;  s++; break; }
    case '+': s++;  // no break 
    default: sign = 1;
  }
  /* get digits before decimal point */
  while (n=(*s++)-'0', n>=0 && n<10) result = 10.0*result+n;
  s--;
  if (*s == '.')   /* get digits after decimal point */ {
    s++;
    while (n=(*s++)-'0', n>=0 && n<10) { result = 10.0*result+n; --ex; }
    s--;
  }
  if (*s == 'e' || *s == 'E') /* get exponent */ ex+=atoi(++s);
  /* exponent adjustment */
  // cerr<<"atoq: result="<<endl; result.dump(""); cerr<<"atoq: ex="<<ex<<endl; 
  while (ex-- > 0) result *= 10;
  while (++ex < 0) result /= 10;
  return (sign>=0) ? result : -result;
}

void Quad::dump(char* s=""){  // debugging use only
  cerr<<s<<setprecision(16)<<"Quad("<<hi<<","<<lo<<")="
         <<setprecision(32)<<Quad(hi,lo)<<endl;
}

// Constructor from string
Quad::Quad(char* s){ Quad y=atoq(s); hi=y.hi; lo=y.lo; }

// Quad = string
Quad& Quad::operator=(char* s){ Quad y=atoq(s); hi=y.hi; lo=y.lo; return *this; }

// improved 96Mar7...
istream& operator >> (istream& s, Quad& x){
  Quad result=0.0; int n, sign=1, ex=0; char ch;
  s>>ws; // eat whitespace
  s>>ch;
  if (ch=='-') { sign = -1;  s>>ch; } else if (ch=='+') { s>>ch; }
  // get digits before decimal point
  n=ch-'0';
  while (n>=0 && n<10) { 
    result = 10*result+n; 
    s.get(ch);  // cannot use s>>ch; as that will eat whitespace
    if (ch == ' ' || ch == '\t' || ch == '\n') goto fixup;
    n=ch-'0'; 
  }
  if (ch=='.') {  // get digits after decimal point
     s>>ch; n=ch-'0';
     while (n>=0 && n<10) { 
       result = 10*result+n; 
       s.get(ch);
       n=ch-'0'; --ex; 
       if (ch == ' ' || ch == '\t' || ch == '\n') goto fixup;
     }
  }
  n=0;
  if (ch == 'e' || ch == 'E') { s>>n; ex+=n; } // get exponent
    else s.putback(ch);
  fixup:
  if (sign<0) result = -result;
  // exponent adjustment
  while (ex-- > 0) result *= 10;
  while (++ex < 0) result /= 10;
  x = result;
  return s;
}

// output
ostream& operator << (ostream& s, const Quad& x){
  if (x.h()==0.0) { s << " 0.0" << " "; return s; }
  int Digits=s.precision();
  Quad ten=10.0,y=fabs(x); double q=log10(y.h());
  int m,n=int(floor(q)); 
  if (n<0) n++;
  Quad l = powint(ten,n);
  y = y/l;
  if (sign(x)<0) s<<"-"; else s<<" ";
  int d = Digits>34 ? 34 : Digits;
  d = d<3 ? 3 : d;
  for (int i=1; i<=d; i++) {
    if (2==i) s<<".";
    m = int(floor(y.h())); 
    s<<m;
    y = (y-Quad(m))*ten;
    if (y.h()<0.0) break; // x must be an integer
    //if (y.h()==0.0) break; // ???????????
  }
  if (n!=0) s<<"e"<<n; else s<<"";
  s << "";
  return s;
}

// rint (round to nearest int)
Quad rint(const Quad& x){ return floor(x+Quad(0.5)); }

// Floor.  Changed 95Nov24.  See Graham, Knuth and Patashnik page 70.
// Greatest integer <= x
Quad floor(const Quad& x) { 
  double fh=floor(x.h()), fl=floor(x.l()); 
  Quad t=(x.h()-fh)+(x.l()-fl);  // t={hi}+{lo}
  if (t<1.0) return Quad(fh,fl);
        else return Quad(fh,fl)+1.0;
}

// ceil.  Added 95Nov24.
// Least integer >= x
// ceil(x)=-floor(-x)
Quad ceil(const Quad& x) { 
  double fh=floor(-x.h()), fl=floor(-x.l()); 
  Quad t=(x.h()-fh)+(x.l()-fl);  // t={hi}+{lo}
  if (t<1.0) return -Quad(fh,fl);
        else return -Quad(fh,fl)-1.0;
}

// trunc.  Changed 95Nov24.
// Round towards zero.
Quad trunc(const Quad& x) { 
  if (x>=0.0) return floor(x); else return -floor(-x);
}

// Modulo
Quad fmod(const Quad& x, const int n) { return x-n*floor(x/n); }

// Signum
int sign(const Quad& x){
  if (x.h()>0.0) return 1; else if (x.h()<0.0) return -1; else return 0;
}

// Comparison
bool operator> (const Quad& x, const Quad& y) { 
   return (x.h()> y.h()) || (x.h()==y.h() && x.l()> y.l()); }
bool operator>=(const Quad& x, const Quad& y) { 
   return (x.h()>y.h()) || (x.h()==y.h() && x.l()>=y.l()); }
bool operator< (const Quad& x, const Quad& y) { 
   return (x.h()< y.h()) || (x.h()==y.h() && x.l()< y.l()); }
bool operator<=(const Quad& x, const Quad& y) { 
   return (x.h()<y.h()) || (x.h()==y.h() && x.l()<=y.l()); }
bool operator==(const Quad& x, const Quad& y) 
   { return x.h()==y.h() && x.l()==y.l(); }
bool operator!=(const Quad& x, const Quad& y) 
   { return x.h()!=y.h() || x.l()!=y.l(); }

#ifndef QUAD_INLINE

// Absolute value
Quad fabs(const Quad& x) { if (x.h()>=0.0) return x; else return -x; }

// Addition
Quad operator +(const Quad& x,const Quad& y ) {
  Quad z; double H, h, T, t, S, s, e, f;
  S = x.hi+y.hi;
  T = x.lo+y.lo;
  e = S-x.hi;
  f = T-x.lo;
  s = S-e; 
  t = T-f; 
  s = (y.hi-e)+(x.hi-s);
  t = (y.lo-f)+(x.lo-t); 
  H = S+(s+T);
  h = (s+T)+(S-H);
  z.hi = H+(t+h);
  z.lo = (t+h)+(H-z.hi);
  return z;
}

// Subtraction
Quad operator -(const Quad& x,const Quad& y) { return x+Quad(-y.hi,-y.lo); }

// Multiplication
Quad operator *(const Quad& x,const Quad& y ) {
  Quad z; double hx, tx, hy, ty, C, c;
  C = Split*x.hi;
  hx = C-x.hi;
  c = Split*y.hi;
  hx = C-hx; 
  tx = x.hi-hx;
  hy = c-y.hi; 
  C = x.hi*y.hi;
  hy = c-hy;
  ty = y.hi-hy;
  c = ((((hx*hy-C)+hx*ty)+tx*hy)+tx*ty)+(x.hi*y.lo+x.lo*y.hi);
  hx = C+c;
  z.hi = hx;
  z.lo = c+double(C-hx);
  return z;
}

// double*Quad
Quad operator *(const double& x,const Quad& y ) {
  Quad z; double hx, tx, hy, ty, C, c;
  C = Split*x;
  hx = C-x;
  c = Split*y.hi;
  hx = C-hx ; 
  tx = x-hx;
  hy = c-y.hi; 
  C = x*y.hi;
  hy = c-hy; 
  ty = y.hi - hy;
  c = ((((hx*hy-C)+hx*ty)+tx*hy)+tx*ty)+x*y.lo;
  hx = C+c;
  z.hi = hx;
  z.lo = c+double(C-hx); 
  return z;
}

// int*Quad
Quad operator *(const int& x,const Quad& y ) {
  if (x==-1) return -y;
  if (x==0) return Quad(0.0);
  if (x==1) return y;
  if (x==2) return y+y;
  Quad z; double hx, tx, hy, ty, C, c;
  C = Split*x;
  hx = C-x;
  c = Split*y.hi;
  hx = C-hx ; 
  tx = x-hx;
  hy = c-y.hi; 
  C = x*y.hi;
  hy = c-hy; 
  ty = y.hi-hy;
  c = ((((hx*hy-C)+hx*ty)+tx*hy)+tx*ty)+(x*y.lo);
  hx = C+c;
  z.hi = hx;
  z.lo = c+double(C-hx); 
  return z;
}

// Divide by Newton's method (EXPERIMENTAL)
/*
Quad operator |(const Quad& a,const Quad& b ) {
  double dx=1.0/b.hi;
  Quad x=dx+(1.0-dx*b.hi)*dx;  // Newton
  Quad y=a.hi*x.hi;
  return y+(a-y*b).hi*x.hi;
}
*/

// Division
Quad operator /(const Quad& x,const Quad& y ) {
  Quad z; double hc, tc, hy, ty, C, c, U, u;
  C = x.hi/y.hi;
  c = Split*C;
  hc =c-C;  
  u = Split*y.hi;
  hc = c-hc;
  tc = C-hc;
  hy = u-y.hi; 
  U = C * y.hi;
  hy = u-hy;
  ty = y.hi-hy;
  u = (((hc*hy-U)+hc*ty)+tc*hy)+tc*ty;
  c = ((((x.hi-U)-u)+x.lo)-C*y.lo)/y.hi;
  u = C+c;
  z.hi = u;
  z.lo = double(C-u)+c;
  return z;
}

// double/Quad
Quad operator /(const double& x,const Quad& y ) {
  Quad z;
  double  hc, tc, hy, ty, C, c, U, u;
  C = x/y.hi;
  c = Split*C;
  hc =c-C;  
  u = Split*y.hi;
  hc = c-hc; 
  tc = C-hc;
  hy = u-y.hi; 
  U = C*y.hi;
  hy = u-hy;
  ty = y.hi-hy;
  u = (((hc*hy-U)+hc*ty)+tc*hy)+tc*ty;
  c = ((((x-U)-u))-C*y.lo)/y.hi;
  u = C+c;
  z.hi = u;
  z.lo = double(C-u)+c;
  return z;
}

// Quad/double
Quad operator /(const Quad& x,const double& y ) {
  Quad z;  double hc, tc, hy, ty, C, c, U, u;
  C = x.hi/y;
  c = Split*C;
  hc = c-C;  
  u = Split*y;
  hc = c-hc; 
  tc = C-hc;
  hy = u-y; 
  U = C*y;
  hy = u-hy; 
  ty = y-hy;
  u = (((hc*hy-U)+hc*ty)+tc*hy)+tc*ty;
  c = (((x.hi-U)-u)+x.lo)/y;
  u = C+c;
  z.hi = u; 
  z.lo = double(C-u)+c; 
  return z;
}

// Quad/int
Quad operator /(const Quad& x,const int& y ) {
  Quad z; double hc, tc, hy, ty, C, c, U, u;
  C = x.hi/y;
  c = Split*C;
  hc =c-C;  
  u = Split*y;
  hc = c-hc; 
  tc = C-hc;
  hy = u-y; 
  U = C*y;
  hy = u-hy; 
  ty = y-hy;
  u = (((hc*hy-U)+hc*ty)+tc*hy)+tc*ty;
  c = ((((x.hi-U)-u)+x.lo))/y;
  u = C+c;
  z.hi = u;
  z.lo = double(C-u)+c; 
  return z;
}

// quad square of a double
Quad sqr_double(const double& x) {
  double hx, tx, hy, ty, C, c;
  C = Split*x; hx = C-double(C-x); tx = x-hx;
  C = x*x;
  c = (((hx*hx-C)+hx*tx)+tx*hx)+tx*tx;
  hx = C+c;
  return Quad(hx,c+double(C-hx));
}

#endif // ifndef QUAD_INLINE

// Square  (faster than x*x)
Quad sqr(const Quad& x) {
  Quad z; double hx, tx, C, c;
  C = Split*x.h(); hx = C-x.h(); hx = C-hx; tx = x.h()-hx;
  C = x.h()*x.h();
  c = ((((hx*hx-C)+2.0*hx*tx))+tx*tx)+2.0*x.h()*x.l();
  hx = C+c;
  return Quad(hx,c+double(C-hx));
}

// cube
Quad cub(const Quad& x) { Quad z=x*sqr(x); return z; }

// hypot - crude version.  Should be improved to avoid overflow!!
Quad hypot(const Quad& x, const Quad& y) { return sqrt(sqr(x)+sqr(y)); }

// square root
Quad sqrt(const Quad& y) {
  double c,p,q,hx,tx,u,uu,cc,hi=y.h();
  if (hi<0.0) {
    cerr << "\nQuad: Attempt to take sqrt of " << hi << endl;
    DOMAIN_ERROR;
  }
  if (0.0==hi) return Quad(0.0,0.0);
  c = sqrt(hi);
  p = Split*c; hx = double(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;
  return Quad(u,cc+(c-u));
}

/*
Quad sqrt(const Quad& y) { // Newton method for sqrt
  Quad x, p, half=0.5;
  x = 1.0/sqrt(y.h());
  p = x*y;
  return p+half*(y-p*p)*x;
}
*/

// Revised exponential function KMB 96 Mar 7
Quad exp(const Quad& t) {
  if (fabs(t.h())<1.0e-11) return 1.0+t+sqr(t)/2;
  int i, n=int(rint((t/Log2))); Quad r, s;
  r=(t-n*Log2)/256;  // Gives -log(2)/512 <= r <= m=log(2)/512 ~= 0.0013538030870
  // Use Chebyshev series on -m..m; max abs error 1e-32. 
  // Typical abs error on -m..m is about 1.5e-34
  // Power series would have error 4e-32 at m with terms up to x^8 included.
  s = 1.0+(
     "0.9999999999999999999999999999989068405105" + (
     "0.4999999999999999999999999999996963445872" + (
     "0.1666666666666666666666746193055972034254" + (
     "0.04166666666666666666666799210648336186960" + (
     "0.008333333333333317712551616807268910060964" + (
     "0.001388888888888886863972747322975862458609" + (
     "0.0001984127097766793471547633384037452748308" + 
     "0.00002480158856425184664607683931594069104335"*r)*r)*r)*r)*r)*r)*r)*r;
  r=s=powint(s,256);
  for (i=1; i<=abs(n); i++) s = s+s;
  if (n<0) return r*r*recip(s); else return s;
}

/*
// Exponential function using power series
// See Bailey, MPFUN  ACM TOMS vol 19, Algorithm 719.
Quad exp(const Quad& t) {
  int n=int(rint((t/Log2))); // minimizes |r|
  Quad r=(t-n*Log2)/256;
  Quad s=1.0+r; Quad term=0.5*r*r; double f=3.0;
  do { s += term; term *= r/f; f++; } while (fabs(term.h()/s.h())>1.0e-35);
  s=sqr(s); s=sqr(s); s=sqr(s); s=sqr(s); 
  s=sqr(s); s=sqr(s); s=sqr(s); s=sqr(s);
  r = s;
  for (int i=1; i<=abs(n); i++) s = s+s;
  if (n<0) { r=r*r*recip(s); return r; } else return s;
}
*/
// Natural logarithm
Quad log(const Quad& t) { // Newton method. See Bailey, MPFUN
  if (t.h() <= 0.0) {
    cerr << "Attempt to take Quad::log(<=0), quitting!\n";
    DOMAIN_ERROR;
  }
  Quad e, s=log(t.h()); // s=double approx to result
  e=exp(s);
  return s+(t-e)/e;  // Newton step
}

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

// Reciprocal
Quad recip(const Quad& y) {
  double  hc, tc, hy, ty, C, c, U, u;
  C = 1.0/y.h();
  c = Split*C;
  hc =c-C;  
  u = Split*y.h();
  hc = c-hc;
  tc = C-hc;
  hy = u-y.h(); 
  U = C*y.h();
  hy = u-hy;
  ty = y.h()-hy;
  u = (((hc*hy-U)+hc*ty)+tc*hy)+tc*ty;
  c = ((((1.0-U)-u))-C*y.l())/y.h();
  u = C+c;
  hy = double(C-u)+c; 
  return Quad(u,hy);
}

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

// Quad^int
Quad powint(const Quad& u, int c){
  Quad w; int i;  
  switch (c){
    case -2: return recip(u*u);  
    case -1: return recip(u);  
    case  0: return Quad(1.0);  // O^0 = NaN ??
    case  1: return u;  
    case  2: return sqr(u);
    case  3: return sqr(u)*u;
    default: { // binary method
      int n=c, m; Quad y(1.0), z=u; if (n<0) n=-n;
      do {
        m=n; n=n/2;
        if (n+n!=m) { // if m odd
           y = z*y; 
           if (n==0) { if (c>0) return y; else return recip(y); }
        }
        z = sqr(z); 
      } while (1);
    }
  }
}


//  sin by Bailey's method.    See mpfun article in TOMS 1993.
Quad sin(const Quad& x) { 
  const Quad 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
  };
  if (fabs(x.h())<1.0e-7) return x*(1.0-sqr(x)/3);
  int a,b; Quad sins, coss, k1, k3, t2, s, s2, sinb, cosb;
  // reduce x: -Pi < x <= Pi
  k1 = x/TwoPi; k3=k1-rint(k1);
  // determine integers a and b to minimise |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=Pi*(k3+k3-(8*a+b)/16.0); // s is now reduced argument. -Pi/32 < s < Pi/32
  s2=sqr(s); 
  /*
  // Old way, use power series...
  Quad term=s, sum=term, fac=6.0; int k=1;
  do { term*=-s2/fac; sum+=term; k++; fac=2*k*(k+k+1); } 
  while (fabs(term.h())>fabs(sum.h())*1.0e-35); Quad sins=sum;
  */
  // New way 96 Mar 4
  // 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*("0.9999999999999999999999999999993767021096"+
         ("-0.1666666666666666666666666602899977158461"+
         ("8333333333333333333322459353395394180616.0e-42"+
         ("-1984126984126984056685882073709830240680.0e-43"+
         ("2755731922396443936999523827282063607870.0e-45"+
         ("-2505210805220830174499424295197047025509.0e-47"+
          "1605649194713006247696761143618673476113.0e-49"*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.
void sincos(const Quad x, Quad& sinx, Quad& cosx) { 
  const Quad 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
  };
  if (fabs(x.h())<1.0e-11) { sinx=x; cosx=1.0-0.5*sqr(x); return; }
  int a,b; Quad sins, coss, k1, k3, t2, s, s2, sinb, cosb;
  k1=x/TwoPi; k3=k1-rint(k1);
  t2=4*k3;
  a=int(rint(t2));
  b=int(rint((8*(t2-a))));
  s=Pi*(k3+k3-(8*a+b)/16.0);
  s2=sqr(s);
  sins=s*("0.9999999999999999999999999999993767021096"+
         ("-0.1666666666666666666666666602899977158461"+
         ("8333333333333333333322459353395394180616.0e-42"+
         ("-1984126984126984056685882073709830240680.0e-43"+
         ("2755731922396443936999523827282063607870.0e-45"+
         ("-2505210805220830174499424295197047025509.0e-47"+
          "1605649194713006247696761143618673476113.0e-49"*s2)*s2)*s2)*s2)*s2)*s2);
  coss=sqrt(1.0-sqr(sins));
  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
Quad cos(const Quad& x){ return sin(Pion2-x); }

// hyperbolic
Quad sinh(const Quad& x){ Quad t=exp(x); return 0.5*(t-recip(t)); }
Quad cosh(const Quad& x){ Quad t=exp(x); return 0.5*(t+recip(t)); }
Quad tanh(const Quad z) {
  Quad e;
  if (z>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); }
}

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

/* Old way.  Power series converges too slowly.
// arctan
Quad atan(const Quad& x){
  if (x.h()==0.0) return Quad(0.0);
  Quad s=x; int fac=1;
  // use identity to reduce argument...
  while (fabs(s)>0.5) { s=s/(1.0+sqrt(1.0+s*s)); fac = fac+fac; }
  Quad  s2=s*s, term=s, sum=s, k=3.0;
  do {
    term *= -s2;
    sum += term/k;
    k = k+2.0;
    //if (k>100) { cerr << "No convergence in quad::atan !\n"; DOMAIN_ERROR;}
  } while (fabs(term.h())>fabs(sum.h())*1.0e-35);
  return fac*sum;
}
*/

Quad atan2(const Quad& qy, const Quad& 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(Pi, signy);
  if (x == zero) return Qcopysign(Pion2, signy);
/* __isinf not defined on landau...
  if (__isinf(x)) {
      if (__isinf(y)) return Qcopysign(signx==one ? Pion4 : 3.0*Pion4, signy);
      else            return Qcopysign(signx==one ? Quad(0.0) : Pi, signy);
  }
  if (__isinf(y)) return Qcopysign(Pion2, signy);
*/
  Quad aqy = fabs(qy);
  if (x<0.0) // X is negative.
    return Qcopysign(Pi-atan(aqy/(-qx)), signy);
  return Qcopysign(atan(aqy/qx), signy);
}

// arcsin
Quad asin(const Quad& x){
  if (fabs(x)>1.0) { cerr<<"|Argument|>1 in quad::asin!\n"; DOMAIN_ERROR;}
  return atan2(x,sqrt(1.0-x*x));
}
// end of quad.cc