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