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