/*Copyright (C) 1997, 1998, 1999, 2000 Victor ShoupThis program is free software; you can redistribute it and/ormodify it under the terms of the GNU General Public Licenseas published by the Free Software Foundation; either version 2of 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 ofMERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See theGNU General Public License for more details.You should have received a copy of the GNU General Public Licensealong with this program; if not, write to the Free SoftwareFoundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.*****************************************************The quad_float package is derived from the doubledouble package ofKeith Briggs.  However, the version employed in NTL has been extensively modified.  Below, I attach the copyright notice from the originaldoubledouble package, which is currently available at    http://www.labs.bt.com/people/briggsk2*****************************************************Copyright (C) 1997 Keith Martin BriggsExcept where otherwise indicated,this program is free software; you can redistribute it and/ormodify it under the terms of the GNU General Public Licenseas published by the Free Software Foundation; either version 2of 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 ofMERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See theGNU General Public License for more details.You should have received a copy of the GNU General Public Licensealong with this program; if not, write to the Free SoftwareFoundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.*/#include <NTL/quad_float.h>#include <NTL/RR.h>#include <float.h>#include <NTL/new.h>NTL_START_IMPL#if (defined(__GNUC__) && NTL_EXT_DOUBLE && (defined(__i386__) || defined(__i486__) || defined(__i586__)))#if (!defined(NTL_X86_FIX) && !defined(NTL_NO_X86_FIX))#define NTL_X86_FIX#endif#endif#if (NTL_EXT_DOUBLE && !defined(NTL_X86_FIX))#define DOUBLE volatile double#else#define DOUBLE double#endif#ifdef NTL_X86_FIX#define START_FIX \  unsigned short __old_cw, __new_cw; \  asm volatile ("fnstcw %0":"=m" (__old_cw)); \  __new_cw = (__old_cw & ~0x300) | 0x200; \  asm volatile ("fldcw %0": :"m" (__new_cw));#define END_FIX  asm volatile ("fldcw %0": :"m" (__old_cw));#else#define START_FIX#define END_FIX#endifstatic void normalize(quad_float& z, const double& xhi, const double& xlo){START_FIX   DOUBLE u, v;   u = xhi + xlo;   v = xhi - u;   v = v + xlo;   z.hi = u;   z.lo = v;END_FIX}#if (NTL_BITS_PER_LONG >= NTL_DOUBLE_PRECISION)quad_float to_quad_float(long n){START_FIX   DOUBLE xhi, xlo;   DOUBLE u, v;   xhi = double(n);   // Because we are assuming 2's compliment integer   // arithmetic, the following prevents long(xhi) from overflowing.   if (n > 0)      xlo = double(n+long(-xhi));   else      xlo = double(n-long(xhi));   // renormalize...just to be safe     u = xhi + xlo;   v = xhi - u;   v = v + xlo;END_FIX   return quad_float(u, v);}quad_float to_quad_float(unsigned long n){START_FIX   DOUBLE xhi, xlo, t;   DOUBLE u, v;   const double bnd = double(1L << (NTL_BITS_PER_LONG-2))*4.0;   xhi = double(n);      if (xhi >= bnd)      t = xhi - bnd;   else      t = xhi;   long llo = long(n - (unsigned long)(t));   xlo = double(llo);   // renormalize...just to be safe   u = xhi + xlo;   v = xhi - u;   v = v + xlo;END_FIX   return quad_float(u, v);}#endiflong quad_float::oprec = 10;void quad_float::SetOutputPrecision(long p){   if (p < 1) p = 1;   if (p >= (1L << (NTL_BITS_PER_LONG-4)))       Error("quad_float: output precision too big");   oprec = p;}quad_float operator +(const quad_float& x, const quad_float& y ) {START_FIX        DOUBLE    H, h, T, t, S, s, e, f;        DOUBLE    t1;        S = x.hi + y.hi;        T = x.lo + y.lo;        e = S - x.hi;        f = T - x.lo;        t1 = S-e;        t1 = x.hi-t1;        s = y.hi-e;        s = s + t1;                t1 = T-f;        t1 = x.lo-t1;        t = y.lo-f;        t = t + t1;        s = s + T;        H = S + s;        h = S - H;        h = h + s;        h = h + t;        e = H + h;         f = H - e;        f = f + h;END_FIX        return quad_float(e, f);}quad_float& operator +=(quad_float& x, const quad_float& y ) {START_FIX        DOUBLE    H, h, T, t, S, s, e, f;        DOUBLE    t1;        S = x.hi + y.hi;        T = x.lo + y.lo;        e = S - x.hi;        f = T - x.lo;        t1 = S-e;        t1 = x.hi-t1;        s = y.hi-e;        s = s + t1;                t1 = T-f;        t1 = x.lo-t1;        t = y.lo-f;        t = t + t1;        s = s + T;        H = S + s;        h = S - H;        h = h + s;        h = h + t;        e = H + h;         f = H - e;        f = f + h;        x.hi = e;        x.lo = f;END_FIX        return x;}quad_float operator -(const quad_float& x, const quad_float& y ) {START_FIX        DOUBLE    H, h, T, t, S, s, e, f;        DOUBLE    t1, yhi, ylo;        yhi = -y.hi;        ylo = -y.lo;        S = x.hi + yhi;        T = x.lo + ylo;        e = S - x.hi;        f = T - x.lo;        t1 = S-e;        t1 = x.hi-t1;        s = yhi-e;        s = s + t1;                t1 = T-f;        t1 = x.lo-t1;        t = ylo-f;        t = t + t1;        s = s + T;        H = S + s;        h = S - H;        h = h + s;        h = h + t;        e = H + h;         f = H - e;        f = f + h;END_FIX        return quad_float(e, f);}quad_float& operator -=(quad_float& x, const quad_float& y ) {START_FIX	DOUBLE    H, h, T, t, S, s, e, f;        DOUBLE    t1, yhi, ylo;        yhi = -y.hi;        ylo = -y.lo;        S = x.hi + yhi;        T = x.lo + ylo;        e = S - x.hi;        f = T - x.lo;        t1 = S-e;        t1 = x.hi-t1;        s = yhi-e;        s = s + t1;                t1 = T-f;        t1 = x.lo-t1;        t = ylo-f;        t = t + t1;        s = s + T;        H = S + s;        h = S - H;        h = h + s;        h = h + t;        e = H + h;         f = H - e;        f = f + h;        x.hi = e;        x.lo = f;END_FIX        return x;}quad_float operator -(const quad_float& x){START_FIX   DOUBLE xhi, xlo, u, v;   xhi = -x.hi;   xlo = -x.lo;   // it is a good idea to renormalize here, just in case   // the rounding rule depends on sign, and thus we will   // maintain the "normal form" for quad_float's.     u = xhi + xlo;   v = xhi - u;   v = v + xlo;END_FIX   return quad_float(u, v);}quad_float operator *(const quad_float& x,const quad_float& y ) {START_FIX  DOUBLE hx, tx, hy, ty, C, c;  DOUBLE t1, t2;  C = NTL_QUAD_FLOAT_SPLIT*x.hi;  hx = C-x.hi;  c = NTL_QUAD_FLOAT_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);    t1 = hx*hy;  t1 = t1-C;  t2 = hx*ty;  t1 = t1+t2;  t2 = tx*hy;  t1 = t1+t2;  t2 = tx*ty;  c = t1+t2;  t1 = x.hi*y.lo;  t2 = x.lo*y.hi;  t1 = t1+t2;  c = c + t1;  hx = C+c;  tx = C-hx;  tx = tx+c;END_FIX  return quad_float(hx, tx);}quad_float& operator *=(quad_float& x,const quad_float& y ) {START_FIX  DOUBLE hx, tx, hy, ty, C, c;  DOUBLE t1, t2;  C = NTL_QUAD_FLOAT_SPLIT*x.hi;  hx = C-x.hi;  c = NTL_QUAD_FLOAT_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);    t1 = hx*hy;  t1 = t1-C;  t2 = hx*ty;  t1 = t1+t2;  t2 = tx*hy;  t1 = t1+t2;  t2 = tx*ty;  c = t1+t2;  t1 = x.hi*y.lo;  t2 = x.lo*y.hi;  t1 = t1+t2;  c = c + t1;  hx = C+c;  tx = C-hx;  tx = tx+c;  x.hi = hx;  x.lo = tx;END_FIX  return x;}quad_float operator /(const quad_float& x, const quad_float& y ) {START_FIX  DOUBLE hc, tc, hy, ty, C, c, U, u;  DOUBLE t1;  C = x.hi/y.hi;  c = NTL_QUAD_FLOAT_SPLIT*C;  hc = c-C;  u = NTL_QUAD_FLOAT_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;  u = hc*hy;  u = u-U;  t1 = hc*ty;  u = u+t1;  t1 = tc*hy;  u = u+t1;  t1 = tc*ty;  u = u+t1;