]> git.cameronkatri.com Git - bsdgames-darwin.git/blobdiff - factor/factor.c
Follow the Fundamental Theory of Algebra. Disallow factorising of
[bsdgames-darwin.git] / factor / factor.c
index 23fce2b7d198d31d46f5ff1b45154aa9079581e7..10de005d0fa3efe21b46fc487578c3eb95613692 100644 (file)
@@ -1,4 +1,4 @@
-/*     $NetBSD: factor.c,v 1.13 2002/06/18 23:07:36 simonb Exp $       */
+/*     $NetBSD: factor.c,v 1.24 2010/05/15 21:22:39 joerg Exp $        */
 
 /*
  * Copyright (c) 1989, 1993
  * 2. Redistributions in binary form must reproduce the above copyright
  *    notice, this list of conditions and the following disclaimer in the
  *    documentation and/or other materials provided with the distribution.
- * 3. All advertising materials mentioning features or use of this software
- *    must display the following acknowledgement:
- *     This product includes software developed by the University of
- *     California, Berkeley and its contributors.
- * 4. Neither the name of the University nor the names of its contributors
+ * 3. Neither the name of the University nor the names of its contributors
  *    may be used to endorse or promote products derived from this software
  *    without specific prior written permission.
  *
 
 #include <sys/cdefs.h>
 #ifndef lint
-__COPYRIGHT("@(#) Copyright (c) 1989, 1993\n\
      The Regents of the University of California.  All rights reserved.\n");
+__COPYRIGHT("@(#) Copyright (c) 1989, 1993\
The Regents of the University of California.  All rights reserved.");
 #endif /* not lint */
 
 #ifndef lint
 #if 0
 static char sccsid[] = "@(#)factor.c   8.4 (Berkeley) 5/4/95";
 #else
-__RCSID("$NetBSD: factor.c,v 1.13 2002/06/18 23:07:36 simonb Exp $");
+__RCSID("$NetBSD: factor.c,v 1.24 2010/05/15 21:22:39 joerg Exp $");
 #endif
 #endif /* not lint */
 
@@ -64,7 +60,7 @@ __RCSID("$NetBSD: factor.c,v 1.13 2002/06/18 23:07:36 simonb Exp $");
  *
  *     number: factor1 factor1 factor2 factor3 factor3 factor3 ...
  *
- * where factor1 < factor2 < factor3 < ...
+ * where factor1 <= factor2 <= factor3 <= ...
  *
  * If no args are given, the list of numbers are read from stdin.
  */
@@ -82,7 +78,7 @@ __RCSID("$NetBSD: factor.c,v 1.13 2002/06/18 23:07:36 simonb Exp $");
 #else
 typedef long   BIGNUM;
 typedef u_long BN_ULONG;
-int    BN_dec2bn(BIGNUM **a, const char *str);
+static int BN_dec2bn(BIGNUM **a, const char *str);
 #define BN_new()               ((BIGNUM *)calloc(sizeof(BIGNUM), 1))
 #define BN_is_zero(v)          (*(v) == 0)
 #define BN_is_one(v)           (*(v) == 1)
@@ -102,27 +98,29 @@ int        BN_dec2bn(BIGNUM **a, const char *str);
  */
 extern const ubig prime[];
 extern const ubig *pr_limit;           /* largest prime in the prime array */
+#if 0 /* debugging: limit table use to stress the "pollard" code */
+#define pr_limit &prime[0]
+#endif
 
 #define        PRIME_CHECKS    5
 
 #ifdef HAVE_OPENSSL 
-BN_CTX *ctx;                           /* just use a global context */
+static BN_CTX *ctx;                    /* just use a global context */
 #endif
 
-int    main(int, char *[]);
-void   pr_fact(BIGNUM *);              /* print factors of a value */
-void   BN_print_dec_fp(FILE *, const BIGNUM *);
-void   usage(void) __attribute__((__noreturn__));
+static void pr_fact(BIGNUM *);         /* print factors of a value */
+static void BN_print_dec_fp(FILE *, const BIGNUM *);
+static void usage(void) __dead;
 #ifdef HAVE_OPENSSL
-void   pollard_pminus1(BIGNUM *);      /* print factors for big numbers */
+static void pollard_rho(BIGNUM *);     /* print factors for big numbers */
 #else
-char   *BN_bn2dec(const BIGNUM *);
-BN_ULONG BN_div_word(BIGNUM *, BN_ULONG);
+static char *BN_bn2dec(const BIGNUM *);
+static BN_ULONG BN_div_word(BIGNUM *, BN_ULONG);
 #endif
 
 
 #ifndef HAVE_OPENSSL
-int
+static int
 BN_dec2bn(BIGNUM **a, const char *str)
 {
        char *p;
@@ -166,7 +164,7 @@ main(int argc, char *argv[])
                                        err(1, "stdin");
                                exit (0);
                        }
-                       for (p = buf; isblank(*p); ++p);
+                       for (p = buf; isblank((unsigned char)*p); ++p);
                        if (*p == '\n' || *p == '\0')
                                continue;
                        if (*p == '-')
