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:

void bit_and(ZZ& x, const ZZ& a, const ZZ& b); // x = |a| AND |b|
void bit_or(ZZ& x, const ZZ& a, const ZZ& b); // x = |a| OR |b|
void bit_xor(ZZ& x, const ZZ& a, const ZZ& b); // x = |a| XOR |b|

// bit-wise Boolean operations, operator notation:

ZZ operator&(const ZZ& a, const ZZ& b);
ZZ operator|(const ZZ& a, const ZZ& b);
ZZ operator^(const ZZ& a, const ZZ& b);

// PROMOTIONS: the above bit-wise operations (both procedural 
// and operator forms) provide promotions from long to ZZ on (a, b).

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

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

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



// conversions between byte sequences and ZZ's

void ZZFromBytes(ZZ& x, const unsigned char *p, long n);
ZZ ZZFromBytes(const unsigned char *p, long n);
// x = sum(p[i]*256^i, i=0..n-1). 
// NOTE: in the unusual event that a char is more than 8 bits, 
//       only the low order 8 bits of p[i] are used

void BytesFromZZ(unsigned char *p, const ZZ& a, long n);
// Computes p[0..n-1] such that abs(a) == sum(p[i]*256^i, i=0..n-1) mod 256^n.

long NumBytes(const ZZ& a);
long NumBytes(long a);
// returns # of base 256 digits needed to represent abs(a).
// NumBytes(0) == 0.



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

                            Pseudo-Random Numbers

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


// Routines for generating pseudo-random numbers.

// These routines generate high qualtity, cryptographically strong
// pseudo-random numbers.  They are implemented so that their behaviour
// is completely independent of the underlying hardware and long 
// integer implementation.  Note, however, that other routines 
// throughout NTL use pseudo-random numbers, and because of this,
// the word size of the machine can impact the sequence of numbers
// seen by a client program.


void SetSeed(const ZZ& s); 
// Initializes generator with a "seed" s.
// s is first hashed to generate the initial state, so it is
// not necessary that s itself looks random, just that 
// it has a lot of "entropy".
// If SetSeed is not called before using the routines below,
// a default initial state is used.
// Calling SetSeed with s == 0, e.g. SetSeed(ZZ::zero()), 
// has the effect of re-setting the state to the default initial state.
// Routine ZZFromBytes (above) may be useful.


void RandomBnd(ZZ& x, const ZZ& n);
ZZ RandomBnd(const ZZ& n);
long RandomBnd(long n);
// x = pseudo-random number in the range 0..n-1, or 0 if n <= 0

void RandomBits(ZZ& x, long l);
ZZ RandomBits_ZZ(long l);
long RandomBits_long(long l);
// x = pseudo-random number in the range 0..2^l-1.

void RandomLen(ZZ& x, long l);
ZZ RandomLen_ZZ(long l);
long RandomLen_long(long l);
// x = psuedo-random number with precisely l bits,
// or 0 of l <= 0.

unsigned long RandomBits_ulong(long l);
// returns a pseudo-random number in the range 0..2^l-1

unsigned long RandomWord();
// returns a word filled with pseudo-random bits.
// Equivalent to RandomBits_ulong(NTL_BITS_PER_LONG).


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

             Incremental Chinese Remaindering

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

long CRT(ZZ& a, ZZ& p, const ZZ& A, const ZZ& P);
long CRT(ZZ& a, ZZ& p, long A, long P);

// 0 <= A < P, (p, P) = 1; computes a' such that a' = a mod p, 
// a' = A mod P, and -p*P/2 < a' <= p*P/2; sets a := a', p := p*P, and
// returns 1 if a's value has changed, otherwise 0


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

                  Rational Reconstruction

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

long ReconstructRational(ZZ& a, ZZ& b, const ZZ& x, const ZZ& m, 
                         const ZZ& a_bound, const ZZ& b_bound);

// 0 <= x < m, m > 2 * a_bound * b_bound,
// a_bound >= 0, b_bound > 0

// This routine either returns 0, leaving a and b unchanged, 
// or returns 1 and sets a and b so that
//   (1) a = b x (mod m),
//   (2) |a| <= a_bound, 0 < b <= b_bound, and
//   (3) gcd(m, b) = gcd(a, b).

// If there exist a, b satisfying (1), (2), and 
//   (3') gcd(m, b) = 1,
// then a, b are uniquely determined if we impose the additional
// condition that gcd(a, b) = 1;  moreover, if such a, b exist,
// then these values are returned by the routine.

