g_lip_impl.h

数值算法库for Windows

#include <NTL/lip.h>
#include <NTL/IsFinite.h>


#include <stdlib.h>
#include <stdio.h>
#include <math.h>



#include <gmp.h>

typedef mp_limb_t *_ntl_limb_t_ptr;

int _ntl_gmp_hack = 0;

#if (__GNU_MP_VERSION < 3)

#error "You have to use GNP version >= 3.1"

#endif

#if ((__GNU_MP_VERSION == 3) && (__GNU_MP_VERSION_MINOR < 1))

#error "You have to use GNP version >= 3.1"

#endif



/* v 3.1 is supposed  mpn_tdiv_qr defined, but it doesn't.
   Here's a workaround */

#if ((__GNU_MP_VERSION == 3) && (__GNU_MP_VERSION_MINOR == 1) && (__GNU_MP_VERSION_PATCHLEVEL == 0))

#define mpn_tdiv_qr __MPN(tdiv_qr)


#ifdef __cplusplus
extern "C" 
#endif
void mpn_tdiv_qr(mp_ptr, mp_ptr, mp_size_t, mp_srcptr, mp_size_t, 
                 mp_srcptr, mp_size_t);

#endif


#if (defined(NTL_CXX_ONLY) && !defined(__cplusplus))
#error "CXX_ONLY flag set...must use C++ compiler"
#endif


union gbigint_header {

   long info[2];
   mp_limb_t alignment;

};

/* A bigint is represented as two long's, ALLOC and SIZE, followed by a 
 * vector DATA of mp_limb_t's.  
 * 
 * ALLOC is of the form
 *    (alloc << 2) | continue_flag | frozen_flag
 * where 
 *    - alloc is the number of allocated mp_limb_t's,
 *    - continue flag is either 2 or 0,
 *    - frozen_flag is either 1 or 0.
 * If frozen_flag is set, then the space for this bigint is *not*
 * managed by the _ntl_gsetlength and _ntl_gfree routines,
 * but are instead managed by the vec_ZZ_p and ZZVec routines.
 * The continue_flag is only set when the frozen_flag is set.
 * 
 * SIZE is the number of mp_limb_t's actually
 * used by the bigint, with the sign of SIZE having
 * the sign of the bigint.
 * Note that the zero bigint is represented as SIZE=0.
 * 
 * Bigint's are accessed through a handle, which is pointer to void.
 * A null handle logically represents the bigint zero.
 * This is done so that the interface presented to higher level
 * routines is essentially the same as that of NTL's traditional
 * long integer package.
 * 
 * The components ALLOC, SIZE, and DATA are all accessed through
 * macros using pointer casts.  While all of may seem a bit dirty, 
 * it should be quite portable: objects are never referenced
 * through pointers of different types, and no alignmement
 * problems should arise.
 * 
 * Actually, mp_limb_t is usually the type unsigned long.
 * However, on some 64-bit platforms, the type long is only 32 bits,
 * and gmp makes mp_limb_t unsigned long long in this case.
 * This is fairly rare, as the industry standard for Unix is to
 * have 64-bit longs on 64-bit machines.
 */ 

#if 1

#define ALLOC(p) (((long *) (p))[0])
#define SIZE(p) (((long *) (p))[1])
#define DATA(p) ((mp_limb_t *) (((long *) (p)) + 2))

#define STORAGE(len) ((long)(2*sizeof(long) + (len)*sizeof(mp_limb_t)))

/* STORAGE computes the number of bytes to allocate for a bigint
 * of maximal SIZE len.  This should be computed so that one
 * can store several such bigints in a contiguous array
 * of memory without breaking any alignment requirements.
 * Currently, it is assumed (and explicitly checked in the NTL installation
 * script) that sizeof(mp_limb_t) is either sizeof(long) or
 * 2*sizeof(long), and therfore, nothing special needs to
 * be done to enfoce alignment requirements.  If this assumption
 * should change, then the storage layout for bigints must be
 * re-designed.   However, this is very unlikely to ever happen.
 * 
 */

#define MustAlloc(c, len)  (!(c) || (ALLOC(c) >> 2) < (len))



