// quad.h
// Keith Briggs.   Last revised 96 Mar 6
// Version 1.9 for g++ 2.7.1
// Seems to work with g++ -O3
// However, g++ 2.5.8 needs g++ -O3 -ffloat-store to work properly
// Change log:  See quad.cc
// 
// TO DO
//  New name DoubledDouble?
//  Do something about global constants.   Make them static data members of Quad?
//  As it is, quad.h cannot be included unless quad.o is also linked.
//  Perhaps name them Quad_Pi etc.
//  Perhaps define them by a struc of two doubles.
//  All mixed ops Quad, double.   Double args do not need to be const double&.
//  double modf    (double x, double *n);
//  double trunc   (double x);

#define QUAD_INLINE
#define DEBUG_QUAD 0

#ifdef __GNUC__
  #if (__GNUC_MINOR__<=6)
    #define bool int
    #define explicit
  #endif
#endif

/*
C++ functions for quad (i.e. double+double) precision.

These functions use techniques due to Dekker, Linnainmaa, Kahan, Knuth 
and Priest.  I credit Kahan with the addition algorithm (the simplification
which permits the elimination of the tests and branches is due to Knuth);
Dekker and Linnainmaa with the multiply, divide, and square root routines,
and Priest for the initial transcription into C++.

See:

T. J. Dekker
   Point Technique for Extending the Available Precision,
   Numer. Math. 18 (1971), pp. 224-242.
S. Linnainmaa
   Software for doubled-precision floating point computations
   ACM TOMS 7, 172-283 (1081).
D. Priest 
  On properties of floating point arithmetics: numerical stability
  and the cost of accurate computations.   Ph.D. Diss, Berkeley 1992.

NB.   These functions assume IEEE double!
Note that all functions returns values such that |x.lo|<|x.hi|.

*/

#ifndef __QUAD_H__
#define __QUAD_H__
#include <iostream.h>
#include <iomanip.h>
#include <math.h>
#include <stdlib.h>

static const double Split=134217729.0L;   // 2^27 + 1, appropriate for IEEE double

class Quad {
protected:
  double hi, lo;
public:
  // Public access to hi and lo
  inline double h() const;
  inline double l() const;

  // Constructors
  inline Quad();
  inline Quad(int n) { hi=double(n); lo=0.0; };
  inline Quad(double x) { hi=x; lo=0.0; };
  inline Quad(double, double);
  inline Quad(const Quad&);
  Quad(char*);

  // Operators
  Quad& operator=(const Quad&);
  Quad& operator=(const double&);
  Quad& operator=(const int&);  // Get funny errors without this
  Quad& operator=(char*);
  Quad normalize(const Quad&) { 
    double h=hi+lo; return Quad(h,lo+double(hi-h)); };
  // Type conversion operator.  Not really necessary...
  operator double() { return hi+lo; };
  friend Quad operator +(const Quad&, const Quad& );
  friend Quad operator +(const double&, const Quad& );
  friend Quad operator -(const Quad&, const Quad& );
  friend Quad operator -(const double&, const Quad& );
  inline friend Quad operator -(const Quad& x);  // Unary -
  friend Quad operator *(const Quad&, const Quad& );
  friend Quad operator *(const double&, const Quad& );
  friend Quad operator *(const int&, const Quad& );
  friend Quad operator /(const Quad&, const Quad& );
  friend Quad operator /(const Quad&, const double& );
  friend Quad operator /(const double&, const Quad& );
  friend Quad operator /(const Quad&, const int& );
  Quad& operator +=(const Quad&);
  Quad& operator -=(const Quad&);
  Quad& operator *=(const Quad&);
  Quad& operator /=(const Quad&);
  //friend Quad operator %(const Quad&, const Quad& );
  friend Quad operator |(const Quad&, const Quad& );
  friend double dnorm(const Quad&);
  friend int intq(const Quad&);
  // member functions
  void dump(char*);  // debugging use only
};  // end class Quad

void base_and_prec(void);
Quad atoq(const char *);  // string to Quad conversion
bool operator> (const Quad&, const Quad&);
bool operator>=(const Quad&, const Quad&);
bool operator< (const Quad&, const Quad&);
bool operator<=(const Quad&, const Quad&);
bool operator==(const Quad&, const Quad&);
bool operator!=(const Quad&, const Quad&);