// Unless the calling routine can *a priori* guarantee the existence
// of a, b satisfying (1), (2), and (3'),
// then to ensure correctness, the calling routine should check
// that gcd(m, b) = 1, or equivalently, gcd(a, b) = 1.

// This is implemented using a variant of Lehmer's extended
// Euclidean algorithm.

// Literature:  see G. Collins and M. Encarnacion, J. Symb. Comp. 20:287-297, 
// 1995; P. Wang, M. Guy, and J. Davenport, SIGSAM Bulletin 16:2-3, 1982. 


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

                                Primality Testing 
                           and Prime Number Generation

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

void GenPrime(ZZ& n, long l, long err = 80);
ZZ GenPrime_ZZ(long l, long err = 80);
long GenPrime_long(long l, long err = 80);

// GenPrime generates a random prime n of length l so that the
// probability that the resulting n is composite is bounded by 2^(-err).
// This calls the routine RandomPrime below, and uses results of 
// Damgard, Landrock, Pomerance to "optimize" 
// the number of Miller-Rabin trials at the end.

void GenGermainPrime(ZZ& n, long l, long err = 80);
ZZ GenGermainPrime_ZZ(long l, long err = 80);
long GenGermainPrime_long(long l, long err = 80);

// A (Sophie) Germain prime is a prime p such that p' = 2*p+1 is also a prime.
// Such primes are useful for cryptographic applications...cryptographers
// sometimes call p' a "strong" or "safe" prime.
// GenGermainPrime generates a random Germain prime n of length l
// so that the probability that either n or 2*n+1 is not a prime
// is bounded by 2^(-err).


long ProbPrime(const ZZ& n, long NumTrials = 10);
long ProbPrime(long n, long NumTrials = 10);
// performs up to NumTrials Miller-witness tests (after some trial division).

void RandomPrime(ZZ& n, long l, long NumTrials=10);
ZZ RandomPrime_ZZ(long l, long NumTrials=10);
long RandomPrime_long(long l, long NumTrials=10);
// n = random l-bit prime.  Uses ProbPrime with NumTrials.

void NextPrime(ZZ& n, const ZZ& m, long NumTrials=10);
ZZ NextPrime(const ZZ& m, long NumTrials=10);
// n = smallest prime >= m.  Uses ProbPrime with NumTrials.

long NextPrime(long m, long NumTrials=10);
// Single precision version of the above.
// Result will always be bounded by NTL_ZZ_SP_BOUND, and an
// error is raised if this cannot be satisfied.

long MillerWitness(const ZZ& n, const ZZ& w);
// Tests if w is a witness to compositeness a la Miller.  Assumption: n is
// odd and positive, 0 <= w < n.
// Return value of 1 implies n is composite.
// Return value of 0 indicates n might be prime.


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

                               Exponentiation

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


void power(ZZ& x, const ZZ& a, long e); // x = a^e (e >= 0)
ZZ power(const ZZ& a, long e); 

void power(ZZ& x, long a, long e);

// two functional variants:
ZZ power_ZZ(long a, long e);
long power_long(long a, long e);

void power2(ZZ& x, long e); // x = 2^e (e >= 0)
ZZ power2_ZZ(long e);


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

                               Square Roots

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


void SqrRoot(ZZ& x, const ZZ& a); // x = floor(a^{1/2}) (a >= 0)
ZZ SqrRoot(const ZZ& a); 

long SqrRoot(long a); 




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

                    Jacobi symbol and modular square roots

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


long Jacobi(const ZZ& a, const ZZ& n);
//  compute Jacobi symbol of a and n; assumes 0 <= a < n, n odd

void SqrRootMod(ZZ& x, const ZZ& a, const ZZ& n);
ZZ SqrRootMod(const ZZ& a, const ZZ& n); 
//  computes square root of a mod n; assumes n is an odd prime, and
//  that a is a square mod n, with 0 <= a < n.




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

                             Input/Output

I/O Format:

Numbers are written in base 10, with an optional minus sign.

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

istream& operator>>(istream& s, ZZ& x);  
ostream& operator<<(ostream& s, const ZZ& a); 



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

                            Miscellany

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


// The following macros are defined:

#define NTL_ZZ_NBITS (...)  // number of bits in a zzigit;
                            // a ZZ is represented as a sequence of zzigits.