#define GET_SIZE_NEG(sz, neg, p)  \
do  \
{   \
   long _s;   \
   _s = SIZE(p);   \
   if (_s < 0) {  \
      sz = -_s;  \
      neg = 1;  \
   }  \
   else {  \
      sz = _s;  \
      neg = 0;  \
   }  \
}  \
while (0)

#define STRIP(sz, p)  \
do  \
{  \
   long _i;  \
   _i = sz - 1;  \
   while (_i >= 0 && p[_i] == 0) _i--;  \
   sz = _i + 1;  \
}  \
while (0) 

#define ZEROP(p) (!p || !SIZE(p))

#define ONEP(p) (p && SIZE(p) == 1 && DATA(p)[0] == 1)

#define SWAP_BIGINT(a, b)  \
do  \
{  \
   _ntl_gbigint _t;  \
   _t = a;  \
   a = b;  \
   b = _t;  \
}  \
while (0) 

#define SWAP_LONG(a, b)  \
do  \
{  \
   long _t;  \
   _t = a;  \
   a = b;  \
   b = _t;  \
}  \
while (0) 

#define SWAP_LIMB_PTR(a, b)  \
do  \
{  \
   _ntl_limb_t_ptr _t;  \
   _t = a;  \
   a = b;  \
   b = _t;  \
}  \
while (0) 

#define COUNT_BITS(cnt, a)  \
do  \
{  \
   long _i = 0;  \
   mp_limb_t _a = (a);  \
  \
   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++;  \
  \
   cnt = _i;  \
}  \
while (0) 

#else

/* These are C++ inline functions that are equivalent to the above
 * macros.  They are mainly intended as a debugging aid.
 */


inline long& ALLOC(_ntl_gbigint p) 
   { return (((long *) p)[0]); }

inline long& SIZE(_ntl_gbigint p) 
   { return (((long *) p)[1]); }

inline mp_limb_t * DATA(_ntl_gbigint p) 
   { return ((mp_limb_t *) (((long *) (p)) + 2)); }

inline long STORAGE(long len)
   { return ((long)(2*sizeof(long) + (len)*sizeof(mp_limb_t))); }

inline long MustAlloc(_ntl_gbigint c, long len)  
   { return (!(c) || (ALLOC(c) >> 2) < (len)); }


inline void GET_SIZE_NEG(long& sz, long& neg, _ntl_gbigint p)
{ 
   long s; 
   s = SIZE(p); 
   if (s < 0) {
      sz = -s;
      neg = 1;
   }
   else {
      sz = s;
      neg = 0;
   }
}

inline void STRIP(long& sz, mp_limb_t *p)
{
   long i;
   i = sz - 1;
   while (i >= 0 && p[i] == 0) i--;
   sz = i + 1;
}

inline long ZEROP(_ntl_gbigint p)
{
   return !p || !SIZE(p);
}

inline long ONEP(_ntl_gbigint p)
{
   return p && SIZE(p) == 1 && DATA(p)[0] == 1;
}

inline void SWAP_BIGINT(_ntl_gbigint& a, _ntl_gbigint& b)
{
   _ntl_gbigint t;
   t = a;
   a = b;
   b = t;
}

inline void SWAP_LONG(long& a, long& b)
{
   long t;
   t = a;
   a = b;
   b = t;
}

inline void SWAP_LIMB_PTR(_ntl_limb_t_ptr& a, _ntl_limb_t_ptr& b)
{
   _ntl_limb_t_ptr t;
   t = a;
   a = b;
   b = t;
}


inline void COUNT_BITS(long& cnt, mp_limb_t a)
{
   long i = 0;

   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++;

   cnt = i;
}

#endif



static 
void ghalt(char *c)
{
   fprintf(stderr,"fatal error:\n   %s\nexit...\n",c);
   fflush(stderr);
   abort();
}


#define MIN_SETL	(4)
   /* _ntl_gsetlength allocates a multiple of MIN_SETL digits */



