Real algebraic numbers

00001 /*********************************************************************
00002 *      This file is part of the source code of SYNAPS kernel.        *
00003 *   Author(s): B. Mourrain, GALAAD, INRIA                            *
00004 **********************************************************************/
00005 #ifndef SYNAPS_ARITHM_ZP_H
00006 #define SYNAPS_ARITHM_ZP_H
00007 //--------------------------------------------------------------------
00008 #include <iostream>
00009 #include <cstdlib>
00010 #include <synaps/init.h>
00011 #include <synaps/arithm/let.h>
00012 //--------------------------------------------------------------------
00013 __BEGIN_NAMESPACE_SYNAPS
00014 //--------------------------------------------------------------------
00022 template <unsigned int p, class T=int>
00023 class Z
00024 {
00025 private:
00026   T rep_;
00027 
00028 public:
00029   Z(): rep_(0) {}
00030   Z(int n): rep_(n%p){ if(rep_<0) rep_+=p; }
00031   Z(const char* s);
00032   Z(const Z<p,T> & z): rep_(z.rep_) {}
00033 
00034   T rep() const {return rep_;}
00035   Z<p,T>&  operator= (const Z<p,T>  z)
00036     {rep_=z.rep_; return *this;}
00037   Z<p,T>&  operator+=(const Z<p,T>  z)
00038     {rep_ += z.rep_; if(rep_>=(int)p) rep_-=p; return *this;}
00039   Z<p,T>&  operator-=(const Z<p,T>  z)
00040     {rep_ -= z.rep_; if(rep_<0) rep_+=p; return *this;}
00041   Z<p,T>&  operator*=(const Z<p,T>  z)
00042     {(rep_ *= z.rep_)%= p; if(rep_<0) rep_+=p; return *this;}
00043   Z<p,T>&  operator/=(const Z<p,T> &  z);
00044   bool   operator==(const Z<p,T> & z) const {return (rep_==z.rep_);}
00045   bool   operator!=(const Z<p,T> & z) const {return !(*this==z);}
00046 
00047   friend inline
00048   Z<p,T> operator -(const Z<p,T> & a) {return Z<p,T>(-a.rep_);}
00049   friend inline
00050   Z<p,T> operator+(const Z<p,T>  & a, const  Z<p,T> & b){return (Z<p,T>(a.rep_+b.rep_));}
00051   friend inline
00052   Z<p,T> operator-(const Z<p,T> & a, const  Z<p,T> & b){return (Z<p,T>(a.rep_-b.rep_));}
00053   friend inline
00054   Z<p,T> operator*(const Z<p,T> & a, const  Z<p,T> & b){return (Z<p,T>(a.rep_*b.rep_));}
00055   friend inline
00056   Z<p,T> operator/(const Z<p,T> & a, const  Z<p,T> & b){return (Z<p,T>(a)/=b);}
00057 
00058   friend
00059   std::ostream &  operator<<(std::ostream&  os, const Z<p,T> & z) {return (os << z.rep_);}
00060   friend
00061   std::istream & operator>>(std::istream&  is, Z<p,T> & z) {return (is >> z.rep_);}
00062 
00063 };
00064 
00065 
00066 //----------------------------------------------------------------------
00068 //----------------------------------------------------------------------
00069 template <unsigned int p, class T>
00070 Z<p,T>::Z(const char* s) :rep_(0)
00071 {
00072   int i=0; while(s[i]!='\0') i++;
00073   unsigned b=1;
00074   while(i>0){
00075     rep_+=(s[i-1]-'0')*b; rep_%=p;
00076     b*=10; b%=p;
00077     i--;
00078   }
00079 }
00080 //----------------------------------------------------------------------
00082 //----------------------------------------------------------------------
00083 template <unsigned int p, class T>
00084 Z<p,T> & Z<p,T>::operator/=(const Z<p,T> & z)
00085 {
00086   T r=z.rep_, b=z.rep_, a=p,u1=rep_,u2=rep_, q, s=1;
00087   if(z==0){ std::cerr<< "xxxx division by zero"<<std::endl;}
00088   rep_=0;
00089   if(a<0) {s=-1; a*=s;}
00090   while(r!=0)
00091     {
00092       q = a/b;
00093       r = a-q*b; a=b; b=r;
00094       u2= rep_-q*u1; rep_=u1; u1=u2;
00095     }
00096   rep_ =(rep_ %p + p)%p;  rep_ *=s;
00097   return *this;
00098 }
00099 //----------------------------------------------------------------------
00101 //----------------------------------------------------------------------
00102 template <unsigned int p, class T>
00103 Z<p,T> operator^(const Z<p,T> a, const  int n)
00104 {
00105     int m;
00106     Z<p,T> r=1, e;
00107     if(n <0){e=1/a;m=-n;}else {e=a;m=n;}
00108     while(m !=0) {
00109       if(m%2){
00110         m -=1; m/=2; r*=e; e*=e;
00111       }else{
00112         m/=2; e*=e;
00113       }
00114     }
00115     return r;
00116 }
00117 //--------------------------------------------------------------------
00118 template <unsigned int p, class T>
00119 bool with_plus_sign(const Z<p,T>& z)
00120 { 
00121   return (z.rep()>0);
00122 }
00123 //--------------------------------------------------------------------
00124 __END_NAMESPACE_SYNAPS
00125 //--------------------------------------------------------------------
00126 #endif // SYNAPS_ARITHM_ZP_H
00127
synaps/arithm/Zp.h