#define NTL_SP_NBITS (...)  // max number of bits in a "single-precision" number

#define NTL_WSP_NBITS (...)  // max number of bits in a "wide single-precision"
                             // number

// The following relations hold:
//    NTL_SP_NBITS <= NTL_WSP_NBITS <= NTL_ZZ_NBITS
//    26 <= NTL_SP_NBITS <= min(NTL_BITS_PER_LONG-2, NTL_DOUBLE_PRECISION-3)
//    NTL_WSP_NBITS <= NTL_BITS_PER_LONG-2
//
// Note that NTL_ZZ_NBITS may be less than, equal to, or greater than
// NTL_BITS_PER_LONG  -- no particular relationship should be assumed to hold.
// In particular, expressions like (1L << NTL_ZZ_BITS) might overflow.
//
// "single-precision" numbers are meant to be used in conjunction with the
//  single-precision modular arithmetic routines.
//
// "wide single-precision" numbers are meant to be used in conjunction
//  with the ZZ arithmetic routines for optimal efficiency.

// The following auxilliary macros are also defined

#define NTL_FRADIX (...) // double-precision value of 2^NTL_ZZ_NBITS

#define NTL_SP_BOUND (1L << NTL_SP_NBITS)
#define NTL_WSP_BOUND (1L << NTL_WSP_NBITS)


// Backward compatability note:
// Prior to version 5.0, the macro NTL_NBITS was defined,
// along with the macro NTL_RADIX defined to be (1L << NTL_NBITS).
// While these macros are still available when using NTL's traditional 
// long integer package (i.e., when NTL_GMP_LIP is not set), 
// they are not available when using the GMP as the primary long integer 
// package (i.e., when NTL_GMP_LIP is set).
// Furthermore, when writing portable programs, one should avoid these macros.
// Note that when using traditional long integer arithmetic, we have
//    NTL_ZZ_NBITS = NTL_SP_NBITS = NTL_WSP_NBITS = NTL_NBITS.


// Here are some additional functions.

void clear(ZZ& x); // x = 0
void set(ZZ& x);   // x = 1

void swap(ZZ& x, ZZ& y);
// swap x and y (done by "pointer swapping", if possible).

double log(const ZZ& a);
// returns double precision approximation to log(a)

long NextPowerOfTwo(long m);
// returns least nonnegative k such that 2^k >= m

long ZZ::size() const;
// a.size() returns the number of zzigits of |a|; the
// size of 0 is 0.

void ZZ::SetSize(long k)
// a.SetSize(k) does not change the value of a, but simply pre-allocates
// space for k zzigits.

long ZZ::SinglePrecision() const;
// a.SinglePrecision() is a predicate that tests if abs(a) < NTL_SP_BOUND

long ZZ::WideSinglePrecision() const;
// a.WideSinglePrecision() is a predicate that tests if abs(a) < NTL_WSP_BOUND

long digit(const ZZ& a, long k);
// returns k-th zzigit of |a|, position 0 being the low-order
// zzigit.
// NOTE: this routine is only available when using NTL's traditional
// long integer arithmetic, and should not be used in programs
// that are meant to be portable.

void ZZ::kill();
// a.kill() sets a to zero and frees the space held by a.

ZZ::ZZ(INIT_SIZE_TYPE, long k);
// ZZ(INIT_SIZE, k) initializes to 0, but space is pre-allocated so
// that numbers x with x.size() <= k can be stored without
// re-allocation.

static const ZZ& ZZ::zero();
// ZZ::zero() yields a read-only reference to zero, if you need it.




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

                    Small Prime Generation

primes are generated in sequence, starting at 2, and up to a maximum
that is no more than min(NTL_SP_BOUND, 2^30).

Example: print the primes up to 1000

#include <NTL/ZZ.h>

main()
{
   PrimeSeq s;
   long p;

   p = s.next();
   while (p <= 1000) {
      cout << p << "\n";
      p = s.next();
   }
}

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



class PrimeSeq {
public:
   PrimeSeq();
   ~PrimeSeq();

   long next();
   // returns next prime in the sequence.  returns 0 if list of small
   // primes is exhausted.

   void reset(long b);
   // resets generator so that the next prime in the sequence is the
   // smallest prime >= b.

private:
   PrimeSeq(const PrimeSeq&);        // disabled
   void operator=(const PrimeSeq&);  // disabled

};