Changeset c4f4fd in git for factory/cf_gcd_smallp.cc
 Timestamp:
 Sep 16, 2010, 5:34:13 PM (13 years ago)
 Branches:
 (u'spielwiese', 'e7cc1ebecb61be8b9ca6c18016352af89940b21a')
 Children:
 8e7db4cc8c23b1a45ca71f3217a02819e04d9c84
 Parents:
 7b8a4166957e36a8fb67484b68a9449846e1995d
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

factory/cf_gcd_smallp.cc
r7b8a416 rc4f4fd 30 30 #include "ftmpl_functions.h" 31 31 #include "cf_random.h" 32 #include "cf_algorithm.h" 33 #include "ffreval.h" 34 #include "fac_iterfor.h" 35 #include "cf_binom.h" 36 #include "facHensel.h" 37 #include "facFqFactorizeUtil.h" 38 39 // inline helper functions: 32 33 // iinline helper functions: 40 34 #include "cf_map_ext.h" 41 35 … … 52 46 TIMING_DEFINE_PRINT(newton_interpolation); 53 47 54 CanonicalForm55 gcd_univar_ntlp( const CanonicalForm & F, const CanonicalForm & G )56 {57 if (fac_NTL_char!=getCharacteristic())58 {59 fac_NTL_char=getCharacteristic();60 #ifdef NTL_ZZ61 ZZ r;62 r=getCharacteristic();63 ZZ_pContext ccc(r);64 #else65 zz_pContext ccc(getCharacteristic());66 #endif67 ccc.restore();68 #ifdef NTL_ZZ69 ZZ_p::init(r);70 #else71 zz_p::init(getCharacteristic());72 #endif73 }74 #ifdef NTL_ZZ75 ZZ_pX F1=convertFacCF2NTLZZpX(F);76 ZZ_pX G1=convertFacCF2NTLZZpX(G);77 ZZ_pX R=GCD(F1,G1);78 return convertNTLZZpX2CF(R,F.mvar());79 #else80 zz_pX F1=convertFacCF2NTLzzpX(F);81 zz_pX G1=convertFacCF2NTLzzpX(G);82 zz_pX R=GCD(F1,G1);83 return convertNTLzzpX2CF(R,F.mvar());84 #endif85 }86 87 #endif88 89 CanonicalForm90 gcd_poly_p( const CanonicalForm & f, const CanonicalForm & g)91 {92 CanonicalForm pi, pi1;93 CanonicalForm C, Ci, Ci1, Hi, bi, pi2;94 bool bpure;95 int delta = degree( f )  degree( g );96 97 if ( delta >= 0 )98 {99 pi = f; pi1 = g;100 }101 else102 {103 pi = g; pi1 = f; delta = delta;104 }105 Ci = content( pi ); Ci1 = content( pi1 );106 pi1 = pi1 / Ci1; pi = pi / Ci;107 C = gcd( Ci, Ci1 );108 if ( !( pi.isUnivariate() && pi1.isUnivariate() ) )109 {110 //out_cf("F:",f,"\n");111 //out_cf("G:",g,"\n");112 //out_cf("newGCD:",newGCD(f,g),"\n");113 if ( gcd_test_one( pi1, pi, true ) )114 {115 C=abs(C);116 //out_cf("GCD:",C,"\n");117 return C;118 }119 bpure = false;120 }121 else122 {123 bpure = isPurePoly(pi) && isPurePoly(pi1);124 #ifdef HAVE_NTL125 if ( isOn(SW_USE_NTL_GCD_P) && bpure && (CFFactory::gettype() != GaloisFieldDomain))126 return gcd_univar_ntlp(pi, pi1 ) * C;127 #endif128 }129 Variable v = f.mvar();130 Hi = power( LC( pi1, v ), delta );131 if ( (delta+1) % 2 )132 bi = 1;133 else134 bi = 1;135 while ( degree( pi1, v ) > 0 )136 {137 pi2 = psr( pi, pi1, v );138 pi2 = pi2 / bi;139 pi = pi1; pi1 = pi2;140 if ( degree( pi1, v ) > 0 )141 {142 delta = degree( pi, v )  degree( pi1, v );143 if ( (delta+1) % 2 )144 bi = LC( pi, v ) * power( Hi, delta );145 else146 bi = LC( pi, v ) * power( Hi, delta );147 Hi = power( LC( pi1, v ), delta ) / power( Hi, delta1 );148 }149 }150 if ( degree( pi1, v ) == 0 )151 {152 C=abs(C);153 //out_cf("GCD:",C,"\n");154 return C;155 }156 pi /= content( pi );157 if ( bpure )158 pi /= pi.lc();159 C=abs(C*pi);160 //out_cf("GCD:",C,"\n");161 return C;162 }163 164 #ifdef HAVE_NTL165 48 /// compressing two polynomials F and G, M is used for compressing, 166 49 /// N to reverse the compression 167 50 static inline 168 51 int myCompress (const CanonicalForm& F, const CanonicalForm& G, CFMap & M, 169 CFMap & N, const bool& topLevel)52 CFMap & N, bool& top_level) 170 53 { 171 54 int n= tmax (F.level(), G.level()); … … 184 67 int both_zero= 0; 185 68 186 if (top Level)69 if (top_level) 187 70 { 188 71 for (int i= 1; i <= n; i++) … … 315 198 delete [] degsg; 316 199 317 if (topLevel) 318 return both_non_zero; 319 else 320 return 1; 200 return 1; 321 201 } 322 323 CanonicalForm324 uniCoeff (const CanonicalForm& F)325 {326 if (F.inCoeffDomain())327 return F;328 if (F.isUnivariate() && F.level() == 1)329 return F;330 else if (F.isUnivariate () && F.level() != 1)331 return 0;332 CanonicalForm result= 0;333 if (F.level() == 2)334 {335 for (CFIterator i= F; i.hasTerms(); i++)336 result += i.coeff();337 }338 else339 {340 for (CFIterator i= F; i.hasTerms(); i++)341 result += uniCoeff (i.coeff());342 }343 return result;344 }345 346 int347 substituteCheck (const CanonicalForm& F, const CanonicalForm& G)348 {349 if (F.inBaseDomain()  G.inBaseDomain())350 return 0;351 Variable x= Variable (1);352 CanonicalForm f= swapvar (F, F.mvar(), x); //TODO swapping seems to be pretty expensive353 CanonicalForm g= swapvar (G, G.mvar(), x);354 int sizef= 0;355 int sizeg= 0;356 for (CFIterator i= f; i.hasTerms(); i++, sizef++)357 {358 if (i.exp() == 1)359 return 0;360 }361 for (CFIterator i= g; i.hasTerms(); i++, sizeg++)362 {363 if (i.exp() == 1)364 return 0;365 }366 int * expf= new int [sizef];367 int * expg= new int [sizeg];368 int j= 0;369 for (CFIterator i= f; i.hasTerms(); i++, j++)370 {371 expf [j]= i.exp();372 }373 j= 0;374 for (CFIterator i= g; i.hasTerms(); i++, j++)375 {376 expg [j]= i.exp();377 }378 379 int indf= sizef  1;380 int indg= sizeg  1;381 if (expf[indf] == 0)382 indf;383 if (expg[indg] == 0)384 indg;385 386 if (expg[indg] != expf [indf]  (expg[indg] == 1 && expf[indf] == 1))387 return 0;388 389 // expf[indg] == expf[indf] after here390 int result= expg[indg];391 for (int i= indf  1; i >= 0; i)392 {393 if (expf [i]%result != 0)394 return 0;395 }396 397 for (int i= indg  1; i >= 0; i)398 {399 if (expg [i]%result != 0)400 return 0;401 }402 403 delete [] expg;404 delete [] expf;405 return result;406 }407 408 // substiution409 void410 subst (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& A,411 CanonicalForm& B, const int d412 )413 {414 if (d == 1)415 {416 A= F;417 B= G;418 return;419 }420 421 CanonicalForm C= 0;422 CanonicalForm D= 0;423 Variable x= Variable (1);424 CanonicalForm f= swapvar (F, x, F.mvar());425 CanonicalForm g= swapvar (G, x, G.mvar());426 for (CFIterator i= f; i.hasTerms(); i++)427 C += i.coeff()*power (f.mvar(), i.exp()/ d);428 for (CFIterator i= g; i.hasTerms(); i++)429 D += i.coeff()*power (g.mvar(), i.exp()/ d);430 A= swapvar (C, x, F.mvar());431 B= swapvar (D, x, G.mvar());432 }433 434 CanonicalForm435 reverseSubst (const CanonicalForm& F, const int d)436 {437 if (d == 1)438 return F;439 440 Variable x= Variable (1);441 CanonicalForm f= swapvar (F, x, F.mvar());442 CanonicalForm result= 0;443 for (CFIterator i= f; i.hasTerms(); i++)444 result += i.coeff()*power (f.mvar(), d*i.exp());445 return swapvar (result, x, F.mvar());446 447 }448 449 static inline CanonicalForm450 uni_content (const CanonicalForm & F);451 452 CanonicalForm453 uni_content (const CanonicalForm& F, const Variable& x)454 {455 if (F.inBaseDomain())456 return F.genOne();457 if (F.level() == x.level() && F.isUnivariate())458 return F;459 if (F.level() != x.level() && F.isUnivariate())460 return F.genOne();461 462 if (x.level() != 1)463 {464 CanonicalForm f= swapvar (F, x, Variable (1));465 CanonicalForm result= uni_content (f);466 return swapvar (result, x, Variable (1));467 }468 else469 return uni_content (F);470 }471 202 472 203 /// compute the content of F, where F is considered as an element of … … 481 212 if (F.level() != 1 && F.isUnivariate()) 482 213 return F.genOne(); 483 if (degree (F, 214 if (degree (F,1) == 0) return F.genOne(); 484 215 485 216 int l= F.level(); … … 502 233 } 503 234 504 CanonicalForm505 extractContents (const CanonicalForm& F, const CanonicalForm& G,506 CanonicalForm& contentF, CanonicalForm& contentG,507 CanonicalForm& ppF, CanonicalForm& ppG, const int d)508 {509 CanonicalForm uniContentF, uniContentG, gcdcFcG;510 contentF= 1;511 contentG= 1;512 ppF= F;513 ppG= G;514 CanonicalForm result= 1;515 for (int i= 1; i <= d; i++)516 {517 uniContentF= uni_content (F, Variable (i));518 uniContentG= uni_content (G, Variable (i));519 gcdcFcG= gcd (uniContentF, uniContentG);520 contentF *= uniContentF;521 contentG *= uniContentG;522 ppF /= uniContentF;523 ppG /= uniContentG;524 result *= gcdcFcG;525 }526 return result;527 }528 529 235 /// compute the leading coefficient of F, where F is considered as an element 530 236 /// of \f$ R[x_{1}][x_{2},\ldots ,x_{n}] \f$, order on Mon (x_{2},\ldots ,x_{n}) 531 /// is Dp.237 /// is dp. 532 238 static inline 533 239 CanonicalForm uni_lcoeff (const CanonicalForm& F) … … 580 286 int p= getCharacteristic (); 581 287 int d= degree (mipo); 582 double bound= pow ((double)p, d);288 int bound= ipower (p, d); 583 289 do 584 290 { … … 620 326 ZZ_pX NTLirredpoly; 621 327 //TODO: replace d by max_{i} (deg_x{i}(f)) 622 int i= (int) (log ((double) pow ((double) d + 1, (double) num_vars))/log ((double) p));328 int i= (int) (log ((double) ipower (d + 1, num_vars))/log ((double) p)); 623 329 int m= degree (getMipo (beta)); 624 330 if (i <= 1) … … 630 336 631 337 /// GCD of F and G over \f$ F_{p}(\alpha ) \f$ , 632 /// l and top Level are only used internally, output is monic338 /// l and top_level are only used internally, output is monic 633 339 /// based on Alg. 7.2. as described in "Algorithms for 634 340 /// Computer Algebra" by Geddes, Czapor, Labahn 635 341 CanonicalForm 636 342 GCD_Fp_extension (const CanonicalForm& F, const CanonicalForm& G, 637 Variable & alpha, CFList& l, bool& top Level)343 Variable & alpha, CFList& l, bool& top_level) 638 344 { 639 345 CanonicalForm A= F; … … 648 354 649 355 CFMap M,N; 650 int best_level= myCompress (A, B, M, N, top Level);356 int best_level= myCompress (A, B, M, N, top_level); 651 357 652 358 if (best_level == 0) return B.genOne(); … … 660 366 if (A.isUnivariate() && B.isUnivariate()) 661 367 return N (gcd(A,B)); 662 663 int substitute= substituteCheck (A, B);664 368 665 if (substitute > 1)666 subst (A, B, A, B, substitute);667 668 369 CanonicalForm cA, cB; // content of A and B 669 370 CanonicalForm ppA, ppB; // primitive part of A and B 670 371 CanonicalForm gcdcAcB; 671 if (topLevel) 672 { 673 if (best_level <= 2) 674 gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, best_level); 675 else 676 gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, 2); // TODO needs some more tuning 677 } 678 else 679 { 680 cA = uni_content (A); 681 cB = uni_content (B); 682 gcdcAcB= gcd (cA, cB); 683 ppA= A/cA; 684 ppB= B/cB; 685 } 372 373 cA = uni_content (A); 374 cB = uni_content (B); 375 376 gcdcAcB= gcd (cA, cB); 377 378 ppA= A/cA; 379 ppB= B/cB; 686 380 687 381 CanonicalForm lcA, lcB; // leading coefficients of A and B … … 705 399 int d= totaldegree (ppA, Variable(2), Variable (ppA.level())); 706 400 707 if (d == 0) 708 { 709 if (substitute > 1) 710 return N (reverseSubst (gcdcAcB, substitute)); 711 else 712 return N(gcdcAcB); 713 } 401 if (d == 0) return N(gcdcAcB); 714 402 int d0= totaldegree (ppB, Variable(2), Variable (ppB.level())); 715 403 if (d0 < d) 716 404 d= d0; 717 if (d == 0) 718 { 719 if (substitute > 1) 720 return N (reverseSubst (gcdcAcB, substitute)); 721 else 722 return N(gcdcAcB); 723 } 405 if (d == 0) return N(gcdcAcB); 724 406 725 407 CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH; … … 731 413 H= 0; 732 414 bool fail= false; 733 top Level= false;415 top_level= false; 734 416 bool inextension= false; 735 417 Variable V_buf= alpha; … … 745 427 int num_vars= tmin (getNumVars(A), getNumVars(B)); 746 428 num_vars; 747 429 748 430 choose_extension (d, num_vars, V_buf); 749 431 bool prim_fail= false; 750 432 Variable V_buf2; 751 433 prim_elem= primitiveElement (alpha, V_buf2, prim_fail); 434 752 435 ASSERT (!prim_fail, "failure in integer factorizer"); 753 436 if (prim_fail) … … 755 438 else 756 439 im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf); 757 440 758 441 DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha)); 759 442 DEBOUTLN (cerr, "getMipo (V_buf2)= " << getMipo (V_buf2)); … … 778 461 G_random_element= 779 462 GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), V_buf, 780 list, top Level);463 list, top_level); 781 464 TIMING_END_AND_PRINT (gcd_recursion, 782 465 "time for recursive call: "); … … 789 472 G_random_element= 790 473 GCD_Fp_extension (ppA(random_element, x), ppB(random_element, x), V_buf, 791 list, top Level);474 list, top_level); 792 475 TIMING_END_AND_PRINT (gcd_recursion, 793 476 "time for recursive call: "); … … 797 480 d0= totaldegree (G_random_element, Variable(2), 798 481 Variable (G_random_element.level())); 799 if (d0 == 0) 800 { 801 if (substitute > 1) 802 return N (reverseSubst (gcdcAcB, substitute)); 803 else 804 return N(gcdcAcB); 805 } 482 if (d0 == 0) return N(gcdcAcB); 806 483 if (d0 > d) 807 484 { … … 844 521 ppH /= Lc(ppH); 845 522 DEBOUTLN (cerr, "ppH after mapDown= " << ppH); 846 if (fdivides (ppH, A) && fdivides (ppH, B)) 847 { 848 if (substitute > 1) 849 { 850 ppH= reverseSubst (ppH, substitute); 851 gcdcAcB= reverseSubst (gcdcAcB, substitute); 852 } 523 if (fdivides (ppH, A) && fdivides (ppH, B)) 853 524 return N(gcdcAcB*ppH); 854 } 855 } 856 else if (fdivides (ppH, A) && fdivides (ppH, B)) 857 { 858 if (substitute > 1) 859 { 860 ppH= reverseSubst (ppH, substitute); 861 gcdcAcB= reverseSubst (gcdcAcB, substitute); 862 } 863 return N(gcdcAcB*ppH); 864 } 525 } 526 else if (fdivides (ppH, A) && fdivides (ppH, B)) 527 return N(gcdcAcB*(ppH/Lc(ppH))); 865 528 } 866 529 … … 886 549 int p= getCharacteristic(); 887 550 int d= getGFDegree(); 888 double bound= pow ((double) p, (double)d);551 int bound= ipower (p, d); 889 552 do 890 553 { … … 917 580 /// the size of the base field is bounded by 2^16, output is monic 918 581 CanonicalForm GCD_GF (const CanonicalForm& F, const CanonicalForm& G, 919 CFList& l, bool& top Level)582 CFList& l, bool& top_level) 920 583 { 921 584 CanonicalForm A= F; … … 930 593 931 594 CFMap M,N; 932 int best_level= myCompress (A, B, M, N, top Level);595 int best_level= myCompress (A, B, M, N, top_level); 933 596 934 597 if (best_level == 0) return B.genOne(); … … 989 652 H= 0; 990 653 bool fail= false; 991 top Level= false;654 top_level= false; 992 655 bool inextension= false; 993 656 int p; … … 1004 667 num_vars; 1005 668 p= getCharacteristic(); 1006 expon= (int) floor ((log ((double) pow ((double)d + 1, num_vars))/log ((double) p)));669 expon= (int) floor ((log ((double) ipower (d + 1, num_vars))/log ((double) p))); 1007 670 if (expon < 2) 1008 671 expon= 2; 1009 672 kk= getGFDegree(); 1010 if ( pow ( (double)p, kk*expon) < (1 << 16))673 if (ipower (p, kk*expon) < (1 << 16)) 1011 674 setCharacteristic (p, kk*(int)expon, 'b'); 1012 675 else … … 1031 694 TIMING_START (gcd_recursion); 1032 695 G_random_element= GCD_GF (ppA(random_element, x), ppB(random_element, x), 1033 list, top Level);696 list, top_level); 1034 697 TIMING_END_AND_PRINT (gcd_recursion, 1035 698 "time for recursive call: "); … … 1041 704 TIMING_START (gcd_recursion); 1042 705 G_random_element= GCD_GF (ppA(random_element, x), ppB(random_element, x), 1043 list, top Level);706 list, top_level); 1044 707 TIMING_END_AND_PRINT (gcd_recursion, 1045 708 "time for recursive call: "); … … 1166 829 1167 830 CanonicalForm GCD_small_p (const CanonicalForm& F, const CanonicalForm& G, 1168 bool& top Level, CFList& l)831 bool& top_level, CFList& l) 1169 832 { 1170 833 CanonicalForm A= F; … … 1179 842 1180 843 CFMap M,N; 1181 int best_level= myCompress (A, B, M, N, top Level);844 int best_level= myCompress (A, B, M, N, top_level); 1182 845 1183 846 if (best_level == 0) return B.genOne(); … … 1185 848 A= M(A); 1186 849 B= M(B); 850 851 Variable x= Variable (1); 1187 852 1188 853 //univariate case 1189 854 if (A.isUnivariate() && B.isUnivariate()) 1190 855 return N (gcd (A, B)); 1191 1192 int substitute= substituteCheck (A, B);1193 1194 if (substitute > 1)1195 subst (A, B, A, B, substitute);1196 1197 Variable x= Variable (1);1198 856 1199 857 CanonicalForm cA, cB; // content of A and B 1200 858 CanonicalForm ppA, ppB; // primitive part of A and B 1201 859 CanonicalForm gcdcAcB; 1202 if (topLevel) 1203 { 1204 if (best_level <= 2) 1205 gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, best_level); 1206 else 1207 gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, 2); 1208 } 1209 else 1210 { 1211 cA = uni_content (A); 1212 cB = uni_content (B); 1213 gcdcAcB= gcd (cA, cB); 1214 ppA= A/cA; 1215 ppB= B/cB; 1216 } 860 cA = uni_content (A); 861 cB = uni_content (B); 862 gcdcAcB= gcd (cA, cB); 863 ppA= A/cA; 864 ppB= B/cB; 1217 865 1218 866 CanonicalForm lcA, lcB; // leading coefficients of A and B … … 1231 879 return G/Lc(G); 1232 880 } 1233 881 1234 882 gcdlcAlcB= gcd (lcA, lcB); 1235 883 … … 1237 885 int d0; 1238 886 1239 if (d == 0) 1240 { 1241 if (substitute > 1) 1242 return N (reverseSubst (gcdcAcB, substitute)); 1243 else 1244 return N(gcdcAcB); 1245 } 887 if (d == 0) return N(gcdcAcB); 1246 888 d0= totaldegree (ppB, Variable (2), Variable (ppB.level())); 1247 889 … … 1249 891 d= d0; 1250 892 1251 if (d == 0) 1252 { 1253 if (substitute > 1) 1254 return N (reverseSubst (gcdcAcB, substitute)); 1255 else 1256 return N(gcdcAcB); 1257 } 893 if (d == 0) return N(gcdcAcB); 1258 894 1259 895 CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH; … … 1263 899 G_m= 0; 1264 900 Variable alpha, V_buf; 901 int p, k; 1265 902 bool fail= false; 1266 903 bool inextension= false; 1267 904 bool inextensionextension= false; 1268 topLevel= false; 905 top_level= false; 906 CanonicalForm CF_buf; 1269 907 CFList source, dest; 908 CanonicalForm gcdcheck; 1270 909 do 1271 910 { … … 1280 919 TIMING_START (gcd_recursion); 1281 920 G_random_element= 1282 GCD_small_p (ppA (random_element,x), ppB (random_element,x), top Level,921 GCD_small_p (ppA (random_element,x), ppB (random_element,x), top_level, 1283 922 list); 1284 923 TIMING_END_AND_PRINT (gcd_recursion, … … 1292 931 G_random_element= 1293 932 GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), alpha, 1294 list, top Level);933 list, top_level); 1295 934 TIMING_END_AND_PRINT (gcd_recursion, 1296 935 "time for recursive call: "); … … 1303 942 CFList list; 1304 943 CanonicalForm mipo; 1305 int deg= 2;944 int deg= 3; 1306 945 do { 1307 l= CFList();1308 946 mipo= randomIrredpoly (deg, x); 1309 947 alpha= rootOf (mipo); … … 1312 950 random_element= randomElement (m, alpha, l, fail); 1313 951 deg++; 1314 } while (fail); 952 } while (fail); 1315 953 list= CFList(); 1316 954 TIMING_START (gcd_recursion); 1317 955 G_random_element= 1318 956 GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), alpha, 1319 list, top Level);957 list, top_level); 1320 958 TIMING_END_AND_PRINT (gcd_recursion, 1321 959 "time for recursive call: "); … … 1334 972 CanonicalForm prim_elem, im_prim_elem; 1335 973 prim_elem= primitiveElement (alpha, V_buf2, prim_fail); 974 1336 975 ASSERT (!prim_fail, "failure in integer factorizer"); 1337 976 if (prim_fail) … … 1362 1001 G_random_element= 1363 1002 GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), V_buf, 1364 list, top Level);1003 list, top_level); 1365 1004 TIMING_END_AND_PRINT (gcd_recursion, 1366 1005 "time for recursive call: "); … … 1374 1013 if (d0 == 0) 1375 1014 { 1376 if (substitute > 1) 1377 return N (reverseSubst (gcdcAcB, substitute)); 1378 else 1379 return N(gcdcAcB); 1015 return N(gcdcAcB); 1380 1016 } 1381 1017 if (d0 > d) … … 1413 1049 ppH /= Lc (ppH); 1414 1050 DEBOUTLN (cerr, "ppH= " << ppH); 1415 if (fdivides (ppH, A) && fdivides (ppH, B)) 1416 { 1417 if (substitute > 1) 1418 { 1419 ppH= reverseSubst (ppH, substitute); 1420 gcdcAcB= reverseSubst (gcdcAcB, substitute); 1421 } 1051 if (fdivides (ppH, A) && fdivides (ppH, B)) 1422 1052 return N(gcdcAcB*ppH); 1423 }1424 1053 } 1425 1054 … … 1432 1061 } 1433 1062 1434 CFArray1435 solveVandermonde (const CFArray& M, const CFArray& A)1436 {1437 int r= M.size();1438 ASSERT (A.size() == r, "vector does not have right size");1439 1440 if (r == 1)1441 {1442 CFArray result= CFArray (1);1443 result [0]= A [0] / M [0];1444 return result;1445 }1446 // check solvability1447 bool notDistinct= false;1448 for (int i= 0; i < r  1; i++)1449 {1450 for (int j= i + 1; j < r; j++)1451 {1452 if (M [i] == M [j])1453 {1454 notDistinct= true;1455 break;1456 }1457 }1458 }1459 if (notDistinct)1460 return CFArray();1461 1462 CanonicalForm master= 1;1463 Variable x= Variable (1);1464 for (int i= 0; i < r; i++)1465 master *= x  M [i];1466 CFList Pj;1467 CanonicalForm tmp;1468 for (int i= 0; i < r; i++)1469 {1470 tmp= master/(x  M [i]);1471 tmp /= tmp (M [i], 1);1472 Pj.append (tmp);1473 }1474 CFArray result= CFArray (r);1475 1476 CFListIterator j= Pj;1477 for (int i= 1; i <= r; i++, j++)1478 {1479 tmp= 0;1480 for (int l= 0; l < A.size(); l++)1481 tmp += A[l]*j.getItem()[l];1482 result[i  1]= tmp;1483 }1484 return result;1485 }1486 1487 /// solve M*x=L, where M is a Vandermonde matrix1488 CFArray1489 solveVandermonde (const CFMatrix& M, const CFArray& A)1490 {1491 int r= M.rows();1492 int c= M.columns();1493 ASSERT (c == r, "number of columns and rows not equal");1494 ASSERT (A.size() == r, "vector does not have right size");1495 1496 if (c == 1)1497 {1498 CFArray result= CFArray (1);1499 result [0]= A [0] / M (1,1);1500 return result;1501 }1502 // check solvability1503 CFMatrix secondRow= M (2, 2, 1, c);1504 bool notDistinct= false;1505 for (int i= 1; i < r; i++)1506 {1507 for (int j= i + 1; j <= r; j++)1508 {1509 if (secondRow (1, i) == secondRow (1, j))1510 {1511 notDistinct= true;1512 break;1513 }1514 }1515 }1516 if (notDistinct)1517 return CFArray();1518 1519 CanonicalForm master= 1;1520 Variable x= Variable (1);1521 for (int i= 1; i <= r; i++)1522 master *= x  secondRow (1, i);1523 CFList Pj;1524 CanonicalForm tmp;1525 for (int i= 1; i <= r; i++)1526 {1527 tmp= master/(x  secondRow (1, i));1528 tmp /= tmp (secondRow (1, i), 1);1529 Pj.append (tmp);1530 }1531 CFArray result= CFArray (r);1532 1533 CFListIterator j= Pj;1534 for (int i= 1; i <= r; i++, j++)1535 {1536 tmp= 0;1537 for (int l= 0; l < A.size(); l++)1538 tmp += A[l]*j.getItem()[l];1539 result[i  1]= tmp;1540 }1541 return result;1542 }1543 1544 CFArray1545 solveGeneralVandermonde (const CFArray& M, const CFArray& A)1546 {1547 int r= M.size();1548 ASSERT (c == r, "number of columns and rows not equal");1549 ASSERT (A.size() == r, "vector does not have right size");1550 if (r == 1)1551 {1552 CFArray result= CFArray (1);1553 result [0]= A[0] / M [0];1554 return result;1555 }1556 // check solvability1557 bool notDistinct= false;1558 for (int i= 0; i < r  1; i++)1559 {1560 for (int j= i + 1; j < r; j++)1561 {1562 if (M [i] == M [j])1563 {1564 notDistinct= true;1565 break;1566 }1567 }1568 }1569 if (notDistinct)1570 return CFArray();1571 1572 CanonicalForm master= 1;1573 Variable x= Variable (1);1574 for (int i= 0; i < r; i++)1575 master *= x  M [i];1576 master *= x;1577 CFList Pj;1578 CanonicalForm tmp;1579 for (int i= 0; i < r; i++)1580 {1581 tmp= master/(x  M [i]);1582 tmp /= tmp (M [i], 1);1583 Pj.append (tmp);1584 }1585 1586 CFArray result= CFArray (r);1587 1588 CFListIterator j= Pj;1589 for (int i= 1; i <= r; i++, j++)1590 {1591 tmp= 0;1592 1593 for (int l= 1; l <= A.size(); l++)1594 tmp += A[l  1]*j.getItem()[l];1595 result[i  1]= tmp;1596 }1597 return result;1598 }1599 1600 CFArray1601 solveGeneralVandermonde (const CFMatrix& M, const CFArray& A)1602 {1603 int r= M.rows();1604 int c= M.columns();1605 ASSERT (c == r, "number of columns and rows not equal");1606 ASSERT (A.size() == r, "vector does not have right size");1607 if (c == 1)1608 {1609 CFArray result= CFArray (1);1610 result [0]= A[0] / M (1,1);1611 return result;1612 }1613 // check solvability1614 CFMatrix firstRow= M (1, 1, 1, c);1615 bool notDistinct= false;1616 for (int i= 1; i < r; i++)1617 {1618 for (int j= i + 1; j <= r; j++)1619 {1620 if (firstRow (1, i) == firstRow (1, j))1621 {1622 notDistinct= true;1623 break;1624 }1625 }1626 }1627 if (notDistinct)1628 return CFArray();1629 1630 CanonicalForm master= 1;1631 Variable x= Variable (1);1632 for (int i= 1; i <= r; i++)1633 master *= x  firstRow (1, i);1634 master *= x;1635 CFList Pj;1636 CanonicalForm tmp;1637 for (int i= 1; i <= r; i++)1638 {1639 tmp= master/(x  firstRow (1, i));1640 tmp /= tmp (firstRow (1, i), 1);1641 Pj.append (tmp);1642 }1643 1644 CFArray result= CFArray (r);1645 CFListIterator j= Pj;1646 for (int i= 1; i <= r; i++, j++)1647 {1648 tmp= 0;1649 for (int l= 1; l <= A.size(); l++)1650 tmp += A[l  1]*j.getItem()[l];1651 result[i  1]= tmp;1652 }1653 return result;1654 }1655 1656 /// M in row echolon form, rk rank of M1657 CFArray1658 readOffSolution (const CFMatrix& M, const long rk)1659 {1660 CFArray result= CFArray (rk);1661 CanonicalForm tmp1, tmp2, tmp3;1662 for (int i= rk; i >= 1; i)1663 {1664 tmp3= 0;1665 tmp1= M (i, M.columns());1666 for (int j= M.columns()  1; j >= 1; j)1667 {1668 tmp2= M (i, j);1669 if (j == i)1670 break;1671 else1672 tmp3 += tmp2*result[j  1];1673 }1674 result[i  1]= (tmp1  tmp3)/tmp2;1675 }1676 return result;1677 }1678 1679 CFArray1680 readOffSolution (const CFMatrix& M, const CFArray& L, const CFArray& partialSol)1681 {1682 CFArray result= CFArray (M.rows());1683 CanonicalForm tmp1, tmp2, tmp3;1684 int k;1685 for (int i= M.rows(); i >= 1; i)1686 {1687 tmp3= 0;1688 tmp1= L[i  1];1689 k= 0;1690 for (int j= M.columns(); j >= 1; j, k++)1691 {1692 tmp2= M (i, j);1693 if (j == i)1694 break;1695 else1696 {1697 if (k > partialSol.size()  1)1698 tmp3 += tmp2*result[j  1];1699 else1700 tmp3 += tmp2*partialSol[partialSol.size()  k  1];1701 }1702 }1703 result [i  1]= (tmp1  tmp3)/tmp2;1704 }1705 return result;1706 }1707 1708 long1709 gaussianElimFp (CFMatrix& M, CFArray& L)1710 {1711 ASSERT (L.size() <= M.rows(), "dimension exceeded");1712 CFMatrix *N;1713 N= new CFMatrix (M.rows(), M.columns() + 1);1714 1715 for (int i= 1; i <= M.rows(); i++)1716 for (int j= 1; j <= M.columns(); j++)1717 (*N) (i, j)= M (i, j);1718 1719 int j= 1;1720 for (int i= 0; i < L.size(); i++, j++)1721 (*N) (j, M.columns() + 1)= L[i];1722 int p= getCharacteristic ();1723 zz_p::init (p);1724 mat_zz_p *NTLN= convertFacCFMatrix2NTLmat_zz_p(*N);1725 long rk= gauss (*NTLN);1726 1727 N= convertNTLmat_zz_p2FacCFMatrix (*NTLN);1728 1729 L= CFArray (M.rows());1730 for (int i= 0; i < M.rows(); i++)1731 L[i]= (*N) (i + 1, M.columns() + 1);1732 M= (*N) (1, M.rows(), 1, M.columns());1733 return rk;1734 }1735 1736 long1737 gaussianElimFq (CFMatrix& M, CFArray& L, const Variable& alpha)1738 {1739 ASSERT (L.size() <= M.rows(), "dimension exceeded");1740 CFMatrix *N;1741 N= new CFMatrix (M.rows(), M.columns() + 1);1742 1743 for (int i= 1; i <= M.rows(); i++)1744 for (int j= 1; j <= M.columns(); j++)1745 (*N) (i, j)= M (i, j);1746 1747 int j= 1;1748 for (int i= 0; i < L.size(); i++, j++)1749 (*N) (j, M.columns() + 1)= L[i];1750 int p= getCharacteristic ();1751 zz_p::init (p);1752 zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));1753 zz_pE::init (NTLMipo);1754 mat_zz_pE *NTLN= convertFacCFMatrix2NTLmat_zz_pE(*N);1755 long rk= gauss (*NTLN);1756 1757 N= convertNTLmat_zz_pE2FacCFMatrix (*NTLN, alpha);1758 1759 M= (*N) (1, M.rows(), 1, M.columns());1760 L= CFArray (M.rows());1761 for (int i= 0; i < M.rows(); i++)1762 L[i]= (*N) (i + 1, M.columns() + 1);1763 return rk;1764 }1765 1766 CFArray1767 solveSystemFp (const CFMatrix& M, const CFArray& L)1768 {1769 ASSERT (L.size() <= M.rows(), "dimension exceeded");1770 CFMatrix *N;1771 N= new CFMatrix (M.rows(), M.columns() + 1);1772 1773 for (int i= 1; i <= M.rows(); i++)1774 for (int j= 1; j <= M.columns(); j++)1775 (*N) (i, j)= M (i, j);1776 1777 int j= 1;1778 for (int i= 0; i < L.size(); i++, j++)1779 (*N) (j, M.columns() + 1)= L[i];1780 int p= getCharacteristic ();1781 zz_p::init (p);1782 mat_zz_p *NTLN= convertFacCFMatrix2NTLmat_zz_p(*N);1783 long rk= gauss (*NTLN);1784 if (rk != M.columns())1785 return CFArray();1786 1787 N= convertNTLmat_zz_p2FacCFMatrix (*NTLN);1788 1789 CFArray A= readOffSolution (*N, rk);1790 1791 return A;1792 }1793 1794 CFArray1795 solveSystemFq (const CFMatrix& M, const CFArray& L, const Variable& alpha)1796 {1797 ASSERT (L.size() <= M.rows(), "dimension exceeded");1798 CFMatrix *N;1799 N= new CFMatrix (M.rows(), M.columns() + 1);1800 1801 for (int i= 1; i <= M.rows(); i++)1802 for (int j= 1; j <= M.columns(); j++)1803 (*N) (i, j)= M (i, j);1804 int j= 1;1805 for (int i= 0; i < L.size(); i++, j++)1806 (*N) (j, M.columns() + 1)= L[i];1807 int p= getCharacteristic ();1808 zz_p::init (p);1809 zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));1810 zz_pE::init (NTLMipo);1811 mat_zz_pE *NTLN= convertFacCFMatrix2NTLmat_zz_pE(*N);1812 long rk= gauss (*NTLN);1813 if (rk != M.columns())1814 return CFArray();1815 1816 N= convertNTLmat_zz_pE2FacCFMatrix (*NTLN, alpha);1817 1818 CFArray A= readOffSolution (*N, rk);1819 1820 return A;1821 }1822 1823 CFArray1824 getMonoms (const CanonicalForm& F)1825 {1826 if (F.inCoeffDomain())1827 {1828 CFArray result= CFArray (1);1829 result [0]= 1;1830 return result;1831 }1832 if (F.isUnivariate())1833 {1834 CFArray result= CFArray (size(F));1835 int j= 0;1836 for (CFIterator i= F; i.hasTerms(); i++, j++)1837 result[j]= power (F.mvar(), i.exp());1838 return result;1839 }1840 int numMon= size (F);1841 CFArray result= CFArray (numMon);1842 int j= 0;1843 CFArray recResult;1844 Variable x= F.mvar();1845 CanonicalForm powX;1846 for (CFIterator i= F; i.hasTerms(); i++)1847 {1848 powX= power (x, i.exp());1849 recResult= getMonoms (i.coeff());1850 for (int k= 0; k < recResult.size(); k++)1851 result[j+k]= powX*recResult[k];1852 j += recResult.size();1853 }1854 return result;1855 }1856 1857 CFArray1858 evaluateMonom (const CanonicalForm& F, const CFList& evalPoints)1859 {1860 if (F.inCoeffDomain())1861 {1862 CFArray result= CFArray (1);1863 result [0]= F;1864 return result;1865 }1866 if (F.isUnivariate())1867 {1868 ASSERT (evalPoints.length() == 1,1869 "expected an eval point with only one component");1870 CFArray result= CFArray (size(F));1871 int j= 0;1872 CanonicalForm evalPoint= evalPoints.getLast();1873 for (CFIterator i= F; i.hasTerms(); i++, j++)1874 result[j]= power (evalPoint, i.exp());1875 return result;1876 }1877 int numMon= size (F);1878 CFArray result= CFArray (numMon);1879 int j= 0;1880 CanonicalForm evalPoint= evalPoints.getLast();1881 CFList buf= evalPoints;1882 buf.removeLast();1883 CFArray recResult;1884 CanonicalForm powEvalPoint;1885 for (CFIterator i= F; i.hasTerms(); i++)1886 {1887 powEvalPoint= power (evalPoint, i.exp());1888 recResult= evaluateMonom (i.coeff(), buf);1889 for (int k= 0; k < recResult.size(); k++)1890 result[j+k]= powEvalPoint*recResult[k];1891 j += recResult.size();1892 }1893 return result;1894 }1895 1896 CFArray1897 evaluate (const CFArray& A, const CFList& evalPoints)1898 {1899 CFArray result= A.size();1900 CanonicalForm tmp;1901 int k;1902 for (int i= 0; i < A.size(); i++)1903 {1904 tmp= A[i];1905 k= 1;1906 for (CFListIterator j= evalPoints; j.hasItem(); j++, k++)1907 tmp= tmp (j.getItem(), k);1908 result[i]= tmp;1909 }1910 return result;1911 }1912 1913 CFList1914 evaluationPoints (const CanonicalForm& F, const CanonicalForm& G,1915 CanonicalForm& Feval, CanonicalForm& Geval,1916 const CanonicalForm& LCF, const bool& GF,1917 const Variable& alpha, bool& fail, CFList& list1918 )1919 {1920 int k= tmax (F.level(), G.level())  1;1921 Variable x= Variable (1);1922 CFList result;1923 FFRandom genFF;1924 GFRandom genGF;1925 int p= getCharacteristic ();1926 double bound;1927 if (alpha != Variable (1))1928 {1929 bound= pow ((double) p, (double) degree (getMipo(alpha)));1930 bound= pow ((double) bound, (double) k);1931 }1932 else if (GF)1933 {1934 bound= pow ((double) p, (double) getGFDegree());1935 bound= pow ((double) bound, (double) k);1936 }1937 else1938 bound= pow ((double) p, (double) k);1939 1940 CanonicalForm random;1941 int j;1942 bool zeroOneOccured= false;1943 bool allEqual= false;1944 CanonicalForm buf;1945 do1946 {1947 random= 0;1948 // possible overflow if list.length() does not fit into a int1949 if (list.length() >= bound)1950 {1951 fail= true;1952 break;1953 }1954 for (int i= 0; i < k; i++)1955 {1956 if (GF)1957 {1958 result.append (genGF.generate());1959 random += result.getLast()*power (x, i);1960 }1961 else if (alpha != Variable(1))1962 {1963 AlgExtRandomF genAlgExt (alpha);1964 result.append (genAlgExt.generate());1965 random += result.getLast()*power (x, i);1966 }1967 else1968 {1969 result.append (genFF.generate());1970 random += result.getLast()*power (x, i);1971 }1972 if (result.getLast().isOne()  result.getLast().isZero())1973 zeroOneOccured= true;1974 }1975 if (find (list, random))1976 {1977 zeroOneOccured= false;1978 allEqual= false;1979 result= CFList();1980 continue;1981 }1982 if (zeroOneOccured)1983 {1984 list.append (random);1985 zeroOneOccured= false;1986 allEqual= false;1987 result= CFList();1988 continue;1989 }1990 // no zero at this point1991 if (k > 1)1992 {1993 allEqual= true;1994 CFIterator iter= random;1995 buf= iter.coeff();1996 iter++;1997 for (; iter.hasTerms(); iter++)1998 if (buf != iter.coeff())1999 allEqual= false;2000 }2001 if (allEqual)2002 {2003 list.append (random);2004 allEqual= false;2005 zeroOneOccured= false;2006 result= CFList();2007 continue;2008 }2009 2010 Feval= F;2011 Geval= G;2012 CanonicalForm LCeval= LCF;2013 j= 1;2014 for (CFListIterator i= result; i.hasItem(); i++, j++)2015 {2016 Feval= Feval (i.getItem(), j);2017 Geval= Geval (i.getItem(), j);2018 LCeval= LCeval (i.getItem(), j);2019 }2020 2021 if (LCeval.isZero())2022 {2023 if (!find (list, random))2024 list.append (random);2025 zeroOneOccured= false;2026 allEqual= false;2027 result= CFList();2028 continue;2029 }2030 2031 if (list.length() >= bound)2032 {2033 fail= true;2034 break;2035 }2036 } while (find (list, random));2037 2038 return result;2039 }2040 2041 /*CFList2042 evalPointsSmallChar (const CanonicalForm& F, const CanonicalForm& G,2043 CanonicalForm& Feval, CanonicalForm& Geval,2044 const CanonicalForm& LCF, const bool& GF,2045 const Variable& alpha, bool& fail, CFList& list2046 )2047 {2048 int k= tmax (F.level(), G.level())  1;2049 Variable x= Variable (1);2050 CFList result;2051 FFRandom genFF;2052 GFRandom genGF;2053 int p= getCharacteristic ();2054 double bound;2055 if (alpha != Variable (1))2056 {2057 bound= pow ((double) p, (double) degree (getMipo(alpha)));2058 bound= pow ((double) bound, (double) k);2059 }2060 else if (GF)2061 {2062 bound= pow ((double) p, (double) getGFDegree());2063 bound= pow ((double) bound, (double) k);2064 }2065 else2066 bound= pow ((double) p, (double) k);2067 2068 CanonicalForm random;2069 int j;2070 do2071 {2072 random= 0;2073 // possible overflow if list.length() does not fit into a int2074 if (list.length() >= bound)2075 {2076 fail= true;2077 break;2078 }2079 for (int i= 0; i < k; i++)2080 {2081 if (GF)2082 {2083 result.append (genGF.generate());2084 random += result.getLast()*power (x, i);2085 }2086 else if (alpha != Variable(1))2087 {2088 AlgExtRandomF genAlgExt (alpha);2089 result.append (genAlgExt.generate());2090 random += result.getLast()*power (x, i);2091 }2092 else2093 {2094 result.append (genFF.generate());2095 random += result.getLast()*power (x, i);2096 }2097 }2098 if (find (list, random))2099 {2100 result= CFList();2101 continue;2102 }2103 2104 Feval= F;2105 Geval= G;2106 CanonicalForm LCeval= LCF;2107 j= 1;2108 for (CFListIterator i= result; i.hasItem(); i++, j++)2109 {2110 Feval= Feval (i.getItem(), j);2111 Geval= Geval (i.getItem(), j);2112 LCeval= LCeval (i.getItem(), j);2113 }2114 2115 if (LCeval.isZero())2116 {2117 if (!find (list, random))2118 list.append (random);2119 result= CFList();2120 continue;2121 }2122 2123 if (list.length() >= bound)2124 {2125 fail= true;2126 break;2127 }2128 } while (find (list, random));2129 2130 return result;2131 }*/2132 2133 2134 /// multiply two lists componentwise2135 void mult (CFList& L1, const CFList& L2)2136 {2137 ASSERT (L1.length() == L2.length(), "lists of the same size expected");2138 2139 CFListIterator j= L2;2140 for (CFListIterator i= L1; i.hasItem(); i++, j++)2141 i.getItem() *= j.getItem();2142 }2143 2144 void eval (const CanonicalForm& A, const CanonicalForm& B, CanonicalForm& Aeval,2145 CanonicalForm& Beval, const CFList& L)2146 {2147 Aeval= A;2148 Beval= B;2149 int j= 1;2150 for (CFListIterator i= L; i.hasItem(); i++, j++)2151 {2152 Aeval= Aeval (i.getItem(), j);2153 Beval= Beval (i.getItem(), j);2154 }2155 }2156 2157 CanonicalForm2158 monicSparseInterpol (const CanonicalForm& F, const CanonicalForm& G,2159 const CanonicalForm& skeleton, const Variable& alpha,2160 bool& fail, CFArray*& coeffMonoms, CFArray& Monoms2161 )2162 {2163 CanonicalForm A= F;2164 CanonicalForm B= G;2165 if (F.isZero() && degree(G) > 0) return G/Lc(G);2166 else if (G.isZero() && degree (F) > 0) return F/Lc(F);2167 else if (F.isZero() && G.isZero()) return F.genOne();2168 if (F.inBaseDomain()  G.inBaseDomain()) return F.genOne();2169 if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);2170 if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);2171 if (F == G) return F/Lc(F);2172 2173 CFMap M,N;2174 int best_level= myCompress (A, B, M, N, false);2175 2176 if (best_level == 0)2177 return B.genOne();2178 2179 A= M(A);2180 B= M(B);2181 2182 Variable x= Variable (1);2183 ASSERT (degree (A, y) == 0, "expected degree (F, 1) == 0");2184 ASSERT (degree (B, y) == 0, "expected degree (G, 1) == 0");2185 2186 //univariate case2187 if (A.isUnivariate() && B.isUnivariate())2188 return N (gcd (A, B));2189 2190 CanonicalForm skel= M(skeleton);2191 CanonicalForm cA, cB; // content of A and B2192 CanonicalForm ppA, ppB; // primitive part of A and B2193 CanonicalForm gcdcAcB;2194 cA = uni_content (A);2195 cB = uni_content (B);2196 gcdcAcB= gcd (cA, cB);2197 ppA= A/cA;2198 ppB= B/cB;2199 2200 CanonicalForm lcA, lcB; // leading coefficients of A and B2201 CanonicalForm gcdlcAlcB;2202 lcA= uni_lcoeff (ppA);2203 lcB= uni_lcoeff (ppB);2204 2205 if (fdivides (lcA, lcB))2206 {2207 if (fdivides (A, B))2208 return F/Lc(F);2209 }2210 if (fdivides (lcB, lcA))2211 {2212 if (fdivides (B, A))2213 return G/Lc(G);2214 }2215 2216 gcdlcAlcB= gcd (lcA, lcB);2217 int skelSize= size (skel, skel.mvar());2218 2219 int j= 0;2220 int biggestSize= 0;2221 2222 for (CFIterator i= skel; i.hasTerms(); i++, j++)2223 biggestSize= tmax (biggestSize, size (i.coeff()));2224 2225 CanonicalForm g, Aeval, Beval;2226 2227 CFList evalPoints;2228 bool evalFail= false;2229 CFList list;2230 bool GF= false;2231 CanonicalForm LCA= LC (A);2232 CanonicalForm tmp;2233 CFArray gcds= CFArray (biggestSize);2234 CFList * pEvalPoints= new CFList [biggestSize];2235 Variable V_buf= alpha;2236 CFList source, dest;2237 CanonicalForm prim_elem, im_prim_elem;2238 for (int i= 0; i < biggestSize; i++)2239 {2240 if (i == 0)2241 evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf, evalFail,2242 list);2243 else2244 {2245 mult (evalPoints, pEvalPoints [0]);2246 eval (A, B, Aeval, Beval, evalPoints);2247 }2248 2249 if (evalFail)2250 {2251 if (V_buf != Variable (1))2252 {2253 do2254 {2255 int num_vars= tmin (getNumVars(A), getNumVars(B));2256 int d= totaldegree (A, Variable(2), Variable (A.level()));2257 d= tmin (d, totaldegree (B, Variable(2), Variable (B.level())));2258 Variable V_buf2= V_buf;2259 choose_extension (d, num_vars, V_buf2);2260 source= CFList();2261 dest= CFList();2262 2263 bool prim_fail= false;2264 Variable V_buf3;2265 prim_elem= primitiveElement (V_buf, V_buf3, prim_fail);2266 2267 ASSERT (!prim_fail, "failure in integer factorizer");2268 if (prim_fail)2269 ; //ERROR2270 else2271 im_prim_elem= mapPrimElem (prim_elem, V_buf, V_buf2);2272 2273 DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf));2274 DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));2275 2276 for (CFListIterator j= list; j.hasItem(); j++)2277 j.getItem()= mapUp (j.getItem(), V_buf, V_buf2, prim_elem,2278 im_prim_elem, source, dest);2279 for (int k= 0; k < i; k++)2280 {2281 for (CFListIterator j= pEvalPoints[k]; j.hasItem(); j++)2282 j.getItem()= mapUp (j.getItem(), V_buf, V_buf2, prim_elem,2283 im_prim_elem, source, dest);2284 gcds[k]= mapUp (gcds[k], V_buf, V_buf2, prim_elem, im_prim_elem,2285 source, dest);2286 }2287 2288 if (alpha != Variable (1))2289 {2290 A= mapUp (A, V_buf, V_buf2, prim_elem, im_prim_elem, source,dest);2291 B= mapUp (B, V_buf, V_buf2, prim_elem, im_prim_elem, source,dest);2292 }2293 evalFail= false;2294 evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,2295 evalFail, list);2296 } while (evalFail);2297 }2298 else2299 {2300 CanonicalForm mipo;2301 int deg= 2;2302 do {2303 mipo= randomIrredpoly (deg, x);2304 V_buf= rootOf (mipo);2305 evalFail= false;2306 evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,2307 evalFail, list);2308 deg++;2309 } while (evalFail);2310 }2311 }2312 2313 g= gcd (Aeval, Beval);2314 g /= Lc (g);2315 2316 if (degree (g) != degree (skel)  (skelSize != size (g)))2317 {2318 delete[] pEvalPoints;2319 fail= true;2320 return 0;2321 }2322 CFIterator l= skel;2323 for (CFIterator k= g; k.hasTerms(); k++, l++)2324 {2325 if (k.exp() != l.exp())2326 {2327 delete[] pEvalPoints;2328 fail= true;2329 return 0;2330 }2331 }2332 pEvalPoints[i]= evalPoints;2333 gcds[i]= g;2334 2335 tmp= 0;2336 int j= 0;2337 for (CFListIterator k= evalPoints; k.hasItem(); k++, j++)2338 tmp += k.getItem()*power (x, j);2339 list.append (tmp);2340 }2341 2342 if (Monoms.size() == 0)2343 Monoms= getMonoms (skel);2344 if (coeffMonoms == NULL)2345 coeffMonoms= new CFArray [skelSize];2346 j= 0;2347 for (CFIterator i= skel; i.hasTerms(); i++, j++)2348 coeffMonoms[j]= getMonoms (i.coeff());2349 2350 CFArray* pL= new CFArray [skelSize];2351 CFArray* pM= new CFArray [skelSize];2352 for (int i= 0; i < biggestSize; i++)2353 {2354 CFIterator l= gcds [i];2355 evalPoints= pEvalPoints [i];2356 for (int k= 0; k < skelSize; k++, l++)2357 {2358 if (i == 0)2359 pL[k]= CFArray (biggestSize);2360 pL[k] [i]= l.coeff();2361 2362 if (i == 0)2363 pM[k]= evaluate (coeffMonoms [k], evalPoints);2364 }2365 }2366 2367 CFArray solution;2368 CanonicalForm result= 0;2369 int ind= 0;2370 CFArray bufArray;2371 CFMatrix Mat;2372 for (int k= 0; k < skelSize; k++)2373 {2374 if (biggestSize != coeffMonoms[k].size())2375 {2376 bufArray= CFArray (coeffMonoms[k].size());2377 for (int i= 0; i < coeffMonoms[k].size(); i++)2378 bufArray [i]= pL[k] [i];2379 solution= solveGeneralVandermonde (pM [k], bufArray);2380 }2381 else2382 solution= solveGeneralVandermonde (pM [k], pL[k]);2383 2384 if (solution.size() == 0)2385 {2386 delete[] pEvalPoints;2387 delete[] pM;2388 delete[] pL;2389 delete[] coeffMonoms;2390 fail= true;2391 return 0;2392 }2393 for (int l= 0; l < solution.size(); l++)2394 result += solution[l]*Monoms [ind + l];2395 ind += solution.size();2396 }2397 2398 delete[] pEvalPoints;2399 delete[] pM;2400 delete[] pL;2401 2402 if (alpha != Variable (1) && V_buf != alpha)2403 {2404 CFList u, v;2405 result= mapDown (result, prim_elem, im_prim_elem, alpha, u, v);2406 }2407 2408 result= N(result);2409 if (fdivides (result, F) && fdivides (result, G))2410 return result;2411 else2412 {2413 delete[] coeffMonoms;2414 fail= true;2415 return 0;2416 }2417 }2418 2419 CanonicalForm2420 nonMonicSparseInterpol (const CanonicalForm& F, const CanonicalForm& G,2421 const CanonicalForm& skeleton, const Variable& alpha,2422 bool& fail, CFArray*& coeffMonoms, CFArray& Monoms2423 )2424 {2425 CanonicalForm A= F;2426 CanonicalForm B= G;2427 if (F.isZero() && degree(G) > 0) return G/Lc(G);2428 else if (G.isZero() && degree (F) > 0) return F/Lc(F);2429 else if (F.isZero() && G.isZero()) return F.genOne();2430 if (F.inBaseDomain()  G.inBaseDomain()) return F.genOne();2431 if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);2432 if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);2433 if (F == G) return F/Lc(F);2434 2435 CFMap M,N;2436 int best_level= myCompress (A, B, M, N, false);2437 2438 if (best_level == 0)2439 return B.genOne();2440 2441 A= M(A);2442 B= M(B);2443 2444 Variable x= Variable (1);2445 ASSERT (degree (A, y) == 0, "expected degree (F, 1) == 0");2446 ASSERT (degree (B, y) == 0, "expected degree (G, 1) == 0");2447 2448 //univariate case2449 if (A.isUnivariate() && B.isUnivariate())2450 return N (gcd (A, B));2451 2452 CanonicalForm skel= M(skeleton);2453 2454 CanonicalForm cA, cB; // content of A and B2455 CanonicalForm ppA, ppB; // primitive part of A and B2456 CanonicalForm gcdcAcB;2457 cA = uni_content (A);2458 cB = uni_content (B);2459 gcdcAcB= gcd (cA, cB);2460 ppA= A/cA;2461 ppB= B/cB;2462 2463 CanonicalForm lcA, lcB; // leading coefficients of A and B2464 CanonicalForm gcdlcAlcB;2465 lcA= uni_lcoeff (ppA);2466 lcB= uni_lcoeff (ppB);2467 2468 if (fdivides (lcA, lcB))2469 {2470 if (fdivides (A, B))2471 return F/Lc(F);2472 }2473 if (fdivides (lcB, lcA))2474 {2475 if (fdivides (B, A))2476 return G/Lc(G);2477 }2478 2479 gcdlcAlcB= gcd (lcA, lcB);2480 int skelSize= size (skel, skel.mvar());2481 2482 int j= 0;2483 int biggestSize= 0;2484 int bufSize;2485 int numberUni= 0;2486 for (CFIterator i= skel; i.hasTerms(); i++, j++)2487 {2488 bufSize= size (i.coeff());2489 biggestSize= tmax (biggestSize, bufSize);2490 numberUni += bufSize;2491 }2492 numberUni;2493 numberUni= (int) ceil ( (double) numberUni / (skelSize  1));2494 biggestSize= tmax (biggestSize , numberUni);2495 2496 numberUni= biggestSize + size (LC(skel))  1;2497 int biggestSize2= tmax (numberUni, biggestSize);2498 2499 CanonicalForm g, Aeval, Beval;2500 2501 CFList evalPoints;2502 CFArray coeffEval;2503 bool evalFail= false;2504 CFList list;2505 bool GF= false;2506 CanonicalForm LCA= LC (A);2507 CanonicalForm tmp;2508 CFArray gcds= CFArray (biggestSize);2509 CFList * pEvalPoints= new CFList [biggestSize];2510 Variable V_buf= alpha;2511 CFList source, dest;2512 CanonicalForm prim_elem, im_prim_elem;2513 for (int i= 0; i < biggestSize; i++)2514 {2515 if (i == 0)2516 {2517 if (getCharacteristic() > 3)2518 evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,2519 evalFail, list);2520 else2521 evalFail= true;2522 2523 if (evalFail)2524 {2525 if (V_buf != Variable (1))2526 {2527 do2528 {2529 int num_vars= tmin (getNumVars(A), getNumVars(B));2530 int d= totaldegree (A, Variable(2), Variable (A.level()));2531 d= tmin (d, totaldegree (B, Variable(2), Variable (B.level())));2532 Variable V_buf2= V_buf;2533 choose_extension (d, num_vars, V_buf2);2534 source= CFList();2535 dest= CFList();2536 2537 bool prim_fail= false;2538 Variable V_buf3;2539 prim_elem= primitiveElement (V_buf, V_buf3, prim_fail);2540 2541 ASSERT (!prim_fail, "failure in integer factorizer");2542 if (prim_fail)2543 ; //ERROR2544 else2545 im_prim_elem= mapPrimElem (prim_elem, V_buf, V_buf2);2546 2547 DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf));2548 DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));2549 2550 for (CFListIterator i= list; i.hasItem(); i++)2551 i.getItem()= mapUp (i.getItem(), V_buf, V_buf2, prim_elem,2552 im_prim_elem, source, dest);2553 evalFail= false;2554 evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,2555 evalFail, list);2556 } while (evalFail);2557 }2558 else2559 {2560 CanonicalForm mipo;2561 int deg= 2;2562 do {2563 mipo= randomIrredpoly (deg, x);2564 V_buf= rootOf (mipo);2565 evalFail= false;2566 evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,2567 evalFail, list);2568 deg++;2569 } while (evalFail);2570 }2571 }2572 }2573 else2574 {2575 mult (evalPoints, pEvalPoints[0]);2576 eval (A, B, Aeval, Beval, evalPoints);2577 }2578 2579 g= gcd (Aeval, Beval);2580 g /= Lc (g);2581 2582 if (degree (g) != degree (skel)  (skelSize != size (g)))2583 {2584 delete[] pEvalPoints;2585 fail= true;2586 return 0;2587 }2588 CFIterator l= skel;2589 for (CFIterator k= g; k.hasTerms(); k++, l++)2590 {2591 if (k.exp() != l.exp())2592 {2593 delete[] pEvalPoints;2594 fail= true;2595 return 0;2596 }2597 }2598 pEvalPoints[i]= evalPoints;2599 gcds[i]= g;2600 2601 tmp= 0;2602 int j= 0;2603 for (CFListIterator k= evalPoints; k.hasItem(); k++, j++)2604 tmp += k.getItem()*power (x, j);2605 list.append (tmp);2606 }2607 2608 if (Monoms.size() == 0)2609 Monoms= getMonoms (skel);2610 2611 if (coeffMonoms == NULL)2612 coeffMonoms= new CFArray [skelSize];2613 2614 j= 0;2615 for (CFIterator i= skel; i.hasTerms(); i++, j++)2616 coeffMonoms[j]= getMonoms (i.coeff());2617 2618 int minimalColumnsIndex;2619 if (skelSize > 1)2620 minimalColumnsIndex= 1;2621 else2622 minimalColumnsIndex= 0;2623 int minimalColumns;2624 2625 CFArray* pM= new CFArray [skelSize];2626 CFMatrix Mat;2627 // find the Matrix with minimal number of columns2628 for (int i= 0; i < skelSize; i++)2629 {2630 pM[i]= CFArray (coeffMonoms[i].size());2631 if (i == 1)2632 minimalColumns= coeffMonoms[i].size();2633 if (i > 1)2634 {2635 if (minimalColumns > coeffMonoms[i].size())2636 {2637 minimalColumns= coeffMonoms[i].size();2638 minimalColumnsIndex= i;2639 }2640 }2641 }2642 CFMatrix* pMat= new CFMatrix [2];2643 pMat[0]= CFMatrix (biggestSize, coeffMonoms[0].size());2644 pMat[1]= CFMatrix (biggestSize, coeffMonoms[minimalColumnsIndex].size());2645 CFArray* pL= new CFArray [skelSize];2646 for (int i= 0; i < biggestSize; i++)2647 {2648 CFIterator l= gcds [i];2649 evalPoints= pEvalPoints [i];2650 for (int k= 0; k < skelSize; k++, l++)2651 {2652 if (i == 0)2653 pL[k]= CFArray (biggestSize);2654 pL[k] [i]= l.coeff();2655 2656 if (i == 0 && k != 0 && k != minimalColumnsIndex)2657 pM[k]= evaluate (coeffMonoms[k], evalPoints);2658 else if (i == 0 && (k == 0  k == minimalColumnsIndex))2659 coeffEval= evaluate (coeffMonoms[k], evalPoints);2660 if (i > 0 && (k == 0  k == minimalColumnsIndex))2661 coeffEval= evaluate (coeffMonoms[k], evalPoints);2662 2663 if (k == 0)2664 {2665 if (pMat[k].rows() >= i + 1)2666 {2667 for (int ind= 1; ind < coeffEval.size() + 1; ind++)2668 pMat[k] (i + 1, ind)= coeffEval[ind  1];2669 }2670 }2671 if (k == minimalColumnsIndex)2672 {2673 if (pMat[1].rows() >= i + 1)2674 {2675 for (int ind= 1; ind < coeffEval.size() + 1; ind++)2676 pMat[1] (i + 1, ind)= coeffEval[ind  1];2677 }2678 }2679 }2680 }2681 2682 CFArray solution;2683 CanonicalForm result= 0;2684 int ind= 1;2685 int matRows, matColumns;2686 matRows= pMat[1].rows();2687 matColumns= pMat[0].columns()  1;2688 matColumns += pMat[1].columns();2689 2690 Mat= CFMatrix (matRows, matColumns);2691 for (int i= 1; i <= matRows; i++)2692 for (int j= 1; j <= pMat[1].columns(); j++)2693 Mat (i, j)= pMat[1] (i, j);2694 2695 for (int j= pMat[1].columns() + 1; j <= matColumns;2696 j++, ind++)2697 {2698 for (int i= 1; i <= matRows; i++)2699 Mat (i, j)= (pMat [0] (i, ind + 1))*pL[minimalColumnsIndex] [i  1];2700 }2701 2702 CFArray firstColumn= CFArray (pMat[0].rows());2703 for (int i= 0; i < pMat[0].rows(); i++)2704 firstColumn [i]= pMat[0] (i + 1, 1);2705 2706 CFArray bufArray= pL[minimalColumnsIndex];2707 2708 for (int i= 0; i < biggestSize; i++)2709 pL[minimalColumnsIndex] [i] *= firstColumn[i];2710 2711 CFMatrix bufMat= pMat[1];2712 pMat[1]= Mat;2713 2714 if (V_buf != x)2715 solution= solveSystemFq (pMat[1],2716 pL[minimalColumnsIndex], V_buf);2717 else2718 solution= solveSystemFp (pMat[1],2719 pL[minimalColumnsIndex]);2720 2721 if (solution.size() == 0) //Javadi's method failed, try deKleine, Monagan, Wittkopf instead2722 {2723 CFMatrix bufMat0= pMat[0];2724 delete [] pMat;2725 pMat= new CFMatrix [skelSize];2726 pL[minimalColumnsIndex]= bufArray;2727 CFList* bufpEvalPoints;2728 CFArray bufGcds;2729 if (biggestSize != biggestSize2)2730 {2731 bufpEvalPoints= pEvalPoints;2732 pEvalPoints= new CFList [biggestSize2];2733 bufGcds= gcds;2734 gcds= CFArray (biggestSize2);2735 for (int i= 0; i < biggestSize; i++)2736 {2737 pEvalPoints[i]= bufpEvalPoints [i];2738 gcds[i]= bufGcds[i];2739 }2740 for (int i= 0; i < biggestSize2  biggestSize; i++)2741 {2742 mult (evalPoints, pEvalPoints[0]);2743 eval (A, B, Aeval, Beval, evalPoints);2744 g= gcd (Aeval, Beval);2745 g /= Lc (g);2746 if (degree (g) != degree (skel)  (skelSize != size (g)))2747 {2748 delete[] pEvalPoints;2749 delete[] pMat;2750 delete[] pL;2751 delete[] coeffMonoms;2752 delete[] pM;2753 fail= true;2754 return 0;2755 }2756 CFIterator l= skel;2757 for (CFIterator k= g; k.hasTerms(); k++, l++)2758 {2759 if (k.exp() != l.exp())2760 {2761 delete[] pEvalPoints;2762 delete[] pMat;2763 delete[] pL;2764 delete[] coeffMonoms;2765 delete[] pM;2766 fail= true;2767 return 0;2768 }2769 }2770 pEvalPoints[i + biggestSize]= evalPoints;2771 gcds[i + biggestSize]= g;2772 }2773 }2774 for (int i= 0; i < biggestSize; i++)2775 {2776 CFIterator l= gcds [i];2777 evalPoints= pEvalPoints [i];2778 for (int k= 1; k < skelSize; k++, l++)2779 {2780 if (i == 0)2781 pMat[k]= CFMatrix (biggestSize2, coeffMonoms[k].size() + biggestSize2  1);2782 if (k == minimalColumnsIndex)2783 continue;2784 coeffEval= evaluate (coeffMonoms[k], evalPoints);2785 if (pMat[k].rows() >= i + 1)2786 {2787 for (int ind= 1; ind < coeffEval.size() + 1; ind++)2788 pMat[k] (i + 1, ind)= coeffEval[ind  1];2789 }2790 }2791 }2792 Mat= bufMat0;2793 pMat[0]= CFMatrix (biggestSize2, coeffMonoms[0].size() + biggestSize2  1);2794 for (int j= 1; j <= Mat.rows(); j++)2795 for (int k= 1; k <= Mat.columns(); k++)2796 pMat [0] (j,k)= Mat (j,k);2797 Mat= bufMat;2798 for (int j= 1; j <= Mat.rows(); j++)2799 for (int k= 1; k <= Mat.columns(); k++)2800 pMat [minimalColumnsIndex] (j,k)= Mat (j,k);2801 // write old matrix entries into new matrices2802 for (int i= 0; i < skelSize; i++)2803 {2804 bufArray= pL[i];2805 pL[i]= CFArray (biggestSize2);2806 for (int j= 0; j < bufArray.size(); j++)2807 pL[i] [j]= bufArray [j];2808 }2809 //write old vector entries into new and add new entries to old matrices2810 for (int i= 0; i < biggestSize2  biggestSize; i++)2811 {2812 CFIterator l= gcds [i + biggestSize];2813 evalPoints= pEvalPoints [i + biggestSize];2814 for (int k= 0; k < skelSize; k++, l++)2815 {2816 pL[k] [i + biggestSize]= l.coeff();2817 coeffEval= evaluate (coeffMonoms[k], evalPoints);2818 if (pMat[k].rows() >= i + biggestSize + 1)2819 {2820 for (int ind= 1; ind < coeffEval.size() + 1; ind++)2821 pMat[k] (i + biggestSize + 1, ind)= coeffEval[ind  1];2822 }2823 }2824 }2825 // begin new2826 for (int i= 0; i < skelSize; i++)2827 {2828 if (pL[i].size() > 1)2829 {2830 for (int j= 2; j <= pMat[i].rows(); j++)2831 pMat[i] (j, coeffMonoms[i].size() + j  1)=2832 pL[i] [j  1];2833 }2834 }2835 2836 2837 long rk;2838 matColumns= biggestSize2  1;2839 matRows= 0;2840 for (int i= 0; i < skelSize; i++)2841 {2842 if (V_buf == x)2843 rk= gaussianElimFp (pMat[i], pL[i]);2844 else2845 rk= gaussianElimFq (pMat[i], pL[i], V_buf);2846 2847 if (pMat[i] (coeffMonoms[i].size(), coeffMonoms[i].size()) == 0)2848 {2849 delete[] pEvalPoints;2850 delete[] pMat;2851 delete[] pL;2852 delete[] coeffMonoms;2853 delete[] pM;2854 fail= true;2855 return 0;2856 }2857 matRows += pMat[i].rows()  coeffMonoms[i].size();2858 }2859 2860 CFMatrix bufMat;2861 Mat= CFMatrix (matRows, matColumns);2862 ind= 0;2863 bufArray= CFArray (matRows);2864 CFArray bufArray2;2865 for (int i= 0; i < skelSize; i++)2866 {2867 bufMat= pMat[i] (coeffMonoms[i].size() + 1, pMat[i].rows(),2868 coeffMonoms[i].size() + 1, pMat[i].columns());2869 2870 for (int j= 1; j <= bufMat.rows(); j++)2871 for (int k= 1; k <= bufMat.columns(); k++)2872 Mat (j + ind, k)= bufMat(j, k);2873 bufArray2= coeffMonoms[i].size();2874 for (int j= 1; j <= pMat[i].rows(); j++)2875 {2876 if (j > coeffMonoms[i].size())2877 bufArray [jcoeffMonoms[i].size() + ind  1]= pL[i] [j  1];2878 else2879 bufArray2 [j  1]= pL[i] [j  1];2880 }2881 pL[i]= bufArray2;2882 ind += bufMat.rows();2883 pMat[i]= pMat[i] (1, coeffMonoms[i].size(), 1, pMat[i].columns());2884 }2885 2886 if (V_buf != x)2887 solution= solveSystemFq (Mat, bufArray, V_buf);2888 else2889 solution= solveSystemFp (Mat, bufArray);2890 2891 if (solution.size() == 0)2892 {2893 delete[] pEvalPoints;2894 delete[] pMat;2895 delete[] pL;2896 delete[] coeffMonoms;2897 delete[] pM;2898 fail= true;2899 return 0;2900 }2901 2902 ind= 0;2903 result= 0;2904 CFArray bufSolution;2905 for (int i= 0; i < skelSize; i++)2906 {2907 bufSolution= readOffSolution (pMat[i], pL[i], solution);2908 for (int i= 0; i < bufSolution.size(); i++, ind++)2909 result += Monoms [ind]*bufSolution[i];2910 }2911 if (alpha != Variable (1) && V_buf != alpha)2912 {2913 CFList u, v;2914 result= mapDown (result, prim_elem, im_prim_elem, alpha, u, v);2915 }2916 result= N(result);2917 if (fdivides (result, F) && fdivides (result, G))2918 return result;2919 else2920 {2921 fail= true;2922 return 0;2923 }2924 } // end of deKleine, Monagan & Wittkopf2925 2926 result += Monoms[0];2927 int ind2= 0, ind3= 2;2928 ind= 0;2929 for (int l= 0; l < minimalColumnsIndex; l++)2930 ind += coeffMonoms[l].size();2931 for (int l= solution.size()  pMat[0].columns() + 1; l < solution.size();2932 l++, ind2++, ind3++)2933 {2934 result += solution[l]*Monoms [1 + ind2];2935 for (int i= 0; i < pMat[0].rows(); i++)2936 firstColumn[i] += solution[l]*pMat[0] (i + 1, ind3);2937 }2938 for (int l= 0; l < solution.size() + 1  pMat[0].columns(); l++)2939 result += solution[l]*Monoms [ind + l];2940 ind= coeffMonoms[0].size();2941 for (int k= 1; k < skelSize; k++)2942 {2943 if (k == minimalColumnsIndex)2944 {2945 ind += coeffMonoms[k].size();2946 continue;2947 }2948 if (k != minimalColumnsIndex)2949 {2950 for (int i= 0; i < biggestSize; i++)2951 pL[k] [i] *= firstColumn [i];2952 }2953 2954 if (biggestSize != coeffMonoms[k].size() && k != minimalColumnsIndex)2955 {2956 bufArray= CFArray (coeffMonoms[k].size());2957 for (int i= 0; i < bufArray.size(); i++)2958 bufArray [i]= pL[k] [i];2959 solution= solveGeneralVandermonde (pM[k], bufArray);2960 }2961 else2962 solution= solveGeneralVandermonde (pM[k], pL[k]);2963 2964 if (solution.size() == 0)2965 {2966 delete[] pEvalPoints;2967 delete[] pMat;2968 delete[] pL;2969 delete[] coeffMonoms;2970 delete[] pM;2971 fail= true;2972 return 0;2973 }2974 if (k != minimalColumnsIndex)2975 {2976 for (int l= 0; l < solution.size(); l++)2977 result += solution[l]*Monoms [ind + l];2978 ind += solution.size();2979 }2980 }2981 2982 delete[] pEvalPoints;2983 delete[] pMat;2984 delete[] pL;2985 delete[] pM;2986 2987 if (alpha != Variable (1) && V_buf != alpha)2988 {2989 CFList u, v;2990 result= mapDown (result, prim_elem, im_prim_elem, alpha, u, v);2991 }2992 result= N(result);2993 2994 if (fdivides (result, F) && fdivides (result, G))2995 return result;2996 else2997 {2998 delete[] coeffMonoms;2999 fail= true;3000 return 0;3001 }3002 }3003 3004 CanonicalForm sparseGCDFq (const CanonicalForm& F, const CanonicalForm& G,3005 const Variable & alpha, CFList& l, bool& topLevel)3006 {3007 CanonicalForm A= F;3008 CanonicalForm B= G;3009 if (F.isZero() && degree(G) > 0) return G/Lc(G);3010 else if (G.isZero() && degree (F) > 0) return F/Lc(F);3011 else if (F.isZero() && G.isZero()) return F.genOne();3012 if (F.inBaseDomain()  G.inBaseDomain()) return F.genOne();3013 if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);3014 if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);3015 if (F == G) return F/Lc(F);3016 3017 CFMap M,N;3018 int best_level= myCompress (A, B, M, N, topLevel);3019 3020 if (best_level == 0) return B.genOne();3021 3022 A= M(A);3023 B= M(B);3024 3025 Variable x= Variable (1);3026 3027 //univariate case3028 if (A.isUnivariate() && B.isUnivariate())3029 return N (gcd (A, B));3030 3031 int substitute= substituteCheck (A, B);3032 3033 if (substitute > 1)3034 subst (A, B, A, B, substitute);3035 3036 CanonicalForm cA, cB; // content of A and B3037 CanonicalForm ppA, ppB; // primitive part of A and B3038 CanonicalForm gcdcAcB;3039 if (topLevel)3040 {3041 if (best_level <= 2)3042 gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, best_level);3043 else3044 gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, 2);3045 }3046 else3047 {3048 cA = uni_content (A);3049 cB = uni_content (B);3050 gcdcAcB= gcd (cA, cB);3051 ppA= A/cA;3052 ppB= B/cB;3053 }3054 3055 CanonicalForm lcA, lcB; // leading coefficients of A and B3056 CanonicalForm gcdlcAlcB;3057 lcA= uni_lcoeff (ppA);3058 lcB= uni_lcoeff (ppB);3059 3060 if (fdivides (lcA, lcB))3061 {3062 if (fdivides (A, B))3063 return F/Lc(F);3064 }3065 if (fdivides (lcB, lcA))3066 {3067 if (fdivides (B, A))3068 return G/Lc(G);3069 }3070 3071 gcdlcAlcB= gcd (lcA, lcB);3072 3073 int d= totaldegree (ppA, Variable (2), Variable (ppA.level()));3074 int d0;3075 3076 if (d == 0)3077 {3078 if (substitute > 1)3079 return N(reverseSubst (gcdcAcB, substitute));3080 else3081 return N(gcdcAcB);3082 }3083 d0= totaldegree (ppB, Variable (2), Variable (ppB.level()));3084 3085 if (d0 < d)3086 d= d0;3087 3088 if (d == 0)3089 {3090 if (substitute > 1)3091 return N(reverseSubst (gcdcAcB, substitute));3092 else3093 return N(gcdcAcB);3094 }3095 3096 CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH, skeleton;3097 CanonicalForm newtonPoly= 1;3098 m= gcdlcAlcB;3099 G_m= 0;3100 H= 0;3101 bool fail= false;3102 topLevel= false;3103 bool inextension= false;3104 Variable V_buf= alpha;3105 CanonicalForm prim_elem, im_prim_elem;3106 CFList source, dest;3107 do // first do3108 {3109 random_element= randomElement (m, V_buf, l, fail);3110 if (random_element == 0 && !fail)3111 {3112 if (!find (l, random_element))3113 l.append (random_element);3114 continue;3115 }3116 if (fail)3117 {3118 source= CFList();3119 dest= CFList();3120 int num_vars= tmin (getNumVars(A), getNumVars(B));3121 num_vars;3122 3123 choose_extension (d, num_vars, V_buf);3124 bool prim_fail= false;3125 Variable V_buf2;3126 prim_elem= primitiveElement (alpha, V_buf2, prim_fail);3127 3128 ASSERT (!prim_fail, "failure in integer factorizer");3129 if (prim_fail)3130 ; //ERROR3131 else3132 im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);3133 3134 DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));3135 DEBOUTLN (cerr, "getMipo (V_buf2)= " << getMipo (V_buf2));3136 inextension= true;3137 for (CFListIterator i= l; i.hasItem(); i++)3138 i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,3139 im_prim_elem, source, dest);3140 m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);3141 G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);3142 newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,3143 source, dest);3144 ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);3145 ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);3146 gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,3147 source, dest);3148 3149 fail= false;3150 random_element= randomElement (m, V_buf, l, fail );3151 DEBOUTLN (cerr, "fail= " << fail);3152 CFList list;3153 TIMING_START (gcd_recursion);3154 G_random_element=3155 sparseGCDFq (ppA (random_element, x), ppB (random_element, x), V_buf,3156 list, topLevel);3157 TIMING_END_AND_PRINT (gcd_recursion,3158 "time for recursive call: ");3159 DEBOUTLN (cerr, "G_random_element= " << G_random_element);3160 }3161 else3162 {3163 CFList list;3164 TIMING_START (gcd_recursion);3165 G_random_element=3166 sparseGCDFq (ppA(random_element, x), ppB(random_element, x), V_buf,3167 list, topLevel);3168 TIMING_END_AND_PRINT (gcd_recursion,3169 "time for recursive call: ");3170 DEBOUTLN (cerr, "G_random_element= " << G_random_element);3171 }3172 3173 d0= totaldegree (G_random_element, Variable(2),3174 Variable (G_random_element.level()));3175 if (d0 == 0)3176 {3177 if (substitute > 1)3178 return N(reverseSubst (gcdcAcB, substitute));3179 else3180 return N(gcdcAcB);3181 }3182 if (d0 > d)3183 {3184 if (!find (l, random_element))3185 l.append (random_element);3186 continue;3187 }3188 3189 G_random_element=3190 (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))3191 * G_random_element;3192 3193 skeleton= G_random_element;3194 d0= totaldegree (G_random_element, Variable(2),3195 Variable(G_random_element.level()));3196 if (d0 < d)3197 {3198 m= gcdlcAlcB;3199 newtonPoly= 1;3200 G_m= 0;3201 d= d0;3202 }3203 3204 H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);3205 if (uni_lcoeff (H) == gcdlcAlcB)3206 {3207 cH= uni_content (H);3208 ppH= H/cH;3209 if (inextension)3210 {3211 CFList u, v;3212 //maybe it's better to test if ppH is an element of F(\alpha) before3213 //mapping down3214 DEBOUTLN (cerr, "ppH before mapDown= " << ppH);3215 ppH= mapDown (ppH, prim_elem, im_prim_elem, alpha, u, v);3216 ppH /= Lc(ppH);3217 DEBOUTLN (cerr, "ppH after mapDown= " << ppH);3218 if (fdivides (ppH, A) && fdivides (ppH, B))3219 {3220 if (substitute > 1)3221 {3222 ppH= reverseSubst (ppH, substitute);3223 gcdcAcB= reverseSubst (gcdcAcB, substitute);3224 }3225 return N(gcdcAcB*ppH);3226 }3227 }3228 else if (fdivides (ppH, A) && fdivides (ppH, B))3229 {3230 if (substitute > 1)3231 {3232 ppH= reverseSubst (ppH, substitute);3233 gcdcAcB= reverseSubst (gcdcAcB, substitute);3234 }3235 return N(gcdcAcB*ppH);3236 }3237 }3238 G_m= H;3239 newtonPoly= newtonPoly*(x  random_element);3240 m= m*(x  random_element);3241 if (!find (l, random_element))3242 l.append (random_element);3243 3244 if (getCharacteristic () > 3 && size (skeleton) < 100) // TODO consider a ratio of characteristic and degree?3245 {3246 CFArray Monoms;3247 CFArray *coeffMonoms= NULL;3248 do //second do3249 {3250 random_element= randomElement (m, V_buf, l, fail);3251 if (random_element == 0 && !fail)3252 {3253 if (!find (l, random_element))3254 l.append (random_element);3255 continue;3256 }3257 if (fail)3258 {3259 source= CFList();3260 dest= CFList();3261 int num_vars= tmin (getNumVars(A), getNumVars(B));3262 num_vars;3263 3264 choose_extension (d, num_vars, V_buf);3265 bool prim_fail= false;3266 Variable V_buf2;3267 prim_elem= primitiveElement (alpha, V_buf2, prim_fail);3268 3269 ASSERT (!prim_fail, "failure in integer factorizer");3270 if (prim_fail)3271 ; //ERROR3272 else3273 im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);3274 3275 DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));3276 DEBOUTLN (cerr, "getMipo (V_buf2)= " << getMipo (V_buf2));3277 inextension= true;3278 for (CFListIterator i= l; i.hasItem(); i++)3279 i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,3280 im_prim_elem, source, dest);3281 m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);3282 G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);3283 newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,3284 source, dest);3285 ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);3286 ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);3287 3288 gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,3289 source, dest);3290 3291 fail= false;3292 random_element= randomElement (m, V_buf, l, fail);3293 DEBOUTLN (cerr, "fail= " << fail);3294 CFList list;3295 TIMING_START (gcd_recursion);3296 3297 //sparseInterpolation3298 bool sparseFail= false;3299 if (LC (skeleton).inCoeffDomain())3300 G_random_element=3301 monicSparseInterpol (ppA (random_element, x), ppB (random_element, x),3302 skeleton,V_buf, sparseFail, coeffMonoms, Monoms);3303 else3304 G_random_element=3305 nonMonicSparseInterpol (ppA(random_element,x),ppB (random_element,x),3306 skeleton, V_buf, sparseFail, coeffMonoms,3307 Monoms);3308 if (sparseFail)3309 break;3310 3311 TIMING_END_AND_PRINT (gcd_recursion,3312 "time for recursive call: ");3313 DEBOUTLN (cerr, "G_random_element= " << G_random_element);3314 }3315 else3316 {3317 CFList list;3318 TIMING_START (gcd_recursion);3319 bool sparseFail= false;3320 if (LC (skeleton).inCoeffDomain())3321 G_random_element=3322 monicSparseInterpol (ppA (random_element, x), ppB (random_element, x),3323 skeleton,V_buf, sparseFail, coeffMonoms, Monoms);3324 else3325 G_random_element=3326 nonMonicSparseInterpol (ppA(random_element,x), ppB(random_element,x),3327 skeleton, V_buf, sparseFail, coeffMonoms,3328 Monoms);3329 if (sparseFail)3330 break;3331 3332 TIMING_END_AND_PRINT (gcd_recursion,3333 "time for recursive call: ");3334 DEBOUTLN (cerr, "G_random_element= " << G_random_element);3335 }3336 3337 d0= totaldegree (G_random_element, Variable(2),3338 Variable (G_random_element.level()));3339 if (d0 == 0)3340 {3341 if (substitute > 1)3342 return N(reverseSubst (gcdcAcB, substitute));3343 else3344 return N(gcdcAcB);3345 }3346 if (d0 > d)3347 {3348 if (!find (l, random_element))3349 l.append (random_element);3350 continue;3351 }3352 3353 G_random_element=3354 (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))3355 * G_random_element;3356 3357 d0= totaldegree (G_random_element, Variable(2),3358 Variable(G_random_element.level()));3359 if (d0 < d)3360 {3361 m= gcdlcAlcB;3362 newtonPoly= 1;3363 G_m= 0;3364 d= d0;3365 }3366 3367 TIMING_START (newton_interpolation);3368 H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);3369 TIMING_END_AND_PRINT (newton_interpolation,3370 "time for newton interpolation: ");3371 3372 //termination test3373 if (uni_lcoeff (H) == gcdlcAlcB)3374 {3375 cH= uni_content (H);3376 ppH= H/cH;3377 if (inextension)3378 {3379 CFList u, v;3380 //maybe it's better to test if ppH is an element of F(\alpha) before3381 //mapping down3382 DEBOUTLN (cerr, "ppH before mapDown= " << ppH);3383 ppH= mapDown (ppH, prim_elem, im_prim_elem, alpha, u, v);3384 ppH /= Lc(ppH);3385 DEBOUTLN (cerr, "ppH after mapDown= " << ppH);3386 if (fdivides (ppH, A) && fdivides (ppH, B))3387 {3388 if (substitute > 1)3389 {3390 ppH= reverseSubst (ppH, substitute);3391 gcdcAcB= reverseSubst (gcdcAcB, substitute);3392 }3393 return N(gcdcAcB*ppH);3394 }3395 }3396 else if (fdivides (ppH, A) && fdivides (ppH, B))3397 {3398 if (substitute > 1)3399 {3400 ppH= reverseSubst (ppH, substitute);3401 gcdcAcB= reverseSubst (gcdcAcB, substitute);3402 }3403 return N(gcdcAcB*ppH);3404 }3405 }3406 3407 G_m= H;3408 newtonPoly= newtonPoly*(x  random_element);3409 m= m*(x  random_element);3410 if (!find (l, random_element))3411 l.append (random_element);3412 3413 } while (1);3414 }3415 } while (1);3416 }3417 3418 CanonicalForm sparseGCDFp (const CanonicalForm& F, const CanonicalForm& G,3419 bool& topLevel, CFList& l)3420 {3421 CanonicalForm A= F;3422 CanonicalForm B= G;3423 if (F.isZero() && degree(G) > 0) return G/Lc(G);3424 else if (G.isZero() && degree (F) > 0) return F/Lc(F);3425 else if (F.isZero() && G.isZero()) return F.genOne();3426 if (F.inBaseDomain()  G.inBaseDomain()) return F.genOne();3427 if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);3428 if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);3429 if (F == G) return F/Lc(F);3430 3431 CFMap M,N;3432 int best_level= myCompress (A, B, M, N, topLevel);3433 3434 if (best_level == 0) return B.genOne();3435 3436 A= M(A);3437 B= M(B);3438 3439 Variable x= Variable (1);3440 3441 //univariate case3442 if (A.isUnivariate() && B.isUnivariate())3443 return N (gcd (A, B));3444 3445 int substitute= substituteCheck (A, B);3446 3447 if (substitute > 1)3448 subst (A, B, A, B, substitute);3449 3450 CanonicalForm cA, cB; // content of A and B3451 CanonicalForm ppA, ppB; // primitive part of A and B3452 CanonicalForm gcdcAcB;3453 if (topLevel)3454 {3455 if (best_level <= 2)3456 gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, best_level);3457 else3458 gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, 2);3459 }3460 else3461 {3462 cA = uni_content (A);3463 cB = uni_content (B);3464 gcdcAcB= gcd (cA, cB);3465 ppA= A/cA;3466 ppB= B/cB;3467 }3468 3469 CanonicalForm lcA, lcB; // leading coefficients of A and B3470 CanonicalForm gcdlcAlcB;3471 lcA= uni_lcoeff (ppA);3472 lcB= uni_lcoeff (ppB);3473 3474 if (fdivides (lcA, lcB))3475 {3476 if (fdivides (A, B))3477 return F/Lc(F);3478 }3479 if (fdivides (lcB, lcA))3480 {3481 if (fdivides (B, A))3482 return G/Lc(G);3483 }3484 3485 gcdlcAlcB= gcd (lcA, lcB);3486 3487 int d= totaldegree (ppA, Variable (2), Variable (ppA.level()));3488 int d0;3489 3490 if (d == 0)3491 {3492 if (substitute > 1)3493 return N(reverseSubst (gcdcAcB, substitute));3494 else3495 return N(gcdcAcB);3496 }3497 d0= totaldegree (ppB, Variable (2), Variable (ppB.level()));3498 3499 if (d0 < d)3500 d= d0;3501 3502 if (d == 0)3503 {3504 if (substitute > 1)3505 return N(reverseSubst (gcdcAcB, substitute));3506 else3507 return N(gcdcAcB);3508 }3509 3510 CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH, skeleton;3511 CanonicalForm newtonPoly= 1;3512 m= gcdlcAlcB;3513 G_m= 0;3514 H= 0;3515 bool fail= false;3516 bool topLevel2= topLevel;3517 int loops= 0;3518 topLevel= false;3519 bool inextension= false;3520 bool inextensionextension= false;3521 Variable V_buf, alpha;3522 CanonicalForm prim_elem, im_prim_elem;3523 CFList source, dest;3524 do //first do3525 {3526 if (inextension)3527 random_element= randomElement (m, alpha, l, fail);3528 else3529 random_element= FpRandomElement (m, l, fail);3530 if (random_element == 0 && !fail)3531 {3532 if (!find (l, random_element))3533 l.append (random_element);3534 continue;3535 }3536 3537 if (!fail && !inextension)3538 {3539 CFList list;3540 TIMING_START (gcd_recursion);3541 G_random_element=3542 sparseGCDFp (ppA (random_element,x), ppB (random_element,x), topLevel,3543 list);3544 TIMING_END_AND_PRINT (gcd_recursion,3545 "time for recursive call: ");3546 DEBOUTLN (cerr, "G_random_element= " << G_random_element);3547 }3548 else if (!fail && inextension)3549 {3550 CFList list;3551 TIMING_START (gcd_recursion);3552 G_random_element=3553 sparseGCDFq (ppA (random_element, x), ppB (random_element, x), alpha,3554 list, topLevel);3555 TIMING_END_AND_PRINT (gcd_recursion,3556 "time for recursive call: ");3557 DEBOUTLN (cerr, "G_random_element= " << G_random_element);3558 }3559 else if (fail && !inextension)3560 {3561 source= CFList();3562 dest= CFList();3563 CFList list;3564 CanonicalForm mipo;3565 int deg= 3;3566 do3567 {3568 mipo= randomIrredpoly (deg, x);3569 alpha= rootOf (mipo);3570 inextension= true;3571 fail= false;3572 random_element= randomElement (m, alpha, l, fail);3573 deg++;3574 } while (fail);3575 list= CFList();3576 TIMING_START (gcd_recursion);3577 G_random_element=3578 sparseGCDFq (ppA (random_element, x), ppB (random_element, x), alpha,3579 list, topLevel);3580 TIMING_END_AND_PRINT (gcd_recursion,3581 "time for recursive call: ");3582 DEBOUTLN (cerr, "G_random_element= " << G_random_element);3583 }3584 else if (fail && inextension)3585 {3586 source= CFList();3587 dest= CFList();3588 int num_vars= tmin (getNumVars(A), getNumVars(B));3589 num_vars;3590 V_buf= alpha;3591 choose_extension (d, num_vars, V_buf);3592 bool prim_fail= false;3593 Variable V_buf2;3594 CanonicalForm prim_elem, im_prim_elem;3595 prim_elem= primitiveElement (alpha, V_buf2, prim_fail);3596 3597 ASSERT (!prim_fail, "failure in integer factorizer");3598 if (prim_fail)3599 ; //ERROR3600 else3601 im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);3602 3603 DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));3604 DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));3605 3606 inextensionextension= true;3607 for (CFListIterator i= l; i.hasItem(); i++)3608 i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,3609 im_prim_elem, source, dest);3610 m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);3611 G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);3612 newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,3613 source, dest);3614 ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);3615 ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);3616 gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,3617 source, dest);3618 fail= false;3619 random_element= randomElement (m, V_buf, l, fail );3620 DEBOUTLN (cerr, "fail= " << fail);3621 CFList list;3622 TIMING_START (gcd_recursion);3623 G_random_element=3624 sparseGCDFq (ppA (random_element, x), ppB (random_element, x), V_buf,3625 list, topLevel);3626 TIMING_END_AND_PRINT (gcd_recursion,3627 "time for recursive call: ");3628 DEBOUTLN (cerr, "G_random_element= " << G_random_element);3629 alpha= V_buf;3630 }3631 3632 d0= totaldegree (G_random_element, Variable(2),3633 Variable (G_random_element.level()));3634 if (d0 == 0)3635 {3636 if (substitute > 1)3637 return N(reverseSubst (gcdcAcB, substitute));3638 else3639 return N(gcdcAcB);3640 }3641 if (d0 > d)3642 {3643 if (!find (l, random_element))3644 l.append (random_element);3645 continue;3646 }3647 3648 G_random_element=3649 (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))3650 * G_random_element;3651 3652 skeleton= G_random_element;3653 3654 if (size (skeleton) < 3 && getCharacteristic() <= 3)3655 {3656 return N (gcdcAcB*gcd_poly_p (ppA, ppB));3657 }3658 3659 d0= totaldegree (G_random_element, Variable(2),3660 Variable(G_random_element.level()));3661 if (d0 < d)3662 {3663 m= gcdlcAlcB;3664 newtonPoly= 1;3665 G_m= 0;3666 d= d0;3667 }3668 3669 H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);3670 3671 if (uni_lcoeff (H) == gcdlcAlcB)3672 {3673 cH= uni_content (H);3674 ppH= H/cH;3675 ppH /= Lc (ppH);3676 DEBOUTLN (cerr, "ppH= " << ppH);3677 3678 if (fdivides (ppH, A) && fdivides (ppH, B))3679 {3680 if (substitute > 1)3681 {3682 ppH= reverseSubst (ppH, substitute);3683 gcdcAcB= reverseSubst (gcdcAcB, substitute);3684 }3685 return N(gcdcAcB*ppH);3686 }3687 }3688 G_m= H;3689 newtonPoly= newtonPoly*(x  random_element);3690 m= m*(x  random_element);3691 if (!find (l, random_element))3692 l.append (random_element);3693 3694 if (getCharacteristic () > 3 && size (skeleton) < 100) // TODO consider a ratio of characteristic and degree?3695 {3696 CFArray Monoms;3697 CFArray* coeffMonoms= NULL;3698 3699 do //second do3700 {3701 if (inextension)3702 random_element= randomElement (m, alpha, l, fail);3703 else3704 random_element= FpRandomElement (m, l, fail);3705 if (random_element == 0 && !fail)3706 {3707 if (!find (l, random_element))3708 l.append (random_element);3709 continue;3710 }3711 3712 bool sparseFail= false;3713 if (!fail && !inextension)3714 {3715 CFList list;3716 TIMING_START (gcd_recursion);3717 if (LC (skeleton).inCoeffDomain())3718 G_random_element=3719 monicSparseInterpol (ppA (random_element, x), ppB (random_element, x),3720 skeleton, Variable(1), sparseFail, coeffMonoms,3721 Monoms);3722 else3723 G_random_element=3724 nonMonicSparseInterpol (ppA(random_element,x), ppB(random_element,x),3725 skeleton, Variable (1), sparseFail,3726 coeffMonoms, Monoms);3727 TIMING_END_AND_PRINT (gcd_recursion,3728 "time for recursive call: ");3729 DEBOUTLN (cerr, "G_random_element= " << G_random_element);3730 }3731 else if (!fail && inextension)3732 {3733 CFList list;3734 TIMING_START (gcd_recursion);3735 if (LC (skeleton).inCoeffDomain())3736 G_random_element=3737 monicSparseInterpol (ppA (random_element, x), ppB (random_element, x),3738 skeleton, alpha, sparseFail, coeffMonoms,3739 Monoms);3740 else3741 G_random_element=3742 nonMonicSparseInterpol (ppA(random_element,x), ppB(random_element,x),3743 skeleton, alpha, sparseFail, coeffMonoms,3744 Monoms);3745 TIMING_END_AND_PRINT (gcd_recursion,3746 "time for recursive call: ");3747 DEBOUTLN (cerr, "G_random_element= " << G_random_element);3748 }3749 else if (fail && !inextension)3750 {3751 source= CFList();3752 dest= CFList();3753 CFList list;3754 CanonicalForm mipo;3755 int deg= 3;3756 do3757 {3758 mipo= randomIrredpoly (deg, x);3759 alpha= rootOf (mipo);3760 inextension= true;3761 fail= false;3762 random_element= randomElement (m, alpha, l, fail);3763 deg++;3764 } while (fail);3765 list= CFList();3766 TIMING_START (gcd_recursion);3767 if (LC (skeleton).inCoeffDomain())3768 G_random_element=3769 monicSparseInterpol (ppA (random_element, x), ppB (random_element, x),3770 skeleton, alpha, sparseFail, coeffMonoms,3771 Monoms);3772 else3773 G_random_element=3774 nonMonicSparseInterpol (ppA(random_element,x), ppB(random_element,x),3775 skeleton, alpha, sparseFail, coeffMonoms,3776 Monoms);3777 TIMING_END_AND_PRINT (gcd_recursion,3778 "time for recursive call: ");3779 DEBOUTLN (cerr, "G_random_element= " << G_random_element);3780 }3781 else if (fail && inextension)3782 {3783 source= CFList();3784 dest= CFList();3785 int num_vars= tmin (getNumVars(A), getNumVars(B));3786 num_vars;3787 V_buf= alpha;3788 choose_extension (d, num_vars, V_buf);3789 bool prim_fail= false;3790 Variable V_buf2;3791 CanonicalForm prim_elem, im_prim_elem;3792 prim_elem= primitiveElement (alpha, V_buf2, prim_fail);3793 3794 ASSERT (!prim_fail, "failure in integer factorizer");3795 if (prim_fail)3796 ; //ERROR3797 else3798 im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);3799 3800 DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));3801 DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));3802 3803 inextensionextension= true;3804 for (CFListIterator i= l; i.hasItem(); i++)3805 i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,3806 im_prim_elem, source, dest);3807 m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);3808 G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);3809 newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,3810 source, dest);3811 ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);3812 ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);3813 gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,3814 source, dest);3815 fail= false;3816 random_element= randomElement (m, V_buf, l, fail );3817 DEBOUTLN (cerr, "fail= " << fail);3818 CFList list;3819 TIMING_START (gcd_recursion);3820 if (LC (skeleton).inCoeffDomain())3821 G_random_element=3822 monicSparseInterpol (ppA (random_element, x), ppB (random_element, x),3823 skeleton, V_buf, sparseFail, coeffMonoms,3824 Monoms);3825 else3826 G_random_element=3827 nonMonicSparseInterpol (ppA(random_element,x), ppB(random_element,x),3828 skeleton, V_buf, sparseFail,3829 coeffMonoms, Monoms);3830 TIMING_END_AND_PRINT (gcd_recursion,3831 "time for recursive call: ");3832 DEBOUTLN (cerr, "G_random_element= " << G_random_element);3833 alpha= V_buf;3834 }3835 3836 if (sparseFail)3837 break;3838 3839 d0= totaldegree (G_random_element, Variable(2),3840 Variable (G_random_element.level()));3841 if (d0 == 0)3842 {3843 if (substitute > 1)3844 return N(reverseSubst (gcdcAcB, substitute));3845 else3846 return N(gcdcAcB);3847 }3848 if (d0 > d)3849 {3850 if (!find (l, random_element))3851 l.append (random_element);3852 continue;3853 }3854 3855 G_random_element=3856 (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))3857 * G_random_element;3858 3859 d0= totaldegree (G_random_element, Variable(2),3860 Variable(G_random_element.level()));3861 if (d0 < d)3862 {3863 m= gcdlcAlcB;3864 newtonPoly= 1;3865 G_m= 0;3866 d= d0;3867 }3868 3869 TIMING_START (newton_interpolation);3870 H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);3871 TIMING_END_AND_PRINT (newton_interpolation,3872 "time for newton interpolation: ");3873 3874 //termination test3875 if (uni_lcoeff (H) == gcdlcAlcB)3876 {3877 cH= uni_content (H);3878 ppH= H/cH;3879 ppH /= Lc (ppH);3880 DEBOUTLN (cerr, "ppH= " << ppH);3881 if (fdivides (ppH, A) && fdivides (ppH, B))3882 {3883 if (substitute > 1)3884 {3885 ppH= reverseSubst (ppH, substitute);3886 gcdcAcB= reverseSubst (gcdcAcB, substitute);3887 }3888 return N(gcdcAcB*ppH);3889 }3890 }3891 3892 G_m= H;3893 newtonPoly= newtonPoly*(x  random_element);3894 m= m*(x  random_element);3895 if (!find (l, random_element))3896 l.append (random_element);3897 3898 } while (1); //end of second do3899 }3900 } while (1); //end of first do3901 }3902 3903 static inline3904 int compress4EZGCD (const CanonicalForm& F, const CanonicalForm& G, CFMap & M,3905 CFMap & N)3906 {3907 int n= tmax (F.level(), G.level());3908 int * degsf= new int [n + 1];3909 int * degsg= new int [n + 1];3910 3911 for (int i = 0; i <= n; i++)3912 degsf[i]= degsg[i]= 0;3913 3914 degsf= degrees (F, degsf);3915 degsg= degrees (G, degsg);3916 3917 int both_non_zero= 0;3918 int f_zero= 0;3919 int g_zero= 0;3920 int both_zero= 0;3921 3922 for (int i= 1; i <= n; i++)3923 {3924 if (degsf[i] != 0 && degsg[i] != 0)3925 {3926 both_non_zero++;3927 continue;3928 }3929 if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level())3930 {3931 f_zero++;3932 continue;3933 }3934 if (degsg[i] == 0 && degsf[i] && i <= F.level())3935 {3936 g_zero++;3937 continue;3938 }3939 }3940 3941 if (both_non_zero == 0) return 0;3942 3943 // map Variables which do not occur in both polynomials to higher levels3944 int k= 1;3945 int l= 1;3946 for (int i= 1; i <= n; i++)3947 {3948 if (degsf[i] != 0 && degsg[i] == 0 && i <= F.level())3949 {3950 if (k + both_non_zero != i)3951 {3952 M.newpair (Variable (i), Variable (k + both_non_zero));3953 N.newpair (Variable (k + both_non_zero), Variable (i));3954 }3955 k++;3956 }3957 if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level())3958 {3959 if (l + g_zero + both_non_zero != i)3960 {3961 M.newpair (Variable (i), Variable (l + g_zero + both_non_zero));3962 N.newpair (Variable (l + g_zero + both_non_zero), Variable (i));3963 }3964 l++;3965 }3966 }3967 3968 // sort Variables x_{i} in decreasing order of3969 // min(deg_{x_{i}}(f),deg_{x_{i}}(g))3970 int m= tmin (F.level(), G.level());3971 int max_min_deg;3972 k= both_non_zero;3973 l= 0;3974 int i= 1;3975 while (k > 0)3976 {3977 max_min_deg= tmin (degsf[i], degsg[i]);3978 while (max_min_deg == 0)3979 {3980 i++;3981 max_min_deg= tmin (degsf[i], degsg[i]);3982 }3983 for (int j= i + 1; j <= m; j++)3984 {3985 if (tmin (degsf[j],degsg[j]) < max_min_deg)3986 {3987 max_min_deg= tmin (degsf[j], degsg[j]);3988 l= j;3989 }3990 }3991 if (l != 0)3992 {3993 if (l != k)3994 {3995 M.newpair (Variable (l), Variable(k));3996 N.newpair (Variable (k), Variable(l));3997 degsf[l]= 0;3998 degsg[l]= 0;3999 l= 0;4000 }4001 else4002 {4003 degsf[l]= 0;4004 degsg[l]= 0;4005 l= 0;4006 }4007 }4008 else if (l == 0)4009 {4010 if (i != k)4011 {4012 M.newpair (Variable (i), Variable (k));4013 N.newpair (Variable (k), Variable (i));4014 degsf[i]= 0;4015 degsg[i]= 0;4016 }4017 else4018 {4019 degsf[i]= 0;4020 degsg[i]= 0;4021 }4022 i++;4023 }4024 k;4025 }4026 4027 delete [] degsf;4028 delete [] degsg;4029 4030 return both_non_zero;4031 }4032 4033 static4034 CanonicalForm myShift2Zero (const CanonicalForm& F, CFList& Feval,4035 const CFList& evaluation)4036 {4037 CanonicalForm A= F;4038 int k= 2;4039 for (CFListIterator i= evaluation; i.hasItem(); i++, k++)4040 A= A (Variable (k) + i.getItem(), k);4041 4042 CanonicalForm buf= A;4043 Feval= CFList();4044 Feval.append (buf);4045 for (k= evaluation.length() + 1; k > 2; k)4046 {4047 buf= mod (buf, Variable (k));4048 Feval.insert (buf);4049 }4050 return A;4051 }4052 4053 static CanonicalForm myReverseShift (const CanonicalForm& F,4054 const CFList& evaluation)4055 {4056 int l= evaluation.length() + 1;4057 CanonicalForm result= F;4058 CFListIterator j= evaluation;4059 for (int i= 2; i < l + 1; i++, j++)4060 {4061 if (F.level() < i)4062 continue;4063 result= result (Variable (i)  j.getItem(), i);4064 }4065 return result;4066 }4067 4068 static4069 bool Hensel_P (const CanonicalForm & U, CFArray & G, const Evaluation & A,4070 const Variable & x, const CFArray& LCs )4071 {4072 CFList factors;4073 factors.append (G[1]);4074 factors.append (G[2]);4075 CFList evaluation;4076 for (int i= A.min(); i <= A.max(); i++)4077 evaluation.append (A [i]);4078 CFList UEval;4079 CanonicalForm shiftedU= myShift2Zero (U, UEval, evaluation);4080 CFArray shiftedLCs= CFArray (2);4081 CFList shiftedLCsEval1, shiftedLCsEval2;4082 shiftedLCs[0]= myShift2Zero (LCs[1], shiftedLCsEval1, evaluation);4083 shiftedLCs[1]= myShift2Zero (LCs[2], shiftedLCsEval2, evaluation);4084 CanonicalForm LcUEval= LC (UEval.getFirst(), x);4085 factors.insert (1);4086 int liftBound= degree (UEval.getFirst(), 2) + 1;4087 CFArray Pi;4088 CFMatrix M= CFMatrix (liftBound, factors.length()  1);4089 CFList diophant;4090 CFArray lcs= CFArray (2);4091 lcs [0]= shiftedLCsEval1.getFirst();4092 lcs [1]= shiftedLCsEval2.getFirst();4093 henselLift122 (UEval.getFirst(), factors, liftBound, Pi, diophant, M,4094 lcs, false);4095 4096 bool noChange= true;4097 for (CFListIterator i= factors; i.hasItem(); i++)4098 {4099 if (degree (i.getItem(), 2) != 0)4100 noChange= false;4101 }4102 if (noChange)4103 return !noChange;4104 int * liftBounds;4105 if (U.level() > 2)4106 {4107 liftBounds= new int [U.level()  1];4108 liftBounds[0]= liftBound;4109 for (int i= 1; i < U.level()  1; i++)4110 liftBounds[i]= degree (shiftedU, Variable (i + 2)) + 1;4111 factors= henselLift2 (UEval, factors, liftBounds, U.level()  1, false,4112 shiftedLCsEval1, shiftedLCsEval2, Pi, diophant);4113 }4114 G[1]= factors.getFirst();4115 G[2]= factors.getLast();4116 G[1]= myReverseShift (G[1], evaluation);4117 G[2]= myReverseShift (G[2], evaluation);4118 return true;4119 }4120 4121 static4122 bool findeval_P (const CanonicalForm & F, const CanonicalForm & G,4123 CanonicalForm & Fb, CanonicalForm & Gb, CanonicalForm & Db,4124 FFREvaluation & b, int delta, int degF, int degG, int maxeval,4125 int & count )4126 {4127 if( count == 0 && delta != 0)4128 {4129 if( count++ > maxeval )4130 return false;4131 }4132 if (count > 0)4133 {4134 b.nextpoint();4135 if (count++ > maxeval)4136 return false;4137 }4138 while( true )4139 {4140 Fb = b( F );4141 if( degree( Fb, 1 ) == degF )4142 {4143 Gb = b( G );4144 if( degree( Gb, 1 ) == degG )4145 {4146 Db = gcd( Fb, Gb );4147 if( delta > 0 )4148 {4149 if( degree( Db, 1 ) <= delta )4150 return true;4151 }4152 else4153 return true;4154 }4155 }4156 b.nextpoint();4157 if( count++ > maxeval )4158 return false;4159 }4160 }4161 4162 /// Extended Zassenhaus GCD for finite fields4163 CanonicalForm EZGCD_P( const CanonicalForm & FF, const CanonicalForm & GG )4164 {4165 if (FF.isZero() && degree(GG) > 0) return GG/Lc(GG);4166 else if (GG.isZero() && degree (FF) > 0) return FF/Lc(FF);4167 else if (FF.isZero() && GG.isZero()) return FF.genOne();4168 if (FF.inBaseDomain()  GG.inBaseDomain()) return FF.genOne();4169 if (FF.isUnivariate() && fdivides(FF, GG)) return FF/Lc(FF);4170 if (GG.isUnivariate() && fdivides(GG, FF)) return GG/Lc(GG);4171 if (FF == GG) return FF/Lc(FF);4172 4173 CanonicalForm F, G, f, g, d, Fb, Gb, Db, Fbt, Gbt, Dbt, B0, B, D0, lcF, lcG,4174 lcD;4175 CFArray DD( 1, 2 ), lcDD( 1, 2 );4176 int degF, degG, delta, count;4177 int maxeval;4178 // bound on the number of eval. to use4179 maxeval = (getCharacteristic()/ilog2 (getCharacteristic()))*2;4180 count = 0; // number of eval. used4181 FFREvaluation b, bt;4182 bool gcdfound = false;4183 Variable x = Variable(1);4184 4185 F= FF;4186 G= GG;4187 4188 CFMap M,N;4189 int best_level= compress4EZGCD (F, G, M, N);4190 4191 if (best_level == 0) return G.genOne();4192 4193 F= M (F);4194 G= M (G);4195 4196 f = content( F, x ); g = content( G, x ); d = gcd( f, g );4197 F /= f; G /= g;4198 4199 if( F.isUnivariate() && G.isUnivariate() )4200 {4201 if( F.mvar() == G.mvar() )4202 d *= gcd( F, G );4203 return N (d);4204 }4205 4206 int maxNumVars= tmax (getNumVars (F), getNumVars (G));4207 Variable a,bv;4208 if (size (F)/maxNumVars > 500 && size (G)/maxNumVars > 500) //TODO find a4209 //better bound?4210 {4211 if (hasFirstAlgVar (F, a)  hasFirstAlgVar (G, a))4212 return N (d*GCD_Fp_extension (F, G, a));4213 else if (CFFactory::gettype() == GaloisFieldDomain)4214 return N (d*GCD_GF (F, G));4215 else4216 return N (d*GCD_small_p (F, G));4217 }4218 4219 if( gcd_test_one( F, G, false ) )4220 {4221 return N (d);4222 }4223 4224 lcF = LC( F, x ); lcG = LC( G, x );4225 lcD = gcd( lcF, lcG );4226 4227 if (fdivides (lcF, lcG)) {4228 if (fdivides (F, G))4229 return FF/Lc(FF);4230 }4231 if (fdivides (lcG, lcF))4232 {4233 if (fdivides (G, F))4234 return GG/Lc(GG);4235 }4236 4237 delta = 0;4238 degF = degree( F, x ); degG = degree( G, x );4239 if(hasFirstAlgVar(F,a))4240 {4241 if(hasFirstAlgVar(G,bv))4242 {4243 if(bv.level() > a.level())4244 a = bv;4245 }4246 b = FFREvaluation( 2, tmax(F.level(), G.level()), AlgExtRandomF( a ) );4247 }4248 else // F not in extension4249 {4250 if(hasFirstAlgVar(G,a))4251 b = FFREvaluation( 2, tmax(F.level(), G.level()), AlgExtRandomF( a ) );4252 else4253 { // both not in extension given by algebraic variable4254 if (CFFactory::gettype() != GaloisFieldDomain)4255 b = FFREvaluation( 2, tmax(F.level(), G.level()), FFRandom() );4256 else4257 b = FFREvaluation( 2, tmax(F.level(), G.level()), GFRandom() );4258 }4259 }4260 CanonicalForm cand;4261 while( !gcdfound )4262 {4263 if( !findeval_P( F, G, Fb, Gb, Db, b, delta, degF, degG, maxeval, count ))4264 { // too many eval. used > try another method4265 if (hasFirstAlgVar (F, a)  hasFirstAlgVar (G, a))4266 return N (d*sparseGCDFq (F, G, a));4267 if (CFFactory::gettype() == GaloisFieldDomain)4268 return N (d*GCD_GF (F, G));4269 else4270 return N (d*sparseGCDFp (F,G));4271 }4272 delta = degree( Db );4273 if( delta == 0 )4274 return N (d);4275 while( true )4276 {4277 bt = b;4278 if( !findeval_P( F, G, Fbt, Gbt, Dbt, bt, delta, degF, degG, maxeval, count))4279 { // too many eval. used > try another method4280 if (hasFirstAlgVar (F, a)  hasFirstAlgVar (G, a))4281 return N (d*sparseGCDFq (F, G, a));4282 if (CFFactory::gettype() == GaloisFieldDomain)4283 return N (d*GCD_GF (F, G));4284 else4285 return N (d*sparseGCDFp (F,G));4286 }4287 int dd = degree( Dbt );4288 if( dd == 0 )4289 return N (d);4290 if( dd == delta )4291 break;4292 if( dd < delta )4293 {4294 delta = dd;4295 b = bt;4296 Db = Dbt; Fb = Fbt; Gb = Gbt;4297 }4298 }4299 if( degF <= degG && delta == degF && fdivides( F, G ) )4300 return N (d*F);4301 if( degG < degF && delta == degG && fdivides( G, F ) )4302 return N (d*G);4303 if( delta != degF && delta != degG )4304 {4305 bool B_is_F;4306 CanonicalForm xxx1, xxx2;4307 CanonicalForm buf;4308 DD[1] = Fb / Db;4309 buf= Gb/Db;4310 xxx1 = gcd( DD[1], Db );4311 xxx2 = gcd( buf, Db );4312 if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) &&4313 (size (F) <= size (G)))4314  (xxx1.inCoeffDomain() && !xxx2.inCoeffDomain()))4315 {4316 B = F;4317 DD[2] = Db;4318 lcDD[1] = lcF;4319 lcDD[2] = lcD;4320 B_is_F = true;4321 }4322 else if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) &&4323 (size (G) < size (F)))4324  (!xxx1.inCoeffDomain() && xxx2.inCoeffDomain()))4325 {4326 DD[1] = buf;4327 B = G;4328 DD[2] = Db;4329 lcDD[1] = lcG;4330 lcDD[2] = lcD;4331 B_is_F = false;4332 }4333 else // special case handling4334 {4335 if (hasFirstAlgVar (F, a)  hasFirstAlgVar (G, a))4336 return N (d*sparseGCDFq (F, G, a));4337 if (CFFactory::gettype() == GaloisFieldDomain)4338 return N (d*GCD_GF (F, G));4339 else4340 return N (d*sparseGCDFp (F,G));4341 }4342 DD[2] = DD[2] * ( b( lcDD[2] ) / lc( DD[2] ) );4343 DD[1] = DD[1] * ( b( lcDD[1] ) / lc( DD[1] ) );4344 4345 gcdfound= Hensel_P (B*lcD, DD, b, x, lcDD);4346 4347 if( gcdfound )4348 {4349 cand = DD[2] / content( DD[2], Variable(1) );4350 if( B_is_F )4351 gcdfound = fdivides( cand, G );4352 else4353 gcdfound = fdivides( cand, F );4354 }4355 }4356 delta= 0;4357 }4358 return N (d*cand);4359 }4360 4361 1063 #endif
Note: See TracChangeset
for help on using the changeset viewer.