lzz_px.h

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

// x = a % f, no restrictions on deg(a);  makes repeated calls to rem21

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

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


void DivRem(zz_pX& q, zz_pX& r, const zz_pX& a, const zz_pXModulus& F);

void div(zz_pX& q, const zz_pX& a, const zz_pXModulus& F);

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

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



void MulMod(zz_pX& x, const zz_pX& a, const zz_pX& b, const zz_pXModulus& F);
// x = (a * b) % f
// deg(a), deg(b) < n

inline zz_pX MulMod(const zz_pX& a, const zz_pX& b, const zz_pXModulus& 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_pXModulus& F);
// x = a^2 % f
// deg(a) < n


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


void PowerMod(zz_pX& x, const zz_pX& a, const ZZ& e, const zz_pXModulus& F);
// x = a^e % f, e >= 0

inline zz_pX PowerMod(const zz_pX& a, const ZZ& e, const zz_pXModulus& F)
   { zz_pX x; PowerMod(x, a, e, F); NTL_OPT_RETURN(zz_pX, x); }

inline void PowerMod(zz_pX& x, const zz_pX& a, long e, const zz_pXModulus& F)
   { PowerMod(x, a, ZZ_expo(e), F); }

inline zz_pX PowerMod(const zz_pX& a, long e, const zz_pXModulus& F)
   { zz_pX x; PowerMod(x, a, e, F); NTL_OPT_RETURN(zz_pX, x); }



void PowerXMod(zz_pX& x, const ZZ& e, const zz_pXModulus& F);
// x = X^e % f, e >= 0

inline zz_pX PowerXMod(const ZZ& e, const zz_pXModulus& F)
   { zz_pX x; PowerXMod(x, e, F); NTL_OPT_RETURN(zz_pX, x); }

inline void PowerXMod(zz_pX& x, long e, const zz_pXModulus& F)
   { PowerXMod(x, ZZ_expo(e), F); }

inline zz_pX PowerXMod(long e, const zz_pXModulus& F)
   { zz_pX x; PowerXMod(x, e, F); NTL_OPT_RETURN(zz_pX, x); }

void PowerXPlusAMod(zz_pX& x, zz_p a, const ZZ& e, const zz_pXModulus& F);
// x = (X + a)^e % f, e >= 0

inline zz_pX PowerXPlusAMod(zz_p a, const ZZ& e, const zz_pXModulus& F)
   { zz_pX x; PowerXPlusAMod(x, a, e, F); NTL_OPT_RETURN(zz_pX, x); }

inline void PowerXPlusAMod(zz_pX& x, zz_p a, long e, const zz_pXModulus& F)
   { PowerXPlusAMod(x, a, ZZ_expo(e), F); }


inline zz_pX PowerXPlusAMod(zz_p a, long e, const zz_pXModulus& F)
   { zz_pX x; PowerXPlusAMod(x, a, e, F); NTL_OPT_RETURN(zz_pX, x); }

// If you need to compute a * b % f for a fixed b, but for many a's
// (for example, computing powers of b modulo f), it is
// much more efficient to first build a zz_pXMultiplier B for b,
// and then use the routine below.

class zz_pXMultiplier {
public:
   zz_pXMultiplier() : UseFFT(0)  { }
   zz_pXMultiplier(const zz_pX& b, const zz_pXModulus& F);

   ~zz_pXMultiplier() { }

   zz_pX b;   
   long UseFFT;
   fftRep B1; 
   fftRep B2; 

   const zz_pX& val() const { return b; }

};

void build(zz_pXMultiplier& B, const zz_pX& b, const zz_pXModulus& F);



void MulMod(zz_pX& x, const zz_pX& a, const zz_pXMultiplier& B,
                                      const zz_pXModulus& F);

// x = (a * b) % f

inline zz_pX MulMod(const zz_pX& a, const zz_pXMultiplier& B,
                                          const zz_pXModulus& F)
   { zz_pX x; MulMod(x, a, B, F); NTL_OPT_RETURN(zz_pX, x); }




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

              Evaluation and related problems

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


void BuildFromRoots(zz_pX& x, const vec_zz_p& a);
// computes the polynomial (X-a[0]) ... (X-a[n-1]), where n = a.length()

inline zz_pX BuildFromRoots(const vec_zz_p& a)
   { zz_pX x; BuildFromRoots(x, a); NTL_OPT_RETURN(zz_pX, x); }



void eval(zz_p& b, const zz_pX& f, zz_p a);
// b = f(a)


inline zz_p eval(const zz_pX& f, zz_p a)
   { zz_p x; eval(x, f, a); return x; }


void eval(vec_zz_p& b, const zz_pX& f, const vec_zz_p& a);
//  b[i] = f(a[i])

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


