c_lip_impl.h

数值算法库for Windows

   if ((sa = *a) < 0)
      sa = -sa;

   d = *(a += sa);
   while (--sa) {
      d <<= NTL_NBITS;
      d += *(--a);
   }

   if ((*(--a)) < 0)
      return (-d);
   return (d);
}

unsigned long _ntl_ztouint(_ntl_verylong a)
{
   unsigned long d;
   long sa;

   if (!a)
      return (0);

   if ((sa = *a) < 0)
      sa = -sa;

   d = (unsigned long) (*(a += sa));
   while (--sa) {
      d <<= NTL_NBITS;
      d += (unsigned long) (*(--a));
   }

   if ((*(--a)) < 0)
      return (-d);
   return (d);
}


long _ntl_zcompare(_ntl_verylong a, _ntl_verylong b)
{
   long sa;
   long sb;

   if (!a) {
      if (!b)
         return (0);
      if (b[0] < 0)
         return (1);
      if (b[0] > 1)
         return (-1);
      if (b[1])
         return (-1);
      return (0);
   }
   if (!b) {
      if (a[0] < 0)
         return (-1);
      if (a[0] > 1)
         return (1);
      if (a[1])
         return (1);
      return (0);
   }

   if ((sa = *a) > (sb = *b))
      return (1);
   if (sa < sb)
      return (-1);
   if (sa < 0)
      sa = (-sa);
   a += sa;
   b += sa;
   for (; sa; sa--) {
      if (*a > *b) {
         if (sb < 0)
            return (-1);
         return (1);
      }
      if (*a-- < *b--) {
         if (sb < 0)
            return (1);
         return (-1);
      }
   }
   return (0);
}

void _ntl_znegate(_ntl_verylong *aa)
{
   _ntl_verylong a = *aa;

   if (!a)
      return;
   if (a[1] || a[0] != 1)
      a[0] = (-a[0]);
}

void _ntl_zsadd(_ntl_verylong a, long d, _ntl_verylong *b)
{
   static _ntl_verylong x = 0;

   _ntl_zintoz(d, &x);
   _ntl_zadd(a, x, b);
}


void
_ntl_zadd(_ntl_verylong a, _ntl_verylong b, _ntl_verylong *cc)
{
   long sa;
   long sb;
   long anegative;

   if (!a) {
      if (b)
         _ntl_zcopy(b, cc);
      else
         _ntl_zzero(cc);
      return;
   }

   if (!b) {
      _ntl_zcopy(a, cc);
      return;
   }

   if ((anegative = ((sa = a[0]) < 0)) == ((sb = b[0]) < 0)) {
      /* signs a and b are the same */
      _ntl_verylong pc, c;
      long carry;
      long i;
      long maxab;

      if (anegative) {
         sa = -sa;
         sb = -sb;
      }

      if (sa < sb) {
         i = sa;
         maxab = sb;
      }
      else {
         i = sb;
         maxab = sa;
      }

      c = *cc;
      if (MustAlloc(c, maxab+1)) {
         _ntl_verylong c1 = c;
         _ntl_zsetlength(&c1, maxab + 1);
         if (a == c) a = c1;
         if (b == c) b = c1;
         *cc = c = c1;
      }

      pc = c;
      carry = 0;

      do {
         long t = (*(++a)) + (*(++b)) + carry;
         carry = t >> NTL_NBITS;
         *(++pc) = t & NTL_RADIXM;
         i--;
      } while (i);

      i = sa-sb;
      if (!i) goto i_exit;

      if (i < 0) {
         i = -i;
         a = b;
      }

      if (!carry) goto carry_exit;

      for (;;) {
         long t = (*(++a)) + 1;
         carry = t >> NTL_NBITS;
         *(++pc) = t & NTL_RADIXM;
         i--;
         if (!i) goto i_exit;
         if (!carry) goto carry_exit;
      }

      i_exit:
      if (carry) {
         *(++pc) = 1;
         maxab++;
      }
      *c = anegative ? -maxab : maxab;
      return;

      carry_exit:
      if (pc != a) {
         do {
            *(++pc) = *(++a);
            i--;
         } while (i);
      }
      *c = anegative ? -maxab : maxab;
   }
   else {
      /* signs a and b are different...use _ntl_zsub */

      _ntl_verylong c = *cc;

      if (anegative) {
         a[0] = -sa;
         _ntl_zsub(b, a, cc);
         if (a != c) a[0] = sa;
      }
      else {
         b[0] = -sb;
         _ntl_zsub(a, b, cc);
         if (b != c) b[0] = sb;
      }
   }
}

