1 // ========================================================================== 2 // Copyright(c)'1994-2009 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 // file: gfqext.h 8 // Time-stamp: <12 Jun 15 16:28:43 Jean-Guillaume.Dumas@imag.fr> 9 // date: 2007 10 // version: 11 // author: Jean-Guillaume.Dumas 12 13 /*! @file gfqext.h 14 * @ingroup zpz 15 * @brief Arithmetic on GF(p^k), with p a prime number less than 2^15. 16 * Specialized for fast conversions to floating point numbers. 17 * Main difference in interface is init/convert. 18 * @bib 19 * - JG Dumas, <i>Q-adic Transform Revisited</i>, ISSAC 2008 20 * @warning k strictly greater than 1 21 */ 22 23 #ifndef __GIVARO_gfq_extension_H 24 #define __GIVARO_gfq_extension_H 25 26 #include "givaro/gfq.h" 27 #include "givaro/givpower.h" 28 29 #include <limits> 30 #include <deque> 31 32 namespace Givaro { 33 34 // init with preconditions, bad entry could segfault 35 template<class TT> class GFqExtFast; 36 37 // init defensive, bad entry are transformed, to the cost of slowdown 38 template<class TT> class GFqExt; 39 40 41 //! GFq Ext 42 template<class TT> class GFqExtFast : public GFqDom<TT> { 43 protected: 44 typedef typename Signed_Trait<TT>::unsigned_type UTT; 45 typedef TT Rep; 46 typedef GFqDom<TT> Father_t; 47 48 UTT _BITS; // the q-adic transform will be with q=2^_BITS 49 UTT _BASE; // this is 2^_BITS 50 UTT _MASK; // 2^_BITS - 1 51 UTT _maxn; // Worst case Maximal number of multiplications 52 // without reduction 53 UTT _degree;// exponent-1 54 UTT _pceil; // smallest such that characteristic<2^_pceil, 55 // used for fast for table indexing 56 UTT _MODOUT;// Largest accepted double for init 57 58 // Conversion tables from exponent to double z-adic representation 59 std::vector<double> _log2dbl; // Exponent to double 60 std::vector<UTT> _high2log; // Half double to exponent 61 std::vector<UTT> _low2log; // Other half in Time-Memory Trade-Off 62 63 64 public: 65 typedef GFqExtFast<TT> Self_t; 66 67 typedef Rep Element; 68 typedef Element* Element_ptr ; 69 typedef const Element* ConstElement_ptr; 70 71 72 typedef UTT Residu_t; 73 74 typedef Rep* Array; 75 typedef const Rep* constArray; 76 GFqExtFast()77 GFqExtFast(): 78 Father_t() 79 // , balanced(false) 80 { 81 } 82 83 // Extension MUST be a parameter of the constructor GFqExtFast(const UTT P,const UTT e)84 GFqExtFast( const UTT P, const UTT e) : 85 Father_t(P,e), 86 _BITS( std::numeric_limits< double >::digits/( (e<<1)-1) ), 87 _BASE(1 << _BITS), 88 _MASK( _BASE - 1), 89 _maxn( _BASE/(P-1)/(P-1)/e), 90 _degree( e-1 ) 91 // , balanced(false) 92 { 93 94 GIVARO_ASSERT(_maxn>0 , "[GFqExtFast]: field too large"); 95 builddoubletables(); 96 97 } 98 99 // Extension MUST be a parameter of the constructor 100 template<typename Vector> GFqExtFast(const UTT P,const UTT e,const Vector & modPoly)101 GFqExtFast( const UTT P, const UTT e, const Vector& modPoly) : 102 Father_t(P,e, modPoly), 103 _BITS( std::numeric_limits< double >::digits/( (e<<1)-1) ), 104 _BASE(1 << _BITS), 105 _MASK( _BASE - 1), 106 _maxn( _BASE/(P-1)/(P-1)/e), 107 _degree( e-1 ) 108 // , balanced(false) 109 { 110 GIVARO_ASSERT(_maxn>0 , "[GFqExtFast]: field too large"); 111 builddoubletables(); 112 113 } 114 ~GFqExtFast()115 virtual ~GFqExtFast() {}; 116 117 Self_t operator=( const Self_t& F) 118 { 119 this->zero = F.zero; 120 this->one = F.one; 121 this->mOne = F.mOne; 122 this->_characteristic = F._characteristic; 123 this->_dcharacteristic = F._dcharacteristic; 124 this->_exponent = F._exponent; 125 this->_q = F._q; 126 this->_qm1 = F._qm1; 127 this->_log2pol = F._log2pol; 128 this->_pol2log = F._pol2log; 129 this->_plus1 = F._plus1; 130 this->_BITS = F._BITS; 131 this->_BASE = F._BASE; 132 this->_MASK = F._MASK; 133 this->_maxn = F._maxn; 134 this->_degree = F._degree; 135 this->_log2dbl = F._log2dbl; 136 this->_low2log = F._low2log; 137 this->_high2log = F._high2log; 138 return *this; 139 } 140 GFqExtFast(const GFqDom<TT> & F)141 GFqExtFast( const GFqDom<TT>& F) : 142 Father_t(F), 143 _BITS( F._BITS ), _BASE( F._BASE ),_MASK( F._MASK ), 144 _maxn( F._maxn ),_degree( F._degree ), 145 _log2dbl ( F._log2dbl ), _low2log( F._low2log ), 146 _high2log (F._high2log ) 147 // , balanced(false) 148 { 149 } 150 151 // Accesses 152 bits()153 UTT bits() const 154 { return _BITS;} 155 base()156 UTT base() const 157 { return _BASE;} 158 mask()159 UTT mask() const 160 { return _MASK;} 161 maxdot()162 UTT maxdot() const 163 { return _maxn; } 164 characteristic(UTT & a)165 UTT& characteristic(UTT& a) const 166 { return a=this->_characteristic; } 167 characteristic()168 UTT characteristic() const 169 { return this->_characteristic; } 170 171 // const bool balanced; 172 173 // Init/Convert 174 175 using Father_t::init; 176 init(Rep & r,const unsigned long l)177 Rep& init( Rep& r, const unsigned long l) const 178 { 179 return Father_t::init(r,l); 180 } 181 init(Rep & pad,const double d)182 virtual Rep& init(Rep& pad, const double d) const 183 { 184 GIVARO_ASSERT(d>=0.0 , "[GFqExtFast]: init from a negative number"); 185 GIVARO_ASSERT(d<_MODOUT, "[GFqExtFast]: init from a too large number"); 186 // WARNING WARNING WARNING WARNING 187 // Precondition : 0 <= d < _MODOUT 188 // Can segfault if d is too large 189 // WARNING WARNING WARNING WARNING 190 uint64_t rll( static_cast<uint64_t>(d) ); 191 uint64_t tll( static_cast<uint64_t>(d/this->_dcharacteristic) ); 192 UTT prec(0); 193 UTT padl = (UTT)(rll - tll*this->_characteristic); 194 195 if (padl == this->_characteristic) { 196 padl -= this->_characteristic; 197 tll += 1; 198 } 199 for(size_t j = 0;j<_degree;++j) { 200 rll >>= _BITS; 201 tll >>= _BITS; 202 prec = (UTT)(rll-tll*this->_characteristic); 203 204 padl <<= _pceil; 205 padl ^= prec; 206 } 207 208 pad = (Rep)prec; 209 for(size_t j = 0;j<_degree;++j) { 210 rll >>= _BITS; 211 tll >>= _BITS; 212 prec = (UTT)(rll-tll*this->_characteristic); 213 214 pad <<= _pceil; 215 pad ^= prec; 216 } 217 218 padl = this->_low2log[(size_t)padl]; 219 pad = (Rep)this->_high2log[(size_t)pad]; 220 return this->addin(pad,(Rep)padl); 221 } 222 init(Rep & pad,const float d)223 virtual Rep& init(Rep& pad, const float d) const 224 { 225 return init(pad, (double)d); 226 } 227 228 using Father_t::convert; 229 convert(double & d,const Rep a)230 virtual double& convert(double& d, const Rep a) const 231 { 232 return d=_log2dbl[(size_t)a]; 233 } 234 convert(float & d,const Rep a)235 virtual float& convert(float& d, const Rep a) const 236 { 237 return d=(float)_log2dbl[(size_t)a]; 238 } 239 random(RandIter & g,Rep & r)240 template<class RandIter> Rep& random(RandIter& g, Rep& r) const 241 { 242 return init(r, static_cast<double>( (UTT)g() % _MODOUT)); 243 } 244 245 protected: 246 builddoubletables()247 void builddoubletables() 248 { 249 _log2dbl.resize(this->_log2pol.size()); 250 _pceil = 1; 251 for(unsigned long ppow = 2; ppow < this->_characteristic; ppow <<= 1,++_pceil) {} 252 unsigned long powersize = 1<<(_pceil * this->_exponent); 253 _MODOUT = UTT(powersize - 1); 254 _high2log.resize(powersize); 255 _low2log.resize(powersize); 256 257 258 259 typedef typename Father_t::Element ZElem; 260 Father_t Zp(this->_characteristic,1); 261 ZElem q,mq; Zp.init(q,2); 262 dom_power(q,q,_BITS,Zp); 263 Zp.neg(mq,q); 264 265 Element xkmu; 266 // This is X 267 xkmu = (Element)this->_pol2log[(size_t)this->_characteristic]; 268 // This is X^{e-1} 269 dom_power(xkmu,xkmu,this->_exponent-1,*this); 270 271 272 273 typedef Poly1FactorDom< Father_t, Dense > PolDom; 274 PolDom Pdom( Zp ); 275 276 typedef Poly1PadicDom< Father_t, Dense > PadicDom; 277 PadicDom PAD(Pdom); 278 279 Father_t Z2B(2,_BITS); 280 PolDom P2dom( Z2B ); 281 PadicDom P2AD( P2dom ); 282 283 std::vector<double>::iterator dblit = _log2dbl.begin(); 284 typename std::vector<UTT>::const_iterator polit = this->_log2pol.begin(); 285 286 287 for( ; polit != this->_log2pol.end(); ++polit, ++dblit) { 288 289 std::vector<double> vect; 290 std::deque<ZElem> low_ui; 291 292 P2AD.evaldirect( *dblit, 293 PAD.radixdirect( 294 vect, 295 (double)(*polit), 296 this->_exponent) 297 ); 298 299 unsigned long binpolit = static_cast<unsigned long>(vect[0]); 300 for(size_t i =1; i<this->_exponent; ++i) { 301 binpolit <<= _pceil; 302 binpolit += static_cast<unsigned long>(vect[i]); 303 } 304 305 ZElem tmp, prec, cour; Zp.init(cour); 306 Zp.init(prec, vect[0]); 307 for(size_t i = 1; i<this->_exponent; ++i) { 308 Zp.init(cour, vect[i]); 309 Zp.axpy(tmp, mq, cour, prec); 310 low_ui.push_back(tmp); 311 prec = cour; 312 } 313 314 PAD.eval(tmp , low_ui ); 315 _low2log[binpolit] = this->_pol2log[(size_t)tmp]; 316 317 low_ui.push_back(cour); 318 PAD.eval( tmp, low_ui); 319 Father_t::mul((Element&)_high2log[(size_t)binpolit],(Element)this->_pol2log[(size_t)tmp], xkmu); 320 } 321 } 322 }; 323 324 325 326 //! GFq Ext (other) 327 template<class TT> class GFqExt : public GFqExtFast<TT> { 328 protected: 329 typedef typename Signed_Trait<TT>::unsigned_type UTT; 330 typedef TT Rep; 331 typedef GFqDom<TT> Father_t; 332 typedef GFqExtFast<TT> DirectFather_t; 333 334 double _fMODOUT; 335 336 public: 337 typedef GFqExt<TT> Self_t; 338 339 typedef Rep Element; 340 typedef UTT Residu_t; 341 342 typedef Rep* Array; 343 typedef const Rep* constArray; 344 GFqExt()345 GFqExt(): DirectFather_t(), 346 _fMODOUT(static_cast<double>(this->_MODOUT)) {} 347 GFqExt(const UTT P,const UTT e)348 GFqExt( const UTT P, const UTT e) : 349 DirectFather_t(P,e), 350 _fMODOUT(static_cast<double>(this->_MODOUT)) {} 351 GFqExt(const GFqDom<TT> & F)352 GFqExt( const GFqDom<TT>& F) : 353 DirectFather_t(F), 354 _fMODOUT(static_cast<double>(this->_MODOUT)) {} 355 ~GFqExt()356 ~GFqExt() {} 357 358 using Father_t::init; 359 init(Rep & pad,const double d)360 virtual Rep& init(Rep& pad, const double d) const 361 { 362 // Defensive init 363 const double tmp(fmod(d,this->_fMODOUT)); 364 return DirectFather_t::init(pad, (tmp>0.0)?tmp:(tmp+_fMODOUT) ); 365 } 366 }; 367 368 } // namespace Givaro 369 370 #endif // __GIVARO_gfq_extension_H 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