[Previous] [Up] [Next]

A Tour of NTL: Examples: Modular Arithmetic

NTL also supports modular integer arithmetic. The class ZZ_p represents the integers mod p. Despite the notation, p need not in general be prime, except in situations where this is mathematically required. The classes Vec<ZZ_p> (a.k.a., vec_ZZ_p), Mat<ZZ_p> (a.k.a., mat_ZZ_p), and ZZ_pX represent vectors, matrices, and polynomials mod p, and work much the same way as the corresponding classes for ZZ.

Here is a program that reads a prime number p, and a polynomial f modulo p, and factors it.

#include <NTL/ZZ_pXFactoring.h>

using namespace std;
using namespace NTL;

int main()
{
   ZZ p;
   cin >> p;
   ZZ_p::init(p);

   ZZ_pX f;
   cin >> f;

   Vec< Pair< ZZ_pX, long > > factors;

   CanZass(factors, f);  // calls "Cantor/Zassenhaus" algorithm

   cout << factors << "\n";
    
}

As a program is running, NTL keeps track of a "current modulus" for the class ZZ_p, which can be initialized or changed using ZZ_p::init. This must be done before any variables are declared or computations are done that depend on this modulus.

Please note that for efficiency reasons, NTL does not make any attempt to ensure that variables declared under one modulus are not used under a different one. If that happens, the behavior of a program is completely unpredictable.

Here are two more examples that illustrate the ZZ_p-related classes. The first is a vector addition routine (already supplied by NTL):

#include <NTL/ZZ_p.h>

using namespace std;
using namespace NTL;

void add(Vec<ZZ_p>& x, const Vec<ZZ_p>& a, const Vec<ZZ_p>& b)
{
   long n = a.length();
   if (b.length() != n) Error("vector add: dimension mismatch");

   x.SetLength(n);
   long i;
   for (i = 0; i < n; i++)
      add(x[i], a[i], b[i]);
}

The second example is an inner product routine (also supplied by NTL):

#include <NTL/ZZ_p.h>

using namespace std;
using namespace NTL;

void InnerProduct(ZZ_p& x, const Vec<ZZ_p>& a, const Vec<ZZ_p>& b)
{
   long n = min(a.length(), b.length());
   long i;
   ZZ accum, t;

   accum = 0;
   for (i = 0; i < n; i++) {
      mul(t, rep(a[i]), rep(b[i]));
      add(accum, accum, t);
   }

   conv(x, accum);
}

This second example illustrates two things. First, it illustrates the use of function rep which returns a read-only reference to the representation of a ZZ_p as a ZZ between 0 and p-1. Second, it illustrates a useful algorithmic technique, whereby one computes over ZZ, reducing mod p only when necessary. This reduces the number of divisions that need to be performed significantly, leading to much faster execution.

The class ZZ_p supports all the basic arithmetic operations in both operator and procedural form. All of the basic operations support a "promotion logic", promoting long to ZZ_p.

Note that the class ZZ_p is mainly useful only when you want to work with vectors, matrices, or polynomials mod p. If you just want to do some simple modular arithemtic, it is probably easier to just work with ZZ's directly. This is especially true if you want to work with many different moduli: modulus switching is supported, but it is a bit awkward.

The class ZZ_pX supports all the basic arithmetic operations in both operator and procedural form. All of the basic operations support a "promotion logic", promoting both long and ZZ_p to ZZ_pX.

See ZZ_p.txt for details on ZZ_p; see ZZ_pX.txt for details on ZZ_pX; see ZZ_pXFactoring.txt for details on the routines for factoring polynomials over ZZ_p; see vec_ZZ_p.txt for details on mathematical operations on Vec<ZZ_p>'s; see mat_ZZ_p.txt for details on mathematical operations on Mat<ZZ_p>'s.

There is a mechanism for saving and restoring a modulus, which the following example illustrates. This routine takes as input an integer polynomial and a prime, and tests if the polynomial is irreducible modulo the prime.

#include <NTL/ZZX.h>
#include <NTL/ZZ_pXFactoring.h>

using namespace std;
using namespace NTL;

long IrredTestMod(const ZZX& f, const ZZ& p)
{
   ZZ_pPush push(p); // save current modulus and install p
                     // as current modulus

   return DetIrredTest(conv<ZZ_pX>(f));

   // old modulus is restored automatically when push is destroyed
   // upon return
}

The modulus switching mechanism is actually quite a bit more general and flexible than this example illustrates.

Note the use of the conversion function conv<ZZ_pX>. We could also have used the equivalent procedural form:

   ZZ_pX f1;
   conv(f1, f);
   return DetIrredTest(f1);

Suppose in the above example that p is known in advance to be a small, single-precision prime. In this case, NTL provides a class zz_p, that acts just like ZZ_p, along with corresponding classes Vec<zz_p>, Mat<zz_p>, and zz_pX. The interfaces to all of the routines are generally identical to those for ZZ_p. However, the routines are much more efficient, in both time and space.

For small primes, the routine in the previous example could be coded as follows.

#include <NTL/ZZX.h>
#include <NTL/lzz_pXFactoring.h>

