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