void
_ntl_zsub(_ntl_verylong a, _ntl_verylong b, _ntl_verylong *cc)
{
   long sa;
   long sb;
   long anegative;

   if (!b) {
      if (a)
         _ntl_zcopy(a, cc);
      else
         _ntl_zzero(cc);
      return;
   }

   if (!a) {
      _ntl_zcopy(b, cc);
      _ntl_znegate(cc);
      return;
   }

   if ((anegative = ((sa = a[0]) < 0)) == ((sb = b[0]) < 0)) {
      /* signs agree */

      long i, carry, *pc, *c;

      if (anegative) {
         sa = -sa;
         sb = -sb;
      }

      carry = sa - sb;
      if (!carry) {
         long *aa = a + sa;
         long *bb = b + sa;
         
         i = sa;
         while (i && !(carry = (*aa - *bb))) {
            aa--; bb--; i--;
         }
      }

      if (!carry) {
         _ntl_zzero(cc);
         return;
      }

      if (carry < 0) {
         { long t = sa; sa = sb; sb = t; }
         { long *t = a; a = b; b = t; }
         anegative = !anegative;
      }

      c = *cc;
      if (MustAlloc(c, sa)) {
         _ntl_verylong c1 = c;
         _ntl_zsetlength(&c1, sa);
         if (b == c) b = c1;
         *cc = c = c1;
      }

      i = sb;
      carry = 0;
      pc = c;

      do {
#if (!NTL_ARITH_RIGHT_SHIFT)
         long t = (*(++a)) - (*(++b)) - carry;
         carry = (t < 0);
#else
         long t = (*(++a)) - (*(++b)) + carry;
         carry = t >> NTL_NBITS;
#endif
         *(++pc) = t & NTL_RADIXM;
         i--;
      } while (i);

      i = sa-sb;
      while (carry) {
         long t = (*(++a)) - 1;
#if (!NTL_ARITH_RIGHT_SHIFT)
         carry = (t < 0);
#else
         carry = t >> NTL_NBITS;
#endif
         *(++pc) = t & NTL_RADIXM;
         i--;
      }

      if (i) {
         if (pc != a) {
            do {
               *(++pc) = *(++a);
               i--;
            } while (i);
         }
      }
      else {
         while (sa > 1 && *pc == 0) { sa--; pc--; } 
      }

      if (anegative) sa = -sa;
      *c = sa;
   }
   else {
      /* signs of a and b are different...use _ntl_zadd */

      _ntl_verylong c = *cc;

      if (anegative) {
         a[0] = -sa;
         _ntl_zadd(a, b, cc);
         if (a != c) a[0] = sa;
         c = *cc;
         c[0] = -c[0];
      }
      else {
         b[0] = -sb;
         _ntl_zadd(a, b, cc);
         if (b != c) b[0] = sb;
      }
   }
}


void
_ntl_zsmul(_ntl_verylong a, long d, _ntl_verylong *bb)
{
   long sa;
   long anegative, bnegative;
   _ntl_verylong b;


   if (d == 2) {
      _ntl_z2mul(a, bb);
      return;
   }


   if ((d >= NTL_RADIX) || (d <= -NTL_RADIX)) {
      static _ntl_verylong x = 0;
      _ntl_zintoz(d,&x);
      _ntl_zmul(a, x, bb);
      return;
   }

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

   if (!d) {
      _ntl_zzero(bb);
      return;
   }


   anegative = 0;
   bnegative = 0;

   if ((sa = a[0]) < 0) {
      anegative = 1;
      a[0] = sa = (-sa);
      if (d < 0)
         d = (-d);
      else
         bnegative = 1;
   }
   else if (bnegative = (d < 0))
      d = (-d);

   b = *bb;

   if (MustAlloc(b, sa + 1)) {
      _ntl_zsetlength(&b, sa + 1);
      if (a == *bb) a = b;
      *bb = b;
   }

   zsxmul(d, b+1, a);

   sa++;
   while ((sa > 1) && (!(b[sa])))
      sa--;
   b[0] = sa;

   if (bnegative && (b[1] || b[0] != 1))
      b[0] = (-b[0]);

   if (anegative && a != b)
      a[0] = -a[0];
}

