1 /* Generate perfect square testing data.
3 Copyright 2002-2004, 2012, 2014 Free Software Foundation, Inc.
5 This file is part of the GNU MP Library.
7 The GNU MP Library is free software; you can redistribute it and/or modify
8 it under the terms of either:
10 * the GNU Lesser General Public License as published by the Free
11 Software Foundation; either version 3 of the License, or (at your
12 option) any later version.
16 * the GNU General Public License as published by the Free Software
17 Foundation; either version 2 of the License, or (at your option) any
20 or both in parallel, as here.
22 The GNU MP Library is distributed in the hope that it will be useful, but
23 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
24 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
27 You should have received copies of the GNU General Public License and the
28 GNU Lesser General Public License along with the GNU MP Library. If not,
29 see https://www.gnu.org/licenses/. */
34 #include "bootstrap.c"
37 /* The aim of this program is to choose either mpn_mod_34lsub1 or mpn_mod_1
38 (plus a PERFSQR_PP modulus), and generate tables indicating quadratic
39 residues and non-residues modulo small factors of that modulus.
41 For the usual 32 or 64 bit cases mpn_mod_34lsub1 gets used. That
42 function exists specifically because 2^24-1 and 2^48-1 have nice sets of
43 prime factors. For other limb sizes it's considered, but if it doesn't
44 have good factors then mpn_mod_1 will be used instead.
46 When mpn_mod_1 is used, the modulus PERFSQR_PP is created from a
47 selection of small primes, chosen to fill PERFSQR_MOD_BITS of a limb,
48 with that bit count chosen so (2*GMP_LIMB_BITS)*2^PERFSQR_MOD_BITS <=
49 GMP_LIMB_MAX, allowing PERFSQR_MOD_IDX in mpn/generic/perfsqr.c to do its
50 calculation within a single limb.
52 In either case primes can be combined to make divisors. The table data
53 then effectively indicates remainders which are quadratic residues mod
54 all the primes. This sort of combining reduces the number of steps
55 needed after mpn_mod_34lsub1 or mpn_mod_1, saving code size and time.
56 Nothing is gained or lost in terms of detections, the same total fraction
57 of non-residues will be identified.
59 Nothing particularly sophisticated is attempted for combining factors to
60 make divisors. This is probably a kind of knapsack problem so it'd be
61 too hard to attempt anything completely general. For the usual 32 and 64
62 bit limbs we get a good enough result just pairing the biggest and
63 smallest which fit together, repeatedly.
65 Another aim is to get powerful combinations, ie. divisors which identify
66 biggest fraction of non-residues, and have those run first. Again for
67 the usual 32 and 64 bits it seems good enough just to pair for big
68 divisors then sort according to the resulting fraction of non-residues
71 Also in this program, a table sq_res_0x100 of residues modulo 256 is
72 generated. This simply fills bits into limbs of the appropriate
73 build-time GMP_LIMB_BITS each.
78 /* Normally we aren't using const in gen*.c programs, so as not to have to
79 bother figuring out if it works, but using it with f_cmp_divisor and
80 f_cmp_fraction avoids warnings from the qsort calls. */
82 /* Same tests as gmp.h. */
83 #if defined (__STDC__) \
84 || defined (__cplusplus) \
87 || (defined (__mips) && defined (_SYSTYPE_SVR4)) \
88 || defined (_MSC_VER) \
98 mpz_t
*sq_res_0x100
; /* table of limbs */
99 int nsq_res_0x100
; /* elements in sq_res_0x100 array */
100 int sq_res_0x100_num
; /* squares in sq_res_0x100 */
101 double sq_res_0x100_fraction
; /* sq_res_0x100_num / 256 */
103 int mod34_bits
; /* 3*GMP_NUMB_BITS/4 */
104 int mod_bits
; /* bits from PERFSQR_MOD_34 or MOD_PP */
105 int max_divisor
; /* all divisors <= max_divisor */
106 int max_divisor_bits
; /* ceil(log2(max_divisor)) */
107 double total_fraction
; /* of squares */
108 mpz_t pp
; /* product of primes, or 0 if mod_34lsub1 used */
109 mpz_t pp_norm
; /* pp shifted so NUMB high bit set */
110 mpz_t pp_inverted
; /* invert_limb style inverse */
111 mpz_t mod_mask
; /* 2^mod_bits-1 */
112 char mod34_excuse
[128]; /* why mod_34lsub1 not used (if it's not) */
114 /* raw list of divisors of 2^mod34_bits-1 or pp, just to show in a comment */
119 struct rawfactor_t
*rawfactor
;
122 /* factors of 2^mod34_bits-1 or pp and associated data, after combining etc */
125 mpz_t inverse
; /* 1/divisor mod 2^mod_bits */
126 mpz_t mask
; /* indicating squares mod divisor */
127 double fraction
; /* squares/total */
129 struct factor_t
*factor
;
130 int nfactor
; /* entries in use in factor array */
131 int factor_alloc
; /* entries allocated to factor array */
135 f_cmp_divisor (const void *parg
, const void *qarg
)
137 const struct factor_t
*p
, *q
;
138 p
= (const struct factor_t
*) parg
;
139 q
= (const struct factor_t
*) qarg
;
140 if (p
->divisor
> q
->divisor
)
142 else if (p
->divisor
< q
->divisor
)
149 f_cmp_fraction (const void *parg
, const void *qarg
)
151 const struct factor_t
*p
, *q
;
152 p
= (const struct factor_t
*) parg
;
153 q
= (const struct factor_t
*) qarg
;
154 if (p
->fraction
> q
->fraction
)
156 else if (p
->fraction
< q
->fraction
)
162 /* Remove array[idx] by copying the remainder down, and adjust narray
164 #define COLLAPSE_ELEMENT(array, idx, narray) \
166 memmove (&(array)[idx], \
168 ((narray)-((idx)+1)) * sizeof (array[0])); \
173 /* return n*2^p mod m */
175 mul_2exp_mod (int n
, int p
, int m
)
182 /* return -n mod m */
184 neg_mod (int n
, int m
)
186 assert (n
>= 0 && n
< m
);
187 return (n
== 0 ? 0 : m
-n
);
190 /* Set "mask" to a value such that "mask & (1<<idx)" is non-zero if
191 "-(idx<<mod_bits)" can be a square modulo m. */
193 square_mask (mpz_t mask
, int m
)
197 p
= mul_2exp_mod (1, mod_bits
, m
);
200 mpz_set_ui (mask
, 0L);
201 for (i
= 0; i
< m
; i
++)
205 mpz_setbit (mask
, (unsigned long) idx
);
210 generate_sq_res_0x100 (int limb_bits
)
214 nsq_res_0x100
= (0x100 + limb_bits
- 1) / limb_bits
;
215 sq_res_0x100
= (mpz_t
*) xmalloc (nsq_res_0x100
* sizeof (*sq_res_0x100
));
217 for (i
= 0; i
< nsq_res_0x100
; i
++)
218 mpz_init_set_ui (sq_res_0x100
[i
], 0L);
220 for (i
= 0; i
< 0x100; i
++)
222 res
= (i
* i
) % 0x100;
223 mpz_setbit (sq_res_0x100
[res
/ limb_bits
],
224 (unsigned long) (res
% limb_bits
));
227 sq_res_0x100_num
= 0;
228 for (i
= 0; i
< nsq_res_0x100
; i
++)
229 sq_res_0x100_num
+= mpz_popcount (sq_res_0x100
[i
]);
230 sq_res_0x100_fraction
= (double) sq_res_0x100_num
/ 256.0;
234 generate_mod (int limb_bits
, int nail_bits
)
236 int numb_bits
= limb_bits
- nail_bits
;
239 mpz_init_set_ui (pp
, 0L);
240 mpz_init_set_ui (pp_norm
, 0L);
241 mpz_init_set_ui (pp_inverted
, 0L);
243 /* no more than limb_bits many factors in a one limb modulus (and of
244 course in reality nothing like that many) */
245 factor_alloc
= limb_bits
;
246 factor
= (struct factor_t
*) xmalloc (factor_alloc
* sizeof (*factor
));
247 rawfactor
= (struct rawfactor_t
*) xmalloc (factor_alloc
* sizeof (*rawfactor
));
249 if (numb_bits
% 4 != 0)
251 strcpy (mod34_excuse
, "GMP_NUMB_BITS % 4 != 0");
255 max_divisor
= 2*limb_bits
;
256 max_divisor_bits
= log2_ceil (max_divisor
);
258 if (numb_bits
/ 4 < max_divisor_bits
)
260 /* Wind back to one limb worth of max_divisor, if that will let us use
262 max_divisor
= limb_bits
;
263 max_divisor_bits
= log2_ceil (max_divisor
);
265 if (numb_bits
/ 4 < max_divisor_bits
)
267 strcpy (mod34_excuse
, "GMP_NUMB_BITS / 4 too small");
273 /* Can use mpn_mod_34lsub1, find small factors of 2^mod34_bits-1. */
277 mod34_bits
= (numb_bits
/ 4) * 3;
279 /* mpn_mod_34lsub1 returns a full limb value, PERFSQR_MOD_34 folds it at
280 the mod34_bits mark, adding the two halves for a remainder of at most
281 mod34_bits+1 many bits */
282 mod_bits
= mod34_bits
+ 1;
284 mpz_init_set_ui (m
, 1L);
285 mpz_mul_2exp (m
, m
, mod34_bits
);
286 mpz_sub_ui (m
, m
, 1L);
291 for (i
= 3; i
<= max_divisor
; i
+=2)
296 mpz_tdiv_qr_ui (q
, r
, m
, (unsigned long) i
);
297 if (mpz_sgn (r
) != 0)
300 /* if a repeated prime is found it's used as an i^n in one factor */
305 if (divisor
> max_divisor
/ i
)
309 mpz_tdiv_qr_ui (q
, r
, m
, (unsigned long) i
);
311 while (mpz_sgn (r
) == 0);
313 assert (nrawfactor
< factor_alloc
);
314 rawfactor
[nrawfactor
].divisor
= i
;
315 rawfactor
[nrawfactor
].multiplicity
= multiplicity
;
328 sprintf (mod34_excuse
, "only %d small factor%s",
329 nrawfactor
, nrawfactor
== 1 ? "" : "s");
332 /* reset to two limbs of max_divisor, in case the mpn_mod_34lsub1 code
333 tried with just one */
334 max_divisor
= 2*limb_bits
;
335 max_divisor_bits
= log2_ceil (max_divisor
);
339 mod_bits
= MIN (numb_bits
, limb_bits
- max_divisor_bits
);
341 /* one copy of each small prime */
343 for (i
= 3; i
<= max_divisor
; i
+=2)
348 mpz_mul_ui (new_pp
, pp
, (unsigned long) i
);
349 if (mpz_sizeinbase (new_pp
, 2) > mod_bits
)
351 mpz_set (pp
, new_pp
);
353 assert (nrawfactor
< factor_alloc
);
354 rawfactor
[nrawfactor
].divisor
= i
;
355 rawfactor
[nrawfactor
].multiplicity
= 1;
359 /* Plus an extra copy of one or more of the primes selected, if that
360 still fits in max_divisor and the total in mod_bits. Usually only
361 3 or 5 will be candidates */
362 for (i
= nrawfactor
-1; i
>= 0; i
--)
364 if (rawfactor
[i
].divisor
> max_divisor
/ rawfactor
[i
].divisor
)
366 mpz_mul_ui (new_pp
, pp
, (unsigned long) rawfactor
[i
].divisor
);
367 if (mpz_sizeinbase (new_pp
, 2) > mod_bits
)
369 mpz_set (pp
, new_pp
);
371 rawfactor
[i
].multiplicity
++;
374 mod_bits
= mpz_sizeinbase (pp
, 2);
376 mpz_set (pp_norm
, pp
);
377 while (mpz_sizeinbase (pp_norm
, 2) < numb_bits
)
378 mpz_add (pp_norm
, pp_norm
, pp_norm
);
380 mpz_preinv_invert (pp_inverted
, pp_norm
, numb_bits
);
385 /* start the factor array */
386 for (i
= 0; i
< nrawfactor
; i
++)
389 assert (nfactor
< factor_alloc
);
390 factor
[nfactor
].divisor
= 1;
391 for (j
= 0; j
< rawfactor
[i
].multiplicity
; j
++)
392 factor
[nfactor
].divisor
*= rawfactor
[i
].divisor
;
397 /* Combine entries in the factor array. Combine the smallest entry with
398 the biggest one that will fit with it (ie. under max_divisor), then
399 repeat that with the new smallest entry. */
400 qsort (factor
, nfactor
, sizeof (factor
[0]), f_cmp_divisor
);
401 for (i
= nfactor
-1; i
>= 1; i
--)
403 if (factor
[i
].divisor
<= max_divisor
/ factor
[0].divisor
)
405 factor
[0].divisor
*= factor
[i
].divisor
;
406 COLLAPSE_ELEMENT (factor
, i
, nfactor
);
411 total_fraction
= 1.0;
412 for (i
= 0; i
< nfactor
; i
++)
414 mpz_init (factor
[i
].inverse
);
415 mpz_invert_ui_2exp (factor
[i
].inverse
,
416 (unsigned long) factor
[i
].divisor
,
417 (unsigned long) mod_bits
);
419 mpz_init (factor
[i
].mask
);
420 square_mask (factor
[i
].mask
, factor
[i
].divisor
);
422 /* fraction of possible squares */
423 factor
[i
].fraction
= (double) mpz_popcount (factor
[i
].mask
)
426 /* total fraction of possible squares */
427 total_fraction
*= factor
[i
].fraction
;
430 /* best tests first (ie. smallest fraction) */
431 qsort (factor
, nfactor
, sizeof (factor
[0]), f_cmp_fraction
);
435 print (int limb_bits
, int nail_bits
)
440 printf ("/* This file generated by gen-psqr.c - DO NOT EDIT. */\n");
443 printf ("#if GMP_LIMB_BITS != %d || GMP_NAIL_BITS != %d\n",
444 limb_bits
, nail_bits
);
445 printf ("Error, error, this data is for %d bit limb and %d bit nail\n",
446 limb_bits
, nail_bits
);
450 printf ("/* Non-zero bit indicates a quadratic residue mod 0x100.\n");
451 printf (" This test identifies %.2f%% as non-squares (%d/256). */\n",
452 (1.0 - sq_res_0x100_fraction
) * 100.0,
453 0x100 - sq_res_0x100_num
);
454 printf ("static const mp_limb_t\n");
455 printf ("sq_res_0x100[%d] = {\n", nsq_res_0x100
);
456 for (i
= 0; i
< nsq_res_0x100
; i
++)
458 printf (" CNST_LIMB(0x");
459 mpz_out_str (stdout
, 16, sq_res_0x100
[i
]);
465 if (mpz_sgn (pp
) != 0)
467 printf ("/* mpn_mod_34lsub1 not used due to %s */\n", mod34_excuse
);
468 printf ("/* PERFSQR_PP = ");
471 printf ("/* 2^%d-1 = ", mod34_bits
);
472 for (i
= 0; i
< nrawfactor
; i
++)
476 printf ("%d", rawfactor
[i
].divisor
);
477 if (rawfactor
[i
].multiplicity
!= 1)
478 printf ("^%d", rawfactor
[i
].multiplicity
);
480 printf (" %s*/\n", mpz_sgn (pp
) == 0 ? "... " : "");
482 printf ("#define PERFSQR_MOD_BITS %d\n", mod_bits
);
483 if (mpz_sgn (pp
) != 0)
485 printf ("#define PERFSQR_PP CNST_LIMB(0x");
486 mpz_out_str (stdout
, 16, pp
);
488 printf ("#define PERFSQR_PP_NORM CNST_LIMB(0x");
489 mpz_out_str (stdout
, 16, pp_norm
);
491 printf ("#define PERFSQR_PP_INVERTED CNST_LIMB(0x");
492 mpz_out_str (stdout
, 16, pp_inverted
);
500 printf ("/* This test identifies %.2f%% as non-squares. */\n",
501 (1.0 - total_fraction
) * 100.0);
502 printf ("#define PERFSQR_MOD_TEST(up, usize) \\\n");
503 printf (" do { \\\n");
504 printf (" mp_limb_t r; \\\n");
505 if (mpz_sgn (pp
) != 0)
506 printf (" PERFSQR_MOD_PP (r, up, usize); \\\n");
508 printf (" PERFSQR_MOD_34 (r, up, usize); \\\n");
510 for (i
= 0; i
< nfactor
; i
++)
513 printf (" /* %5.2f%% */ \\\n",
514 (1.0 - factor
[i
].fraction
) * 100.0);
516 printf (" PERFSQR_MOD_%d (r, CNST_LIMB(%2d), CNST_LIMB(0x",
517 factor
[i
].divisor
<= limb_bits
? 1 : 2,
519 mpz_out_str (stdout
, 16, factor
[i
].inverse
);
521 printf (" CNST_LIMB(0x");
523 if ( factor
[i
].divisor
<= limb_bits
)
525 mpz_out_str (stdout
, 16, factor
[i
].mask
);
529 mpz_tdiv_r_2exp (mlo
, factor
[i
].mask
, (unsigned long) limb_bits
);
530 mpz_tdiv_q_2exp (mhi
, factor
[i
].mask
, (unsigned long) limb_bits
);
531 mpz_out_str (stdout
, 16, mhi
);
532 printf ("), CNST_LIMB(0x");
533 mpz_out_str (stdout
, 16, mlo
);
538 printf (" } while (0)\n");
541 printf ("/* Grand total sq_res_0x100 and PERFSQR_MOD_TEST, %.2f%% non-squares. */\n",
542 (1.0 - (total_fraction
* 44.0/256.0)) * 100.0);
545 printf ("/* helper for tests/mpz/t-perfsqr.c */\n");
546 printf ("#define PERFSQR_DIVISORS { 256,");
547 for (i
= 0; i
< nfactor
; i
++)
548 printf (" %d,", factor
[i
].divisor
);
557 main (int argc
, char *argv
[])
559 int limb_bits
, nail_bits
;
563 fprintf (stderr
, "Usage: gen-psqr <limbbits> <nailbits>\n");
567 limb_bits
= atoi (argv
[1]);
568 nail_bits
= atoi (argv
[2]);
572 || nail_bits
>= limb_bits
)
574 fprintf (stderr
, "Invalid limb/nail bits: %d %d\n",
575 limb_bits
, nail_bits
);
579 generate_sq_res_0x100 (limb_bits
);
580 generate_mod (limb_bits
, nail_bits
);
582 print (limb_bits
, nail_bits
);