core/vnl/vnl_finite.h
Go to the documentation of this file.
00001 // This is core/vnl/vnl_finite.h
00002 #ifndef vnl_finite_h_
00003 #define vnl_finite_h_
00004 //:
00005 // \file
00006 // \brief modulo-N arithmetic (finite ring Z_N and Z_N[X])
00007 //
00008 // The templated vnl_finite_int<N> provides arithmetic "modulo N", i.e.,
00009 // arithmetic in the finite (Galois) field GF(N) in case N is a prime
00010 // or just in the finite ring (or semi-field) of integers modulo N otherwise.
00011 // In that case division makes no sense (unless no zero divisor is involved),
00012 // but all other operations remain valid.
00013 //
00014 // Note that this does not cover finite fields with non-prime sizes (4,8,9,...).
00015 // These are covered by the vnl_finite_int_poly<N,M> class, which implements
00016 // arithmetic with polynomials of degree < M over vnl_finite_int<N>.
00017 // Multiplication is defined modulo a degree M polynomial.
00018 //
00019 // For N prime, and when the "modulo" polynomial is irreducible,
00020 // vnl_finite_int_poly<N,M> implements the finite field GF(N^M).
00021 //
00022 // \author
00023 //  Peter Vanroose, K.U.Leuven, ESAT/PSI.
00024 // \date 5 May 2002.
00025 //
00026 // \verbatim
00027 // Modifications
00028 //  1 June 2002 - Peter Vanroose - added totient(), decompose(), is_unit(), order(), log(), exp().
00029 //  4 June 2002 - Peter Vanroose - renamed class and file name
00030 //  8 June 2002 - Peter Vanroose - added vnl_finite_int_poly<N,M>
00031 //  16 Dec 2007 - Peter Vanroose - more efficient implementation of Ntothe()
00032 // \endverbatim
00033 
00034 #include <vcl_iostream.h>
00035 #include <vcl_cassert.h>
00036 #include <vcl_vector.h>
00037 #include <vcl_cstddef.h>
00038 
00039 //: finite modulo-N arithmetic
00040 //
00041 // The templated vnl_finite_int<N> provides arithmetic "modulo N", i.e.,
00042 // arithmetic in the finite (Galois) field GF(N) in case N is a prime
00043 // or just in the finite ring (or semi-field) of integers modulo N otherwise.
00044 // In that case division makes no sense (unless no zero divisor is involved),
00045 // but all other operations remain valid.
00046 //
00047 template <int N>
00048 class vnl_finite_int
00049 {
00050   int val_; //!< value of this number (smallest nonnegative representation)
00051 
00052   typedef vnl_finite_int<N> Base;
00053  public:
00054   //: The number of different finite_int numbers of this type
00055   static unsigned int cardinality() { return N; }
00056 
00057   //: Creates a finite int element.
00058   //  Default constructor gives 0.
00059   //  Also serves as automatic cast from int to vnl_finite_int.
00060   inline vnl_finite_int(int x = 0) : val_((x%=N)<0?N+x:x), mo_(0), lp1_(0) {assert(N>1);}
00061   //  Copy constructor
00062   inline vnl_finite_int(Base const& x) : val_(int(x)), mo_(x.mo_), lp1_(x.lp1_) {}
00063   //  Destructor
00064   inline ~vnl_finite_int() {}
00065   // Implicit conversions
00066   inline operator int() const { return val_; }
00067   inline operator int() { return val_; }
00068   inline operator long() const { return val_; }
00069   inline operator long() { return val_; }
00070   inline operator short() const { short r = (short)val_; assert(r == val_); return r; }
00071   inline operator short() { short r = (short)val_; assert(r == val_); return r; }
00072 
00073   //: Assignment
00074   inline Base& operator=(Base const& x) { val_ = int(x); mo_=x.mo_; lp1_=x.lp1_; return *this; }
00075   inline Base& operator=(int x) { x%=N; val_ = x<0 ? N+x : x; mo_=lp1_=0; return *this; }
00076 
00077   //: Comparison of finite int numbers.
00078   // Note that finite ints have no order, so < and > make no sense.
00079   inline bool operator==(Base const& x) const { return val_ == int(x); }
00080   inline bool operator!=(Base const& x) const { return val_ != int(x); }
00081   inline bool operator==(int x) const { return operator==(Base(x)); }
00082   inline bool operator!=(int x) const { return !operator==(x); }
00083 
00084   //: Unary minus - returns the additive inverse
00085   inline Base operator-() const { return Base(N-val_); }
00086   //: Unary plus - returns the current number
00087   inline Base operator+() const { return *this; }
00088   //: Unary not - returns true if finite int number is equal to zero.
00089   inline bool operator!() const { return val_ == 0; }
00090 
00091   //: Plus&assign: replace lhs by lhs + rhs
00092   inline Base& operator+=(Base const& r) { mo_=lp1_=0; val_ += int(r); if (val_ >= int(N)) val_ -= N; return *this; }
00093   inline Base& operator+=(int r) { return operator=(val_+r); }
00094   //: Minus&assign: replace lhs by lhs - rhs
00095   inline Base& operator-=(Base const& r) { mo_=lp1_=0; val_ -= int(r); if (val_ < 0) val_ += N; return *this; }
00096   inline Base& operator-=(int r) { return operator=(val_-r); }
00097   //: Multiply&assign: replace lhs by lhs * rhs
00098   inline Base& operator*=(int r) {
00099     r %=N; if (r<0) r=N+r;
00100     // This rather complicated implementation is necessary to avoid integer overflow
00101     if (N<=0x7fff || (val_<=0x7fff && r<=0x7fff)) { val_ *= r; val_ %= N; return *this; }
00102     else { int v=val_; operator+=(v); operator*=(r/2); if (r%2) operator+=(v); return *this; }
00103   }
00104   //: Multiply&assign: replace lhs by lhs * rhs
00105   inline Base& operator*=(Base const& r) {
00106     mo_=0;
00107     if (lp1_!=0 && r.lp1_!=0) set_log(lp1_+r.lp1_-2); else lp1_=0;
00108     // This rather complicated implementation is necessary to avoid integer overflow
00109     unsigned int s=int(r);
00110     if (N<=0x7fff || (val_<=0x7fff && s<=0x7fff)) { val_ *= s; val_ %= N; return *this; }
00111     else { int v=val_; operator+=(v); operator*=(s/2); if (s%2) operator+=(v); return *this; }
00112   }
00113 
00114   //: Return the Euler totient function, i.e., the number of units of this ring
00115   //  This number also equals the periodicity of the exponent: every unit,
00116   //  when raised to this power, yields 1.
00117   static unsigned int totient() {
00118     static unsigned int t_ = 0; // cached value
00119     if (t_ != 0) return t_;
00120     vcl_vector<unsigned int> d = decompose();
00121     t_ = 1; unsigned int p = 1;
00122     for (unsigned int i=0; i<d.size(); ++i)
00123     {
00124       if (p != d[i]) t_ *= d[i]-1;
00125       else           t_ *= d[i];
00126       p = d[i];
00127     }
00128     return t_;
00129   }
00130 
00131   //: Multiplicative inverse.
00132   //  Uses exp() and log() for efficient computation, unless log() is not defined.
00133   Base reciproc() const {
00134     assert(is_unit());
00135     if (val_==1) return *this;
00136     Base z = smallest_generator();
00137     if (z!=1) return exp(Base::totient()-log());
00138     for (unsigned int r=1; r<=N/2; ++r) {
00139       unsigned int t=int(*this*int(r));
00140       if (t==1) return r; else if (t==N-1) return N-r;
00141     }
00142     return 0; // This will never be executed
00143   }
00144 
00145   //: Divide&assign.  Uses r.reciproc() for efficient computation.
00146   inline Base& operator/=(Base const& r) {
00147     assert(r.is_unit());
00148     return val_ == 0 ? operator=(0) : operator*=(r.reciproc());
00149   }
00150 
00151   //: Pre-increment (++r).
00152   inline Base& operator++() { mo_=lp1_=0; ++val_; if (val_==N) val_=0; return *this; }
00153   //: Pre-decrement (--r).
00154   inline Base& operator--() { mo_=lp1_=0; if (val_==0) val_=N; --val_; return *this; }
00155   //: Post-increment (r++).
00156   inline Base operator++(int){Base b=*this; mo_=lp1_=0; ++val_; if (val_==N) val_=0; return b; }
00157   //: Post-decrement (r--).
00158   inline Base operator--(int){Base b=*this; mo_=lp1_=0; if (val_==0) val_=N; --val_; return b;}
00159 
00160   //: Write N as the unique product of prime factors.
00161   static vcl_vector<unsigned int> decompose() {
00162     static vcl_vector<unsigned int> decomposition_ = vcl_vector<unsigned int>(); // cached value
00163     if (decomposition_.size() > 0) return decomposition_;
00164     unsigned int r = N;
00165     for (unsigned int d=2; d*d<=r; ++d)
00166       while (r%d == 0) { decomposition_.push_back(d); r /= d; }
00167     if (r > 1) decomposition_.push_back(r);
00168     return decomposition_;
00169   }
00170 
00171   //: Return true when N is a prime number, i.e., when this ring is a field
00172   static inline bool is_field() {
00173     vcl_vector<unsigned int> d = Base::decompose();
00174     return d.size() == 1;
00175   }
00176 
00177   //: Return true only when x is a unit in this ring.
00178   //  In a field, all numbers except 0 are units.
00179   //  The total number of units is given by the Euler totient function.
00180   inline bool is_unit() const { return gcd(val_) == 1; }
00181 
00182   //: Return true only when x is a zero divisor, i.e., is not a unit.
00183   inline bool is_zero_divisor() const { return gcd(val_) != 1; }
00184 
00185   //: The additive order of x is the smallest nonnegative r such that r*x == 0.
00186   inline unsigned int additive_order() const { if (val_ == 0) return 1; return N/gcd(val_); }
00187 
00188   //: The multiplicative order of x is the smallest r (>0) such that x^r == 1.
00189   unsigned int multiplicative_order() const {
00190     if (mo_ != 0) return mo_;
00191     if (gcd(val_) != 1) return -1; // should actually return infinity
00192     Base y = val_;
00193     for (int r=1; r<N; ++r) { if (y==1) return mo_=r; y *= val_; }
00194     return 0; // this should not happen
00195   }
00196 
00197   //: Return the smallest multiplicative generator of the units in this ring.
00198   //  This is only possible if the units form a cyclic group for multiplication.
00199   //  If not, smallest_generator() returns 1 to indicate this fact.
00200   //  Note that the multiplication group of a finite \e field is always cyclic.
00201   static Base smallest_generator() {
00202     static Base gen_ = 0; // cached value
00203     if (gen_ != 0) return gen_;
00204     if (N==2) return gen_=1;
00205     unsigned int h = Base::totient() / 2; // note that totient() is always even
00206     for (gen_=2; gen_!=0; ++gen_) {
00207       // calculate gen_^h
00208       unsigned int g=h; Base y = 1, z = gen_; while (g>0) { if (g%2) y *= z; g/=2; z*=z; }
00209       // quick elimination of non-generator:
00210       if (y == 1) continue;
00211       // calculate gen_.multiplicative_order() only if really necessary:
00212       if (gen_.multiplicative_order() == Base::totient()) { gen_.set_log(1); return gen_; }
00213     }
00214     assert(!Base::is_field()); // can only reach this point when N is not prime
00215     return gen_=1;
00216   }
00217 
00218   //: Return the r-th power of this number.
00219   inline Base pow(int r) {
00220     r %= Base::totient(); if (r<0) r += Base::totient();
00221     if (r==0) return 1; if (r==1) return *this;
00222     Base y = 1, z = *this; int s=r; while (s!=0) { if (s%2) y*=z; s/=2; z*=z; }
00223     if (lp1_ != 0) y.set_log(r*(lp1_-1));
00224     return y;
00225   }
00226 
00227   //: Return the smallest nonnegative exponent r for which x=g^r, where g is the smallest generator.
00228   unsigned int log() const {
00229     if (is_zero_divisor()) return -1; // should actually return minus infinity
00230     if (lp1_ != 0) return lp1_-1;
00231     Base z = smallest_generator();
00232     assert(N==2||z!=1); // otherwise, the units of this ring do not form a cyclic group
00233     Base y = 1;
00234     for (lp1_=1; lp1_<=N; ++lp1_) {
00235       if (y == *this) return lp1_-1;
00236       y *= z;
00237     }
00238     return -1; // should never reach this point
00239   }
00240 
00241   //: Return the inverse of log(), i.e., return g^r where g is the smallest generator.
00242   static inline Base exp(int r) {
00243     Base z = smallest_generator();
00244     assert(N==2||z!=1); // otherwise, the units of this ring do not form a cyclic group
00245     return z.pow(r);
00246   }
00247 
00248   //: Calculate the greatest common divisor of l and N.
00249   static inline unsigned int gcd(unsigned int l, unsigned int n) {
00250     unsigned int l1 = n;
00251     while (l!=0) { unsigned int t = l; l = l1 % l; l1 = t; }
00252     return l1;
00253   }
00254   static inline unsigned int gcd(unsigned int l) { return gcd(l, N); }
00255 
00256  private:
00257   //: private function to set cached value of lp1_ when available
00258   void set_log(unsigned int r) const { r %= Base::totient(); lp1_ = r+1; }
00259 
00260   mutable unsigned int mo_; //!< cached value for multiplicative order
00261   mutable unsigned int lp1_; //!< cached value for 1+log()
00262 };
00263 
00264 //: formatted output
00265 // \relatesalso vnl_finite_int
00266 template <int N>
00267 inline vcl_ostream& operator<< (vcl_ostream& s, vnl_finite_int<N> const& r)
00268 {
00269   return s << int(r);
00270 }
00271 
00272 //: simple input
00273 // \relatesalso vnl_finite_int
00274 template <int N>
00275 inline vcl_istream& operator>> (vcl_istream& s, vnl_finite_int<N>& r)
00276 {
00277   int n; s >> n; r=n; return s;
00278 }
00279 
00280 //: Returns the sum of two finite int numbers.
00281 // \relatesalso vnl_finite_int
00282 template <int N>
00283 inline vnl_finite_int<N> operator+ (vnl_finite_int<N> const& r1, vnl_finite_int<N> const& r2)
00284 {
00285   vnl_finite_int<N> result(r1); return result += r2;
00286 }
00287 
00288 //: Returns the sum of two finite int numbers.
00289 // \relatesalso vnl_finite_int
00290 template <int N>
00291 inline vnl_finite_int<N> operator+ (vnl_finite_int<N> const& r1, int r2)
00292 {
00293   vnl_finite_int<N> result(r1); return result += r2;
00294 }
00295 
00296 //: Returns the sum of two finite int numbers.
00297 // \relatesalso vnl_finite_int
00298 template <int N>
00299 inline vnl_finite_int<N> operator+ (int r2, vnl_finite_int<N> const& r1)
00300 {
00301   vnl_finite_int<N> result(r1); return result += r2;
00302 }
00303 
00304 //: Returns the difference of two finite int numbers.
00305 // \relatesalso vnl_finite_int
00306 template <int N>
00307 inline vnl_finite_int<N> operator- (vnl_finite_int<N> const& r1, vnl_finite_int<N> const& r2)
00308 {
00309   vnl_finite_int<N> result(r1); return result -= r2;
00310 }
00311 
00312 //: Returns the difference of two finite int numbers.
00313 // \relatesalso vnl_finite_int
00314 template <int N>
00315 inline vnl_finite_int<N> operator- (vnl_finite_int<N> const& r1, int r2)
00316 {
00317   vnl_finite_int<N> result(r1); return result -= r2;
00318 }
00319 
00320 //: Returns the difference of two finite int numbers.
00321 // \relatesalso vnl_finite_int
00322 template <int N>
00323 inline vnl_finite_int<N> operator- (int r2, vnl_finite_int<N> const& r1)
00324 {
00325   vnl_finite_int<N> result(-r1); return result += r2;
00326 }
00327 
00328 //: Returns the product of two finite int numbers.
00329 // \relatesalso vnl_finite_int
00330 template <int N>
00331 inline vnl_finite_int<N> operator* (vnl_finite_int<N> const& r1, vnl_finite_int<N> const& r2)
00332 {
00333   vnl_finite_int<N> result(r1); return result *= r2;
00334 }
00335 
00336 //: Returns the product of two finite int numbers.
00337 // \relatesalso vnl_finite_int
00338 template <int N>
00339 inline vnl_finite_int<N> operator* (vnl_finite_int<N> const& r1, int r2)
00340 {
00341   vnl_finite_int<N> result(r1); return result *= r2;
00342 }
00343 
00344 //: Returns the product of two finite int numbers.
00345 // \relatesalso vnl_finite_int
00346 template <int N>
00347 inline vnl_finite_int<N> operator* (int r2, vnl_finite_int<N> const& r1)
00348 {
00349   vnl_finite_int<N> result(r1); return result *= r2;
00350 }
00351 
00352 //: Returns the quotient of two finite int numbers.
00353 //  Uses r2.reciproc() for efficient computation.
00354 // \relatesalso vnl_finite_int
00355 template <int N>
00356 inline vnl_finite_int<N> operator/(vnl_finite_int<N> const& r1, vnl_finite_int<N> const& r2)
00357 {
00358   assert(r2.is_unit()); return r1 == 0 ? vnl_finite_int<N>(0) : r1*r2.reciproc();
00359 }
00360 
00361 //: Returns the quotient of two finite int numbers.
00362 // \relatesalso vnl_finite_int
00363 template <int N>
00364 inline vnl_finite_int<N> operator/ (vnl_finite_int<N> const& r1, int r2)
00365 {
00366   vnl_finite_int<N> result(r1); return result /= r2;
00367 }
00368 
00369 //: Returns the quotient of two finite int numbers.
00370 // \relatesalso vnl_finite_int
00371 template <int N>
00372 inline vnl_finite_int<N> operator/ (int r1, vnl_finite_int<N> const& r2)
00373 {
00374   vnl_finite_int<N> result(r1); return result /= r2;
00375 }
00376 
00377 //:
00378 // \relatesalso vnl_finite_int
00379 template <int N>
00380 inline bool operator== (int  r1, vnl_finite_int<N> const& r2) { return r2==r1; }
00381 
00382 //:
00383 // \relatesalso vnl_finite_int
00384 template <int N>
00385 inline bool operator!= (int  r1, vnl_finite_int<N> const& r2) { return r2!=r1; }
00386 
00387 //:
00388 // \relatesalso vnl_finite_int
00389 template <int N>
00390 inline vnl_finite_int<N> vnl_math_squared_magnitude(vnl_finite_int<N> const& x) { return x*x; }
00391 
00392 //:
00393 // \relatesalso vnl_finite_int
00394 template <int N>
00395 inline vnl_finite_int<N> vnl_math_sqr(vnl_finite_int<N> const& x) { return x*x; }
00396 
00397 //:
00398 // \relatesalso vnl_finite_int
00399 template <int N>
00400 inline bool vnl_math_isnan(vnl_finite_int<N> const& ){return false;}
00401 
00402 //:
00403 // \relatesalso vnl_finite_int
00404 template <int N>
00405 inline bool vnl_math_isfinite(vnl_finite_int<N> const& x){return true;}
00406 
00407 //: finite modulo-N arithmetic with polynomials of degree < M
00408 //
00409 // This class implements arithmetic with polynomials of degree < M over
00410 // vnl_finite_int<N>. Multiplication is defined modulo a polynomial of degree M.
00411 //
00412 // For N prime, and when the "modulo" polynomial is irreducible,
00413 // vnl_finite_int_poly<N,M> implements the finite field GF(N^M).
00414 //
00415 // Addition, subtraction and scalar multiplication are already defined without
00416 // the presence of a "modulo" polynomial.  Restricted to these operations,
00417 // vnl_finite_int_poly<N,M> forms an M-dimensional vector space over
00418 // vnl_finite_int<N>.  The current implementation does not yet implement
00419 // anything more than that.
00420 //
00421 template <int N, int M>
00422 class vnl_finite_int_poly
00423 {
00424   typedef vnl_finite_int_poly<N,M> Base;
00425   typedef vnl_finite_int<N> Scalar;
00426 
00427   vcl_vector<Scalar> val_; //!< M-tuple (or degree M-1 polynomial) representing this
00428 
00429   // This essentially implements std::pow(N,int) which is not always available
00430   static unsigned int Ntothe(unsigned int m) { return m==0?1:m==1?N:Ntothe(m/2)*Ntothe((m+1)/2); }
00431  public:
00432   //: The number of different finite_int polynomials of this type
00433   static unsigned int cardinality() { return Ntothe(M); }
00434 
00435   //: Creates a general finite_int_poly.
00436   inline vnl_finite_int_poly(vcl_vector<Scalar> const& p) : val_(p) { assert(N>1); assert(M>0); assert(p.size()<=M); }
00437   //: Creates a degree 0 finite_int_poly.
00438   inline vnl_finite_int_poly(Scalar const& n) : val_(vcl_vector<Scalar>(1)) { assert(N>1); assert(M>0); val_[0]=n; }
00439   //: Default constructor. Creates a degree 0 finite_int_poly equal to 0.
00440   inline vnl_finite_int_poly() : val_(vcl_vector<Scalar>(1)) { assert(N>1); assert(M>0); val_[0]=0; }
00441   //  Copy constructor
00442   inline vnl_finite_int_poly(Base const& x) : val_(x.val_) {}
00443   //  Destructor
00444   inline ~vnl_finite_int_poly() {}
00445 
00446   //: Formal degree of this polynomial
00447   inline vcl_size_t deg() const { return val_.size() - 1; }
00448 
00449   //: Effective degree of this polynomial; equals -1 when this polynomial is 0.
00450   int degree() const { for (int i=deg(); i>=0; --i) if (val_[i]!=0) return i; return -1; }
00451 
00452   //: Access to individual coefficients
00453   inline Scalar operator[](unsigned int i) const { assert(i<M); return i<=deg() ? val_[i] : Scalar(0); }
00454 
00455   //: Assignment
00456   inline Base& operator=(Base const& x) { val_ = x.val_; return *this; }
00457   inline Base& operator=(Scalar const& n) { val_ = vcl_vector<Scalar>(1); val_[0] = n; return *this; }
00458 
00459   //: Comparison of finite int polys.
00460   bool operator==(Base const& x) const {
00461     for (unsigned int i=0; i<=deg(); ++i)
00462       if (val_[i] != x[i]) return false;
00463     for (unsigned int i=deg()+1; i<=x.deg(); ++i)
00464       if (x[i] != 0) return false;
00465     return true;
00466   }
00467   inline bool operator!=(Base const& x) const { return !operator==(x); }
00468   bool operator==(Scalar const& x) const {
00469     if (x!=val_[0]) return false;
00470     for (unsigned int i=1; i<=deg(); ++i) if (val_[i] != 0) return false;
00471     return true;
00472   }
00473   inline bool operator!=(Scalar const& x) const { return !operator==(x); }
00474 
00475   //: Unary minus - returns the additive inverse
00476   Base operator-() const { vcl_vector<Scalar> p = val_; for (unsigned int i=0; i<p.size(); ++i) p[i]=-p[i]; return p; }
00477   //: Unary plus - returns the current polynomial
00478   inline Base operator+() const { return *this; }
00479   //: Unary not - returns true if finite int poly is equal to zero.
00480   bool operator!() const { for (unsigned int i=0; i<=deg(); ++i) if (val_[i] != 0) return false; return true; }
00481 
00482   //: Plus&assign: replace lhs by lhs + rhs
00483   Base& operator+=(Base const& r) {
00484     for (unsigned int i=0; i<=r.deg(); ++i)
00485       if (i<=deg()) val_[i] += r[i];
00486       else          val_.push_back(r[i]);
00487     return *this;
00488   }
00489   //: Minus&assign: replace lhs by lhs - rhs
00490   Base& operator-=(Base const& r) {
00491     for (unsigned int i=0; i<=r.deg(); ++i)
00492       if (i<=deg()) val_[i] -= r[i];
00493       else          val_.push_back(-r[i]);
00494     return *this;
00495   }
00496 
00497   //: Scalar multiple of this
00498   Base& operator*=(Scalar const& n) { for (unsigned int i=0; i<=deg(); ++i) val_[i] *= n; return *this; }
00499 
00500   //: The additive order of x is the smallest positive r such that r*x == 0.
00501   unsigned int additive_order() const {
00502     unsigned int r = N;
00503     for (unsigned int i=0; i<=deg(); ++i)
00504       if (val_[i] != 0) r=Scalar::gcd(val_[i],r);
00505     return N/r;
00506   }
00507 
00508   //: get/set the (irreducible) modulo polynomial of degree M
00509   //  Note that this polynomial has to be set only once, i.e., once set,
00510   //  it applies to all vnl_finite_int_polys with the same N and M.
00511   static vcl_vector<Scalar>& modulo_polynomial(vcl_vector<Scalar> p = vcl_vector<Scalar>())
00512   {
00513     static vcl_vector<Scalar> poly_(M+1, Scalar(0));
00514     if (p.size() == 0) { // retrieval
00515       assert(poly_[M] != 0); // cannot retrieve before having set
00516     }
00517     else
00518     {
00519       assert(p.size() == M+1 && p[M].is_unit());// must be of effective degree M
00520       // Now set poly_, thereby making the coefficient poly_[M] equal to -1.
00521       Scalar f = -1/p[M];
00522       for (int m=0; m<=M; ++m) poly_[m] = f*p[m];
00523     }
00524     return poly_;
00525   }
00526 
00527   //: Multiply&assign: replace lhs by lhs * rhs, modulo the "modulo" polynomial
00528   Base& operator*=(Base const& r) {
00529     Base x = *this; *this = r; *this *= x[0];
00530     while (val_.size() < M) val_.push_back(0);
00531     for (int i=1; i<=x.degree(); ++i)
00532       for (unsigned int j=0; j<=r.deg(); ++j)
00533         add_modulo_poly(i+j, x[i]*r[j]);
00534     return *this;
00535   }
00536 
00537   //: Return the multiplicative order of this polynomial.
00538   inline unsigned int multiplicative_order() const {
00539     Base t = Scalar(1);
00540     unsigned int order = 0;
00541     do t *= *this, ++order; while (t != Scalar(1) && t != Scalar(0));
00542     return t==Scalar(1) ? order : -1;
00543   }
00544 
00545   //: Return true when this ring is a field.
00546   //  Note that this requires that N is a prime, but that condition is not
00547   //  sufficient: also the "modulo" polynomial must be irreducible.
00548   static inline bool is_field() {
00549     if (!Scalar::is_field()) return false;
00550 
00551     vcl_vector<Scalar> mod_p = Base::modulo_polynomial();
00552     mod_p.pop_back(); // remove the "-1" coefficient of X^M
00553     return Base(mod_p).multiplicative_order()+1 == Base::cardinality();
00554   }
00555 
00556   //: Return the smallest multiplicative generator of the units in this ring.
00557   //  This is only possible if the units form a cyclic group for multiplication.
00558   //  If not, smallest_generator() returns 1 to indicate this fact.
00559   //  Note that the multiplication group of a finite \e field is always cyclic.
00560   static Base smallest_generator() {
00561     if (!Base::is_field()) return Scalar(1);
00562     vcl_vector<Scalar> mod_p = Base::modulo_polynomial();
00563     mod_p.pop_back(); // remove the "-1" coefficient of X^M
00564     return Base(mod_p);
00565   }
00566 
00567  private:
00568   //: Add x to the i-th degree coefficient of val_, possibly reducing modulo the "modulo" poly.
00569   void add_modulo_poly(unsigned int m, Scalar const& x)
00570   {
00571     if (m < M) val_[m] += x;
00572     else {
00573       vcl_vector<Scalar> p = modulo_polynomial(); // where p[M] == -1
00574       for (int k=0; k<M; ++k) add_modulo_poly(m-M+k, x*p[k]); // recursive call
00575     }
00576   }
00577 };
00578 
00579 //: Returns the sum of two finite int polynomials.
00580 // \relatesalso vnl_finite_int_poly
00581 template <int N, int M>
00582 inline vnl_finite_int_poly<N,M> operator+ (vnl_finite_int_poly<N,M> const& r1, vnl_finite_int_poly<N,M> const& r2)
00583 {
00584   vnl_finite_int_poly<N,M> result=r1; return result += r2;
00585 }
00586 
00587 //: Returns the difference of two finite int polynomials.
00588 // \relatesalso vnl_finite_int_poly
00589 template <int N, int M>
00590 inline vnl_finite_int_poly<N,M> operator- (vnl_finite_int_poly<N,M> const& r1, vnl_finite_int_poly<N,M> const& r2)
00591 {
00592   vnl_finite_int_poly<N,M> result=r1; return result -= r2;
00593 }
00594 
00595 //: Returns a scalar multiple of a finite int polynomial.
00596 // \relatesalso vnl_finite_int
00597 // \relatesalso vnl_finite_int_poly
00598 template <int N, int M>
00599 inline vnl_finite_int_poly<N,M> operator* (vnl_finite_int_poly<N,M> const& r1, vnl_finite_int<N> const& r2)
00600 {
00601   vnl_finite_int_poly<N,M> result(r1); return result *= r2;
00602 }
00603 
00604 //: Returns a scalar multiple of a finite int polynomial.
00605 // \relatesalso vnl_finite_int
00606 // \relatesalso vnl_finite_int_poly
00607 template <int N, int M>
00608 inline vnl_finite_int_poly<N,M> operator* (vnl_finite_int<N> const& r2, vnl_finite_int_poly<N,M> const& r1)
00609 {
00610   vnl_finite_int_poly<N,M> result(r1); return result *= r2;
00611 }
00612 
00613 //: Multiplies two finite int polynomials.
00614 //  NOTE: this requires the "modulo" polynomial to be set.
00615 //  Do this by calling modulo_polynomial(p), where p is a vector of length M+1.
00616 // \relatesalso vnl_finite_int_poly
00617 template <int N, int M>
00618 inline vnl_finite_int_poly<N,M> operator* (vnl_finite_int_poly<N,M> const& r1, vnl_finite_int_poly<N,M> const& r2)
00619 {
00620   vnl_finite_int_poly<N,M> result(r1); return result *= r2;
00621 }
00622 
00623 //: formatted output
00624 // \relatesalso vnl_finite_int_poly
00625 template <int N, int M>
00626 inline vcl_ostream& operator<< (vcl_ostream& s, vnl_finite_int_poly<N,M> const& r)
00627 {
00628   bool out = false;
00629   for (unsigned int i=0; i<=r.deg(); ++i) {
00630     if (r[i] == 0) continue;
00631     if (out) s << '+';
00632     out = true;
00633     if (r[i] != 1 || i==0) s << r[i];
00634     if (i>0) s << 'X';
00635     if (i>1) s << '^' << i;
00636   }
00637   if (!out) s << '0';
00638   return s;
00639 }
00640 
00641 #endif // vnl_finite_h_