1 /* Reads integers "N", or "p N", or "p O N" from stdin and checks N for
2 smoothness, where smooth is defined as
3 For given B1, B2, k, d,
4 N is smooth if c=N/gcd(N,E) is 1 or is prime and B1 < c <= B2,
5 or c is coprime to d and <d/2, or c divides d,
6 where E = k * lcm(1, 2, 3, ..., B1).
7 By default, k=1 and d=1, so that they have no effect.
8 The purpose of the value d is including those N where ECM would find
9 the factor during stage 2 initialisation.
10 With -pmin and -pmax, test only those input where pmin <= p <= pmax,
11 (or pmin <= N <= pmax if input has no p value).
12 Prints the number of smooth numbers, and with -v option each "p O N"
13 where N is smooth.
14 Several B2 can be specified on the command line (but only one B1 value).
15 For each B2, the number of smooth N are counted.
16 If an "m" value can be specified, the number of input numbers and smooth
17 (w.r.t. the largest B2) numbers with p==r (mod m) for the possible r
18 are also counted and printed.
19 Also prints the average exponent of primes up to 19 in N and in O/N
20 (which requires N|O).
21 */
22
23 #include "cado.h" // IWYU pragma: keep
24 #include <stdio.h>
25 #include <stdlib.h>
26 #include <string.h>
27 #include <limits.h>
28 #include <gmp.h>
29 #include "macros.h"
30 #include "gcd.h"
31
32 #define NR_EXPONENTS 8
33 #define EXP_PRIMES {2,3,5,7,11,13,17,19}
34
35 static const unsigned long exp_primes[NR_EXPONENTS] = EXP_PRIMES;
36
37 static int
isprime(const unsigned long N,mpz_t t)38 isprime (const unsigned long N, mpz_t t)
39 {
40 mpz_set_ui (t, N);
41 return mpz_probab_prime_p (t, 1);
42 }
43
44
45 static unsigned int
valuation(const unsigned long n,const unsigned long p)46 valuation (const unsigned long n, const unsigned long p)
47 {
48 unsigned long c = n;
49 unsigned int i = 0;
50 while (c % p == 0)
51 {
52 c /= p;
53 i++;
54 }
55 return i;
56 }
57
58 typedef struct {
59 unsigned long input, *smooth;
60 unsigned int exp_order[NR_EXPONENTS];
61 unsigned int exp_index[NR_EXPONENTS];
62 } stats_t;
63
64
65 int
main(int argc,char ** argv)66 main (int argc, char **argv)
67 {
68 mpz_t E, m_c;
69 unsigned long N, O, c, c_gcdiv_d, p, i, j, imin = 0, imax = 0;
70 /* Parameters */
71 unsigned long B1, *B2, maxB2, k = 1, d = 1, pmin = 0, pmax = ULONG_MAX,
72 m = 1;
73 int quiet = 0, verbose = 0;
74 unsigned int nr_B2 = 0;
75 stats_t *stats;
76 int scanf_code;
77 char *d_coprime;
78 char buf[256];
79
80 /* Parse arguments */
81
82 while (argc > 1 && argv[1][0] == '-')
83 {
84 if (argv[1][1] == 'v')
85 {
86 verbose = 1;
87 argc--;
88 argv++;
89 }
90 else if (argv[1][1] == 'q')
91 {
92 quiet = 1;
93 argc--;
94 argv++;
95 }
96 else if (argc > 2 && argv[1][1] == 'k')
97 {
98 k = strtoul (argv[2], NULL, 10);
99 argc -= 2;
100 argv += 2;
101 }
102 else if (argc > 2 && argv[1][1] == 'd')
103 {
104 d = strtoul (argv[2], NULL, 10);
105 if (d == 0) {
106 fprintf (stderr, "d == 0 does not make sense\n");
107 exit (EXIT_FAILURE);
108 }
109 argc -= 2;
110 argv += 2;
111 }
112 else if (argc > 2 && argv[1][1] == 'm')
113 {
114 m = strtoul (argv[2], NULL, 10);
115 if (m == 0) {
116 fprintf (stderr, "m == 0 does not make sense\n");
117 exit (EXIT_FAILURE);
118 }
119 #ifdef __COVERITY__
120 /* any non-zero value will do */
121 __coverity_mark_pointee_as_sanitized__(&m, GENERIC); // allocation, divisor, loop bound
122 #endif
123 argc -= 2;
124 argv += 2;
125 }
126 else if (argc > 2 && strcmp (argv[1], "-pmin") == 0)
127 {
128 pmin = strtoul (argv[2], NULL, 10);
129 argc -= 2;
130 argv += 2;
131 }
132 else if (argc > 2 && strcmp (argv[1], "-pmax") == 0)
133 {
134 pmax = strtoul (argv[2], NULL, 10);
135 argc -= 2;
136 argv += 2;
137 }
138 else
139 {
140 fprintf (stderr, "Unknown option %s\n", argv[1]);
141 exit (EXIT_FAILURE);
142 }
143 }
144
145 if (argc < 3)
146 {
147 printf ("Please provide B1 and B2 values\n");
148 exit (EXIT_FAILURE);
149 }
150
151 /* Get B1 and all the B2 values */
152 B1 = strtoul (argv[1], NULL, 10);
153 nr_B2 = argc - 2;
154 B2 = malloc(nr_B2 * sizeof(unsigned long));
155 maxB2 = 0;
156 for (i = 0; i < nr_B2; i++)
157 {
158 B2[i] = strtoul (argv[2 + i], NULL, 10);
159 maxB2 = MAX(B2[i], maxB2);
160 }
161
162 d_coprime = malloc (d * sizeof(char));
163 for (i = 1; i < d; i++)
164 d_coprime[i] = (gcd_ul(i,d) == 1) ? 1 : 0;
165 if (d > 1)
166 {
167 imin = (B1 + d / 2) / d;
168 imax = (maxB2 - d / 2) / d;
169 }
170
171 /* Allocate a stat_t for each residue class (mod m) */
172 stats = malloc (m * sizeof (stats_t));
173 for (i = 0; i < m; i++)
174 {
175 stats[i].input = 0;
176 for (j = 0; j < NR_EXPONENTS; j++)
177 stats[i].exp_order[j] = stats[i].exp_index[j] = 0;
178
179 stats[i].smooth = malloc (nr_B2 * sizeof(unsigned long));
180 for (j = 0; j < nr_B2; j++)
181 stats[i].smooth[j] = 0;
182 }
183
184
185 /* Compute E */
186 mpz_init (E);
187 mpz_set_ui (E, 1UL);
188 for (i = 2UL; i <= B1; i = (i + 1UL) | 1UL)
189 if (mpz_gcd_ui (NULL, E, i) == 1UL) /* A prime? */
190 {
191 for (j = 1UL; j <= B1 / i; j *= i);
192 mpz_mul_ui (E, E, j);
193 }
194 mpz_mul_ui (E, E, k);
195
196 if (!quiet)
197 {
198 fprintf (stderr, "B1 = %lu, ", B1);
199 for (i = 0; i < nr_B2; i++)
200 fprintf (stderr, "B2[%lu] = %lu, ", i, B2[i]);
201 gmp_fprintf (stderr, "k = %lu, pmin = %lu, pmax = %lu, E = %Zd\n",
202 k, pmin, pmax, E);
203 }
204
205
206 /* Test input numbers */
207 mpz_init (m_c);
208 while (!feof(stdin))
209 {
210 if (fgets (buf, 256, stdin) == NULL)
211 break;
212 scanf_code = sscanf (buf, "%lu %lu %lu\n", &p, &O, &N);
213 if (pmin <= p && p <= pmax)
214 {
215 unsigned long r;
216 int was_smooth = 0, ispr = 0;
217 /* input line may contain only "N", or "p N", or "p O N" */
218 if (scanf_code == 1)
219 N = p;
220 else if (scanf_code == 2)
221 N = O;
222
223 r = p % m;
224 stats[r].input++;
225
226 c = N / mpz_gcd_ui (NULL, E, N);
227 c_gcdiv_d = c / gcd_ul (c, d);
228 ispr = (B1 < c && c <= maxB2 && isprime(c, m_c));
229
230 for (i = 0; i < nr_B2; i++)
231 if (c_gcdiv_d == 1UL || /* c == 1 (factor found after stage 1) || c | d */
232 (imin <= c_gcdiv_d && c_gcdiv_d <= imax) || /* factor appears in list of idP */
233 (c < d / 2 && d_coprime[c] == 1) || /* factor appears in list of jP */
234 (ispr && c <= B2[i])) /* factor appears normally in stage 2 */
235 {
236 stats[r].smooth[i]++;
237 was_smooth = 1;
238 }
239
240 /* If order was smooth for any B2 and verbose is set, print it */
241 if (was_smooth && verbose)
242 printf ("%lu %lu\n", p, N);
243
244 if (scanf_code == 3)
245 {
246 ASSERT (O % N == 0);
247 for (i = 0; i < NR_EXPONENTS; i++)
248 {
249 stats[r].exp_order[i] += valuation (O, exp_primes[i]);
250 stats[r].exp_index[i] += valuation (O/N, exp_primes[i]);
251 }
252 }
253 }
254 }
255 mpz_clear (E);
256 mpz_clear (m_c);
257
258 if (!quiet)
259 {
260 unsigned long r;
261 for (r = 0; r < m; r++)
262 {
263 if (stats[r].input > 0)
264 {
265 printf ("%lu (mod %lu): %lu input numbers\n",
266 r, m, stats[r].input);
267 for (i = 0; i < nr_B2; i++)
268 printf ("%lu (mod %lu): %lu %lu,%lu-smooth numbers, ratio %f\n",
269 r, m, stats[r].smooth[i], B1, B2[i],
270 (double) stats[r].smooth[i] / (double) stats[r].input);
271 }
272
273 if (stats[r].exp_order[0] > 0)
274 {
275 printf ("Sum of valuation of q in O:\n");
276 for (i = 0; i < NR_EXPONENTS; i++)
277 if (stats[r].exp_order[i])
278 printf ("q = %lu: %u ", exp_primes[i], stats[r].exp_order[i]);
279 printf ("\n");
280 }
281 if (stats[r].exp_index[0] > 0)
282 {
283 printf ("Sum of valuation of q in N/O (index of starting element):\n");
284 for (i = 0; i < NR_EXPONENTS; i++)
285 if (stats[r].exp_index[i])
286 printf ("q = %lu: %u ", exp_primes[i], stats[r].exp_index[i]);
287 printf ("\n");
288 }
289 }
290 }
291
292 /* Free memory */
293 free (B2);
294 B2 = NULL;
295 free (d_coprime);
296 d_coprime = NULL;
297 for (i = 0; i < m; i++)
298 {
299 free(stats[i].smooth);
300 stats[i].smooth = NULL;
301 }
302 free (stats);
303 stats = NULL;
304
305 return 0;
306 }
307