@@ -179,7 +177,7 @@ main(int argc, char *argv[])
        else
                for (; *argv != NULL; ++argv) {
                        if (argv[0][0] == '-')
-                               errx(1, "negative numbers aren't permitted.");
+                               errx(1, "numbers <= 1 aren't permitted.");
                        if (BN_dec2bn(&val, argv[0]) == 0)
                                errx(1, "%s: illegal numeric format.", argv[0]);
                        pr_fact(val);
@@ -200,18 +198,14 @@ main(int argc, char *argv[])
  *
  * Factors are printed with leading tabs.
  */
-void
+static void
 pr_fact(BIGNUM *val)
 {
        const ubig *fact;               /* The factor found. */
 
        /* Firewall - catch 0 and 1. */
-       if (BN_is_zero(val))    /* Historical practice; 0 just exits. */
-               exit(0);
-       if (BN_is_one(val)) {
-               printf("1: 1\n");
-               return;
-       }
+       if (BN_is_zero(val) || BN_is_one(val))
+               errx(1, "numbers <= 1 aren't permitted.");
 
        /* Factor value. */
 
@@ -219,10 +213,11 @@ pr_fact(BIGNUM *val)
        putchar(':');
        for (fact = &prime[0]; !BN_is_one(val); ++fact) {
                /* Look for the smallest factor. */
-               do {
+               while (fact <= pr_limit) {
                        if (BN_mod_word(val, (BN_ULONG)*fact) == 0)
                                break;
-               } while (++fact <= pr_limit);
+                       fact++;
+               }
 
                /* Watch for primes larger than the table. */
                if (fact > pr_limit) {
@@ -230,13 +225,15 @@ pr_fact(BIGNUM *val)
                        BIGNUM *bnfact;
 
                        bnfact = BN_new();
-                       BN_set_word(bnfact, *(fact - 1));
+                       BN_set_word(bnfact, (BN_ULONG)*(fact - 1));
                        BN_sqr(bnfact, bnfact, ctx);
-                       if (BN_cmp(bnfact, val) > 0) {
+                       if (BN_cmp(bnfact, val) > 0
+                           || BN_is_prime(val, PRIME_CHECKS, NULL, NULL,
+                                          NULL) == 1) {
                                putchar(' ');
                                BN_print_dec_fp(stdout, val);
                        } else
-                               pollard_pminus1(val);
+                               pollard_rho(val);
 #else
                        printf(" %s", BN_bn2dec(val));
 #endif
@@ -258,7 +255,7 @@ pr_fact(BIGNUM *val)
 /*
  * Sigh..  No _decimal_ output to file functions in BN.
  */
-void
+static void
 BN_print_dec_fp(FILE *fp, const BIGNUM *num)
 {
        char *buf;
@@ -281,38 +278,57 @@ usage(void)
 
 
 #ifdef HAVE_OPENSSL
-/* pollard rho, algorithm from Jim Gillogly, May 2000 */
-
-void
-pollard_pminus1(BIGNUM *val)
+static void
+pollard_rho(BIGNUM *val)
 {
-       BIGNUM *base, *num, *i, *x;
+       BIGNUM *x, *y, *tmp, *num;
+       BN_ULONG a;
+       unsigned int steps_taken, steps_limit;
 
-       base = BN_new();
-       num = BN_new();
-       i = BN_new();
        x = BN_new();
-
-       BN_set_word(i, 2);
-       BN_set_word(base, 2);
+       y = BN_new();
+       tmp = BN_new();
+       num = BN_new();
+       a = 1;
+restart:
+       steps_taken = 0;
+       steps_limit = 2;
+       BN_set_word(x, 1);
+       BN_copy(y, x);
 
        for (;;) {
-               BN_mod_exp(base, base, i, val, ctx);
-
-               BN_copy(x, base);
-               BN_sub_word(x, 1);
-               BN_gcd(x, x, val, ctx);
+               BN_sqr(tmp, x, ctx);
+               BN_add_word(tmp, a);
+               BN_mod(x, tmp, val, ctx);
+               BN_sub(tmp, x, y);
+               if (BN_is_zero(tmp)) {
+#ifdef DEBUG
+                       printf(" (loop)");
+#endif
+                       a++;
+                       goto restart;
+               }
+               BN_gcd(tmp, tmp, val, ctx);
 
-               if (!BN_is_one(x)) {
-                       if (BN_is_prime(x, PRIME_CHECKS, NULL, NULL,
+               if (!BN_is_one(tmp)) {
+                       if (BN_is_prime(tmp, PRIME_CHECKS, NULL, NULL,
                            NULL) == 1) {
                                putchar(' ');
-                               BN_print_dec_fp(stdout, x);
-                       } else
-                               pollard_pminus1(x);
+                               BN_print_dec_fp(stdout, tmp);
+                       } else {
+#ifdef DEBUG
+                               printf(" (recurse for ");
+                               BN_print_dec_fp(stdout, tmp);
+                               putchar(')');
+#endif
+                               pollard_rho(BN_dup(tmp));
+#ifdef DEBUG
+                               printf(" (back)");
+#endif
+                       }
                        fflush(stdout);
 
-                       BN_div(num, NULL, val, x, ctx);
+                       BN_div(num, NULL, val, tmp, ctx);
                        if (BN_is_one(num))
                                return;
                        if (BN_is_prime(num, PRIME_CHECKS, NULL, NULL,
@@ -323,8 +339,21 @@ pollard_pminus1(BIGNUM *val)
                                return;
                        }
                        BN_copy(val, num);
+                       goto restart;
+               }
+               steps_taken++;
+               if (steps_taken == steps_limit) {
+                       BN_copy(y, x); /* teleport the turtle */
+                       steps_taken = 0;
+                       steps_limit *= 2;
+                       if (steps_limit == 0) {
+#ifdef DEBUG
+                               printf(" (overflow)");
+#endif
+                               a++;
+                               goto restart;
+                       }
                }
-               BN_add_word(i, 1);
        }
 }
 #else
@@ -340,7 +369,7 @@ BN_bn2dec(const BIGNUM *val)
        return buf;
 }
 
-BN_ULONG
+static BN_ULONG
 BN_div_word(BIGNUM *a, BN_ULONG b)
 {
        BN_ULONG mod;