zz.txt

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

/**************************************************************************\

MODULE: ZZ

SUMMARY:

The class ZZ is used to represent signed, arbitrary length integers.

Routines are provided for all of the basic arithmetic operations, as
well as for some more advanced operations such as primality testing.
Space is automatically managed by the constructors and destructors.

This module also provides routines for generating small primes, and
fast routines for performing modular arithmetic on single-precision
numbers.


\**************************************************************************/

#include <NTL/tools.h>


class ZZ {
public:

   ZZ(); // initial value is 0

   ZZ& operator=(const ZZ& a);  // assignment operator
   ZZ& operator=(long a);  

   ZZ(const ZZ& a);  // copy constructor
   ~ZZ(); // destructor


   // ...

};


// NOTE: A ZZ is represented as a sequence of "zzigits",
// where each zzigit is between 0 and 2^{NTL_ZZ_NBITS-1}.

// NTL_ZZ_NBITS is  macros defined in <NTL/ZZ.h>.

// SIZE INVARIANT: the number of bits in a ZZ is always less than
// 2^(NTL_BITS_PER_LONG-4).


/**************************************************************************\

                                 Comparison

\**************************************************************************/



// The usual comparison operators: 
   
long operator==(const ZZ& a, const ZZ& b);
long operator!=(const ZZ& a, const ZZ& b);
long operator<(const ZZ& a, const ZZ& b);
long operator>(const ZZ& a, const ZZ& b);
long operator<=(const ZZ& a, const ZZ& b);
long operator>=(const ZZ& a, const ZZ& b);

// other stuff:

long sign(const ZZ& a); // returns sign of a (-1, 0, +1)
long IsZero(const ZZ& a); // test for 0
long IsOne(const ZZ& a); // test for 1

long compare(const ZZ& a, const ZZ& b); // returns sign of a-b (-1, 0, or 1).

// PROMOTIONS: the comparison operators and the function compare
// support promotion from long to ZZ on (a, b).


/**************************************************************************\

                                 Addition

\**************************************************************************/


// operator notation:

ZZ operator+(const ZZ& a, const ZZ& b);
ZZ operator-(const ZZ& a, const ZZ& b);
ZZ operator-(const ZZ& a); // unary -

ZZ& operator+=(ZZ& x, const ZZ& a); 
ZZ& operator+=(ZZ& x, long a); 

ZZ& operator-=(ZZ& x, const ZZ& a); 
ZZ& operator-=(ZZ& x, long a); 

ZZ& operator++(ZZ& x);  // prefix
void operator++(ZZ& x, int);  // postfix

ZZ& operator--(ZZ& x);  // prefix
void operator--(ZZ& x, int);  // postfix



// procedural versions:

void add(ZZ& x, const ZZ& a, const ZZ& b); // x = a + b
void sub(ZZ& x, const ZZ& a, const ZZ& b); // x = a - b
void SubPos(ZZ& x, const ZZ& a, const ZZ& b); // x = a-b; assumes a >= b >= 0.
void negate(ZZ& x, const ZZ& a); // x = -a

void abs(ZZ& x, const ZZ& a); // x = |a|
ZZ abs(const ZZ& a);

// PROMOTIONS: binary +, -, as well as the procedural versions add, sub
// support promotions from long to ZZ on (a, b).


/**************************************************************************\

                             Multiplication

\**************************************************************************/

// operator notation:

ZZ operator*(const ZZ& a, const ZZ& b);

ZZ& operator*=(ZZ& x, const ZZ& a);
ZZ& operator*=(ZZ& x, long a);

// procedural versions:

void mul(ZZ& x, const ZZ& a, const ZZ& b); // x = a * b

void sqr(ZZ& x, const ZZ& a); // x = a*a
ZZ sqr(const ZZ& a); 

// PROMOTIONS: operator * and procedure mul support promotion
// from long to ZZ on (a, b).


/**************************************************************************\

                                 Division

\**************************************************************************/


// operator notation:

ZZ operator/(const ZZ& a, const ZZ& b);
ZZ operator/(const ZZ& a, long  b);

ZZ operator%(const ZZ& a, const ZZ& b);
long operator%(const ZZ& a, long b);

