1 // Copyright (c) 1997-2006 John Abbott
2
3 // This file is part of the source of CoCoALib, the CoCoA Library.
4
5 // CoCoALib is free software: you can redistribute it and/or modify
6 // it under the terms of the GNU General Public License as published by
7 // the Free Software Foundation, either version 3 of the License, or
8 // (at your option) any later version.
9
10 // CoCoALib is distributed in the hope that it will be useful,
11 // but WITHOUT ANY WARRANTY; without even the implied warranty of
12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 // GNU General Public License for more details.
14
15 // You should have received a copy of the GNU General Public License
16 // along with CoCoALib. If not, see <http://www.gnu.org/licenses/>.
17
18
19 #include "DUPZfactor_combine.h"
20 #include "mpz_log.h"
21 /* We use a random number in DUPZfactor_combine_init, so include stdlib.h */
22 #include <stdlib.h>
23 #include <stddef.h>
24 /* Use float.h just to get DBL_EPSILON and DBL_MANT_DIG, used in DUPZfactor_combiner_set_n1coeffd_bounds */
25 #include <float.h>
26
27 #include <math.h>
28 #include "jalloc.h"
29 #include "logi.h"
30 #include "WGD.h"
31 #include "DUPZevaluate.h"
32 #include "DUPZfactor_bound.h"
33 #include "DUPFF.h"
34 #include "DUPZ_DUPFF.h"
35 #include "DUPZfactor_refine_fds.h"
36
37 /* This is defined in DUPZfactor.c. */
38 extern void DUPZfactor_finish(DUPZfactor_info THIS);
39
40 /***************************************************************************/
41 /* These static variables are initialised by DUPZfactor_combine_init. */
42
43 static mpz_t num, den, tmp, lcm;
44
45 /* for accounting/debugging */
46 static int divisions;
47 static int count1, count2, count3, count4, count5, count6;
48
49 /***************************************************************************/
50 /* This function sets values for P and Q the numerator and denominator */
51 /* bounds to be used during rational reconstruction. During early */
52 /* searching the value for Q is "guessed" (approx 1/(2*sqrt(lc(f)))). */
53
DUPZfactor_combiner_set_PQ(DUPZfactor_combiner THIS)54 static void DUPZfactor_combiner_set_PQ(DUPZfactor_combiner THIS)
55 {
56 /* numerator/denominator bounds for rational reconstruction */
57 /* Updated whenever non-monic factors are found */
58 if (mpz_log(THIS->modulus) > logi(2) + 2*mpz_log(DUPZlc(THIS->info->f)))
59 mpz_abs(THIS->Q, DUPZlc(THIS->info->f));
60 else
61 {
62 mpz_set_ui(THIS->Q, 1);
63 mpz_mul_2exp(THIS->Q, THIS->Q, mpz_sizeinbase(THIS->modulus, 4)-1);
64 }
65
66 mpz_fdiv_q(THIS->P, THIS->modulus, THIS->Q);
67 mpz_fdiv_q_ui(THIS->P, THIS->P, 2);
68 }
69
70 #if 0
71 static void DUPZfactor_combiner_set_n1coeff_bounds(DUPZfactor_combiner THIS)
72 {
73 mpz_t tmp;
74 int i;
75 double log_bound;
76
77 mpz_init(tmp);
78 for (i=0; i < nfactors; i++)
79 {
80 fi = *THIS->factors[i];
81 mpz_mul(tmp, fi->coeff[degfi - delta], DUPZlc(THIS->info->f));
82 mpz_fdiv_r(tmp, tmp, modulus);
83 mpz_mul_2exp(tmp, tmp, 32); /* assumes unsigned is 32 bits */
84 mpz_fdiv_q(tmp, tmp, modulus);
85 n1coeff[i] = mpz_get_ui(tmp);
86 }
87
88 log_bound = THIS->info->bounds->root_bound + THIS->info->bounds->leading_coeff - THIS->height * logi(THIS->p);
89 if (log_bound < -32*logi(2)) n1bound = 1;
90 else n1bound = 1 + exp(32*logi(2) + log_bound);
91 mpz_clear(tmp);
92 }
93 #endif
94
DUPZfactor_combiner_set_n1coeffd_bounds(DUPZfactor_combiner THIS)95 static void DUPZfactor_combiner_set_n1coeffd_bounds(DUPZfactor_combiner THIS)
96 {
97 /* const int nfactors = DUPFFlist_length(THIS->info->pfactors);*/
98 int log2_modulus, shift;
99 double tmpd, log_modulus;
100 const int float_bits = DBL_MANT_DIG;
101 const double eps = DBL_EPSILON;
102 DUPZ fi;
103 mpz_t n1coeff, half_modulus; /* workspace */
104 int i, delta;
105
106 log_modulus = mpz_log(THIS->modulus);
107 mpz_init(n1coeff);
108 mpz_init_set(half_modulus, THIS->modulus);
109 mpz_fdiv_q_ui(half_modulus, half_modulus, 2);
110 log2_modulus = mpz_sizeinbase(THIS->modulus, 2);
111 shift = log2_modulus - float_bits - 10; /* 10 is rather arbitrary */
112 if (shift < 0) shift = 0;
113 delta=1;
114 for (i=0; i < THIS->nfactors_orig; i++)
115 {
116 if (THIS->used[i]) continue;
117 fi = *THIS->factors[i];
118 mpz_set(n1coeff, fi->coeffs[DUPZdeg(fi)-delta]);
119 /*printf("d-1 coeff of f[%d]=",i);mpz_out_str(stdout,10,n1coeff);printf("\n");*/
120 mpz_mul(n1coeff, n1coeff, DUPZlc(THIS->info->f));
121 mpz_fdiv_r(n1coeff, n1coeff, THIS->modulus);
122 if (mpz_cmp(n1coeff, half_modulus) > 0)
123 mpz_sub(n1coeff, n1coeff, THIS->modulus);
124 mpz_tdiv_q_2exp(n1coeff, n1coeff, shift);
125 THIS->n1coeffd[i] = mpz_get_d(n1coeff);
126 mpz_tdiv_q_2exp(n1coeff, THIS->modulus, shift);
127 THIS->n1coeffd[i] /= mpz_get_d(n1coeff);
128 /*DUPZprint(fi);*/
129 /*printf("THIS->n1coeffd[%d]=%g\n",i,THIS->n1coeffd[i]);*/
130 }
131 mpz_clear(n1coeff);
132 mpz_clear(half_modulus);
133 THIS->n1coeffd_bound = exp(THIS->info->bounds->leading_coeff + THIS->info->bounds->root_bound - mpz_log(THIS->modulus));
134 THIS->n1coeffd_bound += 2*eps;
135 THIS->n1coeffd_sum[0] = 0;
136 /*
137 *printf("mpz_log(THIS->modulus)=%g\n",mpz_log(THIS->modulus));
138 *printf("THIS->info->bounds->root_bound=%g\n",THIS->info->bounds->root_bound);
139 *printf("THIS->info->bounds->leading_coeff=%g\n",THIS->info->bounds->leading_coeff);
140 *printf("THIS->n1coeffd_bound=%g\n",THIS->n1coeffd_bound);
141 */
142 THIS->n1coeffd_flag = 0;
143 tmpd = THIS->info->bounds->leading_coeff + THIS->info->bounds->root_bound;
144 if (log_modulus < logi(10) + logi(DUPZdeg(THIS->info->f)) + tmpd) return;
145 THIS->n1coeffd_flag = 1;
146 }
147
148
149
DUPZfactor_combiner_ctor(DUPZfactor_info info)150 DUPZfactor_combiner DUPZfactor_combiner_ctor(DUPZfactor_info info)
151 {
152 DUPZ f = info->f; /* aliasing for convenience */
153 DUPZfactor_combiner THIS;
154 int i;
155 double rt_bd;
156 const int nfactors = DUPFFlist_length(info->pfactors);
157 const int range = 10; /* should be the same as range in DUPZfactor_sqfr */
158
159 THIS = (DUPZfactor_combiner)MALLOC(sizeof(struct DUPZfactor_combine_struct));
160 THIS->info = info;
161 THIS->tuple_deg = 0;
162
163 mpz_init_set(THIS->modulus, info->Q);
164 THIS->nfactors_orig = nfactors; /* this never changes */
165 THIS->nfactors = nfactors; /* this counts down as factors are used */
166 THIS->factors = (DUPZ**)MALLOC(nfactors*sizeof(DUPZ*));
167 DUPZfactor_lift_output(THIS->factors, info->lifter, 0);
168 /*
169 *printf("There are %d modular factors:\n", nfactors);
170 *for(i=0;i<nfactors;i++)DUPZprint(*THIS->factors[i]);
171 */
172 THIS->combination = (int*)MALLOC(nfactors*sizeof(int));
173 THIS->used = (int*)MALLOC(nfactors*sizeof(int));
174 for (i=0; i < nfactors; i++) THIS->used[i] = 0;
175
176 THIS->n1coeffd = (double*)MALLOC(nfactors*sizeof(double));
177 THIS->n1coeffd_sum = (double*)MALLOC(nfactors*sizeof(double));
178 DUPZfactor_combiner_set_n1coeffd_bounds(THIS);
179
180 /* pick random value from [-range, range] avoiding 0 */
181 THIS->evaluation_point = rand()%(2*range) - range;
182 if (THIS->evaluation_point == 0) THIS->evaluation_point = range;
183
184 /* the two lines below are odd to avoid risk of overflow */
185 rt_bd = info->bounds->root_bound;
186 THIS->log_ev_pt = rt_bd + log(1+exp(logi(abs(THIS->evaluation_point)) - rt_bd));
187 mpz_init(THIS->f_value);
188 DUPZevaluate(THIS->f_value, f, THIS->evaluation_point);
189 /* we have removed small linear factors from f, so f_value is non-zero */
190 THIS->factor_value = (mpz_t*)MALLOC(nfactors*sizeof(mpz_t));
191 for (i=0; i < nfactors; i++)
192 {
193 mpz_init(THIS->factor_value[i]);
194 DUPZevaluate(THIS->factor_value[i], *THIS->factors[i], THIS->evaluation_point);
195 }
196
197 mpz_init(THIS->P);
198 mpz_init(THIS->Q);
199 DUPZfactor_combiner_set_PQ(THIS);
200 THIS->factor = DUPZnew(DUPZdeg(f)); /* Must be big enough to contain either the */
201 THIS->quot = DUPZnew(DUPZdeg(f)); /* putative factor or its complement. */
202 THIS->rem = DUPZnew(DUPZdeg(f));
203 mpz_init(THIS->value_at_1);
204 DUPZevaluate(THIS->value_at_1, f, 1);
205 mpz_init(THIS->value_at_m1);
206 DUPZevaluate(THIS->value_at_m1, f, -1);
207
208 THIS->tc_combination = (int*)MALLOC(nfactors*sizeof(int));
209 THIS->tc_combination[0] = -1;
210 THIS->tc_product = (mpz_t*)MALLOC(nfactors*sizeof(mpz_t));
211 for (i=0; i < nfactors; i++) mpz_init(THIS->tc_product[i]);
212
213 return THIS;
214 }
215
216
DUPZfactor_combiner_dtor(DUPZfactor_combiner THIS)217 void DUPZfactor_combiner_dtor(DUPZfactor_combiner THIS)
218 {
219 int i;
220
221 DUPZfree(THIS->factor);
222 DUPZfree(THIS->quot);
223 DUPZfree(THIS->rem);
224 /*****for (i=0; i < THIS->nfactors_orig; i++) DUPZfree(THIS->factors[i]);*****/
225 FREE(THIS->factors);
226 mpz_clear(THIS->modulus);
227 mpz_clear(THIS->f_value);
228 mpz_clear(THIS->value_at_1);
229 mpz_clear(THIS->value_at_m1);
230 mpz_clear(THIS->P);
231 mpz_clear(THIS->Q);
232 for(i=0; i < THIS->nfactors_orig; i++) mpz_clear(THIS->factor_value[i]);
233 FREE(THIS->factor_value);
234 FREE(THIS->used);
235 FREE(THIS->n1coeffd);
236 FREE(THIS->n1coeffd_sum);
237
238 FREE(THIS->combination);
239 FREE(THIS->tc_combination);
240 for (i=0; i < THIS->nfactors_orig; i++) mpz_clear(THIS->tc_product[i]);
241 FREE(THIS->tc_product);
242
243 FREE(THIS);
244 }
245
246
247
248 /***************************************************************************/
249 /* These are not constructor/destructor; "_init" must be called before */
250 /* searching combinations, and "_done" must be called after all searching */
251 /* has been finished. These functions deal with globals. */
252
DUPZfactor_combine_init()253 void DUPZfactor_combine_init()
254 {
255 divisions=count1=count2=count3=count4=count5=count6=0;
256 mpz_init(num);
257 mpz_init(den);
258 mpz_init(tmp);
259 mpz_init(lcm);
260 }
261
DUPZfactor_combine_done()262 void DUPZfactor_combine_done()
263 {
264 mpz_clear(num); mpz_clear(den); mpz_clear(tmp); mpz_clear(lcm);
265 #ifdef FACTOR_DEBUG
266 printf("\nSearch complete. count1=%d, count2=%d, count3=%d, count4=%d, count5=%d, count6=%d\n", count1, count2, count3, count4, count5, count6);
267 printf("Total number of divisions: %d\n", divisions);
268 #endif
269 }
270
271
272 /***************************************************************************/
273
274 static int DUPZfactor_combine1(DUPZfactor_combiner THIS, int n, int r, int i);
275
DUPZfactor_combine(DUPZfactor_combiner THIS)276 static void DUPZfactor_combine(DUPZfactor_combiner THIS)
277 {
278 int tuple_size, max_tuple_size;
279 double log_modulus, log_lc, factor_bound;
280
281 THIS->single_factor = 0;
282 log_modulus = mpz_log(THIS->modulus);
283 if (THIS->early)
284 {
285 max_tuple_size = (int)(10/logi(THIS->nfactors));
286 if (max_tuple_size > THIS->nfactors/2) max_tuple_size = THIS->nfactors/2;
287 }
288 else
289 {
290 log_lc = THIS->info->bounds->leading_coeff;
291 factor_bound = THIS->info->bounds->all_factors;
292 /* Search through only half the combinations if modulus is big enough */
293 if (log_modulus > 2*log_lc + factor_bound + 0.7) max_tuple_size = THIS->nfactors/2;
294 else { THIS->single_factor = 1; max_tuple_size = THIS->nfactors - 1; }
295 }
296 for (tuple_size=1; tuple_size <= max_tuple_size; tuple_size++)
297 {
298 if (tuple_size >= THIS->nfactors) break; /* only relevant during early searching */
299 #ifdef FACTOR_DEBUG
300 printf("\nInterim. count1=%d, count2=%d, count3=%d, count4=%d, count5=%d, count6=%d", count1, count2, count3, count4, count5, count6);
301 printf("\n%d-tuples: ", tuple_size);
302 #endif
303 THIS->tuple_size = tuple_size;
304 THIS->tc_combination[tuple_size-1] = -1; /* so first call to tctest works properly */
305 DUPZfactor_combine1(THIS, THIS->nfactors, 0, 0);
306 if (THIS->early) continue;
307 /* Check if we can reduce to a search through half the combinations. */
308 log_lc = THIS->info->bounds->leading_coeff;
309 factor_bound = THIS->info->bounds->all_factors;
310 THIS->single_factor = 0;
311 if (log_modulus > 2*log_lc + factor_bound + 0.7) max_tuple_size = THIS->nfactors/2;
312 else { THIS->single_factor = 1; max_tuple_size = THIS->nfactors - 1; }
313 }
314 if (THIS-> early || DUPZdeg(THIS->info->f) <= 0) return;
315 DUPZfactors_add(THIS->info->irreds, THIS->info->f);
316 THIS->info->f->deg = 0;
317 }
318
319
DUPZfactor_combine_early(DUPZfactor_combiner THIS)320 void DUPZfactor_combine_early(DUPZfactor_combiner THIS)
321 {
322 THIS->early = 1; /* flag that we're doing early searching */
323 DUPZfactor_combine(THIS);
324 }
325
326
DUPZfactor_combine_final(DUPZfactor_combiner THIS)327 void DUPZfactor_combine_final(DUPZfactor_combiner THIS)
328 {
329 THIS->early = 0; /* flag that this is the final search */
330 DUPZfactor_combine(THIS);
331 }
332
333
334
335 /***************************************************************************/
336 /***************************************************************************/
337 /* n-1 coeff test */
338
DUPZfactor_combine_n1test(DUPZfactor_combiner THIS)339 static int DUPZfactor_combine_n1test(DUPZfactor_combiner THIS)
340 {
341 int i, deg, df, D;
342 const DUPZ f = THIS->info->f;
343 DUPZ fi; /* just a convenient alias */
344 double n1_coeff_bound;
345
346 deg = THIS->tuple_deg;
347
348 df = DUPZdeg(f);
349
350 /* Compute in tmp the n-1 coeff of the current tuple (or its complement). */
351 /* Cannot just subtract from n-1 coeff in f because that value isn't to hand. */
352 mpz_set_ui(tmp, 0);
353 if (THIS->single_factor || deg <= df/2) /* tuple has low degree, so compute n-1 coeff directly */
354 {
355 D = deg;
356 for (i=0; i < THIS->tuple_size; i++)
357 {
358 fi = *THIS->factors[THIS->combination[i]];
359 mpz_add(tmp, tmp, fi->coeffs[DUPZdeg(fi)-1]);
360 }
361 }
362 else /* tuple has high degree, so compute n-1 coeff of its complement */
363 {
364 D = df - deg;
365 for (i=0; i < THIS->nfactors_orig; i++)
366 {
367 if (THIS->used[i]) continue;
368 fi = *THIS->factors[i];
369 mpz_add(tmp, tmp, fi->coeffs[DUPZdeg(fi)-1]);
370 }
371 }
372 mpz_fdiv_r(tmp, tmp, THIS->modulus);
373 /* Now compute the rational corresponding to the modular value. */
374 modular_to_rational(num, den, THIS->Q, tmp, THIS->modulus);
375
376 /* check magnitude of numerator */
377 if (mpz_sgn(num) == 0) return 1;
378 n1_coeff_bound = THIS->info->bounds->root_bound + logi(D) + 0.001;
379 if (mpz_log(num) - mpz_log(den) > n1_coeff_bound) return 0;
380
381 /* check denominator divides leading coeff */
382 mpz_fdiv_r(num, DUPZlc(f), den);
383 if (mpz_sgn(num) != 0) return 0;
384 return 1;
385 }
386
387
388 /***************************************************************************/
389 /* Trailing coeff test */
390 /* SIDE EFFECT: sets THIS->log_tc (used in coefftest below) */
391
DUPZfactor_combine_tctest(DUPZfactor_combiner THIS)392 static int DUPZfactor_combine_tctest(DUPZfactor_combiner THIS)
393 {
394 DUPZ f = THIS->info->f; /* convenient alias */
395 int i, deg, df;
396 DUPZ fi; /* just a convenient alias */
397
398 deg = THIS->tuple_deg;
399 df = DUPZdeg(f);
400 /* compute in tmp the trailing coeff of the current combination */
401 mpz_set_ui(tmp, 1);
402 if (THIS->single_factor || deg <= df/2)
403 {
404 for (i=0; i < THIS->tuple_size; i++)
405 {
406 if (THIS->combination[i] == THIS->tc_combination[i]) continue;
407 THIS->tc_combination[i] = THIS->combination[i];
408 THIS->tc_combination[i+1] = -1;
409 fi = *THIS->factors[THIS->combination[i]];
410 if (i == 0) mpz_set(THIS->tc_product[i], fi->coeffs[0]);
411 else
412 {
413 mpz_mul(THIS->tc_product[i], THIS->tc_product[i-1], fi->coeffs[0]);
414 mpz_fdiv_r(THIS->tc_product[i], THIS->tc_product[i], THIS->modulus);
415 }
416 }
417 mpz_set(tmp, THIS->tc_product[i-1]);
418 }
419 else
420 {
421 for (i=0; i < THIS->nfactors_orig; i++)
422 {
423 if (THIS->used[i]) continue;
424 fi = *THIS->factors[i];
425 mpz_mul(tmp, tmp, fi->coeffs[0]);
426 mpz_fdiv_r(tmp, tmp, THIS->modulus);
427 }
428 }
429
430 /* NB tmp != 0 because we excluded primes dividing the trailing coeff */
431 modular_to_rational(num, den, THIS->Q, tmp, THIS->modulus);
432
433 /* check numerator divides trailing coeff of original polynomial */
434 mpz_fdiv_r(tmp, f->coeffs[0], num);
435 if (mpz_sgn(tmp) != 0) return 0;
436
437 /* check denominator divides leading coeff */
438 mpz_fdiv_r(tmp, DUPZlc(f), den);
439 if (mpz_sgn(tmp) != 0) return 0;
440
441 THIS->log_tc = mpz_log(num) - mpz_log(den); /* used in coefftest below */
442 return 1;
443 }
444
445
DUPZfactor_combine_evaltest(DUPZfactor_combiner THIS)446 static int DUPZfactor_combine_evaltest(DUPZfactor_combiner THIS)
447 {
448 int i, deg, df, D;
449 DUPZ f = THIS->info->f; /* convenient alias */
450
451 df = DUPZdeg(f);
452 deg = THIS->tuple_deg;
453 if (THIS->single_factor || deg <= df/2) D = deg; else D = df - deg;
454 /* Do eval at r test only if P/Q > (r + root_bound)^D */
455 if (mpz_log(THIS->P) - mpz_log(THIS->Q) < D*THIS->log_ev_pt + 0.001) return 1;
456 /* compute in tmp the value at r of the current combination (or its complement) */
457 mpz_set_ui(tmp, 1);
458 if (THIS->single_factor || deg <= df/2) /* tuple has low degree, so compute its value directly */
459 {
460 for (i=0; i < THIS->tuple_size; i++)
461 {
462 mpz_mul(tmp, tmp, THIS->factor_value[THIS->combination[i]]);
463 mpz_fdiv_r(tmp, tmp, THIS->modulus);
464 }
465 }
466 else /* tuple has high degree, so compute value of its complement instead */
467 {
468 for (i=0; i < THIS->nfactors_orig; i++)
469 {
470 if (THIS->used[i]) continue;
471 mpz_mul(tmp, tmp, THIS->factor_value[i]);
472 mpz_fdiv_r(tmp, tmp, THIS->modulus);
473 }
474 }
475
476 modular_to_rational(num, den, THIS->Q, tmp, THIS->modulus);
477
478 if (mpz_sgn(num) == 0) return 0;
479
480 /* check numerator divides the value of f */
481 mpz_fdiv_r(num, THIS->f_value, num);
482 if (mpz_sgn(num) != 0) return 0;
483
484 /* check denominator divides lc(f) */
485 mpz_fdiv_r(den, DUPZlc(f), den);
486 if (mpz_sgn(den) != 0) return 0;
487
488 return 1;
489 }
490
491
492 /***************************************************************************/
493 /* Coefficient test: the modular factors are multiplied, and the coeffs of */
494 /* the product converted to rationals checking numerator and denominator */
495 /* have possible values. This test is relatively expensive. */
496
DUPZfactor_combine_coefftest(DUPZfactor_combiner THIS)497 static int DUPZfactor_combine_coefftest(DUPZfactor_combiner THIS)
498 {
499 int i, j, deg, df, D;
500 DUPZ f = THIS->info->f; /* alias */
501 DUPZ g = THIS->factor; /* alias */
502 double log_coeff, forward, reverse; /* forward and reverse bounds */
503
504 deg = THIS->tuple_deg;
505 df = DUPZdeg(f);
506 /* if deg <= df/2 then we multiply out the combination, */
507 /* otherwise we multiply out the complementary combination. */
508 if (THIS->single_factor || deg <= df/2)
509 {
510 D = deg;
511 DUPZcopy2(g, *THIS->factors[THIS->combination[0]]);
512 for (i=1; i < THIS->tuple_size; i++)
513 {
514 DUPZmul3(g, g, *THIS->factors[THIS->combination[i]]);
515 for (j=DUPZdeg(g); j >= 0; j--)
516 mpz_fdiv_r(g->coeffs[j], g->coeffs[j], THIS->modulus);
517 }
518 }
519 else
520 {
521 D = df - deg;
522 mpz_set_ui(g->coeffs[0], 1);
523 g->deg = 0;
524 for (i=0; i < THIS->nfactors_orig; i++)
525 {
526 if (THIS->used[i]) continue;
527 DUPZmul3(g, g, *THIS->factors[i]);
528 for (j=DUPZdeg(g); j >= 0; j--)
529 mpz_fdiv_r(g->coeffs[j], g->coeffs[j], THIS->modulus);
530 }
531 }
532
533
534 forward = 0.001; /* add 0.001 to allow for rounding error */
535 reverse = THIS->log_tc + D*THIS->info->bounds->reverse_root_bound + 0.001;
536 mpz_set_ui(lcm, 1);
537 for (i=D-1; i >= 0; i--)
538 {
539 forward += logi(i+1) - logi(D-i) + THIS->info->bounds->root_bound;
540 reverse += logi(i+1) - logi(D-i) - THIS->info->bounds->reverse_root_bound;
541 if (mpz_sgn(g->coeffs[i]) == 0) continue;
542 modular_to_rational(num, den, THIS->Q, g->coeffs[i], THIS->modulus);
543 mpz_fdiv_r(tmp, DUPZlc(f), den); /* check denom divides lc(f) */
544 if (mpz_sgn(tmp) != 0) return 0; /* if not, factor is false */
545
546 /* check coeff is within permitted range */
547 log_coeff = mpz_log(num) - mpz_log(den);
548 if (log_coeff > forward || log_coeff > reverse) return 0;
549
550 /* store rational in THIS->factor and THIS->quot */
551 mpz_set(g->coeffs[i], num);
552 mpz_set(THIS->quot->coeffs[i], den);
553
554 /* update the common denominator */
555 mpz_gcd(tmp, lcm, den);
556 mpz_divexact(den, den, tmp);
557 mpz_mul(lcm, lcm, den);
558 }
559
560 /* make THIS->factor have a common denominator */
561 for (i=D-1; i >= 0; i--)
562 {
563 if (mpz_sgn(g->coeffs[i]) == 0) continue;
564 if (mpz_cmp(lcm, THIS->quot->coeffs[i]) == 0) continue;
565 mpz_divexact(den, lcm, THIS->quot->coeffs[i]);
566 mpz_mul(g->coeffs[i], g->coeffs[i], den);
567 }
568 mpz_set(g->coeffs[D], lcm);
569
570 return 1;
571 }
572
573
574 /***************************************************************************/
575 /* Division test: we divide the putative factor into the polynomial. */
576 /* We initially check that the values at 1 and -1 divide those of f. */
577 /* Then we just do long division. This test is potentially very expensive */
578 /* if the putative factor is not a true one -- should hardly ever happen. */
579
DUPZfactor_combine_dividetest(DUPZfactor_combiner THIS)580 static int DUPZfactor_combine_dividetest(DUPZfactor_combiner THIS)
581 {
582 DUPZ f = THIS->info->f; /* convenient alias */
583
584 /* These two preliminary checks are most relevant during early detection */
585 /* check value of factor at 1 divides the value of f at 1 */
586 DUPZevaluate(tmp, THIS->factor, 1);
587 if (mpz_sgn(tmp) == 0) return 0;
588 mpz_fdiv_r(tmp, THIS->value_at_1, tmp);
589 if (mpz_sgn(tmp) != 0) return 0;
590
591 /* check value of factor at -1 divides the value of f at -1 */
592 DUPZevaluate(tmp, THIS->factor, -1);
593 if (mpz_sgn(tmp) == 0) return 0;
594 mpz_fdiv_r(tmp, THIS->value_at_m1, tmp);
595 if (mpz_sgn(tmp) != 0) return 0;
596
597 /* we expect this division almost always to be exact */
598 DUPZdiv4(THIS->quot, THIS->rem, f, THIS->factor);
599 divisions++;
600 if (DUPZdeg(THIS->rem) >= 0) return 0;
601
602 /* Swap factor and quot if we have built the complement in factor. */
603 if (!THIS->single_factor && THIS->tuple_deg > DUPZdeg(THIS->info->f)/2)
604 DUPZswap(THIS->quot, THIS->factor);
605 return 1;
606 }
607
608
609 /***************************************************************************/
610 /* This function is called if we cannot immediately show that the factor */
611 /* is irreducible. A more stringent irreducibility test is conducted, and */
612 /* if necessary, a subsidiary factorization is performed. */
613
DUPZfactor_combine_fork(DUPZfactor_combiner old,DUPZ factor)614 static void DUPZfactor_combine_fork(DUPZfactor_combiner old, DUPZ factor)
615 {
616 int i, irred, dmax;
617 DUPZfactor_info fork;
618 DUPFF fp, rem;
619 DUPFFlist iter;
620 bound_info bd;
621 double lift_height = mpz_log(old->modulus) - 0.001; /* allow for rounding */
622
623 /* initial quick tests for "obvious" irreducibility */
624 if (old->tuple_size == 1 || !old->early)
625 {
626 DUPZfactors_add(old->info->irreds, factor);
627 return;
628 }
629 dmax = DUPZdeg(factor)/2;
630 while (!old->info->fds[dmax]) dmax--;
631 bd = DUPZcopy_refine_bound(old->info->bounds, factor, dmax);
632 if (bd->lift_bound < lift_height)
633 {
634 FREE(bd);
635 DUPZfactors_add(old->info->irreds, factor);
636 return;
637 }
638 fork = (DUPZfactor_info)MALLOC(sizeof(struct DUPZfactor_info_struct));
639 fork->f = DUPZcopy(factor);
640 fork->irreds = old->info->irreds;
641 fork->p = old->info->p;
642
643 fork->nprimes = old->info->nprimes;
644 fork->FFq = (FF*)MALLOC(fork->nprimes*sizeof(FF));
645 fork->qfactors = (DUPFFlist*)MALLOC(fork->nprimes*sizeof(DUPFFlist));
646 fork->fds = (int*)MALLOC((1+DUPZdeg(fork->f))*sizeof(int));
647 for (i=0; i <= DUPZdeg(fork->f); i++) fork->fds[i] = old->info->fds[i];
648 fork->bounds = NULL;
649 mpz_init_set(fork->recip_lcf, old->info->recip_lcf);
650 mpz_mul(fork->recip_lcf, fork->recip_lcf, DUPZlc(old->info->f));
651 mpz_divexact(fork->recip_lcf, fork->recip_lcf, DUPZlc(factor));
652 mpz_init_set(fork->Q, old->info->Q);
653 fork->lifter = NULL;
654
655 /* determine new factorizations modulo the various primes */
656 fp = DUPFFnew(DUPZdeg(fork->f));
657 rem = DUPFFnew(DUPZdeg(fork->f));
658 irred = 0;
659 for (i=0; i < old->info->nprimes; i++)
660 {
661 FFselect(old->info->FFq[i]);
662 fork->FFq[i] = FFctor(CurrentFF.prime); /* wasteful copy */
663 fork->qfactors[i] = NULL;
664 DUPZ_to_DUPFF2(fp, fork->f);
665 for (iter = old->info->qfactors[i]; iter; iter = iter->next)
666 {
667 DUPFFcopy2(rem, fp);
668 DUPFFrem2(rem, iter->poly);
669 if (DUPFFdeg(rem) >= 0) continue;
670 fork->qfactors[i] = DUPFFlist_append(DUPFFlist_ctor(DUPFFcopy(iter->poly), 1), fork->qfactors[i]);
671 }
672 irred = DUPZfactor_refine_fds(fork->fds, fork->qfactors[i]);
673 /* if (irred) break;*/ /* including this line makes deallocation trickier */
674 }
675 /* select correct prime for rest of computation */
676 fork->p_index = old->info->p_index;
677 fork->pfactors = fork->qfactors[fork->p_index];
678 FFselect(fork->FFq[fork->p_index]);
679 DUPFFfree(fp);
680 DUPFFfree(rem);
681 if (irred)
682 {
683 FREE(bd);
684 DUPZfactors_add(fork->irreds, fork->f);
685 DUPZfactor_info_dtor(fork);
686 return;
687 }
688
689 /* update bounds */
690 fork->dmax = old->info->dmax;
691 if (fork->dmax > DUPZdeg(fork->f)/2) fork->dmax = DUPZdeg(fork->f)/2;
692
693 while (!fork->fds[dmax]) dmax--;
694 fork->dmax = dmax;
695 /* fork->dmin = old->info->dmin; */
696 /* while (!fork->fds[fork->dmin]) fork->dmin++; */
697 DUPZrefine_bound(bd, fork->f, dmax);
698 fork->bounds = bd;
699 if (bd->lift_bound < lift_height)
700 {
701 DUPZfactors_add(fork->irreds, fork->f);
702 DUPZfactor_info_dtor(fork);
703 return;
704 }
705 DUPZfactor_info_set_target_height(fork);
706
707
708 /************************************************************* DUPZfactor_lift_revise1(fork, old->info->lifter->r); **********************/
709 #ifdef FACTOR_DEBUG
710 printf("\nFactor may be reducible: starting subsidiary factorization\n");
711 #endif
712 DUPZfactor_finish(fork);
713 #ifdef FACTOR_DEBUG
714 printf("Subsidiary factorization complete; resuming previous factorization\n");
715 #endif
716 DUPZfactor_info_dtor(fork);
717 }
718
719 /***************************************************************************/
720 /* Upon entry to this function THIS->factor contains the factor built from */
721 /* the current tuple, while THIS->quot contains its complement (i.e. the */
722 /* quotient of f by the factor found. */
723
DUPZfactor_combine_found(DUPZfactor_combiner THIS)724 static void DUPZfactor_combine_found(DUPZfactor_combiner THIS)
725 {
726 const DUPZ f = THIS->info->f; /* convenient alias */
727 DUPZ fi;
728 double rt_bd;
729 int i, irred, df;
730 DUPFF rem, fp;
731 DUPFFlist iter, next, *prev;
732
733 DUPZfactor_combine_fork(THIS, THIS->factor); /* recurse on factor */
734 for(i=0;i<THIS->tuple_size;i++){fi=*THIS->factors[THIS->combination[i]];fi->deg=-fi->deg;}
735 DUPZcopy2(f, THIS->quot);
736 mpz_mul(THIS->info->recip_lcf, THIS->info->recip_lcf, DUPZlc(THIS->factor));
737 df = DUPZdeg(f);
738 THIS->nfactors -= THIS->tuple_size;
739
740 /* ---- we may be able to prove that what's left is irreducible ----- */
741 /* determine new factorizations modulo the various primes */
742 fp = DUPFFnew(df);
743 rem = DUPFFnew(df);
744 irred = 0;
745 for (i=0; i < THIS->info->nprimes; i++)
746 {
747 FFselect(THIS->info->FFq[i]);
748 prev = &THIS->info->qfactors[i];
749 DUPZ_to_DUPFF2(fp, f);
750 for (iter = THIS->info->qfactors[i]; iter; iter = next)
751 {
752 next = iter->next;
753 DUPFFcopy2(rem, fp);
754 DUPFFrem2(rem, iter->poly);
755 if (DUPFFdeg(rem) >= 0) DUPFFlist_elem_dtor(iter);
756 else { *prev = iter; prev = &iter->next; }
757 }
758 *prev = NULL;
759 irred = DUPZfactor_refine_fds(THIS->info->fds, THIS->info->qfactors[i]);
760 if (irred) break;
761 }
762
763 THIS->info->pfactors = THIS->info->qfactors[THIS->info->p_index];
764 FFselect(THIS->info->FFq[THIS->info->p_index]);
765 DUPFFfree(fp);
766 DUPFFfree(rem);
767 if (irred)
768 {
769 DUPZfactors_add(THIS->info->irreds, f);
770 for (i=0; i < THIS->nfactors_orig; i++) THIS->used[i] = 1;
771 for (i=0; i < THIS->nfactors_orig; i++) {fi=*THIS->factors[i];if (DUPZdeg(fi)>0) fi->deg = -fi->deg; }
772 THIS->nfactors = 0;
773 f->deg = 0;
774 return;
775 }
776 if (THIS->info->dmax > df/2) THIS->info->dmax = df/2;
777 while (!THIS->info->fds[THIS->info->dmax]) THIS->info->dmax--;
778 /* while (!THIS->info->fds[THIS->info->dmin]) THIS->info->dmin++; */
779
780 /* update bounds */
781 DUPZrefine_bound(THIS->info->bounds, f, THIS->info->dmax);
782 DUPZfactor_info_set_target_height(THIS->info);
783
784 if (THIS->early && THIS->info->target_height <= THIS->info->current_height)
785 {
786 THIS->early = 0;
787 #ifdef FACTOR_DEBUG
788 printf("[early search has now become final search]");
789 #endif
790 }
791
792 /* update log_ev_pt */
793 rt_bd = THIS->info->bounds->root_bound;
794 THIS->log_ev_pt = rt_bd + log(1+exp(logi(abs(THIS->evaluation_point)) - rt_bd));
795
796 /* update f_value */
797 DUPZevaluate(THIS->f_value, f, THIS->evaluation_point);
798
799 /* update P, Q */
800 DUPZfactor_combiner_set_PQ(THIS);
801
802 /* update n1coeff_bounds */
803 DUPZfactor_combiner_set_n1coeffd_bounds(THIS);
804 }
805
806
807
DUPZfactor_combine1(DUPZfactor_combiner THIS,int n,int r,int i)808 static int DUPZfactor_combine1(DUPZfactor_combiner THIS, int n, int r, int i)
809 /* any true factors found will be added to ans. */
810 /* n is number of factors still to be considered, */
811 /* r is number of factors already picked, */
812 /* i is index of next factor to consider */
813 {
814 int deg, found;
815
816 if (n < THIS->tuple_size - r) return 0;
817 if (r < THIS->tuple_size)
818 {
819 /* if (deg > dmax) return;*/
820 if (THIS->used[i]) return DUPZfactor_combine1(THIS, n, r, i+1);
821 THIS->combination[r] = i;
822 THIS->used[i] = 2;
823 deg = DUPZdeg(*THIS->factors[i]);
824 THIS->tuple_deg += deg;
825 THIS->n1coeffd_sum[r+1] = THIS->n1coeffd_sum[r] + THIS->n1coeffd[i];
826 /* first try with this factor */
827 found = DUPZfactor_combine1(THIS, n-1, r+1, i+1);
828 THIS->tuple_deg -= deg;
829 if (!found) THIS->used[i] = 0;
830 /* special case if n = nfactors/2 */
831 if (!found && !THIS->single_factor && (r == 0) && (THIS->nfactors == 2*THIS->tuple_size)) return 0;
832 if (!found) return DUPZfactor_combine1(THIS, n-1, r, i+1);
833 if (r > 0) return 1; /* unwind to "top level" */
834 /* now we must be at "top level" */
835 if (!THIS->single_factor && THIS->nfactors < 2*THIS->tuple_size) return 0;
836 if (THIS->nfactors <= THIS->tuple_size) return 0;
837
838 /* now try without this factor */
839 return DUPZfactor_combine1(THIS, n - THIS->tuple_size, r, i+1);
840 }
841
842 /*
843 *printf("Trying tuple: \t"); for (i=0; i < THIS->tuple_size; i++) printf("%d ", THIS->combination[i]); printf("\n");
844 */
845 count1++;
846
847 /* check degree is possible */
848 if (THIS->info->fds[THIS->tuple_deg] == 0) return 0;
849 count2++;
850
851
852 if (THIS->n1coeffd_flag)
853 {
854 /* Compute in deg the nearest integer to THIS->n1coeffd_sum */
855 if (THIS->n1coeffd_sum[THIS->tuple_size] >= 0)
856 deg = (int)(THIS->n1coeffd_sum[THIS->tuple_size] + 0.5);
857 else deg = (int)(THIS->n1coeffd_sum[THIS->tuple_size] - 0.5);
858
859 if (fabs(THIS->n1coeffd_sum[THIS->tuple_size] - deg) >= THIS->tuple_deg * THIS->n1coeffd_bound) return 0;
860 }
861 else
862 if (!DUPZfactor_combine_n1test(THIS)) return 0;
863
864 /*
865 * if (!DUPZfactor_combine_1test(THIS)) {printf("!");fflush(stdout);return 0;}
866 */
867 count3++;
868
869 if (!DUPZfactor_combine_tctest(THIS)) return 0;
870 count4++;
871
872 if (!DUPZfactor_combine_evaltest(THIS)) return 0;
873 count5++;
874
875 if (!DUPZfactor_combine_coefftest(THIS)) return 0;
876 count6++;
877
878 if (!DUPZfactor_combine_dividetest(THIS)) return 0;
879
880 DUPZfactor_combine_found(THIS);
881
882 return 1;
883 }
884