1 // ========================================================================== 2 // Copyright(c)'1994-2015 by The Givaro group 3 // This file is part of Givaro. 4 // Givaro is governed by the CeCILL-B license under French law 5 // and abiding by the rules of distribution of free software. 6 // see the COPYRIGHT file for more details. 7 // Authors: W. J. Turner <wjturner@acm.org> 8 // Bradford Hovinen <hovinen@cis.udel.edu> 9 // Clement Pernet <clement.pernet@gmail.com> (inserted into FFLAS-FFPACK) 10 // A. Breust (taken from FFLAS-FFPACK) 11 // ========================================================================== 12 13 14 /*! @file field/zring.h 15 * @ingroup field 16 * @brief representation of a field of characteristic 0. 17 */ 18 19 #ifndef __GIVARO_ring_zring_H 20 #define __GIVARO_ring_zring_H 21 22 #include <algorithm> 23 #include <typeinfo> 24 #include <math.h> 25 26 #include "givaro/unparametric-operations.h" 27 #include "givaro/givranditer.h" 28 #include "givaro/givinteger.h" 29 30 namespace Givaro 31 { 32 template<class _Element> class ZRing; 33 34 template<typename Domain> struct DomainRandIter { 35 typedef GeneralRingRandIter<Domain> RandIter; 36 }; 37 38 template<> struct DomainRandIter<ZRing<Integer>> { 39 typedef IntegerDom::RandIter RandIter; 40 }; 41 42 43 /** Class ZRing. 44 * Ring of integers, using the _Element base type. 45 */ 46 template<class _Element> 47 class ZRing : public UnparametricOperations<_Element> 48 { 49 public: 50 51 // ----- Exported Types and constantes 52 using Element = _Element; 53 using Rep = _Element; 54 using Self_t = ZRing<Element>; 55 using Residu_t = Element; 56 using Element_ptr = Element*; 57 using ConstElement_ptr = const Element*; 58 enum { size_rep = sizeof(Residu_t) }; 59 60 const Element one = 1; 61 const Element zero = 0; 62 const Element mOne = -1; 63 64 //----- Constructors 65 ZRing() {} 66 ZRing(const ZRing& F) {} 67 // Needed in FFLAS, when ZRing is used as delayed field. 68 template<class T> ZRing(const T&) {} 69 70 //----- Access 71 Residu_t residu() const { return static_cast<Residu_t>(0); } 72 Residu_t size() const { return static_cast<Residu_t>(0); } 73 Residu_t cardinality() const { return static_cast<Residu_t>(0); } 74 Residu_t characteristic() const { return static_cast<Residu_t>(0); } 75 template<typename T> T& cardinality(T& c) const { return c = static_cast<T>(0); } 76 template<typename T> T& characteristic(T& c) const { return c = static_cast<T>(0); } 77 78 static inline Residu_t maxCardinality() { return -1; } 79 static inline Residu_t minCardinality() { return 2; } 80 81 //----- Ring-wise operations 82 inline bool operator==(const Self_t& F) const { return true; } 83 inline bool operator!=(const Self_t& F) const { return false; } 84 inline ZRing<Element>& operator=(const ZRing<Element>&) { return *this; } 85 // Ring tests 86 bool isZero(const Element& a) const { return a == zero; } 87 bool isOne (const Element& a) const { return a == one; } 88 bool isMOne(const Element& a) const { return a == mOne; } 89 bool isUnit(const Element& a) const { return isOne(a) || isMOne(a); } 90 91 Element& abs(Element& x, const Element& a) const {return x= (a>0)? a: -a;} 92 Element abs(const Element& a) const {return (a>0)? a: -a;} 93 long compare (const Element& a, const Element& b) const {return (a>b)? 1: ((a<b)? -1 : 0);} 94 Element& gcd (Element& g, const Element& a, const Element& b) const {return Givaro::gcd(g,a,b);} 95 Element& gcdin (Element& g, const Element& b) const{return gcd(g, g, b);} 96 Element& gcd (Element& g, Element& s, Element& t, const Element& a, const Element& b) const{return Givaro::gcd(g,s,t,a,b);} 97 Element &dxgcd(Element &g, Element &s, Element &t, Element &u, Element &v, const Element &a, const Element &b) const 98 { 99 gcd(g,s,t,a,b); 100 div(u,a,g); 101 div(v,b,g); 102 return g; 103 } 104 Element& lcm (Element& c, const Element& a, const Element& b) const 105 { 106 if ((a==Element(0)) || (b==Element(0))) return c = Element(0); 107 else { 108 Element g; 109 gcd (g, a, b); 110 c= a*b; 111 c /= g; 112 c=abs (c); 113 return c; 114 } 115 } 116 117 Element& lcmin (Element& l, const Element& b) const 118 { 119 if ((l==Element(0)) || (b==Element(0))) return l = Element(0); 120 else { 121 Element g; 122 gcd (g, l, b); 123 l*= b; 124 l/= g; 125 l=abs (l); 126 return l; 127 } 128 } 129 130 void reconstructRational (Element& a, Element& b, const Element& x, const Element& m) const 131 {this->RationalReconstruction(a,b, x, m, Givaro::sqrt(m), true, true);} 132 void reconstructRational (Element& a, Element& b, const Element& x, const Element& m, const Element& bound) const 133 {this->RationalReconstruction(a,b, x, m, bound, true, true);} 134 bool reconstructRational (Element& a, Element& b, const Element& x, const Element& m, const Element& a_bound, const Element& b_bound) const 135 { 136 Element bound = x/b_bound; 137 this->RationalReconstruction(a,b,x,m, (bound>a_bound?bound:a_bound), true, false); 138 return b <= b_bound; 139 } 140 Element& quo (Element& q, const Element& a, const Element& b) const {return Integer::floor(q, a, b);} 141 Element& rem (Element& r, const Element& a, const Element& b) const {return Integer::mod(r,a,b);} 142 Element& quoin (Element& a, const Element& b) const{return quo(a,a,b);} 143 Element& remin (Element& a, const Element& b) const {return rem(a,a,b);} 144 void quoRem (Element& q, Element& r, const Element& a, const Element& b) const{Integer::divmod(q,r,a,b);} 145 bool isDivisor (const Element& a, const Element& b) const 146 { 147 Element r; 148 return rem(r,a,b)==Element(0); 149 } 150 151 inline Element& sqrt(Element& x, const Element& y) const{return Givaro::sqrt(x,y);} 152 inline Element powtwo(Element& z, const Element& x) const 153 { 154 z = 1; 155 Element max; init(max, (int64_t)(1<<30)); 156 if (x < 0) return z; 157 //if (x < (Element)max-1) { 158 if (x < max) { 159 z<<=(int32_t)x; 160 return z; 161 } 162 else { 163 Element n,m; 164 quoRem(n,m,x,max); 165 //quoRem(n,m,x,(Element)(LONG_MAX-1)); 166 for (int i=0; i < n; ++i) { 167 z <<=(int32_t)max; 168 //z <<=(long int)(LONG_MAX-1); 169 } 170 z <<= (int32_t)m; 171 return z; 172 } 173 //for (Element i=0; i < x; ++i) { 174 // z <<= 1; 175 //} 176 //return z; // BB peut pas ! 177 } 178 179 inline Element logtwo(Element& z, const Element& x) const {return z = x.bitsize() - 1;} 180 181 182 //----- Initialisation 183 Element& init(Element& x) const { return x; } 184 template <typename T> Element& init(Element& x, const T& s) const 185 { return Caster(x,s); } 186 187 Element& assign(Element& x, const Element& y) const { return x = y; } 188 189 //----- Convert 190 template <typename T> T& convert(T& x, const Element& y) const 191 { return Caster(x,y); } 192 193 Element& reduce (Element& x, const Element& y) const { return x = y; } 194 Element& reduce (Element& x) const { return x; } 195 196 // To ensure interface consistency 197 Element minElement() const { return 0; } 198 Element maxElement() const { return 0; } 199 200 201 // ----- Random generators 202 typedef typename DomainRandIter<Self_t>::RandIter RandIter; 203 typedef GeneralRingNonZeroRandIter<Self_t> NonZeroRandIter; 204 template< class Random > Element& random(const Random& g, Element& r) const 205 { return init(r, g()); } 206 template< class Random > Element& nonzerorandom(const Random& g, Element& a) const 207 { while (isZero(init(a, g()))) 208 ; 209 return a; } 210 211 //----- IO 212 std::ostream& write(std::ostream &os) const 213 { 214 return os << "ZRing<" << typeid(Element).name() << ">"; 215 } 216 std::ostream& write(std::ostream &os, const Element& a) const 217 { 218 return os << a; 219 } 220 std::istream& read(std::istream &is, Element& a) const 221 { 222 return is >> a; 223 } 224 protected: 225 /*! Rational number reconstruction. 226 * \f$\frac{n}{d} \equiv f \mod m\f$, with \f$\vert n 227 \vert <k\f$ and \f$0 < \vert d \vert \leq \frac{f}{k}\f$. 228 * @bib 229 * - von zur Gathen & Gerhard, <i>Modern Computer Algebra</i>, 230 * 5.10, Cambridge Univ. Press 1999 231 */ 232 inline void RationalReconstruction( Element& a, Element& b, 233 const Element& f, const Element& m, 234 const Element& k, 235 bool forcereduce, bool recursive ) const 236 { 237 Element x(f); 238 if (x<0) { 239 if ((-x)>m) 240 x %= m; 241 if (x<0) 242 x += m; 243 } 244 else { 245 if (x>m) 246 x %= m; 247 } 248 249 if (x == 0) { 250 a = 0; 251 b = 1; 252 } 253 else { 254 bool res = ratrecon(a,b,x,m,k, forcereduce, recursive); 255 if (recursive) 256 for( Element newk = k + 1; (!res) && (newk<f) ; ++newk) 257 res = ratrecon(a,b,x,m,newk,forcereduce, true); 258 } 259 } 260 261 // Precondition f is suppposed strictly positive and strictly less than m 262 inline bool ratrecon( Element& num, Element& den, 263 const Element& f, const Element& m, 264 const Element& k, 265 bool forcereduce, bool recursive ) const 266 { 267 268 //std::cerr << "RatRecon : " << f << " " << m << " " << k << std::endl; 269 Element r0, t0, q, u; 270 r0=m; 271 t0=0; 272 num=f; 273 den=1; 274 while(num>=k) 275 { 276 q = r0; 277 q /= num; // r0/num 278 u = num; 279 num = r0; // num <-- r0 280 r0 = u; // r0 <-- num 281 //maxpyin(num,u,q); 282 Integer::maxpyin(num,u,q); 283 if (num == 0) return false; 284 285 u = den; 286 den = t0; // num <-- r0 287 t0 = u; // r0 <-- num 288 //maxpyin(den,u,q); 289 Integer::maxpyin(den,u,q); 290 } 291 292 if (forcereduce) { 293 // [GG, MCA, 1999] Theorem 5.26 294 295 // (ii) 296 Element gg; 297 if (gcd(gg,num,den) != 1) { 298 299 Element ganum, gar2; 300 for( q = 1, ganum = r0-num, gar2 = r0 ; (ganum < k) && (gar2>=k); ++q ) { 301 ganum -= num; 302 gar2 -= num; 303 } 304 305 //maxpyin(r0,q,num); 306 Integer::maxpyin(r0,q,num); 307 //maxpyin(t0,q,den); 308 Integer::maxpyin(t0,q,den); 309 310 if (t0 < 0) { 311 num = -r0; 312 den = -t0; 313 } 314 else { 315 num = r0; 316 den = t0; 317 } 318 319 // if (t0 > m/k) 320 if (den > m/k) { 321 if (!recursive) 322 std::cerr 323 << "*** Error *** No rational reconstruction of " 324 << f 325 << " modulo " 326 << m 327 << " with denominator <= " 328 << (m/k) 329 << std::endl; 330 } 331 if (gcd(gg,num,den) != 1) { 332 if (!recursive) 333 std::cerr 334 << "*** Error *** There exists no rational reconstruction of " 335 << f 336 << " modulo " 337 << m 338 << " with |numerator| < " 339 << k 340 << std::endl 341 << "*** Error *** But " 342 << num 343 << " = " 344 << den 345 << " * " 346 << f 347 << " modulo " 348 << m 349 << std::endl; 350 return false; 351 } 352 } 353 } 354 // (i) 355 if (den < 0) { 356 Integer::negin(num); 357 Integer::negin(den); 358 } 359 360 // std::cerr << "RatRecon End " << num << "/" << den << std::endl; 361 return true; 362 } 363 }; 364 365 typedef ZRing<float> FloatDomain; 366 typedef ZRing<double> DoubleDomain; 367 typedef ZRing<Integer> IntegerDomain; 368 } 369 370 #endif 371 /* -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ 372 // vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s 373