void _ntl_zsubpos(_ntl_verylong a, _ntl_verylong b, _ntl_verylong *cc)
{
   long sa;
   long sb;

   long *c, *pc;
   long i, carry;

   if (!b) {
      if (a)
         _ntl_zcopy(a, cc);
      else
         _ntl_zzero(cc);
      return;
   }

   if (!a) {
      _ntl_zzero(cc);
      return;
   }

   sa = a[0];
   sb = b[0];

   c = *cc;
   if (MustAlloc(c, sa)) {
      _ntl_verylong c1 = c;
      _ntl_zsetlength(&c1, sa);
      if (b == c) b = c1;
      *cc = c = c1;
   }

   i = sb;
   carry = 0;
   pc = c;

   while (i) {
#if (!NTL_ARITH_RIGHT_SHIFT)
      long t = (*(++a)) - (*(++b)) - carry;
      carry = (t < 0);
#else
      long t = (*(++a)) - (*(++b)) + carry;
      carry = t >> NTL_NBITS;
#endif
      *(++pc) = t & NTL_RADIXM;
      i--;
   }

   i = sa-sb;
   while (carry) {
      long t = (*(++a)) - 1;
#if (!NTL_ARITH_RIGHT_SHIFT)
      carry = (t < 0);
#else
      carry = t >> NTL_NBITS;
#endif
      *(++pc) = t & NTL_RADIXM;
      i--;
   }

   if (i) {
      if (pc != a) {
         do {
            *(++pc) = *(++a);
            i--;
         } while (i);
      }
   }
   else {
      while (sa > 1 && *pc == 0) { sa--; pc--; } 
   }

   *c = sa;
}



static long *kmem = 0;     /* globals for Karatsuba */
static long max_kmem = 0;


/* These cross-over points were estimated using
   a Sparc-10, a Sparc-20, and a Pentium-90. */

#if (defined(NTL_SINGLE_MUL))
#define KARX (18)
#else
#define KARX (16) 
#endif

/* Auxilliary routines for Karatsuba */


static
void kar_fold(long *T, long *b, long hsa)
{
   long sb, *p2, *p3, i, carry;

   sb = *b;
   p2 = b + hsa;
   p3 = T;
   carry = 0;

   for (i = sb-hsa; i>0; i--) {
      long t = (*(++b)) + (*(++p2)) + carry;
      carry = t >> NTL_NBITS;
      *(++p3) = t & NTL_RADIXM;
   }

   for (i = (hsa << 1) - sb; i>0; i--) {
      long t = (*(++b)) + carry;
      carry = t >> NTL_NBITS;
      *(++p3) = t & NTL_RADIXM;
   }

   if (carry) {
      *(++p3) = carry;
      *T = hsa + 1;
   }
   else
      *T = hsa;
}

static
void kar_sub(long *T, long *c)
{
   long i, carry;

   i = *c;
   carry = 0;

   while (i>0) {
#if (!NTL_ARITH_RIGHT_SHIFT)
      long t = (*(++T)) - (*(++c)) - carry;
      carry = (t < 0);
#else
      long t = (*(++T)) - (*(++c)) + carry;
      carry = t >> NTL_NBITS;
#endif
      *T = t & NTL_RADIXM;
      i--;
   }

   while (carry) {
      long t = (*(++T)) - 1;
#if (!NTL_ARITH_RIGHT_SHIFT)
      carry = (t < 0);
#else
      carry = t >> NTL_NBITS;
#endif
      *T = t & NTL_RADIXM;
   }
}


static
void kar_add(long *c, long *T, long hsa)
{
   long i, carry;

   c += hsa;
   i = *T;
   while (T[i] == 0 && i > 0) i--;
   carry = 0;

   while (i>0) {
      long t = (*(++c)) + (*(++T)) + carry;
      carry = t >> NTL_NBITS;
      *c = t & NTL_RADIXM;
      i--;
   }

   while (carry) {
      long t = (*(++c)) + 1;
      carry = t >> NTL_NBITS;
      *c = t & NTL_RADIXM;
   }
}

static
void kar_fix(long *c, long *T, long hsa)
{
   long i, carry, s;

   s = *T;

   i = hsa;
   while (i>0) {
      *(++c) = *(++T);
      i--;
   }


   i = s - hsa;
   carry = 0;

   while (i > 0) {
      long t = (*(++c)) + (*(++T)) + carry;
      carry = t >> NTL_NBITS;
      *c = t & NTL_RADIXM;
      i--;
   }

   while (carry) {
      long t = (*(++c)) + 1;
      carry = t >> NTL_NBITS;
      *c = t & NTL_RADIXM;
   }
}
      
   

static
void kar_mul(long *c, long *a, long *b, long *stk)
{
   long sa, sb, sc;