1 // Copyright (c) 2008 Max-Planck-Institute Saarbruecken (Germany). 2 // All rights reserved. 3 // 4 // This file is part of CGAL (www.cgal.org) 5 // 6 // $URL: https://github.com/CGAL/cgal/blob/v5.3/Polynomial/include/CGAL/Polynomial/subresultants.h $ 7 // $Id: subresultants.h 0779373 2020-03-26T13:31:46+01:00 Sébastien Loriot 8 // SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial 9 // 10 // Author(s) : Michael Kerber <mkerber@mpi-inf.mpg.de> 11 // 12 // ============================================================================ 13 #ifndef CGAL_POLYNOMIAL_SUBRESULTANTS_H 14 #define CGAL_POLYNOMIAL_SUBRESULTANTS_H 15 16 #include <list> 17 18 #include <CGAL/Polynomial_traits_d.h> 19 #include <CGAL/Polynomial/bezout_matrix.h> 20 21 namespace CGAL { 22 23 24 namespace internal { 25 26 // Intern function needed for Ducos algorithm 27 lazard_optimization(typename Polynomial_traits_d::Coefficient_type y,double n,typename Polynomial_traits_d::Polynomial_d B,typename Polynomial_traits_d::Polynomial_d & C)28 template<typename Polynomial_traits_d> void lazard_optimization 29 (typename Polynomial_traits_d::Coefficient_type y, 30 double n, 31 typename Polynomial_traits_d::Polynomial_d B, 32 typename Polynomial_traits_d::Polynomial_d& C) { 33 34 typedef typename Polynomial_traits_d::Coefficient_type NT; 35 typename CGAL::Algebraic_structure_traits<NT>::Integral_division idiv; 36 37 CGAL_precondition(n>0); 38 NT x = typename Polynomial_traits_d::Leading_coefficient() (B); 39 double a = pow(2.,std::floor(log(n)/log(2.))); 40 NT c = x; 41 n -= a; 42 while(a!=1) { 43 a/=2; 44 c=idiv(c*c,y); 45 if(n>=a) { 46 c=idiv(c*x,y); 47 n-=a; 48 } 49 } 50 C=c*B/y; 51 } 52 53 template<typename Polynomial_traits_d> lickteig_roy_optimization(typename Polynomial_traits_d::Polynomial_d A,typename Polynomial_traits_d::Polynomial_d B,typename Polynomial_traits_d::Polynomial_d C,typename Polynomial_traits_d::Coefficient_type s,typename Polynomial_traits_d::Polynomial_d & D)54 void lickteig_roy_optimization 55 (typename Polynomial_traits_d::Polynomial_d A, 56 typename Polynomial_traits_d::Polynomial_d B, 57 typename Polynomial_traits_d::Polynomial_d C, 58 typename Polynomial_traits_d::Coefficient_type s, 59 typename Polynomial_traits_d::Polynomial_d& D) { 60 61 typedef typename Polynomial_traits_d::Polynomial_d Poly; 62 typedef typename Polynomial_traits_d::Coefficient_type NT; 63 typename Polynomial_traits_d::Degree degree; 64 typename Polynomial_traits_d::Leading_coefficient lcoeff; 65 typename Polynomial_traits_d::Construct_polynomial construct; 66 typename Polynomial_traits_d::Get_coefficient coeff; 67 68 int d = degree(A), e = degree(B); 69 CGAL_precondition(d>=e); 70 std::vector<Poly> H(d+1); 71 std::list<NT> initial; 72 initial.push_front(lcoeff(C)); 73 for(int i=0;i<e;i++) { 74 H[i] = construct(initial.begin(),initial.end()); 75 initial.push_front(NT(0)); 76 } 77 H[e]=construct(initial.begin(),initial.end())-C; 78 CGAL_assertion(degree(H[e])<e); 79 initial.clear(); 80 std::copy(H[e].begin(),H[e].end(),std::back_inserter(initial)); 81 initial.push_front(NT(0)); 82 for(int i=e+1;i<d;i++) { 83 H[i]=construct(initial.begin(),initial.end()); 84 NT h_i_e=H[i].degree()>=e ? coeff(H[i],e) : NT(0); 85 H[i]-=(h_i_e*B)/lcoeff(B); 86 initial.clear(); 87 std::copy(H[i].begin(),H[i].end(),std::back_inserter(initial)); 88 initial.push_front(NT(0)); 89 } 90 H[d]=construct(initial.begin(),initial.end()); 91 D=construct(0); 92 for(int i=0;i<d;i++) { 93 D+=A[i]*H[i]; 94 } 95 D/=lcoeff(A); 96 NT Hde = degree(H[d])>=e ? coeff(H[d],e) : NT(0); 97 D=(lcoeff(B)*(H[d]+D)-Hde*B)/s; 98 if((d-e)%2==0) { 99 D=-D; 100 } 101 return; 102 } 103 104 template<typename Polynomial_traits_d> 105 typename Polynomial_traits_d::Coefficient_type resultant_for_constant_polynomial(typename Polynomial_traits_d::Polynomial_d P,typename Polynomial_traits_d::Polynomial_d Q)106 resultant_for_constant_polynomial 107 (typename Polynomial_traits_d::Polynomial_d P, 108 typename Polynomial_traits_d::Polynomial_d Q) { 109 110 typedef typename Polynomial_traits_d::Polynomial_d Polynomial; 111 typedef typename Polynomial_traits_d::Coefficient_type NT; 112 typename Polynomial_traits_d::Leading_coefficient lcoeff; 113 typename Polynomial_traits_d::Degree degree; 114 typename CGAL::Algebraic_structure_traits<Polynomial>::Is_zero is_zero; 115 CGAL_assertion(degree(P) < 1 || degree(Q) < 1); 116 if(is_zero(P) || is_zero(Q) ) { 117 return NT(0); 118 } 119 if(degree(P)==0) { 120 return CGAL::ipower(lcoeff(P),degree(Q)); 121 } else { 122 return CGAL::ipower(lcoeff(Q),degree(P)); 123 } 124 } 125 126 127 /*! 128 * \brief Compute the sequence of subresultants with pseudo-division 129 * 130 * This is an implementation of Ducos' algorithm. It improves on the 131 * classical methods for subresultant computation by reducing the 132 * swell-up of intermediate results. For all details, see 133 * L.Ducos: Optimazations of the Subresultant algorithm. <i>Journal of Pure 134 * and Applied Algebra</i> <b>145</b> (2000) 149--163 135 */ 136 template <typename Polynomial_traits_d,typename OutputIterator> inline prs_polynomial_subresultants(typename Polynomial_traits_d::Polynomial_d P,typename Polynomial_traits_d::Polynomial_d Q,OutputIterator out)137 OutputIterator prs_polynomial_subresultants 138 (typename Polynomial_traits_d::Polynomial_d P, 139 typename Polynomial_traits_d::Polynomial_d Q, 140 OutputIterator out) { 141 142 typedef typename Polynomial_traits_d::Polynomial_d Polynomial; 143 typedef typename Polynomial_traits_d::Coefficient_type NT; 144 typename Polynomial_traits_d::Leading_coefficient lcoeff; 145 typename Polynomial_traits_d::Degree degree; 146 typename Polynomial_traits_d::Construct_polynomial construct; 147 typename CGAL::Algebraic_structure_traits<Polynomial>::Is_zero is_zero; 148 149 if(degree(P) < 1 || degree(Q) < 1) { 150 *out++ = Polynomial(CGAL::internal::resultant_for_constant_polynomial 151 <Polynomial_traits_d> (P,Q)); 152 return out; 153 } 154 155 bool poly_swapped = (degree(P) < degree(Q)); 156 157 if(poly_swapped) { 158 std::swap(P,Q); 159 } 160 161 Polynomial zero_pol = construct(NT(0)); 162 std::vector<Polynomial> sres; 163 164 int deg_diff=degree(P)-degree(Q); 165 166 if(deg_diff==0) { 167 sres.push_back(Q); 168 } else { 169 sres.push_back(CGAL::ipower(lcoeff(Q),deg_diff-1)*Q); 170 } 171 172 Polynomial A,B,C,D,dummy_pol; 173 NT s,dummy_nt; 174 int delta,d,e; 175 176 A=Q; 177 178 s=CGAL::ipower(lcoeff(Q),deg_diff); 179 180 typename Polynomial_traits_d::Pseudo_division() 181 (P, -Q, dummy_pol, B, dummy_nt); 182 183 while(true) { 184 d=degree(A); 185 e=degree(B); 186 if(is_zero(B)) { 187 for(int i=0;i<d;i++) { 188 sres.push_back(zero_pol); 189 } 190 break; 191 } 192 sres.push_back(B); 193 delta=d-e; 194 if(delta>1) { 195 CGAL::internal::lazard_optimization<Polynomial_traits_d> 196 (s,double(delta-1),B,C); 197 //C=CGAL::ipower(CGAL::integral_division(lcoeff(B),s),delta-1)*B; 198 for(int i=0;i<delta-2;i++) { 199 sres.push_back(zero_pol); 200 } 201 sres.push_back(C); 202 } 203 else { 204 C=B; 205 } 206 if(e==0) { 207 break; 208 } 209 CGAL::internal::lickteig_roy_optimization<Polynomial_traits_d>(A,B,C,s,D); 210 B=D; 211 //typename Polynomial_traits_d::Pseudo_division() 212 // (A, -B, dummy_pol, D, dummy_nt); 213 //B= D / (CGAL::ipower(s,delta)*lcoeff(A)); 214 A=C; 215 s=lcoeff(A); 216 } 217 218 CGAL_assertion(static_cast<int>(sres.size()) 219 == degree(Q)+1); 220 221 // If P and Q were swapped, correct the signs 222 if(poly_swapped) { 223 int p = degree(P); 224 int q = degree(Q); 225 for(int i=0;i<=q;i++) { 226 if((p-i)*(q-i) % 2 == 1) { 227 sres[q-i]=-sres[q-i]; 228 } 229 } 230 } 231 232 // Now, reverse the entries 233 return std::copy(sres.rbegin(),sres.rend(),out); 234 } 235 236 237 /*! 238 * \brief Computes the polynomial subresultants 239 * as minors of the Bezout matrix 240 */ 241 template <typename Polynomial_traits_d,typename OutputIterator> inline bezout_polynomial_subresultants(typename Polynomial_traits_d::Polynomial_d P,typename Polynomial_traits_d::Polynomial_d Q,OutputIterator out)242 OutputIterator bezout_polynomial_subresultants 243 (typename Polynomial_traits_d::Polynomial_d P, 244 typename Polynomial_traits_d::Polynomial_d Q, 245 OutputIterator out) { 246 247 typedef typename Polynomial_traits_d::Polynomial_d Polynomial; 248 typedef typename Polynomial_traits_d::Coefficient_type NT; 249 typename Polynomial_traits_d::Leading_coefficient lcoeff; 250 typename Polynomial_traits_d::Degree degree; 251 typename Polynomial_traits_d::Construct_polynomial construct; 252 253 if(degree(P) < 1 || degree(Q) < 1) { 254 *out++ = Polynomial(CGAL::internal::resultant_for_constant_polynomial 255 <Polynomial_traits_d> (P,Q)); 256 return out; 257 } 258 259 typedef CGAL::internal::Simple_matrix<NT> Matrix; 260 Matrix M = CGAL::internal::polynomial_subresultant_matrix 261 <Polynomial_traits_d> (P,Q); 262 263 int r = static_cast<int>(M.row_dimension()); 264 265 for(int i = 0;i < r; i++) { 266 std::vector<NT> curr_row; 267 std::copy(M[r-1-i].begin(), 268 M[r-1-i].end(), 269 std::back_inserter(curr_row)); 270 //std::reverse(curr_row.begin(),curr_row.end()); 271 *out++ = construct(curr_row.rbegin(),curr_row.rend()); 272 } 273 int deg_diff=degree(P)-degree(Q); 274 if(deg_diff==0) { 275 *out++=Q; 276 } else if(deg_diff>0) { 277 *out++=CGAL::ipower(lcoeff(Q),deg_diff-1)*Q; 278 } else { 279 *out++=CGAL::ipower(lcoeff(P),-deg_diff-1)*P; 280 } 281 282 return out; 283 284 } 285 286 /*! 287 * \brief Compute the sequence of principal subresultants 288 * with pseudo-division 289 * 290 * Uses Ducos algorithm for the polynomial subresultant, and 291 * returns the formal leading coefficients. 292 */ 293 template <typename Polynomial_traits_d,typename OutputIterator> inline prs_principal_subresultants(typename Polynomial_traits_d::Polynomial_d P,typename Polynomial_traits_d::Polynomial_d Q,OutputIterator out)294 OutputIterator prs_principal_subresultants 295 (typename Polynomial_traits_d::Polynomial_d P, 296 typename Polynomial_traits_d::Polynomial_d Q, 297 OutputIterator out) { 298 299 typedef typename Polynomial_traits_d::Polynomial_d Polynomial; 300 typedef typename Polynomial_traits_d::Coefficient_type NT; 301 typename Polynomial_traits_d::Degree degree; 302 typename Polynomial_traits_d::Get_coefficient coeff; 303 304 std::vector<Polynomial> sres; 305 int q = (std::min)(degree(Q),degree(P)); 306 307 CGAL::internal::prs_polynomial_subresultants<Polynomial_traits_d> 308 (P,Q,std::back_inserter(sres)); 309 CGAL_assertion(static_cast<int>(sres.size()) == q+1); 310 for(int i=0; i <= q; i++) { 311 int d = degree(sres[i]); 312 CGAL_assertion(d<=i); 313 if(d<i) { 314 *out++ = NT(0); 315 } else { 316 *out++ = coeff(sres[i],i); 317 } 318 } 319 return out; 320 } 321 322 /*! 323 * \brief Compute the sequence of principal subresultants 324 * with minors of the Bezout matrix 325 * 326 */ 327 template <typename Polynomial_traits_d,typename OutputIterator> inline bezout_principal_subresultants(typename Polynomial_traits_d::Polynomial_d P,typename Polynomial_traits_d::Polynomial_d Q,OutputIterator out)328 OutputIterator bezout_principal_subresultants 329 (typename Polynomial_traits_d::Polynomial_d P, 330 typename Polynomial_traits_d::Polynomial_d Q, 331 OutputIterator out) { 332 333 typedef typename Polynomial_traits_d::Coefficient_type NT; 334 typename Polynomial_traits_d::Leading_coefficient lcoeff; 335 typename Polynomial_traits_d::Degree degree; 336 337 if(degree(P) < 1 || degree(Q) < 1) { 338 *out++ = CGAL::internal::resultant_for_constant_polynomial 339 <Polynomial_traits_d> (P,Q); 340 return out; 341 } 342 343 typedef CGAL::internal::Simple_matrix<NT> Matrix; 344 Matrix M = CGAL::internal::polynomial_subresultant_matrix 345 <Polynomial_traits_d> (P,Q,1); 346 347 int r = static_cast<int>(M.row_dimension()); 348 349 for(int i = r - 1;i >=0; i--) { 350 *out++=M[i][i]; 351 } 352 int deg_diff=degree(P)-degree(Q); 353 if(deg_diff==0) { 354 *out++=NT(1); 355 } else if(deg_diff>0) { 356 *out++=CGAL::ipower(lcoeff(Q),deg_diff); 357 } else { 358 *out++=CGAL::ipower(lcoeff(P),-deg_diff); 359 } 360 return out; 361 362 } 363 364 /*! 365 * \brief Computes the subresultants together with the according cofactors 366 * 367 * For details, see S.Basu, R.Pollack, M.-F.Roy: Algorithms in Real 368 * Algebraic Geometry, Second edition, Alg.8.22 369 */ 370 template<typename Polynomial_traits_d, 371 typename OutputIterator1, 372 typename OutputIterator2, 373 typename OutputIterator3> prs_subresultants_with_cofactors(typename Polynomial_traits_d::Polynomial_d P,typename Polynomial_traits_d::Polynomial_d Q,OutputIterator1 sres_out,OutputIterator2 coP_out,OutputIterator3 coQ_out)374 OutputIterator1 prs_subresultants_with_cofactors 375 (typename Polynomial_traits_d::Polynomial_d P, 376 typename Polynomial_traits_d::Polynomial_d Q, 377 OutputIterator1 sres_out, 378 OutputIterator2 coP_out, 379 OutputIterator3 coQ_out) { 380 381 typedef typename Polynomial_traits_d::Polynomial_d Polynomial; 382 typedef typename Polynomial_traits_d::Coefficient_type NT; 383 typename Polynomial_traits_d::Leading_coefficient lcoeff; 384 typename Polynomial_traits_d::Degree degree; 385 typename Polynomial_traits_d::Construct_polynomial construct; 386 387 if(degree(P) < 1 || degree(Q) < 1) { 388 *sres_out++ = Polynomial(CGAL::internal::resultant_for_constant_polynomial 389 <Polynomial_traits_d> (P,Q)); 390 *coP_out++ = Polynomial(lcoeff(Q)); 391 *coQ_out++ = Polynomial(lcoeff(P)); 392 return sres_out; 393 } 394 395 bool poly_swapped = (degree(P) < degree(Q)); 396 397 if(poly_swapped) { 398 std::swap(P,Q); 399 } 400 401 Polynomial zero_pol = construct(NT(0)); 402 std::vector<Polynomial> sres, coP, coQ; 403 404 #if 0 // old algorithm, there is some problem when deg_diff>1 405 406 int deg_diff=degree(P)-degree(Q); 407 408 409 410 if(deg_diff==0) { 411 sres.push_back(Q); 412 } else { 413 sres.push_back(CGAL::ipower(lcoeff(Q),deg_diff-1)*Q); 414 } 415 416 417 418 Polynomial A,B,C,D,Quo, coPA, coPB, coQA, coQB, coPC, coQC; 419 NT s,m; 420 int delta,d,e; 421 422 coPA = construct(NT(0)); 423 if(deg_diff==0) { 424 coQA = construct(NT(1)); 425 } else { 426 coQA = construct(CGAL::ipower(lcoeff(Q),deg_diff-1)); 427 } 428 429 coP.push_back(coPA); 430 coQ.push_back(coQA); 431 432 A=Q; 433 434 s=CGAL::ipower(lcoeff(Q),deg_diff); 435 //s=CGAL::ipower(lcoeff(Q),1); 436 437 typename Polynomial_traits_d::Pseudo_division() (P, -Q, Quo, B, m); 438 439 coPB = construct(m); 440 coQB = Quo; 441 442 //CGAL_assertion(m*P+Quo*Q==B); 443 //CGAL_assertion(CGAL::degree(B)<CGAL::degree(-Q)); 444 445 while(true) { 446 d=degree(A); 447 e=degree(B); 448 if(B.is_zero()) { 449 for(int i=0;i<d;i++) { 450 sres.push_back(zero_pol); 451 coP.push_back(zero_pol); 452 coQ.push_back(zero_pol); 453 } 454 break; 455 } 456 457 sres.push_back(B); 458 coP.push_back(coPB); 459 coQ.push_back(coQB); 460 461 //CGAL_assertion(coPB*P+coQB*Q==B); 462 463 delta=d-e; 464 if(delta>1) { 465 C=CGAL::ipower(lcoeff(B),delta-1)*B / CGAL::ipower(s,delta-1); 466 467 coPC = CGAL::ipower(lcoeff(B),delta-1)*coPB / 468 CGAL::ipower(s,delta-1); 469 coQC = CGAL::ipower(lcoeff(B),delta-1)*coQB / 470 CGAL::ipower(s,delta-1); 471 for(int i=0;i<delta-2;i++) { 472 sres.push_back(zero_pol); 473 coP.push_back(zero_pol); 474 coQ.push_back(zero_pol); 475 } 476 sres.push_back(C); 477 coP.push_back(coPC); 478 coQ.push_back(coQC); 479 480 } 481 else { 482 C=B; 483 coPC = coPB; 484 coQC = coQB; 485 } 486 if(e==0) { 487 break; 488 } 489 NT denominator = CGAL::ipower(s,delta)*lcoeff(A); 490 typename Polynomial_traits_d::Pseudo_division() (A, -B, Quo, D, m); 491 coPB = (m*coPA + Quo*coPB) / denominator; 492 coQB = (m*coQA + Quo*coQB) / denominator; 493 B = D / denominator; 494 A = C; 495 coPA = coPC; 496 coQA = coQC; 497 s = lcoeff(A); 498 } 499 500 501 #endif 502 503 int p = degree(P); 504 int q = degree(Q); 505 506 bool same_degree = (p==q); 507 508 if(same_degree) { 509 p++; 510 } 511 512 std::vector<Polynomial> sResP,sResU,sResV,C; 513 std::vector<NT> s,t; 514 515 for(int i=0;i<p+1;i++) { 516 sResP.push_back(construct(NT(0))); 517 sResU.push_back(construct(NT(0))); 518 sResV.push_back(construct(NT(0))); 519 C.push_back(construct(NT(0))); 520 s.push_back(NT(0)); 521 t.push_back(NT(0)); 522 } 523 sResP[p]=P; 524 s[p]=t[p]=(CGAL::sign(lcoeff(P))==CGAL::POSITIVE) ? NT(1) : NT(-1); 525 sResP[p-1]=Q; 526 t[p-1]=lcoeff(Q); 527 sResU[p]=sResV[p-1]=construct(NT(1)); 528 sResV[p]=sResU[p-1]=construct(NT(0)); 529 if(p-q>1) { 530 NT eps_p_minus_1 = ((p-q)%4==0 || (p-q)%4==1) ? NT(1) : NT(-1); 531 sResP[q]=eps_p_minus_1*CGAL::ipower(lcoeff(Q),p-q-1)*Q; 532 s[q]=eps_p_minus_1*CGAL::ipower(lcoeff(Q),p-q); 533 sResU[q]=construct(NT(0)); 534 sResV[q]=construct(eps_p_minus_1*CGAL::ipower(lcoeff(Q),p-q-1)); 535 for(int i=q+1;i<=p-2;i++) { 536 sResP[i]=sResU[i]=sResV[i]=construct(NT(0)); 537 s[i]=NT(0); 538 } 539 } 540 int i = p+1; 541 int j = p; 542 int k = 0; 543 while(!CGAL::is_zero(sResP[j-1])) { 544 k=degree(sResP[j-1]); 545 if(k>=j-1) { 546 if(k==0) { 547 break; 548 } 549 s[j-1]=t[j-1]; 550 NT prefac=CGAL::ipower(s[j-1],2); 551 NT denom=s[j]*t[i-1]; 552 Polynomial Quo,Rem; 553 NT D; 554 CGAL::pseudo_division(prefac*sResP[i-1],sResP[j-1],Quo,Rem,D); 555 C[k-1]=CGAL::integral_division(Quo,D); 556 sResP[k-1]=CGAL::integral_division 557 (-prefac*sResP[i-1]+C[k-1]*sResP[j-1], 558 denom); 559 sResU[k-1]=CGAL::integral_division 560 (-prefac*sResU[i-1]+C[k-1]*sResU[j-1], 561 denom); 562 sResV[k-1]=CGAL::integral_division 563 (-prefac*sResV[i-1]+C[k-1]*sResV[j-1], 564 denom); 565 } else { // k < j-1 566 567 s[j-1]=NT(0); 568 for(int delta=1;delta<=j-k-1;delta++) { 569 t[j-delta-1]=CGAL::ipower(NT(-1),delta)*CGAL::integral_division 570 (t[j-1]*t[j-delta],s[j]); 571 } 572 s[k]=t[k]; 573 sResP[k]=CGAL::integral_division(s[k]*sResP[j-1],t[j-1]); 574 sResU[k]=CGAL::integral_division(s[k]*sResU[j-1],t[j-1]); 575 sResV[k]=CGAL::integral_division(s[k]*sResV[j-1],t[j-1]); 576 for(int ell=k+1;ell<=j-2;ell++) { 577 sResP[ell]=sResU[ell]=sResV[ell]=construct(NT(0)); 578 s[ell]=NT(0); 579 } 580 if(k==0) { 581 break; 582 } 583 Polynomial Quo,Rem; 584 NT D; 585 NT prefac=s[k]*t[j-1]; 586 CGAL::pseudo_division(prefac*sResP[i-1],sResP[j-1],Quo,Rem,D); 587 C[k-1]=CGAL::integral_division(Quo,D); 588 589 NT denom = s[j]*t[i-1]; 590 sResP[k-1]=CGAL::integral_division 591 (-prefac*sResP[i-1]+C[k-1]*sResP[j-1],denom); 592 sResU[k-1]=CGAL::integral_division 593 (-prefac*sResU[i-1]+C[k-1]*sResU[j-1],denom); 594 sResV[k-1]=CGAL::integral_division 595 (-prefac*sResV[i-1]+C[k-1]*sResV[j-1],denom); 596 } 597 t[k-1]=lcoeff(sResP[k-1]); 598 i=j; 599 j=k; 600 } 601 if(k>0) { 602 for(int ell=0;ell<=j-2;ell++) { 603 sResP[ell]=sResU[ell]=sResV[ell]=construct(NT(0)); 604 s[ell]=NT(0); 605 } 606 } 607 608 // Correct factors for same degree (hack) 609 610 if(same_degree) { 611 for(int i = q-1;i>=0;i--) { 612 NT d = lcoeff(Q); 613 CGAL_assertion(CGAL::divides(d,sResP[i])); 614 sResP[i]=CGAL::integral_division(sResP[i],d); 615 CGAL_assertion(CGAL::divides(d,sResU[i])); 616 sResU[i]=CGAL::integral_division(sResU[i],d); 617 CGAL_assertion(CGAL::divides(d,sResV[i])); 618 sResV[i]=CGAL::integral_division(sResV[i],d); 619 } 620 } 621 622 623 // Correct the signs (the algorithm computes the signed subresultants) 624 625 if(degree(P)==degree(Q)) { 626 p--; 627 CGAL_assertion(p==q); 628 } 629 630 for(int i = q;i>=0;i--) { 631 if((p-i)%4==0 || (p-i)%4==1) { 632 sres.push_back(sResP[i]); 633 coP.push_back(sResU[i]); 634 coQ.push_back(sResV[i]); 635 } else { 636 sres.push_back(-sResP[i]); 637 coP.push_back(-sResU[i]); 638 coQ.push_back(-sResV[i]); 639 } 640 } 641 642 CGAL_assertion(static_cast<int>(sres.size()) 643 == degree(Q)+1); 644 645 // If P and Q were swapped, correct the signs 646 if(poly_swapped) { 647 int p = degree(P); 648 int q = degree(Q); 649 for(int i=0;i<=q;i++) { 650 if((p-i)*(q-i) % 2 == 1) { 651 sres[q-i] = -sres[q-i]; 652 coP[q-i] = -coP[q-i]; 653 coQ[q-i] = -coQ[q-i]; 654 } 655 } 656 for(int i=0;i<=q;i++) { 657 // Swap coP and coQ: 658 Polynomial help = coP[i]; 659 coP[i] = coQ[i]; 660 coQ[i] = help; 661 } 662 } 663 664 // Now, reverse the entries 665 std::copy(coP.rbegin(),coP.rend(),coP_out); 666 std::copy(coQ.rbegin(),coQ.rend(),coQ_out); 667 return std::copy(sres.rbegin(),sres.rend(),sres_out); 668 669 } 670 671 // the general function for CGAL::Integral_domain_without_division_tag 672 template <typename Polynomial_traits_d,typename OutputIterator> inline 673 OutputIterator polynomial_subresultants_(typename Polynomial_traits_d::Polynomial_d A,typename Polynomial_traits_d::Polynomial_d B,OutputIterator out,CGAL::Integral_domain_without_division_tag)674 polynomial_subresultants_(typename Polynomial_traits_d::Polynomial_d A, 675 typename Polynomial_traits_d::Polynomial_d B, 676 OutputIterator out, 677 CGAL::Integral_domain_without_division_tag){ 678 679 return bezout_polynomial_subresultants<Polynomial_traits_d>(A,B,out); 680 681 } 682 683 684 // the specialization for CGAL::Integral_domain_tag 685 template <typename Polynomial_traits_d,typename OutputIterator> inline 686 OutputIterator polynomial_subresultants_(typename Polynomial_traits_d::Polynomial_d A,typename Polynomial_traits_d::Polynomial_d B,OutputIterator out,CGAL::Integral_domain_tag)687 polynomial_subresultants_(typename Polynomial_traits_d::Polynomial_d A, 688 typename Polynomial_traits_d::Polynomial_d B, 689 OutputIterator out, 690 CGAL::Integral_domain_tag){ 691 692 return prs_polynomial_subresultants<Polynomial_traits_d>(A,B,out); 693 694 } 695 696 template <typename Polynomial_traits_d,typename OutputIterator> inline polynomial_subresultants_(typename Polynomial_traits_d::Polynomial_d A,typename Polynomial_traits_d::Polynomial_d B,OutputIterator out)697 OutputIterator polynomial_subresultants_ 698 (typename Polynomial_traits_d::Polynomial_d A, 699 typename Polynomial_traits_d::Polynomial_d B, 700 OutputIterator out) { 701 702 typedef typename Polynomial_traits_d::Coefficient_type NT; 703 704 typedef typename 705 CGAL::Algebraic_structure_traits<NT>::Algebraic_category 706 Algebraic_category; 707 return polynomial_subresultants_<Polynomial_traits_d> 708 (A,B,out,Algebraic_category()); 709 } 710 711 // the general function for CGAL::Integral_domain_without_division_tag 712 template <typename Polynomial_traits_d,typename OutputIterator> inline 713 OutputIterator principal_subresultants_(typename Polynomial_traits_d::Polynomial_d A,typename Polynomial_traits_d::Polynomial_d B,OutputIterator out,CGAL::Integral_domain_without_division_tag)714 principal_subresultants_(typename Polynomial_traits_d::Polynomial_d A, 715 typename Polynomial_traits_d::Polynomial_d B, 716 OutputIterator out, 717 CGAL::Integral_domain_without_division_tag){ 718 719 return bezout_principal_subresultants<Polynomial_traits_d>(A,B,out); 720 721 } 722 723 // the specialization for CGAL::Integral_domain_tag 724 template <typename Polynomial_traits_d,typename OutputIterator> inline 725 OutputIterator principal_subresultants_(typename Polynomial_traits_d::Polynomial_d A,typename Polynomial_traits_d::Polynomial_d B,OutputIterator out,CGAL::Integral_domain_tag)726 principal_subresultants_(typename Polynomial_traits_d::Polynomial_d A, 727 typename Polynomial_traits_d::Polynomial_d B, 728 OutputIterator out, 729 CGAL::Integral_domain_tag){ 730 731 return prs_principal_subresultants<Polynomial_traits_d>(A,B,out); 732 733 } 734 735 template <typename Polynomial_traits_d,typename OutputIterator> inline principal_subresultants_(typename Polynomial_traits_d::Polynomial_d A,typename Polynomial_traits_d::Polynomial_d B,OutputIterator out)736 OutputIterator principal_subresultants_ 737 (typename Polynomial_traits_d::Polynomial_d A, 738 typename Polynomial_traits_d::Polynomial_d B, 739 OutputIterator out) { 740 741 typedef typename Polynomial_traits_d::Coefficient_type NT; 742 743 typedef typename 744 CGAL::Algebraic_structure_traits<NT>::Algebraic_category 745 Algebraic_category; 746 return principal_subresultants_<Polynomial_traits_d> 747 (A,B,out,Algebraic_category()); 748 } 749 750 template<typename Polynomial_traits_d, 751 typename OutputIterator1, 752 typename OutputIterator2, 753 typename OutputIterator3> polynomial_subresultants_with_cofactors_(typename Polynomial_traits_d::Polynomial_d P,typename Polynomial_traits_d::Polynomial_d Q,OutputIterator1 sres_out,OutputIterator2 coP_out,OutputIterator3 coQ_out,CGAL::Integral_domain_tag)754 OutputIterator1 polynomial_subresultants_with_cofactors_ 755 (typename Polynomial_traits_d::Polynomial_d P, 756 typename Polynomial_traits_d::Polynomial_d Q, 757 OutputIterator1 sres_out, 758 OutputIterator2 coP_out, 759 OutputIterator3 coQ_out, 760 CGAL::Integral_domain_tag) { 761 return prs_subresultants_with_cofactors<Polynomial_traits_d> 762 (P,Q,sres_out,coP_out,coQ_out); 763 } 764 765 template<typename Polynomial_traits_d, 766 typename OutputIterator1, 767 typename OutputIterator2, 768 typename OutputIterator3> polynomial_subresultants_with_cofactors_(typename Polynomial_traits_d::Polynomial_d,typename Polynomial_traits_d::Polynomial_d,OutputIterator1 sres_out,OutputIterator2,OutputIterator3,CGAL::Integral_domain_without_division_tag)769 OutputIterator1 polynomial_subresultants_with_cofactors_ 770 (typename Polynomial_traits_d::Polynomial_d /* P */, 771 typename Polynomial_traits_d::Polynomial_d /* Q */, 772 OutputIterator1 sres_out, 773 OutputIterator2 /* coP_out */, 774 OutputIterator3 /* coQ_out */, 775 CGAL::Integral_domain_without_division_tag) { 776 // polynomial_subresultants_with_cofactors requires 777 // a model of IntegralDomain as coefficient type; 778 CGAL_static_assertion(sizeof(Polynomial_traits_d)==0); 779 return sres_out; 780 } 781 782 template<typename Polynomial_traits_d, 783 typename OutputIterator1, 784 typename OutputIterator2, 785 typename OutputIterator3> polynomial_subresultants_with_cofactors_(typename Polynomial_traits_d::Polynomial_d P,typename Polynomial_traits_d::Polynomial_d Q,OutputIterator1 sres_out,OutputIterator2 coP_out,OutputIterator3 coQ_out)786 OutputIterator1 polynomial_subresultants_with_cofactors_ 787 (typename Polynomial_traits_d::Polynomial_d P, 788 typename Polynomial_traits_d::Polynomial_d Q, 789 OutputIterator1 sres_out, 790 OutputIterator2 coP_out, 791 OutputIterator3 coQ_out) { 792 793 typedef typename Polynomial_traits_d::Coefficient_type NT; 794 795 typedef typename 796 CGAL::Algebraic_structure_traits<NT>::Algebraic_category 797 Algebraic_category; 798 return polynomial_subresultants_with_cofactors_<Polynomial_traits_d> 799 (P,Q,sres_out,coP_out,coQ_out,Algebraic_category()); 800 } 801 802 /*! \relates CGAL::Polynomial 803 * \brief compute the polynomial subresultants of the polynomials 804 * \c A and \c B 805 * 806 * If \c n and \c m are the degrees of p and q, 807 * the routine returns a sequence 808 * of length min(n,m)+1, the (polynomial) subresultants of \c p and \c q. 809 * It starts with the resultant of \c p and \c q. 810 * The <tt>i</tt>th polynomial has degree at most i. 811 * 812 * The way the subresultants are computed depends on the Algebra_type. 813 * In general the subresultant will be computed by the function 814 * CGAL::bezout_polynomial_subresultants, but if possible the function 815 * CGAL::prs_polynomial_subresultants is used. 816 */ 817 template <typename Polynomial_traits_d,typename OutputIterator> inline polynomial_subresultants(typename Polynomial_traits_d::Polynomial_d p,typename Polynomial_traits_d::Polynomial_d q,OutputIterator out)818 OutputIterator polynomial_subresultants 819 (typename Polynomial_traits_d::Polynomial_d p, 820 typename Polynomial_traits_d::Polynomial_d q, 821 OutputIterator out) { 822 return CGAL::internal::polynomial_subresultants_<Polynomial_traits_d> 823 (p, q, out); 824 } 825 826 /*! \relates CGAL::Polynomial 827 * \brief compute the principal subresultants of the polynomials 828 * \c p and \c q 829 * 830 * If \c n and \c m are the degrees of A and B, 831 * the routine returns a sequence 832 * of length min(n,m)+1, the (principal) subresultants of \c p and \c q, 833 * which starts with the resultant of \c p and \c q. 834 * 835 * The way the subresultants are computed depends on the Algebra_type. 836 * In general the subresultant will be computed by the function 837 * CGAL::bezout_principal_subresultants, but if possible the function 838 * CGAL::prs_principal_subresultants is used. 839 */ 840 template <typename Polynomial_traits_d,typename OutputIterator> inline principal_subresultants(typename Polynomial_traits_d::Polynomial_d p,typename Polynomial_traits_d::Polynomial_d q,OutputIterator out)841 OutputIterator principal_subresultants 842 (typename Polynomial_traits_d::Polynomial_d p, 843 typename Polynomial_traits_d::Polynomial_d q, 844 OutputIterator out) { 845 return CGAL::internal::principal_subresultants_<Polynomial_traits_d> 846 (p, q, out); 847 } 848 849 template<typename Polynomial_traits_d, 850 typename OutputIterator1, 851 typename OutputIterator2, 852 typename OutputIterator3> polynomial_subresultants_with_cofactors(typename Polynomial_traits_d::Polynomial_d p,typename Polynomial_traits_d::Polynomial_d q,OutputIterator1 sres_out,OutputIterator2 coP_out,OutputIterator3 coQ_out)853 OutputIterator1 polynomial_subresultants_with_cofactors 854 (typename Polynomial_traits_d::Polynomial_d p, 855 typename Polynomial_traits_d::Polynomial_d q, 856 OutputIterator1 sres_out, 857 OutputIterator2 coP_out, 858 OutputIterator3 coQ_out) { 859 return CGAL::internal::polynomial_subresultants_with_cofactors_ 860 <Polynomial_traits_d> (p,q,sres_out,coP_out,coQ_out); 861 } 862 863 } // namespace internal 864 865 } //namespace CGAL 866 867 #endif// CGAL_POLYNOMIAL_SUBRESULTANTS_H 868