g_lip_impl.h

NTL is a high-performance, portable C++ library providing data structures and algorithms for manipul

      if (direction == 0 && residual != 0) {
         if (residual == sgn)
            direction = 1;
         else 
            direction = -1;
      }

      if (direction == 0) {
         /* round to even */

         wh = wh << 1;

         /*
          * DIRT: if GMP has non-empty "nails", this won't work.
          */

         if (wh == 0) {
            wh = 1;
            bl++;
         }

         if (adata[bl] & wh)
            direction = 1;
         else
            direction = -1;
      }
   }
   else
      direction = -1;

   if (direction == 1)
      return sgn;

   return 0;
}




double _ntl_gdoub(_ntl_gbigint n)
{
   static _ntl_gbigint tmp = 0;

   long s;
   long shamt;
   long correction;
   double x;

   s = _ntl_g2log(n);
   shamt = s - NTL_DOUBLE_PRECISION;

   if (shamt <= 0)
      return _ntl_gdoub_aux(n);

   _ntl_grshift(n, shamt, &tmp);

   correction = _ntl_ground_correction(n, shamt, 0);

   if (correction) _ntl_gsadd(tmp, correction, &tmp);

   x = _ntl_gdoub_aux(tmp);

   x = _ntl_ldexp(x, shamt);

   return x;
}


double _ntl_glog(_ntl_gbigint n)
{
   static _ntl_gbigint tmp = 0;

   static double log_2;
   static long init = 0;

   long s;
   long shamt;
   long correction;
   double x;

   if (!init) {
      log_2 = log(2.0);
      init = 1;
   }

   if (_ntl_gsign(n) <= 0)
      ghalt("log argument <= 0");

   s = _ntl_g2log(n);
   shamt = s - NTL_DOUBLE_PRECISION;

   if (shamt <= 0)
      return log(_ntl_gdoub_aux(n));

   _ntl_grshift(n, shamt, &tmp);

   correction = _ntl_ground_correction(n, shamt, 0);

   if (correction) _ntl_gsadd(tmp, correction, &tmp);

   x = _ntl_gdoub_aux(tmp);

   return log(x) + shamt*log_2;
}





/* To implement _ntl_gdoubtoz, I've implemented essentially the
 * same algorithm as in LIP, processing in blocks of
 * NTL_SP_NBITS bits, rather than NTL_ZZ_NBITS.
 * This is conversion is rather delicate, and I don't want to
 * make any new assumptions about the underlying arithmetic.
 * This implementation should be quite portable. */

void _ntl_gdoubtoz(double a, _ntl_gbigint *xx)
{
   static _ntl_gbigint x = 0;

   long neg, i, t, sz;

   a = floor(a);

   if (!_ntl_IsFinite(&a))
      ghalt("_ntl_gdoubtoz: attempt to convert non-finite value");

   if (a < 0) {
      a = -a;
      neg = 1;
   }
   else
      neg = 0;

   if (a == 0) {
      _ntl_gzero(xx);
      return;
   }

   sz = 0;
   while (a >= 1) {
      a = a*(1.0/NTL_SP_FBOUND);
      sz++;
   }

   i = 0;
   _ntl_gzero(&x);

   while (a != 0) {
      i++;
      a = a*NTL_SP_FBOUND;
      t = (long) a;
      a = a - t;

      if (i == 1) {
         _ntl_gintoz(t, &x);
      }
      else {
         _ntl_glshift(x, NTL_SP_NBITS, &x);
         _ntl_gsadd(x, t, &x);
      }
   }

   if (i > sz) ghalt("bug in _ntl_gdoubtoz");

   _ntl_glshift(x, (sz-i)*NTL_SP_NBITS, xx);
   if (neg) _ntl_gnegate(xx);
}



/* I've adapted LIP's extended euclidean algorithm to
 * do rational reconstruction.  -- VJS.
 */


