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