Line data Source code
1 : /* Functions for sets of long ints, to factor (Z/NZ)* into a set of sums
2 : as described in section 5 of "Improved Stage 2 to $P\pm{}1$ Factoring
3 : Algorithms" by Peter L. Montgomery and Alexander Kruppa, ANTS 2008
4 : (8th Algorithmic Number Theory Symposium).
5 :
6 : Copyright 2007, 2008, 2009, 2012 Alexander Kruppa, Paul Zimmermann.
7 :
8 : This file is part of the ECM Library.
9 :
10 : The ECM Library is free software; you can redistribute it and/or modify
11 : it under the terms of the GNU Lesser General Public License as published by
12 : the Free Software Foundation; either version 3 of the License, or (at your
13 : option) any later version.
14 :
15 : The ECM Library is distributed in the hope that it will be useful, but
16 : WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
17 : or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
18 : License for more details.
19 :
20 : You should have received a copy of the GNU Lesser General Public License
21 : along with the ECM Library; see the file COPYING.LIB. If not, see
22 : http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
23 : 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
24 :
25 : #include "config.h"
26 : #include "ecm-impl.h"
27 : #include <stdlib.h>
28 : #ifdef HAVE_ALLOCA_H
29 : #include <alloca.h>
30 : #endif
31 : #ifdef TESTDRIVE
32 : #include <stdio.h>
33 : FILE *ECM_STDOUT, *ECM_STDERR;
34 : #endif
35 :
36 : /*****************************************************************
37 :
38 : Functions for processing sets
39 :
40 : A set is a cardinality of unsigned long type and an array of
41 : long ints.
42 : A set of sets is an unsigned long telling the number of sets,
43 : an array that has several sets stored back-to-back.
44 :
45 : *****************************************************************/
46 :
47 :
48 : /* Copy a set from "*S" to "*T". Assumes that the sets do not overlap,
49 : or that T < S. */
50 :
51 : static void
52 7572 : set_copy (set_long_t *T, set_long_t *S)
53 : {
54 : unsigned long i;
55 7572 : const unsigned long c = S->card; /* We might overwrite S->card */
56 :
57 7572 : T->card = c;
58 33397 : for (i = 0UL; i < c; i++)
59 25825 : T->elem[i] = S->elem[i];
60 7572 : }
61 :
62 :
63 : /* Exchange two adjacent sets in memory. Since all "elem" arrays are stored
64 : in the same chunk of allocated memory, and not in different chunks, we
65 : cannot simply swap the "elem" pointers.
66 : If the set T has size c and the next has size d, after the swap the set T
67 : will have size d and the next will have size c.
68 : */
69 : static void
70 1840 : set_swap (set_long_t *T)
71 : {
72 : set_long_t *next, *tmp;
73 :
74 1840 : next = sets_nextset (T);
75 1840 : tmp = alloca (set_sizeof (T->card));
76 : ASSERT(tmp != NULL);
77 1840 : set_copy (tmp, T);
78 1840 : set_copy (T, next);
79 : /* warning: sets_nextset(T) might differ from next, if T and next had
80 : different sizes */
81 1840 : set_copy (sets_nextset(T), tmp);
82 1840 : }
83 :
84 :
85 : /* Functions for sorting an array of longs */
86 :
87 : static inline void
88 1 : swap_long (long *a, long *b)
89 : {
90 : long t;
91 1 : t = *a;
92 1 : *a = *b;
93 1 : *b = t;
94 1 : }
95 :
96 : static inline void
97 596 : swapsort_long (long *a, long *b)
98 : {
99 596 : if (*a > *b)
100 1 : swap_long (a, b);
101 596 : }
102 :
103 : void
104 554 : quicksort_long (long *a, unsigned long l)
105 : {
106 : unsigned long i, j;
107 : long pivot;
108 :
109 554 : if (l < 2)
110 182 : return;
111 :
112 372 : j = l - 1;
113 372 : swapsort_long (a, a+j);
114 372 : if (l == 2)
115 260 : return;
116 :
117 112 : i = j / 2;
118 112 : swapsort_long (a, a+i);
119 112 : swapsort_long (a+i, a+j);
120 112 : if (l == 3)
121 75 : return;
122 :
123 37 : pivot = a[i]; /* Median of three */
124 :
125 : /* Stuff <= pivot goes in first list */
126 :
127 : /* Invariant: a[0 ... i-1] <= pivot, a[j+1 ... l-1] > pivot */
128 126 : for (i = 1; i < j;)
129 89 : if (a[i] > pivot)
130 : {
131 113 : for (; a[j] > pivot; j--);
132 36 : if (i < j)
133 0 : swap_long (a+(i++), a+j);
134 : }
135 : else
136 53 : i++;
137 :
138 : #ifdef WANT_ASSERT
139 : for (j = 0; j < i; j++)
140 : ASSERT (a[j] <= pivot);
141 : for (j = i; j < l; j++)
142 : ASSERT(a[j] > pivot);
143 : #endif
144 :
145 37 : quicksort_long (a, i);
146 37 : quicksort_long (a + i, l - i);
147 :
148 : #ifdef WANT_ASSERT
149 : for (j = 0; i < l - 1; i++)
150 : ASSERT (a[j] <= a[j + 1]);
151 : #endif
152 : }
153 :
154 :
155 : /* Returns max(S), where S == (Z/\beta Z)* as chosen by
156 : sets_get_factored_sorted() */
157 : /* Assumes that S == 0 at recursion entry */
158 : static void
159 4546523 : sets_max_recurse (mpz_t S, const unsigned long beta)
160 : {
161 4546523 : unsigned long P = beta, p, pk;
162 : unsigned int k;
163 :
164 4546523 : if (beta == 1UL)
165 1051016 : return;
166 :
167 3495507 : p = find_factor (P);
168 3495507 : k = 1; pk = p; P /= p;
169 4360850 : while (P % p == 0)
170 : {
171 865343 : k++;
172 865343 : pk *= p;
173 865343 : P /= p; /* P*pk=beta is invariant */
174 : }
175 3495507 : sets_max_recurse (S, P);
176 3495507 : mpz_mul_ui (S, S, pk);
177 3495507 : if (p == 2UL && k == 1)
178 0 : mpz_add_ui (S, S, P);
179 3495507 : else if (p == 2UL)
180 0 : mpz_add_ui (S, S, P * (pk / 2UL - 1UL));
181 3495507 : else if (p % 4UL == 1UL)
182 1352828 : mpz_add_ui (S, S, P * ((pk + p) / 2UL - 2UL));
183 2142679 : else if (p % 4UL == 3UL)
184 2142679 : mpz_add_ui (S, S, P * ((pk - 1UL) / 2UL));
185 : else
186 0 : abort();
187 : }
188 :
189 : void
190 1051016 : sets_max (mpz_t S, const unsigned long beta)
191 : {
192 1051016 : mpz_set_ui (S, 0UL);
193 1051016 : sets_max_recurse (S, beta);
194 1051016 : }
195 :
196 :
197 : /* Compute the set of sums over the "nr_sets" different sets in "*sets".
198 : The value of "add" is added to each element of the set of sums.
199 : "*sum" will have {\prod_{S \in "*sets"} #S} entries and must have
200 : enough memory allocated. This number of elements in the set of sums
201 : is the return value. In case of nr_sets == 0, "add" is written to *sets
202 : and 1 is returned. The sets in "*sets" are assumed to be non-empty.
203 : If "*sum" is NULL, nothing is written, but the return value is computed
204 : correctly. */
205 :
206 : static unsigned long
207 1461 : sets_sumset_recurse (long *sum, const set_long_t *sets,
208 : const unsigned long nr_sets, const long add)
209 : {
210 1461 : unsigned long i, j = 0UL;
211 :
212 1461 : if (nr_sets == 0UL)
213 : {
214 927 : if (sum != NULL)
215 927 : sum[0] = add;
216 927 : return 1UL;
217 : }
218 :
219 : ASSERT (sets->card > 0UL);
220 1515 : for (i = 0UL; i < sets->card; i++)
221 : {
222 : /* Test for overflow */
223 981 : ASSERT_ALWAYS (add <= 0 || add + sets->elem[i] > sets->elem[i]);
224 981 : ASSERT_ALWAYS (add >= 0 || add + sets->elem[i] < sets->elem[i]);
225 981 : j += sets_sumset_recurse (sum + j, sets_nextset(sets), nr_sets - 1UL,
226 981 : add + sets->elem[i]);
227 : }
228 :
229 534 : return j;
230 : }
231 :
232 :
233 : void
234 480 : sets_sumset (set_long_t *sum, const sets_long_t *sets)
235 : {
236 480 : sum->card = sets_sumset_recurse (sum->elem, sets->sets, sets->nr, 0L);
237 480 : }
238 :
239 :
240 : /* Returns the minimal (if minmax == -1) or maximal (minmax == 1) value
241 : in the set of sums over the sets in "*sets". */
242 :
243 : void
244 480 : sets_sumset_minmax (mpz_t sum, const sets_long_t *sets, const int minmax)
245 : {
246 : unsigned long i, nr;
247 480 : const set_long_t *set = sets->sets;
248 : long extremum;
249 :
250 : ASSERT (minmax == 1 || minmax == -1);
251 480 : mpz_set_ui (sum, 0UL);
252 :
253 3615 : for (nr = 0; nr < sets->nr; nr++)
254 : {
255 : ASSERT (set->card > 0UL);
256 3135 : extremum = set->elem[0];
257 :
258 8301 : for (i = 1UL; i < set->card; i++)
259 5166 : if ((minmax == -1 && set->elem[i] < extremum) ||
260 5166 : (minmax == 1 && set->elem[i] > extremum))
261 5166 : extremum = set->elem[i];
262 :
263 3135 : if (extremum >= 0)
264 3135 : mpz_add_ui (sum, sum, extremum);
265 : else
266 0 : mpz_sub_ui (sum, sum, -extremum);
267 3135 : set = sets_nextset (set);
268 : }
269 :
270 480 : return;
271 : }
272 :
273 :
274 : /* Store in (**L) arithmetic progressions of prime length whose sumset is
275 : k/2*R_n, an arithmetic progression centered at 0 of common difference k
276 : and cardinality n. If n is even, k must be as well to ensure integer
277 : results.
278 : I.e. n = 1: k/2*R_n = {0},
279 : n = 2: k/2*R_n = k/2 * {1, -1},
280 : n = 3: k/2*R_n = k * {-1, 0, 1},
281 : n = 4: k/2*R_n = k/2 * {-3, -1, 1, 3},
282 : n = 5: k/2*R_n = k * {-2, -1, 0, 1, 2} etc.
283 : _ADDS_ the size in bytes of the set to "*sets_size"
284 : */
285 :
286 : static unsigned long
287 5748 : sets_factored_Rn2 (set_long_t **L, size_t *sets_size, const long n,
288 : const long k)
289 : {
290 5748 : unsigned long nr = 0UL;
291 : long i, m, q, r;
292 5748 : size_t size = 0;
293 :
294 : /* n must be odd, or n and k both even */
295 5748 : ASSERT_ALWAYS(n % 2L == 1L || k % 2L == 0L);
296 : ASSERT(L != NULL);
297 5748 : m = k; /* The multiplier accumulated so far, init to k */
298 5748 : r = n; /* The remaining cofactor of n */
299 13512 : for (q = 2L; r > 1L; q = (q + 1L) | 1L) /* Find prime factors of n */
300 : {
301 : ASSERT (q <= r);
302 14034 : while (r % q == 0L)
303 : {
304 6270 : if (*L != NULL)
305 : {
306 : /* Add m*R_q/2 to list */
307 3135 : (*L)->card = q;
308 11436 : for (i = 0L; i < q; i++)
309 : {
310 8301 : const long t = m * (2L * i - q + 1L);
311 : ASSERT(t % 2L == 0L);
312 8301 : (*L)->elem[i] = t / 2L;
313 : }
314 3135 : *L = sets_nextset (*L);
315 3135 : nr++;
316 : }
317 6270 : size += set_sizeof ((unsigned long) q);
318 : /* Multiply this t to multiplier and treat remaining
319 : factors of the set */
320 6270 : m *= q;
321 6270 : r /= q;
322 : }
323 : }
324 5748 : if (sets_size != NULL)
325 5748 : *sets_size += size;
326 5748 : return nr;
327 : }
328 :
329 :
330 : /* Return a set L of sets M_i so that M_1 + ... + M_k is congruent to
331 : (Z/nZ)*, which is the set of residue classes coprime to n. The M_i all
332 : have prime cardinality.
333 : The size of the set of sets "*L" in bytes is computed and stored in
334 : "*sets_size" unless "*sets_size" is NULL.
335 : Return the number of sets in L.
336 : If L is the NULL pointer, nothing will be stored in L. The correct
337 : return value (number of set in L) and "*sets_size" value will still
338 : be computed, for example so that the correct amount of space can be
339 : allocated and factor_coprimeset() be called again.
340 : */
341 :
342 : static unsigned long
343 960 : sets_factor_coprime (sets_long_t *sets, size_t *sets_size,
344 : const unsigned long n)
345 : {
346 960 : unsigned long r, k, nr = 0UL;
347 : long p, np;
348 960 : size_t size = sizeof (unsigned long);
349 960 : set_long_t *set = NULL;
350 :
351 : ASSERT (n > 0UL);
352 960 : if (sets != NULL)
353 480 : set = sets->sets;
354 :
355 960 : r = n;
356 4040 : while (r > 1UL)
357 : {
358 18856 : for (p = 2L; r % p > 0L; p++); /* Find smallest prime p that divides r */
359 6740 : for (k = 0UL; r % p == 0UL; k++, r /= p); /* Find p^k || r */
360 3080 : np = n/p;
361 :
362 3080 : if (p == 2L && k == 1UL) /* Case 2^1. Deal with it before the */
363 : { /* while loop below decreases k. */
364 0 : if (set != NULL)
365 : {
366 0 : set->card = 1UL;
367 0 : set->elem[0] = np;
368 0 : set = sets_nextset (set);
369 : }
370 0 : size += set_sizeof (1UL);
371 0 : nr++;
372 : }
373 :
374 : /* If k > 1, do the \sum_{i=1}^{k-1} p^i (Z/pZ) part here.
375 : (Z/pZ) is represented by an arithmetic progression of
376 : common difference 1 and length p. */
377 :
378 3660 : while (k-- > 1UL)
379 : {
380 580 : nr += sets_factored_Rn2 (&set, &size, p, np);
381 580 : np /= p;
382 : }
383 :
384 3080 : if (p % 4L == 3L)
385 : {
386 : /* We can use \hat{S}_p. Factor as
387 : {-(p+1)/4, (p+1)/4} + C_{(p-1)/2} */
388 :
389 : /* Add the {-(p+1)/4, (p+1)/4} set to L */
390 2088 : nr += sets_factored_Rn2 (&set, &size, 2L, (p + 1L) / 2L * np);
391 :
392 : /* Add the np / 2 * R_{(p-1)/2} set to L */
393 2088 : nr += sets_factored_Rn2 (&set, &size, (p - 1L) / 2L, np);
394 : }
395 992 : else if (p % 4L == 1L)
396 : {
397 : /* Factor into arithmetic progressions of prime length.
398 : R_{p} = {-p+1, -p+3, ..., p-3, p+1}, i.e.
399 : R_2 = {-1, 1}, R_3 = {-2, 0, 2}, R_4 = {-3, -1, 1, 3}
400 : We have R_{sq} = R_q + q*R_s */
401 :
402 992 : nr += sets_factored_Rn2 (&set, &size, p - 1L, 2L * np);
403 : }
404 : }
405 :
406 960 : if (sets_size != NULL)
407 480 : *sets_size = size;
408 :
409 960 : if (sets != NULL)
410 480 : sets->nr = nr;
411 :
412 960 : return nr;
413 : }
414 :
415 :
416 : /* Sort the sets in F into order of ascending cardinality. Uses a simple
417 : Bubble sort. */
418 :
419 : static void
420 480 : sets_sort (sets_long_t *sets)
421 : {
422 : unsigned long i, nr_unsorted, highest_swap;
423 : set_long_t *set;
424 :
425 : /* The last sets->nr - nr_unsorted sets in "*sets" are known to be
426 : sorted and each one larger than any of the first nr_unsorted sets
427 : in "*sets". */
428 480 : nr_unsorted = sets->nr;
429 1517 : while (nr_unsorted > 1UL)
430 : {
431 1037 : outputf (OUTPUT_TRACE, "nr_unsorted = %lu. ", nr_unsorted);
432 1037 : sets_print (OUTPUT_TRACE, sets);
433 1037 : set = sets->sets;
434 1037 : highest_swap = 1UL;
435 6783 : for (i = 1UL; i < nr_unsorted; i++)
436 : {
437 5746 : if (set->card > sets_nextset(set)->card)
438 : {
439 1840 : outputf (OUTPUT_TRACE, "sets_sort: swapping %lu and %lu\n",
440 : i - 1, i);
441 1840 : set_swap (set);
442 1840 : highest_swap = i;
443 : }
444 5746 : set = sets_nextset (set);
445 : }
446 1037 : nr_unsorted = highest_swap;
447 : }
448 :
449 : #ifdef WANT_ASSERT
450 : set = sets->sets;
451 : for (i = 0UL; i + 1UL < sets->nr; i++)
452 : {
453 : ASSERT(set->card <= sets_nextset (set)->card);
454 : set = sets_nextset (set);
455 : }
456 : #endif
457 480 : }
458 :
459 : /* Print all the sets in "*sets", formatted as a sum of sets */
460 :
461 : void
462 1092 : sets_print (const int verbosity, sets_long_t *sets)
463 : {
464 : unsigned long i, j;
465 1092 : set_long_t *set = sets->sets;
466 :
467 9470 : for (i = 0UL; i < sets->nr; i++)
468 : {
469 8378 : if (i == 0UL)
470 1092 : outputf (verbosity, "{");
471 : else
472 7286 : outputf (verbosity, " + {");
473 :
474 : ASSERT(set->card > 0UL);
475 8378 : outputf (verbosity, "%ld", set->elem[0]);
476 22442 : for (j = 1UL; j < set->card; j++)
477 14064 : outputf (verbosity, ", %ld", set->elem[j]);
478 8378 : outputf (verbosity, "}");
479 8378 : set = sets_nextset (set);
480 : }
481 1092 : outputf (verbosity, "\n");
482 1092 : }
483 :
484 :
485 : /* Extract sets whose set of sums has cardinality "d". We expect that
486 : "d" divides the cardinality of the set of sums of "sets" and that
487 : the cardinalities of the sets in "sets" are all prime.
488 : The amount of memory in bytes needed to store the extracted sets in
489 : "*extracted" is stored at "*extr_size". The number of sets extracted
490 : is returned. (If d = p_1 * ... * p_k, the return value is k and
491 : "*extr_size" is set_sizeof(p_1) + ... + set_sizeof(p_k).)
492 : If "*extracted" is NULL, nothing is written and no sets are removed
493 : from "*sets", but "*extr_size" is computed as if they were. */
494 :
495 : void
496 960 : sets_extract (sets_long_t *extracted, size_t *extr_size, sets_long_t *sets,
497 : const unsigned long d)
498 : {
499 960 : unsigned long i, c, remaining_d = d;
500 960 : set_long_t *readfrom, *readnext, *moveto, *extractto = NULL;
501 960 : size_t extracted_size = sizeof (unsigned long);
502 :
503 960 : ASSERT_ALWAYS (d > 0UL);
504 :
505 960 : if (d == 1UL)
506 : {
507 : /* d == 1 means we need to extract a set of cardinality 1, which we
508 : most likely don't have in "*sets". (FIXME: check for set of
509 : cardinality 1?) We return the set containing only zero, which
510 : can be added to any set of sets without changing the set of sums */
511 362 : if (extracted != NULL)
512 : {
513 181 : extracted->nr = 1;
514 181 : extractto = extracted->sets;
515 181 : extractto->card = 1UL;
516 181 : extractto->elem[0] = 0L;
517 : }
518 362 : if (extr_size != NULL)
519 181 : *extr_size = sizeof (unsigned long) + set_sizeof (1UL);
520 362 : return;
521 : }
522 :
523 598 : if (extracted != NULL)
524 : {
525 299 : extracted->nr = 0UL;
526 299 : extractto = extracted->sets;
527 : }
528 : /* All sets from *sets are read via *readfrom, and (assuming we actually
529 : extract them) are either copied to *extractto to *moveto */
530 598 : readfrom = moveto = sets->sets;
531 4702 : for (i = 0UL; i < sets->nr; i++)
532 : {
533 4104 : c = readfrom->card; /* readfrom->card may get garbled */
534 4104 : readnext = sets_nextset (readfrom);
535 4104 : if (remaining_d % c == 0UL)
536 : {
537 652 : if (extracted != NULL)
538 : {
539 : /* Copy this set to extractto */
540 326 : set_copy (extractto, readfrom);
541 326 : extractto = sets_nextset (extractto);
542 326 : extracted->nr++;
543 : }
544 652 : remaining_d /= c;
545 652 : extracted_size += set_sizeof (c);
546 : } else {
547 3452 : if (extracted != NULL)
548 : {
549 : /* Move this set within "*sets", filling the gaps left by
550 : extracted sets */
551 1726 : set_copy (moveto, readfrom);
552 1726 : moveto = sets_nextset (moveto);
553 : }
554 : }
555 4104 : readfrom = readnext;
556 : }
557 :
558 598 : ASSERT_ALWAYS (remaining_d == 1UL);
559 598 : if (extr_size != NULL)
560 299 : *extr_size = extracted_size;
561 598 : if (extracted != NULL)
562 299 : sets->nr -= extracted->nr;
563 : }
564 :
565 : sets_long_t *
566 480 : sets_get_factored_sorted (const unsigned long beta)
567 : {
568 : sets_long_t *sets;
569 : size_t size;
570 :
571 480 : sets_factor_coprime (NULL, &size, beta);
572 480 : sets = malloc (size);
573 480 : if (sets == NULL)
574 0 : return NULL;
575 480 : sets_factor_coprime (sets, NULL, beta);
576 :
577 480 : if (test_verbose (OUTPUT_TRACE))
578 : {
579 10 : outputf (OUTPUT_TRACE,
580 : "sets_get_factored_sorted: Factored sets before sorting are ");
581 10 : sets_print (OUTPUT_TRACE, sets);
582 : }
583 :
584 480 : sets_sort (sets);
585 :
586 480 : if (test_verbose (OUTPUT_TRACE))
587 : {
588 10 : outputf (OUTPUT_TRACE, "Factored sets after sorting are ");
589 10 : sets_print (OUTPUT_TRACE, sets);
590 : }
591 :
592 480 : return sets;
593 : }
594 :
595 : #ifdef TESTDRIVE
596 : static void
597 : selftest (const unsigned long beta)
598 : {
599 : sets_long_t *sets;
600 : set_long_t *sumset;
601 : unsigned long i, j, phibeta;
602 : mpz_t max;
603 :
604 : ASSERT_ALWAYS (beta > 0);
605 : sets = sets_get_factored_sorted (beta);
606 :
607 : /* Test that the sumset % beta is equal to (Z/betaZ)* % beta */
608 : phibeta = eulerphi (beta);
609 : sumset = malloc (set_sizeof (phibeta));
610 : if (sumset == NULL)
611 : {
612 : fprintf (stderr, "Cannot allocate memory in selftest\n");
613 : exit (1);
614 : }
615 : sets_sumset (sumset, sets);
616 : ASSERT_ALWAYS (sumset->card = phibeta);
617 : /* Also test that max (sumset) == sets_max (beta) */
618 : mpz_init (max);
619 : sets_max (max, beta);
620 : if (phibeta > 0)
621 : {
622 : long maxelem;
623 : maxelem = sumset->elem[0];
624 : for (i = 1; i < phibeta; i++)
625 : if (maxelem < sumset->elem[i])
626 : maxelem = sumset->elem[i];
627 : ASSERT_ALWAYS (mpz_cmp_si (max, maxelem) == 0);
628 : }
629 : else
630 : {
631 : ASSERT_ALWAYS (mpz_cmp_ui (max, 0UL) == 0);
632 : }
633 : mpz_clear (max);
634 :
635 : /* printf ("sumset, before reduction: ");
636 : for (i = 0; i < phibeta; i++)
637 : printf ("%ld%s", sumset->elem[i], i < phibeta-1 ? ", " : "\n"); */
638 : for (i = 0; i < phibeta; i++)
639 : {
640 : sumset->elem[i] = (sumset->elem[i] < 0L) ?
641 : beta - (long) ((unsigned long) (-sumset->elem[i]) % beta)
642 : : (unsigned long) sumset->elem[i] % beta;
643 : ASSERT_ALWAYS (sumset->elem[i] >= 0L);
644 : ASSERT_ALWAYS (sumset->elem[i] < (long) beta);
645 : }
646 : /* printf ("sumset, after reduction: ");
647 : for (i = 0; i < phibeta; i++)
648 : printf ("%ld%s", sumset->elem[i], i < phibeta-1 ? ", " : "\n"); */
649 :
650 : quicksort_long (sumset->elem, sumset->card);
651 : /* printf ("sumset, after sorting: ");
652 : for (i = 0; i < phibeta; i++)
653 : printf ("%ld%s", sumset->elem[i], i < phibeta-1 ? ", " : "\n"); */
654 :
655 : j = 0;
656 : for (i = 1; i < beta; i++)
657 : {
658 : if (gcd (i, beta) == 1)
659 : {
660 : if (sumset->elem[j] != (long) i)
661 : {
662 : printf ("sumset->elem[%ld] = %ld != %ld\n",
663 : j, sumset->elem[j], i);
664 : abort();
665 : }
666 : j++;
667 : }
668 : }
669 : free (sumset);
670 : free (sets);
671 : }
672 :
673 : int
674 : main (int argc, char **argv)
675 : {
676 : unsigned long beta;
677 : const unsigned long selftest_max = 1000;
678 : int loop = 1;
679 :
680 : ECM_STDOUT = stdout;
681 : ECM_STDERR = stderr;
682 :
683 : if (argc > 1)
684 : {
685 : beta = atol (argv[1]);
686 : loop = 0;
687 : }
688 :
689 : if (!loop)
690 : set_verbose (OUTPUT_TRACE);
691 :
692 : if (!loop)
693 : selftest (beta);
694 : else
695 : {
696 : printf ("Testing beta = 1, ..., %lu\n", selftest_max);
697 : for (beta = 1; beta < selftest_max; beta++)
698 : selftest (beta);
699 : }
700 :
701 : return 0;
702 : }
703 : #endif
|