ZZ& operator/=(ZZ& x, const ZZ& b);
ZZ& operator/=(ZZ& x, long b);

ZZ& operator%=(ZZ& x, const ZZ& b);


// procedural versions:

void DivRem(ZZ& q, ZZ& r, const ZZ& a, const ZZ& b);
// q = floor(a/b), r = a - b*q.
// This implies that:
//    |r| < |b|, and if r != 0, sign(r) = sign(b)

void div(ZZ& q, const ZZ& a, const ZZ& b);
// q = floor(a/b)

void rem(ZZ& r, const ZZ& a, const ZZ& b);
// q = floor(a/b), r = a - b*q


// single-precision variants:

long DivRem(ZZ& q, const ZZ& a, long b);
// q = floor(a/b), r = a - b*q, return value is r.

long rem(const ZZ& a, long b);
// q = floor(a/b), r = a - b*q, return value is r.


// divisibility testing:

long divide(ZZ& q, const ZZ& a, const ZZ& b);
long divide(ZZ& q, const ZZ& a, long b);
// if b | a, sets q = a/b and returns 1; otherwise returns 0.

long divide(const ZZ& a, const ZZ& b);
long divide(const ZZ& a, long b);
// if b | a, returns 1; otherwise returns 0.


/**************************************************************************\

                                    GCD's

\**************************************************************************/


void GCD(ZZ& d, const ZZ& a, const ZZ& b);
ZZ GCD(const ZZ& a, const ZZ& b); 

// d = gcd(a, b) (which is always non-negative).  Uses a binary GCD
// algorithm.



void XGCD(ZZ& d, ZZ& s, ZZ& t, const ZZ& a, const ZZ& b);

//  d = gcd(a, b) = a*s + b*t.

// The coefficients s and t are defined according to the standard
// Euclidean algorithm applied to |a| and |b|, with the signs then
// adjusted according to the signs of a and b.

// Uses a variant of Lehmer's algorithm (see Knuth, The Art of Computer
// Programming, vol. 2).




// special-purpose single-precision variants:

long GCD(long a, long b);
// return value is gcd(a, b) (which is always non-negative)

void XGCD(long& d, long& s, long& t, long a, long b);
//  d = gcd(a, b) = a*s + b*t.

//  The coefficients s and t are defined according to the standard
//  Euclidean algorithm applied to |a| and |b|, with the signs then
//  adjusted according to the signs of a and b.



/**************************************************************************\

                             Modular Arithmetic

The following routines perform arithmetic mod n, where n > 1.

All arguments (other than exponents) are assumed to be in the range
0..n-1.  Some routines may check this and raise an error if this
does not hold.  Others may not, and the behaviour is unpredictable
in this case.

\**************************************************************************/



void AddMod(ZZ& x, const ZZ& a, const ZZ& b, const ZZ& n); // x = (a+b)%n
ZZ AddMod(const ZZ& a, const ZZ& b, const ZZ& n);

void SubMod(ZZ& x, const ZZ& a, const ZZ& b, const ZZ& n); // x = (a-b)%n
ZZ SubMod(const ZZ& a, const ZZ& b, const ZZ& n);

void NegateMod(ZZ& x, const ZZ& a, const ZZ& n); // x = -a % n
ZZ NegateMod(const ZZ& a, const ZZ& n);

void MulMod(ZZ& x, const ZZ& a, const ZZ& b, const ZZ& n); // x = (a*b)%n
ZZ MulMod(const ZZ& a, const ZZ& b, const ZZ& n);

void SqrMod(ZZ& x, const ZZ& a, const ZZ& n); // x = a^2 % n
ZZ SqrMod(const ZZ& a, const ZZ& n);

void InvMod(ZZ& x, const ZZ& a, const ZZ& n);
ZZ InvMod(const ZZ& a, const ZZ& n);
// x = a^{-1} mod n (0 <= x < n); error is raised occurs if inverse
// not defined

long InvModStatus(ZZ& x, const ZZ& a, const ZZ& n);
// if gcd(a,b) = 1, then return-value = 0, x = a^{-1} mod n;
// otherwise, return-value = 1, x = gcd(a, n)