void interpolate(zz_pX& f, const vec_zz_p& a, const vec_zz_p& b);
// computes f such that f(a[i]) = b[i]

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



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

                       vectors of zz_pX's

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

NTL_vector_decl(zz_pX,vec_zz_pX)

NTL_eq_vector_decl(zz_pX,vec_zz_pX)

NTL_io_vector_decl(zz_pX,vec_zz_pX)



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

         Modular Composition and Minimal Polynomials

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


// algorithms for computing g(h) mod f




void CompMod(zz_pX& x, const zz_pX& g, const zz_pX& h, const zz_pXModulus& F);
// x = g(h) mod f

inline zz_pX CompMod(const zz_pX& g, const zz_pX& h,
                           const zz_pXModulus& F)
   { zz_pX x; CompMod(x, g, h, F); NTL_OPT_RETURN(zz_pX, x); }


void Comp2Mod(zz_pX& x1, zz_pX& x2, const zz_pX& g1, const zz_pX& g2,
              const zz_pX& h, const zz_pXModulus& F);
// xi = gi(h) mod f (i=1,2)

void Comp3Mod(zz_pX& x1, zz_pX& x2, zz_pX& x3, 
              const zz_pX& g1, const zz_pX& g2, const zz_pX& g3,
              const zz_pX& h, const zz_pXModulus& F);
// xi = gi(h) mod f (i=1..3)



// The routine build (see below) which is implicitly called
// by the various compose and UpdateMap routines builds a table
// of polynomials.  
// If zz_pXArgBound > 0, then the table is limited in
// size to approximamtely that many KB.
// If zz_pXArgBound <= 0, then it is ignored, and space is allocated
// so as to maximize speed.
// Initially, zz_pXArgBound = 0.


// If a single h is going to be used with many g's
// then you should build a zz_pXArgument for h,
// and then use the compose routine below.
// build computes and stores h, h^2, ..., h^m mod f.
// After this pre-computation, composing a polynomial of degree 
// roughly n with h takes n/m multiplies mod f, plus n^2
// scalar multiplies.
// Thus, increasing m increases the space requirement and the pre-computation
// time, but reduces the composition time.
// If zz_pXArgBound > 0, a table of size less than m may be built.

struct zz_pXArgument {
   vec_zz_pX H;
};

extern long zz_pXArgBound;


void build(zz_pXArgument& H, const zz_pX& h, const zz_pXModulus& F, long m);

// m must be > 0, otherwise an error is raised

void CompMod(zz_pX& x, const zz_pX& g, const zz_pXArgument& H, 
             const zz_pXModulus& F);

inline zz_pX
CompMod(const zz_pX& g, const zz_pXArgument& H, const zz_pXModulus& F)
   { zz_pX x; CompMod(x, g, H, F); NTL_OPT_RETURN(zz_pX, x); }



#ifndef NTL_TRANSITION

void UpdateMap(vec_zz_p& x, const vec_zz_p& a,
               const zz_pXMultiplier& B, const zz_pXModulus& F);

inline vec_zz_p
UpdateMap(const vec_zz_p& a,
          const zz_pXMultiplier& B, const zz_pXModulus& F)
   { vec_zz_p x; UpdateMap(x, a, B, F);
     NTL_OPT_RETURN(vec_zz_p, x); }

#endif


/* computes (a, b), (a, (b*X)%f), ..., (a, (b*X^{n-1})%f),
   where ( , ) denotes the vector inner product.

   This is really a "transposed" MulMod by B.
*/

void PlainUpdateMap(vec_zz_p& x, const vec_zz_p& a,
                    const zz_pX& b, const zz_pX& f);


// same as above, but uses only classical arithmetic


void ProjectPowers(vec_zz_p& x, const vec_zz_p& a, long k,
                   const zz_pX& h, const zz_pXModulus& F);

// computes (a, 1), (a, h), ..., (a, h^{k-1} % f)
// this is really a "transposed" compose.

inline vec_zz_p ProjectPowers(const vec_zz_p& a, long k,
                   const zz_pX& h, const zz_pXModulus& F)
{
   vec_zz_p x;
   ProjectPowers(x, a, k, h, F);
   NTL_OPT_RETURN(vec_zz_p, x);
}


void ProjectPowers(vec_zz_p& x, const vec_zz_p& a, long k,
                   const zz_pXArgument& H, const zz_pXModulus& F);

inline vec_zz_p ProjectPowers(const vec_zz_p& a, long k,
                   const zz_pXArgument& H, const zz_pXModulus& F)
{
   vec_zz_p x;
   ProjectPowers(x, a, k, H, F);
   NTL_OPT_RETURN(vec_zz_p, x);
}

