1 /* Copyright (C) 2007 LinBox 2 * written by JG Dumas 3 * 4 * 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 /*! @file algorithms/cra-builder-single.h 26 * @ingroup algorithms 27 * @brief Chinese remaindering of a single value 28 */ 29 30 31 #ifndef __LINBOX_cra_single_H 32 #define __LINBOX_cra_single_H 33 34 #include "linbox/util/timer.h" 35 #include <stdlib.h> 36 #include "linbox/integer.h" 37 #include "linbox/solutions/methods.h" 38 #include "linbox/algorithms/cra-domain.h" 39 #include "linbox/algorithms/cra-builder-full-multip.h" 40 #include <vector> 41 #include <array> 42 #include <utility> 43 44 namespace LinBox 45 { 46 /** @brief Lower bound on number of b-bit primes 47 * @ingroup CRA 48 */ primes_count(size_t pbits)49 inline uint64_t primes_count(size_t pbits) 50 { 51 static std::array<uint64_t, 36> pctable = {{ 52 0, 0, 2, 2, 2, 5, 7, 13, 23, 43, 75, 137, 255, 53 464, 872, 1612, 3030, 5709, 10749, 20390, 38635, 54 73586, 140336, 268216, 513708, 985818, 1894120, 55 3645744, 7027290, 13561907, 26207278, 50697537, 56 98182656, 190335585, 369323305, 717267168, 57 }}; // source: http://oeis.org/A162145/list 58 static double C = 3. / (10. * log(2.)); 59 if (pbits < pctable.size()) { 60 return pctable[pbits]; 61 } 62 else if (pbits <= 71) { 63 // source: Rosser & Schoenfeld 64 return ceil((C / (pbits - 1)) * pow(2., pbits)); 65 } 66 else { 67 return std::numeric_limits<uint64_t>::max(); 68 } 69 } 70 71 /** @brief Abstract base class for CRA builders 72 * 73 * Subclasses must implement the termination functionality. 74 * 75 * @ingroup CRA 76 * 77 */ 78 template<class Domain_Type> 79 struct CRABuilderSingleBase { 80 typedef Domain_Type Domain; 81 typedef typename Domain::Element DomainElement; 82 typedef CRABuilderSingleBase<Domain> Self_t; 83 84 protected: 85 // PrimeProd*nextM_ is the modulus 86 Integer primeProd_; 87 Integer nextM_; 88 Integer residue_; // remainder to be reconstructed 89 90 #ifdef _LB_CRATIMING 91 mutable Timer tInit, tIteration, tImaging, tIRecon, tOther; 92 mutable CRATimer totalTime; 93 mutable size_t IterCounter_; 94 #endif 95 96 /** @brief Update the residue and check whether it changed 97 * 98 * The eventually-recovered number will be congruent to e modulo D. 99 * 100 * The initialize function should be called at least once before 101 * calling this one. 102 * 103 * @param D The modulus of the new image 104 * @param e The residue modulo D 105 * @returns true iff the residue changed with this update 106 */ progress_checkCRABuilderSingleBase107 bool progress_check (const Integer& D, const Integer& e) 108 { 109 // Precondition : initialize has been called once before 110 #ifdef _LB_CRATIMING 111 tIRecon.clear(); 112 tIRecon.start(); 113 ++IterCounter_; 114 #endif 115 bool was_updated = false; 116 primeProd_ *= nextM_; 117 nextM_ =D; 118 Integer u0 = residue_ % D;//0 119 Integer u1 = e % D;//e 120 Integer m0 = primeProd_;//1 121 if (u0 != u1) { 122 was_updated = true; 123 inv(m0, m0, D); // res <-- m0^{-1} mod m1//1 124 u0 = u1 - u0; // tmp <-- (u1-u0)//e 125 u0 *= m0; // res <-- (u1-u0)( m0^{-1} mod m1 )//e 126 u0 %= D; 127 Integer tmp(u0);//e 128 if (u0 < 0) 129 tmp += D;//e+D 130 else 131 tmp -= D;//e-D 132 if (absCompare(u0,tmp) > 0) u0 = tmp; 133 u0 *= primeProd_; // res <-- (u1-u0)( m0^{-1} mod m1 ) m0 and res <= (m0m1-m0) 134 residue_ += u0; // res <-- u0 + (u1-u0)( m0^{-1} mod m1 ) m0 and res < m0m1 135 } 136 #ifdef _LB_CRATIMING 137 tIRecon.stop(); 138 totalTime.ttIRecon += tIRecon; 139 #endif 140 return was_updated; 141 } 142 143 /** @brief Update the residue and check whether it changed 144 * 145 * The eventually-recovered number will be congruent to e modulo D. 146 * 147 * The initialize function should be called at least once before 148 * calling this one. 149 * 150 * @param D The modulus of the new image 151 * @param e The residue modulo D 152 * @returns true iff the residue changed with this update 153 */ progress_checkCRABuilderSingleBase154 bool progress_check (const Domain& D, const DomainElement& e) 155 { 156 // Precondition : initialize has been called once before 157 #ifdef _LB_CRATIMING 158 tIRecon.clear(); 159 tIRecon.start(); 160 ++IterCounter_; 161 #endif 162 bool was_updated = false; 163 primeProd_ *= nextM_; 164 D.characteristic( nextM_ ); 165 166 DomainElement u0; 167 if (! D.areEqual( D.init(u0, residue_), e)) { 168 was_updated = true; 169 170 D.negin(u0); // u0 <-- -u0 171 D.addin(u0, e); // u0 <-- e-u0 172 173 DomainElement m0; 174 D.init(m0, primeProd_); 175 D.invin(m0); // m0 <-- m0^{-1} mod nextM_ 176 D.mulin(u0, m0); // u0 <-- (e-u0)( m0^{-1} mod nextM_ ) 177 178 Integer res; 179 D.convert(res, u0); // res <-- (e-u0)( m0^{-1} mod nextM_ ) 180 // and res < nextM_ 181 182 Integer tmp(res); 183 tmp -= nextM_; 184 if (absCompare(res,tmp)>0) res = tmp; // Normalize 185 186 res *= primeProd_; // res <-- (e-u0)( m0^{-1} mod nextM_ ) m0 187 // and res <= (m0.nextM_-m0) 188 189 residue_ += res; // <-- u0 + (e-u0)( m0^{-1} mod nextM_ ) m0 190 // and res < m0.nextM_ 191 } 192 #ifdef _LB_CRATIMING 193 tIRecon.stop(); 194 totalTime.ttIRecon += tIRecon; 195 #endif 196 return was_updated; 197 } 198 199 public: 200 CRABuilderSingleBaseCRABuilderSingleBase201 CRABuilderSingleBase() : 202 primeProd_(1U), 203 nextM_(1U) 204 #ifdef _LB_CRATIMING 205 , IterCounter_(0) 206 #endif 207 { 208 #ifdef _LB_CRATIMING 209 clearTimers(); 210 totalTime.clear(); 211 #endif 212 } 213 214 /** @brief Initialize the CRA with the first residue. 215 * 216 * The eventually-recovered number will be congruent to e modulo D. 217 * This function must be called just once. Subsequent calls 218 * should be made to the progress() function. 219 * 220 * @param D The modulus 221 * @param e The residue 222 */ initializeCRABuilderSingleBase223 void initialize (const Integer& D, const Integer& e) 224 { 225 #ifdef _LB_CRATIMING 226 tInit.clear(); 227 tInit.start(); 228 IterCounter_=1; 229 #endif 230 primeProd_ = D; 231 nextM_ = 1U; 232 residue_ = e; 233 #ifdef _LB_CRATIMING 234 tInit.stop(); 235 totalTime.ttInit += tInit; 236 #endif 237 } 238 239 /** @brief Initialize the CRA with the first residue. 240 * 241 * The eventually-recovered number will be congruent to e modulo D. 242 * This function must be called just once. Subsequent calls 243 * should be made to the progress() function. 244 * 245 * @param D The modulus 246 * @param e The residue 247 */ initializeCRABuilderSingleBase248 void initialize (const Domain& D, const DomainElement& e) 249 { 250 #ifdef _LB_CRATIMING 251 tInit.clear(); 252 tInit.start(); 253 IterCounter_=1; 254 #endif 255 D.characteristic( primeProd_ ); 256 nextM_ = 1U; 257 D.convert( residue_, e); 258 #ifdef _LB_CRATIMING 259 tInit.stop(); 260 totalTime.ttInit += tInit; 261 #endif 262 } 263 264 /** @brief Gets the result recovered so far. 265 * 266 * (This is the same as getResidue.) 267 */ resultCRABuilderSingleBase268 Integer& result(Integer& d) 269 { 270 return d=residue_; 271 } 272 273 /** @brief Gets the result recovered so far. 274 */ getResidueCRABuilderSingleBase275 Integer& getResidue(Integer& r ) 276 { 277 return r= residue_; 278 } 279 280 /** @brief Gets the modulus of the result recovered so far. 281 */ getModulusCRABuilderSingleBase282 Integer& getModulus(Integer& m) 283 { 284 285 #ifdef _LB_CRATIMING 286 tOther.clear(); 287 tOther.start(); 288 #endif 289 m = primeProd_ * nextM_; 290 #ifdef _LB_CRATIMING 291 tOther.stop(); 292 totalTime.ttOther += tOther; 293 #endif 294 return m; 295 } 296 297 /** @brief Checks whether i is co-prime to the modulus so far. 298 * 299 * @return true iff i shares a common factor with the modulus 300 */ noncoprimeCRABuilderSingleBase301 bool noncoprime(const Integer& i) const 302 { 303 Integer g; 304 return ( (gcd(g, i, nextM_) != 1) || (gcd(g, i, primeProd_) != 1) ); 305 } 306 307 /** @brief Returns a lower bound on the number of bits in the modulus. 308 */ IntegerCRABuilderSingleBase309 decltype(Integer().bitsize()) modbits() const 310 { 311 return primeProd_.bitsize() + nextM_.bitsize() - 1; 312 } 313 ~CRABuilderSingleBaseCRABuilderSingleBase314 virtual ~CRABuilderSingleBase() {} 315 316 #ifdef _LB_CRATIMING clearTimersCRABuilderSingleBase317 void clearTimers() const 318 { 319 tInit.clear(); 320 //tIteration.clear(); 321 //tImaging.clear(); 322 tIRecon.clear(); 323 tOther.clear(); 324 } 325 public: 326 inline std::ostream& printTime(const Timer& timer, const char* title, std::ostream& os, const char* pref = "") const 327 { 328 if (timer.count() > 0) { 329 os << pref << title; 330 for (int i=strlen(title)+strlen(pref); i<28; i++) 331 os << ' '; 332 return os << timer << std::endl; 333 } 334 else 335 return os; 336 } 337 printCRATimeCRABuilderSingleBase338 inline std::ostream& printCRATime(const CRATimer& timer, const char* title, std::ostream& os) const 339 { 340 printTime(timer.ttInit, " Init: ", os, title); 341 //printTime(timer.ttImaging, "Imaging", os, title); 342 //printTime(timer.ttIteration, "Iteration", os, title); 343 printTime(timer.ttIRecon, " Integer reconstruction: ", os, title); 344 printTime(timer.ttOther, " Other: ", os, title); 345 return os; 346 } 347 reportTimesCRABuilderSingleBase348 std::ostream& reportTimes(std::ostream& os) const 349 { 350 os << "CRA Iterations:" << IterCounter_ << "\n"; 351 printCRATime(totalTime, "CRA Time", os); 352 return os; 353 } 354 #endif 355 356 }; 357 358 359 /** @brief Heuristic Chinese Remaindering with early termination 360 * 361 * This approach stops as soon as the modulus doesn't changed for some 362 * fixed number of steps in a row. 363 * 364 * @ingroup CRA 365 * 366 */ 367 template<class Domain_Type> 368 struct CRABuilderEarlySingle :public CRABuilderSingleBase<Domain_Type> { 369 typedef CRABuilderSingleBase<Domain_Type> Base; 370 typedef Domain_Type Domain; 371 typedef typename Domain::Element DomainElement; 372 typedef CRABuilderEarlySingle<Domain> Self_t; 373 374 const unsigned int EARLY_TERM_THRESHOLD; 375 376 protected: 377 unsigned int occurency_; // number of equalities 378 379 public: 380 381 /** @brief Creates a new heuristic CRA object. 382 * 383 * @param EARLY how many unchanging iterations until termination. 384 */ 385 CRABuilderEarlySingle(const size_t EARLY=LINBOX_DEFAULT_EARLY_TERMINATION_THRESHOLD) : 386 EARLY_TERM_THRESHOLD((unsigned)EARLY-1), 387 occurency_(0U) 388 { } 389 390 /** @brief Initialize the CRA with the first residue. 391 * 392 * The eventually-recovered number will be congruent to e modulo D. 393 * This function must be called just once. Subsequent calls 394 * should be made to the progress() function. 395 * 396 * Either the types of D and e should both be Integer, 397 * or D is the domain type (e.g., Modular<double>) and 398 * e is the element type (e.g., double) 399 * 400 * @param D The modulus 401 * @param e The residue 402 */ 403 template <typename ModType, typename ResType> initializeCRABuilderEarlySingle404 void initialize (const ModType& D, const ResType& e) 405 { 406 Base::initialize(D,e); 407 occurency_ = 1; 408 } 409 410 /** @brief Update the residue and termination condition. 411 * 412 * The eventually-recovered number will be congruent to e modulo D. 413 * 414 * The initialize function must be called at least once before 415 * calling this one. 416 * 417 * Either the types of D and e should both be Integer, 418 * or D is the domain type (e.g., Modular<double>) and 419 * e is the element type (e.g., double) 420 * 421 * @param D The modulus of the new image 422 * @param e The residue modulo D 423 */ 424 template <typename ModType, typename ResType> progressCRABuilderEarlySingle425 void progress (const ModType& D, const ResType& e) 426 { 427 // Precondition : initialize has been called once before 428 // linbox_check(occurency_ > 0); 429 if (Base::progress_check(D,e)) 430 occurency_ = 1; 431 else 432 occurency_ ++; 433 } 434 435 /** @brief Checks whether the CRA is (heuristically) finished. 436 * 437 * @return true iff the early termination condition has been reached. 438 */ terminatedCRABuilderEarlySingle439 bool terminated() const 440 { 441 return occurency_ > EARLY_TERM_THRESHOLD; 442 } 443 }; 444 445 446 /** @brief Chinese Remaindering with guaranteed probability bound and early termination. 447 * 448 * This strategy terminates when the probability of failure is below a 449 * certain threshold. 450 * 451 * @ingroup CRA 452 * 453 */ 454 template<class Domain_Type> 455 struct CRABuilderProbSingle :public CRABuilderSingleBase<Domain_Type> { 456 typedef CRABuilderSingleBase<Domain_Type> Base; 457 typedef Domain_Type Domain; 458 typedef typename Domain::Element DomainElement; 459 typedef CRABuilderProbSingle<Domain> Self_t; 460 461 protected: 462 const size_t bitbound_; 463 const double failbound_; 464 double curfailprob_; // the probability the result right now is incorrect 465 mod_bitsizeCRABuilderProbSingle466 size_t mod_bitsize(const Integer& D) const { return D.bitsize(); } 467 468 template <typename Field> mod_bitsizeCRABuilderProbSingle469 size_t mod_bitsize(const Field& D) const { 470 Integer p; 471 D.characteristic(p); 472 return p.bitsize(); 473 } 474 475 public: 476 /** @brief Creates a new probabilistic CRA object. 477 * 478 * @param bitbound An upper bound on the number of bits in the result. 479 * @param failprob An upper bound on the probability of failure. 480 */ 481 CRABuilderProbSingle(const size_t bitbound, const double failprob = LINBOX_DEFAULT_FAILURE_PROBABILITY) : bitbound_CRABuilderProbSingle482 bitbound_(bitbound), 483 failbound_(failprob), 484 curfailprob_(-1.) 485 { 486 } 487 488 /** @brief Initialize the CRA with the first residue. 489 * 490 * The eventually-recovered number will be congruent to e modulo D. 491 * This function must be called just once. Subsequent calls 492 * should be made to the progress() function. 493 * 494 * Either the types of D and e should both be Integer, 495 * or D is the domain type (e.g., Modular<double>) and 496 * e is the element type (e.g., double) 497 * 498 * @param D The modulus 499 * @param e The residue 500 */ 501 template <typename ModType, typename ResType> initializeCRABuilderProbSingle502 void initialize (const ModType& D, const ResType& e) 503 { 504 Base::initialize(D,e); 505 curfailprob_ = 1.; 506 } 507 508 /** @brief Update the residue and termination condition. 509 * 510 * The eventually-recovered number will be congruent to e modulo D. 511 * 512 * The initialize function must be called at least once before 513 * calling this one. 514 * 515 * Either the types of D and e should both be Integer, 516 * or D is the domain type (e.g., Modular<double>) and 517 * e is the element type (e.g., double) 518 * 519 * @param D The modulus of the new image 520 * @param e The residue modulo D 521 */ 522 template <typename ModType, typename ResType> progressCRABuilderProbSingle523 void progress (const ModType& D, const ResType& e) 524 { 525 // Precondition : initialize has been called once before 526 // linbox_check(curfailprob_ >= 0.); 527 if (Base::progress_check(D,e)) { 528 curfailprob_ = (Base::modbits() > bitbound_) ? 0. : 1.; 529 } 530 else { 531 // failure iff (failed so far) AND (this image fooled you) 532 // = curfailprob * (# bad primes in range) / (total # primes in range) 533 // <= curfailprob * (floor((bitbound_ - modbits + 1)/(pbits - 1)) / (# pbits primes) 534 auto mbits = Base::modbits(); 535 if (mbits > bitbound_) 536 curfailprob_ = 0.; 537 else { 538 auto pbits = mod_bitsize(D); 539 double failupdate = double((bitbound_ + 1 - Base::modbits()) / (pbits - 1)) 540 / primes_count(pbits); 541 if (failupdate < 1.) 542 curfailprob_ *= failupdate; 543 } 544 } 545 } 546 547 /** @brief Checks whether the CRA is (heuristically) finished. 548 * 549 * @return true iff the early termination condition has been reached. 550 */ terminatedCRABuilderProbSingle551 bool terminated() const 552 { 553 return curfailprob_ <= failbound_; 554 } 555 }; 556 557 558 /** @brief Chinese Remaindering with full precision and no chance of failure. 559 * 560 * @ingroup CRA 561 * 562 */ 563 template<class Domain_Type> 564 struct CRABuilderFullSingle :public CRABuilderFullMultip<Domain_Type>{ 565 typedef CRABuilderSingleBase<Domain_Type> Base; 566 typedef Domain_Type Domain; 567 typedef typename Domain::Element DomainElement; 568 typedef CRABuilderFullSingle<Domain> Self_t; 569 using MultiParent = CRABuilderFullMultip<Domain>; 570 571 public: 572 /** @brief Creates a new deterministic CRA object. 573 * 574 * @param bitbound An upper bound on the number of bits in the result. 575 * @param failprob An upper bound on the probability of failure. 576 */ CRABuilderFullSingleCRABuilderFullSingle577 CRABuilderFullSingle(const size_t bitbound) : 578 MultiParent((double)bitbound) 579 { 580 } 581 582 /** @brief Initialize the CRA with the first residue. 583 * 584 * The eventually-recovered number will be congruent to e modulo D. 585 * This function must be called just once. Subsequent calls 586 * should be made to the progress() function. 587 * 588 * Either the types of D and e should both be Integer, 589 * or D is the domain type (e.g., Modular<double>) and 590 * e is the element type (e.g., double) 591 * 592 * @param D The modulus 593 * @param e The residue 594 */ 595 template <typename ModType, typename ResType> initializeCRABuilderFullSingle596 void initialize (const ModType& D, const ResType& e) 597 { 598 std::array<ResType,1> v {e}; 599 MultiParent::initialize(D,v); 600 } 601 602 /** @brief Update the residue and termination condition. 603 * 604 * The eventually-recovered number will be congruent to e modulo D. 605 * 606 * The initialize function must be called at least once before 607 * calling this one. 608 * 609 * Either the types of D and e should both be Integer, 610 * or D is the domain type (e.g., Modular<double>) and 611 * e is the element type (e.g., double) 612 * 613 * @param D The modulus of the new image 614 * @param e The residue modulo D 615 */ 616 template <typename ModType, typename ResType> progressCRABuilderFullSingle617 void progress (const ModType& D, const ResType& e) 618 { 619 // Precondition : initialize has been called once before 620 std::array<ResType,1> v {e}; 621 MultiParent::progress(D, v); 622 } 623 624 /** @brief Gets the result recovered so far. 625 */ getResidueCRABuilderFullSingle626 inline const Integer& getResidue() const 627 { 628 return MultiParent::getResidue().front(); 629 } 630 631 /** @brief Gets the result recovered so far. 632 */ getResidueCRABuilderFullSingle633 inline Integer& getResidue(Integer& r) const 634 { 635 return r = getResidue(); 636 } 637 638 /** @brief alias for getResidue 639 */ resultCRABuilderFullSingle640 inline Integer& result(Integer& r) const 641 { 642 return r = getResidue(); 643 } 644 }; 645 646 } 647 648 #endif //__LINBOX_cra_single_H 649 650 // Local Variables: 651 // mode: C++ 652 // tab-width: 4 653 // indent-tabs-mode: nil 654 // c-basic-offset: 4 655 // End: 656 // vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s 657