vnl_finite.h

InsightToolkit-1.4.0(有大量的优化算法程序)

}

//: Returns the difference of two finite int numbers.
// \relates vnl_finite_int
template <int N>
inline vnl_finite_int<N> operator- (vnl_finite_int<N> const& r1, vnl_finite_int<N> const& r2) {
  vnl_finite_int<N> result(r1); return result -= r2;
}
template <int N>
inline vnl_finite_int<N> operator- (vnl_finite_int<N> const& r1, int r2) {
  vnl_finite_int<N> result(r1); return result -= r2;
}
template <int N>
inline vnl_finite_int<N> operator- (int r2, vnl_finite_int<N> const& r1) {
  vnl_finite_int<N> result(-r1); return result += r2;
}

//: Returns the product of two finite int numbers.
// \relates vnl_finite_int
template <int N>
inline vnl_finite_int<N> operator* (vnl_finite_int<N> const& r1, vnl_finite_int<N> const& r2) {
  vnl_finite_int<N> result(r1); return result *= r2;
}
template <int N>
inline vnl_finite_int<N> operator* (vnl_finite_int<N> const& r1, int r2) {
  vnl_finite_int<N> result(r1); return result *= r2;
}
template <int N>
inline vnl_finite_int<N> operator* (int r2, vnl_finite_int<N> const& r1) {
  vnl_finite_int<N> result(r1); return result *= r2;
}
//: Returns the quotient of two finite int numbers.
//  Uses r2.reciproc() for efficient computation.
// \relates vnl_finite_int
template <int N>
inline vnl_finite_int<N> operator/(vnl_finite_int<N> const& r1, vnl_finite_int<N> const& r2) {
  assert(r2.is_unit()); return r1 == 0 ? vnl_finite_int<N>(0) : r1*r2.reciproc();
}
template <int N>
inline vnl_finite_int<N> operator/ (vnl_finite_int<N> const& r1, int r2) {
  vnl_finite_int<N> result(r1); return result /= r2;
}
template <int N>
inline vnl_finite_int<N> operator/ (int r1, vnl_finite_int<N> const& r2) {
  vnl_finite_int<N> result(r1); return result /= r2;
}

template <int N>
inline bool operator== (int  r1, vnl_finite_int<N> const& r2) { return r2==r1; }
template <int N>
inline bool operator!= (int  r1, vnl_finite_int<N> const& r2) { return r2!=r1; }

//:
// \relates vnl_finite_int
template <int N>
inline vnl_finite_int<N> vnl_math_squared_magnitude(vnl_finite_int<N> const& x) { return x*x; }
template <int N>
inline vnl_finite_int<N> vnl_math_sqr(vnl_finite_int<N> const& x) { return x*x; }
template <int N>
inline bool vnl_math_isnan(vnl_finite_int<N> const& ){return false;}
template <int N>
inline bool vnl_math_isfinite(vnl_finite_int<N> const& x){return true;} 

//: finite modulo-N arithmetic with polynomials of degree < M
//
// This class implements arithmetic with polynomials of degree < M over
// vnl_finite_int<N>. Multiplication is defined modulo a polynomial of degree M.
//
// For N prime, and when the "modulo" polynomial is irreducible,
// vnl_finite_int_poly<N,M> implements the finite field GF(N^M).
//
// Addition, subtraction and scalar multiplication are already defined without
// the presence of a "modulo" polynomial.  Restricted to these operations,
// vnl_finite_int_poly<N,M> forms an M-dimensional vector space over
// vnl_finite_int<N>.  The current implementation does not yet implement
// anything more than that.
//
template <int N, int M>
class vnl_finite_int_poly {
  typedef vnl_finite_int_poly<N,M> Base;
  typedef vnl_finite_int<N> Scalar;

  vcl_vector<Scalar> val_; //!< M-tuple (or degree M-1 polynomial) representing this

  // This essentially implements std::pow(int,int) which is not always available
  static unsigned int Ntothe(unsigned int m) { return m==0?1:N*Ntothe(m-1); }
 public:
  //: The number of different finite_int polynomials of this type
  static unsigned int cardinality() { return Ntothe(M); }

