g_lip.c

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

      sc = sa;
      if (MustAlloc(c, sc)) {
         _ntl_gsetlength(&c, sc);
         if (a == *cc) a = c; 
         if (b == *cc) b = c;
         *cc = c;
      }

      adata = DATA(a);
      bdata = DATA(b);
      cdata = DATA(c);

      if (sa > sb) 
         cmp = 1;
      else
         cmp = mpn_cmp(adata, bdata, sa);

      if (cmp == 0) {
         SIZE(c) = 0;
      }
      else {
         if (cmp < 0) cmp = 0;
         if (cmp > 0) cmp = 1;
         /* abs(a) != abs(b) && (abs(a) > abs(b) <=> cmp) */

         if (cmp)
            mpn_sub(cdata, adata, sa, bdata, sb);
         else
            mpn_sub(cdata, bdata, sb, adata, sa); /* sa == sb */

         STRIP(sc, cdata);
         if ((aneg == cmp) ^ rev) sc = -sc;
         SIZE(c) = sc;
      }
   }
}

void
_ntl_gsubpos(_ntl_gbigint a, _ntl_gbigint b, _ntl_gbigint *cc)
{
   long sa, sb, sc;
   mp_limb_t *adata, *bdata, *cdata;
   _ntl_gbigint c;

   if (ZEROP(a)) {
      _ntl_gzero(cc);
      return;
   }

   if (ZEROP(b)) {
      _ntl_gcopy(a, cc);
      return;
   }

   sa = SIZE(a);
   sb = SIZE(b);

   c = *cc;

   sc = sa;
   if (MustAlloc(c, sc)) {
      _ntl_gsetlength(&c, sc);
      if (a == *cc) a = c; 
      if (b == *cc) b = c;
      *cc = c;
   }

   adata = DATA(a);
   bdata = DATA(b);
   cdata = DATA(c);

   mpn_sub(cdata, adata, sa, bdata, sb);

   STRIP(sc, cdata);
   SIZE(c) = sc;
}

void _ntl_gmul(_ntl_gbigint a, _ntl_gbigint b, _ntl_gbigint *cc)
{
   static _ntl_gbigint mem = 0;

   long sa, aneg, sb, bneg, alias, sc;
   mp_limb_t *adata, *bdata, *cdata, msl;
   _ntl_gbigint c;

   if (ZEROP(a) || ZEROP(b)) {
      _ntl_gzero(cc);
      return;
   }

   GET_SIZE_NEG(sa, aneg, a);
   GET_SIZE_NEG(sb, bneg, b);

   if (a == *cc || b == *cc) {
      c = mem;
      alias = 1;
   }
   else {
      c = *cc;
      alias = 0;
   }

   sc = sa + sb;
   if (MustAlloc(c, sc))
      _ntl_gsetlength(&c, sc);

   if (alias)
      mem = c;
   else
      *cc = c;

   adata = DATA(a);
   bdata = DATA(b);
   cdata = DATA(c);

   if (sa >= sb)
      msl = mpn_mul(cdata, adata, sa, bdata, sb);
   else
      msl = mpn_mul(cdata, bdata, sb, adata, sa);

   if (!msl) sc--;
   if (aneg != bneg) sc = -sc;
   SIZE(c) = sc;

   if (alias) _ntl_gcopy(mem, cc);
}

void _ntl_gsq(_ntl_gbigint a, _ntl_gbigint *cc)
{
   _ntl_gmul(a, a, cc);
   /* this is good enough...eventually, mpn_sqr_n will be called */
}


void
_ntl_gsmul(_ntl_gbigint a, long d, _ntl_gbigint *bb)
{
   long sa, sb;
   long anegative, bnegative;
   _ntl_gbigint b;
   mp_limb_t *adata, *bdata;
   mp_limb_t dd, carry;

   if (ZEROP(a) || !d) {
      _ntl_gzero(bb);
      return;
   }

   GET_SIZE_NEG(sa, anegative, a);

   if (d < 0) {
      dd = - ((mp_limb_t) d); /* careful ! */
      bnegative = 1-anegative;
   }
   else {
      dd = (mp_limb_t) d;
      bnegative = anegative;
   }

   sb = sa + 1;

   b = *bb;
   if (MustAlloc(b, sb)) {
      _ntl_gsetlength(&b, sb);
      if (a == *bb) a = b;
      *bb = b;
   }

   adata = DATA(a);
   bdata = DATA(b);

   if (dd == 2) 
      carry = mpn_lshift(bdata, adata, sa, 1);
   else
      carry = mpn_mul_1(bdata, adata, sa, dd);

   if (carry) 
      bdata[sa] = carry;
   else
      sb--;

   if (bnegative) sb = -sb;
   SIZE(b) = sb;
}

