00001
00002 #ifndef vnl_finite_h_
00003 #define vnl_finite_h_
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034 #include <vcl_iostream.h>
00035 #include <vcl_cassert.h>
00036 #include <vcl_vector.h>
00037 #include <vcl_cstddef.h>
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047 template <int N>
00048 class vnl_finite_int
00049 {
00050 int val_;
00051
00052 typedef vnl_finite_int<N> Base;
00053 public:
00054
00055 static unsigned int cardinality() { return N; }
00056
00057
00058
00059
00060 inline vnl_finite_int(int x = 0) : val_((x%=N)<0?N+x:x), mo_(0), lp1_(0) {assert(N>1);}
00061
00062 inline vnl_finite_int(Base const& x) : val_(int(x)), mo_(x.mo_), lp1_(x.lp1_) {}
00063
00064 inline ~vnl_finite_int() {}
00065
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
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
00078
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
00085 inline Base operator-() const { return Base(N-val_); }
00086
00087 inline Base operator+() const { return *this; }
00088
00089 inline bool operator!() const { return val_ == 0; }
00090
00091
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
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
00098 inline Base& operator*=(int r) {
00099 r %=N; if (r<0) r=N+r;
00100
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
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
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
00115
00116
00117 static unsigned int totient() {
00118 static unsigned int t_ = 0;
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
00132
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;
00143 }
00144
00145
00146 inline Base& operator/=(Base const& r) {
00147 assert(r.is_unit());
00148 return val_ == 0 ? operator=(0) : operator*=(r.reciproc());
00149 }
00150
00151
00152 inline Base& operator++() { mo_=lp1_=0; ++val_; if (val_==N) val_=0; return *this; }
00153
00154 inline Base& operator--() { mo_=lp1_=0; if (val_==0) val_=N; --val_; return *this; }
00155
00156 inline Base operator++(int){Base b=*this; mo_=lp1_=0; ++val_; if (val_==N) val_=0; return b; }
00157
00158 inline Base operator--(int){Base b=*this; mo_=lp1_=0; if (val_==0) val_=N; --val_; return b;}
00159
00160
00161 static vcl_vector<unsigned int> decompose() {
00162 static vcl_vector<unsigned int> decomposition_ = vcl_vector<unsigned int>();
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
00172 static inline bool is_field() {
00173 vcl_vector<unsigned int> d = Base::decompose();
00174 return d.size() == 1;
00175 }
00176
00177
00178
00179
00180 inline bool is_unit() const { return gcd(val_) == 1; }
00181
00182
00183 inline bool is_zero_divisor() const { return gcd(val_) != 1; }
00184
00185
00186 inline unsigned int additive_order() const { if (val_ == 0) return 1; return N/gcd(val_); }
00187
00188
00189 unsigned int multiplicative_order() const {
00190 if (mo_ != 0) return mo_;
00191 if (gcd(val_) != 1) return -1;
00192 Base y = val_;
00193 for (int r=1; r<N; ++r) { if (y==1) return mo_=r; y *= val_; }
00194 return 0;
00195 }
00196
00197
00198
00199
00200
00201 static Base smallest_generator() {
00202 static Base gen_ = 0;
00203 if (gen_ != 0) return gen_;
00204 if (N==2) return gen_=1;
00205 unsigned int h = Base::totient() / 2;
00206 for (gen_=2; gen_!=0; ++gen_) {
00207
00208 unsigned int g=h; Base y = 1, z = gen_; while (g>0) { if (g%2) y *= z; g/=2; z*=z; }
00209
00210 if (y == 1) continue;
00211
00212 if (gen_.multiplicative_order() == Base::totient()) { gen_.set_log(1); return gen_; }
00213 }
00214 assert(!Base::is_field());
00215 return gen_=1;
00216 }
00217
00218
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
00228 unsigned int log() const {
00229 if (is_zero_divisor()) return -1;
00230 if (lp1_ != 0) return lp1_-1;
00231 Base z = smallest_generator();
00232 assert(N==2||z!=1);
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;
00239 }
00240
00241
00242 static inline Base exp(int r) {
00243 Base z = smallest_generator();
00244 assert(N==2||z!=1);
00245 return z.pow(r);
00246 }
00247
00248
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
00258 void set_log(unsigned int r) const { r %= Base::totient(); lp1_ = r+1; }
00259
00260 mutable unsigned int mo_;
00261 mutable unsigned int lp1_;
00262 };
00263
00264
00265
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
00273
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
00281
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
00289
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
00297
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
00305
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
00313
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
00321
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
00329
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
00337
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
00345
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
00353
00354
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
00362
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
00370
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
00379 template <int N>
00380 inline bool operator== (int r1, vnl_finite_int<N> const& r2) { return r2==r1; }
00381
00382
00383
00384 template <int N>
00385 inline bool operator!= (int r1, vnl_finite_int<N> const& r2) { return r2!=r1; }
00386
00387
00388
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
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
00399 template <int N>
00400 inline bool vnl_math_isnan(vnl_finite_int<N> const& ){return false;}
00401
00402
00403
00404 template <int N>
00405 inline bool vnl_math_isfinite(vnl_finite_int<N> const& x){return true;}
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419
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_;
00428
00429
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
00433 static unsigned int cardinality() { return Ntothe(M); }
00434
00435
00436 inline vnl_finite_int_poly(vcl_vector<Scalar> const& p) : val_(p) { assert(N>1); assert(M>0); assert(p.size()<=M); }
00437
00438 inline vnl_finite_int_poly(Scalar const& n) : val_(vcl_vector<Scalar>(1)) { assert(N>1); assert(M>0); val_[0]=n; }
00439
00440 inline vnl_finite_int_poly() : val_(vcl_vector<Scalar>(1)) { assert(N>1); assert(M>0); val_[0]=0; }
00441
00442 inline vnl_finite_int_poly(Base const& x) : val_(x.val_) {}
00443
00444 inline ~vnl_finite_int_poly() {}
00445
00446
00447 inline vcl_size_t deg() const { return val_.size() - 1; }
00448
00449
00450 int degree() const { for (int i=deg(); i>=0; --i) if (val_[i]!=0) return i; return -1; }
00451
00452
00453 inline Scalar operator[](unsigned int i) const { assert(i<M); return i<=deg() ? val_[i] : Scalar(0); }
00454
00455
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
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
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
00478 inline Base operator+() const { return *this; }
00479
00480 bool operator!() const { for (unsigned int i=0; i<=deg(); ++i) if (val_[i] != 0) return false; return true; }
00481
00482
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
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
00498 Base& operator*=(Scalar const& n) { for (unsigned int i=0; i<=deg(); ++i) val_[i] *= n; return *this; }
00499
00500
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
00509
00510
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) {
00515 assert(poly_[M] != 0);
00516 }
00517 else
00518 {
00519 assert(p.size() == M+1 && p[M].is_unit());
00520
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
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
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
00546
00547
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();
00553 return Base(mod_p).multiplicative_order()+1 == Base::cardinality();
00554 }
00555
00556
00557
00558
00559
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();
00564 return Base(mod_p);
00565 }
00566
00567 private:
00568
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();
00574 for (int k=0; k<M; ++k) add_modulo_poly(m-M+k, x*p[k]);
00575 }
00576 }
00577 };
00578
00579
00580
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
00588
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
00596
00597
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
00605
00606
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
00614
00615
00616
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
00624
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_