  //: Creates a general finite_int_poly.
  inline vnl_finite_int_poly(vcl_vector<Scalar> const& p) : val_(p) { assert(N>1); assert(M>0); assert(p.size()<=M); }
  //: Creates a degree 0 finite_int_poly.
  inline vnl_finite_int_poly(Scalar const& n) : val_(vcl_vector<Scalar>(1)) { assert(N>1); assert(M>0); val_[0]=n; }
  //: Default constructor. Creates a degree 0 finite_int_poly equal to 0.
  inline vnl_finite_int_poly() : val_(vcl_vector<Scalar>(1)) { assert(N>1); assert(M>0); val_[0]=0; }
  //  Copy constructor
  inline vnl_finite_int_poly(Base const& x) : val_(x.val_) {}
  //  Destructor
  inline ~vnl_finite_int_poly() {}

  //: Formal degree of this polynomial
  inline unsigned int deg() const { return val_.size() - 1; }

  //: Effective degree of this polynomial; equals -1 when this polynomial is 0.
  inline int degree() const { for (int i=deg(); i>=0; --i) if (val_[i]!=0) return i; return -1; }

  //: Access to individual coefficients
  inline Scalar operator[](unsigned int i) const { assert(i<M); return i<=deg() ? val_[i] : Scalar(0); }

  //: Assignment
  inline Base& operator=(Base const& x) { val_ = x.val_; return *this; }
  inline Base& operator=(Scalar const& n) { val_ = vcl_vector<Scalar>(1); val_[0] = n; return *this; }

  //: Comparison of finite int polys.
  inline bool operator==(Base const& x) const {
    for (unsigned int i=0; i<=deg(); ++i)
      if (val_[i] != x[i]) return false;
    for (unsigned int i=deg()+1; i<=x.deg(); ++i)
      if (x[i] != 0) return false;
    return true;
  }
  inline bool operator!=(Base const& x) const { return !operator==(x); }
  inline bool operator==(Scalar const& x) const {
    if (x!=val_[0]) return false;
    for (unsigned int i=1; i<=deg(); ++i) if (val_[i] != 0) return false;
    return true;
  }
  inline bool operator!=(Scalar const& x) const { return !operator==(x); }

  //: Unary minus - returns the additive inverse
  inline Base operator-() const { vcl_vector<Scalar> p = val_; for (unsigned int i=0; i<p.size(); ++i) p[i]=-p[i]; return p; }
  //: Unary plus - returns the current polynomial
  inline Base operator+() const { return *this; }
  //: Unary not - returns true if finite int poly is equal to zero.
  inline bool operator!() const { for (unsigned int i=0; i<=deg(); ++i) if (val_[i] != 0) return false; return true; }

  //: Plus&assign: replace lhs by lhs + rhs
  inline Base& operator+=(Base const& r) {
    for (unsigned int i=0; i<=r.deg(); ++i)
      if (i<=deg()) val_[i] += r[i];
      else          val_.push_back(r[i]);
    return *this;
  }
  //: Minus&assign: replace lhs by lhs - rhs
  inline Base& operator-=(Base const& r) {
    for (unsigned int i=0; i<=r.deg(); ++i)
      if (i<=deg()) val_[i] -= r[i];
      else          val_.push_back(-r[i]);
    return *this;
  }

  //: Scalar multiple of this
  inline Base& operator*=(Scalar const& n) { for (unsigned int i=0; i<=deg(); ++i) val_[i] *= n; return *this; }

  //: The additive order of x is the smallest positive r such that r*x == 0.
  inline unsigned int additive_order() const {
    unsigned int r = N;
    for (unsigned int i=0; i<=deg(); ++i)
      if (val_[i] != 0) r=Scalar::gcd(val_[i],r);
    return N/r;
  }

  //: get/set the (irreducible) modulo polynomial of degree M
  //  Note that this polynomial has to be set only once, i.e., once set,
  //  it applies to all vnl_finite_int_polys with the same N and M.
  static vcl_vector<Scalar>& modulo_polynomial(vcl_vector<Scalar> p = vcl_vector<Scalar>())
  {
    static vcl_vector<Scalar> poly_(M+1, Scalar(0));
    if (p.size() == 0) { // retrieval
      assert(poly_[M] != 0); // cannot retrieve before having set
    }
    else
    {
      assert(p.size() == M+1 && p[M].is_unit());// must be of effective degree M
      // Now set poly_, thereby making the coefficient poly_[M] equal to -1.
      Scalar f = -1/p[M];
      for (int m=0; m<=M; ++m) poly_[m] = f*p[m];
    }
    return poly_;
  }

