g_lip_impl.h

数值算法库for Windows

}

static long
gweights_mp_limb(
	mp_limb_t a
	)
{
	long res = 0;
   
	while (a) {
		if (a & 1) res ++;
		a >>= 1;
	}
	return (res);
}

long
_ntl_gweight(
        _ntl_gbigint a
        )
{
	long i;
	long sa;
	mp_limb_t *adata;
	long res;

	if (!a) return (0);

	sa = SIZE(a);
	if (sa < 0) sa = -sa;
	adata = DATA(a);

	res = 0;
	for (i = 0; i < sa; i++)
		res += gweights_mp_limb(adata[i]);

	return (res);
}

long 
_ntl_g2logs(
	long aa
	)
{
	long i = 0;
	unsigned long a;

	if (aa < 0)
		a = - ((unsigned long) aa);
	else
		a = aa;

	while (a>=256)
		i += 8, a >>= 8;
	if (a >=16)
		i += 4, a >>= 4;
	if (a >= 4)
		i += 2, a >>= 2;
	if (a >= 2)
		i += 2;
	else if (a >= 1)
		i++;
	return (i);
}

long _ntl_g2log(_ntl_gbigint a)
{
   long la;
   long t;

   if (!a) return 0;
   la = SIZE(a);
   if (la == 0) return 0;
   if (la < 0) la = -la;
   COUNT_BITS(t, DATA(a)[la-1]);
   return NTL_ZZ_NBITS*(la - 1) + t;
}


long _ntl_gmakeodd(_ntl_gbigint *nn)
{
   _ntl_gbigint n = *nn;
   long shift;
   mp_limb_t *ndata;
   mp_limb_t i;

   if (ZEROP(n))
      return (0);

   shift = 0;
   ndata = DATA(n);

   while (ndata[shift] == 0)
      shift++;

   i = ndata[shift];

   shift = NTL_ZZ_NBITS * shift;

   while ((i & 1) == 0) {
      shift++;
      i >>= 1;
   }
   _ntl_grshift(n, shift, &n);
   return shift;
}

long _ntl_gnumtwos(_ntl_gbigint n)
{
   long shift;
   mp_limb_t *ndata;
   mp_limb_t i;

   if (ZEROP(n))
      return (0);

   shift = 0;
   ndata = DATA(n);

   while (ndata[shift] == 0)
      shift++;

   i = ndata[shift];

   shift = NTL_ZZ_NBITS * shift;

   while ((i & 1) == 0) {
      shift++;
      i >>= 1;
   }

   return shift;
}

