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