lzz_px.h

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

inline zz_pX operator-(const zz_pX& a)
   { zz_pX x; negate(x, a); NTL_OPT_RETURN(zz_pX, x); }

inline zz_pX& operator++(zz_pX& x) { add(x, x, 1); return x; }
inline void operator++(zz_pX& x, int) { add(x, x, 1); }
inline zz_pX& operator--(zz_pX& x) { sub(x, x, 1); return x; }
inline void operator--(zz_pX& x, int) { sub(x, x, 1); }



/*****************************************************************

                        Multiplication

******************************************************************/


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

void sqr(zz_pX& x, const zz_pX& a);
inline zz_pX sqr(const zz_pX& a)
   { zz_pX x; sqr(x, a); NTL_OPT_RETURN(zz_pX, x); }
// x = a^2

void mul(zz_pX& x, const zz_pX& a, zz_p b);
inline void mul(zz_pX& x, const zz_pX& a, long b) { mul(x, a, to_zz_p(b)); }

inline void mul(zz_pX& x, zz_p a, const zz_pX& b) { mul(x, b, a); }
inline void mul(zz_pX& x, long a, const zz_pX& b) { mul(x, b, a); }


inline zz_pX operator*(const zz_pX& a, const zz_pX& b)
   { zz_pX x; mul(x, a, b); NTL_OPT_RETURN(zz_pX, x); }

inline zz_pX operator*(const zz_pX& a, zz_p b)
   { zz_pX x; mul(x, a, b); NTL_OPT_RETURN(zz_pX, x); }

inline zz_pX operator*(const zz_pX& a, long b)
   { zz_pX x; mul(x, a, b); NTL_OPT_RETURN(zz_pX, x); }

inline zz_pX operator*(zz_p a, const zz_pX& b)
   { zz_pX x; mul(x, a, b); NTL_OPT_RETURN(zz_pX, x); }

inline zz_pX operator*(long a, const zz_pX& b)
   { zz_pX x; mul(x, a, b); NTL_OPT_RETURN(zz_pX, x); }

inline zz_pX& operator*=(zz_pX& x, const zz_pX& b)
   { mul(x, x, b); return x; }

inline zz_pX& operator*=(zz_pX& x, zz_p b)
   { mul(x, x, b); return x; }

inline zz_pX& operator*=(zz_pX& x, long b)
   { mul(x, x, b); return x; }


void PlainMul(zz_pX& x, const zz_pX& a, const zz_pX& b);
// always uses the "classical" algorithm

void PlainSqr(zz_pX& x, const zz_pX& a);
// always uses the "classical" algorithm


void FFTMul(zz_pX& x, const zz_pX& a, const zz_pX& b);
// always uses the FFT

void FFTSqr(zz_pX& x, const zz_pX& a);
// always uses the FFT

void MulTrunc(zz_pX& x, const zz_pX& a, const zz_pX& b, long n);
// x = a * b % X^n

inline zz_pX MulTrunc(const zz_pX& a, const zz_pX& b, long n)
   { zz_pX x; MulTrunc(x, a, b, n); NTL_OPT_RETURN(zz_pX, x); }

void PlainMulTrunc(zz_pX& x, const zz_pX& a, const zz_pX& b, long n);
void FFTMulTrunc(zz_pX& x, const zz_pX& a, const zz_pX& b, long n);

void SqrTrunc(zz_pX& x, const zz_pX& a, long n);
// x = a^2 % X^n

inline zz_pX SqrTrunc(const zz_pX& a, long n)
   { zz_pX x; SqrTrunc(x, a, n); NTL_OPT_RETURN(zz_pX, x); }

void PlainSqrTrunc(zz_pX& x, const zz_pX& a, long n);
void FFTSqrTrunc(zz_pX& x, const zz_pX& a, long n);

void power(zz_pX& x, const zz_pX& a, long e);
inline zz_pX power(const zz_pX& a, long e)
   { zz_pX x; power(x, a, e); NTL_OPT_RETURN(zz_pX, x); }






// The following data structures and routines allow one
// to hand-craft various algorithms, using the FFT convolution
// algorithms directly.
// Look in the file zz_pX.c for examples.




// FFT representation of polynomials

class fftRep {

public:
   long k;                // a 2^k point representation
   long MaxK;             // maximum space allocated
   long *tbl[4];
   long NumPrimes;

   fftRep(const fftRep&); 
   fftRep& operator=(const fftRep&); 

   void SetSize(long NewK);

   fftRep() { k = MaxK = -1; NumPrimes = zz_pInfo->NumPrimes; }
   fftRep(INIT_SIZE_TYPE, long InitK) 
   { k = MaxK = -1; NumPrimes = zz_pInfo->NumPrimes; SetSize(InitK); }
   ~fftRep();
};


void TofftRep(fftRep& y, const zz_pX& x, long k, long lo, long hi);
// computes an n = 2^k point convolution of x[lo..hi].