void _ntl_gand(_ntl_gbigint a, _ntl_gbigint b, _ntl_gbigint *cc)
{
   _ntl_gbigint c = *cc;
   long sa;
   long sb;
   long sm;
   long i;
   mp_limb_t *adata, *bdata, *cdata;

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

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

   sb = SIZE(b);
   if (sb < 0) sb = -sb;

   sm = (sa > sb ? sb : sa);

   _ntl_gsetlength(&c, sm);
   if (a == *cc) a = c;
   if (b == *cc) b = c;
   *cc = c;

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

   for (i = 0; i < sm; i++)
      cdata[i] = adata[i] & bdata[i];

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

void _ntl_gxor(_ntl_gbigint a, _ntl_gbigint b, _ntl_gbigint *cc)
{
   _ntl_gbigint c = *cc;
   long sa;
   long sb;
   long sm;
   long la;
   long i;
   mp_limb_t *adata, *bdata, *cdata;

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

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

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

   sb = SIZE(b);
   if (sb < 0) sb = -sb;

   if (sa > sb) {
      la = sa;
      sm = sb;
   } 
   else {
      la = sb;
      sm = sa;
   }

   _ntl_gsetlength(&c, la);
   if (a == *cc) a = c;
   if (b == *cc) b = c;
   *cc = c;

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

   for (i = 0; i < sm; i ++)
      cdata[i] = adata[i] ^ bdata[i];

   if (sa > sb)
      for (;i < la; i++) cdata[i] = adata[i];
   else
      for (;i < la; i++) cdata[i] = bdata[i];

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

void _ntl_gor(_ntl_gbigint a, _ntl_gbigint b, _ntl_gbigint *cc)
{
   _ntl_gbigint c = *cc;
   long sa;
   long sb;
   long sm;
   long la;
   long i;
   mp_limb_t *adata, *bdata, *cdata;

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

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

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

   sb = SIZE(b);
   if (sb < 0) sb = -sb;

   if (sa > sb) {
      la = sa;
      sm = sb;
   } 
   else {
      la = sb;
      sm = sa;
   }

   _ntl_gsetlength(&c, la);
   if (a == *cc) a = c;
   if (b == *cc) b = c;
   *cc = c;

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

   for (i = 0; i < sm; i ++)
      cdata[i] = adata[i] | bdata[i];

   if (sa > sb)
      for (;i < la; i++) cdata[i] = adata[i];
   else
      for (;i < la; i++) cdata[i] = bdata[i];

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


void _ntl_gnegate(_ntl_gbigint *aa)
{
   _ntl_gbigint a = *aa;
   if (a) SIZE(a) = -SIZE(a);
}

void _ntl_gintoz(long d, _ntl_gbigint *aa)
{
   _ntl_gbigint a = *aa;

   if (d == 0) {
      if (a) SIZE(a) = 0;
   }
   else if (d > 0) {
      if (!a) {
         _ntl_gsetlength(&a, 1);
         *aa = a;
      }
   
      SIZE(a) = 1;
      DATA(a)[0] = d;
   }
   else {
      if (!a) {
         _ntl_gsetlength(&a, 1);
         *aa = a;
      }
   
      SIZE(a) = -1;
      DATA(a)[0] = -((mp_limb_t) d); /* careful! */
   }
}

void _ntl_guintoz(unsigned long d, _ntl_gbigint *aa)
{
   _ntl_gbigint a = *aa;

   if (d == 0) {
      if (a) SIZE(a) = 0;
   }
   else {
      if (!a) {
         _ntl_gsetlength(&a, 1);
         *aa = a;
      }
   
      SIZE(a) = 1;
      DATA(a)[0] = d;
   }
}

long _ntl_gtoint(_ntl_gbigint a)
{
   if (ZEROP(a)) 
      return 0;

   if (SIZE(a) > 0) 
      return DATA(a)[0];

   return -DATA(a)[0];
}

unsigned long _ntl_gtouint(_ntl_gbigint a)
{
   if (ZEROP(a)) 
      return 0;

   if (SIZE(a) > 0) 
      return DATA(a)[0];

   return -DATA(a)[0];
}

long _ntl_gcompare(_ntl_gbigint a, _ntl_gbigint b)
{
   long sa, sb, cmp;
   mp_limb_t *adata, *bdata;

   if (!a) 
      sa = 0;
   else
      sa = SIZE(a);

   if (!b) 
      sb = 0;
   else
      sb = SIZE(b);

   if (sa != sb) {
      if (sa > sb)
         return 1;
      else
         return -1;
   }

   if (sa == 0)
      return 0;

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

   if (sa > 0) {
      cmp = mpn_cmp(adata, bdata, sa);

      if (cmp > 0)
         return 1;
      else if (cmp < 0) 
         return -1;
      else
         return 0;
   }
   else {
      cmp = mpn_cmp(adata, bdata, -sa);

      if (cmp > 0)
         return -1;
      else if (cmp < 0) 
         return 1;
      else
         return 0;
   }
}


long _ntl_gsign(_ntl_gbigint a)
{
   long sa;

   if (!a) return 0;

   sa = SIZE(a);
   if (sa > 0) return 1;
   if (sa == 0) return 0;
   return -1;
}

void _ntl_gabs(_ntl_gbigint *pa)
{
   _ntl_gbigint a = *pa;

   if (!a) return;
   if (SIZE(a) < 0) SIZE(a) = -SIZE(a);
}

long _ntl_gscompare(_ntl_gbigint a, long b)
{
   if (b == 0) {
      long sa;
      if (!a) return 0;
      sa = SIZE(a);
      if (sa > 0) return 1;
      if (sa == 0) return 0;
      return -1;
   }
   else {
      static _ntl_gbigint B = 0;
      _ntl_gintoz(b, &B);
      return _ntl_gcompare(a, B);
   }
}


void _ntl_glshift(_ntl_gbigint n, long k, _ntl_gbigint *rres)
{
   _ntl_gbigint res;
   mp_limb_t *ndata, *resdata, *resdata1;
   long limb_cnt, i, sn, nneg, sres;

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

   if (!k) {
      if (n != *rres)
         _ntl_gcopy(n, rres);
      return;
   }

   if (k < 0) {
      if (k < -NTL_MAX_LONG) ghalt("overflow in _ntl_glshift");
      _ntl_grshift(n, -k, rres);
      return;
   }

   GET_SIZE_NEG(sn, nneg, n);

   limb_cnt = k/NTL_ZZ_NBITS;
   sres = sn + limb_cnt + 1; 

   res = *rres;
   if (MustAlloc(res, sres)) {
      _ntl_gsetlength(&res, sres);
      if (n == *rres) n = res;
      *rres = res;
   }

   ndata = DATA(n);
   resdata = DATA(res);
   resdata1 = resdata + limb_cnt;
   k %= NTL_ZZ_NBITS;
   sres--;

   if (k != 0) {
      mp_limb_t t = mpn_lshift(resdata1, ndata, sn, k);
      if (t != 0) {
         resdata[sres] = t;
         sres++;
      }
   }
   else {
      for (i = sn-1; i >= 0; i--)
         resdata1[i] = ndata[i];
   }

   for (i = 0; i < limb_cnt; i++)
      resdata[i] = 0;

   if (nneg) sres = -sres;
   SIZE(res) = sres;
}

void _ntl_grshift(_ntl_gbigint n, long k, _ntl_gbigint *rres)
{
   _ntl_gbigint res;
   mp_limb_t *ndata, *resdata, *ndata1;
   long limb_cnt, i, sn, nneg, sres;

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

   if (!k) {
      if (n != *rres)
         _ntl_gcopy(n, rres);
      return;
   }

   if (k < 0) {
      if (k < -NTL_MAX_LONG) ghalt("overflow in _ntl_glshift");
      _ntl_glshift(n, -k, rres);
      return;
   }

   GET_SIZE_NEG(sn, nneg, n);

   limb_cnt = k/NTL_ZZ_NBITS;

   sres = sn - limb_cnt;

   if (sres <= 0) {
      _ntl_gzero(rres);
      return;
   }

   res = *rres;
   if (MustAlloc(res, sres)) {
      _ntl_gsetlength(&res, sres);
      if (n == *rres) n = res;
      *rres = res;
   }

   ndata = DATA(n);
   resdata = DATA(res);
   ndata1 = ndata + limb_cnt;
   k %= NTL_ZZ_NBITS;

   if (k != 0) {
      mpn_rshift(resdata, ndata1, sres, k);
      if (resdata[sres-1] == 0)
         sres--;
   }
   else {
      for (i = 0; i < sres; i++)
         resdata[i] = ndata1[i];
   }

   if (nneg) sres = -sres;
   SIZE(res) = sres;
}
   
   

void
_ntl_gadd(_ntl_gbigint a, _ntl_gbigint b, _ntl_gbigint *cc)
{
   long sa, aneg, sb, bneg, sc, cmp;
   mp_limb_t *adata, *bdata, *cdata, carry;
   _ntl_gbigint c;

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

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

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

   if (sa < sb) {
      SWAP_BIGINT(a, b);
      SWAP_LONG(sa, sb);
      SWAP_LONG(aneg, bneg);
   }

   /* sa >= sb */

   c = *cc;

   if (aneg == bneg) {
      /* same sign => addition */

      sc = sa + 1;
      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);

      carry = mpn_add(cdata, adata, sa, bdata, sb);
      if (carry) 
         cdata[sc-1] = carry;
      else
         sc--;

      if (aneg) sc = -sc;
      SIZE(c) = sc;
   }
   else {
      /* opposite sign => subtraction */

      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) sc = -sc;
         SIZE(c) = sc;
      }
   }
}

void
_ntl_gsadd(_ntl_gbigint a, long b, _ntl_gbigint *cc)
{
   static _ntl_gbigint B = 0;
   _ntl_gintoz(b, &B);
   _ntl_gadd(a, B, cc);