Changeset 9ff686 in git
 Timestamp:
 Jun 14, 2011, 11:29:25 AM (12 years ago)
 Branches:
 (u'jengelhdatetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', 'a800fe4b3e9d37a38c5a10cc0ae9dfa0c15a4ee6')
 Children:
 d3b70ae172ff67aecbfd2f47077e21aa319554fd
 Parents:
 3ffd6d2e233ebf7267b3b68dc7d6d7d506305516
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

factory/cf_gcd_smallp.cc
r3ffd6d r9ff686 52 52 /// compressing two polynomials F and G, M is used for compressing, 53 53 /// N to reverse the compression 54 static inline55 54 int myCompress (const CanonicalForm& F, const CanonicalForm& G, CFMap & M, 56 55 CFMap & N, bool topLevel) … … 126 125 // sort Variables x_{i} in increasing order of 127 126 // min(deg_{x_{i}}(f),deg_{x_{i}}(g)) 128 int m= tm in(F.level(), G.level());129 int m ax_min_deg;127 int m= tmax (F.level(), G.level()); 128 int min_max_deg; 130 129 k= both_non_zero; 131 130 l= 0; … … 133 132 while (k > 0) 134 133 { 135 max_min_deg= tmin (degsf[i], degsg[i]); 136 while (max_min_deg == 0) 134 if (degsf [i] != 0 && degsg [i] != 0) 135 min_max_deg= tmax (degsf[i], degsg[i]); 136 else 137 min_max_deg= 0; 138 while (min_max_deg == 0) 137 139 { 138 140 i++; 139 max_min_deg= tmin (degsf[i], degsg[i]); 141 if (degsf [i] != 0 && degsg [i] != 0) 142 min_max_deg= tmax (degsf[i], degsg[i]); 143 else 144 min_max_deg= 0; 140 145 } 141 146 for (int j= i + 1; j <= m; j++) 142 147 { 143 if (tmin (degsf[j],degsg[j]) >= max_min_deg) 144 { 145 max_min_deg= tmin (degsf[j], degsg[j]); 148 if (degsf[j] != 0 && degsg [j] != 0 && 149 tmax (degsf[j],degsg[j]) <= min_max_deg) 150 { 151 min_max_deg= tmax (degsf[j], degsg[j]); 146 152 l= j; 147 153 } … … 209 215 return 1; 210 216 } 211 212 217 213 218 int … … 219 224 if (degree (F, x) <= 1  degree (G, x) <= 1) 220 225 return 0; 221 CanonicalForm f= swapvar (F, F.mvar(), x); //TODO swapping seems to bepretty expensive226 CanonicalForm f= swapvar (F, F.mvar(), x); //TODO swapping is pretty expensive 222 227 CanonicalForm g= swapvar (G, G.mvar(), x); 223 228 int sizef= 0; … … 437 442 /// newtonPoly must be equal to (x  alpha_1) * ... * (x  alpha_{n1}) 438 443 static inline CanonicalForm 439 newtonInterp(const CanonicalForm alpha, const CanonicalForm u, const CanonicalForm newtonPoly, const CanonicalForm oldInterPoly, const Variable & x) 444 newtonInterp(const CanonicalForm alpha, const CanonicalForm u, 445 const CanonicalForm newtonPoly, const CanonicalForm oldInterPoly, 446 const Variable & x) 440 447 { 441 448 CanonicalForm interPoly; 442 449 443 interPoly = oldInterPoly + ((u  oldInterPoly(alpha, x)) / newtonPoly(alpha, x)) * newtonPoly; 450 interPoly= oldInterPoly + ((u  oldInterPoly(alpha, x))/newtonPoly(alpha, x)) 451 *newtonPoly; 444 452 return interPoly; 445 453 } … … 487 495 } while (find (list, random)); 488 496 return random; 497 } 498 499 static inline 500 Variable chooseExtension (const Variable & alpha) 501 { 502 zz_p::init (getCharacteristic()); 503 zz_pX NTLIrredpoly; 504 int i, m; 505 // extension of F_p needed 506 if (alpha.level() == 1) 507 { 508 i= 1; 509 m= 2; 510 } //extension of F_p(alpha) 511 if (alpha.level() != 1) 512 { 513 i= 4; 514 m= degree (getMipo (alpha)); 515 } 516 BuildIrred (NTLIrredpoly, i*m); 517 CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1)); 518 return rootOf (newMipo); 489 519 } 490 520 … … 625 655 source= CFList(); 626 656 dest= CFList(); 627 int num_vars= tmin (getNumVars(A), getNumVars(B)); 628 num_vars; 629 630 choose_extension (d, num_vars, V_buf); 657 658 Variable V_buf3= V_buf; 659 V_buf= chooseExtension (V_buf); 631 660 bool prim_fail= false; 632 661 Variable V_buf2; 633 662 prim_elem= primitiveElement (alpha, V_buf2, prim_fail); 663 664 if (V_buf3 != alpha) 665 { 666 m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest); 667 G_m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest); 668 newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha, 669 source, dest); 670 ppA= mapDown (ppA, prim_elem, im_prim_elem, alpha, source, dest); 671 ppB= mapDown (ppB, prim_elem, im_prim_elem, alpha, source, dest); 672 gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha, 673 source, dest); 674 for (CFListIterator i= l; i.hasItem(); i++) 675 i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha, 676 source, dest); 677 } 634 678 635 679 ASSERT (!prim_fail, "failure in integer factorizer"); … … 909 953 if (fail) 910 954 { 911 int num_vars= tmin (getNumVars(A), getNumVars(B));912 num_vars;913 955 p= getCharacteristic(); 914 expon= (int) floor ((log ((double) ipower (d + 1, num_vars))/log ((double) p))); 915 if (expon < 2) 916 expon= 2; 956 expon= 2; 917 957 kk= getGFDegree(); 918 958 if (ipower (p, kk*expon) < (1 << 16)) … … 1202 1242 { 1203 1243 if (inextension) 1204 random_element= randomElement (m, alpha, l, fail);1244 random_element= randomElement (m, V_buf, l, fail); 1205 1245 else 1206 1246 random_element= FpRandomElement (m, l, fail); … … 1244 1284 } while (fail); 1245 1285 list= CFList(); 1286 V_buf= alpha; 1246 1287 TIMING_START (gcd_recursion); 1247 1288 G_random_element= … … 1256 1297 source= CFList(); 1257 1298 dest= CFList(); 1258 int num_vars= tmin (getNumVars(A), getNumVars(B)); 1259 num_vars; 1260 V_buf= alpha; 1261 choose_extension (d, num_vars, V_buf); 1299 1300 Variable V_buf3= V_buf; 1301 V_buf= chooseExtension (V_buf); 1262 1302 bool prim_fail= false; 1263 1303 Variable V_buf2; 1264 1304 CanonicalForm prim_elem, im_prim_elem; 1265 1305 prim_elem= primitiveElement (alpha, V_buf2, prim_fail); 1306 1307 if (V_buf3 != alpha) 1308 { 1309 m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest); 1310 G_m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest); 1311 newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha, 1312 source, dest); 1313 ppA= mapDown (ppA, prim_elem, im_prim_elem, alpha, source, dest); 1314 ppB= mapDown (ppB, prim_elem, im_prim_elem, alpha, source, dest); 1315 gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha, source, 1316 dest); 1317 for (CFListIterator i= l; i.hasItem(); i++) 1318 i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha, 1319 source, dest); 1320 } 1266 1321 1267 1322 ASSERT (!prim_fail, "failure in integer factorizer"); … … 1297 1352 "time for recursive call: "); 1298 1353 DEBOUTLN (cerr, "G_random_element= " << G_random_element); 1299 alpha= V_buf;1300 1354 } 1301 1355 … … 1983 2037 if (evalFail) 1984 2038 { 1985 if (V_buf != Variable (1))2039 if (V_buf.level() != 1) 1986 2040 { 1987 2041 do 1988 2042 { 1989 int num_vars= tmin (getNumVars(A), getNumVars(B)); 1990 int d= totaldegree (A, Variable(2), Variable (A.level())); 1991 d= tmin (d, totaldegree (B, Variable(2), Variable (B.level()))); 1992 Variable V_buf2= V_buf; 1993 choose_extension (d, num_vars, V_buf2); 2043 Variable V_buf2= chooseExtension (V_buf); 1994 2044 source= CFList(); 1995 2045 dest= CFList(); … … 2020 2070 } 2021 2071 2022 if (alpha != Variable (1))2072 if (alpha.level() != 1) 2023 2073 { 2024 2074 A= mapUp (A, V_buf, V_buf2, prim_elem, im_prim_elem, source,dest); … … 2134 2184 delete[] pL; 2135 2185 2136 if (alpha != Variable (1)&& V_buf != alpha)2186 if (alpha.level() != 1 && V_buf != alpha) 2137 2187 { 2138 2188 CFList u, v; … … 2257 2307 if (evalFail) 2258 2308 { 2259 if (V_buf != Variable (1))2309 if (V_buf.level() != 1) 2260 2310 { 2261 2311 do 2262 2312 { 2263 int num_vars= tmin (getNumVars(A), getNumVars(B)); 2264 int d= totaldegree (A, Variable(2), Variable (A.level())); 2265 d= tmin (d, totaldegree (B, Variable(2), Variable (B.level()))); 2266 Variable V_buf2= V_buf; 2267 choose_extension (d, num_vars, V_buf2); 2313 Variable V_buf2= chooseExtension (V_buf); 2268 2314 source= CFList(); 2269 2315 dest= CFList(); … … 2446 2492 pMat[1]= Mat; 2447 2493 2448 if (V_buf != x)2494 if (V_buf.level() != 1) 2449 2495 solution= solveSystemFq (pMat[1], 2450 2496 pL[minimalColumnsIndex], V_buf); … … 2577 2623 for (int i= 0; i < skelSize; i++) 2578 2624 { 2579 if (V_buf == x)2625 if (V_buf.level() == 1) 2580 2626 rk= gaussianElimFp (pMat[i], pL[i]); 2581 2627 else … … 2623 2669 } 2624 2670 2625 if (V_buf != x)2671 if (V_buf.level() != 1) 2626 2672 solution= solveSystemFq (Mat, bufArray, V_buf); 2627 2673 else … … 2650 2696 result += Monoms [ind]*bufSolution[i]; 2651 2697 } 2652 if (alpha != Variable (1)&& V_buf != alpha)2698 if (alpha.level() != 1 && V_buf != alpha) 2653 2699 { 2654 2700 CFList u, v; … … 2741 2787 delete[] pM; 2742 2788 2743 if (alpha != Variable (1)&& V_buf != alpha)2789 if (alpha.level() != 1 && V_buf != alpha) 2744 2790 { 2745 2791 CFList u, v; … … 2874 2920 source= CFList(); 2875 2921 dest= CFList(); 2876 int num_vars= tmin (getNumVars(A), getNumVars(B)); 2877 num_vars; 2878 2879 choose_extension (d, num_vars, V_buf); 2922 2923 Variable V_buf3= V_buf; 2924 V_buf= chooseExtension (V_buf); 2880 2925 bool prim_fail= false; 2881 2926 Variable V_buf2; 2882 2927 prim_elem= primitiveElement (alpha, V_buf2, prim_fail); 2928 2929 if (V_buf3 != alpha) 2930 { 2931 m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest); 2932 G_m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest); 2933 newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha, 2934 source, dest); 2935 ppA= mapDown (ppA, prim_elem, im_prim_elem, alpha, source, dest); 2936 ppB= mapDown (ppB, prim_elem, im_prim_elem, alpha, source, dest); 2937 gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha, source, 2938 dest); 2939 for (CFListIterator i= l; i.hasItem(); i++) 2940 i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha, 2941 source, dest); 2942 } 2883 2943 2884 2944 ASSERT (!prim_fail, "failure in integer factorizer"); … … 3015 3075 source= CFList(); 3016 3076 dest= CFList(); 3017 int num_vars= tmin (getNumVars(A), getNumVars(B)); 3018 num_vars; 3019 3020 choose_extension (d, num_vars, V_buf); 3077 3078 Variable V_buf3= V_buf; 3079 V_buf= chooseExtension (V_buf); 3021 3080 bool prim_fail= false; 3022 3081 Variable V_buf2; 3023 3082 prim_elem= primitiveElement (alpha, V_buf2, prim_fail); 3083 3084 if (V_buf3 != alpha) 3085 { 3086 m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest); 3087 G_m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest); 3088 newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha, 3089 source, dest); 3090 ppA= mapDown (ppA, prim_elem, im_prim_elem, alpha, source, dest); 3091 ppB= mapDown (ppB, prim_elem, im_prim_elem, alpha, source, dest); 3092 gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha, 3093 source, dest); 3094 for (CFListIterator i= l; i.hasItem(); i++) 3095 i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha, 3096 source, dest); 3097 } 3024 3098 3025 3099 ASSERT (!prim_fail, "failure in integer factorizer"); … … 3279 3353 { 3280 3354 if (inextension) 3281 random_element= randomElement (m, alpha, l, fail);3355 random_element= randomElement (m, V_buf, l, fail); 3282 3356 else 3283 3357 random_element= FpRandomElement (m, l, fail); … … 3327 3401 deg++; 3328 3402 } while (fail); 3403 V_buf= alpha; 3329 3404 list= CFList(); 3330 3405 TIMING_START (gcd_recursion); … … 3340 3415 source= CFList(); 3341 3416 dest= CFList(); 3342 int num_vars= tmin (getNumVars(A), getNumVars(B)); 3343 num_vars; 3344 V_buf= alpha; 3345 choose_extension (d, num_vars, V_buf); 3417 3418 Variable V_buf3= V_buf; 3419 V_buf= chooseExtension (V_buf); 3346 3420 bool prim_fail= false; 3347 3421 Variable V_buf2; 3348 3422 CanonicalForm prim_elem, im_prim_elem; 3349 3423 prim_elem= primitiveElement (alpha, V_buf2, prim_fail); 3424 3425 if (V_buf3 != alpha) 3426 { 3427 m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest); 3428 G_m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest); 3429 newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha, source, 3430 dest); 3431 ppA= mapDown (ppA, prim_elem, im_prim_elem, alpha, source, dest); 3432 ppB= mapDown (ppB, prim_elem, im_prim_elem, alpha, source, dest); 3433 gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha, source, 3434 dest); 3435 for (CFListIterator i= l; i.hasItem(); i++) 3436 i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha, 3437 source, dest); 3438 } 3350 3439 3351 3440 ASSERT (!prim_fail, "failure in integer factorizer"); … … 3381 3470 "time for recursive call: "); 3382 3471 DEBOUTLN (cerr, "G_random_element= " << G_random_element); 3383 alpha= V_buf;3384 3472 } 3385 3473 … … 3512 3600 deg++; 3513 3601 } while (fail); 3602 V_buf= alpha; 3514 3603 list= CFList(); 3515 3604 TIMING_START (gcd_recursion); … … 3532 3621 source= CFList(); 3533 3622 dest= CFList(); 3534 int num_vars= tmin (getNumVars(A), getNumVars(B)); 3535 num_vars; 3536 V_buf= alpha; 3537 choose_extension (d, num_vars, V_buf); 3623 3624 Variable V_buf3= V_buf; 3625 V_buf= chooseExtension (V_buf); 3538 3626 bool prim_fail= false; 3539 3627 Variable V_buf2; 3540 3628 CanonicalForm prim_elem, im_prim_elem; 3541 3629 prim_elem= primitiveElement (alpha, V_buf2, prim_fail); 3630 3631 if (V_buf3 != alpha) 3632 { 3633 m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest); 3634 G_m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest); 3635 newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha, 3636 source, dest); 3637 ppA= mapDown (ppA, prim_elem, im_prim_elem, alpha, source, dest); 3638 ppB= mapDown (ppB, prim_elem, im_prim_elem, alpha, source, dest); 3639 gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha, 3640 source, dest); 3641 for (CFListIterator i= l; i.hasItem(); i++) 3642 i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha, 3643 source, dest); 3644 } 3542 3645 3543 3646 ASSERT (!prim_fail, "failure in integer factorizer"); … … 3580 3683 "time for recursive call: "); 3581 3684 DEBOUTLN (cerr, "G_random_element= " << G_random_element); 3582 alpha= V_buf;3583 3685 } 3584 3686 … … 3822 3924 3823 3925 static inline 3824 bool Hensel_P (const CanonicalForm & U, CFArray & G, const Evaluation & A, 3825 const Variable & x, const CFArray& LCs ) 3926 Evaluation optimize4Lift (const CanonicalForm& F, CFMap & M, 3927 CFMap & N, const Evaluation& A) 3928 { 3929 int n= F.level(); 3930 int * degsf= new int [n + 1]; 3931 3932 for (int i = 0; i <= n; i++) 3933 degsf[i]= 0; 3934 3935 degsf= degrees (F, degsf); 3936 3937 Evaluation result= Evaluation (A.min(), A.max()); 3938 ASSERT (A.min() == 2, "expected A.min() == 2"); 3939 ASSERT (A.max() == n, "expected A.max() == n"); 3940 int max_deg; 3941 int k= n; 3942 int l= 1; 3943 int i= 2; 3944 int pos= 2; 3945 while (k > 1) 3946 { 3947 max_deg= degsf [i]; 3948 while (max_deg == 0) 3949 { 3950 i++; 3951 max_deg= degsf [i]; 3952 } 3953 l= i; 3954 for (int j= i + 1; j <= n; j++) 3955 { 3956 if (degsf[j] > max_deg) 3957 { 3958 max_deg= degsf[j]; 3959 l= j; 3960 } 3961 } 3962 3963 if (l <= n) 3964 { 3965 if (l != pos) 3966 { 3967 result.setValue (pos, A [l]); 3968 M.newpair (Variable (l), Variable (pos)); 3969 N.newpair (Variable (pos), Variable (l)); 3970 degsf[l]= 0; 3971 l= 2; 3972 if (k == 2 && n == 3) 3973 { 3974 result.setValue (l, A [pos]); 3975 M.newpair (Variable (pos), Variable (l)); 3976 N.newpair (Variable (l), Variable (pos)); 3977 degsf[pos]= 0; 3978 } 3979 } 3980 else 3981 { 3982 result.setValue (l, A [l]); 3983 degsf [l]= 0; 3984 } 3985 } 3986 pos++; 3987 k; 3988 l= 2; 3989 } 3990 3991 delete [] degsf; 3992 3993 return result; 3994 } 3995 3996 static inline 3997 int Hensel_P (const CanonicalForm & UU, CFArray & G, const Evaluation & AA, 3998 const Variable & x, const CFArray& LeadCoeffs ) 3826 3999 { 3827 4000 CFList factors; 3828 4001 factors.append (G[1]); 3829 4002 factors.append (G[2]); 4003 4004 CFMap NN, MM; 4005 Evaluation A= optimize4Lift (UU, MM, NN, AA); 4006 4007 CanonicalForm U= MM (UU); 4008 CFArray LCs= CFArray (1,2); 4009 LCs [1]= MM (LeadCoeffs [1]); 4010 LCs [2]= MM (LeadCoeffs [2]); 4011 3830 4012 CFList evaluation; 3831 4013 for (int i= A.min(); i <= A.max(); i++) … … 3833 4015 CFList UEval; 3834 4016 CanonicalForm shiftedU= myShift2Zero (U, UEval, evaluation); 4017 4018 if (size (shiftedU)/getNumVars (U) > 500) 4019 return 1; 4020 3835 4021 CFArray shiftedLCs= CFArray (2); 3836 4022 CFList shiftedLCsEval1, shiftedLCsEval2; 3837 4023 shiftedLCs[0]= myShift2Zero (LCs[1], shiftedLCsEval1, evaluation); 3838 4024 shiftedLCs[1]= myShift2Zero (LCs[2], shiftedLCsEval2, evaluation); 3839 CanonicalForm LcUEval= LC (UEval.getFirst(), x);3840 4025 factors.insert (1); 3841 4026 int liftBound= degree (UEval.getLast(), 2) + 1; … … 3849 4034 lcs, false); 3850 4035 3851 bool noChange= true;3852 4036 for (CFListIterator i= factors; i.hasItem(); i++) 3853 4037 { 3854 if (degree (i.getItem(), 2) != 0) 3855 noChange= false; 3856 } 3857 if (noChange) 3858 return !noChange; 4038 if (!fdivides (i.getItem(), UEval.getFirst())) 4039 return 0; 4040 } 4041 3859 4042 int * liftBounds; 3860 noChange= false;4043 bool noOneToOne; 3861 4044 if (U.level() > 2) 3862 4045 { … … 3866 4049 liftBounds[i]= degree (shiftedU, Variable (i + 2)) + 1; 3867 4050 factors= henselLift2 (UEval, factors, liftBounds, U.level()  1, false, 3868 shiftedLCsEval1, shiftedLCsEval2, Pi, diophant, noChange); 3869 delete [] liftBounds; 3870 if (noChange) 3871 return !noChange; 4051 shiftedLCsEval1, shiftedLCsEval2, Pi, diophant, 4052 noOneToOne); 4053 delete [] liftBounds; 4054 if (noOneToOne) 4055 return 0; 3872 4056 } 3873 4057 G[1]= factors.getFirst(); … … 3875 4059 G[1]= myReverseShift (G[1], evaluation); 3876 4060 G[2]= myReverseShift (G[2], evaluation); 3877 return true; 4061 G[1]= NN (G[1]); 4062 G[2]= NN (G[2]); 4063 return 1; 3878 4064 } 3879 4065 … … 3881 4067 bool findeval_P (const CanonicalForm & F, const CanonicalForm & G, 3882 4068 CanonicalForm & Fb, CanonicalForm & Gb, CanonicalForm & Db, 3883 FFREvaluation & b, int delta, int degF, int degG, int maxeval,3884 int & count )4069 REvaluation & b, int delta, int degF, int degG, int maxeval, 4070 int & count, int& k, int bound, int& l) 3885 4071 { 3886 4072 if( count == 0 && delta != 0) … … 3889 4075 return false; 3890 4076 } 4077 CanonicalForm evalPoint; 3891 4078 if (count > 0) 3892 4079 { 3893 b.nextpoint(); 4080 b.nextpoint(k); 4081 l++; 4082 if (l > bound) 4083 { 4084 l= 1; 4085 k++; 4086 if (k > tmax (F.level(), G.level())  1) 4087 return false; 4088 b.nextpoint (k); 4089 } 3894 4090 if (count++ > maxeval) 3895 4091 return false; … … 3913 4109 } 3914 4110 } 3915 b.nextpoint(); 4111 if (k == 0) 4112 k++; 4113 b.nextpoint(k); 4114 l++; 4115 if (l > bound) 4116 { 4117 l= 1; 4118 k++; 4119 if (k > tmax (F.level(), G.level())  1) 4120 return false; 4121 b.nextpoint (k); 4122 } 3916 4123 if( count++ > maxeval ) 3917 4124 return false; … … 3922 4129 static int maxNumEval= 200; 3923 4130 static int sizePerVars1= 500; //try dense gcd if size/#variables is bigger 3924 static int sizePerVars2= 290; //try psr gcd if size/#variables is less3925 4131 3926 4132 /// Extended Zassenhaus GCD for finite fields … … 3943 4149 (int)(ilog2(getCharacteristic())*log2exp))*2, maxNumEval); 3944 4150 count= 0; // number of eval. used 3945 FFREvaluation b, bt;3946 bool gcdfound = false;4151 REvaluation b, bt; 4152 int gcdfound = 0; 3947 4153 Variable x = Variable(1); 3948 4154 … … 3971 4177 3972 4178 int maxNumVars= tmax (getNumVars (F), getNumVars (G)); 3973 Variable a, bv;4179 Variable a, oldA; 3974 4180 int sizeF= size (F); 3975 4181 int sizeG= size (G); … … 3990 4196 } 3991 4197 4198 bool passToGF= false; 4199 bool extOfExt= false; 4200 int p= getCharacteristic(); 4201 bool algExtension= (hasFirstAlgVar(F,a)  hasFirstAlgVar(G,a)); 4202 int k= 1; 4203 CanonicalForm primElem, imPrimElem; 4204 CFList source, dest; 4205 if (p < 50 && CFFactory::gettype() != GaloisFieldDomain && !algExtension) 4206 { 4207 if (p == 2) 4208 setCharacteristic (2, 6, 'Z'); 4209 else if (p == 3) 4210 setCharacteristic (3, 4, 'Z'); 4211 else if (p == 5  p == 7) 4212 setCharacteristic (p, 3, 'Z'); 4213 else 4214 setCharacteristic (p, 2, 'Z'); 4215 passToGF= true; 4216 F= F.mapinto(); 4217 G= G.mapinto(); 4218 maxeval= 2*ipower (p, getGFDegree()); 4219 } 4220 else if (CFFactory::gettype() == GaloisFieldDomain && 4221 ipower (p , getGFDegree()) < 50) 4222 { 4223 k= getGFDegree(); 4224 if (ipower (p, 2*k) > 50) 4225 setCharacteristic (p, 2*k, gf_name); 4226 else 4227 setCharacteristic (p, 3*k, gf_name); 4228 F= GFMapUp (F, k); 4229 G= GFMapUp (G, k); 4230 maxeval= tmin (2*ipower (p, getGFDegree()), maxNumEval); 4231 } 4232 else if (p < 50 && algExtension && !CFFactory::gettype() == GaloisFieldDomain) 4233 { 4234 int d= degree (getMipo (a)); 4235 oldA= a; 4236 Variable v2; 4237 if (p == 2 && d < 6) 4238 { 4239 zz_p::init (p); 4240 bool primFail= false; 4241 Variable vBuf; 4242 primElem= primitiveElement (a, vBuf, primFail); 4243 ASSERT (!primFail, "failure in integer factorizer"); 4244 if (d < 3) 4245 { 4246 zz_pX NTLIrredpoly; 4247 BuildIrred (NTLIrredpoly, d*3); 4248 CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1)); 4249 v2= rootOf (newMipo); 4250 } 4251 else 4252 { 4253 zz_pX NTLIrredpoly; 4254 BuildIrred (NTLIrredpoly, d*2); 4255 CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1)); 4256 v2= rootOf (newMipo); 4257 } 4258 imPrimElem= mapPrimElem (primElem, a, v2); 4259 extOfExt= true; 4260 } 4261 else if ((p == 3 && d < 4)  ((p == 5  p == 7) && d < 3)) 4262 { 4263 zz_p::init (p); 4264 bool primFail= false; 4265 Variable vBuf; 4266 primElem= primitiveElement (a, vBuf, primFail); 4267 ASSERT (!primFail, "failure in integer factorizer"); 4268 zz_pX NTLIrredpoly; 4269 BuildIrred (NTLIrredpoly, d*2); 4270 CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1)); 4271 v2= rootOf (newMipo); 4272 imPrimElem= mapPrimElem (primElem, a, v2); 4273 extOfExt= true; 4274 } 4275 if (extOfExt) 4276 { 4277 maxeval= tmin (2*ipower (p, degree (getMipo (v2))), maxNumEval); 4278 F= mapUp (F, a, v2, primElem, imPrimElem, source, dest); 4279 G= mapUp (G, a, v2, primElem, imPrimElem, source, dest); 4280 a= v2; 4281 } 4282 } 4283 3992 4284 lcF = LC( F, x ); lcG = LC( G, x ); 3993 4285 lcD = gcd( lcF, lcG ); … … 3995 4287 delta = 0; 3996 4288 degF = degree( F, x ); degG = degree( G, x ); 3997 if(hasFirstAlgVar(F,a)) 3998 { 3999 if(hasFirstAlgVar(G,bv)) 4000 { 4001 if(bv.level() > a.level()) 4002 a = bv; 4003 } 4004 b = FFREvaluation( 2, tmax(F.level(), G.level()), AlgExtRandomF( a ) ); 4005 } 4006 else // F not in extension 4007 { 4008 if(hasFirstAlgVar(G,a)) 4009 b = FFREvaluation( 2, tmax(F.level(), G.level()), AlgExtRandomF( a ) ); 4289 4290 if(hasFirstAlgVar(G,a)) 4291 b = REvaluation( 2, tmax(F.level(), G.level()), AlgExtRandomF( a ) ); 4292 else 4293 { // both not in extension given by algebraic variable 4294 if (CFFactory::gettype() != GaloisFieldDomain) 4295 b = REvaluation( 2, tmax(F.level(), G.level()), FFRandom() ); 4010 4296 else 4011 { // both not in extension given by algebraic variable 4012 if (CFFactory::gettype() != GaloisFieldDomain) 4013 b = FFREvaluation( 2, tmax(F.level(), G.level()), FFRandom() ); 4014 else 4015 b = FFREvaluation( 2, tmax(F.level(), G.level()), GFRandom() ); 4016 } 4017 } 4297 b = REvaluation( 2, tmax(F.level(), G.level()), GFRandom() ); 4298 } 4299 4018 4300 CanonicalForm cand; 4301 CanonicalForm result; 4302 int o, t; 4303 o= 0; 4304 t= 1; 4305 int goodPointCount= 0; 4019 4306 while( !gcdfound ) 4020 4307 { 4021 if( !findeval_P( F, G, Fb, Gb, Db, b, delta, degF, degG, maxeval, count )) 4308 if( !findeval_P( F, G, Fb, Gb, Db, b, delta, degF, degG, maxeval, count, o, 4309 maxeval/maxNumVars, t )) 4022 4310 { // too many eval. used > try another method 4023 if (sizeF/maxNumVars < sizePerVars2 && sizeG/maxNumVars < sizePerVars2) 4024 { 4025 CanonicalForm dummy1, dummy2; 4026 Variable y= Variable (tmax (F.level(), G.level())); 4027 Variable z= Variable (smallestDegLev); 4028 dummy1= swapvar (F, x, y); 4029 dummy2= swapvar (G, x, y); 4030 if (size (psr (dummy1, dummy2, y)) <= tmin (size (F), size (G))) 4031 { 4032 Off (SW_USE_EZGCD_P); 4033 CanonicalForm result= swapvar (gcd (dummy1,dummy2), x, y); 4034 On (SW_USE_EZGCD_P); 4035 return N (d*result); 4036 } 4037 } 4038 if (hasFirstAlgVar (F, a)  hasFirstAlgVar (G, a)) 4039 return N (d*sparseGCDFq (F, G, a)); 4040 if (CFFactory::gettype() == GaloisFieldDomain) 4041 return N (d*GCD_GF (F, G)); 4042 else 4043 return N (d*sparseGCDFp (F,G)); 4311 Off (SW_USE_EZGCD_P); 4312 result= gcd (F,G); 4313 On (SW_USE_EZGCD_P); 4314 if (passToGF) 4315 { 4316 Variable alpha= rootOf (gf_mipo); 4317 setCharacteristic (p); 4318 result= GF2FalphaRep (result, alpha); 4319 } 4320 if (k > 1) 4321 { 4322 result= GFMapDown (result, k); 4323 setCharacteristic (p, k, gf_name); 4324 } 4325 if (extOfExt) 4326 result= mapDown (result, primElem, imPrimElem, oldA, dest, source); 4327 return N (d*result); 4044 4328 } 4045 4329 delta = degree( Db ); 4046 4330 if( delta == 0 ) 4331 { 4332 if (passToGF) 4333 setCharacteristic (p); 4334 if (k > 1) 4335 setCharacteristic (p, k, gf_name); 4047 4336 return N (d); 4337 } 4048 4338 while( true ) 4049 4339 { 4050 4340 bt = b; 4051 if( !findeval_P(F,G,Fbt,Gbt,Dbt, bt, delta, degF, degG, maxeval, count )) 4341 if( !findeval_P(F,G,Fbt,Gbt,Dbt, bt, delta, degF, degG, maxeval, count, o, 4342 maxeval/maxNumVars, t )) 4052 4343 { // too many eval. used > try another method 4053 if (sizeF/maxNumVars < sizePerVars2 && sizeG/maxNumVars < sizePerVars2) 4054 { 4055 CanonicalForm dummy1, dummy2; 4056 Variable y= Variable (tmax (F.level(), G.level())); 4057 Variable z= Variable (smallestDegLev); 4058 dummy1= swapvar (F, x, y); 4059 dummy2= swapvar (G, x, y); 4060 if (size (psr (dummy1, dummy2, y)) <= tmin (size (F), size (G))) 4061 { 4062 Off (SW_USE_EZGCD_P); 4063 CanonicalForm result= swapvar (gcd (dummy1,dummy2), x, y); 4064 On (SW_USE_EZGCD_P); 4065 return N (d*result); 4066 } 4067 } 4068 if (hasFirstAlgVar (F, a)  hasFirstAlgVar (G, a)) 4069 return N (d*sparseGCDFq (F, G, a)); 4070 if (CFFactory::gettype() == GaloisFieldDomain) 4071 return N (d*GCD_GF (F, G)); 4072 else 4073 return N (d*sparseGCDFp (F,G)); 4344 Off (SW_USE_EZGCD_P); 4345 result= gcd (F,G); 4346 On (SW_USE_EZGCD_P); 4347 if (passToGF) 4348 { 4349 Variable alpha= rootOf (gf_mipo); 4350 setCharacteristic (p); 4351 result= GF2FalphaRep (result, alpha); 4352 } 4353 if (k > 1) 4354 { 4355 result= GFMapDown (result, k); 4356 setCharacteristic (p, k, gf_name); 4357 } 4358 if (extOfExt) 4359 result= mapDown (result, primElem, imPrimElem, oldA, dest, source); 4360 return N (d*result); 4074 4361 } 4075 4362 int dd = degree( Dbt ); 4076 4363 if( dd == 0 ) 4364 { 4365 if (passToGF) 4366 setCharacteristic (p); 4367 if (k > 1) 4368 setCharacteristic (p, k, gf_name); 4077 4369 return N (d); 4370 } 4078 4371 if( dd == delta ) 4079 break; 4372 { 4373 goodPointCount++; 4374 if (goodPointCount == 5) 4375 break; 4376 } 4080 4377 if( dd < delta ) 4081 4378 { 4379 goodPointCount= 0; 4082 4380 delta = dd; 4083 4381 b = bt; 4084 4382 Db = Dbt; Fb = Fbt; Gb = Gbt; 4085 4383 } 4086 } 4087 if( degF <= degG && delta == degF && fdivides( F, G ) ) 4088 return N (d*F); 4089 if( degG < degF && delta == degG && fdivides( G, F ) ) 4090 return N (d*G); 4384 if (delta == degF) 4385 { 4386 if (degF <= degG && fdivides (F, G)) 4387 { 4388 if (passToGF) 4389 { 4390 Variable alpha= rootOf (gf_mipo); 4391 setCharacteristic (p); 4392 F= GF2FalphaRep (F, alpha); 4393 } 4394 if (k > 1) 4395 { 4396 F= GFMapDown (F, k); 4397 setCharacteristic (p, k, gf_name); 4398 } 4399 if (extOfExt) 4400 F= mapDown (F, primElem, imPrimElem, oldA, dest, source); 4401 return N (d*F); 4402 } 4403 else 4404 delta; 4405 } 4406 else if (delta == degG) 4407 { 4408 if (degG <= degF && fdivides (G, F)) 4409 { 4410 if (passToGF) 4411 { 4412 Variable alpha= rootOf (gf_mipo); 4413 setCharacteristic (p); 4414 G= GF2FalphaRep (G, alpha); 4415 } 4416 if (k > 1) 4417 { 4418 G= GFMapDown (G, k); 4419 setCharacteristic (p, k, gf_name); 4420 } 4421 if (extOfExt) 4422 G= mapDown (G, primElem, imPrimElem, oldA, dest, source); 4423 return N (d*G); 4424 } 4425 else 4426 delta; 4427 } 4428 } 4091 4429 if( delta != degF && delta != degG ) 4092 4430 { … … 4121 4459 else // special case handling 4122 4460 { 4123 if (sizeF/maxNumVars < sizePerVars2 && sizeG/maxNumVars < sizePerVars2) 4124 { 4125 CanonicalForm dummy1, dummy2; 4126 Variable y= Variable (tmax (F.level(), G.level())); 4127 Variable z= Variable (smallestDegLev); 4128 dummy1= swapvar (F, x, y); 4129 dummy2= swapvar (G, x, y); 4130 if (size (psr (dummy1, dummy2, y)) <= tmin (size (F), size (G))) 4131 { 4132 Off (SW_USE_EZGCD_P); 4133 CanonicalForm result= swapvar (gcd (dummy1,dummy2), x, y); 4134 On (SW_USE_EZGCD_P); 4135 return N (d*result); 4136 } 4137 } 4138 if (hasFirstAlgVar (F, a)  hasFirstAlgVar (G, a)) 4139 return N (d*sparseGCDFq (F, G, a)); 4140 if (CFFactory::gettype() == GaloisFieldDomain) 4141 return N (d*GCD_GF (F, G)); 4142 else 4143 return N (d*sparseGCDFp (F,G)); 4461 Off (SW_USE_EZGCD_P); 4462 result= gcd (F,G); 4463 On (SW_USE_EZGCD_P); 4464 if (passToGF) 4465 { 4466 Variable alpha= rootOf (gf_mipo); 4467 setCharacteristic (p); 4468 result= GF2FalphaRep (result, alpha); 4469 } 4470 if (k > 1) 4471 { 4472 result= GFMapDown (result, k); 4473 setCharacteristic (p, k, gf_name); 4474 } 4475 if (extOfExt) 4476 result= mapDown (result, primElem, imPrimElem, oldA, dest, source); 4477 return N (d*result); 4144 4478 } 4145 4479 DD[2] = DD[2] * ( b( lcDD[2] ) / lc( DD[2] ) ); 4146 4480 DD[1] = DD[1] * ( b( lcDD[1] ) / lc( DD[1] ) ); 4147 4481 4482 if (size (B*lcDD[2])/maxNumVars > sizePerVars1) 4483 { 4484 if (algExtension) 4485 { 4486 result= GCD_Fp_extension (F, G, a); 4487 if (extOfExt) 4488 result= mapDown (result, primElem, imPrimElem, oldA, dest, source); 4489 return N (d*result); 4490 } 4491 if (CFFactory::gettype() == GaloisFieldDomain) 4492 { 4493 result= GCD_GF (F, G); 4494 if (passToGF) 4495 { 4496 Variable alpha= rootOf (gf_mipo); 4497 setCharacteristic (p); 4498 result= GF2FalphaRep (result, alpha); 4499 } 4500 if (k > 1) 4501 { 4502 result= GFMapDown (result, k); 4503 setCharacteristic (p, k, gf_name); 4504 } 4505 return N (d*result); 4506 } 4507 else 4508 return N (d*GCD_small_p (F,G)); 4509 } 4510 4148 4511 gcdfound= Hensel_P (B*lcD, DD, b, x, lcDD); 4149 4512 4150 if( gcdfound ) 4513 if (gcdfound == 1) 4514 { 4515 Off (SW_USE_EZGCD_P); 4516 result= gcd (F,G); 4517 On (SW_USE_EZGCD_P); 4518 if (passToGF) 4519 { 4520 Variable alpha= rootOf (gf_mipo); 4521 setCharacteristic (p); 4522 result= GF2FalphaRep (result, alpha); 4523 } 4524 if (k > 1) 4525 { 4526 result= GFMapDown (result, k); 4527 setCharacteristic (p, k, gf_name); 4528 } 4529 if (extOfExt) 4530 result= mapDown (result, primElem, imPrimElem, oldA, dest, source); 4531 return N (d*result); 4532 } 4533 4534 if (gcdfound == 1) 4151 4535 { 4152 4536 cand = DD[2] / content( DD[2], Variable(1) ); 4153 if( B_is_F ) 4154 gcdfound = fdivides( cand, G ); 4155 else 4156 gcdfound = fdivides( cand, F ); 4157 } 4158 } 4159 delta= 0; 4537 gcdfound = fdivides( cand, G ) && fdivides ( cand, F ); 4538 4539 if (passToGF && gcdfound) 4540 { 4541 Variable alpha= rootOf (gf_mipo); 4542 setCharacteristic (p); 4543 cand= GF2FalphaRep (cand, alpha); 4544 } 4545 if (k > 1 && gcdfound) 4546 { 4547 cand= GFMapDown (cand, k); 4548 setCharacteristic (p, k, gf_name); 4549 } 4550 if (extOfExt && gcdfound) 4551 cand= mapDown (cand, primElem, imPrimElem, oldA, dest, source); 4552 } 4553 } 4554 delta; 4555 goodPointCount= 0; 4160 4556 } 4161 4557 return N (d*cand);
Note: See TracChangeset
for help on using the changeset viewer.