void _ntl_gsetlength(_ntl_gbigint *v, long len)
{
   _ntl_gbigint x = *v;

   if (len < 0)
      ghalt("negative size allocation in _ntl_zgetlength");

   if (len >= (1L << (NTL_BITS_PER_LONG-4))/NTL_ZZ_NBITS)
      ghalt("size too big in _ntl_gsetlength");

#ifdef NTL_SMALL_MP_SIZE_T
   /* this makes sure that numbers don't get too big for GMP */
   if (len >= (1L << (NTL_BITS_PER_INT-4)))
      ghalt("size too big for GMP");
#endif


   if (x) {
      long oldlen = ALLOC(x);
      long fixed = oldlen & 1;
      oldlen = oldlen >> 2;

      if (fixed) {
         if (len > oldlen) 
            ghalt("internal error: can't grow this _ntl_gbigint");
         else
            return;
      }

      if (len <= oldlen) return;

      len++;  /* always allocate at least one more than requested */

      oldlen = (long) (oldlen * 1.2); /* always increase by at least 20% */
      if (len < oldlen)
         len = oldlen;

      /* round up to multiple of MIN_SETL */
      len = ((len+(MIN_SETL-1))/MIN_SETL)*MIN_SETL; 

      /* test len again */
      if (len >= (1L << (NTL_BITS_PER_LONG-4))/NTL_ZZ_NBITS)
         ghalt("size too big in _ntl_gsetlength");

      ALLOC(x) = len << 2;
      if (!(x = (_ntl_gbigint)realloc((void*) x, STORAGE(len)))) {
         ghalt("reallocation failed in _ntl_gsetlength");
      }
   }
   else {
      len++;
      len = ((len+(MIN_SETL-1))/MIN_SETL)*MIN_SETL;

      /* test len again */
      if (len >= (1L << (NTL_BITS_PER_LONG-4))/NTL_ZZ_NBITS)
         ghalt("size too big in _ntl_gsetlength");

      if (!(x = (_ntl_gbigint)malloc(STORAGE(len)))) {
         ghalt("allocation failed in _ntl_gsetlength");
      }
      ALLOC(x) = len << 2;
      SIZE(x) = 0;
   }

   *v = x;
}

void _ntl_gfree(_ntl_gbigint *xx)
{
   _ntl_gbigint x = *xx;

   if (!x)
      return;

   if (ALLOC(x) & 1)
      ghalt("Internal error: can't free this _ntl_gbigint");

   free((void*) x);
   *xx = 0;
   return;
}

void
_ntl_gswap(_ntl_gbigint *a, _ntl_gbigint *b)
{
   _ntl_gbigint c;

   if ((*a && (ALLOC(*a) & 1)) || (*b && (ALLOC(*b) & 1))) {
      static _ntl_gbigint t = 0; 
      _ntl_gcopy(*a, &t);
      _ntl_gcopy(*b, a);
      _ntl_gcopy(t, b);
      return;
   }

   c = *a;
   *a = *b;
   *b = c;
}


void _ntl_gcopy(_ntl_gbigint a, _ntl_gbigint *bb)
{
   _ntl_gbigint b;
   long sa, abs_sa, i;
   mp_limb_t *adata, *bdata;

   b = *bb;

   if (!a || (sa = SIZE(a)) == 0) {
      if (b) SIZE(b) = 0;
   }
   else {
      if (a != b) {
         if (sa >= 0)
            abs_sa = sa;
         else
            abs_sa = -sa;

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

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

         for (i = 0; i < abs_sa; i++)
            bdata[i] = adata[i];

         SIZE(b) = sa;
      }
   }
}


void _ntl_gzero(_ntl_gbigint *aa) 
{
   _ntl_gbigint a = *aa;

   if (a) SIZE(a) = 0;
}

void _ntl_gone(_ntl_gbigint *aa)
{
   _ntl_gbigint a = *aa;
   if (!a) {
      _ntl_gsetlength(&a, 1);
      *aa = a;
   }

   SIZE(a) = 1;
   DATA(a)[0] = 1;
}

long _ntl_giszero(_ntl_gbigint a)
{
   return ZEROP(a);
}

long _ntl_godd(_ntl_gbigint a)
{
   if (ZEROP(a)) 
      return 0;
   else
      return DATA(a)[0]&1;
}