inline void TofftRep(fftRep& y, const zz_pX& x, long k)

   { TofftRep(y, x, k, 0, deg(x)); }

void RevTofftRep(fftRep& y, const vec_zz_p& x,
                 long k, long lo, long hi, long offset);
// computes an n = 2^k point convolution of X^offset*x[lo..hi] mod X^n-1
// using "inverted" evaluation points.



void FromfftRep(zz_pX& x, fftRep& y, long lo, long hi);
// converts from FFT-representation to coefficient representation
// only the coefficients lo..hi are computed
// NOTE: this version destroys the data in y

// non-destructive versions of the above

void NDFromfftRep(zz_pX& x, const fftRep& y, long lo, long hi, fftRep& temp);
void NDFromfftRep(zz_pX& x, const fftRep& y, long lo, long hi);

void RevFromfftRep(vec_zz_p& x, fftRep& y, long lo, long hi);

   // converts from FFT-representation to coefficient representation
   // using "inverted" evaluation points.
   // only the coefficients lo..hi are computed




void FromfftRep(zz_p* x, fftRep& y, long lo, long hi);
// convert out coefficients lo..hi of y, store result in x.
// no normalization is done.


// direct manipulation of FFT reps

void mul(fftRep& z, const fftRep& x, const fftRep& y);
void sub(fftRep& z, const fftRep& x, const fftRep& y);
void add(fftRep& z, const fftRep& x, const fftRep& y);

void reduce(fftRep& x, const fftRep& a, long k);
// reduces a 2^l point FFT-rep to a 2^k point FFT-rep

void AddExpand(fftRep& x, const fftRep& a);
//  x = x + (an "expanded" version of a)







/*************************************************************

                      Division

**************************************************************/

void DivRem(zz_pX& q, zz_pX& r, const zz_pX& a, const zz_pX& b);
// q = a/b, r = a%b

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


void div(zz_pX& q, const zz_pX& a, zz_p b);
inline void div(zz_pX& q, const zz_pX& a, long b)
   { div(q, a, to_zz_p(b)); }

void rem(zz_pX& r, const zz_pX& a, const zz_pX& b);
// r = a%b

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

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


void InvTrunc(zz_pX& x, const zz_pX& a, long m);
// computes x = a^{-1} % X^m 
// constant term must be non-zero

inline zz_pX InvTrunc(const zz_pX& a, long m)
   { zz_pX x; InvTrunc(x, a, m); NTL_OPT_RETURN(zz_pX, x); }




// These always use "classical" arithmetic
void PlainDivRem(zz_pX& q, zz_pX& r, const zz_pX& a, const zz_pX& b);
void PlainDiv(zz_pX& q, const zz_pX& a, const zz_pX& b);
void PlainRem(zz_pX& r, const zz_pX& a, const zz_pX& b);


// These always use FFT arithmetic
void FFTDivRem(zz_pX& q, zz_pX& r, const zz_pX& a, const zz_pX& b);
void FFTDiv(zz_pX& q, const zz_pX& a, const zz_pX& b);
void FFTRem(zz_pX& r, const zz_pX& a, const zz_pX& b);

void PlainInvTrunc(zz_pX& x, const zz_pX& a, long m);
// always uses "classical" algorithm
// ALIAS RESTRICTION: input may not alias output

void NewtonInvTrunc(zz_pX& x, const zz_pX& a, long m);
// uses a Newton Iteration with the FFT.
// ALIAS RESTRICTION: input may not alias output


inline zz_pX operator/(const zz_pX& a, const zz_pX& b)
   { zz_pX x; div(x, a, b); NTL_OPT_RETURN(zz_pX, x); }

inline zz_pX operator/(const zz_pX& a, zz_p b)
   { zz_pX x; div(x, a, b); NTL_OPT_RETURN(zz_pX, x); }

inline zz_pX operator/(const zz_pX& a, long b)
   { zz_pX x; div(x, a, b); NTL_OPT_RETURN(zz_pX, x); }

inline zz_pX& operator/=(zz_pX& x, zz_p b)
   { div(x, x, b); return x; }

inline zz_pX& operator/=(zz_pX& x, long b)
   { div(x, x, b); return x; }

inline zz_pX& operator/=(zz_pX& x, const zz_pX& b)
   { div(x, x, b); return x; }


inline zz_pX operator%(const zz_pX& a, const zz_pX& b)
   { zz_pX x; rem(x, a, b); NTL_OPT_RETURN(zz_pX, x); }

inline zz_pX& operator%=(zz_pX& x, const zz_pX& b)
   { rem(x, x, b); return x; }




/***********************************************************

                         GCD's

************************************************************/


void GCD(zz_pX& x, const zz_pX& a, const zz_pX& b);
// x = GCD(a, b),  x is always monic (or zero if a==b==0).

inline zz_pX GCD(const zz_pX& a, const zz_pX& b)
   { zz_pX x; GCD(x, a, b); NTL_OPT_RETURN(zz_pX, x); }