  //: Multiply&assign: replace lhs by lhs * rhs, modulo the "modulo" polynomial
  inline Base& operator*=(Base const& r) {
    Base x = *this; *this = r; *this *= x[0];
    while (val_.size() < M) val_.push_back(0);
    for (int i=1; i<=x.degree(); ++i)
      for (unsigned int j=0; j<=r.deg(); ++j)
        add_modulo_poly(i+j, x[i]*r[j]);
    return *this;
  }

  //: Return the multiplicative order of this polynomial.
  inline unsigned int multiplicative_order() const {
    Base t = Scalar(1);
    unsigned int order = 0;
    do t *= *this, ++order; while (t != Scalar(1) && t != Scalar(0));
    return t==Scalar(1) ? order : -1;
  }

  //: Return true when this ring is a field.
  //  Note that this requires that N is a prime, but that condition is not
  //  sufficient: also the "modulo" polynomial must be irreducible.
  static inline bool is_field() {
    if (!Scalar::is_field()) return false;

    vcl_vector<Scalar> mod_p = Base::modulo_polynomial();
    mod_p.pop_back(); // remove the "-1" coefficient of X^M
    return Base(mod_p).multiplicative_order()+1 == Base::cardinality();
  }

  //: Return the smallest multiplicative generator of the units in this ring.
  //  This is only possible if the units form a cyclic group for multiplication.
  //  If not, smallest_generator() returns 1 to indicate this fact.
  //  Note that the multiplication group of a finite *field* is always cyclic.
  static Base smallest_generator() {
    if (!Base::is_field()) return Scalar(1);
    vcl_vector<Scalar> mod_p = Base::modulo_polynomial();
    mod_p.pop_back(); // remove the "-1" coefficient of X^M
    return Base(mod_p);
  }

 private:
  //: Add x to the i-th degree coefficient of val_, possibly reducing modulo the "modulo" poly.
  inline void add_modulo_poly(unsigned int m, Scalar const& x)
  {
    if (m < M) val_[m] += x;
    else {
      vcl_vector<Scalar> p = modulo_polynomial(); // where p[M] == -1
      for (int k=0; k<M; ++k) add_modulo_poly(m-M+k, x*p[k]); // recursive call
    }
  }
};

//: Returns the sum of two finite int polynomials.
// \relates vnl_finite_int_poly
template <int N, int M>
inline vnl_finite_int_poly<N,M> operator+ (vnl_finite_int_poly<N,M> const& r1, vnl_finite_int_poly<N,M> const& r2) {
  vnl_finite_int_poly<N,M> result=r1; return result += r2;
}

//: Returns the difference of two finite int polynomials.
// \relates vnl_finite_int_poly
template <int N, int M>
inline vnl_finite_int_poly<N,M> operator- (vnl_finite_int_poly<N,M> const& r1, vnl_finite_int_poly<N,M> const& r2) {
  vnl_finite_int_poly<N,M> result=r1; return result -= r2;
}

//: Returns a scalar multiple of a finite int polynomial.
// \relates vnl_finite_int
// \relates vnl_finite_int_poly
template <int N, int M>
inline vnl_finite_int_poly<N,M> operator* (vnl_finite_int_poly<N,M> const& r1, vnl_finite_int<N> const& r2) {
  vnl_finite_int_poly<N,M> result(r1); return result *= r2;
}
template <int N, int M>
inline vnl_finite_int_poly<N,M> operator* (vnl_finite_int<N> const& r2, vnl_finite_int_poly<N,M> const& r1) {
  vnl_finite_int_poly<N,M> result(r1); return result *= r2;
}

//: Multiplies two finite int polynomials.
//  NOTE: this requires the "modulo" polynomial to be set.
//  Do this by calling modulo_polynomial(p), where p is a vector of length M+1.
// \relates vnl_finite_int_poly
template <int N, int M>
inline vnl_finite_int_poly<N,M> operator* (vnl_finite_int_poly<N,M> const& r1, vnl_finite_int_poly<N,M> const& r2) {
  vnl_finite_int_poly<N,M> result(r1); return result *= r2;
}

//: formatted output
// \relates vnl_finite_int_poly
template <int N, int M>
inline vcl_ostream& operator<< (vcl_ostream& s, vnl_finite_int_poly<N,M> const& r) {
  bool out = false;
  for (unsigned int i=0; i<=r.deg(); ++i) {
    if (r[i] == 0) continue;
    if (out) s << '+';
    out = true;
    if (r[i] != 1 || i==0) s << r[i];
    if (i>0) s << 'X';
    if (i>1) s << '^' << i;
  }
  if (!out) s << '0';
  return s;
}

#endif // vnl_finite_h_