1 | #include "config.h" |
---|
2 | #ifdef HAVE_CSTDIO |
---|
3 | #include <cstdio> |
---|
4 | #else |
---|
5 | #include <stdio.h> |
---|
6 | #endif |
---|
7 | #ifndef NOSTREAMIO |
---|
8 | #ifdef HAVE_IOSTREAM |
---|
9 | #include <iostream> |
---|
10 | #define ISTREAM std::istream |
---|
11 | #define OSTREAM std::ostream |
---|
12 | #define CERR std::cerr |
---|
13 | #elif defined(HAVE_IOSTREAM_H) |
---|
14 | #include <iostream.h> |
---|
15 | #define ISTREAM istream |
---|
16 | #define OSTREAM ostream |
---|
17 | #define CERR cerr |
---|
18 | #endif |
---|
19 | #endif |
---|
20 | #include <math.h> |
---|
21 | #include <cf_factory.h> |
---|
22 | #include <variable.h> |
---|
23 | #include <canonicalform.h> |
---|
24 | #include <cf_iter.h> |
---|
25 | #include <cf_random.h> |
---|
26 | #include <cf_algorithm.h> |
---|
27 | #include <cf_map.h> |
---|
28 | #include <cf_util.h> |
---|
29 | |
---|
30 | static CanonicalForm contentWRT(const CanonicalForm & F, const int lev); |
---|
31 | static int degWRT(CanonicalForm F, const int lev); |
---|
32 | static CanonicalForm lcoefWRT( const CanonicalForm & F, const int lev); |
---|
33 | static CanonicalForm simpleGCD(const CanonicalForm & A, const CanonicalForm & B); |
---|
34 | CanonicalForm newGCD(CanonicalForm A, CanonicalForm B); |
---|
35 | static CanonicalForm GFPowDown(const CanonicalForm & A, int k); |
---|
36 | static CanonicalForm GFPowUp(const CanonicalForm & A, int k); |
---|
37 | static CanonicalForm GFMapUp(const CanonicalForm & A, int k); |
---|
38 | static CanonicalForm GFMapDown(const CanonicalForm & A, int k); |
---|
39 | |
---|
40 | |
---|
41 | |
---|
42 | |
---|
43 | static CanonicalForm |
---|
44 | contentWRT0(const CanonicalForm & F, const int lev, const int lev0) |
---|
45 | // Computes the content of a polynomial, considering the variables of level |
---|
46 | // larger than lev as as parameters. |
---|
47 | { |
---|
48 | if(F.inBaseDomain() || (F.mvar().level() <= lev)) |
---|
49 | { |
---|
50 | return(1); |
---|
51 | } |
---|
52 | else |
---|
53 | { |
---|
54 | if (F.isUnivariate()) |
---|
55 | { |
---|
56 | return(F); |
---|
57 | } |
---|
58 | else |
---|
59 | { |
---|
60 | CanonicalForm pol; |
---|
61 | CFIterator i = CFIterator(F); |
---|
62 | CanonicalForm c = 0; |
---|
63 | for(; i.hasTerms(); i++) |
---|
64 | { |
---|
65 | CanonicalForm cc=i.coeff(); |
---|
66 | if (cc.level() > lev0) |
---|
67 | pol = contentWRT0(cc, lev, lev0-1 ); |
---|
68 | else |
---|
69 | pol = contentWRT(cc, lev -1); |
---|
70 | c = gcd(c, pol); |
---|
71 | if(c.isOne()) |
---|
72 | { |
---|
73 | return c; |
---|
74 | } |
---|
75 | } |
---|
76 | return c; |
---|
77 | } |
---|
78 | } |
---|
79 | } |
---|
80 | |
---|
81 | static CanonicalForm |
---|
82 | contentWRT(const CanonicalForm & F, const int lev) |
---|
83 | // Computes the content of a polynomial, considering the variables of level |
---|
84 | // larger than lev as as parameters. |
---|
85 | { |
---|
86 | //if (lev <0) |
---|
87 | // cout << "contentWRT:" <<F << " " << lev <<"\n";//printf("level=%d\n",lev); |
---|
88 | //int debug = 0; |
---|
89 | //if(debug == 1)cout << "contentWRT input - F: " << F << " lev: " << lev << endl; |
---|
90 | if(F.inBaseDomain() || (F.mvar().level() <= lev)) |
---|
91 | { |
---|
92 | //if(debug == 1)cout << "output 1:" << F << endl; |
---|
93 | return(1); |
---|
94 | } |
---|
95 | else |
---|
96 | { |
---|
97 | if (F.isUnivariate()) |
---|
98 | { |
---|
99 | //if(debug == 1)cout << "output 2:" << F << endl; |
---|
100 | return(F); |
---|
101 | } |
---|
102 | else |
---|
103 | { |
---|
104 | CanonicalForm pol; |
---|
105 | CFIterator i = CFIterator(F); |
---|
106 | CanonicalForm c = 0; |
---|
107 | for(; i.hasTerms(); i++) |
---|
108 | { |
---|
109 | CanonicalForm cc=i.coeff(); |
---|
110 | if (cc.level() > lev) |
---|
111 | pol = contentWRT0(cc, lev, cc.level()); |
---|
112 | else |
---|
113 | pol = contentWRT(i.coeff(), lev - 1); |
---|
114 | //if(debug == 1)cout << "c: " << c << " - pol : " << pol << endl; |
---|
115 | c = gcd(c, pol); |
---|
116 | //if(debug == 1)cout << "c temp:" << c << endl; |
---|
117 | if(c.isOne()) |
---|
118 | { |
---|
119 | //if(debug == 1)cout << "output 3:" << c << endl; |
---|
120 | return c; |
---|
121 | } |
---|
122 | } |
---|
123 | //if(debug == 1)cout << "output 4:" << c << endl; |
---|
124 | return c; |
---|
125 | } |
---|
126 | } |
---|
127 | } |
---|
128 | |
---|
129 | // Degree of F wrt all variables of level smaller than or equal to lev |
---|
130 | static int degWRT(CanonicalForm F, const int lev) |
---|
131 | { |
---|
132 | int deg = F.degree(lev); |
---|
133 | for(int i = lev; i >= 1; i--) |
---|
134 | { |
---|
135 | F = F.LC(i); |
---|
136 | deg = deg + F.degree(i - 1); |
---|
137 | } |
---|
138 | return deg; |
---|
139 | } |
---|
140 | |
---|
141 | // Leading coefficient of F wrt all variables of level smaller than or equal to lev. |
---|
142 | static CanonicalForm lcoefWRT(const CanonicalForm & F, const int lev) |
---|
143 | { |
---|
144 | CanonicalForm FL; |
---|
145 | FL = F; |
---|
146 | |
---|
147 | for(int i = lev; i >= 1; i--) |
---|
148 | { |
---|
149 | FL = FL.LC(i); |
---|
150 | } |
---|
151 | return FL; |
---|
152 | } |
---|
153 | |
---|
154 | |
---|
155 | static CanonicalForm |
---|
156 | newtonInterp(const CanonicalForm alpha, const CanonicalForm u, const CanonicalForm newtonPoly, const CanonicalForm oldInterPoly, const Variable & x) |
---|
157 | // Newton interpolation - Incremental algorithm. |
---|
158 | // Given a list of values alpha_i and a list of polynomials u_i, 1 <= i <= n, |
---|
159 | // computes the interpolation polynomial assuming that |
---|
160 | // the polynomials in u are the results of evaluating the variabe x |
---|
161 | // of the unknown polynomial at the alpha values. |
---|
162 | // This incremental version receives only the values of alpha_n and u_n and |
---|
163 | // the previous interpolation polynomial for points 1 <= i <= n-1, and computes |
---|
164 | // the polynomial interpolating in all the points. |
---|
165 | // newtonPoly must be equal to (x - alpha_1) * ... * (x - alpha_{n-1}) |
---|
166 | { |
---|
167 | /* |
---|
168 | int debug = 0; |
---|
169 | if(debug) |
---|
170 | { |
---|
171 | cout << "newtonInterpIncremental input - variable " << x <<endl; |
---|
172 | cout << "newtonInterpIncremental input - alpha:" << alpha << " - u: " << u << endl; |
---|
173 | cout << "newtonInterpIncremental input - newtonPoly: " << newtonPoly << " - oldInterPoly: " << oldInterPoly << endl; |
---|
174 | } |
---|
175 | */ |
---|
176 | |
---|
177 | CanonicalForm interPoly; |
---|
178 | |
---|
179 | interPoly = oldInterPoly + ((u - oldInterPoly(alpha, x)) / newtonPoly(alpha, x)) * newtonPoly; |
---|
180 | /*cout << "newtonPoly (alpha,x)= " << newtonPoly (alpha,x) << "\n"; |
---|
181 | cout << "1/newtonPoly (alpha,x)= " << 1/newtonPoly (alpha,x) << "\n"; |
---|
182 | cout << "newtonPoly (alpha,x)*(1/newtonPoly (alpha,x))= " << newtonPoly (alpha,x)*(1/newtonPoly (alpha,x)) << "\n";*/ |
---|
183 | //if(debug)cout << "newtonInterpIncremental output:" << interPoly << endl; |
---|
184 | return interPoly; |
---|
185 | } |
---|
186 | |
---|
187 | |
---|
188 | // Univariate GCD |
---|
189 | // Naive implementation of Euclidean Algorithm. |
---|
190 | // Should be used only for GCD of univariate polys over finite fields, |
---|
191 | // where there is no coefficient growth. |
---|
192 | static CanonicalForm simpleGCD(const CanonicalForm & A, const CanonicalForm & B) |
---|
193 | { |
---|
194 | //int debug = 0; |
---|
195 | CanonicalForm P1, P2; |
---|
196 | /* |
---|
197 | if(debug == 1) |
---|
198 | { |
---|
199 | cout << "simpleGCD input" << endl; |
---|
200 | cout << "A: " << A << endl; |
---|
201 | cout << "B: " << B << endl; |
---|
202 | } |
---|
203 | */ |
---|
204 | if (A.degree() > B.degree()) |
---|
205 | { |
---|
206 | P1 = A; |
---|
207 | P2 = B; |
---|
208 | } |
---|
209 | else |
---|
210 | { |
---|
211 | P1 = B; |
---|
212 | P2 = A; |
---|
213 | } |
---|
214 | CanonicalForm rem; |
---|
215 | while(!P2.isZero()) |
---|
216 | { |
---|
217 | // cout << "P1: " << P1 << " - P2: " << P2 << endl; |
---|
218 | rem = P1 % P2; |
---|
219 | // cout << "rem: " << rem << endl; |
---|
220 | P1 = P2; |
---|
221 | P2 = rem; |
---|
222 | } |
---|
223 | //if(debug == 1)cout << "simpleGCD output: " << P1 << endl; |
---|
224 | return(P1/lc(P1)); |
---|
225 | } |
---|
226 | |
---|
227 | |
---|
228 | static CanonicalForm GFPowDown(const CanonicalForm & A, int k) |
---|
229 | // Replaces all coefficients of A in the ground field by their k-th roots. |
---|
230 | // It assumes that the k-th roots exist. |
---|
231 | // This procedure is used to map down a polynomial from an extension field. |
---|
232 | { |
---|
233 | CanonicalForm result = 0; |
---|
234 | int i, j; |
---|
235 | int fieldSize = ipower(getCharacteristic(), getGFDegree()); |
---|
236 | CanonicalForm g; |
---|
237 | for(i = 0; i <= degree(A); i++) |
---|
238 | { |
---|
239 | if(!A[i].isZero()) |
---|
240 | { |
---|
241 | //cout << "i: " << i << endl; |
---|
242 | //cout << "A[i]: " << A[i] << endl; |
---|
243 | if(A[i].inBaseDomain()) |
---|
244 | { |
---|
245 | //cout << "GFPowDown - inBaseDomain" << endl; |
---|
246 | g = getGFGenerator(); |
---|
247 | j = 0; |
---|
248 | while(((power(g, j * k)*1) != (A[i]*1)) && (j <= fieldSize)) |
---|
249 | { |
---|
250 | //cout << "power(g, j * k)*1: " << power(g, j * k)*1 << " - A[i]*1: " << A[i]*1 << endl; |
---|
251 | j++; |
---|
252 | } |
---|
253 | if(j == fieldSize + 1) return(-1); |
---|
254 | result = result + power(g, j) * power(A.mvar(), i); |
---|
255 | //cout << "power(A[i], k): " << power(A[i], k) << endl; |
---|
256 | //cout << "result: " << result << endl; |
---|
257 | } |
---|
258 | else |
---|
259 | { |
---|
260 | //cout << "not inBaseDomain" << endl; |
---|
261 | result = result + GFPowDown(A[i], k) * power(A.mvar(), i);; |
---|
262 | // power(pol[i], k) * power(A.mvar(), i); |
---|
263 | //pol[i] = GFMap(pol[i], k); |
---|
264 | } |
---|
265 | } |
---|
266 | } |
---|
267 | return(result); |
---|
268 | } |
---|
269 | |
---|
270 | static CanonicalForm GFPowUp(const CanonicalForm & A, int k) |
---|
271 | // Raises all coefficients of A in the ground field to the power k. |
---|
272 | // Used for maping a polynomial to an extension field. |
---|
273 | { |
---|
274 | //cout << "A:" << A <<"\n"; |
---|
275 | CanonicalForm result = 0; |
---|
276 | int i; |
---|
277 | for(i = 0; i <= degree(A); i++) |
---|
278 | { |
---|
279 | if(!A[i].isZero()) |
---|
280 | { |
---|
281 | //cout << "i: " << i << endl; |
---|
282 | //cout << "A[i]: " << A[i] << endl; |
---|
283 | if(A[i].inBaseDomain()) |
---|
284 | { |
---|
285 | //cout << "inBaseDomain" << endl; |
---|
286 | result = result + power(A[i], k) * power(A.mvar(), i); |
---|
287 | //cout << "power(A[i], k): " << power(A[i], k) << endl; |
---|
288 | //cout << "result: " << result << endl; |
---|
289 | } |
---|
290 | else |
---|
291 | { |
---|
292 | //cout << "not inBaseDomain" << endl; |
---|
293 | result = result + GFPowUp(A[i], k) * power(A.mvar(), i);; |
---|
294 | } |
---|
295 | } |
---|
296 | } |
---|
297 | //cout << "result:" << result <<"\n"; |
---|
298 | return(result); |
---|
299 | } |
---|
300 | |
---|
301 | CanonicalForm GFMapUp(const CanonicalForm & A, int k) |
---|
302 | // Maps all coefficients of A to the base domain, asumming A is in GF(p^k). |
---|
303 | // The current base domain must be GF(p^j), with j a multiple of k. |
---|
304 | { |
---|
305 | //int p = getCharacteristic(); |
---|
306 | int expon = getGFDegree(); |
---|
307 | int extExp = expon / k; |
---|
308 | // Assumes that we are using Conway polynomials |
---|
309 | return(GFPowUp(A, extExp)); |
---|
310 | } |
---|
311 | |
---|
312 | CanonicalForm GFMapDown(const CanonicalForm & A, int k) |
---|
313 | // Maps all coefficients of A from the base domain to GF(p^k). |
---|
314 | // The current base domain must be GF(p^j), with j a multiple of k. |
---|
315 | { |
---|
316 | //int p = getCharacteristic(); |
---|
317 | int expon = getGFDegree(); |
---|
318 | //cout << "Expon: " << expon << endl; |
---|
319 | int extExp = expon / k; |
---|
320 | // Assumes that we are using Conway polynomials |
---|
321 | return(GFPowDown(A, extExp)); |
---|
322 | } |
---|
323 | |
---|
324 | CanonicalForm newGCD(CanonicalForm A, CanonicalForm B) |
---|
325 | // Computes the GCD of two polynomials over a prime field. |
---|
326 | // Based on Algorithm 7.2 from "Algorithms for Computer Algebra" |
---|
327 | { |
---|
328 | //int debug = 1; |
---|
329 | //if(debug) |
---|
330 | //{ |
---|
331 | // cout << "newGCD input" << endl; |
---|
332 | // cout << "A: " << A << endl; |
---|
333 | // cout << "B: " << B << endl; |
---|
334 | //} |
---|
335 | |
---|
336 | if(A*1 == 0) |
---|
337 | { |
---|
338 | //if(debug)cout << "A Zero!" << endl; |
---|
339 | return(abs(B)); |
---|
340 | } |
---|
341 | if(B*1 == 0) |
---|
342 | { |
---|
343 | //if(debug)cout << "B Zero!" << endl; |
---|
344 | return(abs(A)); |
---|
345 | } |
---|
346 | if (A.inBaseDomain()) return A.genOne(); |
---|
347 | if (B.inBaseDomain()) return B.genOne(); |
---|
348 | |
---|
349 | CanonicalForm pol; |
---|
350 | int p = getCharacteristic(); |
---|
351 | int k = 1; |
---|
352 | int fieldSize = p; |
---|
353 | if (CFFactory::gettype() == GaloisFieldDomain) |
---|
354 | { |
---|
355 | k=getGFDegree(); |
---|
356 | fieldSize = ipower(p, k); |
---|
357 | } |
---|
358 | |
---|
359 | //if(debug) |
---|
360 | //{ |
---|
361 | // cout << "p: " << p << endl; |
---|
362 | // cout << "k: " << k << endl; |
---|
363 | // cout << "fieldSize: " << fieldSize << endl; |
---|
364 | // cout << "pow: " << ipower(p, k) << endl; |
---|
365 | //} |
---|
366 | |
---|
367 | CFMap M,N; |
---|
368 | compress(A,B,M,N); |
---|
369 | A=M(A); |
---|
370 | B=M(B); |
---|
371 | // Compute the main variable (the largest in A and B) |
---|
372 | Variable x = A.mvar(); |
---|
373 | int mv = x.level(); |
---|
374 | if(B.mvar().level() > mv) |
---|
375 | { |
---|
376 | x = B.mvar(); |
---|
377 | mv = x.level(); |
---|
378 | } |
---|
379 | //cout << "main variable " << x << endl << endl; |
---|
380 | |
---|
381 | // Checks for univariate polynomial |
---|
382 | if ( mv == 1 ) |
---|
383 | { |
---|
384 | // This call can be replaced by a faster GCD algorithm |
---|
385 | pol = simpleGCD(A, B); // Univariate GCD |
---|
386 | //cout << "newGCD output 1: " << pol << endl; |
---|
387 | return N(pol); |
---|
388 | } |
---|
389 | |
---|
390 | CanonicalForm b; |
---|
391 | CFArray bArray(0, fieldSize-1); // Stores the bs already used |
---|
392 | int sizebArray; |
---|
393 | |
---|
394 | GFRandom genGF; |
---|
395 | FFRandom genFF; |
---|
396 | int i; |
---|
397 | int used; |
---|
398 | |
---|
399 | CanonicalForm c; // gcd of the contents |
---|
400 | CanonicalForm C; |
---|
401 | CanonicalForm Cb; // gcd of Ab and Bb |
---|
402 | CanonicalForm H = 0; // for Newton Interpolation |
---|
403 | CanonicalForm newtonPoly = 1; // for Newton Interpolation |
---|
404 | CanonicalForm Cblc; |
---|
405 | |
---|
406 | // Contents and primparts of A and B |
---|
407 | CanonicalForm contA, contB; // Contents of A and B |
---|
408 | CanonicalForm Ap, Bp; // primpart of A and B |
---|
409 | contA = contentWRT(A, mv - 1); |
---|
410 | contB = contentWRT(B, mv - 1); |
---|
411 | c = simpleGCD(contA, contB); // gcd of univariate polynomials |
---|
412 | Ap = A / contA; |
---|
413 | Bp = B / contB; |
---|
414 | |
---|
415 | CanonicalForm AL, BL; // leading coefficients of A and B |
---|
416 | CanonicalForm g; // gcd of AL and BL |
---|
417 | AL = lcoefWRT(Ap, mv - 1); |
---|
418 | BL = lcoefWRT(Bp, mv - 1); |
---|
419 | g = simpleGCD(AL, BL); // gcd of univariate polynomials |
---|
420 | |
---|
421 | CanonicalForm Ab, Bb, gb; // A, B and g evaluated at b |
---|
422 | |
---|
423 | int m; // degree of the gcd of Ab and Bb |
---|
424 | int n = degWRT(Ap, mv - 1); |
---|
425 | if(degWRT(Bp, mv - 1) < n) |
---|
426 | { |
---|
427 | n = degWRT(Bp, mv - 1); |
---|
428 | } |
---|
429 | |
---|
430 | int fail = 0; |
---|
431 | int goodb; |
---|
432 | // The main loop |
---|
433 | //cout << "start loop\n"; |
---|
434 | sizebArray = 0; |
---|
435 | do |
---|
436 | { |
---|
437 | // Searches for a good b. |
---|
438 | // If there are no more possible b, the algorithm fails. |
---|
439 | goodb = 0; |
---|
440 | if(sizebArray >= fieldSize-1) |
---|
441 | { |
---|
442 | // An extension field is needed. |
---|
443 | fail = 1; |
---|
444 | } |
---|
445 | else |
---|
446 | { |
---|
447 | do |
---|
448 | { |
---|
449 | // Searches for a new element of the ground field |
---|
450 | do |
---|
451 | { |
---|
452 | // New element of the ground field |
---|
453 | if(k > 1) |
---|
454 | { |
---|
455 | b = genGF.generate(); |
---|
456 | } |
---|
457 | else |
---|
458 | { |
---|
459 | b = genFF.generate(); |
---|
460 | } |
---|
461 | //cout << "try(" << sizebArray << "):" << b; |
---|
462 | // Checks if this element was already used |
---|
463 | used = 0; |
---|
464 | for(i = 0; i < sizebArray; i++) |
---|
465 | { |
---|
466 | // Multiplication by 1 is used to prevent a bug which |
---|
467 | // produces b=a^(q-1) as a random element. |
---|
468 | if((bArray[i]*1) == (b*1)) |
---|
469 | { |
---|
470 | used = 1; |
---|
471 | //cout << " used\n"; |
---|
472 | } |
---|
473 | } |
---|
474 | //if(debug==1)cout << "b: " << b << endl; |
---|
475 | } |
---|
476 | while(used == 1); |
---|
477 | bArray[sizebArray] = b; |
---|
478 | sizebArray++; |
---|
479 | // b must not cancel the gcd of the leading coefficients |
---|
480 | if(g(b, mv) != 0) |
---|
481 | { |
---|
482 | goodb = 1; |
---|
483 | // cout << " good\n"; |
---|
484 | } |
---|
485 | else |
---|
486 | { |
---|
487 | // cout << " bad"; |
---|
488 | if(sizebArray == fieldSize - 1) |
---|
489 | { |
---|
490 | fail = 1; |
---|
491 | //cout << " out of elems " << sizebArray << "tried"; |
---|
492 | } |
---|
493 | //cout << "\n"; |
---|
494 | } |
---|
495 | } |
---|
496 | while((goodb == 0) && (fail == 0)); |
---|
497 | } |
---|
498 | if(fail) |
---|
499 | { |
---|
500 | // Algorithm fails. An extension is needed. |
---|
501 | |
---|
502 | // Computes the exponent of the ring extension so as to have enough interpolation points. |
---|
503 | int degMax; |
---|
504 | if(totaldegree(A) > totaldegree(B)) |
---|
505 | { |
---|
506 | degMax = totaldegree(A); |
---|
507 | } |
---|
508 | else |
---|
509 | { |
---|
510 | degMax = totaldegree(B); |
---|
511 | } |
---|
512 | int expon = 2; // expon <= will not extend the field |
---|
513 | while(ipower(fieldSize, expon) < degMax) |
---|
514 | { |
---|
515 | expon++; |
---|
516 | } |
---|
517 | //if(debug)cout << "Not enough elements in the base field. An extension to " << p << "^" << k*expon << " is needed." << endl; |
---|
518 | if(k > 1) |
---|
519 | { |
---|
520 | if(ipower(p,k * expon) < (1<<16)) |
---|
521 | { |
---|
522 | setCharacteristic(p, k * expon, 'b'); |
---|
523 | CanonicalForm P1 = GFMapUp(A, k); |
---|
524 | CanonicalForm P2 = GFMapUp(B, k); |
---|
525 | //cout << "newGCD(mapped):" << P1 << " , " << P2 <<"\n"; |
---|
526 | pol = newGCD(P1, P2); |
---|
527 | pol = GFMapDown(pol, k); |
---|
528 | setCharacteristic(p, k, 'a'); |
---|
529 | //if(debug)cout << "newGCD output 5: " << pol << endl; |
---|
530 | CanonicalForm temp=N(pol); |
---|
531 | temp/=temp.lc(); |
---|
532 | return temp; |
---|
533 | } |
---|
534 | else |
---|
535 | { |
---|
536 | Off(SW_USE_GCD_P); |
---|
537 | CanonicalForm temp=N(gcd( A, B)); |
---|
538 | On(SW_USE_GCD_P); |
---|
539 | return temp; |
---|
540 | } |
---|
541 | } |
---|
542 | else |
---|
543 | { |
---|
544 | if(ipower(p,k * expon) < (1<<16)) |
---|
545 | { |
---|
546 | setCharacteristic(p, k * expon, 'a'); |
---|
547 | CanonicalForm P1 = A.mapinto(); |
---|
548 | CanonicalForm P2 = B.mapinto(); |
---|
549 | pol = newGCD(P1, P2); |
---|
550 | setCharacteristic(p); |
---|
551 | //if(debug)cout << "newGCD output 4: " << pol << endl; |
---|
552 | CanonicalForm temp=N(pol); |
---|
553 | temp/=temp.lc(); |
---|
554 | return temp; |
---|
555 | } |
---|
556 | else |
---|
557 | { |
---|
558 | Off(SW_USE_GCD_P); |
---|
559 | CanonicalForm temp=N(gcd( A, B)); |
---|
560 | On(SW_USE_GCD_P); |
---|
561 | return temp; |
---|
562 | } |
---|
563 | } |
---|
564 | } |
---|
565 | else |
---|
566 | { |
---|
567 | // Evaluate the polynomials in b |
---|
568 | Ab = Ap(b, mv); |
---|
569 | Bb = Bp(b, mv); |
---|
570 | gb = g(b, mv); |
---|
571 | Cb = newGCD(Ab, Bb); |
---|
572 | //if(debug)cout << "newGCD received: " << Cb << endl; |
---|
573 | m = Cb.degree(); |
---|
574 | |
---|
575 | Cblc = Cb.lc(); |
---|
576 | Cb *= gb; |
---|
577 | Cb /= Cblc; |
---|
578 | // Test for unlucky homomorphism |
---|
579 | if(m < n) |
---|
580 | { |
---|
581 | // The degree decreased, we have to start it all over. |
---|
582 | H = Cb; |
---|
583 | newtonPoly = newtonPoly * (x - b); |
---|
584 | n = m; |
---|
585 | } |
---|
586 | else if(m == n) |
---|
587 | { |
---|
588 | // Good b |
---|
589 | H = newtonInterp(b, Cb, newtonPoly, H, x); |
---|
590 | newtonPoly = newtonPoly * (x - b); |
---|
591 | if(lcoefWRT(H, mv - 1) == g) |
---|
592 | { |
---|
593 | C = H / contentWRT(H, mv - 1); // primpart |
---|
594 | if(fdivides(C, A) && fdivides(C, B)) |
---|
595 | { |
---|
596 | //if(debug)cout << "newGCD output 2: " << c * C<< endl; |
---|
597 | return N(c * C); |
---|
598 | } |
---|
599 | else |
---|
600 | { |
---|
601 | if(m == 0) |
---|
602 | { |
---|
603 | //if(debug)cout << "newGCD output 3: " << c << endl; |
---|
604 | return N(c); |
---|
605 | } |
---|
606 | } |
---|
607 | } |
---|
608 | } |
---|
609 | } |
---|
610 | } while(1); |
---|
611 | //cout << "No way to get here!" << endl; |
---|
612 | //cout << "H:" << H << endl; |
---|
613 | // No way to get here! |
---|
614 | return H; |
---|
615 | } |
---|
616 | |
---|
617 | #if 0 |
---|
618 | main() |
---|
619 | { |
---|
620 | CanonicalForm A; |
---|
621 | CanonicalForm B; |
---|
622 | CanonicalForm pol; |
---|
623 | Variable x('x'), y('y'), z('z'); |
---|
624 | |
---|
625 | cout << "setCharacteristic(2, 4, 'a')" << endl; |
---|
626 | setCharacteristic(2, 4, 'a'); |
---|
627 | |
---|
628 | int k = getGFDegree(); |
---|
629 | int p = getCharacteristic(); |
---|
630 | int fieldSize = ipower(p, k); |
---|
631 | cout << "p: " << p << endl; |
---|
632 | cout << "GFDegree: " << k << endl; |
---|
633 | cout << "fieldSize: " << fieldSize << endl << endl; |
---|
634 | |
---|
635 | CanonicalForm g = getGFGenerator(); |
---|
636 | //CanonicalForm g = 1; |
---|
637 | A = power((x*y*y-power(x, 7)), 42) * power((power(y, 5) + g*y*y*power(x,13) + g*g), 13); |
---|
638 | B = power((x+y+1), 37) * power((power(y, 5) + g*y*y*power(x,13) + g*g), 13); |
---|
639 | cout << "A: " << A << endl; |
---|
640 | cout << "B: " << B << endl << endl; |
---|
641 | |
---|
642 | int i; |
---|
643 | CanonicalForm qa; |
---|
644 | CanonicalForm lco; |
---|
645 | |
---|
646 | for(i = 1; i <= 1; i++){ |
---|
647 | pol = newGCD(A * i, B * i); |
---|
648 | lco = pol.lc(); |
---|
649 | cout << "new GCD: " << (pol / lco) << endl; |
---|
650 | } |
---|
651 | |
---|
652 | for(i = 1; i<=1; i++){ |
---|
653 | pol = gcd(A*i, B*i); |
---|
654 | lco = pol.lc(); |
---|
655 | cout << "old GCD: " << pol / lco << endl; |
---|
656 | } |
---|
657 | |
---|
658 | } |
---|
659 | #endif |
---|
660 | |
---|