zz.h

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

inline ZZ SubMod(const ZZ& a, long b, const ZZ& n)
   { ZZ x; SubMod(x, a, b, n); NTL_OPT_RETURN(ZZ, x); }

void SubMod(ZZ& x, long a, const ZZ& b, const ZZ& n);
inline ZZ SubMod(long a, const ZZ& b, const ZZ& n)
   { ZZ x; SubMod(x, a, b, n); NTL_OPT_RETURN(ZZ, x); }

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

   { NTL_zmulmod(a.rep, b.rep, n.rep, &x.rep); }

inline ZZ MulMod(const ZZ& a, const ZZ& b, const ZZ& n)
   { ZZ x; MulMod(x, a, b, n); NTL_OPT_RETURN(ZZ, x); }

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

   { NTL_zsmulmod(a.rep, b, n.rep, &x.rep); }

inline ZZ MulMod(const ZZ& a, long b, const ZZ& n)
   { ZZ x; MulMod(x, a, b, n); NTL_OPT_RETURN(ZZ, x); }

inline void MulMod(ZZ& x, long a, const ZZ& b, const ZZ& n)
   { MulMod(x, b, a, n); }

inline ZZ MulMod(long a, const ZZ& b, const ZZ& n)
   { ZZ x; MulMod(x, a, b, n); NTL_OPT_RETURN(ZZ, x); }


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

   { NTL_zsqmod(a.rep, n.rep, &x.rep); }

inline ZZ SqrMod(const ZZ& a, const ZZ& n)
   {  ZZ x; SqrMod(x, a, n); NTL_OPT_RETURN(ZZ, x); }

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

   { NTL_zinvmod(a.rep, n.rep, &x.rep); }

inline ZZ InvMod(const ZZ& a, const ZZ& n)
   {  ZZ x; InvMod(x, a, n); NTL_OPT_RETURN(ZZ, x); }


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

  { return NTL_zinv(a.rep, n.rep, &x.rep); }


inline void PowerMod(ZZ& x, const ZZ& a, const ZZ& e, const ZZ& n)
   { NTL_zpowermod(a.rep, e.rep, n.rep, &x.rep); }

inline ZZ PowerMod(const ZZ& a, const ZZ& e, const ZZ& n)
   { ZZ x; PowerMod(x, a, e, n); NTL_OPT_RETURN(ZZ, x); }

inline void PowerMod(ZZ& x, const ZZ& a, long e, const ZZ& n)
   { PowerMod(x, a, ZZ_expo(e), n); }

inline ZZ PowerMod(const ZZ& a, long e, const ZZ& n)
   { ZZ x; PowerMod(x, a, e, n); NTL_OPT_RETURN(ZZ, x); }






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

   Jacobi symbol and modular squre 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);
//  computes square root of a mod n;
//  assumes n is an odd prime, and that a is a square mod n

inline ZZ SqrRootMod(const ZZ& a, const ZZ& n)
   { ZZ x; SqrRootMod(x, a, n); NTL_OPT_RETURN(ZZ, x); }




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


                    Small Prime Generation


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


// primes are generated in sequence, starting at 2, 
// and up until (2*NTL_PRIME_BND+1)^2, which is less than NTL_SP_BOUND.

#if (NTL_SP_NBITS > 30)
#define NTL_PRIME_BND ((1L << 14) - 1)
#else
#define NTL_PRIME_BND ((1L << (NTL_SP_NBITS/2-1)) - 1)
#endif


class PrimeSeq {


char *movesieve;
char *movesieve_mem;
long pindex;
long pshift;
long exhausted;

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

// auxilliary routines

void start();
void shift(long);

};




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

                      Input/Output

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

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




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

    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 1..2^{NTL_SP_NBITS}-1.
*/



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

{
   long res = a + b;
#if (NTL_ARITH_RIGHT_SHIFT && defined(NTL_AVOID_BRANCHING))
   res -= n;
   res += (res >> (NTL_BITS_PER_LONG-1)) & n;
   return res;
#else
   if (res >= n)
      return res - n;
   else
      return res;
#endif
}

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

