1 /* linbox/algorithms/coppersmith.h 2 * evolved from block-wiedemann.h by George Yuhasz 3 * 4 * ========LICENCE======== 5 * This file is part of the library LinBox. 6 * 7 * LinBox is free software: you can redistribute it and/or modify 8 * it under the terms of the GNU Lesser General Public 9 * License as published by the Free Software Foundation; either 10 * version 2.1 of the License, or (at your option) any later version. 11 * 12 * This library is distributed in the hope that it will be useful, 13 * but WITHOUT ANY WARRANTY; without even the implied warranty of 14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 15 * Lesser General Public License for more details. 16 * 17 * You should have received a copy of the GNU Lesser General Public 18 * License along with this library; if not, write to the Free Software 19 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA 20 * ========LICENCE======== 21 */ 22 23 24 #ifndef __LINBOX_coppersmith_H 25 #define __LINBOX_coppersmith_H 26 27 #include <vector> 28 #include <numeric> 29 #include <algorithm> 30 #include <iostream> 31 #include <givaro/givpoly1crt.h> 32 33 34 #include "linbox/integer.h" 35 #include "linbox/util/commentator.h" 36 #include "linbox/algorithms/blackbox-block-container.h" 37 #include "linbox/algorithms/block-coppersmith-domain.h" 38 #include "linbox/solutions/det.h" 39 40 #include "linbox/util/error.h" 41 #include "linbox/util/debug.h" 42 43 namespace LinBox 44 { 45 46 template <class _Domain> 47 class CoppersmithSolver{ 48 49 public: 50 typedef _Domain Domain; 51 typedef typename Domain::Field Field; 52 typedef typename Domain::Element Element; 53 typedef typename Domain::OwnMatrix Block; 54 typedef typename Domain::Matrix Sub; 55 domain()56 inline const Domain & domain() const { return *_MD; } field()57 inline const Field & field() const { return domain().field(); } 58 protected: 59 const Domain *_MD; 60 size_t blocking; 61 62 public: 63 CoppersmithSolver(const Domain &MD, size_t blocking_ = 0) : 64 _MD(&MD), blocking(blocking_) 65 {} 66 67 68 template <class Vector, class Blackbox> solveNonSingular(Vector & x,const Blackbox & B,const Vector & y)69 Vector &solveNonSingular (Vector &x, const Blackbox &B, const Vector &y) const 70 { 71 commentator().start ("Coppersmith solveNonSingular", "solveNonSingular"); 72 #if 1 73 std::ostream& report = commentator().report(Commentator::LEVEL_IMPORTANT, INTERNAL_DESCRIPTION); 74 #endif 75 76 //Set up the projection matrices and their dimensions 77 size_t d = B.coldim(); 78 size_t r,c; 79 integer tmp = uint64_t(d); 80 81 //Set the blocking size, Using Pascal Giorgi's convention 82 if(blocking==0){ 83 r=tmp.bitsize()-1; 84 c=tmp.bitsize()-1; 85 }else 86 r=c=blocking; 87 88 //Create the block 89 Block U(field(),r,d); 90 Block W(field(),d,c-1); 91 Block V(field(),d,c); 92 93 //Pick random entries for U and W. W will become the last c-1 columns of V 94 95 U.random(); 96 W.random(); 97 98 //Multiply W by B on the left and place it in the last c-1 columns of V 99 Sub V2(V,0,1,d,c-1); 100 domain().mul(V2,B,W); 101 102 //Make the first column of V a copy of the right side of the system, y 103 for(size_t i=0; i<d; i++) 104 V.setEntry(i,0,y[i]); 105 106 //Create the sequence container and its iterator that will compute the projection 107 BlackboxBlockContainer<Field, Blackbox > blockseq(&B,field(),U,V); 108 109 //Get the generator of the projection using the Coppersmith algorithm (slightly modified by Yuhasz) 110 BlockCoppersmithDomain<Domain, BlackboxBlockContainer<Field, Blackbox> > BCD(domain(), &blockseq,d); 111 std::vector<Block> gen; 112 std::vector<size_t> deg; 113 deg = BCD.right_minpoly(gen); 114 report << "Size of blockseq " << blockseq.size() << std::endl; 115 report << "Size of gen " << gen.size() << std::endl; 116 for(size_t i = 0; i < gen[0].coldim(); i++) 117 report << "Column " << i << " has degree " << deg[i] << std::endl; 118 119 //Reconstruct the solution 120 //Pick a column of the generator with a nonzero element in the first row of the constant coefficient 121 size_t idx = 0; 122 if(field().isZero(gen[0].getEntry(0,0))){ 123 size_t i = 1; 124 while(i<c && field().isZero(gen[0].getEntry(0,i))) 125 i++; 126 if(i==c) 127 throw LinboxError(" block minpoly: matrix seems to be singular - abort"); 128 else 129 idx=i; 130 } 131 132 //from 1 to the degree of the index column, multiply A^(i-1)V times the idx column of the generator coefficient x^i 133 //Accumulate these results in xm 134 size_t mu = deg[idx]; 135 Block BVo(V); 136 Block BVe(field(),d,c); 137 Block xm(field(),d,1); 138 bool odd = true; 139 for(size_t i = 1; i < mu+1; i++){ 140 Sub gencol(gen[i],0,idx,c,1); // BB changed d,1 to c,1 141 Block BVgencol(field(),d,1); 142 if(odd){ 143 domain().mul(BVgencol,BVo,gencol); 144 domain().addin(xm, BVgencol); 145 domain().mul(BVe,B,BVo); 146 odd=false; 147 } 148 else{ 149 domain().mul(BVgencol,BVe,gencol); 150 domain().addin(xm, BVgencol); 151 domain().mul(BVo,B,BVe); 152 odd=true; 153 } 154 155 } 156 157 //For the constant coefficient, loop over the elements in the idx column except the first row 158 //Multiply the corresponding column of W (the last c-1 columns of V before application of B) by the generator element 159 //Accumulate the results in xm 160 for(size_t i = 1; i < c; i++){ 161 Sub Wcol(W,0,i-1,d,1); 162 Block Wcolgen0(field(),d,1); 163 domain().mul(Wcolgen0, Wcol, gen[0].getEntry(i,idx)); 164 domain().addin(xm,Wcolgen0); 165 } 166 167 //Multiply xm by -1(move to the correct side of the equation) and divide the the 0,idx entry of the generator constant 168 Element gen0inv; 169 field().inv(gen0inv,gen[0].getEntry(0,idx)); 170 field().negin(gen0inv); 171 domain().mulin(xm, gen0inv); 172 173 #if 0 174 //Test to see if the answer works with U 175 Block Bxm(field(),d,1), UBxm(field(),r,1), Uycol(field(), r,1); 176 Sub ycol(V,0,0,d,1); 177 domain().mul(Uycol, U, ycol); 178 domain().mul(Bxm, B, xm); 179 domain().mul(UBxm, U, Bxm); 180 181 if(domain().areEqual(UBxm, Uycol)) 182 report << "The solution matches when projected by U" << endl; 183 else 184 report << "The solution does not match when projected by U" << endl; 185 #endif 186 187 188 //Copy xm into x (Change type from 1 column matrix to Vector) 189 for(size_t i =0; i<d; i++) 190 x[i]=xm.getEntry(i,0); 191 192 commentator().stop ("done", NULL, "solveNonSingular"); 193 return x; 194 } 195 196 197 198 199 }; // end of class CoppersmithSolver 200 201 template <class _Domain> 202 class CoppersmithRank{ 203 204 public: 205 typedef _Domain Domain; 206 typedef typename Domain::Field Field; 207 typedef typename Domain::Element Element; 208 typedef typename Domain::Matrix Block; 209 typedef typename Domain::Submatrix Sub; 210 typedef typename Field::RandIter Random; 211 domain()212 inline const Domain & domain() const { return *_MD; } field()213 inline const Field & field() const { return domain().field(); } 214 protected: 215 const Domain *_MD; 216 Random iter; 217 size_t blocking; 218 219 //Compute the determinant of a polynomial matrix at the given set of evaluation points 220 //Store the results in the vector dets. EvalPolyMat(std::vector<Element> & dets,std::vector<Element> & values,std::vector<Block> & mat)221 void EvalPolyMat(std::vector<Element> &dets, std::vector<Element> &values, std::vector<Block> & mat) const { 222 223 size_t deg = mat.size() -1; 224 size_t numv = values.size(); 225 //Compute the determinant of the evaluation at values[i] for each i 226 for(size_t i = 0; i<numv; i++){ 227 //copy the highest matrix coefficient 228 Block evalmat(field(), mat[0].rowdim(), mat[0].coldim()); 229 domain().copy(evalmat,mat[deg]); 230 //Evaluate using a horner style evaluation 231 typename std::vector<Block>::reverse_iterator addit = mat.rbegin(); 232 addit++; 233 for(addit; addit != mat.rend(); addit++){ 234 domain().mulin(evalmat,values[i]); 235 domain().addin(evalmat,*addit); 236 }//end loop computing horner evaluation 237 //Compute the determinant of the evaluation and store it in dets[i] 238 dets[i] = det(dets[i],evalmat); 239 }//end loop over evaluation points 240 }//end evaluation of polynominal matrix determinant 241 242 public: 243 CoppersmithRank(const Domain &MD, size_t blocking_ = 0) : 244 _MD(&MD), blocking(blocking_), iter(MD.field()) 245 {} 246 247 248 template <class Blackbox> rank(const Blackbox & B)249 size_t rank (const Blackbox &B) const 250 { 251 commentator().start ("Coppersmith rank", "rank"); 252 #if 1 253 std::ostream& report = commentator().report(Commentator::LEVEL_IMPORTANT, INTERNAL_DESCRIPTION); 254 #endif 255 256 //Set up the projection matrices and their dimensions 257 size_t d = B.coldim(); 258 size_t r,c; 259 integer tmp = uint64_t(d); 260 261 //Set the blocking size, Using Pascal Giorgi's convention 262 if(blocking==0){ 263 r=tmp.bitsize()-1; 264 c=tmp.bitsize()-1; 265 }else 266 r=c=blocking; 267 268 //Create the block 269 Block U(field(),r,d); 270 Block V(field(),d,c); 271 272 //Pick random entries for U and W. W will become the last c-1 columns of V 273 274 U.random(); 275 V.random(); 276 277 278 BlackboxBlockContainer<Field, Blackbox > blockseq(&B,field(),U,V); 279 280 //Get the generator of the projection using the Coppersmith algorithm (slightly modified by Yuhasz) 281 BlockCoppersmithDomain<Domain, BlackboxBlockContainer<Field, Blackbox> > BCD(domain(), &blockseq,d); 282 std::vector<Block> gen; 283 std::vector<size_t> deg; 284 deg = BCD.right_minpoly(gen); 285 for(size_t i = 0; i < gen[0].coldim(); i++) 286 report << "Column " << i << " has degree " << deg[i] << std::endl; 287 288 //Compute the rank via the determinant of the generator 289 //Get the sum of column degrees 290 //This is the degree of the determinant via Yuhasz thesis 291 //size_t detdeg = std::accumulate(deg.begin(), deg.end(), 0); 292 size_t detdeg= 0; 293 for(size_t i = 0; i < gen[0].coldim(); i++) 294 detdeg+=deg[i]; 295 //Set up interpolation with one more evaluation point than degree 296 size_t numpoints = d+1; 297 std::vector<Element> evalpoints(numpoints), evaldets(numpoints); 298 for(typename std::vector<Element>::iterator evalit = evalpoints.begin(); evalit != evalpoints.end(); evalit++){ 299 do{ 300 //do iter.random(*evalit); while(field().isZero(*evalit)); 301 iter.random(*evalit); 302 }while ((std::find(evalpoints.begin(), evalit, *evalit) != evalit)); 303 }//end evaluation point construction loop 304 305 //Evaluate the generator determinant at the points 306 EvalPolyMat(evaldets, evalpoints, gen); 307 for(size_t k = 0; k <numpoints; k++) 308 report << evalpoints[k] << " " << evaldets[k] << std::endl; 309 //Construct the polynomial using Givare interpolation 310 //Stolen from Pascal Giorgi, linbox/examples/omp-block-rank.C 311 typedef Givaro::Poly1CRT< Field > PolyCRT; 312 PolyCRT Interpolator(field(), evalpoints, "x"); 313 typename PolyCRT::Element Determinant; 314 Interpolator.RnsToRing(Determinant,evaldets); 315 Givaro::Degree intdetdeg; 316 Interpolator.getpolydom().degree(intdetdeg,Determinant); 317 Givaro::Degree intdetval; 318 Interpolator.getpolydom().val(intdetval,Determinant); 319 if(detdeg != (size_t) intdetdeg.value()){ 320 report << "sum of column degrees " << detdeg << std::endl; 321 report << "interpolation degree " << intdetdeg.value() << std::endl; 322 } 323 report << "sum of column degrees " << detdeg << std::endl; 324 report << "interpolation degree " << intdetdeg.value() << std::endl; 325 report << "valence (trailing degree) " << intdetval.value() << std::endl; 326 for(size_t k = 0; k<gen.size(); k++) 327 domain().write(report, gen[k]) << "x^" << k << std::endl; 328 Interpolator.write(report << "Interpolated determinant: ", Determinant) << std::endl; 329 size_t myrank = size_t(intdetdeg.value() - intdetval.value()); 330 commentator().stop ("done", NULL, "Coppersmith rank"); 331 return myrank; 332 } 333 334 335 336 337 }; // end of class CoppersmithRank 338 339 //Use the coppersmith block wiedemann to compute the determinant 340 template <class _Domain> 341 class CoppersmithDeterminant{ 342 343 public: 344 typedef _Domain Domain; 345 typedef typename Domain::Field Field; 346 typedef typename Domain::Element Element; 347 typedef typename Domain::Matrix Block; 348 typedef typename Domain::Submatrix Sub; 349 typedef typename Field::RandIter Random; 350 domain()351 inline const Domain & domain() const { return *_MD; } field()352 inline const Field & field() const { return domain().field(); } 353 protected: 354 const Domain *_MD; 355 Random iter; 356 size_t blocking; 357 358 //Compute the determinant of a polynomial matrix at the given set of evaluation points 359 //Store the results in the vector dets. EvalPolyMat(std::vector<Element> & dets,std::vector<Element> & values,std::vector<Block> & mat)360 void EvalPolyMat(std::vector<Element> &dets, std::vector<Element> &values, std::vector<Block> & mat) const { 361 362 size_t deg = mat.size() -1; 363 size_t numv = values.size(); 364 //Compute the determinant of the evaluation at values[i] for each i 365 for(size_t i = 0; i<numv; i++){ 366 //copy the highest matrix coefficient 367 Block evalmat(mat[deg]); 368 //Evaluate using a horner style evaluation 369 typename std::vector<Block>::reverse_iterator addit = mat.rbegin(); 370 addit++; 371 for(addit; addit != mat.rend(); addit++){ 372 domain().mulin(evalmat,values[i]); 373 domain().addin(evalmat,*addit); 374 }//end loop computing horner evaluation 375 //Compute the determinant of the evaluation and store it in dets[i] 376 dets[i] = det(dets[i],evalmat); 377 }//end loop over evaluation points 378 }//end evaluation of polynominal matrix determinant 379 380 public: 381 CoppersmithDeterminant(const Domain &MD, size_t blocking_ = 0) : 382 _MD(&MD), blocking(blocking_), iter(MD.field()) 383 {} 384 385 386 template <class Blackbox> det(const Blackbox & B)387 Element det (const Blackbox &B) const 388 { 389 commentator().start ("Coppersmith determinant", "determinant"); 390 #if 1 391 std::ostream& report = commentator().report(Commentator::LEVEL_IMPORTANT, INTERNAL_DESCRIPTION); 392 #endif 393 394 //Set up the projection matrices and their dimensions 395 size_t d = B.coldim(); 396 size_t r,c; 397 integer tmp = uint64_t(d); 398 399 //Use given blocking size, if not given use Pascal Giorgi's convention 400 if(blocking==0){ 401 r=tmp.bitsize()-1; 402 c=tmp.bitsize()-1; 403 }else 404 r=c=blocking; 405 406 //Create the block 407 Block U(field(),r,d); 408 Block V(field(),d,c); 409 410 //Pick random entries for U and W. W will become the last c-1 columns of V 411 412 U.random(); 413 V.random(); 414 415 //Multiply V by B on the left 416 domain().leftMulin(B,V); 417 418 //Create the sequence container and its iterator that will compute the projection 419 BlackboxBlockContainer<Field, Blackbox > blockseq(&B,field(),U,V); 420 421 //Get the generator of the projection using the Coppersmith algorithm (slightly modified by Yuhasz) 422 BlockCoppersmithDomain<Domain, BlackboxBlockContainer<Field, Blackbox> > BCD(domain(), &blockseq,d); 423 std::vector<Block> gen; 424 std::vector<size_t> deg; 425 deg = BCD.right_minpoly(gen); 426 427 //Compute the determinant via the constant coefficient of the determinant of the generator 428 //Get the sum of column degrees 429 //This is the degree of the determinant via Yuhasz thesis 430 //size_t detdeg = std::accumulate(deg.begin(), deg.end(), 0); 431 size_t detdeg= 0; 432 for(size_t i = 0; i < gen[0].coldim(); i++) 433 detdeg+=deg[i]; 434 //Set up interpolation with one more evaluation point than degree 435 size_t numpoints = 2*d; 436 std::vector<Element> evalpoints(numpoints), evaldets(numpoints); 437 for(typename std::vector<Element>::iterator evalit = evalpoints.begin(); evalit != evalpoints.end(); evalit++){ 438 do{ 439 do iter.random(*evalit); while(field().isZero(*evalit)); 440 }while ((std::find(evalpoints.begin(), evalit, *evalit) != evalit)); 441 }//end evaluation point construction loop 442 443 //Evaluate the generator determinant at the points 444 EvalPolyMat(evaldets, evalpoints, gen); 445 //Construct the polynomial using Givare interpolation 446 //Stolen from Pascal Giorgi, linbox/examples/omp-block-rank.C 447 typedef Givaro::Poly1CRT<Field> PolyCRT; 448 PolyCRT Interpolator(field(), evalpoints, "x"); 449 typename PolyCRT::Element Determinant; 450 Interpolator.RnsToRing(Determinant,evaldets); 451 Givaro::Degree intdetdeg; 452 Interpolator.getpolydom().degree(intdetdeg,Determinant); 453 Givaro::Degree intdetval(0); 454 Interpolator.getpolydom().val(intdetval,Determinant); 455 if(d != (size_t)intdetdeg.value()){ 456 report << "The matrix is singular, determinant is zero" << std::endl; 457 return field(0).zero; 458 } 459 Interpolator.write(report << "Interpolated determinant: ", Determinant) << std::endl; 460 Element intdeterminant(field().zero); 461 Interpolator.getpolydom().getEntry(intdeterminant,intdetval,Determinant); 462 commentator().stop ("done", NULL, "Coppersmith determinant"); 463 return intdeterminant; 464 } 465 466 }; // end of class CoppersmithDeterminant 467 468 469 }// end of namespace LinBox 470 471 #endif //__LINBOX_coppersmith_H 472 473 // Local Variables: 474 // mode: C++ 475 // tab-width: 4 476 // indent-tabs-mode: nil 477 // c-basic-offset: 4 478 // End: 479 // vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s 480