long 
_ntl_gxxratrecon(
   _ntl_gbigint ain,
   _ntl_gbigint nin,
   _ntl_gbigint num_bound,
   _ntl_gbigint den_bound,
   _ntl_gbigint *num_out,
   _ntl_gbigint *den_out
   )
{
   static _ntl_gbigint a = 0;
   static _ntl_gbigint n = 0;
   static _ntl_gbigint q = 0;
   static _ntl_gbigint w = 0;
   static _ntl_gbigint x = 0;
   static _ntl_gbigint y = 0;
   static _ntl_gbigint z = 0;
   static _ntl_gbigint inv = 0;
   static _ntl_gbigint u = 0;
   static _ntl_gbigint a_bak = 0;
   static _ntl_gbigint n_bak = 0;
   static _ntl_gbigint inv_bak = 0;
   static _ntl_gbigint w_bak = 0;

   mp_limb_t *p;

   long diff;
   long ilo;
   long sa;
   long sn;
   long snum;
   long sden;
   long e;
   long fast;
   long temp;
   long parity;
   long gotthem;
   long try11;
   long try12;
   long try21;
   long try22;
   long got11;
   long got12;
   long got21;
   long got22;

   double hi;
   double lo;
   double dt;
   double fhi, fhi1;
   double flo, flo1;
   double num;
   double den;
   double dirt;

   if (_ntl_gsign(num_bound) < 0)
      ghalt("rational reconstruction: bad numerator bound");

   if (!num_bound)
      snum = 0;
   else
      snum = SIZE(num_bound);

   if (_ntl_gsign(den_bound) <= 0)
      ghalt("rational reconstruction: bad denominator bound");

   sden = SIZE(den_bound);

   if (_ntl_gsign(nin) <= 0)
      ghalt("rational reconstruction: bad modulus");

   if (_ntl_gsign(ain) < 0 || _ntl_gcompare(ain, nin) >= 0)
      ghalt("rational reconstruction: bad residue");

      
   e = SIZE(nin);

   _ntl_gsetlength(&a, e);
   _ntl_gsetlength(&n, e);
   _ntl_gsetlength(&q, e);
   _ntl_gsetlength(&w, e);
   _ntl_gsetlength(&x, e);
   _ntl_gsetlength(&y, e);
   _ntl_gsetlength(&z, e);
   _ntl_gsetlength(&inv, e);
   _ntl_gsetlength(&u, e);
   _ntl_gsetlength(&a_bak, e);
   _ntl_gsetlength(&n_bak, e);
   _ntl_gsetlength(&inv_bak, e);
   _ntl_gsetlength(&w_bak, e);

   fhi1 = 1.0 + ((double) 32.0)/NTL_FDOUBLE_PRECISION;
   flo1 = 1.0 - ((double) 32.0)/NTL_FDOUBLE_PRECISION;

   fhi = 1.0 + ((double) 8.0)/NTL_FDOUBLE_PRECISION;
   flo = 1.0 - ((double) 8.0)/NTL_FDOUBLE_PRECISION;

   _ntl_gcopy(ain, &a);
   _ntl_gcopy(nin, &n);

   _ntl_gone(&inv);
   _ntl_gzero(&w);

   while (1)
   {
      if (SIZE(w) >= sden && _ntl_gcompare(w, den_bound) > 0) break;
      if (SIZE(n) <= snum && _ntl_gcompare(n, num_bound) <= 0) break;

      _ntl_gcopy(a, &a_bak);
      _ntl_gcopy(n, &n_bak);
      _ntl_gcopy(w, &w_bak);
      _ntl_gcopy(inv, &inv_bak);

      gotthem = 0;
      sa = SIZE(a);
      sn = SIZE(n);
      diff = sa - sn;
      if (!diff || diff == 1)
      {
         sa = SIZE(a);
         p = DATA(a) + (sa-1);
         num = (double) (*p) * NTL_ZZ_FRADIX;
         if (sa > 1)
            num += (*(--p));
         num *= NTL_ZZ_FRADIX;
         if (sa > 2)
            num += (*(p - 1));

         sn = SIZE(n);
         p = DATA(n) + (sn-1);
         den = (double) (*p) * NTL_ZZ_FRADIX;
         if (sn > 1)
            den += (*(--p));
         den *= NTL_ZZ_FRADIX;
         if (sn > 2)
            den += (*(p - 1));

         hi = fhi1 * (num + 1.0) / den;
         lo = flo1 * num / (den + 1.0);
         if (diff > 0)
         {
            hi *= NTL_ZZ_FRADIX;
            lo *= NTL_ZZ_FRADIX;
         }

         try11 = 1;
         try12 = 0;
         try21 = 0;
         try22 = 1;
         parity = 1;
         fast = 1; 
         while (fast > 0)
         {
            parity = 1 - parity;
            if (hi >= NTL_SP_BOUND)
               fast = 0;
            else
            {
               ilo = (long)lo;
               dirt = hi - ilo;
               if (dirt < 1.0/NTL_FDOUBLE_PRECISION || !ilo || ilo < (long)hi)
                  fast = 0;
               else
               {
                  dt = lo-ilo;
                  lo = flo / dirt;
                  if (dt > 1.0/NTL_FDOUBLE_PRECISION)
                     hi = fhi / dt;
                  else
                     hi = NTL_SP_BOUND;
                  temp = try11;
                  try11 = try21;
                  if ((NTL_WSP_BOUND - temp) / ilo < try21)
                     fast = 0;
                  else
                     try21 = temp + ilo * try21;
                  temp = try12;
                  try12 = try22;
                  if ((NTL_WSP_BOUND - temp) / ilo < try22)
                     fast = 0;
                  else
                     try22 = temp + ilo * try22;
                  if ((fast > 0) && (parity > 0))
                  {
                     gotthem = 1;
                     got11 = try11;
                     got12 = try12;
                     got21 = try21;
                     got22 = try22;
                  }
               }
            }
         }
      }
      if (gotthem)
      {
         _ntl_gsmul(inv, got11, &x);
         _ntl_gsmul(w, got12, &y);
         _ntl_gsmul(inv, got21, &z);
         _ntl_gsmul(w, got22, &w);
         _ntl_gadd(x, y, &inv);
         _ntl_gadd(z, w, &w);
         _ntl_gsmul(a, got11, &x);
         _ntl_gsmul(n, got12, &y);
         _ntl_gsmul(a, got21, &z);
         _ntl_gsmul(n, got22, &n);
         _ntl_gsub(x, y, &a);
         _ntl_gsub(n, z, &n);
      }
      else
      {
         _ntl_gdiv(a, n, &q, &a);
         _ntl_gmul(q, w, &x);
         _ntl_gadd(inv, x, &inv);
         if (!ZEROP(a))
         {
            _ntl_gdiv(n, a, &q, &n);
            _ntl_gmul(q, inv, &x);
            _ntl_gadd(w, x, &w);
         }
         else
         {
            break;
         }
      }
   }

   _ntl_gcopy(a_bak, &a);
   _ntl_gcopy(n_bak, &n);
   _ntl_gcopy(w_bak, &w);
   _ntl_gcopy(inv_bak, &inv);

   _ntl_gnegate(&w);

   while (1)
   {
      sa = SIZE(w);
      if (sa < 0) SIZE(w) = -sa;
      if (SIZE(w) >= sden && _ntl_gcompare(w, den_bound) > 0) return 0;
      SIZE(w) = sa;

      if (SIZE(n) <= snum && _ntl_gcompare(n, num_bound) <= 0) break;
      
      fast = 0;
      sa = SIZE(a);
      sn = SIZE(n);
      diff = sa - sn;
      if (!diff || diff == 1)
      {
         sa = SIZE(a);
         p = DATA(a) + (sa-1);
         num = (double) (*p) * NTL_ZZ_FRADIX;
         if (sa > 1)
            num += (*(--p));
         num *= NTL_ZZ_FRADIX;
         if (sa > 2)
            num += (*(p - 1));

         sn = SIZE(n);
         p = DATA(n) + (sn-1);
         den = (double) (*p) * NTL_ZZ_FRADIX;
         if (sn > 1)
            den += (*(--p));
         den *= NTL_ZZ_FRADIX;
         if (sn > 2)
            den += (*(p - 1));

         hi = fhi1 * (num + 1.0) / den;
         lo = flo1 * num / (den + 1.0);
         if (diff > 0)
         {
            hi *= NTL_ZZ_FRADIX;
            lo *= NTL_ZZ_FRADIX;
         }

         if (hi < NTL_SP_BOUND)
         {
            ilo = (long)lo;
            if (ilo == (long)hi)
               fast = 1;
         }
      }

      if (fast) 
      {
         if (ilo != 0) {
            if (ilo == 1) {
               _ntl_gsub(inv, w, &inv);
               _ntl_gsubpos(a, n, &a);
            }
            else {
               _ntl_gsmul(w, ilo, &x);
               _ntl_gsub(inv, x, &inv);
               _ntl_gsmul(n, ilo, &x);
               _ntl_gsubpos(a, x, &a);
            }
         }
      }
      else {
         _ntl_gdiv(a, n, &q, &a);
         _ntl_gmul(q, w, &x);
         _ntl_gsub(inv, x, &inv);
      }

      _ntl_gswap(&a, &n);
      _ntl_gswap(&inv, &w);
   }

   if (_ntl_gsign(w) < 0) {
      _ntl_gnegate(&w);
      _ntl_gnegate(&n);
   }

   _ntl_gcopy(n, num_out);
   _ntl_gcopy(w, den_out);

   return 1;
}