long _ntl_gsdiv(_ntl_gbigint a, long d, _ntl_gbigint *bb)
{
   long sa, aneg, sb, dneg;
   _ntl_gbigint b;
   mp_limb_t dd, *adata, *bdata;
   long r;

   if (!d) {
      ghalt("division by zero in _ntl_gsdiv");
   }

   if (ZEROP(a)) {
      _ntl_gzero(bb);
      return (0);
   }

   GET_SIZE_NEG(sa, aneg, a);

   if (d < 0) {
      dd = - ((mp_limb_t) d); /* careful ! */
      dneg = 1;
   }
   else {
      dd = (mp_limb_t) d;
      dneg = 0;
   }

   sb = sa;
   b = *bb;
   if (MustAlloc(b, sb)) {
      _ntl_gsetlength(&b, sb);
      if (a == *bb) a = b;
      *bb = b;
   }

   adata = DATA(a);
   bdata = DATA(b);

   if (dd == 2) 
      r = mpn_rshift(bdata, adata, sa, 1) >> (NTL_ZZ_NBITS - 1);
   else
      r = mpn_divmod_1(bdata, adata, sa, dd);

   if (bdata[sb-1] == 0)
      sb--;

   SIZE(b) = sb;

   if (aneg || dneg) {
      if (aneg != dneg) {
         if (!r) {
            SIZE(b) = -SIZE(b);
         }
         else {
            _ntl_gsadd(b, 1, &b);
            SIZE(b) = -SIZE(b);
            if (dneg)
               r = r + d;
            else
               r = d - r;
            *bb = b;
         }
      }
      else
         r = -r;
   }

   return r;
}

long _ntl_gsmod(_ntl_gbigint a, long d)
{
   long sa, aneg, dneg;
   mp_limb_t dd, *adata;
   long r;

   if (!d) {
      ghalt("division by zero in _ntl_gsmod");
   }

   if (ZEROP(a)) {
      return (0);
   }

   GET_SIZE_NEG(sa, aneg, a);

   if (d < 0) {
      dd = - ((mp_limb_t) d); /* careful ! */
      dneg = 1;
   }
   else {
      dd = (mp_limb_t) d;
      dneg = 0;
   }

   adata = DATA(a);

   if (dd == 2) 
      r = adata[0] & 1;
   else
      r = mpn_mod_1(adata, sa, dd);

   if (aneg || dneg) {
      if (aneg != dneg) {
         if (r) {
            if (dneg)
               r = r + d;
            else
               r = d - r;
         }
      }
      else
         r = -r;
   }

   return r;
}

void _ntl_gdiv(_ntl_gbigint a, _ntl_gbigint d, 
               _ntl_gbigint *bb, _ntl_gbigint *rr)
{
   static _ntl_gbigint b = 0, rmem = 0;

   _ntl_gbigint r;

   long sa, aneg, sb, sd, dneg, sr, in_place;
   mp_limb_t *adata, *ddata, *bdata, *rdata;

   if (ZEROP(d)) {
      ghalt("division by zero in _ntl_gdiv");
   }

   if (ZEROP(a)) {
      if (bb) _ntl_gzero(bb);
      if (rr) _ntl_gzero(rr);
      return;
   }

   GET_SIZE_NEG(sa, aneg, a);
   GET_SIZE_NEG(sd, dneg, d);

   if (!aneg && !dneg && rr && *rr != a && *rr != d) {
      in_place = 1;
      r = *rr;
   }
   else {
      in_place = 0;
      r = rmem;
   }

   if (sa < sd) {
      _ntl_gzero(&b);
      _ntl_gcopy(a, &r);
      if (aneg) SIZE(r) = -SIZE(r);
      goto done;
   }

   sb = sa-sd+1;
   if (MustAlloc(b, sb))
      _ntl_gsetlength(&b, sb);

   sr = sd;
   if (MustAlloc(r, sr))
      _ntl_gsetlength(&r, sr);

   adata = DATA(a);
   ddata = DATA(d);
   bdata = DATA(b);
   rdata = DATA(r);

   mpn_tdiv_qr(bdata, rdata, 0, adata, sa, ddata, sd);

   if (bdata[sb-1] == 0)
      sb--;
   SIZE(b) = sb;

   STRIP(sr, rdata);
   SIZE(r) = sr;

done:

   if (aneg || dneg) {
      if (aneg != dneg) {
         if (ZEROP(r)) {
            SIZE(b) = -SIZE(b);
         }
         else {
            if (bb) {
               _ntl_gsadd(b, 1, &b);
               SIZE(b) = -SIZE(b);
            }
            if (rr) {
               if (dneg)
                  _ntl_gadd(r, d, &r);
               else
                  _ntl_gsub(d, r, &r);
            }
         }
      }
      else
         SIZE(r) = -SIZE(r);
   }

   if (bb) _ntl_gcopy(b, bb);

   if (in_place)
      *rr = r;
   else {
      if (rr) _ntl_gcopy(r, rr);
      rmem = r;
   }
}


