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 // Time-stamp: <18 Feb 11 16:04:19 Jean-Guillaume.Dumas@imag.fr> 8 // Author: J-G. Dumas 9 // Description: Pieces of polynomials as defined in 10 // [Arne Storjohann, High-Order Lifting 11 // ISSAC'2002, pp. 246-254, ACM Press, July 2002]. 12 // =============================================================== 13 #ifndef __GIVARO_trunc_domain_H 14 #define __GIVARO_trunc_domain_H 15 #include <givaro/givpoly1dense.h> 16 #ifndef __PATHCC__ 17 #include <utility> 18 #endif 19 20 21 namespace Givaro { 22 23 template <class Domain> 24 class TruncDom : public Poly1Dom<Domain,Dense> { 25 public : 26 // -- Self_t 27 typedef TruncDom<Domain> Self_t; 28 // -- Father_t 29 typedef Poly1Dom<Domain,Dense> Father_t; 30 // -- Exported types 31 typedef Domain Domain_t; 32 typedef typename Domain::Element Type_t; 33 typedef typename Father_t::Storage_t Polynomial_t; 34 typedef std::pair<Polynomial_t, Degree> Storage_t; 35 typedef Storage_t Rep; 36 typedef Storage_t Element; 37 38 39 Storage_t zero, one,mOne; 40 Father_t(d,X)41 TruncDom (const Domain& d, const Indeter& X = Indeter() ) : Father_t(d,X) { 42 this->assign(zero,Father_t::zero); 43 this->assign(one,Father_t::one); 44 this->assign(mOne,Father_t::mOne); 45 } TruncDom(const Self_t & t)46 TruncDom (const Self_t& t) : Father_t(static_cast<const Father_t&>(t)) { 47 this->assign(zero,Father_t::zero); 48 this->assign(one,Father_t::one); 49 this->assign(mOne,Father_t::mOne); 50 } TruncDom(const Father_t & t)51 TruncDom (const Father_t& t) : Father_t(t) { 52 this->assign(zero,Father_t::zero); 53 this->assign(one,Father_t::one); 54 this->assign(mOne,Father_t::mOne); 55 } 56 getpoldomain()57 const Father_t& getpoldomain() const 58 { 59 return static_cast<const Father_t&>(*this); 60 } 61 62 init(Rep & p)63 Rep& init(Rep& p) const 64 { Father_t::init(p.first); p.second=0; return p; } 65 66 template<class XXX> init(Rep & p,const XXX & cste)67 Rep& init(Rep& p, const XXX &cste ) const 68 { 69 Father_t::init(p.first,cste); p.second=0; return p; 70 } 71 72 // -- For polynomial = lcoeff X^deg 73 template<class XXX> init(Rep & p,const Degree deg,const XXX & lcoeff)74 Rep& init (Rep& p, const Degree deg , const XXX& lcoeff) const 75 { 76 Father_t::init(p.first,Degree(0),lcoeff); p.second=deg; return p; 77 } 78 convert(Polynomial_t & r,const Rep & P)79 Polynomial_t& convert (Polynomial_t& r, const Rep& P) const 80 { 81 Father_t::assign(r, P.first); 82 Polynomial_t mon; 83 Father_t::init(mon, P.second); 84 return Father_t::mulin(r, mon); 85 } 86 87 template<class XXX> convert(XXX & r,const Rep & P)88 XXX& convert (XXX& r, const Rep& P) const 89 { 90 Polynomial_t eP; Father_t::init(eP); 91 this->convert(eP, P); 92 return Father_t::convert(r, eP); 93 } 94 95 // F.assign(P[deg], lcoeff); assign(Rep & p,const Degree deg,const Type_t & lcoeff)96 Rep& assign (Rep& p, const Degree deg , const Type_t& lcoeff) const 97 { 98 Father_t::assign(p.first,Degree(0),lcoeff); 99 Father_t::setdegree(p.first); 100 p.second=deg; 101 return setval(p); 102 } 103 104 // -- Assign polynomial with field value : F.assign(p[0],cste) assign(Rep & p,const Type_t & cste)105 Rep& assign(Rep& p, const Type_t &cste ) const 106 { 107 return assign(p, Degree(0), cste); 108 } 109 // -- Assignment p = q assign(Rep & p,const Rep & q)110 Rep& assign( Rep& p, const Rep& q) const 111 { 112 Father_t::assign(p.first, q.first); 113 p.second = q.second; 114 return p; 115 } 116 assign(Rep & p,const Polynomial_t & r)117 Rep& assign(Rep& p, const Polynomial_t& r ) const 118 { 119 Father_t::assign(p.first,r); p.second=0; return setval(p); 120 } 121 assign(Rep & p,const Polynomial_t & r,const Degree v,const Degree d)122 Rep& assign(Rep& p, const Polynomial_t& r, const Degree v, const Degree d) const 123 { 124 Father_t::assign(p.first,r); p.second=0; return truncin(p,v,d); 125 } 126 127 mulin(Rep & p,const Degree & s)128 Rep& mulin(Rep& p, const Degree& s) const 129 { 130 p.second += s; return p; 131 } 132 divin(Rep & p,const Degree & s)133 Rep& divin(Rep& p, const Degree& s) const 134 { 135 p.second -= s; p.second=(p.second<0?0:p.second); return p; 136 } 137 138 // -- Compute the degree of P setdegree(Rep & P)139 Rep& setdegree( Rep& P ) const 140 { 141 Father_t::setdegree(P.first); 142 if (P.first.size() <= 0) P.second=0; 143 return P; 144 } 145 146 // -- Compute the valuation of P setval(Rep & P)147 Rep& setval( Rep& P ) const 148 { 149 setdegree(P); 150 if (P.first.size() <= 0) return P; 151 typename Polynomial_t::iterator it = P.first.begin(); 152 if (! this->_domain.isZero(*it)) return P; 153 for(++it,++P.second; it != P.first.end(); ++it,++P.second) { 154 if (! this->_domain.isZero(*it)) { 155 P.first.erase(P.first.begin(),it); 156 return P; 157 } 158 } 159 P.first.resize(0); P.second=0; 160 return P; 161 } 162 163 // -- Returns the degree of polynomial degree(Degree & d,const Rep & P)164 Degree& degree(Degree& d, const Rep& P) const 165 { 166 Father_t::degree(d, P.first); 167 return d+=P.second; 168 } 169 170 // -- Returns the valuation of polynomial val(Degree & d,const Rep & P)171 Degree& val(Degree& d, const Rep& P) const 172 { 173 this->setval(const_cast<Rep&>(P)); 174 return d=P.second; 175 } 176 177 178 // -- Comparaison operator isZero(const Rep & P)179 int isZero ( const Rep& P ) const 180 { return Father_t::isZero(P.first); } isOne(const Rep & P)181 int isOne ( const Rep& P ) const 182 { 183 Degree vP;val(vP,P); 184 return ((vP == 0) && Father_t::isOne(P.first)); 185 } isMOne(const Rep & P)186 int isMOne ( const Rep& P ) const 187 { 188 Degree vP;val(vP,P); 189 return ((vP == 0) && Father_t::isMOne(P.first)); 190 } 191 192 areEqual(const Rep & P,const Rep & Q)193 int areEqual ( const Rep& P, const Rep& Q ) const 194 { 195 Degree vP;val(vP,P); 196 Degree vQ;val(vQ,Q); 197 return (vP == vQ) && Father_t::areEqual(P.first,Q.first); 198 } 199 areNEqual(const Rep & P,const Rep & Q)200 int areNEqual( const Rep& P, const Rep& Q ) const 201 { 202 Degree vP;val(vP,P); 203 Degree vQ;val(vQ,Q); 204 return (vP != vQ) || Father_t::areNEqual(P.first,Q.first); 205 } 206 207 208 Rep& shift(Rep& p, const Degree& s) const; 209 Rep& truncin(Rep& p, const Degree& v, const Degree& d) const; trunc(Rep & p,const Rep & R,const Degree & v,const Degree & d)210 Rep& trunc(Rep& p, const Rep& R, const Degree& v, const Degree& d) const 211 { 212 return this->truncin(this->assign(p,R),v,d); 213 } 214 215 // I/O read(std::istream & i)216 std::istream& read ( std::istream& i ) { 217 char tmp, t[5]; 218 return Father_t::read(i>> tmp)>> t ; 219 } write(std::ostream & o)220 std::ostream& write( std::ostream& o ) const 221 { 222 return Father_t::write(o << '[') << "]_i^j"; 223 } read(std::istream & i,Rep & n)224 std::istream& read ( std::istream& i, Rep& n) const 225 { 226 char tmp; 227 return Father_t::read(i>>tmp,n.first)>>tmp>>tmp>>tmp>>tmp>> n.second; 228 } write(std::ostream & o,const Rep & n)229 std::ostream& write( std::ostream& o, const Rep& n) const 230 { 231 232 return Father_t::write(o<<'(',n.first)<<")*" << this->_x << '^' << n.second; 233 } 234 expand(Rep & P,const Degree & d)235 Rep& expand(Rep& P, const Degree& d) const 236 { 237 Degree vP; val(vP, P); 238 if (vP > d) { 239 P.first.insert(P.first.begin(),(size_t)value(vP-d),this->_domain.zero); 240 P.second = d; 241 } 242 return P; 243 } 244 245 246 // -- Arithmetics operators 247 Rep& addin ( Rep& R, const Rep& P) const; add(Rep & res,const Rep & u,const Rep & v)248 Rep& add ( Rep& res, const Rep& u, const Rep& v ) const 249 { 250 assign(res,u); 251 return addin(res,v); 252 } 253 254 Rep& addin ( Rep& R, const Rep& P, const Degree& v, const Degree& d) const; add(Rep & res,const Rep & u,const Rep & v,const Degree & Val,const Degree & deg)255 Rep& add ( Rep& res, const Rep& u, const Rep& v, const Degree& Val, const Degree& deg) const 256 { 257 assign(res,u); 258 return addin(res,v,Val,deg); 259 } 260 neg(Rep & R,const Rep & P)261 Rep& neg(Rep& R, const Rep& P) const 262 { 263 Father_t::neg(R.first,P.first); 264 R.second = P.second; 265 return R; 266 } negin(Rep & R)267 Rep& negin(Rep& R) const 268 { 269 Father_t::negin(R.first); 270 return R; 271 } 272 sub(Rep & res,const Rep & u,const Rep & v)273 Rep& sub ( Rep& res, const Rep& u, const Rep& v ) const 274 { 275 assign(res,u); 276 return this->subin(res,v); 277 } 278 Rep& subin ( Rep& R, const Rep& P) const ; sub(Rep & R,const Rep & P,const Rep & Q,const Degree & v,const Degree & d)279 Rep& sub ( Rep& R, const Rep& P, const Rep& Q, const Degree& v, const Degree& d) const 280 { 281 return this->addin(this->neg(R,Q),P,v,d); 282 283 } 284 subin(Rep & R,const Rep & P,const Degree & v,const Degree & d)285 Rep& subin ( Rep& R, const Rep& P, const Degree& v, const Degree& d) const 286 { 287 return this->negin(this->addin(this->negin(R),P,v,d)); 288 } mul(Rep & res,const Rep & u,const Rep & v)289 Rep& mul ( Rep& res, const Rep& u, const Rep& v ) const 290 { 291 Father_t::mul(res.first,u.first,v.first); 292 res.second = u.second+v.second; 293 return res; 294 } mulin(Rep & P,const Rep & Q)295 Rep& mulin ( Rep& P, const Rep& Q ) const 296 { 297 Father_t::mulin(P.first,Q.first); 298 P.second += Q.second; 299 return P; 300 } 301 302 303 Rep& mul( Rep& r, const Rep& u, const Rep& v, const Degree& Val, const Degree& deg) const; mulin(Rep & r,const Rep & v,const Degree & Val,const Degree & deg)304 Rep& mulin( Rep& r, const Rep& v, const Degree& Val, const Degree& deg) const 305 { 306 Rep tmp(r); 307 return mul(r,tmp,v,Val,deg); 308 } 309 axpy(Rep & r,const Rep & a,const Rep & x,const Rep & y)310 Rep& axpy (Rep& r, const Rep& a, const Rep& x, const Rep& y) const 311 { 312 return this->addin( this->mul(r,a,x), y); 313 } axpyin(Rep & r,const Rep & a,const Rep & x)314 Rep& axpyin(Rep& r, const Rep& a, const Rep& x) const 315 { 316 Rep tmp; this->init(tmp); 317 return this->addin(r, this->mul(tmp,a,x)); 318 } axpy(Rep & r,const Rep & a,const Rep & x,const Rep & y,const Degree & Val,const Degree & deg)319 Rep& axpy (Rep& r, const Rep& a, const Rep& x, const Rep& y, const Degree& Val, const Degree& deg) const 320 { 321 return this->addin(this->mul(r,a,x,Val,deg), y, Val, deg); 322 } axpyin(Rep & r,const Rep & a,const Rep & x,const Degree & Val,const Degree & deg)323 Rep& axpyin (Rep& r, const Rep& a, const Rep& x, const Degree& Val, const Degree& deg) const 324 { 325 Rep tmp; this->init(tmp); 326 return this->addin(r, this->mul(tmp,a,x,Val,deg), Val, deg); 327 } 328 axmy(Rep & r,const Rep & a,const Rep & x,const Rep & y)329 Rep& axmy (Rep& r, const Rep& a, const Rep& x, const Rep& y) const 330 { 331 return this->subin( this->mul(r,a,x), y); 332 } axmyin(Rep & r,const Rep & a,const Rep & x)333 Rep& axmyin(Rep& r, const Rep& a, const Rep& x) const 334 { 335 return this->negin(this->maxpyin(r,a,x)); 336 } axmy(Rep & r,const Rep & a,const Rep & x,const Rep & y,const Degree & Val,const Degree & deg)337 Rep& axmy (Rep& r, const Rep& a, const Rep& x, const Rep& y, const Degree& Val, const Degree& deg) const 338 { 339 return this->subin(this->mul(r,a,x,Val,deg), y, Val, deg); 340 } axmyin(Rep & r,const Rep & a,const Rep & x,const Degree & Val,const Degree & deg)341 Rep& axmyin (Rep& r, const Rep& a, const Rep& x, const Degree& Val, const Degree& deg) const 342 { 343 return this->negin(this->maxpyin(r,a,x,Val,deg)); 344 } 345 maxpy(Rep & r,const Rep & a,const Rep & x,const Rep & y)346 Rep& maxpy (Rep& r, const Rep& a, const Rep& x, const Rep& y) const 347 { 348 return this->addin( this->negin(this->mul(r,a,x)), y); 349 } maxpyin(Rep & r,const Rep & a,const Rep & x)350 Rep& maxpyin(Rep& r, const Rep& a, const Rep& x) const 351 { 352 Rep tmp; 353 return this->subin(r, this->mul(tmp,a,x)); 354 } maxpy(Rep & r,const Rep & a,const Rep & x,const Rep & y,const Degree & Val,const Degree & deg)355 Rep& maxpy (Rep& r, const Rep& a, const Rep& x, const Rep& y, const Degree& Val, const Degree& deg) const 356 { 357 return this->addin( this->negin(this->mul(r,a,x,Val,deg)), y, Val, deg); 358 } maxpyin(Rep & r,const Rep & a,const Rep & x,const Degree & Val,const Degree & deg)359 Rep& maxpyin (Rep& r, const Rep& a, const Rep& x, const Degree& Val, const Degree& deg) const 360 { 361 Rep tmp; 362 return this->subin(r, this->mul(tmp,a,x,Val,deg), Val, deg); 363 } 364 365 // -- Random generators 366 // -- Random dense polynomial of degree 0 random(RandIter & g,Rep & r)367 template< class RandIter > Rep& random(RandIter& g, Rep& r) const 368 { 369 Father_t::random(g,r.first); 370 r.second = rand(); 371 return r; 372 } 373 374 // -- Random dense polynomial of size s random(RandIter & g,Rep & r,long s)375 template< class RandIter > Rep& random(RandIter& g, Rep& r, long s) const 376 { 377 Father_t::random(g,r.first,s); 378 r.second = rand() % s; 379 return r; 380 } 381 // -- Random dense polynomial of degree d random(RandIter & g,Rep & r,Degree s)382 template< class RandIter > Rep& random(RandIter& g, Rep& r, Degree s) const 383 { 384 Father_t::random(g,r.first,s); 385 r.second = rand() % s.value(); 386 return r; 387 } 388 random(GivRandom & g,Rep & r,Degree s)389 Rep& random(GivRandom& g, Rep& r, Degree s) const 390 { 391 Father_t::random(g,r.first,s); 392 r.second = (Degree)(long)((unsigned long)g() % (unsigned long)((s.value()<<1)|1)); 393 return r; 394 } 395 // -- Random dense polynomial with same size as b. random(RandIter & g,Rep & r,const Rep & b)396 template< class RandIter > Rep& random(RandIter& g, Rep& r, const Rep& b) const 397 { 398 Father_t::random(g,r.first,b); 399 r.second = b.second; 400 return r; 401 } 402 nonzerorandom(RandIter & g,Rep & r)403 template< class RandIter > Rep& nonzerorandom(RandIter& g, Rep& r) const{ 404 Father_t::nonzerorandom(g,r.first); 405 r.second = rand(); 406 return r; 407 } nonzerorandom(RandIter & g,Rep & r,long s)408 template< class RandIter > Rep& nonzerorandom(RandIter& g, Rep& r, long s) const 409 { 410 Father_t::nonzerorandom(g,r.first,s); 411 r.second = rand() % s; 412 return r; 413 } nonzerorandom(RandIter & g,Rep & r,Degree s)414 template< class RandIter > Rep& nonzerorandom(RandIter& g, Rep& r, Degree s) const 415 { 416 Father_t::nonzerorandom(g,r.first,s); 417 r.second = rand() % s.value(); 418 return r; 419 } nonzerorandom(RandIter & g,Rep & r,const Rep & b)420 template< class RandIter > Rep& nonzerorandom(RandIter& g, Rep& r, const Rep& b) const 421 { 422 Father_t::random(g,r.first,b); 423 r.second = b.second; 424 return r; 425 } 426 427 428 }; 429 430 } // Givaro 431 432 #include "givaro/givtruncdomain.inl" 433 434 #endif 435 /* -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ 436 // vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s 437