// same as above, but uses a pre-computed zz_pXArgument

inline void project(zz_p& x, const vec_zz_p& a, const zz_pX& b)
   { InnerProduct(x, a, b.rep); }

inline zz_p project(const vec_zz_p& a, const zz_pX& b)
   {  zz_p x; project(x, a, b); return x; }


void MinPolySeq(zz_pX& h, const vec_zz_p& a, long m);
// computes the minimum polynomial of a linealy generated sequence;
// m is a bound on the degree of the polynomial;
// required: a.length() >= 2*m

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

void ProbMinPolyMod(zz_pX& h, const zz_pX& g, const zz_pXModulus& F, long m);

inline zz_pX ProbMinPolyMod(const zz_pX& g, const zz_pXModulus& F, long m)
   { zz_pX x; ProbMinPolyMod(x, g, F, m); NTL_OPT_RETURN(zz_pX, x); }


inline void ProbMinPolyMod(zz_pX& h, const zz_pX& g, const zz_pXModulus& F)
   { ProbMinPolyMod(h, g, F, F.n); }

inline zz_pX ProbMinPolyMod(const zz_pX& g, const zz_pXModulus& F)
   { zz_pX x; ProbMinPolyMod(x, g, F); NTL_OPT_RETURN(zz_pX, x); }


// computes the monic minimal polynomial if (g mod f). 
// m = a bound on the degree of the minimal polynomial.
// If this argument is not supplied, it defaults to deg(f).
// The algorithm is probabilistic, always returns a divisor of
// the minimal polynomial, and returns a proper divisor with
// probability at most m/p.

void MinPolyMod(zz_pX& h, const zz_pX& g, const zz_pXModulus& F, long m);

inline zz_pX MinPolyMod(const zz_pX& g, const zz_pXModulus& F, long m)
   { zz_pX x; MinPolyMod(x, g, F, m); NTL_OPT_RETURN(zz_pX, x); }


inline void MinPolyMod(zz_pX& h, const zz_pX& g, const zz_pXModulus& F)
   { MinPolyMod(h, g, F, F.n); }

inline zz_pX MinPolyMod(const zz_pX& g, const zz_pXModulus& F)
   { zz_pX x; MinPolyMod(x, g, F); NTL_OPT_RETURN(zz_pX, x); }


// same as above, but guarantees that result is correct

void IrredPolyMod(zz_pX& h, const zz_pX& g, const zz_pXModulus& F, long m);

inline zz_pX IrredPolyMod(const zz_pX& g, const zz_pXModulus& F, long m)
   { zz_pX x; IrredPolyMod(x, g, F, m); NTL_OPT_RETURN(zz_pX, x); }


inline void IrredPolyMod(zz_pX& h, const zz_pX& g, const zz_pXModulus& F)
   { IrredPolyMod(h, g, F, F.n); }

inline zz_pX IrredPolyMod(const zz_pX& g, const zz_pXModulus& F)
   { zz_pX x; IrredPolyMod(x, g, F); NTL_OPT_RETURN(zz_pX, x); }


// same as above, but assumes that f is irreducible, 
// or at least that the minimal poly of g is itself irreducible.
// The algorithm is deterministic (and is always correct).



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

                   Traces, norms, resultants

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

void TraceVec(vec_zz_p& S, const zz_pX& f);

inline vec_zz_p TraceVec(const zz_pX& f)
   { vec_zz_p x; TraceVec(x, f); NTL_OPT_RETURN(vec_zz_p, x); }


void FastTraceVec(vec_zz_p& S, const zz_pX& f);
void PlainTraceVec(vec_zz_p& S, const zz_pX& f);

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

inline zz_p TraceMod(const zz_pX& a, const zz_pXModulus& F)
   { zz_p x; TraceMod(x, a, F); return x; }

void TraceMod(zz_p& x, const zz_pX& a, const zz_pX& f);

inline zz_p TraceMod(const zz_pX& a, const zz_pX& f)
   { zz_p x; TraceMod(x, a, f); return x; }




void ComputeTraceVec(const zz_pXModulus& F);

void NormMod(zz_p& x, const zz_pX& a, const zz_pX& f);


inline zz_p NormMod(const zz_pX& a, const zz_pX& f)
   { zz_p x; NormMod(x, a, f); return x; }

void resultant(zz_p& rres, const zz_pX& a, const zz_pX& b);

inline zz_p resultant(const zz_pX& a, const zz_pX& b)
   { zz_p x; resultant(x, a, b); return x; }

void CharPolyMod(zz_pX& g, const zz_pX& a, const zz_pX& f);
// g = char poly of (a mod f)
// only implemented for p >= deg(f)+1

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


NTL_CLOSE_NNS

#endif