long _ntl_gbit(_ntl_gbigint a, long p)
{
   long bl;
   long sa;
   mp_limb_t wh;

   if (p < 0 || !a) return 0;

   bl = p/NTL_ZZ_NBITS;
   wh = ((mp_limb_t) 1) << (p - NTL_ZZ_NBITS*bl);

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

   if (sa <= bl) return 0;
   if (DATA(a)[bl] & wh) return 1;
   return 0;
}

void _ntl_glowbits(_ntl_gbigint a, long b, _ntl_gbigint *cc)
{
   _ntl_gbigint c = *cc;

   long bl;
   long wh;
   long sa;
   long i;
   mp_limb_t *adata, *cdata;

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

   bl = b/NTL_ZZ_NBITS;
   wh = b - NTL_ZZ_NBITS*bl;
   if (wh != 0) 
      bl++;
   else
      wh = NTL_ZZ_NBITS;

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

   if (sa < bl) {
      _ntl_gcopy(a,cc);
      _ntl_gabs(cc);
      return;
   }

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

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

   for (i = 0; i < bl-1; i++)
      cdata[i] = adata[i];

   if (wh == NTL_ZZ_NBITS)
      cdata[bl-1] = adata[bl-1];
   else
      cdata[bl-1] = adata[bl-1] & ((((mp_limb_t) 1) << wh) - ((mp_limb_t) 1));

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

long _ntl_gslowbits(_ntl_gbigint a, long p)
{
   static _ntl_gbigint x = 0;

   _ntl_glowbits(a, p, &x);

   return _ntl_gtoint(x);
}

long _ntl_gsetbit(_ntl_gbigint *a, long b)
{
   long bl;
   long sa, aneg;
   long i;
   mp_limb_t wh, *adata, tmp;

   if (b<0) ghalt("_ntl_gsetbit: negative index");

   if (ZEROP(*a)) {
      _ntl_gintoz(1, a);
      _ntl_glshift(*a, b, a);
      return 0;
   }

   bl = (b/NTL_ZZ_NBITS);
   wh = ((mp_limb_t) 1) << (b - NTL_ZZ_NBITS*bl);

   GET_SIZE_NEG(sa, aneg, *a);

   if (sa > bl) {
      adata = DATA(*a);
      tmp = adata[bl] & wh;
      adata[bl] |= wh;
      if (tmp) return 1;
      return 0;
   } 
   else {
      _ntl_gsetlength(a, bl+1);
      adata = DATA(*a);
      for (i = sa; i < bl; i++)
         adata[i] = 0;
      adata[bl] = wh;

      sa = bl+1;
      if (aneg) sa = -sa;
      SIZE(*a) = sa;
      return 0;
   }
}

long _ntl_gswitchbit(_ntl_gbigint *a, long b)
{
   long bl;
   long sa, aneg;
   long i;
   mp_limb_t wh, *adata, tmp;

   if (b<0) ghalt("_ntl_gswitchbit: negative index");


   if (ZEROP(*a)) {
      _ntl_gintoz(1, a);
      _ntl_glshift(*a, b, a);
      return 0;
   }

   bl = (b/NTL_ZZ_NBITS);
   wh = ((mp_limb_t) 1) << (b - NTL_ZZ_NBITS*bl);

   GET_SIZE_NEG(sa, aneg, *a);

   if (sa > bl) {
      adata = DATA(*a);
      tmp = adata[bl] & wh;
      adata[bl] ^= wh;

      if (bl == sa-1) {
         STRIP(sa, adata);
         if (aneg) sa = -sa;
         SIZE(*a) = sa;
      }

      if (tmp) return 1;
      return 0;
   } 
   else {
      _ntl_gsetlength(a, bl+1);
      adata = DATA(*a);
      for (i = sa; i < bl; i++)
         adata[i] = 0;
      adata[bl] = wh;

      sa = bl+1;
      if (aneg) sa = -sa;
      SIZE(*a) = sa;
      return 0;
   }
}

long
_ntl_gweights(
	long aa
	)
{
	unsigned long a;
	long res = 0;
	if (aa < 0) 
		a = -aa;
	else
		a = aa;
   
	while (a) {
		if (a & 1) res ++;
		a >>= 1;
	}
	return (res);