1 /* $NetBSD: factor.c,v 1.38 2020/10/12 13:54:51 christos Exp $ */
2 /*
3 * Copyright (c) 1989, 1993
4 * The Regents of the University of California. All rights reserved.
5 *
6 * This code is derived from software contributed to Berkeley by
7 * Landon Curt Noll.
8 *
9 * Redistribution and use in source and binary forms, with or without
10 * modification, are permitted provided that the following conditions
11 * are met:
12 * 1. Redistributions of source code must retain the above copyright
13 * notice, this list of conditions and the following disclaimer.
14 * 2. Redistributions in binary form must reproduce the above copyright
15 * notice, this list of conditions and the following disclaimer in the
16 * documentation and/or other materials provided with the distribution.
17 * 3. Neither the name of the University nor the names of its contributors
18 * may be used to endorse or promote products derived from this software
19 * without specific prior written permission.
20 *
21 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
22 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
23 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
24 * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
25 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
26 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
27 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
28 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
29 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
30 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
31 * SUCH DAMAGE.
32 */
33
34 #ifndef lint
35 #include <sys/cdefs.h>
36 #ifdef __COPYRIGHT
37 __COPYRIGHT("@(#) Copyright (c) 1989, 1993\
38 The Regents of the University of California. All rights reserved.");
39 #endif
40 #ifdef __SCCSID
41 __SCCSID("@(#)factor.c 8.4 (Berkeley) 5/4/95");
42 #endif
43 #ifdef __RCSID
44 __RCSID("$NetBSD: factor.c,v 1.38 2020/10/12 13:54:51 christos Exp $");
45 #endif
46 #ifdef __FBSDID
47 __FBSDID("$FreeBSD: head/usr.bin/factor/factor.c 356666 2020-01-12 20:25:11Z gad $");
48 #endif
49 #endif /* not lint */
50
51 /*
52 * factor - factor a number into primes
53 *
54 * By: Landon Curt Noll chongo@toad.com, ...!{sun,tolsoft}!hoptoad!chongo
55 *
56 * chongo <for a good prime call: 391581 * 2^216193 - 1> /\oo/\
57 *
58 * usage:
59 * factor [-hx] [number] ...
60 *
61 * The form of the output is:
62 *
63 * number: factor1 factor1 factor2 factor3 factor3 factor3 ...
64 *
65 * where factor1 <= factor2 <= factor3 <= ...
66 *
67 * If the -h flag is specified, the output is in "human-readable" format.
68 * Duplicate factors are printed in the form of x^n.
69 *
70 * If the -x flag is specified numbers are printed in hex.
71 *
72 * If no number args are given, the list of numbers are read from stdin.
73 */
74
75 #include <ctype.h>
76 #include <err.h>
77 #include <errno.h>
78 #include <inttypes.h>
79 #include <limits.h>
80 #include <stdbool.h>
81 #include <stdio.h>
82 #include <stdlib.h>
83 #include <unistd.h>
84
85 #include "primes.h"
86
87 #ifdef HAVE_OPENSSL
88
89 #include <openssl/bn.h>
90
91 #define PRIME_CHECKS 5
92
93 /* print factors for big numbers */
94 static void pollard_pminus1(BIGNUM *, int, int);
95
96 #else
97
98 typedef long BIGNUM;
99 typedef u_long BN_ULONG;
100
101 #define BN_CTX int
102 #define BN_CTX_new() NULL
103 #define BN_new() calloc(sizeof(BIGNUM), 1)
104 #define BN_free(a) free(a)
105
106 #define BN_copy(a, b) *(a) = *(b)
107 #define BN_cmp(a, b) (*(a) - *(b))
108 #define BN_is_zero(v) (*(v) == 0)
109 #define BN_is_one(v) (*(v) == 1)
110 #define BN_mod_word(a, b) (*(a) % (b))
111
112 static BIGNUM *BN_dup(const BIGNUM *);
113 static int BN_dec2bn(BIGNUM **, const char *);
114 static int BN_hex2bn(BIGNUM **, const char *);
115 static BN_ULONG BN_div_word(BIGNUM *, BN_ULONG);
116 static void BN_print_fp(FILE *, const BIGNUM *);
117
118 #endif
119
120 static void BN_print_dec_fp(FILE *, const BIGNUM *);
121 static void convert_str2bn(BIGNUM **, char *);
122 static void pr_fact(BIGNUM *, int, int); /* print factors of a value */
123 static void pr_print(BIGNUM *, int, int); /* print a prime */
124 static void usage(void) __dead;
125
126 static BN_CTX *ctx; /* just use a global context */
127
128 int
main(int argc,char * argv[])129 main(int argc, char *argv[])
130 {
131 BIGNUM *val;
132 int ch, hflag, xflag;
133 char *p, buf[LINE_MAX]; /* > max number of digits. */
134
135 hflag = xflag = 0;
136 ctx = BN_CTX_new();
137 val = BN_new();
138 if (val == NULL)
139 errx(1, "can't initialise bignum");
140
141 while ((ch = getopt(argc, argv, "hx")) != -1)
142 switch (ch) {
143 case 'h':
144 hflag++;
145 break;
146 case 'x':
147 xflag++;
148 break;
149 case '?':
150 default:
151 usage();
152 }
153 argc -= optind;
154 argv += optind;
155
156 /* No args supplied, read numbers from stdin. */
157 if (argc == 0)
158 for (;;) {
159 if (fgets(buf, sizeof(buf), stdin) == NULL) {
160 if (ferror(stdin))
161 err(1, "stdin");
162 exit (0);
163 }
164 for (p = buf; isblank((unsigned char)*p); ++p);
165 if (*p == '\n' || *p == '\0')
166 continue;
167 convert_str2bn(&val, p);
168 pr_fact(val, hflag, xflag);
169 }
170 /* Factor the arguments. */
171 else
172 for (p = *argv; p != NULL; p = *++argv) {
173 convert_str2bn(&val, p);
174 pr_fact(val, hflag, xflag);
175 }
176 exit(0);
177 }
178
179 static void
pr_exp(int i,int xflag)180 pr_exp(int i, int xflag)
181 {
182 printf(xflag && i > 9 ? "^0x%x" : "^%d", i);
183 }
184
185 static void
pr_uint64(uint64_t i,int xflag)186 pr_uint64(uint64_t i, int xflag)
187 {
188 printf(xflag ? " 0x%" PRIx64 : " %" PRIu64, i);
189 }
190
191 /*
192 * pr_fact - print the factors of a number
193 *
194 * Print the factors of the number, from the lowest to the highest.
195 * A factor will be printed multiple times if it divides the value
196 * multiple times.
197 *
198 * If hflag is specified, duplicate factors are printed in "human-
199 * readable" form of x^n.
200 *
201 * If the xflag is specified, numbers will be printed in hex.
202 *
203 * Factors are printed with leading tabs.
204 */
205 static void
pr_fact(BIGNUM * val,int hflag,int xflag)206 pr_fact(BIGNUM *val, int hflag, int xflag)
207 {
208 int i;
209 const uint64_t *fact; /* The factor found. */
210
211 /* Firewall - catch 0 and 1. */
212 if (BN_is_zero(val)) /* Historical practice; 0 just exits. */
213 exit(0);
214 if (BN_is_one(val)) {
215 printf("1: 1\n");
216 return;
217 }
218
219 /* Factor value. */
220
221 if (xflag) {
222 fputs("0x", stdout);
223 BN_print_fp(stdout, val);
224 } else
225 BN_print_dec_fp(stdout, val);
226 putchar(':');
227 fflush(stdout);
228 for (fact = &prime[0]; !BN_is_one(val); ++fact) {
229 /* Look for the smallest factor. */
230 do {
231 if (BN_mod_word(val, (BN_ULONG)*fact) == 0)
232 break;
233 } while (++fact <= pr_limit);
234
235 /* Watch for primes larger than the table. */
236 if (fact > pr_limit) {
237 #ifdef HAVE_OPENSSL
238 BIGNUM *bnfact;
239
240 bnfact = BN_new();
241 BN_set_word(bnfact, *(fact - 1));
242 if (!BN_sqr(bnfact, bnfact, ctx))
243 errx(1, "error in BN_sqr()");
244 if (BN_cmp(bnfact, val) > 0 ||
245 BN_is_prime_ex(val, PRIME_CHECKS, NULL, NULL) == 1)
246 pr_print(val, hflag, xflag);
247 else
248 pollard_pminus1(val, hflag, xflag);
249 #else
250 pr_print(val, hflag, xflag);
251 #endif
252 pr_print(NULL, hflag, xflag);
253 break;
254 }
255
256 /* Divide factor out until none are left. */
257 i = 0;
258 do {
259 i++;
260 if (!hflag)
261 pr_uint64(*fact, xflag);
262 BN_div_word(val, (BN_ULONG)*fact);
263 } while (BN_mod_word(val, (BN_ULONG)*fact) == 0);
264
265 if (hflag) {
266 pr_uint64(*fact, xflag);
267 if (i > 1)
268 pr_exp(i, xflag);
269 }
270
271 /* Let the user know we're doing something. */
272 fflush(stdout);
273 }
274 putchar('\n');
275 }
276
277 static void
pr_print(BIGNUM * val,int hflag,int xflag)278 pr_print(BIGNUM *val, int hflag, int xflag)
279 {
280 static BIGNUM *sval;
281 static int ex = 1;
282 BIGNUM *pval;
283
284 if (hflag) {
285 if (sval == NULL) {
286 sval = BN_dup(val);
287 return;
288 }
289
290 if (val != NULL && BN_cmp(val, sval) == 0) {
291 ex++;
292 return;
293 }
294 pval = sval;
295 } else if (val == NULL) {
296 return;
297 } else {
298 pval = val;
299 }
300
301 if (xflag) {
302 fputs(" 0x", stdout);
303 BN_print_fp(stdout, pval);
304 } else {
305 putchar(' ');
306 BN_print_dec_fp(stdout, pval);
307 }
308
309 if (hflag) {
310 if (ex > 1)
311 pr_exp(ex, xflag);
312
313 if (val != NULL) {
314 BN_copy(sval, val);
315 } else {
316 BN_free(sval);
317 sval = NULL;
318 }
319 ex = 1;
320 }
321 }
322
323 static void
usage(void)324 usage(void)
325 {
326 fprintf(stderr, "Usage: %s [-hx] [value ...]\n", getprogname());
327 exit(1);
328 }
329
330 #ifdef HAVE_OPENSSL
331
332 /* pollard p-1, algorithm from Jim Gillogly, May 2000 */
333 static void
pollard_pminus1(BIGNUM * val,int hflag,int xflag)334 pollard_pminus1(BIGNUM *val, int hflag, int xflag)
335 {
336 BIGNUM *base, *rbase, *num, *i, *x;
337
338 base = BN_new();
339 rbase = BN_new();
340 num = BN_new();
341 i = BN_new();
342 x = BN_new();
343
344 BN_set_word(rbase, 1);
345 newbase:
346 if (!BN_add_word(rbase, 1))
347 errx(1, "error in BN_add_word()");
348 BN_set_word(i, 2);
349 BN_copy(base, rbase);
350
351 for (;;) {
352 BN_mod_exp(base, base, i, val, ctx);
353 if (BN_is_one(base))
354 goto newbase;
355
356 BN_copy(x, base);
357 BN_sub_word(x, 1);
358 if (!BN_gcd(x, x, val, ctx))
359 errx(1, "error in BN_gcd()");
360
361 if (!BN_is_one(x)) {
362 if (BN_is_prime_ex(x, PRIME_CHECKS, NULL, NULL) == 1)
363 pr_print(x, hflag, xflag);
364 else
365 pollard_pminus1(x, hflag, xflag);
366 fflush(stdout);
367
368 BN_div(num, NULL, val, x, ctx);
369 if (BN_is_one(num))
370 return;
371 if (BN_is_prime_ex(num, PRIME_CHECKS, NULL,
372 NULL) == 1) {
373 pr_print(num, hflag, xflag);
374 fflush(stdout);
375 return;
376 }
377 BN_copy(val, num);
378 }
379 if (!BN_add_word(i, 1))
380 errx(1, "error in BN_add_word()");
381 }
382 }
383
384 /*
385 * Sigh.. No _decimal_ output to file functions in BN.
386 */
387 static void
BN_print_dec_fp(FILE * fp,const BIGNUM * num)388 BN_print_dec_fp(FILE *fp, const BIGNUM *num)
389 {
390 char *buf;
391
392 buf = BN_bn2dec(num);
393 if (buf == NULL)
394 return; /* XXX do anything here? */
395 fprintf(fp, "%s", buf);
396 free(buf);
397 }
398
399 #else
400
401 static void
BN_print_fp(FILE * fp,const BIGNUM * num)402 BN_print_fp(FILE *fp, const BIGNUM *num)
403 {
404 fprintf(fp, "%lx", (unsigned long)*num);
405 }
406
407 static void
BN_print_dec_fp(FILE * fp,const BIGNUM * num)408 BN_print_dec_fp(FILE *fp, const BIGNUM *num)
409 {
410 fprintf(fp, "%lu", (unsigned long)*num);
411 }
412
413 static int
BN_dec2bn(BIGNUM ** a,const char * str)414 BN_dec2bn(BIGNUM **a, const char *str)
415 {
416 char *p;
417
418 errno = 0;
419 **a = strtoul(str, &p, 10);
420 return (errno == 0 ? 1 : 0); /* OpenSSL returns 0 on error! */
421 }
422
423 static int
BN_hex2bn(BIGNUM ** a,const char * str)424 BN_hex2bn(BIGNUM **a, const char *str)
425 {
426 char *p;
427
428 errno = 0;
429 **a = strtoul(str, &p, 16);
430 return (errno == 0 ? 1 : 0); /* OpenSSL returns 0 on error! */
431 }
432
433 static BN_ULONG
BN_div_word(BIGNUM * a,BN_ULONG b)434 BN_div_word(BIGNUM *a, BN_ULONG b)
435 {
436 BN_ULONG mod;
437
438 mod = *a % b;
439 *a /= b;
440 return mod;
441 }
442
443 static BIGNUM *
BN_dup(const BIGNUM * a)444 BN_dup(const BIGNUM *a)
445 {
446 BIGNUM *b = BN_new();
447 BN_copy(b, a);
448 return b;
449 }
450
451 #endif
452
453 /* Convert string pointed to by *str to a bignum. */
454 static void
convert_str2bn(BIGNUM ** val,char * p)455 convert_str2bn(BIGNUM **val, char *p)
456 {
457 int n = 0;
458
459 if (*p == '+') p++;
460 if (*p == '-')
461 errx(1, "negative numbers aren't permitted.");
462 if (*p == '0' && (p[1] == 'x' || p[1] == 'X')) {
463 n = BN_hex2bn(val, p + 2);
464 } else {
465 n = BN_dec2bn(val, p);
466 }
467 if (n == 0)
468 errx(1, "%s: illegal numeric format.", p);
469 }
470