Initial import from FreeBSD RELENG_4:
[dragonfly.git] / games / factor / factor.c
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.
35  */
36
37 #ifndef lint
38 static const char copyright[] =
39 "@(#) Copyright (c) 1989, 1993\n\
40         The Regents of the University of California.  All rights reserved.\n";
41 #endif /* not lint */
42
43 #ifndef lint
44 #if 0
45 static char sccsid[] = "@(#)factor.c    8.4 (Berkeley) 5/4/95";
46 __RCSID("$NetBSD: factor.c,v 1.13 2002/06/18 23:07:36 simonb Exp $");
47 #endif
48 static const char rcsid[] =
49  "$FreeBSD: src/games/factor/factor.c,v 1.9.2.2 2002/10/23 14:59:14 fanf Exp $";
50 #endif /* not lint */
51
52 /*
53  * factor - factor a number into primes
54  *
55  * By: Landon Curt Noll   chongo@toad.com,   ...!{sun,tolsoft}!hoptoad!chongo
56  *
57  *   chongo <for a good prime call: 391581 * 2^216193 - 1> /\oo/\
58  *
59  * usage:
60  *      factor [-h] [number] ...
61  *
62  * The form of the output is:
63  *
64  *      number: factor1 factor1 factor2 factor3 factor3 factor3 ...
65  *
66  * where factor1 < factor2 < factor3 < ...
67  *
68  * If no args are given, the list of numbers are read from stdin.
69  */
70
71 #include <ctype.h>
72 #include <err.h>
73 #include <errno.h>
74 #include <limits.h>
75 #include <stdio.h>
76 #include <stdlib.h>
77 #include <unistd.h>
78
79 #include "primes.h"
80
81 #ifdef HAVE_OPENSSL
82
83 #include <openssl/bn.h>
84
85 #define PRIME_CHECKS    5
86
87 static void     pollard_pminus1(BIGNUM *); /* print factors for big numbers */
88
89 #else
90
91 typedef ubig    BIGNUM;
92 typedef u_long  BN_ULONG;
93
94 #define BN_CTX                  int
95 #define BN_CTX_new()            NULL
96 #define BN_new()                ((BIGNUM *)calloc(sizeof(BIGNUM), 1))
97 #define BN_is_zero(v)           (*(v) == 0)
98 #define BN_is_one(v)            (*(v) == 1)
99 #define BN_mod_word(a, b)       (*(a) % (b))
100
101 static int      BN_dec2bn(BIGNUM **a, const char *str);
102 static int      BN_hex2bn(BIGNUM **a, const char *str);
103 static BN_ULONG BN_div_word(BIGNUM *, BN_ULONG);
104 static void     BN_print_fp(FILE *, const BIGNUM *);
105
106 #endif
107
108 static void     BN_print_dec_fp(FILE *, const BIGNUM *);
109
110 static void     pr_fact(BIGNUM *);      /* print factors of a value */
111 static void     pr_print(BIGNUM *);     /* print a prime */
112 static void     usage(void);
113
114 static BN_CTX   *ctx;                   /* just use a global context */
115 static int      hflag;
116
117 int
118 main(int argc, char *argv[])
119 {
120         BIGNUM *val;
121         int ch;
122         char *p, buf[LINE_MAX];         /* > max number of digits. */
123
124         ctx = BN_CTX_new();
125         val = BN_new();
126         if (val == NULL)
127                 errx(1, "can't initialise bignum");
128
129         while ((ch = getopt(argc, argv, "h")) != -1)
130                 switch (ch) {
131                 case 'h':
132                         hflag++;
133                         break;
134                 case '?':
135                 default:
136                         usage();
137                 }
138         argc -= optind;
139         argv += optind;
140
141         /* No args supplied, read numbers from stdin. */
142         if (argc == 0)
143                 for (;;) {
144                         if (fgets(buf, sizeof(buf), stdin) == NULL) {
145                                 if (ferror(stdin))
146                                         err(1, "stdin");
147                                 exit (0);
148                         }
149                         for (p = buf; isblank(*p); ++p);
150                         if (*p == '\n' || *p == '\0')
151                                 continue;
152                         if (*p == '-')
153                                 errx(1, "negative numbers aren't permitted.");
154                         if (BN_dec2bn(&val, buf) == 0 &&
155                             BN_hex2bn(&val, buf) == 0)
156                                 errx(1, "%s: illegal numeric format.", buf);
157                         pr_fact(val);
158                 }
159         /* Factor the arguments. */
160         else
161                 for (; *argv != NULL; ++argv) {
162                         if (argv[0][0] == '-')
163                                 errx(1, "negative numbers aren't permitted.");
164                         if (BN_dec2bn(&val, argv[0]) == 0 &&
165                             BN_hex2bn(&val, argv[0]) == 0)
166                                 errx(1, "%s: illegal numeric format.", argv[0]);
167                         pr_fact(val);
168                 }
169         exit(0);
170 }
171
172 /*
173  * pr_fact - print the factors of a number
174  *
175  * Print the factors of the number, from the lowest to the highest.
176  * A factor will be printed multiple times if it divides the value
177  * multiple times.
178  *
179  * Factors are printed with leading tabs.
180  */
181 static void
182 pr_fact(BIGNUM *val)
183 {
184         const ubig *fact;       /* The factor found. */
185
186         /* Firewall - catch 0 and 1. */
187         if (BN_is_zero(val))    /* Historical practice; 0 just exits. */
188                 exit(0);
189         if (BN_is_one(val)) {
190                 printf("1: 1\n");
191                 return;
192         }
193
194         /* Factor value. */
195
196         if (hflag) {
197                 fputs("0x", stdout);
198                 BN_print_fp(stdout, val);
199         } else
200                 BN_print_dec_fp(stdout, val);
201         putchar(':');
202         for (fact = &prime[0]; !BN_is_one(val); ++fact) {
203                 /* Look for the smallest factor. */
204                 do {
205                         if (BN_mod_word(val, (BN_ULONG)*fact) == 0)
206                                 break;
207                 } while (++fact <= pr_limit);
208
209                 /* Watch for primes larger than the table. */
210                 if (fact > pr_limit) {
211 #ifdef HAVE_OPENSSL
212                         BIGNUM *bnfact;
213
214                         bnfact = BN_new();
215                         BN_set_word(bnfact, *(fact - 1));
216                         BN_sqr(bnfact, bnfact, ctx);
217                         if (BN_cmp(bnfact, val) > 0)
218                                 pr_print(val);
219                         else
220                                 pollard_pminus1(val);
221 #else
222                         pr_print(val);
223 #endif
224                         break;
225                 }
226
227                 /* Divide factor out until none are left. */
228                 do {
229                         printf(hflag ? " 0x%lx" : " %lu", *fact);
230                         BN_div_word(val, (BN_ULONG)*fact);
231                 } while (BN_mod_word(val, (BN_ULONG)*fact) == 0);
232
233                 /* Let the user know we're doing something. */
234                 fflush(stdout);
235         }
236         putchar('\n');
237 }
238
239 static void
240 pr_print(BIGNUM *val)
241 {
242         if (hflag) {
243                 fputs(" 0x", stdout);
244                 BN_print_fp(stdout, val);
245         } else {
246                 putchar(' ');
247                 BN_print_dec_fp(stdout, val);
248         }
249 }
250
251 static void
252 usage(void)
253 {
254         fprintf(stderr, "usage: factor [-h] [value ...]\n");
255         exit(1);
256 }
257
258 #ifdef HAVE_OPENSSL
259
260 /* pollard rho, algorithm from Jim Gillogly, May 2000 */
261 static void
262 pollard_pminus1(BIGNUM *val)
263 {
264         BIGNUM *base, *num, *i, *x;
265
266         base = BN_new();
267         num = BN_new();
268         i = BN_new();
269         x = BN_new();
270
271         BN_set_word(i, 2);
272         BN_set_word(base, 2);
273
274         for (;;) {
275                 BN_mod_exp(base, base, i, val, ctx);
276
277                 BN_copy(x, base);
278                 BN_sub_word(x, 1);
279                 BN_gcd(x, x, val, ctx);
280
281                 if (!BN_is_one(x)) {
282                         if (BN_is_prime(x, PRIME_CHECKS, NULL, NULL,
283                             NULL) == 1)
284                                 pr_print(x);
285                         else
286                                 pollard_pminus1(x);
287                         fflush(stdout);
288
289                         BN_div(num, NULL, val, x, ctx);
290                         if (BN_is_one(num))
291                                 return;
292                         if (BN_is_prime(num, PRIME_CHECKS, NULL, NULL,
293                             NULL) == 1) {
294                                 pr_print(num);
295                                 fflush(stdout);
296                                 return;
297                         }
298                         BN_copy(val, num);
299                 }
300                 BN_add_word(i, 1);
301         }
302 }
303
304 /*
305  * Sigh..  No _decimal_ output to file functions in BN.
306  */
307 static void
308 BN_print_dec_fp(FILE *fp, const BIGNUM *num)
309 {
310         char *buf;
311
312         buf = BN_bn2dec(num);
313         if (buf == NULL)
314                 return; /* XXX do anything here? */
315         fprintf(fp, buf);
316         free(buf);
317 }
318
319 #else
320
321 static void
322 BN_print_fp(FILE *fp, const BIGNUM *num)
323 {
324         fprintf(fp, "%lx", (unsigned long)*num);
325 }
326
327 static void
328 BN_print_dec_fp(FILE *fp, const BIGNUM *num)
329 {
330         fprintf(fp, "%lu", (unsigned long)*num);
331 }
332
333 static int
334 BN_dec2bn(BIGNUM **a, const char *str)
335 {
336         char *p;
337
338         errno = 0;
339         **a = strtoul(str, &p, 10);
340         return (errno == 0 && (*p == '\n' || *p == '\0'));
341 }
342
343 static int
344 BN_hex2bn(BIGNUM **a, const char *str)
345 {
346         char *p;
347
348         errno = 0;
349         **a = strtoul(str, &p, 16);
350         return (errno == 0 && (*p == '\n' || *p == '\0'));
351 }
352
353 static BN_ULONG
354 BN_div_word(BIGNUM *a, BN_ULONG b)
355 {
356         BN_ULONG mod;
357
358         mod = *a % b;
359         *a /= b;
360         return mod;
361 }
362
363 #endif