g_lip.c

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

   _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 <= 0 || !ilo || ilo < (long)hi)
                  fast = 0;
               else
               {
                  dt = lo-ilo;
                  lo = flo / dirt;
                  if (dt > 0)
                     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;
}

static
mp_limb_t 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.
 */

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;
      d = mpn_addmul_1(Tdata+i, Ndata, n, q);
      t = Tdata[i+n] + d;
      Tdata[i+n] = t + c;
      if (t < d || (c == 1 && t + c  == 0)) 
         c = 1;
      else
         c = 0;
   }

   if (c) {
      mpn_sub_n(resdata, Tdata + m, Ndata, n);
   }
   else {
      for (i = 0; i < n; i++)
         resdata[i] = Tdata[m + i];
   }

   i = n;
   STRIP(i, resdata);

   SIZE(res) = i;
   SIZE(T) = 0;
}



#define REDC_CROSS (32)

void _ntl_gpowermod(_ntl_gbigint g, _ntl_gbigint e, _ntl_gbigint F,
                    _ntl_gbigint *h)

/* h = g^e mod f using "sliding window" algorithm

   remark: the notation (h, g, e, F) is strange, because I
   copied the code from BB.c.
*/

{
   _ntl_gbigint res, gg, *v, t;
   long n, i, k, val, cnt, m;
   long use_redc, sF;
   mp_limb_t inv;

   if (_ntl_gsign(g) < 0 || _ntl_gcompare(g, F) >= 0 || 
       _ntl_gscompare(F, 1) <= 0) 
      ghalt("PowerMod: bad args");

   if (_ntl_gscompare(e, 0) == 0) {
      _ntl_gone(h);
      return;
   }

   if (_ntl_gscompare(e, 1) == 0) {
      _ntl_gcopy(g, h);
      return;
   }

   if (_ntl_gscompare(e, -1) == 0) {
      _ntl_ginvmod(g, F, h);
      return;
   }

   if (_ntl_gscompare(e, 2) == 0) {
      _ntl_gsqmod(g, F, h);
      return;
   }

   if (_ntl_gscompare(e, -2) == 0) {
      res = 0;
      _ntl_gsqmod(g, F, &res);
      _ntl_ginvmod(res, F, h);
      _ntl_gfree(&res);
      return;
   }

   n = _ntl_g2log(e);

   sF = SIZE(F);

   res = 0;
   _ntl_gsetlength(&res, sF*2);

   t = 0;
   _ntl_gsetlength(&t, sF*2);

   use_redc = (DATA(F)[0] & 1) && sF < REDC_CROSS;

   gg = 0;

   if (use_redc) {
      _ntl_glshift(g, sF*NTL_ZZ_NBITS, &res);
      _ntl_gmod(res, F, &gg);

      inv = - inv_mod_limb(DATA(F)[0]);
   }
   else
      _ntl_gcopy(g, &gg);


   if (_ntl_gscompare(g, 2) == 0) {
      /* plain square-and-multiply algorithm, optimized for g == 2 */

      _ntl_gbigint F1 = 0;

      if (use_redc) {
         long shamt;

         COUNT_BITS(shamt, DATA(F)[sF-1]);
         shamt = NTL_ZZ_NBITS - shamt;
         _ntl_glshift(F, shamt, &F1);
      }

      _ntl_gcopy(gg, &res);

      for (i = n - 2; i >= 0; i--) {
         _ntl_gsq(res, &t);
         if (use_redc) redc(t, F, sF, inv, res); else _ntl_gmod(t, F, &res);

         if (_ntl_gbit(e, i)) {
            _ntl_gadd(res, res, &res);

            if (use_redc) {
               while (SIZE(res) > sF) {
                  _ntl_gsubpos(res, F1, &res);
               }
            }
            else {
               if (_ntl_gcompare(res, F) >= 0)
                  _ntl_gsubpos(res, F, &res);
            }
         }
      }


      if (use_redc) {
         _ntl_gcopy(res, &t);
         redc(t, F, sF, inv, res);
         if (_ntl_gcompare(res, F) >= 0) {
            _ntl_gsub(res, F, &res);
         }
      }

      if (_ntl_gsign(e) < 0) _ntl_ginvmod(res, F, &res);

      _ntl_gcopy(res, h);
      _ntl_gfree(&res);
      _ntl_gfree(&gg);
      _ntl_gfree(&t);
      _ntl_gfree(&F1);
      return;
   }


   if (n < 16) { 
      /* plain square-and-multiply algorithm */

      _ntl_gcopy(gg, &res);

      for (i = n - 2; i >= 0; i--) {
         _ntl_gsq(res, &t);
         if (use_redc) redc(t, F, sF, inv, res); else _ntl_gmod(t, F, &res);

         if (_ntl_gbit(e, i)) {
            _ntl_gmul(res, gg, &t);
            if (use_redc) redc(t, F, sF, inv, res); else _ntl_gmod(t, F, &res);
         }
      }


      if (use_redc) {
         _ntl_gcopy(res, &t);
         redc(t, F, sF, inv, res);
         if (_ntl_gcompare(res, F) >= 0) {
            _ntl_gsub(res, F, &res);
         }
      }

      if (_ntl_gsign(e) < 0) _ntl_ginvmod(res, F, &res);

      _ntl_gcopy(res, h);
      _ntl_gfree(&res);
      _ntl_gfree(&gg);
      _ntl_gfree(&t);
      return;
   }

   k = OptWinSize(n);

   if (k > 5) k = 5;

   v = (_ntl_gbigint *) malloc((1L << (k-1))*sizeof(_ntl_gbigint));
   if (!v) ghalt("out of memory");
   for (i = 0; i < (1L << (k-1)); i++) {
      v[i] = 0; 
      _ntl_gsetlength(&v[i], sF);
   }

   _ntl_gcopy(gg, &v[0]);
 
   if (k > 1) {
      _ntl_gsq(gg, &t);
      if (use_redc) redc(t, F, sF, inv, res); else _ntl_gmod(t, F, &res);

      for (i = 1; i < (1L << (k-1)); i++) {
         _ntl_gmul(v[i-1], res, &t);
         if (use_redc) redc(t, F, sF, inv, v[i]); else _ntl_gmod(t, F, &v[i]);
      }
   }

   _ntl_gcopy(gg, &res);

   val = 0;
   for (i = n-2; i >= 0; i--) {
      val = (val << 1) | _ntl_gbit(e, i); 
      if (val == 0) {
         _ntl_gsq(res, &t);
         if (use_redc) redc(t, F, sF, inv, res); else _ntl_gmod(t, F, &res);
      }
      else if (val >= (1L << (k-1)) || i == 0) {
         cnt = 0;
         while ((val & 1) == 0) {
            val = val >> 1;
            cnt++;
         }

         m = val;
         while (m > 0) {
            _ntl_gsq(res, &t);
            if (use_redc) redc(t, F, sF, inv, res); else _ntl_gmod(t, F, &res);
            m = m >> 1;
         }

         _ntl_gmul(res, v[val >> 1], &t);
         if (use_redc) redc(t, F, sF, inv, res); else _ntl_gmod(t, F, &res);

         while (cnt > 0) {
            _ntl_gsq(res, &t);
            if (use_redc) redc(t, F, sF, inv, res); else _ntl_gmod(t, F, &res);
            cnt--;
         }

         val = 0;
      }
   }

   if (use_redc) {
      _ntl_gcopy(res, &t);
      redc(t, F, sF, inv, res);
      if (_ntl_gcompare(res, F) >= 0) {
         _ntl_gsub(res, F, &res);
      }
   }

   if (_ntl_gsign(e) < 0) _ntl_ginvmod(res, F, &res);

   _ntl_gcopy(res, h);

   _ntl_gfree(&res);
   _ntl_gfree(&gg);
   _ntl_gfree(&t);
   for (i = 0; i < (1L << (k-1)); i++)
      _ntl_gfree(&v[i]);
   free(v);
}