using namespace std;
using namespace NTL;

long IrredTestMod(const ZZX& f, long p)
{
   zz_pPush push(p);
   return DetIrredTest(conv<zz_pX>(f));
}

The following is a routine (essentially the same as implemented in NTL) for computing the GCD of polynomials with integer coefficients. It uses a "modular" approach: the GCDs are computed modulo small primes, and the results are combined using the Chinese Remainder Theorem (CRT). The small primes are specially chosen "FFT primes", which are of a special form that allows for particular fast polynomial arithmetic.

#include <NTL/ZZX.h>

using namespace std;
using namespace NTL;

void GCD(ZZX& d, const ZZX& a, const ZZX& b)
{
   if (a == 0) {
      d = b;
      if (LeadCoeff(d) < 0) negate(d, d);
      return;
   }

   if (b == 0) {
      d = a;
      if (LeadCoeff(d) < 0) negate(d, d);
      return;
   }

   ZZ c1, c2, c;
   ZZX f1, f2;

   content(c1, a);
   divide(f1, a, c1);

   content(c2, b);
   divide(f2, b, c2);

   GCD(c, c1, c2);

   ZZ ld;
   GCD(ld, LeadCoeff(f1), LeadCoeff(f2));

   ZZX g, res;

   ZZ prod;

   zz_pPush push; // save current modulus, restore upon return

   long FirstTime = 1;

   long i;
   for (i = 0; ;i++) {
      zz_p::FFTInit(i);
      long p = zz_p::modulus();

      if (divide(LeadCoeff(f1), p) || divide(LeadCoeff(f2), p)) continue;

      zz_pX G, F1, F2;
      zz_p  LD;

      conv(F1, f1);
      conv(F2, f2);
      conv(LD, ld);

      GCD(G, F1, F2);
      mul(G, G, LD);


      if (deg(G) == 0) {
         res = 1;
         break;
      }

      if (FirstTime || deg(G) < deg(g)) {
         prod = 1;
         g = 0;
         FirstTime = 0;
      }
      else if (deg(G) > deg(g)) {
         continue;
      }

      if (!CRT(g, prod, G)) {
         PrimitivePart(res, g);
         if (divide(f1, res) && divide(f2, res))
            break;
      }

   }

   mul(d, res, c);
   if (LeadCoeff(d) < 0) negate(d, d);
}

See lzz_p.txt for details on zz_p; see lzz_pX.txt for details on zz_pX; see lzz_pXFactoring.txt for details on the routines for factoring polynomials over zz_p; see vec_lzz_p.txt for details on vec_zz_p; see mat_lzz_p.txt for details on mat_zz_p.

NTL provides a number of "residue class" types with a dynamic modulus stored as a global variable: the types ZZ_p and zz_p, discussed above, as well as the types ZZ_pE, zz_pE, and GF2E, discussed later.

Some caution must be used so that a variable constructed under one modulus is not used "out of context", when a different modulus, or perhaps no modulus, is installed as the current modulus. While arithmetic operations should certainly be avoided, NTL does take care to allow for certain operations to be safely performed "out of context". These operations include default and copy constructors, as well as assignment.

Arithmetic mod 2 is such an important special case that NTL provides a class GF2, that acts just like ZZ_p when p == 2, along with corresponding classes Vec<GF2>, Mat<GF2>, and GF2X. The interfaces to all of the routines are generally identical to those for ZZ_p. However, the routines are much more efficient, in both time and space. Note that Vec<GF2> is an explicit specialization of the template class Vec<T>, with a special implementation that packs the coefficients into the bits of a machine word. You need to include the header file <NTL/vec_GF2.h> to use the class Vec<GF2>.

This example illustrates the GF2X and Mat<GF2> classes with a simple routine to test if a polynomial over GF(2) is irreducible using linear algebra. NTL's built-in irreducibility test is to be preferred, however.

#include <NTL/GF2X.h>
#include <NTL/mat_GF2.h>

using namespace std;
using namespace NTL;

long MatIrredTest(const GF2X& f)
{
   long n = deg(f);

   if (n <= 0) return 0;
   if (n == 1) return 1;

   if (GCD(f, diff(f)) != 1) return 0;

   Mat<GF2> M;

   M.SetDims(n, n);

   GF2X x_squared = GF2X(INIT_MONO, 2);

   GF2X g;
   g = 1;

   for (long i = 0; i < n; i++) {
      VectorCopy(M[i], g, n);
      M[i][i] += 1;
      g = (g * x_squared) % f;
   }

   long rank = gauss(M);

   if (rank == n-1)
      return 1;
   else
      return 0;
}

Note that the statement

   g = (g * x_squared) % f;

could be replace d by the more efficient code sequence

   MulByXMod(g, g, f);
   MulByXMod(g, g, f);

but this would not significantly impact the overall running time, since it is the Gaussian elimination that dominates the running time.

See GF2.txt for details on GF2; see GF2X.txt for details on GF2X; see GF2XFactoring.txt for details on the routines for factoring polynomials over GF2; see vec_GF2.txt for details on vec_GF2; see mat_GF2.txt for details on mat_GF2.

[Previous] [Up] [Next]