/* a simplified mod operation:  assumes a >= 0, d > 0 are non-negative,
 * that space for the result has already been allocated,
 * and that inputs do not alias output. */

void gmod_simple(_ntl_gbigint a, _ntl_gbigint d, _ntl_gbigint *rr)
{
   static _ntl_gbigint b = 0;

   long sa, sb, sd, sr;
   mp_limb_t *adata, *ddata, *bdata, *rdata;
   _ntl_gbigint r;

   if (ZEROP(a)) {
      _ntl_gzero(rr);
      return;
   }

   sa = SIZE(a);
   sd = SIZE(d);

   if (sa < sd) {
      _ntl_gcopy(a, rr);
      return;
   }

   sb = sa-sd+1;
   if (MustAlloc(b, sb))
      _ntl_gsetlength(&b, sb);

   sr = sd;
   r = *rr;

   adata = DATA(a);
   ddata = DATA(d);
   bdata = DATA(b);
   rdata = DATA(r);

   mpn_tdiv_qr(bdata, rdata, 0, adata, sa, ddata, sd);

   STRIP(sr, rdata);
   SIZE(r) = sr;
}

void _ntl_gmod(_ntl_gbigint a, _ntl_gbigint d, _ntl_gbigint *rr)
{
   _ntl_gdiv(a, d, 0, rr);
}

void _ntl_gquickmod(_ntl_gbigint *rr, _ntl_gbigint d)
{
   _ntl_gdiv(*rr, d, 0, rr);
}

void _ntl_gsqrt(_ntl_gbigint n, _ntl_gbigint *rr)
{
   static _ntl_gbigint r = 0;

   long sn, sr;
   mp_limb_t *ndata, *rdata;

   if (ZEROP(n)) {
      _ntl_gzero(rr);
      return;
   }

   sn = SIZE(n);
   if (sn < 0) ghalt("negative argument to _ntl_sqrt");

   sr = (sn+1)/2;
   _ntl_gsetlength(&r, sr);

   ndata = DATA(n);
   rdata = DATA(r);

   mpn_sqrtrem(rdata, 0, ndata, sn);

   STRIP(sr, rdata);
   SIZE(r) = sr;

   _ntl_gcopy(r, rr);
}

long _ntl_gsqrts(long n)
{
   mp_limb_t ndata, rdata;

   if (n == 0) {
      return 0;
   }

   if (n < 0) ghalt("negative argument to _ntl_sqrts");

   ndata = n;

   mpn_sqrtrem(&rdata, 0, &ndata, 1);

   return rdata;
}


void _ntl_ggcd(_ntl_gbigint m1, _ntl_gbigint m2, _ntl_gbigint *r)
{
   static _ntl_gbigint s1 = 0, s2 = 0, res = 0;

   long k1, k2, k_min, l1, l2, ss1, ss2, sres;

   _ntl_gcopy(m1, &s1);
   _ntl_gabs(&s1);

   _ntl_gcopy(m2, &s2);
   _ntl_gabs(&s2);

   if (ZEROP(s1)) {
      _ntl_gcopy(s2, r);
      return;
   }

   if (ZEROP(s2)) {
      _ntl_gcopy(s1, r);
      return;
   }

   k1 = _ntl_gmakeodd(&s1);
   k2 = _ntl_gmakeodd(&s2);

   if (k1 <= k2)
      k_min = k1;
   else
      k_min = k2;

   l1 = _ntl_g2log(s1);
   l2 = _ntl_g2log(s2);

   ss1 = SIZE(s1);
   ss2 = SIZE(s2);

   if (ss1 >= ss2)
      sres = ss1;
   else
      sres = ss2;

   /* set to max: gmp documentation is unclear on this point */

   _ntl_gsetlength(&res, sres);
   
   if (l1 >= l2)
      SIZE(res) = mpn_gcd(DATA(res), DATA(s1), ss1, DATA(s2), ss2);
   else
      SIZE(res) = mpn_gcd(DATA(res), DATA(s2), ss2, DATA(s1), ss1);

   _ntl_glshift(res, k_min, &res);

   _ntl_gcopy(res, r);
}

static long 
gxxeucl(
   _ntl_gbigint ain,
   _ntl_gbigint nin,
   _ntl_gbigint *invv,
   _ntl_gbigint *uu
   )
{
   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;
   _ntl_gbigint inv = *invv;
   _ntl_gbigint u = *uu;
   long diff;
   long ilo;
   long sa;
   long sn;
   long temp;
   long e;
   long fast;
   long parity;
   long gotthem;
   mp_limb_t *p;
   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;

   _ntl_gsetlength(&a, (e = (SIZE(ain) > SIZE(nin) ? SIZE(ain) : SIZE(nin))));
   _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);
   *invv = inv;
   _ntl_gsetlength(&u, e);
   *uu = u;

   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 (SIZE(n) > 0)
   {
      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)
      {