1 /* linbox/algorithms/vector-fraction.h 2 * Copyright (C) 2004 David Pritchard 3 * 4 * Written by David Pritchard <daveagp@mit.edu> 5 * 6 * ========LICENCE======== 7 * This file is part of the library LinBox. 8 * 9 * LinBox is free software: you can redistribute it and/or modify 10 * it under the terms of the GNU Lesser General Public 11 * License as published by the Free Software Foundation; either 12 * version 2.1 of the License, or (at your option) any later version. 13 * 14 * This library is distributed in the hope that it will be useful, 15 * but WITHOUT ANY WARRANTY; without even the implied warranty of 16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 17 * Lesser General Public License for more details. 18 * 19 * You should have received a copy of the GNU Lesser General Public 20 * License along with this library; if not, write to the Free Software 21 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA 22 * ========LICENCE======== 23 */ 24 25 #ifndef __LINBOX_vector_fraction_H 26 #define __LINBOX_vector_fraction_H 27 28 #include "linbox/linbox-config.h" 29 #include "linbox/util/debug.h" 30 #include <stdio.h> 31 #include "linbox/vector/blas-vector.h" 32 #include "linbox/vector/vector-traits.h" 33 34 namespace LinBox 35 { 36 37 /** utility function to reduce a rational pair to lowest form */ 38 template<class Domain> reduceIn(Domain & D,std::pair<typename Domain::Element,typename Domain::Element> & frac)39 void reduceIn(Domain& D, std::pair<typename Domain::Element, typename Domain::Element> &frac) 40 { 41 linbox_check(!D.isZero(frac.second)); 42 43 if (D.isZero(frac.first)){ 44 D.init(frac.second, 1); 45 return; 46 } 47 48 typename Domain::Element gcd; 49 D.gcd(gcd, frac.first, frac.second); 50 D.divin(frac.first, gcd); 51 D.divin(frac.second, gcd); 52 } 53 54 /** utility function to gcd-in a vector of elements over a domain */ 55 //this could be replaced by a fancier version that combines elements linearly at random 56 template<class Domain, class Vector> vectorGcdIn(typename Domain::Element & result,Domain & D,Vector & v)57 void vectorGcdIn(typename Domain::Element& result, Domain& D, Vector& v) { 58 for (typename Vector::iterator i = v.begin(); i != v.end(); ++i) 59 D.gcdin(result, *i); 60 } 61 62 /** utility function, returns gcd of a vector of elements over a domain */ 63 // this could be replaced by a fancier version that combines elements linearly at random 64 template<class Domain, class Vector> vectorGcd(Domain & D,Vector & v)65 typename Domain::Element vectorGcd(Domain& D, Vector& v) { 66 typename Domain::Element result; 67 D.assign(result, D.zero); 68 vectorGcdIn(result, D, v); 69 return result; 70 } 71 72 73 /** 74 * \brief VectorFraction<Domain> is a vector of rational elements with common reduced denominator. 75 * Here Domain is a ring supporting the gcd, eg NTL_ZZ or Givaro::ZRing<Integer> 76 * For compatability with the return type of rationalSolver, it allows conversion from/to 77 * std::vector<std::pair<Domain::Element> >. 78 * All functions will return the fraction in reduced form, calling reduce() if necessary. 79 */ 80 81 template<class Domain> 82 class VectorFraction{ 83 public: 84 typedef typename Domain::Element Element; 85 typedef typename std::pair<Element, Element> Fraction; 86 typedef typename std::vector<Fraction> FVector; 87 typedef BlasVector<Domain> Vector; 88 89 Vector numer; 90 Element denom; 91 const Domain& _domain; 92 93 /** 94 * constructor from vector of rational numbers 95 * reduces individual pairs in-place first unless alreadyReduced=true 96 */ VectorFraction(const Domain & D,FVector & frac)97 VectorFraction(const Domain& D, FVector& frac 98 //,bool alreadyReduced = false 99 ) : 100 _domain(D) 101 { 102 bool alreadyReduced = false; 103 typename FVector::iterator i; 104 105 D.assign(denom, D.one); 106 if (!alreadyReduced) 107 for (i=frac.begin(); i!=frac.end(); ++i) 108 reduceIn(D, *i); 109 110 for (i=frac.begin(); i!=frac.end(); ++i) { 111 linbox_check(!D.isZero(i->second)); 112 D.lcmin(denom, i->second); 113 } 114 115 numer = Vector(frac.size()); 116 typename Vector::iterator j; 117 118 for (i=frac.begin(), j=numer.begin(); i!=frac.end(); ++i, ++j){ 119 D.mul(*j, denom, i->first); 120 D.divin(*j, i->second); 121 } 122 } 123 124 /** allocating constructor, returns [0, 0, ... 0]/1 */ VectorFraction(const Domain & D,size_t n)125 VectorFraction(const Domain& D, size_t n) : 126 numer(D,n) 127 ,_domain(D) 128 { 129 typename Vector::iterator j; 130 131 for (j=numer.begin(); j!=numer.end(); ++j) 132 D.assign(*j, D.zero); 133 } 134 135 /** copy constructor */ VectorFraction(const VectorFraction<Domain> & VF)136 VectorFraction(const VectorFraction<Domain>& VF) : 137 numer(VF._domain) 138 ,_domain(VF._domain) 139 { 140 copy(VF); 141 } 142 size()143 size_t size() const { return numer.size(); } 144 145 /** copy without construction */ copy(const VectorFraction<Domain> & VF)146 void copy(const VectorFraction<Domain>& VF) 147 { 148 //assumes _domain = VF._domain 149 denom = VF.denom; 150 numer.resize(VF.numer.size()); 151 typename Vector::iterator i; 152 typename Vector::const_iterator j; 153 154 for (i=numer.begin(), j=VF.numer.begin(); i!=numer.end(); ++i, ++j) 155 _domain.assign(*i, *j); 156 } 157 158 /** clear and resize without construction */ clearAndResize(size_t size)159 void clearAndResize(size_t size) 160 { 161 _domain.init(denom, (int64_t)1); 162 typename Vector::iterator i; 163 numer.resize(size); 164 for (i=numer.begin(); i!=numer.end(); ++i) 165 _domain.assign(*i, _domain.zero); 166 } 167 168 /** 169 * Replaces *this with a linear combination of *this and other 170 * such that the result has denominator == gcd(this->denom, other.denom) 171 * see Mulders+Storjohann : 'Certified Dense Linear System Solving' Lemma 2.1 172 * return value of true means that there was some improvement (ie denom was reduced) 173 */ combineSolution(const VectorFraction<Domain> & other)174 bool combineSolution(const VectorFraction<Domain>& other) 175 { 176 if (_domain.isDivisor(other.denom, denom)) return false; 177 if (_domain.isDivisor(denom, other.denom)) { 178 denom = other.denom; 179 numer = other.numer; 180 return true; 181 } 182 Element s, t, g; 183 _domain.gcd(g, s, t, denom, other.denom); 184 if (_domain.areEqual(g, denom)) ; //do nothing 185 else { 186 denom = g; 187 typename Vector::iterator it=numer.begin(); 188 typename Vector::const_iterator io=other.numer.begin(); 189 for (; it != numer.end(); ++it, ++io) { 190 _domain.mulin(*it, s); 191 _domain.axpyin(*it, t, *io); 192 } 193 return true; 194 } 195 return false; 196 } 197 198 /** 199 * Adds in-place to *this a multiple of other 200 * such that the result has gcd(denominator, denBound) == gcd(this->denom, other.denom, denBound) 201 * see Mulders+Storjohann : 'Certified Dense Linear System Solving' Lemma 6.1 202 * return value of true means that there was some improvement (ie gcd(denom, denBound) was reduced) 203 * g is gcd(denom, denBound), and is updated by this function when there is improvement 204 */ boundedCombineSolution(const VectorFraction<Domain> & other,const Element & denBound,Element & g)205 bool boundedCombineSolution(const VectorFraction<Domain>& other, const Element& denBound, Element& g) 206 { 207 208 //this means that new solution won't reduce g 209 if (_domain.isDivisor(other.denom, g)) return false; 210 211 //short-circuit in case the new solution is completely better than old one 212 Element _dtmp; 213 if (_domain.isDivisor(g, _domain.gcd(_dtmp, denBound, other.denom))) { 214 denom = other.denom; 215 numer = other.numer; 216 g = _dtmp; 217 return true; 218 } 219 220 Element A, g2, lincomb; 221 _domain.gcd(g, other.denom, g); //we know this reduces g 222 223 // find A s.t. gcd(denBound, denom + A*other.denom) = g 224 // strategy: pick random values of A <= d(y_0) 225 uint64_t tmp; 226 _domain.convert(tmp, denBound); 227 typename Domain::RandIter randiter(_domain, tmp); //seed omitted 228 // TODO: I don't think this random iterator has high-quality low order bits, which are needed 229 do { 230 randiter.random(A); 231 _domain.assign(lincomb, denom); 232 _domain.axpyin(lincomb, A, other.denom); 233 _domain.gcd(g2, lincomb, denBound); 234 } 235 while (!_domain.areEqual(g, g2)); 236 237 _domain.assign(denom, lincomb); 238 typename Vector::iterator it=numer.begin(); 239 typename Vector::const_iterator io=other.numer.begin(); 240 for (; it != numer.end(); ++it, ++io) 241 _domain.axpyin(*it, A, *io); 242 return true; 243 } 244 245 /** 246 * Adds in-place to *this a multiple of other to create an improved certificate ("z") 247 * n1/d1 = *this . b, n2/d2 = other . b in reduced form 248 * n1/d1 are updated so that new denominator is lcm(d1, d2); 249 * see Mulders+Storjohann : 'Certified Dense Linear System Solving' Lemma 6.2 250 * return value of true means that there was some improvement (ie d1 was increased) 251 */ combineCertificate(const VectorFraction<Domain> & other,Element & n1,Element & d1,const Element & n2,const Element d2)252 bool combineCertificate(const VectorFraction<Domain>& other, Element& n1, Element& d1, 253 const Element& n2, const Element d2) 254 { 255 //this means that new solution won't reduce g 256 if (_domain.isDivisor(d1, d2)) return false; 257 258 //short-circuit in case the new solution is completely better than old one 259 if (_domain.isDivisor(d2, d1)) { 260 copy(other); 261 n1 = n2; 262 d1 = d2; 263 return true; 264 } 265 266 Element A, g, l, n1d2_g, n2d1_g, lincomb, g2, tmpe; 267 268 _domain.gcd(g, d1, d2); //compute gcd 269 _domain.mul(l, d1, d2); 270 _domain.divin(l, g); //compute lcm 271 272 _domain.div(n1d2_g, d2, g); 273 _domain.mulin(n1d2_g, n1); //compute n1.d2/g 274 _domain.div(n2d1_g, d1, g); 275 _domain.mulin(n2d1_g, n2); //compute n2.d1/g 276 277 // find A s.t. gcd(denBound, denom + A*other.denom) = g 278 // strategy: pick random values of A <= lcm(d(denom), d(other.denom)) 279 uint64_t tmp; 280 _domain.mul(tmpe, denom, other.denom); 281 _domain.convert(tmp, tmpe); 282 typename Domain::RandIter randiter(_domain, tmp); //seed omitted 283 // TODO: I don't think this random iterator has high-quality low order bits, which are needed 284 do { 285 randiter.random(A); 286 _domain.assign(lincomb, n1d2_g); 287 _domain.axpyin(lincomb, A, n2d1_g); 288 _domain.gcd(g2, lincomb, l); 289 } 290 while (!_domain.areEqual(_domain.one, g2)); 291 292 this->axpyin(A, other); 293 _domain.lcmin(d1, d2); 294 295 return true; 296 } 297 298 /** 299 * this += a * x. performs a rational axpy with an integer multiplier 300 * returns (*this) 301 */ axpyin(Element & a,const VectorFraction<Domain> & x)302 VectorFraction<Domain>& axpyin(Element& a, const VectorFraction<Domain>& x) 303 { 304 Element a_prime, gcd_a_xdenom, xdenom_prime; 305 _domain.gcd(gcd_a_xdenom, a, x.denom); 306 _domain.div(a_prime, a, gcd_a_xdenom); 307 _domain.div(xdenom_prime, x.denom, gcd_a_xdenom); 308 309 Element cdf; //common denominator factor; multiply both sides by this and divide at end 310 _domain.gcd(cdf, denom, xdenom_prime); 311 _domain.divin(denom, cdf); 312 _domain.divin(xdenom_prime, cdf); 313 314 // we perform numer[i] = xdenom_prime * numer[i] + a_prime * denom * x.denom[i] 315 // so multiply denom into a_prime and save a multiplication on each entry 316 _domain.mulin(a_prime, denom); 317 318 typename Vector::iterator i = this->numer.begin(); 319 typename Vector::const_iterator j = x.numer.begin(); 320 for (; i != this->numer.end(); ++i, ++j) { 321 _domain.mulin(*i, xdenom_prime); 322 _domain.axpyin(*i, a_prime, *j); 323 } 324 325 _domain.mulin(denom, cdf); 326 _domain.mulin(denom, xdenom_prime); 327 simplify(); 328 return *this; 329 } 330 331 /** write to a stream */ write(std::ostream & os)332 std::ostream& write(std::ostream& os) const 333 { 334 os << "["; 335 for (typename Vector::const_iterator it=numer.begin(); it != numer.end(); ++it) { 336 if (it != numer.begin()) os << " "; 337 os << *it; 338 } 339 return os << "]/" << denom; 340 } 341 342 std::ostream& operator<<(std::ostream& os) const { 343 return write(os); 344 } 345 346 /** convert to 'answer' type of lifting container */ toFVector(FVector & result)347 FVector& toFVector(FVector& result) const 348 { 349 linbox_check(numer.size()==result.size()); 350 typename Vector::const_iterator it=numer.begin(); 351 typename FVector::iterator ir=result.begin(); 352 for (; it != numer.end(); ++it, ++ir) { 353 _domain.assign(ir->first, *it); 354 _domain.assign(ir->second, denom); 355 } 356 return result; 357 } 358 359 /** reduces to simplest form, returns (*this) */ simplify()360 VectorFraction<Domain>& simplify() 361 { 362 typename Vector::iterator i; 363 Element gcd; 364 _domain.init(gcd, denom); 365 vectorGcdIn(gcd, _domain, numer); 366 367 _domain.divin(denom, gcd); 368 for (i=numer.begin(); i!=numer.end(); ++i) 369 _domain.divin(*i, gcd); 370 return (*this); 371 } 372 }; 373 374 } 375 376 #endif //__LINBOX_vector_fraction_H 377 378 // Local Variables: 379 // mode: C++ 380 // tab-width: 4 381 // indent-tabs-mode: nil 382 // c-basic-offset: 4 383 // End: 384 // vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s 385