void
_ntl_gexp(
	_ntl_gbigint a,
	long e,
	_ntl_gbigint *bb
	)
{
	long k;
	long len_a;
	static _ntl_gbigint res = 0;

	if (!e)
	{
		_ntl_gone(bb);
		return;
	}

	if (e < 0)
		ghalt("negative exponent in _ntl_gexp");

	if (_ntl_giszero(a))
	{
		_ntl_gzero(bb);
		return;
	}

	len_a = _ntl_g2log(a);
	if (len_a > (NTL_MAX_LONG-(NTL_ZZ_NBITS-1))/e)
		ghalt("overflow in _ntl_gexp");

	_ntl_gsetlength(&res, (len_a*e+NTL_ZZ_NBITS-1)/NTL_ZZ_NBITS);

	_ntl_gcopy(a, &res);
	k = 1;
	while ((k << 1) <= e)
		k <<= 1;
	while (k >>= 1) {
		_ntl_gsq(res, &res);
		if (e & k)
			_ntl_gmul(a, res, &res);
	}

	_ntl_gcopy(res, bb);
}

void
_ntl_gexps(
	long a,
	long e,
	_ntl_gbigint *bb
	)
{
	long k;
	long len_a;
	static _ntl_gbigint res = 0;

	if (!e)
	{
		_ntl_gone(bb);
		return;
	}

	if (e < 0)
		ghalt("negative exponent in _ntl_zexps");

	if (!a)
	{
		_ntl_gzero(bb);
		return;
	}

	len_a = _ntl_g2logs(a);
	if (len_a > (NTL_MAX_LONG-(NTL_ZZ_NBITS-1))/e)
		ghalt("overflow in _ntl_gexps");

	_ntl_gsetlength(&res, (len_a*e+NTL_ZZ_NBITS-1)/NTL_ZZ_NBITS);

	_ntl_gintoz(a, &res);
	k = 1;
	while ((k << 1) <= e)
		k <<= 1;
	while (k >>= 1) {
		_ntl_gsq(res, &res);
		if (e & k)
			_ntl_gsmul(res, a, &res);
	}

	_ntl_gcopy(res, bb);
}