{
   long res = a - b;
#if (NTL_ARITH_RIGHT_SHIFT && defined(NTL_AVOID_BRANCHING))
   res += (res >> (NTL_BITS_PER_LONG-1)) & n;
   return res;
#else
   if (res < 0)
      return res + n;
   else
      return res;
#endif
}

inline long NegateMod(long a, long n)
{
   return SubMod(0, a, n);
}


#if (defined(NTL_SINGLE_MUL))


#if (!defined(NTL_FAST_INT_MUL))


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

{
   double ab;
   long q, res;

   ab = ((double) a) * ((double) b);
   q  = (long) (ab/((double) n));  // q could be off by (+/-) 1
   res = (long) (ab - ((double) q)*((double) n));
#if (NTL_ARITH_RIGHT_SHIFT && defined(NTL_AVOID_BRANCHING))
   res += (res >> (NTL_BITS_PER_LONG-1)) & n;
   res -= n;
   res += (res >> (NTL_BITS_PER_LONG-1)) & n;
#else
   if (res >= n)
      res -= n;
   else if (res < 0)
      res += n;
#endif
   return res;
}

/*
The following MulMod takes a fourth argument, ninv,
which is assumed to equal 1/((double) n).
It is usually faster than the above.
*/

inline long MulMod(long a, long b, long n, double ninv)
{
   double ab;
   long q, res;

   ab = ((double) a) * ((double) b);
   q  = (long) (ab*ninv);   // q could be off by (+/-) 1
   res = (long) (ab - ((double) q)*((double) n));
#if (NTL_ARITH_RIGHT_SHIFT && defined(NTL_AVOID_BRANCHING))
   res += (res >> (NTL_BITS_PER_LONG-1)) & n;
   res -= n;
   res += (res >> (NTL_BITS_PER_LONG-1)) & n;
#else
   if (res >= n)
      res -= n;
   else if (res < 0)
      res += n;
#endif
   return res;
}

/*
Yet another MulMod.
This time, the 4th argument should be ((double) b)/((double) n).
*/

inline long MulMod2(long a, long b, long n, double bninv)
{
   double ab;
   long q, res;

   ab = ((double) a)*((double) b);
   q = (long) (((double) a)*bninv);
   res = (long) (ab - ((double) q)*((double) n));
#if (NTL_ARITH_RIGHT_SHIFT && defined(NTL_AVOID_BRANCHING))
   res += (res >> (NTL_BITS_PER_LONG-1)) & n;
   res -= n;
   res += (res >> (NTL_BITS_PER_LONG-1)) & n;
#else
   if (res >= n)
      res -= n;
   else if (res < 0)
      res += n;
#endif
   return res;
}

inline long MulDivRem(long& qq, long a, long b, long n, double bninv)
{
   double ab;
   long q, res;

   ab = ((double) a)*((double) b);
   q = (long) (((double) a)*bninv);
   res = (long) (ab - ((double) q)*((double) n));
   if (res >= n) {
      res -= n;
      q++;
   } else if (res < 0) {
      res += n;
      q--;
   }

   qq = q;
   return res;
}

#else

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

{
   double ab, xx;
   long iab, q, res;

   ab = ((double) a) * ((double) b);
   q  = (long) (ab/((double) n));  // q could be off by (+/-) 1

   xx = ab + 4503599627370496.0;
   NTL_FetchLo(iab, xx);

   res = iab - q*n;

#if (NTL_ARITH_RIGHT_SHIFT && defined(NTL_AVOID_BRANCHING))
   res += (res >> (NTL_BITS_PER_LONG-1)) & n;
   res -= n;
   res += (res >> (NTL_BITS_PER_LONG-1)) & n;
#else
   if (res >= n)
      res -= n;
   else if (res < 0)
      res += n;
#endif
   return res;
}