// Useful constants
static const Quad   Log2="0.6931471805599453094172321214581765680755"; 
static const Quad  Log10="2.302585092994045684017991454684364207601"; 
static const Quad     Pi="3.1415926535897932384626433832795028841972";
static const Quad  TwoPi="6.2831853071795864769252867665590057683943";
static const Quad  Pion2="1.5707963267948966192313216916397514420985";
static const Quad  Pion4="0.7853981633974483096156608458198757210493";

// inline members

inline double Quad::h() const {return hi;}
inline double Quad::l() const {return lo;}
inline Quad::Quad():hi(0.0),lo(0.0) {}
inline Quad::Quad(double x, double y):hi(x),lo(y) {}
// or inline Quad::Quad(double x, double y=0.0) { hi=x+y; lo=y+double(x-hi); } //normalize
inline Quad::Quad(const Quad& x):hi(x.h()),lo(x.l()) {}
inline Quad& Quad::operator=(const Quad& x){ hi=x.hi; lo=x.lo; return *this;}
inline Quad& Quad::operator=(const double& x){ hi=x; lo=0.0; return *this;}
inline Quad& Quad::operator=(const int& x){ hi=x; lo=0.0; return *this;}

// inline functions
inline Quad operator-(const Quad& x) {return Quad(-x.hi,-x.lo);};  // Unary -
inline double dnorm(const Quad& x){ return fabs(x.h());}
inline int intq(const Quad& x){      // explicit type conversion Quad -> int
  return int(x.hi); 
} 

// inline functions (defined in quad.cc)
inline Quad Qcopysign(const Quad&, const double);

// non inline functions (defined in quad.cc)

istream& operator >> (istream&, Quad&);
ostream& operator << (ostream&, const Quad&);
int sign(const Quad&);
Quad hypot(const Quad&, const Quad&);
Quad recip(const Quad&);
Quad sqrt(const Quad&);
Quad sqr(const Quad&);
Quad cub(const Quad&);
Quad sqr_double(const double&);
Quad rint(const Quad&);
Quad floor(const Quad&);
Quad trunc(const Quad&);
Quad fmod(const Quad&, const int);
Quad fabs(const Quad&);
Quad exp(const Quad&); 
Quad log(const Quad&);  
Quad log10(const Quad&);  
Quad powint(const Quad&, const int);
Quad pow(const Quad&, const Quad&);
Quad sin(const Quad&);
void sincos(const Quad x, Quad& sinx, Quad& cosx);
Quad cos(const Quad&);
Quad atan(const Quad&);
Quad atan2(const Quad&, const Quad&);
Quad asin(const Quad&);
Quad sinh(const Quad&);
Quad cosh(const Quad&);
Quad tanh(const Quad&);
int  digits(const Quad&,const Quad&);

#ifdef QUAD_INLINE

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

// Addition
/*      (C) W. Kahan 1989
*       NOTICE:
*       Copyrighted programs may not be translated, used, nor
*       reproduced without the author's permission.  Normally that
*       permission is granted freely for academic and scientific
*       purposes subject to the following three requirements:
*       1. This NOTICE and the copyright notices must remain attached
*          to the programs and their translations.
*       2. Users of such a program should regard themselves as voluntary
*          participants in the author's researches and so are obliged
*          to report their experience with the program back to the author.
*       3. Neither the author nor the institution that employs him will
*          be held responsible for the consequences of using a program
*          for which neither has received payment.
*       Would-be commercial users of these programs must correspond
*       with the author to arrange terms of payment and warranty.
*/

//inline 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;
//}

inline Quad operator +(const Quad& x,const Quad& y ) { // KMB 95Nov30
  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); 
  e = s+T;
  H = S+e; h = e+(S-H);
  e = t+h;
  z.hi = H+e; z.lo = e+(H-z.hi); 
  return z;
}

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

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

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

inline Quad& Quad::operator /=(const Quad& y) {
  Quad z; double hc, tc, hy, ty, C, c, U, u;
  C = 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 = ((((hi-U)-u)+lo)-C*y.lo)/y.hi;
  u = C+c; hi = u; lo = double(C-u)+c; return *this;
}

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

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

// Multiplication
inline 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
inline 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
inline 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 using Newton's method (EXPERIMENTAL)
inline 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
inline 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
inline 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 ?? BUGGY? KMB95Apr5
inline 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
inline 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;
}

#endif // QUAD_INLINE

#endif