quad_float.c

密码大家Shoup写的数论算法c语言实现

/*
Copyright (C) 1997, 1998, 1999, 2000 Victor Shoup

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.

*****************************************************

The quad_float package is derived from the doubledouble package of
Keith Briggs.  However, the version employed in NTL has been extensively 
modified.  Below, I attach the copyright notice from the original
doubledouble package, which is currently available at 

   http://www.labs.bt.com/people/briggsk2

*****************************************************

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 <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

#endif


static 
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);
}
#endif


long 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;