/*
The following MulMod takes a fourth argument, ninv,
which is assumed to equal 1/((double) n).
It is usually faster than the above.
*/

inline long MulMod(long a, long b, long n, double ninv)
{
   double ab, xx;
   long iab, q, res;

   ab = ((double) a) * ((double) b);
   q  = (long) (ab*ninv);   // q could be off by (+/-) 1

   xx = ab + 4503599627370496.0;
   NTL_FetchLo(iab, xx);

   res = iab - q*n;
#if (NTL_ARITH_RIGHT_SHIFT && defined(NTL_AVOID_BRANCHING))
   res += (res >> (NTL_BITS_PER_LONG-1)) & n;
   res -= n;
   res += (res >> (NTL_BITS_PER_LONG-1)) & n;
#else
   if (res >= n)
      res -= n;
   else if (res < 0)
      res += n;
#endif
   return res;
}

/*
Yet another MulMod.
This time, the 4th argument should be ((double) b)/((double) n).
*/

inline long MulMod2(long a, long b, long n, double bninv)
{
   long q, res;

   q = (long) (((double) a)*bninv);
   res = a*b - q*n;
#if (NTL_ARITH_RIGHT_SHIFT && defined(NTL_AVOID_BRANCHING))
   res += (res >> (NTL_BITS_PER_LONG-1)) & n;
   res -= n;
   res += (res >> (NTL_BITS_PER_LONG-1)) & n;
#else
   if (res >= n)
      res -= n;
   else if (res < 0)
      res += n;
#endif
   return res;
}


inline long MulDivRem(long& qq, long a, long b, long n, double bninv)
{
   long q, res;

   q = (long) (((double) a)*bninv);
   res = a*b - q*n;
   if (res >= n) {
      res -= n;
      q++;
   } else if (res < 0) {
      res += n;
      q--;
   }

   qq = q;
   return res;
}

#endif




#else

inline long MulMod(long a, long b, long n)
{
   long q, res;

   q  = (long) ((((double) a) * ((double) b)) / ((double) n)); 
   res = a*b - q*n;
#if (NTL_ARITH_RIGHT_SHIFT && defined(NTL_AVOID_BRANCHING))
   res += (res >> (NTL_BITS_PER_LONG-1)) & n;
   res -= n;
   res += (res >> (NTL_BITS_PER_LONG-1)) & n;
#else
   if (res >= n)
      res -= n;
   else if (res < 0)
      res += n;
#endif
   return res;
}

inline long MulMod(long a, long b, long n, double ninv)
{
   long q, res;

   q  = (long) ((((double) a) * ((double) b)) * ninv); 
   res = a*b - q*n;
#if (NTL_ARITH_RIGHT_SHIFT && defined(NTL_AVOID_BRANCHING))
   res += (res >> (NTL_BITS_PER_LONG-1)) & n;
   res -= n;
   res += (res >> (NTL_BITS_PER_LONG-1)) & n;
#else
   if (res >= n)
      res -= n;
   else if (res < 0)
      res += n;
#endif
   return res;
}


inline long MulMod2(long a, long b, long n, double bninv)
{
   long q, res;

   q  = (long) (((double) a) * bninv);
   res = a*b - q*n;
#if (NTL_ARITH_RIGHT_SHIFT && defined(NTL_AVOID_BRANCHING))
   res += (res >> (NTL_BITS_PER_LONG-1)) & n;
   res -= n;
   res += (res >> (NTL_BITS_PER_LONG-1)) & n;
#else
   if (res >= n)
      res -= n;
   else if (res < 0)
      res += n;
#endif
   return res;
}

inline long MulDivRem(long& qq, long a, long b, long n, double bninv)
{
   long q, res;

   q  = (long) (((double) a) * bninv);
   res = a*b - q*n;
   if (res >= n) {
      res -= n;
      q++;
   } else if (res < 0) {
      res += n;
      q--;
   }

   qq = q;
   return res;
}


#endif


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 >= 0


NTL_CLOSE_NNS


#endif