static
long OptWinSize(long n)
/* finds k that minimizes n/(k+1) + 2^{k-1} */

{
   long k;
   double v, v_new;


   v = n/2.0 + 1.0;
   k = 1;

   for (;;) {
      v_new = n/((double)(k+2)) + ((double)(1L << k));
      if (v_new >= v) break;
      v = v_new;
      k++;
   }

   return k;
}



/* DIRT: will not work with non-empty "nails" */

static
mp_limb_t neg_inv_mod_limb(mp_limb_t m0)
{
   mp_limb_t x; 
   long k;

   x = 1; 
   k = 1; 
   while (k < NTL_ZZ_NBITS) {
      x += x * (1 - x * m0);
      k <<= 1;
   }


   return - x;
}


/* Montgomery reduction:
 * This computes res = T/b^m mod N, where b = 2^{NTL_ZZ_NBITS}.
 * It is assumed that N has n limbs, and that T has at most n + m limbs.
 * Also, inv should be set to -N^{-1} mod b.
 * Finally, it is assumed that T has space allocated for n + m limbs,
 * and that res has space allocated for n limbs.  
 * Note: res should not overlap any inputs, and T is destroyed.
 * Note: res will have at most n limbs, but may not be fully reduced
 * mod N.  In general, we will have res < T/b^m + N.
 */

/* DIRT: this routine may not work with non-empty "nails" */

static
void redc(_ntl_gbigint T, _ntl_gbigint N, long m, mp_limb_t inv, 
          _ntl_gbigint res) 
{
   long n, sT, i;
   mp_limb_t *Ndata, *Tdata, *resdata, q, d, t, c;

   n = SIZE(N);
   Ndata = DATA(N);
   sT = SIZE(T);
   Tdata = DATA(T);
   resdata = DATA(res);

   for (i = sT; i < m+n; i++)
      Tdata[i] = 0;

   c = 0;
   for (i = 0; i < m; i++) {
      q = Tdata[i]*inv;