long _ntl_gsize(_ntl_gbigint rep)
{
   if (!rep)
      return 0;
   else if (SIZE(rep) < 0)
      return -SIZE(rep);
   else
      return SIZE(rep);
}

long _ntl_gisone(_ntl_gbigint rep)
{
   return rep != 0 && SIZE(rep) == 1 && DATA(rep)[0] == 1;
}

long _ntl_gsptest(_ntl_gbigint rep)
{
   return !rep || SIZE(rep) == 0 ||
          ((SIZE(rep) == 1 || SIZE(rep) == -1) && 
           DATA(rep)[0] < ((mp_limb_t) NTL_SP_BOUND));
}

long _ntl_gwsptest(_ntl_gbigint rep)
{
   return !rep || SIZE(rep) == 0 ||
          ((SIZE(rep) == 1 || SIZE(rep) == -1) && 
           DATA(rep)[0] < ((mp_limb_t) NTL_WSP_BOUND));
}



long _ntl_gcrtinrange(_ntl_gbigint g, _ntl_gbigint a)
{
   long sa, sg, i; 
   mp_limb_t carry, u, v;
   mp_limb_t *adata, *gdata;

   if (!a || SIZE(a) <= 0) return 0;

   sa = SIZE(a);

   if (!g) return 1;

   sg = SIZE(g);

   if (sg == 0) return 1;

   if (sg < 0) sg = -sg;

   if (sa-sg > 1) return 1;

   if (sa-sg < 0) return 0;

   adata = DATA(a);
   gdata = DATA(g);

   carry=0;

   if (sa-sg ==