Changeset dec1024 in git
 Timestamp:
 Sep 20, 2010, 5:17:36 PM (13 years ago)
 Branches:
 (u'spielwiese', '8d54773d6c9e2f1d2593a28bc68b7eeab54ed529')
 Children:
 a0b5e27a5fb475f020a02a82976ffe87d3c097ae
 Parents:
 3f46965e1ad849ccb6d39d9e57cca2ecbbf10d86
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

factory/cf_gcd_smallp.cc
r3f4696 rdec1024 50 50 static inline 51 51 int myCompress (const CanonicalForm& F, const CanonicalForm& G, CFMap & M, 52 CFMap & N, bool& top _level)52 CFMap & N, bool& topLevel) 53 53 { 54 54 int n= tmax (F.level(), G.level()); … … 67 67 int both_zero= 0; 68 68 69 if (top _level)69 if (topLevel) 70 70 { 71 71 for (int i= 1; i <= n; i++) … … 88 88 } 89 89 90 if (both_non_zero == 0) return 0; 90 if (both_non_zero == 0) 91 { 92 delete [] degsf; 93 delete [] degsg; 94 return 0; 95 } 91 96 92 97 // map Variables which do not occur in both polynomials to higher levels … … 201 206 } 202 207 208 int 209 substituteCheck (const CanonicalForm& F, const CanonicalForm& G) 210 { 211 if (F.inCoeffDomain()  G.inCoeffDomain()) 212 return 0; 213 Variable x= Variable (1); 214 if (degree (F, x) <= 1  degree (G, x) <= 1) 215 return 0; 216 CanonicalForm f= swapvar (F, F.mvar(), x); //TODO swapping seems to be pretty expensive 217 CanonicalForm g= swapvar (G, G.mvar(), x); 218 int sizef= 0; 219 int sizeg= 0; 220 for (CFIterator i= f; i.hasTerms(); i++, sizef++) 221 { 222 if (i.exp() == 1) 223 return 0; 224 } 225 for (CFIterator i= g; i.hasTerms(); i++, sizeg++) 226 { 227 if (i.exp() == 1) 228 return 0; 229 } 230 int * expf= new int [sizef]; 231 int * expg= new int [sizeg]; 232 int j= 0; 233 for (CFIterator i= f; i.hasTerms(); i++, j++) 234 { 235 expf [j]= i.exp(); 236 } 237 j= 0; 238 for (CFIterator i= g; i.hasTerms(); i++, j++) 239 { 240 expg [j]= i.exp(); 241 } 242 243 int indf= sizef  1; 244 int indg= sizeg  1; 245 if (expf[indf] == 0) 246 indf; 247 if (expg[indg] == 0) 248 indg; 249 250 if (expg[indg] != expf [indf]  (expg[indg] == 1 && expf[indf] == 1)) 251 { 252 delete [] expg; 253 delete [] expf; 254 return 0; 255 } 256 // expf[indg] == expf[indf] after here 257 int result= expg[indg]; 258 for (int i= indf  1; i >= 0; i) 259 { 260 if (expf [i]%result != 0) 261 { 262 delete [] expg; 263 delete [] expf; 264 return 0; 265 } 266 } 267 268 for (int i= indg  1; i >= 0; i) 269 { 270 if (expg [i]%result != 0) 271 { 272 delete [] expg; 273 delete [] expf; 274 return 0; 275 } 276 } 277 278 delete [] expg; 279 delete [] expf; 280 return result; 281 } 282 283 // substiution 284 void 285 subst (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& A, 286 CanonicalForm& B, const int d 287 ) 288 { 289 if (d == 1) 290 { 291 A= F; 292 B= G; 293 return; 294 } 295 296 CanonicalForm C= 0; 297 CanonicalForm D= 0; 298 Variable x= Variable (1); 299 CanonicalForm f= swapvar (F, x, F.mvar()); 300 CanonicalForm g= swapvar (G, x, G.mvar()); 301 for (CFIterator i= f; i.hasTerms(); i++) 302 C += i.coeff()*power (f.mvar(), i.exp()/ d); 303 for (CFIterator i= g; i.hasTerms(); i++) 304 D += i.coeff()*power (g.mvar(), i.exp()/ d); 305 A= swapvar (C, x, F.mvar()); 306 B= swapvar (D, x, G.mvar()); 307 } 308 309 CanonicalForm 310 reverseSubst (const CanonicalForm& F, const int d) 311 { 312 if (d == 1) 313 return F; 314 315 Variable x= Variable (1); 316 if (degree (F, x) <= 0) 317 return F; 318 CanonicalForm f= swapvar (F, x, F.mvar()); 319 CanonicalForm result= 0; 320 for (CFIterator i= f; i.hasTerms(); i++) 321 result += i.coeff()*power (f.mvar(), d*i.exp()); 322 return swapvar (result, x, F.mvar()); 323 } 324 325 static inline CanonicalForm 326 uni_content (const CanonicalForm & F); 327 328 CanonicalForm 329 uni_content (const CanonicalForm& F, const Variable& x) 330 { 331 if (F.inCoeffDomain()) 332 return F.genOne(); 333 if (F.level() == x.level() && F.isUnivariate()) 334 return F; 335 if (F.level() != x.level() && F.isUnivariate()) 336 return F.genOne(); 337 338 if (x.level() != 1) 339 { 340 CanonicalForm f= swapvar (F, x, Variable (1)); 341 CanonicalForm result= uni_content (f); 342 return swapvar (result, x, Variable (1)); 343 } 344 else 345 return uni_content (F); 346 } 347 203 348 /// compute the content of F, where F is considered as an element of 204 349 /// \f$ R[x_{1}][x_{2},\ldots ,x_{n}] \f$ … … 231 376 return c; 232 377 } 378 } 379 380 CanonicalForm 381 extractContents (const CanonicalForm& F, const CanonicalForm& G, 382 CanonicalForm& contentF, CanonicalForm& contentG, 383 CanonicalForm& ppF, CanonicalForm& ppG, const int d) 384 { 385 CanonicalForm uniContentF, uniContentG, gcdcFcG; 386 contentF= 1; 387 contentG= 1; 388 ppF= F; 389 ppG= G; 390 CanonicalForm result= 1; 391 for (int i= 1; i <= d; i++) 392 { 393 uniContentF= uni_content (F, Variable (i)); 394 uniContentG= uni_content (G, Variable (i)); 395 gcdcFcG= gcd (uniContentF, uniContentG); 396 contentF *= uniContentF; 397 contentG *= uniContentG; 398 ppF /= uniContentF; 399 ppG /= uniContentG; 400 result *= gcdcFcG; 401 } 402 return result; 233 403 } 234 404 … … 336 506 337 507 /// GCD of F and G over \f$ F_{p}(\alpha ) \f$ , 338 /// l and top _level are only used internally, output is monic508 /// l and topLevel are only used internally, output is monic 339 509 /// based on Alg. 7.2. as described in "Algorithms for 340 510 /// Computer Algebra" by Geddes, Czapor, Labahn 341 511 CanonicalForm 342 512 GCD_Fp_extension (const CanonicalForm& F, const CanonicalForm& G, 343 Variable & alpha, CFList& l, bool& top _level)513 Variable & alpha, CFList& l, bool& topLevel) 344 514 { 345 515 CanonicalForm A= F; … … 354 524 355 525 CFMap M,N; 356 int best_level= myCompress (A, B, M, N, top _level);526 int best_level= myCompress (A, B, M, N, topLevel); 357 527 358 528 if (best_level == 0) return B.genOne(); … … 367 537 return N (gcd(A,B)); 368 538 539 int substitute= substituteCheck (A, B); 540 541 if (substitute > 1) 542 subst (A, B, A, B, substitute); 543 369 544 CanonicalForm cA, cB; // content of A and B 370 545 CanonicalForm ppA, ppB; // primitive part of A and B 371 546 CanonicalForm gcdcAcB; 372 547 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; 548 if (topLevel) 549 { 550 if (best_level <= 2) 551 gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, best_level); 552 else 553 gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, 2); 554 } 555 else 556 { 557 cA = uni_content (A); 558 cB = uni_content (B); 559 gcdcAcB= gcd (cA, cB); 560 ppA= A/cA; 561 ppB= B/cB; 562 } 380 563 381 564 CanonicalForm lcA, lcB; // leading coefficients of A and B … … 385 568 lcB= uni_lcoeff (ppB); 386 569 387 if (fdivides (lcA, lcB)) { 570 if (fdivides (lcA, lcB)) 571 { 388 572 if (fdivides (A, B)) 389 573 return F/Lc(F); … … 399 583 int d= totaldegree (ppA, Variable(2), Variable (ppA.level())); 400 584 401 if (d == 0) return N(gcdcAcB); 585 if (d == 0) 586 { 587 if (substitute > 1) 588 return N (reverseSubst (gcdcAcB, substitute)); 589 else 590 return N(gcdcAcB); 591 } 402 592 int d0= totaldegree (ppB, Variable(2), Variable (ppB.level())); 403 593 if (d0 < d) 404 594 d= d0; 405 if (d == 0) return N(gcdcAcB); 595 if (d == 0) 596 { 597 if (substitute > 1) 598 return N (reverseSubst (gcdcAcB, substitute)); 599 else 600 return N(gcdcAcB); 601 } 406 602 407 603 CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH; … … 413 609 H= 0; 414 610 bool fail= false; 415 top _level= false;611 topLevel= false; 416 612 bool inextension= false; 417 613 Variable V_buf= alpha; … … 461 657 G_random_element= 462 658 GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), V_buf, 463 list, top _level);659 list, topLevel); 464 660 TIMING_END_AND_PRINT (gcd_recursion, 465 661 "time for recursive call: "); … … 472 668 G_random_element= 473 669 GCD_Fp_extension (ppA(random_element, x), ppB(random_element, x), V_buf, 474 list, top _level);670 list, topLevel); 475 671 TIMING_END_AND_PRINT (gcd_recursion, 476 672 "time for recursive call: "); … … 480 676 d0= totaldegree (G_random_element, Variable(2), 481 677 Variable (G_random_element.level())); 482 if (d0 == 0) return N(gcdcAcB); 678 if (d0 == 0) 679 { 680 if (substitute > 1) 681 return N (reverseSubst (gcdcAcB, substitute)); 682 else 683 return N(gcdcAcB); 684 } 483 685 if (d0 > d) 484 686 { … … 522 724 DEBOUTLN (cerr, "ppH after mapDown= " << ppH); 523 725 if (fdivides (ppH, A) && fdivides (ppH, B)) 726 { 727 if (substitute > 1) 728 { 729 ppH= reverseSubst (ppH, substitute); 730 gcdcAcB= reverseSubst (gcdcAcB, substitute); 731 } 524 732 return N(gcdcAcB*ppH); 733 } 525 734 } 526 735 else if (fdivides (ppH, A) && fdivides (ppH, B)) 527 return N(gcdcAcB*(ppH/Lc(ppH))); 736 { 737 if (substitute > 1) 738 { 739 ppH= reverseSubst (ppH, substitute); 740 gcdcAcB= reverseSubst (gcdcAcB, substitute); 741 } 742 return N(gcdcAcB*ppH); 743 } 528 744 } 529 745 … … 580 796 /// the size of the base field is bounded by 2^16, output is monic 581 797 CanonicalForm GCD_GF (const CanonicalForm& F, const CanonicalForm& G, 582 CFList& l, bool& top _level)798 CFList& l, bool& topLevel) 583 799 { 584 800 CanonicalForm A= F; … … 593 809 594 810 CFMap M,N; 595 int best_level= myCompress (A, B, M, N, top _level);811 int best_level= myCompress (A, B, M, N, topLevel); 596 812 597 813 if (best_level == 0) return B.genOne(); … … 606 822 return N (gcd (A, B)); 607 823 824 int substitute= substituteCheck (A, B); 825 826 if (substitute > 1) 827 subst (A, B, A, B, substitute); 828 608 829 CanonicalForm cA, cB; // content of A and B 609 830 CanonicalForm ppA, ppB; // primitive part of A and B 610 831 CanonicalForm gcdcAcB; 611 832 612 cA = uni_content (A); 613 cB = uni_content (B); 614 615 gcdcAcB= gcd (cA, cB); 616 617 ppA= A/cA; 618 ppB= B/cB; 833 if (topLevel) 834 { 835 if (best_level <= 2) 836 gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, best_level); 837 else 838 gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, 2); 839 } 840 else 841 { 842 cA = uni_content (A); 843 cB = uni_content (B); 844 gcdcAcB= gcd (cA, cB); 845 ppA= A/cA; 846 ppB= B/cB; 847 } 619 848 620 849 CanonicalForm lcA, lcB; // leading coefficients of A and B … … 638 867 639 868 int d= totaldegree (ppA, Variable(2), Variable (ppA.level())); 640 if (d == 0) return N(gcdcAcB); 869 if (d == 0) 870 { 871 if (substitute > 1) 872 return N (reverseSubst (gcdcAcB, substitute)); 873 else 874 return N(gcdcAcB); 875 } 641 876 int d0= totaldegree (ppB, Variable(2), Variable (ppB.level())); 642 877 if (d0 < d) 643 878 d= d0; 644 if (d == 0) return N(gcdcAcB); 879 if (d == 0) 880 { 881 if (substitute > 1) 882 return N (reverseSubst (gcdcAcB, substitute)); 883 else 884 return N(gcdcAcB); 885 } 645 886 646 887 CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH; … … 652 893 H= 0; 653 894 bool fail= false; 654 top _level= false;895 topLevel= false; 655 896 bool inextension= false; 656 897 int p; … … 694 935 TIMING_START (gcd_recursion); 695 936 G_random_element= GCD_GF (ppA(random_element, x), ppB(random_element, x), 696 list, top _level);937 list, topLevel); 697 938 TIMING_END_AND_PRINT (gcd_recursion, 698 939 "time for recursive call: "); … … 704 945 TIMING_START (gcd_recursion); 705 946 G_random_element= GCD_GF (ppA(random_element, x), ppB(random_element, x), 706 list, top _level);947 list, topLevel); 707 948 TIMING_END_AND_PRINT (gcd_recursion, 708 949 "time for recursive call: "); … … 717 958 { 718 959 setCharacteristic (p, k, gf_name_buf); 719 return N(gcdcAcB); 960 { 961 if (substitute > 1) 962 return N (reverseSubst (gcdcAcB, substitute)); 963 else 964 return N(gcdcAcB); 965 } 720 966 } 721 967 else 722 return N(gcdcAcB); 968 { 969 if (substitute > 1) 970 return N (reverseSubst (gcdcAcB, substitute)); 971 else 972 return N(gcdcAcB); 973 } 723 974 } 724 975 if (d0 > d) … … 759 1010 ppH= GFMapDown (ppH, k); 760 1011 DEBOUTLN (cerr, "ppH after mapDown= " << ppH); 1012 if (substitute > 1) 1013 { 1014 ppH= reverseSubst (ppH, substitute); 1015 gcdcAcB= reverseSubst (gcdcAcB, substitute); 1016 } 761 1017 setCharacteristic (p, k, gf_name_buf); 762 1018 return N(gcdcAcB*ppH); … … 766 1022 { 767 1023 if (fdivides (ppH, A) && fdivides (ppH, B)) 1024 { 1025 if (substitute > 1) 1026 { 1027 ppH= reverseSubst (ppH, substitute); 1028 gcdcAcB= reverseSubst (gcdcAcB, substitute); 1029 } 768 1030 return N(gcdcAcB*ppH); 1031 } 769 1032 } 770 1033 } … … 829 1092 830 1093 CanonicalForm GCD_small_p (const CanonicalForm& F, const CanonicalForm& G, 831 bool& top _level, CFList& l)1094 bool& topLevel, CFList& l) 832 1095 { 833 1096 CanonicalForm A= F; … … 842 1105 843 1106 CFMap M,N; 844 int best_level= myCompress (A, B, M, N, top _level);1107 int best_level= myCompress (A, B, M, N, topLevel); 845 1108 846 1109 if (best_level == 0) return B.genOne(); … … 855 1118 return N (gcd (A, B)); 856 1119 1120 int substitute= substituteCheck (A, B); 1121 1122 if (substitute > 1) 1123 subst (A, B, A, B, substitute); 1124 857 1125 CanonicalForm cA, cB; // content of A and B 858 1126 CanonicalForm ppA, ppB; // primitive part of A and B 859 1127 CanonicalForm gcdcAcB; 860 cA = uni_content (A); 861 cB = uni_content (B); 862 gcdcAcB= gcd (cA, cB); 863 ppA= A/cA; 864 ppB= B/cB; 1128 1129 if (topLevel) 1130 { 1131 if (best_level <= 2) 1132 gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, best_level); 1133 else 1134 gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, 2); 1135 } 1136 else 1137 { 1138 cA = uni_content (A); 1139 cB = uni_content (B); 1140 gcdcAcB= gcd (cA, cB); 1141 ppA= A/cA; 1142 ppB= B/cB; 1143 } 865 1144 866 1145 CanonicalForm lcA, lcB; // leading coefficients of A and B … … 885 1164 int d0; 886 1165 887 if (d == 0) return N(gcdcAcB); 1166 if (d == 0) 1167 { 1168 if (substitute > 1) 1169 return N (reverseSubst (gcdcAcB, substitute)); 1170 else 1171 return N(gcdcAcB); 1172 } 888 1173 d0= totaldegree (ppB, Variable (2), Variable (ppB.level())); 889 1174 … … 891 1176 d= d0; 892 1177 893 if (d == 0) return N(gcdcAcB); 1178 if (d == 0) 1179 { 1180 if (substitute > 1) 1181 return N (reverseSubst (gcdcAcB, substitute)); 1182 else 1183 return N(gcdcAcB); 1184 } 894 1185 895 1186 CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH; … … 899 1190 G_m= 0; 900 1191 Variable alpha, V_buf; 901 int p, k;902 1192 bool fail= false; 903 1193 bool inextension= false; 904 1194 bool inextensionextension= false; 905 top_level= false; 906 CanonicalForm CF_buf; 1195 topLevel= false; 907 1196 CFList source, dest; 908 CanonicalForm gcdcheck;909 1197 do 910 1198 { … … 919 1207 TIMING_START (gcd_recursion); 920 1208 G_random_element= 921 GCD_small_p (ppA (random_element,x), ppB (random_element,x), top _level,1209 GCD_small_p (ppA (random_element,x), ppB (random_element,x), topLevel, 922 1210 list); 923 1211 TIMING_END_AND_PRINT (gcd_recursion, … … 931 1219 G_random_element= 932 1220 GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), alpha, 933 list, top _level);1221 list, topLevel); 934 1222 TIMING_END_AND_PRINT (gcd_recursion, 935 1223 "time for recursive call: "); … … 942 1230 CFList list; 943 1231 CanonicalForm mipo; 944 int deg= 3;1232 int deg= 2; 945 1233 do { 946 1234 mipo= randomIrredpoly (deg, x); … … 955 1243 G_random_element= 956 1244 GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), alpha, 957 list, top _level);1245 list, topLevel); 958 1246 TIMING_END_AND_PRINT (gcd_recursion, 959 1247 "time for recursive call: "); … … 1001 1289 G_random_element= 1002 1290 GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), V_buf, 1003 list, top _level);1291 list, topLevel); 1004 1292 TIMING_END_AND_PRINT (gcd_recursion, 1005 1293 "time for recursive call: "); … … 1013 1301 if (d0 == 0) 1014 1302 { 1015 return N(gcdcAcB); 1303 if (substitute > 1) 1304 return N (reverseSubst (gcdcAcB, substitute)); 1305 else 1306 return N(gcdcAcB); 1016 1307 } 1017 1308 if (d0 > d) … … 1050 1341 DEBOUTLN (cerr, "ppH= " << ppH); 1051 1342 if (fdivides (ppH, A) && fdivides (ppH, B)) 1343 { 1344 if (substitute > 1) 1345 { 1346 ppH= reverseSubst (ppH, substitute); 1347 gcdcAcB= reverseSubst (gcdcAcB, substitute); 1348 } 1052 1349 return N(gcdcAcB*ppH); 1350 } 1053 1351 } 1054 1352
Note: See TracChangeset
for help on using the changeset viewer.