1 /* Factoring of uintmax_t numbers. Generation of needed tables.
3 Contributed to the GNU project by Torbjörn Granlund and Niels Möller
4 Contains code from GNU MP.
6 Copyright 2012 Free Software Foundation, Inc.
8 This program is free software; you can redistribute it and/or modify it under
9 the terms of the GNU General Public License as published by the Free Software
10 Foundation; either version 3 of the License, or (at your option) any later
13 This program is distributed in the hope that it will be useful, but WITHOUT ANY
14 WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
15 PARTICULAR PURPOSE. See the GNU General Public License for more details.
17 You should have received a copy of the GNU General Public License along with
18 this program. If not, see http://www.gnu.org/licenses/. */
30 /* Deactivate config.h's "rpl_"-prefixed definitions of these symbols. */
35 /* An unsigned type that is no narrower than 32 bits and no narrower
36 than unsigned int. It's best to make it as wide as possible.
37 For GCC 4.6 and later, use a heuristic to guess whether unsigned
38 __int128 works on your platform. If this heuristic does not work
39 for you, please report a bug; in the meantime compile with, e.g.,
40 -Dwide_uint='unsigned __int128' to override the heuristic. */
42 # if 4 < __GNUC__ + (6 <= __GNUC_MINOR__) && ULONG_MAX >> 31 >> 31 >> 1 != 0
43 typedef unsigned __int128 wide_uint
;
45 typedef uintmax_t wide_uint
;
52 wide_uint pinv
; /* Inverse mod b = 2^{bitsize of wide_uint} */
53 wide_uint lim
; /* floor ((wide_uint) -1 / p) */
56 static wide_uint _GL_ATTRIBUTE_CONST
59 wide_uint x
= 0xf5397db1 >> (4*((a
/2) & 0x7));
62 wide_uint y
= 2*x
- x
*x
*a
;
70 process_prime (struct prime
*info
, unsigned p
)
74 info
->pinv
= binvert (p
);
79 print_wide_uint (wide_uint n
, int nesting
, unsigned wide_uint_bits
)
81 /* Number of bits per integer literal. 8 is too many, because
82 uintmax_t is 32 bits on some machines so we cannot shift by 32 bits.
84 int hex_digits_per_literal
= 7;
85 int bits_per_literal
= hex_digits_per_literal
* 4;
87 unsigned remainder
= n
& ((1 << bits_per_literal
) - 1);
91 int needs_parentheses
= n
>> bits_per_literal
>> bits_per_literal
!= 0;
92 if (needs_parentheses
)
94 print_wide_uint (n
>> bits_per_literal
, nesting
+ 1, wide_uint_bits
);
95 if (needs_parentheses
)
96 printf (")\n%*s", nesting
+ 3, "");
97 printf (" << %d | ", bits_per_literal
);
101 printf ("(uintmax_t) ");
102 hex_digits_per_literal
103 = ((wide_uint_bits
- 1) % bits_per_literal
) % 4 + 1;
106 printf ("0x%0*xU", hex_digits_per_literal
, remainder
);
110 output_primes (const struct prime
*primes
, unsigned nprimes
)
116 /* Compute wide_uint_bits by repeated shifting, rather than by
117 multiplying sizeof by CHAR_BIT, as this works even if the
118 wide_uint representation has holes. */
119 unsigned wide_uint_bits
= 0;
121 for (wide_uint_bits
= 0; mask
; wide_uint_bits
++)
124 puts ("/* Generated file -- DO NOT EDIT */\n");
125 printf ("#define WIDE_UINT_BITS %u\n", wide_uint_bits
);
127 for (i
= 0, p
= 2; i
< nprimes
; i
++)
129 unsigned int d8
= i
+ 8 < nprimes
? primes
[i
+ 8].p
- primes
[i
].p
: 0xff;
130 if (255 < d8
) /* this happens at 668221 */
132 printf ("P (%u, %u,\n (", primes
[i
].p
- p
, d8
);
133 print_wide_uint (primes
[i
].pinv
, 0, wide_uint_bits
);
134 printf ("),\n UINTMAX_MAX / %d)\n", primes
[i
].p
);
138 printf ("\n#undef FIRST_OMITTED_PRIME\n");
140 /* Find next prime */
144 for (i
= 0, is_prime
= 1; is_prime
; i
++)
146 if (primes
[i
].p
* primes
[i
].p
> p
)
148 if (p
* primes
[i
].pinv
<= primes
[i
].lim
)
157 printf ("#define FIRST_OMITTED_PRIME %u\n", p
);
163 void *p
= malloc (s
);
167 fprintf (stderr
, "Virtual memory exhausted.\n");
172 main (int argc
, char **argv
)
179 struct prime
*prime_list
;
184 fprintf (stderr
, "Usage: %s LIMIT\n"
185 "Produces a list of odd primes <= LIMIT\n", argv
[0]);
188 limit
= atoi (argv
[1]);
197 /* sieve[i] represents 3+2*i */
198 sieve
= xalloc (size
);
199 memset (sieve
, 1, size
);
201 prime_list
= xalloc (size
* sizeof (*prime_list
));
204 for (i
= 0; i
< size
;)
209 process_prime (&prime_list
[nprimes
++], p
);
211 for (j
= (p
*p
- 3)/2; j
< size
; j
+= p
)
214 while (i
< size
&& sieve
[++i
] == 0)
218 output_primes (prime_list
, nprimes
);
220 if (ferror (stdout
) + fclose (stdout
))
222 fprintf (stderr
, "write error: %s\n", strerror (errno
));