void XGCD(zz_pX& d, zz_pX& s, zz_pX& t, const zz_pX& a, const zz_pX& b);
// d = gcd(a,b), a s + b t = d 

void PlainXGCD(zz_pX& d, zz_pX& s, zz_pX& t, const zz_pX& a, const zz_pX& b);
// same as above, but uses classical algorithm


void PlainGCD(zz_pX& x, const zz_pX& a, const zz_pX& b);
// always uses "cdlassical" arithmetic


class zz_pXMatrix {
private:

   zz_pXMatrix(const zz_pXMatrix&);  // disable
   zz_pX elts[2][2];

public:

   zz_pXMatrix() { }
   ~zz_pXMatrix() { }

   void operator=(const zz_pXMatrix&);
   zz_pX& operator() (long i, long j) { return elts[i][j]; }
   const zz_pX& operator() (long i, long j) const { return elts[i][j]; }
};


void HalfGCD(zz_pXMatrix& M_out, const zz_pX& U, const zz_pX& V, long d_red);
// deg(U) > deg(V),   1 <= d_red <= deg(U)+1.
//
// This computes a 2 x 2 polynomial matrix M_out such that
//    M_out * (U, V)^T = (U', V')^T,
// where U', V' are consecutive polynomials in the Euclidean remainder
// sequence of U, V, and V' is the polynomial of highest degree
// satisfying deg(V') <= deg(U) - d_red.

void XHalfGCD(zz_pXMatrix& M_out, zz_pX& U, zz_pX& V, long d_red);

// same as above, except that U is replaced by U', and V by V'


/*************************************************************

             Modular Arithmetic without pre-conditioning

**************************************************************/

// arithmetic mod f.
// all inputs and outputs are polynomials of degree less than deg(f).
// ASSUMPTION: f is assumed monic, and deg(f) > 0.
// NOTE: if you want to do many computations with a fixed f,
//       use the zz_pXModulus data structure and associated routines below.



void MulMod(zz_pX& x, const zz_pX& a, const zz_pX& b, const zz_pX& f);
// x = (a * b) % f

inline zz_pX MulMod(const zz_pX& a, const zz_pX& b, const zz_pX& f)
   { zz_pX x; MulMod(x, a, b, f); NTL_OPT_RETURN(zz_pX, x); }

void SqrMod(zz_pX& x, const zz_pX& a, const zz_pX& f);
// x = a^2 % f

inline zz_pX SqrMod(const zz_pX& a, const zz_pX& f)
   { zz_pX x; SqrMod(x, a, f); NTL_OPT_RETURN(zz_pX, x); }

void MulByXMod(zz_pX& x, const zz_pX& a, const zz_pX& f);
// x = (a * X) mod f

inline zz_pX MulByXMod(const zz_pX& a, const zz_pX& f)
   { zz_pX x; MulByXMod(x, a, f); NTL_OPT_RETURN(zz_pX, x); }

void InvMod(zz_pX& x, const zz_pX& a, const zz_pX& f);
// x = a^{-1} % f, error is a is not invertible

inline zz_pX InvMod(const zz_pX& a, const zz_pX& f)
   { zz_pX x; InvMod(x, a, f); NTL_OPT_RETURN(zz_pX, x); }

long InvModStatus(zz_pX& x, const zz_pX& a, const zz_pX& f);
// if (a, f) = 1, returns 0 and sets x = a^{-1} % f
// otherwise, returns 1 and sets x = (a, f)



/******************************************************************

        Modular Arithmetic with Pre-conditioning

*******************************************************************/


// If you need to do a lot of arithmetic modulo a fixed f,
// build zz_pXModulus F for f.  This pre-computes information about f
// that speeds up the computation a great deal.


class zz_pXModulus {
public:
   zz_pXModulus() : UseFFT(0), n(-1)  { }
   ~zz_pXModulus() { }

   zz_pX f;   // the modulus
   long UseFFT;// flag indicating whether FFT should be used.
   long n;     // n = deg(f)
   long k;     // least k s/t 2^k >= n
   long l;     // least l s/t 2^l >= 2n-3
   fftRep FRep; // 2^k point rep of f
                // H = rev((rev(f))^{-1} rem X^{n-1})
   fftRep HRep; // 2^l point rep of H
   vec_zz_p tracevec;  // mutable

   zz_pXModulus(const zz_pX& ff);

   operator const zz_pX& () const { return f; }
   const zz_pX& val() const { return f; }

};


inline long deg(const zz_pXModulus& F) { return F.n; }

void build(zz_pXModulus& F, const zz_pX& f);
// deg(f) > 0


void rem21(zz_pX& x, const zz_pX& a, const zz_pXModulus& F);
// x = a % f
// deg(a) <= 2(n-1), where n = F.n = deg(f)

void rem(zz_pX& x, const zz_pX& a, const zz_pXModulus& F);