| 132 | | /* |
| 133 | | template <class Ring> |
| 134 | | class WangRationalReconstruction: public RationalReconstructionBase<Ring> |
| 135 | | { |
| 136 | | protected: |
| 137 | | bool _reduce; |
| 138 | | bool _recursive; |
| 139 | | Ring _Z; |
| 140 | | OpCounter C; |
| 141 | | public: |
| 142 | | Ring _Z; |
| 143 | | typedef typename Ring::Element Element; |
| 144 | | |
| 145 | | OpCounter C; |
| 146 | | |
| 147 | | WangRationalReconstruction(const Ring& Z, const bool reduce = true, const bool recursive = false): RationalReconstructionBase<Ring>(Z) { |
| 148 | | _reduce = reduce; _recursive = recursive; |
| 149 | | } |
| 150 | | |
| 151 | | WangRationalReconstruction<Ring> (const WangClassicRationalReconstruction<Ring>& RR): RationalReconstructionBase<Ring>(RR._Z) { |
| 152 | | _reduce = RR._reduce; _recursive = RR._recursive; |
| 153 | | } |
| 154 | | |
| 155 | | bool reconstructRational(Element& a, Element& b, const Element& x, const Element& m)=0 { |
| 156 | | Element a_bound; _Z.sqrt(a_bound, m/2); |
| 157 | | bool res = rationalReconstruction(a,b,x,m,a_bound); |
| 158 | | res = res && (b <= a_bound); |
| 159 | | return res; |
| 160 | | } |
| 161 | | |
| 162 | | bool reconstructRational(Element& a, Element& b, const Element& x, const Element& m, const Element& a_bound) { |
| 163 | | bool res; |
| 164 | | |
| 165 | | if (x == 0) { |
| 166 | | a = 0; |
| 167 | | b = 1; |
| 168 | | } else { |
| 169 | | bool res = ratrecon(a,b,x,m,a_bound); |
| 170 | | if (_recursive) |
| 171 | | for(Element newbound = a_bound + 1; (!res) && (newbound<x) ; ++newbound) |
| 172 | | res = ratrecon(a,b,x,m,newbound); |
| 173 | | } |
| 174 | | if (!res) { |
| 175 | | a = x> m/2? x-m: x; |
| 176 | | b = 1; |
| 177 | | if (a > 0) res = (a < a_bound); |
| 178 | | else res = (-a < a_bound); |
| 179 | | } |
| 180 | | |
| 181 | | return res; |
| 182 | | } |
| 183 | | |
| 184 | | virtual ~WangRationalReconstruction() {} |
| 185 | | protected: |
| 186 | | |
| 187 | | // m> x> 0 |
| 188 | | bool ratrecon(Element& a,Element& b,const Element& x,const Element& m,const Element& a_bound) { |
| 189 | | |
| 190 | | Element r0, t0, q, u; |
| 191 | | r0=m; |
| 192 | | t0=0; |
| 193 | | a=x; |
| 194 | | b=1; |
| 195 | | //Element s0,s; s0=1,s=0;//test time gcdex; |
| 196 | | while(a>=a_bound) |
| 197 | | { |
| 198 | | |
| 199 | | q = r0; |
| 200 | | _Z.divin(q,a); // r0/num |
| 201 | | ++C.div_counter; |
| 202 | | |
| 203 | | u = a; |
| 204 | | a = r0; |
| 205 | | r0 = u; // r0 <-- num |
| 206 | | |
| 207 | | _Z.axmyin(a,u,q); // num <-- r0-q*num |
| 208 | | ++C.mul_counter; |
| 209 | | //if (a == 0) return false; |
| 210 | | |
| 211 | | u = b; |
| 212 | | b = t0; |
| 213 | | t0 = u; // t0 <-- den |
| 214 | | |
| 215 | | _Z.axmyin(b,u,q); // den <-- t0-q*den |
| 216 | | ++C.mul_counter; |
| 217 | | } |
| 218 | | |
| 219 | | // if (den < 0) { |
| 220 | | // _Z.negin(num); |
| 221 | | // _Z.negin(den); |
| 222 | | // } |
| 223 | | |
| 224 | | if ((a>0) && (_reduce)) { |
| 225 | | |
| 226 | | // [GG, MCA, 1999] Theorem 5.26 |
| 227 | | // (ii) |
| 228 | | Element gg; |
| 229 | | ++C.gcd_counter; |
| 230 | | if (_Z.gcd(gg,a,b) != 1) { |
| 231 | | |
| 232 | | Element ganum, gar2; |
| 233 | | for( q = 1, ganum = r0-a, gar2 = r0 ; (ganum >= a_bound) || (gar2<a_bound); ++q ) { |
| 234 | | ganum -= a; |
| 235 | | gar2 -= a; |
| 236 | | } |
| 237 | | |
| 238 | | //_Z.axmyin(r0,q,a); |
| 239 | | r0 = ganum; |
| 240 | | _Z.axmyin(t0,q,b); |
| 241 | | ++C.mul_counter;++C.mul_counter; |
| 242 | | if (t0 < 0) { |
| 243 | | a = -r0; |
| 244 | | b = -t0; |
| 245 | | } else { |
| 246 | | a = r0; |
| 247 | | b = t0; |
| 248 | | } |
| 249 | | |
| 250 | | // if (t0 > m/k) { |
| 251 | | if ((double)b > (double)m/(double)a_bound) { |
| 252 | | if (!_recursive) { |
| 253 | | std::cerr |
| 254 | | << "*** Error *** No rational reconstruction of " |
| 255 | | << x |
| 256 | | << " modulo " |
| 257 | | << m |
| 258 | | << " with denominator <= " |
| 259 | | << (m/a_bound) |
| 260 | | << std::endl; |
| 261 | | } |
| 262 | | return false; |
| 263 | | } |
| 264 | | if (_Z.gcd(gg,a,b) != 1) { |
| 265 | | if (!_recursive) |
| 266 | | std::cerr |
| 267 | | << "*** Error *** There exists no rational reconstruction of " |
| 268 | | << x |
| 269 | | << " modulo " |
| 270 | | << m |
| 271 | | << " with |numerator| < " |
| 272 | | << a_bound |
| 273 | | << std::endl |
| 274 | | << "*** Error *** But " |
| 275 | | << a |
| 276 | | << " = " |
| 277 | | << b |
| 278 | | << " * " |
| 279 | | << x |
| 280 | | << " modulo " |
| 281 | | << m |
| 282 | | << std::endl; |
| 283 | | return false; |
| 284 | | } |
| 285 | | } |
| 286 | | } |
| 287 | | // (i) |
| 288 | | if (b < 0) { |
| 289 | | _Z.negin(a); |
| 290 | | _Z.negin(b); |
| 291 | | } |
| 292 | | |
| 293 | | // std::cerr << "RatRecon End " << num << "/" << den << std::endl; |
| 294 | | return true; |
| 295 | | } |
| 296 | | }; |
| 297 | | */ |
| 298 | | template <class Ring> |
| 299 | | class QMatrix { |
| 300 | | public: |
| 301 | | Ring _Z; |
| 302 | | typedef typename Ring::Element Element; |
| 303 | | |
| 304 | | Element a,b,c,d; |
| 305 | | Element q; |
| 306 | | |
| 307 | | QMatrix(const Ring Z): _Z(Z) { |
| 308 | | a=1;b=0;c=0;d=1;q=0; |
| 309 | | } |
| 310 | | |
| 311 | | QMatrix(const QMatrix& Q): _Z(Q._Z) { |
| 312 | | a = Q.a; |
| 313 | | b = Q.b; |
| 314 | | c = Q.c; |
| 315 | | d = Q.d; |
| 316 | | q = Q.q; |
| 317 | | } |
| 318 | | |
| 319 | | QMatrix(Ring Z,const Element& ai, const Element& bi, const Element& ci, const Element& di, const Element& qi=0): _Z(Z) { |
| 320 | | a = ai; |
| 321 | | b = bi; |
| 322 | | c = ci; |
| 323 | | d = di; |
| 324 | | q = qi; |
| 325 | | } |
| 326 | | |
| 327 | | QMatrix& operator=(const QMatrix& Q) { |
| 328 | | a = Q.a; |
| 329 | | b = Q.b; |
| 330 | | c = Q.c; |
| 331 | | d = Q.d; |
| 332 | | q = Q.q; |
| 333 | | return *this; |
| 334 | | } |
| 335 | | |
| 336 | | void leftmultiply(const QMatrix Q) { |
| 337 | | leftmultiply(Q.a,Q.b, Q.c,Q.d); |
| 338 | | } |
| 339 | | |
| 340 | | void leftmultiply(const Element& ai,const Element& bi,const Element& ci,const Element& di) { |
| 341 | | Element tmpa,tmpb,tmpc,tmpd; |
| 342 | | _Z.mul(tmpa,ai,a); |
| 343 | | _Z.axpyin(tmpa,bi,c); |
| 344 | | |
| 345 | | _Z.mul(tmpb,ai,b); |
| 346 | | _Z.axpyin(tmpb,bi,d); |
| 347 | | |
| 348 | | _Z.mul(tmpc,ci,a); |
| 349 | | _Z.axpyin(tmpc,di,c); |
| 350 | | |
| 351 | | _Z.mul(tmpd,ci,b); |
| 352 | | _Z.axpyin(tmpd,di,d); |
| 353 | | |
| 354 | | a = tmpa; b = tmpb; c = tmpc; d = tmpd; |
| 355 | | } |
| 356 | | |
| 357 | | QMatrix& max(const QMatrix& max1, const QMatrix max2) { |
| 358 | | if (max1.q >= max2.q) return QMatrix(max1); |
| 359 | | else return QMatrix(max2); |
| 360 | | } |
| 361 | | |
| 362 | | bool maxin(const QMatrix max2) { |
| 363 | | if (q >= max2.q) return true; |
| 364 | | else { |
| 365 | | a = max2.a; |
| 366 | | b = max2.b; |
| 367 | | c = max2.c; |
| 368 | | d = max2.d; |
| 369 | | q = max2.q; |
| 370 | | return false; |
| 371 | | } |
| 372 | | } |
| 373 | | }; |
| 374 | | |
| 375 | | template <class Ring> |
| 376 | | class myQueue: public deque<QMatrix<Ring > > { |
| 377 | | public: |
| 378 | | typedef typename Ring::Element Element; |
| 379 | | typedef QMatrix<Ring> QMatrix; |
| 380 | | |
| 381 | | size_t _maxSize; |
| 382 | | size_t _size; |
| 383 | | |
| 384 | | myQueue(const size_t K=0) { |
| 385 | | _maxSize = K; |
| 386 | | _size = 0; |
| 387 | | } |
| 388 | | |
| 389 | | bool pushpop(QMatrix& top, const QMatrix& bottom) { |
| 390 | | if (_size+1 < _maxSize) { |
| 391 | | push_back(bottom); |
| 392 | | return false; |
| 393 | | ++_size; |
| 394 | | } else { |
| 395 | | if (!this->empty()) { |
| 396 | | top = this->front(); |
| 397 | | this->pop_front(); |
| 398 | | push_back(bottom); |
| 399 | | |
| 400 | | } else { |
| 401 | | top=bottom; |
| 402 | | } |
| 403 | | return true; |
| 404 | | } |
| 405 | | } |
| 406 | | |
| 407 | | QMatrix& clearmax(QMatrix& max1) { |
| 408 | | while (!this->empty()) { |
| 409 | | QMatrix max2(this->front()); |
| 410 | | if (max2.q > max1.q) return max1=max2; |
| 411 | | this->pop_front(); |
| 412 | | } |
| 413 | | } |
| 414 | | }; |
| 415 | | |
| 416 | | template <class Ring> |
| 417 | | class MaxQFastRationalReconstruction: public RationalReconstructionBase<Ring> |
| 418 | | { |
| 419 | | protected: |
| 420 | | size_t _threshold; |
| 421 | | Ring _Z; |
| 422 | | OpCounter C; |
| 423 | | |
| 424 | | typedef typename Ring::Element Element; |
| 425 | | public: |
| 426 | | |
| 427 | | MaxQFastRationalReconstruction(const Ring& Z): RationalReconstructionBase<Ring>(Z), _Z(Z) { |
| 428 | | _threshold = __FASTRR_DEFAULT_THRESHOLD; |
| 429 | | if (_threshold <5) _threshold=5; |
| 430 | | } |
| 431 | | |
| 432 | | ~MaxQFastRationalReconstruction() {} |
| 433 | | |
| 434 | | bool reconstructRational(Element& a, Element& b, const Element& x, const Element& m) { |
| 435 | | //reconstructRational(a,b,x,m,1); |
| 436 | | return fastQMaxReconstructRational(a,b,x,m); |
| 437 | | } |
| 438 | | |
| 439 | | bool reconstructRational(Element& a, Element& b, const Element& x, const Element& m, const Element& a_bound) { |
| 440 | | bool res = reconstructRational(a,b,x,m); |
| 441 | | return res && (a < a_bound); |
| 442 | | } |
| 443 | | |
| 444 | | protected: |
| 445 | | |
| 446 | | Element cur_ri; |
| 447 | | Element cur_rinext; |
| 448 | | Element cur_ainext; |
| 449 | | Element cur_qinext; |
| 450 | | |
| 451 | | //QMatrix max; |
| 452 | | |
| 453 | | Element& powtwo(Element& h, const Element& log_h) { |
| 454 | | h = 1; |
| 455 | | //cout << "max" << ULONG_MAX << "x" << x << "?" << (x < ULONG_MAX); |
| 456 | | if (log_h < ULONG_MAX) { |
| 457 | | h<<=log_h; |
| 458 | | //cout << "z"<< z; |
| 459 | | return h; |
| 460 | | } else { |
| 461 | | Element n,m; |
| 462 | | quoRem(n,m,log_h,(Element)(ULONG_MAX-1)); |
| 463 | | for (int i=0; i < n; ++i) { |
| 464 | | h <<=(long int)(ULONG_MAX-1); |
| 465 | | } |
| 466 | | h <= (long int)m; |
| 467 | | return h; |
| 468 | | } |
| 469 | | } |
| 470 | | |
| 471 | | Element& powtwo(Element& h, const size_t log_h) { |
| 472 | | h = 1; |
| 473 | | h<<=log_h; |
| 474 | | return h; |
| 475 | | } |
| 476 | | |
| 477 | | bool fastQMaxReconstructRational(Element& n, Element& d, const Element& x, const Element& m) { |
| 478 | | |
| 479 | | size_t log_m = m.bitsize()-1; //true unless m = 2^k |
| 480 | | |
| 481 | | Element ai, bi, ci, di; |
| 482 | | ai=1;bi=0;ci=0;di=1; |
| 483 | | |
| 484 | | cur_ri = m; |
| 485 | | cur_rinext = x; |
| 486 | | _Z.quo(cur_ainext, cur_ri, cur_rinext); |
| 487 | | cur_qinext = cur_ainext;// ri/ri_next |
| 488 | | |
| 489 | | Element powh; |
| 490 | | myQueue<Ring > queueMax(0); |
| 491 | | QMatrix<Ring > maxQ(_Z); |
| 492 | | |
| 493 | | if (!fastQMaxEEA (ai,bi,ci,di,m,log_m,x,powtwo(powh,log_m+1), log_m+1,queueMax,maxQ)) { |
| 494 | | //cout << "false\n" << flush; |
| 495 | | return false; |
| 496 | | } |
| 497 | | |
| 498 | | if (cur_rinext != 0) { |
| 499 | | cout << "bad bounds - should not happen\n" << flush; |
| 500 | | } |
| 501 | | |
| 502 | | if (!queueMax.empty()) { |
| 503 | | cout << "Queue is empty\n - sth wrong" << flush; |
| 504 | | } |
| 505 | | |
| 506 | | _Z.mul(n,x, maxQ.a); |
| 507 | | _Z.axmyin(n,m,maxQ.c); |
| 508 | | //n = x_in*ai-m*ci; |
| 509 | | d = maxQ.a; |
| 510 | | |
| 511 | | /* |
| 512 | | _Z.mul(n, maxQ.d, m); |
| 513 | | _Z.axmyin(n, maxQ.b, n); |
| 514 | | d = maxQ.b; |
| 515 | | */ |
| 516 | | Element _gcd; |
| 517 | | if (_Z.gcd(_gcd, n,d) > 1) { |
| 518 | | /* solution not implemented */ |
| 519 | | cout << "Solution pair is not coprime\n"; |
| 520 | | return false; |
| 521 | | } |
| 522 | | |
| 523 | | return true; |
| 524 | | } |
| 525 | | |
| 526 | | bool classicQMaxEEA(Element& ai, Element& bi, Element& ci, Element& di, const Element& r0, const Element& r1,const Element& powh, myQueue<Ring >& queueMax, QMatrix<Ring>& maxQ) { |
| 527 | | Element ri, rinext; |
| 528 | | ri = r0; rinext = r1; |
| 529 | | ai =di = 1; |
| 530 | | bi =ci = 0; |
| 531 | | if (rinext==0) return true; |
| 532 | | Element ainext,cinext; |
| 533 | | Element qinext; |
| 534 | | _Z.quo(qinext,ri,rinext); |
| 535 | | ++C.div_counter; |
| 536 | | |
| 537 | | ainext =qinext; |
| 538 | | //_Z.axpy(ainext, ai,qinext,bi); |
| 539 | | while (ainext <= powh) { |
| 540 | | |
| 541 | | //_Z.quo(qinext,ri,rinext); |
| 542 | | ++C.div_counter; |
| 543 | | QMatrix<Ring> newQ(_Z,ai,bi,ci,di,qinext); |
| 544 | | QMatrix<Ring> top(_Z); |
| 545 | | if (queueMax.pushpop(top, newQ)) { |
| 546 | | //cout << "1new max " << top.q << "\n" << flush; |
| 547 | | if (maxQ.q < top.q) maxQ = top; |
| 548 | | } |
| 549 | | //cout << "EEA" << qinext << ","; |
| 550 | | //++i; |
| 551 | | //ainext = ai*qinext + bi; |
| 552 | | //_Z.axpy(ainext, ai,qinext,bi); |
| 553 | | ++C.mul_counter; |
| 554 | | |
| 555 | | Element tmpbi = bi; |
| 556 | | Element tmpdi = di; |
| 557 | | bi = ai; |
| 558 | | ai = ainext; |
| 559 | | Element temp = ci; |
| 560 | | _Z.axpy(ci, temp,qinext,di); |
| 561 | | ++C.mul_counter; |
| 562 | | di = temp; |
| 563 | | |
| 564 | | temp = ri; |
| 565 | | ri=rinext; |
| 566 | | _Z.axpy(rinext, -qinext,ri,temp);//rinext = ri-qinext*rinext |
| 567 | | |
| 568 | | ++C.mul_counter; |
| 569 | | if (rinext==0) { |
| 570 | | //ainext = infinity |
| 571 | | cur_ri = ri; |
| 572 | | cur_rinext = rinext; |
| 573 | | cur_ainext = r0+1; |
| 574 | | cur_qinext = cur_ainext;//infinity |
| 575 | | //cout << "->1E:" << cur_ri << " " << cur_rinext <<" " << cur_ainext << "\n"<< flush; |
| 576 | | return true; |
| 577 | | } |
| 578 | | |
| 579 | | _Z.quo(qinext,ri,rinext); |
| 580 | | ++C.div_counter; |
| 581 | | _Z.axpy(ainext, ai,qinext,bi); |
| 582 | | |
| 583 | | } |
| 584 | | cur_ri = ri; |
| 585 | | cur_rinext = rinext; |
| 586 | | cur_ainext = ainext; |
| 587 | | cur_qinext = qinext; |
| 588 | | //cout << "->2E:" << cur_ri << " " << cur_rinext <<" " << cur_ainext << "\n"<< flush; |
| 589 | | return true; |
| 590 | | } |
| 591 | | |
| 592 | | bool fastQMaxEEA(Element& ai,Element& bi, Element& ci, Element& di, const Element& m, const size_t d, const Element& n, const Element& powh, const size_t& h, myQueue<Ring >& queueMax, QMatrix<Ring>& maxQ) { |
| 593 | | |
| 594 | | //cout << m << " " << n << " " << powh << "\n" << flush; |
| 595 | | |
| 596 | | ai=Element(1); |
| 597 | | di=Element(1); |
| 598 | | bi=Element(0); |
| 599 | | ci=Element(0); //Q(0)=Id |
| 600 | | |
| 601 | | if (m==n) { |
| 602 | | cur_ri = n; |
| 603 | | cur_rinext = 0; |
| 604 | | cur_ainext = m+1;//infinity |
| 605 | | cur_qinext = m+1; |
| 606 | | ai=1; bi=1; ci=1; di=0; |
| 607 | | QMatrix<Ring> newQ(_Z,1,0,0,1,1); |
| 608 | | //QMatrix<Ring> newQ(_Z,1,1,1,0,1); |
| 609 | | QMatrix<Ring> top(_Z); |
| 610 | | if (queueMax.pushpop(top, newQ)) { |
| 611 | | //cout << "2new max " << top.q << "\n"<< flush; |
| 612 | | if (maxQ.q < top.q) maxQ = top; |
| 613 | | } |
| 614 | | //cout << "m=n1" << ","; |
| 615 | | return true; |
| 616 | | } |
| 617 | | |
| 618 | | if (powh < 1) return false; //should not happen |
| 619 | | if (n < 2) { |
| 620 | | //cout << "n=1\n"<< flush; |
| 621 | | //n==1 -> Q1=(m,1//1,0) |
| 622 | | //n==0 -> Q1=infinity |
| 623 | | cur_ri = m; |
| 624 | | cur_rinext = 1; |
| 625 | | cur_ainext = m;//infinity |
| 626 | | cur_qinext = m; |
| 627 | | if (cur_ainext > powh) { |
| 628 | | //we do not have to treat identity; |
| 629 | | return true; |
| 630 | | } |
| 631 | | else { |
| 632 | | ai=m; bi=1; ci=1; di=0; |
| 633 | | cur_ri = 1; |
| 634 | | cur_rinext = 0; |
| 635 | | cur_ainext = m+1; |
| 636 | | cur_qinext = m+1; |
| 637 | | QMatrix<Ring> newQ(_Z,1,0,0,1,m); |
| 638 | | //QMatrix<Ring> newQ(_Z,m,1,1,0,m); |
| 639 | | QMatrix<Ring> top(_Z); |
| 640 | | if (queueMax.pushpop(top, newQ)) { |
| 641 | | // cout << "3new max " << top.q << "\n"<< flush; |
| 642 | | if (maxQ.q < top.q) maxQ = top; |
| 643 | | } |
| 644 | | //cout << "n=1" << m << ","; |
| 645 | | |
| 646 | | return true; |
| 647 | | } |
| 648 | | return true; |
| 649 | | } |
| 650 | | |
| 651 | | if (h < 1) { |
| 652 | | //cout << "h < 1" << flush; |
| 653 | | //quo(qinext,m,n); |
| 654 | | //if (qinext==1) { return (1,1,1,0) or (1,0,0,1) |
| 655 | | if ((n << 1) > m) { |
| 656 | | cur_ri = n; |
| 657 | | cur_rinext = m-n; |
| 658 | | _Z.quo(cur_qinext, cur_ri, cur_rinext); |
| 659 | | ++C.div_counter; |
| 660 | | cur_ainext = cur_qinext + 1; |
| 661 | | ai=1; bi=1; ci=1; di=0; |
| 662 | | QMatrix<Ring> newQ(_Z,1,0,0,1,1); |
| 663 | | //QMatrix<Ring> newQ(_Z,1,1,1,0,1); |
| 664 | | QMatrix<Ring> top(_Z); |
| 665 | | if (queueMax.pushpop(top, newQ)) { |
| 666 | | //cout << "4new max " << top.q << "\n"<< flush; |
| 667 | | if (maxQ.q < top.q) maxQ = top; |
| 668 | | } |
| 669 | | //cout << "h=01" << ","; |
| 670 | | |
| 671 | | } else { |
| 672 | | //we do not have to treat identity; |
| 673 | | cur_ri = m; |
| 674 | | cur_rinext = n; |
| 675 | | _Z.quo(cur_ainext,m,n); |
| 676 | | ++C.div_counter; |
| 677 | | cur_qinext = cur_ainext; |
| 678 | | } |
| 679 | | return true; |
| 680 | | } |
| 681 | | |
| 682 | | if (n < _threshold) { //what about m? |
| 683 | | return classicQMaxEEA(ai,bi,ci,di,m,n,powh,queueMax,maxQ); |
| 684 | | } |
| 685 | | |
| 686 | | //size_t log_n = n.bitsize()-1; |
| 687 | | |
| 688 | | // cout << d << " " << log_n << " " << h << "\n"<< flush; |
| 689 | | |
| 690 | | if (2*h+1 < d) { |
| 691 | | //cout << "choice1"<< flush; |
| 692 | | //if (2*h+1 <= log_n) { |
| 693 | | if (true) { |
| 694 | | //cout << ".1\n"<< flush; |
| 695 | | size_t lambda = d-2*h-1; |
| 696 | | |
| 697 | | Element aistar, bistar, cistar, distar; |
| 698 | | aistar=distar=1; |
| 699 | | bistar=cistar=0; |
| 700 | | |
| 701 | | Element mstar = m >> (long unsigned int) lambda; |
| 702 | | Element nstar = n >> (long unsigned int) lambda; |
| 703 | | |
| 704 | | size_t log_mstar = 2*h+1; |
| 705 | | |
| 706 | | queueMax._maxSize +=2; |
| 707 | | if (nstar > 0) if (!fastQMaxEEA(aistar, bistar, cistar, distar, mstar, log_mstar,nstar, powh, h, queueMax, maxQ)) return false; |
| 708 | | |
| 709 | | if (queueMax._size > 1) { |
| 710 | | queueMax.pop_back(); |
| 711 | | QMatrix<Ring> Q_i_2 (queueMax.back()); |
| 712 | | queueMax.pop_back();//pop Q*(i-1) q*(i) |
| 713 | | --queueMax._size ; |
| 714 | | --queueMax._size ; |
| 715 | | ai = Q_i_2.a; |
| 716 | | bi = Q_i_2.b; |
| 717 | | ci = Q_i_2.c; |
| 718 | | di = Q_i_2.d;// Q(i-2)= Q*(i-2) q*(i-1) |
| 719 | | } else { |
| 720 | | queueMax.clear(); |
| 721 | | } |
| 722 | | queueMax._maxSize -=2; |
| 723 | | |
| 724 | | _Z.mul(cur_ri,m,di); |
| 725 | | _Z.axmyin(cur_ri,n,bi); |
| 726 | | _Z.mul(cur_rinext,n,ai); |
| 727 | | _Z.axmyin(cur_rinext,m,ci); |
| 728 | | C.mul_counter+=4; |
| 729 | | //ri = m*di - n*bi; |
| 730 | | //rinext = -m*ci + n*ai; |
| 731 | | |
| 732 | | if (cur_ri < 0) { |
| 733 | | cur_ri = -cur_ri; |
| 734 | | |
