1 // ==========================================================================
2 // Copyright(c)'1994-2016 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: Bastien Vialla <bastien.vialla@lirmm.fr>
8 // ==========================================================================
9 
10 #ifndef __GIVARO_MODULAR_EXTENDED_H
11 #define __GIVARO_MODULAR_EXTENDED_H
12 
13 #include "givaro/givconfig.h"
14 
15 #include "givaro/givranditer.h"
16 #include "givaro/ring-interface.h"
17 #include "givaro/modular-general.h"
18 
19 namespace Givaro{
20     /*
21      *
22      * Modular double/float allowing big moduli
23      * !!: RandIter does not works, use your own random
24      *
25      */
26     template<class _Element>
27     class ModularExtended
28     {
29     public:
30 
31         typedef _Element Element;
32         typedef Element* Element_ptr ;
33         typedef const Element ConstElement;
34         typedef const Element* ConstElement_ptr;
35         // ----- Exported Types and constantes
36         typedef ModularExtended<_Element> Self_t;
37         using Compute_t = _Element;
38         typedef uint64_t Residu_t;
39         enum { size_rep = sizeof(Residu_t) };
40 
41     private:
42         // Verkampt Split
split(const Element x,Element & x_h,Element & x_l)43         inline void split(const Element x, Element &x_h, Element &x_l) const {
44             Element c;
45             if(std::is_same<Element, double>::value){
46                 c = (Element)((1 << 27)+1);
47             }else if(std::is_same<Element, float>::value){
48                 c = (Element)((1 << 13)+1);
49             }
50             c = c*x;
51             x_h = c-(c-x);
52             x_l = x - x_h;
53         }
54 
55         // Dekker mult, a * b = s + t
mult_dekker(const Element a,const Element b,Element & s,Element & t)56         inline void mult_dekker(const Element a, const Element b, Element &s, Element &t) const{
57             s = a*b;
58             Element ah, al, bh, bl;
59             split(a, ah, al);
60             split(b, bh, bl);
61             t = (al*bl-(((s-ah*bh)-al*bh)-ah*bl));
62         }
63 
64     public:
65         // ----- Constantes
66         const Element zero = 0.0;
67         const Element one = 1.0;
68         const Element mOne = -1.0;
69 
70         // ----- Constructors
71         ModularExtended() = default;
72 
ModularExtended(const XXX & p)73         template<class XXX> ModularExtended(const XXX& p)
74         : zero(0.0), one(1.0), mOne((Element)p - 1.0), _p((Element)p), _invp((Element)1/(Element)_p), _negp(-_p), _lp((Residu_t)p)
75         {
76             assert(_p >= minCardinality());
77             assert(_p <= maxCardinality());
78         }
79 
80         // ----- Accessors
minElement()81         inline Element minElement() const { return zero; }
maxElement()82         inline Element maxElement() const { return mOne; }
83 
84         // ----- Access to the modulus
residu()85         inline Residu_t residu() const { return _lp; }
size()86         inline Residu_t size() const { return _lp; }
characteristic()87         inline Residu_t characteristic() const { return _lp; }
characteristic(T & p)88         template<class T> inline T& characteristic(T& p) const { return p = _lp; }
cardinality()89         inline Residu_t cardinality() const { return _lp; }
cardinality(T & p)90         template<class T> inline T& cardinality(T& p) const { return p = _lp; }
maxCardinality()91         static inline Residu_t maxCardinality() {
92             if(std::is_same<Element, double>::value)
93                 return 1125899906842623; // 2^(52-2) - 1
94             else if(std::is_same<Element, float>::value)
95                 return 2097151; // 2^(23-2) - 1
96         }
minCardinality()97         static inline Residu_t minCardinality() { return 2; }
98 
99         // ----- Checkers
isZero(const Element & a)100         inline bool isZero(const Element& a) const { return a == zero; }
isOne(const Element & a)101         inline bool isOne (const Element& a) const { return a == one; }
isMOne(const Element & a)102         inline bool isMOne(const Element& a) const { return a == mOne; }
103         inline bool isUnit(const Element& a) const;
areEqual(const Element & a,const Element & b)104         inline bool areEqual(const Element& a, const Element& b) const { return a == b; }
length(const Element a)105         inline size_t length(const Element a) const { return size_rep; }
106 
107         // ----- Ring-wise operators
108         inline bool operator==(const Self_t& F) const { return _p == F._p; }
109         inline bool operator!=(const Self_t& F) const { return _p != F._p; }
110         inline Self_t& operator=(const Self_t& F)
111         {
112             F.assign(const_cast<Element&>(one),  F.one);
113             F.assign(const_cast<Element&>(zero), F.zero);
114             F.assign(const_cast<Element&>(mOne), F.mOne);
115             _p = F._p;
116             _negp = F._negp;
117             _invp = F._invp;
118             _lp= F._lp;
119             return *this;
120         }
121 
122 
123 
124         // ----- Initialisation
init(Element & x)125         inline Element& init (Element& x) const
126         {
127             return x = zero;
128         }
129 
130         template<typename T>
init(Element & r,T a)131         Element& init (Element& r, T a) const{
132             r = Caster<Element>(a);
133             return reduce(r);
134         }
135 
assign(Element & x,const Element & y)136         Element &assign (Element &x, const Element &y) const{
137             return x = y;
138         }
139 
140         // ----- Convert and reduce
convert(T & r,const Element & a)141         template<typename T> T& convert(T& r, const Element& a) const
142         { return r = static_cast<T>(a); }
143 
reduce(Element & x,const Element & y)144         Element& reduce (Element& x, const Element& y) const{
145             x=y;
146             return reduce(x);
147         }
148 
149         Element& reduce (Element& x) const ;
150 
151         // ----- Classic arithmetic
152         Element& mul(Element& r, const Element& a, const Element& b) const;
153 
div(Element & r,const Element & a,const Element & b)154         Element& div(Element& r, const Element& a, const Element& b) const{
155             return mulin(inv(r, a), b);
156         }
add(Element & r,const Element & a,const Element & b)157         Element& add(Element& r, const Element& a, const Element& b) const {
158             r = a + b;
159             if(r >= _p)
160                 r += _negp;
161             return r;
162         }
sub(Element & r,const Element & a,const Element & b)163         Element& sub(Element& r, const Element& a, const Element& b) const {
164             r = a - b;
165             if(r < 0)
166                 r += _p;
167             return r;
168         }
neg(Element & r,const Element & a)169         Element& neg(Element& r, const Element& a) const {
170             r = -a;
171             if(r < 0)
172                 r += _p;
173             return r;
174         }
inv(Element & x,const Element & y)175         Element& inv(Element& x, const Element& y) const{
176             invext(x,y,_p);
177             if (x<0) x += _p;
178             return x;
179         }
180 
mulin(Element & r,const Element & a)181         Element& mulin(Element& r, const Element& a) const {
182             return mul(r, r, a);
183         }
divin(Element & r,const Element & y)184         Element& divin(Element& r, const Element& y) const{
185             Element iy;
186             return mulin(r, inv(iy, y));
187         }
addin(Element & r,const Element & a)188         Element& addin(Element& r, const Element& a) const {
189             return add(r, r, a);
190         }
subin(Element & r,const Element & a)191         Element& subin(Element& r, const Element& a) const {
192             return sub(r, r, a);
193         }
negin(Element & r)194         Element& negin(Element& r) const {
195             return neg(r, r);
196         }
invin(Element & r)197         Element& invin(Element& r) const {
198             return inv(r, r);
199         }
200 
201         // -- axpy:   r <- a * x + y
202         // -- axpyin: r <- a * x + r
axpy(Element & r,const Element & a,const Element & x,const Element & y)203         Element& axpy  (Element& r, const Element& a, const Element& x, const Element& y) const {
204             Element tmp;
205             mul(tmp, a, x);
206             return add(r, tmp, y);
207         }
axpyin(Element & r,const Element & a,const Element & x)208         Element& axpyin(Element& r, const Element& a, const Element& x) const {
209             Element tmp(r);
210             return axpy(r, a, x, tmp);
211         }
212 
213         // -- axmy:   r <- a * x - y
214         // -- axmyin: r <- a * x - r
axmy(Element & r,const Element & a,const Element & x,const Element & y)215         Element& axmy  (Element& r, const Element& a, const Element& x, const Element& y) const {
216             Element tmp;
217             mul(tmp, a, x);
218             return sub(r, tmp, y);
219         }
axmyin(Element & r,const Element & a,const Element & x)220         Element& axmyin(Element& r, const Element& a, const Element& x) const {
221             return axmy(r, a, x, r);
222         }
223 
224         // -- maxpy:   r <- y - a * x
225         // -- maxpyin: r <- r - a * x
maxpy(Element & r,const Element & a,const Element & x,const Element & y)226         Element& maxpy  (Element& r, const Element& a, const Element& x, const Element& y) const {
227             Element tmp;
228             mul(tmp, a, x);
229             return sub(r, y, tmp);
230         }
maxpyin(Element & r,const Element & a,const Element & x)231         Element& maxpyin(Element& r, const Element& a, const Element& x) const {
232             return maxpy(r, a, x, r);
233         }
234 
235         // ----- Random generators
236         typedef ModularRandIter<Self_t> RandIter;
237         typedef GeneralRingNonZeroRandIter<Self_t> NonZeroRandIter;
random(const Random & g,Element & r)238         template< class Random > Element& random(const Random& g, Element& r) const { return init(r, g()); }
nonzerorandom(const Random & g,Element & a)239         template< class Random > Element& nonzerorandom(const Random& g, Element& a) const {
240             while (isZero(init(a, g())));
241             return a;
242         }
243 
244         // --- IO methods
245         std::istream& read (std::istream& s);
246         std::ostream& write(std::ostream& s) const;
247         std::istream& read (std::istream& s, Element& a) const;
248         std::ostream& write(std::ostream& s, const Element& a) const;
249 
250     protected:
251         Element _p = 0;
252         Element _invp = 0;
253         Element _negp = 0;
254         Residu_t _lp = 0;
255 
256     };
257 
258 } // Givaro
259 
260 #include "givaro/modular-extended.inl"
261 
262 #endif // __GIVARO_MODULAR_EXTENDED_H
263 /* -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
264 // vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
265