void PowerMod(ZZ& x, const ZZ& a, const ZZ& e, const ZZ& n);
ZZ PowerMod(const ZZ& a, const ZZ& e, const ZZ& n);

void PowerMod(ZZ& x, const ZZ& a, long e, const ZZ& n);
ZZ PowerMod(const ZZ& a, long e, const ZZ& n);

// x = a^e % n (e may be negative)


// PROMOTIONS: AddMod, SubMod, and MulMod (both procedural and functional
// forms) support promotions from long to ZZ on (a, b).


/**************************************************************************\

                        Single-precision modular arithmetic

These routines implement single-precision modular arithmetic.  If n is
the modulus, all inputs should be in the range 0..n-1.  The number n
itself should be in the range 2..NTL_SP_BOUND-1.

Most of these routines are, of course, implemented as fast inline
functions.  No checking is done that inputs are in range.

\**************************************************************************/




long AddMod(long a, long b, long n); // return (a+b)%n

long SubMod(long a, long b, long n); // return (a-b)%n

long NegateMod(long a, long n); // return (-a)%n

long MulMod(long a, long b, long n); // return (a*b)%n

long MulMod(long a, long b, long n, double ninv);
// return (a*b)%n.  ninv = 1/((double) n).  This is faster if n is
// fixed for many multiplications.


long MulMod2(long a, long b, long n, double bninv);
// return (a*b)%n.  bninv = ((double) b)/((double) n).  This is faster
// if both n and b are fixed for many multiplications.


long MulDivRem(long& q, long a, long b, long n, double bninv);
// return (a*b)%n, set q = (a*b)/n.  bninv = ((double) b)/((double) n)

long InvMod(long a, long n);
// computes a^{-1} mod n.  Error is raised if undefined.

long PowerMod(long a, long e, long n);
// computes a^e mod n (e may be negative)

/**************************************************************************\

                               Shift Operations

LeftShift by n means multiplication by 2^n
RightShift by n means division by 2^n, with truncation toward zero
  (so the sign is preserved).

A negative shift amount reverses the direction of the shift.

\**************************************************************************/

// operator notation:

ZZ operator<<(const ZZ& a, long n);
ZZ operator>>(const ZZ& a, long n);

ZZ& operator<<=(ZZ& x, long n);
ZZ& operator>>=(ZZ& x, long n);

// procedural versions:

void LeftShift(ZZ& x, const ZZ& a, long n); 
ZZ LeftShift(const ZZ& a, long n);

void RightShift(ZZ& x, const ZZ& a, long n); 
ZZ RightShift(const ZZ& a, long n); 



/**************************************************************************\

                              Bits and Bytes

\**************************************************************************/



long MakeOdd(ZZ& x);
// removes factors of 2 from x, returns the number of 2's removed
// returns 0 if x == 0

long NumTwos(const ZZ& x);
// returns max e such that 2^e divides x if x != 0, and returns 0 if x == 0.

long IsOdd(const ZZ& a); // test if a is odd

long NumBits(const ZZ& a);
long NumBits(long a);  
// returns the number of bits in binary represenation of |a|; 
// NumBits(0) = 0


long bit(const ZZ& a, long k);
long bit(long a, long k); 
// returns bit k of a, position 0 being the low-order bit.
// If  k < 0 or k >= NumBits(a), returns 0.


void trunc(ZZ& x, const ZZ& a, long k);
// x = low order k bits of |a|. 
// If k <= 0, x = 0.

// two functional variants:
ZZ trunc_ZZ(const ZZ& a, long k);  
long trunc_long(const ZZ& a, long k);

long SetBit(ZZ& x, long p);
// returns original value of p-th bit of |a|, and replaces p-th bit of
// a by 1 if it was zero; low order bit is bit 0; error if p < 0;
// the sign of x is maintained

long SwitchBit(ZZ& x, long p);
// returns original value of p-th bit of |a|, and switches the value
// of p-th bit of a; low order bit is bit 0; error if p < 0
// the sign of x is maintained

long weight(const ZZ& a); // returns Hamming weight of |a|
long weight(long a); 

// bit-wise Boolean operations, procedural form: