1 /* $NetBSD: factor.c,v 1.27 2014/10/02 21:36:37 ast Exp $ */
2
3 /*
4 * Copyright (c) 1989, 1993
5 * The Regents of the University of California. All rights reserved.
6 *
7 * This code is derived from software contributed to Berkeley by
8 * Landon Curt Noll.
9 *
10 * Redistribution and use in source and binary forms, with or without
11 * modification, are permitted provided that the following conditions
12 * are met:
13 * 1. Redistributions of source code must retain the above copyright
14 * notice, this list of conditions and the following disclaimer.
15 * 2. Redistributions in binary form must reproduce the above copyright
16 * notice, this list of conditions and the following disclaimer in the
17 * documentation and/or other materials provided with the distribution.
18 * 3. Neither the name of the University nor the names of its contributors
19 * may be used to endorse or promote products derived from this software
20 * without specific prior written permission.
21 *
22 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
23 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
24 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
25 * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
26 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
27 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
28 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
29 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
30 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
31 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
32 * SUCH DAMAGE.
33 */
34
35 #include <sys/cdefs.h>
36 #ifndef lint
37 __COPYRIGHT("@(#) Copyright (c) 1989, 1993\
38 The Regents of the University of California. All rights reserved.");
39 #endif /* not lint */
40
41 #ifndef lint
42 #if 0
43 static char sccsid[] = "@(#)factor.c 8.4 (Berkeley) 5/4/95";
44 #else
45 __RCSID("$NetBSD: factor.c,v 1.27 2014/10/02 21:36:37 ast Exp $");
46 #endif
47 #endif /* not lint */
48
49 /*
50 * factor - factor a number into primes
51 *
52 * By Landon Curt Noll, http://www.isthe.com/chongo/index.html /\oo/\
53 *
54 * usage:
55 * factor [number] ...
56 *
57 * The form of the output is:
58 *
59 * number: factor1 factor1 factor2 factor3 factor3 factor3 ...
60 *
61 * where factor1 <= factor2 <= factor3 <= ...
62 *
63 * If no args are given, the list of numbers are read from stdin.
64 */
65
66 #include <ctype.h>
67 #include <err.h>
68 #include <errno.h>
69 #include <limits.h>
70 #include <stdio.h>
71 #include <stdlib.h>
72 #include <unistd.h>
73 #include <inttypes.h>
74
75 #ifdef HAVE_OPENSSL
76 #include <openssl/bn.h>
77 #else
78 typedef long BIGNUM;
79 typedef u_long BN_ULONG;
80 static int BN_dec2bn(BIGNUM **a, const char *str);
81 #define BN_new() ((BIGNUM *)calloc(sizeof(BIGNUM), 1))
82 #define BN_is_zero(v) (*(v) == 0)
83 #define BN_is_one(v) (*(v) == 1)
84 #define BN_mod_word(a, b) (*(a) % (b))
85 #endif
86
87 #include "primes.h"
88
89 /*
90 * prime[i] is the (i-1)th prime.
91 *
92 * We are able to sieve 2^32-1 because this byte table yields all primes
93 * up to 65537 and 65537^2 > 2^32-1.
94 */
95
96 #if 0 /* debugging: limit table use to stress the "pollard" code */
97 #define pr_limit &prime[0]
98 #endif
99
100 #define PRIME_CHECKS 5
101
102 #ifdef HAVE_OPENSSL
103 static BN_CTX *ctx; /* just use a global context */
104 #endif
105
106 static void pr_fact(BIGNUM *); /* print factors of a value */
107 static void BN_print_dec_fp(FILE *, const BIGNUM *);
108 static void usage(void) __dead;
109 #ifdef HAVE_OPENSSL
110 static void pollard_rho(BIGNUM *); /* print factors for big numbers */
111 #else
112 static char *BN_bn2dec(const BIGNUM *);
113 static BN_ULONG BN_div_word(BIGNUM *, BN_ULONG);
114 #endif
115
116
117 #ifndef HAVE_OPENSSL
118 static int
BN_dec2bn(BIGNUM ** a,const char * str)119 BN_dec2bn(BIGNUM **a, const char *str)
120 {
121 char *p;
122
123 errno = 0;
124 **a = strtoul(str, &p, 10);
125 if (errno)
126 err(1, "%s", str);
127 return (*p == '\n' || *p == '\0');
128 }
129 #endif
130
131 int
main(int argc,char * argv[])132 main(int argc, char *argv[])
133 {
134 BIGNUM *val;
135 int ch;
136 char *p, buf[LINE_MAX]; /* > max number of digits. */
137
138 #ifdef HAVE_OPENSSL
139 ctx = BN_CTX_new();
140 #endif
141 val = BN_new();
142 if (val == NULL)
143 errx(1, "can't initialise bignum");
144
145 while ((ch = getopt(argc, argv, "")) != -1)
146 switch (ch) {
147 case '?':
148 default:
149 usage();
150 }
151 argc -= optind;
152 argv += optind;
153
154 /* No args supplied, read numbers from stdin. */
155 if (argc == 0)
156 for (;;) {
157 if (fgets(buf, sizeof(buf), stdin) == NULL) {
158 if (ferror(stdin))
159 err(1, "stdin");
160 exit (0);
161 }
162 for (p = buf; isblank((unsigned char)*p); ++p);
163 if (*p == '\n' || *p == '\0')
164 continue;
165 if (*p == '-')
166 errx(1, "negative numbers aren't permitted.");
167 if (BN_dec2bn(&val, buf) == 0)
168 errx(1, "%s: illegal numeric format.", argv[0]);
169 pr_fact(val);
170 }
171 /* Factor the arguments. */
172 else
173 for (; *argv != NULL; ++argv) {
174 if (argv[0][0] == '-')
175 errx(1, "numbers <= 1 aren't permitted.");
176 if (BN_dec2bn(&val, argv[0]) == 0)
177 errx(1, "%s: illegal numeric format.", argv[0]);
178 pr_fact(val);
179 }
180 exit(0);
181 }
182
183 /*
184 * pr_fact - print the factors of a number
185 *
186 * If the number is 0 or 1, then print the number and return.
187 * If the number is < 0, print -1, negate the number and continue
188 * processing.
189 *
190 * Print the factors of the number, from the lowest to the highest.
191 * A factor will be printed numtiple times if it divides the value
192 * multiple times.
193 *
194 * Factors are printed with leading tabs.
195 */
196 static void
pr_fact(BIGNUM * val)197 pr_fact(BIGNUM *val)
198 {
199 const uint64_t *fact; /* The factor found. */
200
201 /* Firewall - catch 0 and 1. */
202 if (BN_is_zero(val) || BN_is_one(val))
203 errx(1, "numbers <= 1 aren't permitted.");
204
205 /* Factor value. */
206
207 BN_print_dec_fp(stdout, val);
208 putchar(':');
209 for (fact = &prime[0]; !BN_is_one(val); ++fact) {
210 /* Look for the smallest factor. */
211 while (fact <= pr_limit) {
212 if (BN_mod_word(val, (BN_ULONG)*fact) == 0)
213 break;
214 fact++;
215 }
216
217 /* Watch for primes larger than the table. */
218 if (fact > pr_limit) {
219 #ifdef HAVE_OPENSSL
220 BIGNUM *bnfact;
221
222 bnfact = BN_new();
223 BN_set_word(bnfact, (BN_ULONG)*(fact - 1));
224 BN_sqr(bnfact, bnfact, ctx);
225 if (BN_cmp(bnfact, val) > 0
226 || BN_is_prime(val, PRIME_CHECKS, NULL, NULL,
227 NULL) == 1) {
228 putchar(' ');
229 BN_print_dec_fp(stdout, val);
230 } else
231 pollard_rho(val);
232 #else
233 printf(" %s", BN_bn2dec(val));
234 #endif
235 break;
236 }
237
238 /* Divide factor out until none are left. */
239 do {
240 printf(" %" PRIu64, *fact);
241 BN_div_word(val, (BN_ULONG)*fact);
242 } while (BN_mod_word(val, (BN_ULONG)*fact) == 0);
243
244 /* Let the user know we're doing something. */
245 fflush(stdout);
246 }
247 putchar('\n');
248 }
249
250 /*
251 * Sigh.. No _decimal_ output to file functions in BN.
252 */
253 static void
BN_print_dec_fp(FILE * fp,const BIGNUM * num)254 BN_print_dec_fp(FILE *fp, const BIGNUM *num)
255 {
256 char *buf;
257
258 buf = BN_bn2dec(num);
259 if (buf == NULL)
260 return; /* XXX do anything here? */
261 fprintf(fp, "%s", buf);
262 free(buf);
263 }
264
265 void
usage(void)266 usage(void)
267 {
268 fprintf(stderr, "usage: factor [value ...]\n");
269 exit (0);
270 }
271
272
273
274
275 #ifdef HAVE_OPENSSL
276 static void
pollard_rho(BIGNUM * val)277 pollard_rho(BIGNUM *val)
278 {
279 BIGNUM *x, *y, *tmp, *num;
280 BN_ULONG a;
281 unsigned int steps_taken, steps_limit;
282
283 x = BN_new();
284 y = BN_new();
285 tmp = BN_new();
286 num = BN_new();
287 a = 1;
288 restart:
289 steps_taken = 0;
290 steps_limit = 2;
291 BN_set_word(x, 1);
292 BN_copy(y, x);
293
294 for (;;) {
295 BN_sqr(tmp, x, ctx);
296 BN_add_word(tmp, a);
297 BN_mod(x, tmp, val, ctx);
298 BN_sub(tmp, x, y);
299 if (BN_is_zero(tmp)) {
300 #ifdef DEBUG
301 printf(" (loop)");
302 #endif
303 a++;
304 goto restart;
305 }
306 BN_gcd(tmp, tmp, val, ctx);
307
308 if (!BN_is_one(tmp)) {
309 if (BN_is_prime(tmp, PRIME_CHECKS, NULL, NULL,
310 NULL) == 1) {
311 putchar(' ');
312 BN_print_dec_fp(stdout, tmp);
313 } else {
314 #ifdef DEBUG
315 printf(" (recurse for ");
316 BN_print_dec_fp(stdout, tmp);
317 putchar(')');
318 #endif
319 pollard_rho(BN_dup(tmp));
320 #ifdef DEBUG
321 printf(" (back)");
322 #endif
323 }
324 fflush(stdout);
325
326 BN_div(num, NULL, val, tmp, ctx);
327 if (BN_is_one(num))
328 return;
329 if (BN_is_prime(num, PRIME_CHECKS, NULL, NULL,
330 NULL) == 1) {
331 putchar(' ');
332 BN_print_dec_fp(stdout, num);
333 fflush(stdout);
334 return;
335 }
336 BN_copy(val, num);
337 goto restart;
338 }
339 steps_taken++;
340 if (steps_taken == steps_limit) {
341 BN_copy(y, x); /* teleport the turtle */
342 steps_taken = 0;
343 steps_limit *= 2;
344 if (steps_limit == 0) {
345 #ifdef DEBUG
346 printf(" (overflow)");
347 #endif
348 a++;
349 goto restart;
350 }
351 }
352 }
353 }
354 #else
355 char *
BN_bn2dec(const BIGNUM * val)356 BN_bn2dec(const BIGNUM *val)
357 {
358 char *buf;
359
360 buf = malloc(100);
361 if (!buf)
362 return buf;
363 snprintf(buf, 100, "%ld", (long)*val);
364 return buf;
365 }
366
367 static BN_ULONG
BN_div_word(BIGNUM * a,BN_ULONG b)368 BN_div_word(BIGNUM *a, BN_ULONG b)
369 {
370 BN_ULONG mod;
371
372 mod = *a % b;
373 *a /= b;
374 return mod;
375 }
376 #endif
377