diff options
| -rw-r--r-- | .gitignore | 6 | ||||
| -rw-r--r-- | AUTHORS | 2 | ||||
| -rw-r--r-- | cfg.mk | 23 | ||||
| -rw-r--r-- | src/.gitignore | 2 | ||||
| -rw-r--r-- | src/factor-ng.c | 2396 | ||||
| -rw-r--r-- | src/factor.c | 2463 | ||||
| -rw-r--r-- | src/local.mk | 31 | ||||
| -rw-r--r-- | src/longlong.h | 6 | ||||
| -rwxr-xr-x | src/wheel-gen.pl | 114 | ||||
| -rwxr-xr-x | tests/factor/run.sh | 30 | ||||
| -rw-r--r-- | tests/local.mk | 71 |
11 files changed, 2310 insertions, 2834 deletions
diff --git a/.gitignore b/.gitignore index 56a9bf40e..8a423f752 100644 --- a/.gitignore +++ b/.gitignore @@ -37,7 +37,7 @@ /coreutils-*.tar.xz /coreutils-*.tar.xz.sig /gnulib-tests -/lib/.cvsignore +/lib/.dirstamp /lib/.gitignore /lib/alloca.h /lib/arg-nonnull.h @@ -74,6 +74,7 @@ /lib/ref-del.sed /lib/selinux /lib/signal.h +/lib/spawn.h /lib/stamp-h1 /lib/stdalign.h /lib/stdio.h @@ -152,12 +153,15 @@ /po/remove-potcdate.sin /po/stamp-po /src/cu-progs.mk +/src/make-prime-list +/src/primes.h /src/version.c /src/version.h /stamp-h1 /tests/*/*.log /tests/*/*.trs /tests/.built-programs +/tests/factor/[0-9]*.sh /tests/t? /tests/test-suite.log ID @@ -26,7 +26,7 @@ echo: Brian Fox, Chet Ramey env: Richard Mlynarik, David MacKenzie expand: David MacKenzie expr: Mike Parker, James Youngman, Paul Eggert -factor: Paul Rubin +factor: Paul Rubin, Torbjörn Granlund, Niels Möller false: Jim Meyering fmt: Ross Paterson fold: David MacKenzie @@ -111,7 +111,7 @@ sc_tests_list_consistency: for t in $(all_tests); do echo $$t; done; \ cd $(top_srcdir); \ $(SHELL) build-aux/vc-list-files tests \ - | grep -v '^tests/init\.sh$$' \ + | grep -Ev '^tests/(factor/run|init)\.sh$$' \ | $(EGREP) "$$test_extensions_rx\$$"; \ } | sort | uniq -u | grep . && exit 1; : @@ -515,10 +515,11 @@ update-copyright-env = \ # List syntax-check exemptions. exclude_file_name_regexp--sc_space_tab = \ ^(tests/pr/|tests/misc/nl\.sh$$|gl/.*\.diff$$) -exclude_file_name_regexp--sc_bindtextdomain = ^(gl/.*|lib/euidaccess-stat)\.c$$ +exclude_file_name_regexp--sc_bindtextdomain = \ + ^(gl/.*|lib/euidaccess-stat|src/make-prime-list)\.c$$ exclude_file_name_regexp--sc_trailing_blank = ^tests/pr/ exclude_file_name_regexp--sc_system_h_headers = \ - ^src/((system|copy)\.h|libstdbuf\.c)$$ + ^src/((system|copy)\.h|libstdbuf\.c|make-prime-list\.c)$$ _src = (false|lbracket|ls-(dir|ls|vdir)|tac-pipe|uname-(arch|uname)) exclude_file_name_regexp--sc_require_config_h_first = \ @@ -530,7 +531,8 @@ exclude_file_name_regexp--sc_po_check = ^gl/ exclude_file_name_regexp--sc_prohibit_always-defined_macros = \ ^src/(seq|remove)\.c$$ exclude_file_name_regexp--sc_prohibit_empty_lines_at_EOF = ^tests/pr/ -exclude_file_name_regexp--sc_program_name = ^(gl/.*|lib/euidaccess-stat)\.c$$ +exclude_file_name_regexp--sc_program_name = \ + ^(gl/.*|lib/euidaccess-stat|src/make-prime-list)\.c$$ exclude_file_name_regexp--sc_file_system = \ NEWS|^(init\.cfg|src/df\.c|tests/df/df-P\.sh)$$ exclude_file_name_regexp--sc_prohibit_always_true_header_tests = \ @@ -539,14 +541,20 @@ exclude_file_name_regexp--sc_prohibit_fail_0 = \ (^.*/git-hooks/commit-msg|^tests/init\.sh|Makefile\.am|\.mk|.*\.texi)$$ exclude_file_name_regexp--sc_prohibit_atoi_atof = ^lib/euidaccess-stat\.c$$ +# longlong.h is maintained elsewhere. +_ll = ^src/longlong\.h$$ +exclude_file_name_regexp--sc_useless_cpp_parens = $(_ll) +exclude_file_name_regexp--sc_long_lines = $(_ll) +exclude_file_name_regexp--sc_space_before_open_paren = $(_ll) + tbi_1 = ^tests/pr/|(^gl/lib/reg.*\.c\.diff|\.mk|^man/help2man)$$ tbi_2 = ^scripts/git-hooks/(pre-commit|pre-applypatch|applypatch-msg)$$ -tbi_3 = (GNU)?[Mm]akefile(\.am)?$$ +tbi_3 = (GNU)?[Mm]akefile(\.am)?$$|$(_ll) exclude_file_name_regexp--sc_prohibit_tab_based_indentation = \ $(tbi_1)|$(tbi_2)|$(tbi_3) exclude_file_name_regexp--sc_preprocessor_indentation = \ - ^(gl/lib/rand-isaac\.[ch]|gl/tests/test-rand-isaac\.c)$$ + ^(gl/lib/rand-isaac\.[ch]|gl/tests/test-rand-isaac\.c)$$|$(__ll) exclude_file_name_regexp--sc_prohibit_stat_st_blocks = \ ^(src/system\.h|tests/du/2g\.sh)$$ @@ -560,6 +568,9 @@ exclude_file_name_regexp--sc_prohibit_test_backticks = \ exclude_file_name_regexp--sc_prohibit_operator_at_end_of_line = \ ^src/(ptx|test|head)\.c$$ +exclude_file_name_regexp--sc_error_message_uppercase = ^src/factor\.c$$ +exclude_file_name_regexp--sc_prohibit_atoi_atof = ^src/make-prime-list\.c$$ + # Augment AM_CFLAGS to include our per-directory options: AM_CFLAGS += $($(@D)_CFLAGS) diff --git a/src/.gitignore b/src/.gitignore index 1c1edbbe1..18cccc1d7 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -109,8 +109,6 @@ uptime users vdir wc -wheel-size.h -wheel.h who whoami yes diff --git a/src/factor-ng.c b/src/factor-ng.c deleted file mode 100644 index 36be5e6fe..000000000 --- a/src/factor-ng.c +++ /dev/null @@ -1,2396 +0,0 @@ -/* Factoring of uintmax_t numbers. - - Contributed to the GNU project by Torbjörn Granlund and Niels Möller - Contains code from GNU MP. - -Copyright 1995, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2005, 2009, 2012 -Free Software Foundation, Inc. - -This program is free software; you can redistribute it and/or modify it under -the terms of the GNU General Public License as published by the Free Software -Foundation; either version 3 of the License, or (at your option) any later -version. - -This program is distributed in the hope that it will be useful, but WITHOUT ANY -WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A -PARTICULAR PURPOSE. See the GNU General Public License for more details. - -You should have received a copy of the GNU General Public License along with -this program. If not, see http://www.gnu.org/licenses/. */ - -/* Efficiently factor numbers that fit in one or two words (word = uintmax_t), - or, with GMP, numbers of any size. - - Code organisation: - - There are several variants of many functions, for handling one word, two - words, and GMP's mpz_t type. If the one-word variant is called foo, the - two-word variant will be foo2, and the one for mpz_t will be mp_foo. In - some cases, the plain function variants will handle both one-word and - two-word numbers, evidenced by function arguments. - - The factoring code for two words will fall into the code for one word when - progress allows that. - - Using GMP is optional. Define HAVE_GMP to make this code include GMP - factoring code. The GMP factoring code is based on GMP's demos/factorize.c - (last synched 2012-09-07). The GMP-based factoring code will stay in GMP - factoring code even if numbers get small enough for using the two-word - code. - - Algorithm: - - (1) Perform trial division using a small primes table, but without hardware - division since the primes table store inverses modulo the word base. - (The GMP variant of this code doesn't make use of the precomputed - inverses, but instead relies on GMP for fast divisibility testing.) - (2) Check the nature of any non-factored part using Miller-Rabin for - detecting composites, and Lucas for detecting primes. - (3) Factor any remaining composite part using the Pollard-Brent rho - algorithm or the SQUFOF algorithm, checking status of found factors - again using Miller-Rabin and Lucas. - - We prefer using Hensel norm in the divisions, not the more familiar - Euclidian norm, since the former leads to much faster code. In the - Pollard-Brent rho code and the the prime testing code, we use Montgomery's - trick of multiplying all n-residues by the word base, allowing cheap Hensel - reductions mod n. - - Improvements: - - * Use modular inverses also for exact division in the Lucas code, and - elsewhere. A problem is to locate the inverses not from an index, but - from a prime. We might instead compute the inverse on-the-fly. - - * Tune trial division table size (not forgetting that this is a standalone - program where the table will be read from disk for each invocation). - - * Implement less naive powm, using k-ary exponentiation for k = 3 or - perhaps k = 4. - - * Try to speed trial division code for single uintmax_t numbers, i.e., the - code using DIVBLOCK. It currently runs at 2 cycles per prime (Intel SBR, - IBR), 3 cycles per prime (AMD Stars) and 5 cycles per prime (AMD BD) when - using gcc 4.6 and 4.7. Some software pipelining should help; 1, 2, and 4 - respectively cycles ought to be possible. - - * The redcify function could be vastly improved by using (plain Euclidian) - pre-inversion (such as GMP's invert_limb) and udiv_qrnnd_preinv (from - GMP's gmp-impl.h). The redcify2 function could be vastly improved using - similar methoods. These functions currently dominate run time when using - the -w option. -*/ - -#include <config.h> - -#include <stdlib.h> -#include <stdio.h> -#include <inttypes.h> -#include <errno.h> -#include <assert.h> -#include <ctype.h> -#include <string.h> /* for memmove */ -#include <unistd.h> /* for getopt */ -#include <stdbool.h> - -#include "system.h" - -#if HAVE_GMP -# include <gmp.h> -#endif - -#ifndef STAT_SQUFOF -# define STAT_SQUFOF 0 -#endif - -#ifndef USE_LONGLONG_H -# define USE_LONGLONG_H 1 -#endif - -#if USE_LONGLONG_H - -/* Make definitions for longlong.h to make it do what it can do for us */ -# define W_TYPE_SIZE 64 /* bitcount for uintmax_t */ -# define UWtype uintmax_t -# define UHWtype unsigned long int -# undef UDWtype -# if HAVE_ATTRIBUTE_MODE -typedef unsigned int UQItype __attribute__ ((mode (QI))); -typedef int SItype __attribute__ ((mode (SI))); -typedef unsigned int USItype __attribute__ ((mode (SI))); -typedef int DItype __attribute__ ((mode (DI))); -typedef unsigned int UDItype __attribute__ ((mode (DI))); -# else -typedef unsigned char UQItype; -typedef long SItype; -typedef unsigned long int USItype; -# if HAVE_LONG_LONG -typedef long long int DItype; -typedef unsigned long long int UDItype; -# else /* Assume `long' gives us a wide enough type. Needed for hppa2.0w. */ -typedef long int DItype; -typedef unsigned long int UDItype; -# endif -# endif -# define LONGLONG_STANDALONE /* Don't require GMP's longlong.h mdep files */ -# define ASSERT(x) /* FIXME make longlong.h really standalone */ -# define __GMP_DECLSPEC /* FIXME make longlong.h really standalone */ -# define __clz_tab factor_clz_tab /* Rename to avoid glibc collision */ -# ifndef __GMP_GNUC_PREREQ -# define __GMP_GNUC_PREREQ(a,b) 1 -# endif -# if _ARCH_PPC -# define HAVE_HOST_CPU_FAMILY_powerpc 1 -# endif -# include "longlong.h" -# ifdef COUNT_LEADING_ZEROS_NEED_CLZ_TAB -const unsigned char factor_clz_tab[129] = -{ - 1,2,3,3,4,4,4,4,5,5,5,5,5,5,5,5,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6, - 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7, - 8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8, - 8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8, - 9 -}; -# endif - -#else /* not USE_LONGLONG_H */ - -# define W_TYPE_SIZE (8 * sizeof (uintmax_t)) -# define __ll_B ((uintmax_t) 1 << (W_TYPE_SIZE / 2)) -# define __ll_lowpart(t) ((uintmax_t) (t) & (__ll_B - 1)) -# define __ll_highpart(t) ((uintmax_t) (t) >> (W_TYPE_SIZE / 2)) - -#endif - -enum alg_type { ALG_POLLARD_RHO = 1, ALG_SQUFOF = 2 }; - -static enum alg_type alg; - -/* 2*3*5*7*11...*101 is 128 bits, and has 26 prime factors */ -#define MAX_NFACTS 26 - -struct factors -{ - uintmax_t plarge[2]; /* Can have a single large factor */ - uintmax_t p[MAX_NFACTS]; - unsigned char e[MAX_NFACTS]; - unsigned char nfactors; -}; - -#if HAVE_GMP -struct mp_factors -{ - mpz_t *p; - unsigned long int *e; - unsigned long int nfactors; -}; -#endif - -static void factor (uintmax_t, uintmax_t, struct factors *); - -#ifndef umul_ppmm -# define umul_ppmm(w1, w0, u, v) \ - do { \ - uintmax_t __x0, __x1, __x2, __x3; \ - unsigned long int __ul, __vl, __uh, __vh; \ - uintmax_t __u = (u), __v = (v); \ - \ - __ul = __ll_lowpart (__u); \ - __uh = __ll_highpart (__u); \ - __vl = __ll_lowpart (__v); \ - __vh = __ll_highpart (__v); \ - \ - __x0 = (uintmax_t) __ul * __vl; \ - __x1 = (uintmax_t) __ul * __vh; \ - __x2 = (uintmax_t) __uh * __vl; \ - __x3 = (uintmax_t) __uh * __vh; \ - \ - __x1 += __ll_highpart (__x0);/* this can't give carry */ \ - __x1 += __x2; /* but this indeed can */ \ - if (__x1 < __x2) /* did we get it? */ \ - __x3 += __ll_B; /* yes, add it in the proper pos. */ \ - \ - (w1) = __x3 + __ll_highpart (__x1); \ - (w0) = (__x1 << W_TYPE_SIZE / 2) + __ll_lowpart (__x0); \ - } while (0) -#endif - -#if !defined(udiv_qrnnd) || defined(UDIV_NEEDS_NORMALIZATION) -/* Define our own, not needing normalization. This function is - currently not performance critical, so keep it simple. Similar to - the mod macro below. */ -# undef udiv_qrnnd -# define udiv_qrnnd(q, r, n1, n0, d) \ - do { \ - uintmax_t __d1, __d0, __q, __r1, __r0; \ - \ - assert ((n1) < (d)); \ - __d1 = (d); __d0 = 0; \ - __r1 = (n1); __r0 = (n0); \ - __q = 0; \ - for (unsigned int __i = W_TYPE_SIZE; __i > 0; __i--) \ - { \ - rsh2 (__d1, __d0, __d1, __d0, 1); \ - __q <<= 1; \ - if (ge2 (__r1, __r0, __d1, __d0)) \ - { \ - __q++; \ - sub_ddmmss (__r1, __r0, __r1, __r0, __d1, __d0); \ - } \ - } \ - (r) = __r0; \ - (q) = __q; \ - } while (0) -#endif - -#if !defined (add_ssaaaa) -# define add_ssaaaa(sh, sl, ah, al, bh, bl) \ - do { \ - uintmax_t _add_x; \ - _add_x = (al) + (bl); \ - (sh) = (ah) + (bh) + (_add_x < (al)); \ - (sl) = _add_x; \ - } while (0) -#endif - -#define rsh2(rh, rl, ah, al, cnt) \ - do { \ - (rl) = ((ah) << (W_TYPE_SIZE - (cnt))) | ((al) >> (cnt)); \ - (rh) = (ah) >> (cnt); \ - } while (0) - -#define lsh2(rh, rl, ah, al, cnt) \ - do { \ - (rh) = ((ah) << cnt) | ((al) >> (W_TYPE_SIZE - (cnt))); \ - (rl) = (al) << (cnt); \ - } while (0) - -#define ge2(ah, al, bh, bl) \ - ((ah) > (bh) || ((ah) == (bh) && (al) >= (bl))) - -#define gt2(ah, al, bh, bl) \ - ((ah) > (bh) || ((ah) == (bh) && (al) > (bl))) - -#ifndef sub_ddmmss -# define sub_ddmmss(rh, rl, ah, al, bh, bl) \ - do { \ - uintmax_t _cy; \ - _cy = (al) < (bl); \ - (rl) = (al) - (bl); \ - (rh) = (ah) - (bh) - _cy; \ - } while (0) -#endif - -#ifndef count_leading_zeros -# define count_leading_zeros(count, x) do { \ - uintmax_t __clz_x = (x); \ - unsigned int __clz_c; \ - for (__clz_c = 0; \ - (__clz_x & ((uintmax_t) 0xff << (W_TYPE_SIZE - 8))) == 0; \ - __clz_c += 8) \ - __clz_x <<= 8; \ - for (; (intmax_t)__clz_x >= 0; __clz_c++) \ - __clz_x <<= 1; \ - (count) = __clz_c; \ - } while (0) -#endif - -#ifndef count_trailing_zeros -# define count_trailing_zeros(count, x) do { \ - uintmax_t __ctz_x = (x); \ - unsigned int __ctz_c = 0; \ - while ((__ctz_x & 1) == 0) \ - { \ - __ctz_x >>= 1; \ - __ctz_c++; \ - } \ - (count) = __ctz_c; \ - } while (0) -#endif - -/* Requires that a < n and b <= n */ -#define submod(r,a,b,n) \ - do { \ - uintmax_t _t = - (uintmax_t) (a < b); \ - (r) = ((n) & _t) + (a) - (b); \ - } while (0) - -#define addmod(r,a,b,n) \ - submod((r), (a), ((n) - (b)), (n)) - -/* Modular two-word addition and subtraction. For performance reasons, the - most significant bit of n1 must be clear. The destination variables must be - distinct from the mod operand. */ -#define addmod2(r1, r0, a1, a0, b1, b0, n1, n0) \ - do { \ - add_ssaaaa ((r1), (r0), (a1), (a0), (b1), (b0)); \ - if (ge2 ((r1), (r0), (n1), (n0))) \ - sub_ddmmss ((r1), (r0), (r1), (r0), (n1), (n0)); \ - } while (0) -#define submod2(r1, r0, a1, a0, b1, b0, n1, n0) \ - do { \ - sub_ddmmss ((r1), (r0), (a1), (a0), (b1), (b0)); \ - if ((intmax_t) (r1) < 0) \ - add_ssaaaa ((r1), (r0), (r1), (r0), (n1), (n0)); \ - } while (0) - -#define HIGHBIT_TO_MASK(x) \ - (((intmax_t)-1 >> 1) < 0 \ - ? (uintmax_t)((intmax_t)(x) >> (W_TYPE_SIZE - 1)) \ - : ((x) & ((uintmax_t) 1 << (W_TYPE_SIZE - 1)) \ - ? UINTMAX_MAX : (uintmax_t) 0)) - -/* Compute r = a mod d, where r = <*t1,retval>, a = <a1,a0>, d = <d1,d0>. - Requires that d1 != 0. */ -static uintmax_t -mod2 (uintmax_t *r1, uintmax_t a1, uintmax_t a0, uintmax_t d1, uintmax_t d0) -{ - int cntd, cnta; - - assert (d1 != 0); - - if (a1 == 0) - { - *r1 = 0; - return a0; - } - - count_leading_zeros (cntd, d1); - count_leading_zeros (cnta, a1); - int cnt = cntd - cnta; - lsh2 (d1, d0, d1, d0, cnt); - for (int i = 0; i < cnt; i++) - { - if (ge2 (a1, a0, d1, d0)) - sub_ddmmss (a1, a0, a1, a0, d1, d0); - rsh2 (d1, d0, d1, d0, 1); - } - - *r1 = a1; - return a0; -} - -static uintmax_t _GL_ATTRIBUTE_CONST -gcd_odd (uintmax_t a, uintmax_t b) -{ - if ( (b & 1) == 0) - { - uintmax_t t = b; - b = a; - a = t; - } - if (a == 0) - return b; - - /* Take out least significant one bit, to make room for sign */ - b >>= 1; - - for (;;) - { - uintmax_t t; - uintmax_t bgta; - - while ((a & 1) == 0) - a >>= 1; - a >>= 1; - - t = a - b; - if (t == 0) - return (a << 1) + 1; - - bgta = HIGHBIT_TO_MASK (t); - - /* b <-- min (a, b) */ - b += (bgta & t); - - /* a <-- |a - b| */ - a = (t ^ bgta) - bgta; - } -} - -static uintmax_t -gcd2_odd (uintmax_t *r1, uintmax_t a1, uintmax_t a0, uintmax_t b1, uintmax_t b0) -{ - while ((a0 & 1) == 0) - rsh2 (a1, a0, a1, a0, 1); - while ((b0 & 1) == 0) - rsh2 (b1, b0, b1, b0, 1); - - for (;;) - { - if ((b1 | a1) == 0) - { - *r1 = 0; - return gcd_odd (b0, a0); - } - - if (gt2 (a1, a0, b1, b0)) - { - sub_ddmmss (a1, a0, a1, a0, b1, b0); - do - rsh2 (a1, a0, a1, a0, 1); - while ((a0 & 1) == 0); - } - else if (gt2 (b1, b0, a1, a0)) - { - sub_ddmmss (b1, b0, b1, b0, a1, a0); - do - rsh2 (b1, b0, b1, b0, 1); - while ((b0 & 1) == 0); - } - else - break; - } - - *r1 = a1; - return a0; -} - -static void -factor_insert_multiplicity (struct factors *factors, - uintmax_t prime, unsigned int m) -{ - unsigned int nfactors = factors->nfactors; - uintmax_t *p = factors->p; - unsigned char *e = factors->e; - - /* Locate position for insert new or increment e. */ - int i; - for (i = nfactors - 1; i >= 0; i--) - { - if (p[i] <= prime) - break; - } - - if (i < 0 || p[i] != prime) - { - for (int j = nfactors - 1; j > i; j--) - { - p[j + 1] = p[j]; - e[j + 1] = e[j]; - } - p[i + 1] = prime; - e[i + 1] = m; - factors->nfactors = nfactors + 1; - } - else - { - e[i] += m; - } -} - -#define factor_insert(f, p) factor_insert_multiplicity(f, p, 1) - -static void -factor_insert_large (struct factors *factors, - uintmax_t p1, uintmax_t p0) -{ - if (p1 > 0) - { - assert (factors->plarge[1] == 0); - factors->plarge[0] = p0; - factors->plarge[1] = p1; - } - else - factor_insert (factors, p0); -} - -#if HAVE_GMP -static void mp_factor (mpz_t, struct mp_factors *); - -static void -mp_factor_init (struct mp_factors *factors) -{ - factors->p = malloc (1); - factors->e = malloc (1); - factors->nfactors = 0; -} - -static void -mp_factor_clear (struct mp_factors *factors) -{ - for (unsigned int i = 0; i < factors->nfactors; i++) - mpz_clear (factors->p[i]); - - free (factors->p); - free (factors->e); -} - -static void -mp_factor_insert (struct mp_factors *factors, mpz_t prime) -{ - unsigned long int nfactors = factors->nfactors; - mpz_t *p = factors->p; - unsigned long int *e = factors->e; - long i; - - /* Locate position for insert new or increment e. */ - for (i = nfactors - 1; i >= 0; i--) - { - if (mpz_cmp (p[i], prime) <= 0) - break; - } - - if (i < 0 || mpz_cmp (p[i], prime) != 0) - { - p = realloc (p, (nfactors + 1) * sizeof p[0]); - e = realloc (e, (nfactors + 1) * sizeof e[0]); - - mpz_init (p[nfactors]); - for (long j = nfactors - 1; j > i; j--) - { - mpz_set (p[j + 1], p[j]); - e[j + 1] = e[j]; - } - mpz_set (p[i + 1], prime); - e[i + 1] = 1; - - factors->p = p; - factors->e = e; - factors->nfactors = nfactors + 1; - } - else - { - e[i] += 1; - } -} - -static void -mp_factor_insert_ui (struct mp_factors *factors, unsigned long int prime) -{ - mpz_t pz; - - mpz_init_set_ui (pz, prime); - mp_factor_insert (factors, pz); - mpz_clear (pz); -} -#endif /* HAVE_GMP */ - - -#define P(a,b,c,d) a, -static const unsigned char primes_diff[] = { -#include "primes.h" -0,0,0,0,0,0,0 /* 7 sentinels for 8-way loop */ -}; -#undef P - -#define PRIMES_PTAB_ENTRIES (sizeof(primes_diff) / sizeof(primes_diff[0]) - 8 + 1) - -#define P(a,b,c,d) b, -static const unsigned char primes_diff8[] = { -#include "primes.h" -0,0,0,0,0,0,0 /* 7 sentinels for 8-way loop */ -}; -#undef P - -struct primes_dtab -{ - uintmax_t binv, lim; -}; - -#define P(a,b,c,d) {c,d}, -static const struct primes_dtab primes_dtab[] = { -#include "primes.h" -{1,0},{1,0},{1,0},{1,0},{1,0},{1,0},{1,0} /* 7 sentinels for 8-way loop */ -}; -#undef P - -/* This flag is honoured just in the GMP code. */ -static int flag_verbose = 0; - -/* Prove primality or run probabilistic tests. */ -static int flag_prove_primality = 1; - -/* Number of Miller-Rabin tests to run when not proving primality. */ -#define MR_REPS 25 - -#define LIKELY(cond) __builtin_expect ((cond) != 0, 1) -#define UNLIKELY(cond) __builtin_expect ((cond) != 0, 0) - -static void -factor_insert_refind (struct factors *factors, uintmax_t p, unsigned int i, - unsigned int off) -{ - for (unsigned int j = 0; j < off; j++) - p += primes_diff[i + j]; - factor_insert (factors, p); -} - -/* Trial division with odd primes uses the following trick. - - Let p be an odd prime, and B = 2^{W_TYPE_SIZE}. For simplicity, - consider the case t < B (this is the second loop below). - - From our tables we get - - binv = p^{-1} (mod B) - lim = floor ( (B-1) / p ). - - First assume that t is a multiple of p, t = q * p. Then 0 <= q <= - lim (and all quotients in this range occur for some t). - - Then t = q * p is true also (mod B), and p is invertible we get - - q = t * binv (mod B). - - Next, assume that t is *not* divisible by p. Since multiplication - by binv (mod B) is a one-to-one mapping, - - t * binv (mod B) > lim, - - because all the smaller values are already taken. - - This can be summed up by saying that the function - - q(t) = binv * t (mod B) - - is a permutation of the range 0 <= t < B, with the curious property - that it maps the multiples of p onto the range 0 <= q <= lim, in - order, and the non-multiples of p onto the range lim < q < B. - */ - -static uintmax_t -factor_using_division (uintmax_t *t1p, uintmax_t t1, uintmax_t t0, - struct factors *factors) -{ - if (t0 % 2 == 0) - { - unsigned int cnt; - - if (t0 == 0) - { - count_trailing_zeros (cnt, t1); - t0 = t1 >> cnt; - t1 = 0; - cnt += W_TYPE_SIZE; - } - else - { - count_trailing_zeros (cnt, t0); - rsh2 (t1, t0, t1, t0, cnt); - } - - factor_insert_multiplicity (factors, 2, cnt); - } - - uintmax_t p = 3; - unsigned int i; - for (i = 0; t1 > 0 && i < PRIMES_PTAB_ENTRIES; i++) - { - for (;;) - { - uintmax_t q1, q0, hi, lo; - - q0 = t0 * primes_dtab[i].binv; - umul_ppmm (hi, lo, q0, p); - if (hi > t1) - break; - hi = t1 - hi; - q1 = hi * primes_dtab[i].binv; - if (LIKELY (q1 > primes_dtab[i].lim)) - break; - t1 = q1; t0 = q0; - factor_insert (factors, p); - } - p += primes_diff[i + 1]; - } - if (t1p) - *t1p = t1; - -#define DIVBLOCK(I) \ - do { \ - for (;;) \ - { \ - q = t0 * pd[I].binv; \ - if (LIKELY (q > pd[I].lim)) \ - break; \ - t0 = q; \ - factor_insert_refind (factors, p, i + 1, I); \ - } \ - } while (0) - - for (; i < PRIMES_PTAB_ENTRIES; i += 8) - { - uintmax_t q; - const struct primes_dtab *pd = &primes_dtab[i]; - DIVBLOCK (0); - DIVBLOCK (1); - DIVBLOCK (2); - DIVBLOCK (3); - DIVBLOCK (4); - DIVBLOCK (5); - DIVBLOCK (6); - DIVBLOCK (7); - - p += primes_diff8[i]; - if (p * p > t0) - break; - } - - return t0; -} - -#if HAVE_GMP -static void -mp_factor_using_division (mpz_t t, struct mp_factors *factors) -{ - mpz_t q; - unsigned long int p; - - if (flag_verbose > 0) - { - printf ("[trial division] "); - } - - mpz_init (q); - - p = mpz_scan1 (t, 0); - mpz_div_2exp (t, t, p); - while (p) - { - mp_factor_insert_ui (factors, 2); - --p; - } - - p = 3; - for (unsigned int i = 1; i <= PRIMES_PTAB_ENTRIES;) - { - if (! mpz_divisible_ui_p (t, p)) - { - p += primes_diff[i++]; - if (mpz_cmp_ui (t, p * p) < 0) - break; - } - else - { - mpz_tdiv_q_ui (t, t, p); - mp_factor_insert_ui (factors, p); - } - } - - mpz_clear (q); -} -#endif - -/* Entry i contains (2i+1)^(-1) mod 2^8. */ -static const unsigned char binvert_table[128] = -{ - 0x01, 0xAB, 0xCD, 0xB7, 0x39, 0xA3, 0xC5, 0xEF, - 0xF1, 0x1B, 0x3D, 0xA7, 0x29, 0x13, 0x35, 0xDF, - 0xE1, 0x8B, 0xAD, 0x97, 0x19, 0x83, 0xA5, 0xCF, - 0xD1, 0xFB, 0x1D, 0x87, 0x09, 0xF3, 0x15, 0xBF, - 0xC1, 0x6B, 0x8D, 0x77, 0xF9, 0x63, 0x85, 0xAF, - 0xB1, 0xDB, 0xFD, 0x67, 0xE9, 0xD3, 0xF5, 0x9F, - 0xA1, 0x4B, 0x6D, 0x57, 0xD9, 0x43, 0x65, 0x8F, - 0x91, 0xBB, 0xDD, 0x47, 0xC9, 0xB3, 0xD5, 0x7F, - 0x81, 0x2B, 0x4D, 0x37, 0xB9, 0x23, 0x45, 0x6F, - 0x71, 0x9B, 0xBD, 0x27, 0xA9, 0x93, 0xB5, 0x5F, - 0x61, 0x0B, 0x2D, 0x17, 0x99, 0x03, 0x25, 0x4F, - 0x51, 0x7B, 0x9D, 0x07, 0x89, 0x73, 0x95, 0x3F, - 0x41, 0xEB, 0x0D, 0xF7, 0x79, 0xE3, 0x05, 0x2F, - 0x31, 0x5B, 0x7D, 0xE7, 0x69, 0x53, 0x75, 0x1F, - 0x21, 0xCB, 0xED, 0xD7, 0x59, 0xC3, 0xE5, 0x0F, - 0x11, 0x3B, 0x5D, 0xC7, 0x49, 0x33, 0x55, 0xFF -}; - -/* Compute n^(-1) mod B, using a Newton iteration. */ -#define binv(inv,n) \ - do { \ - uintmax_t __n = (n); \ - uintmax_t __inv; \ - \ - __inv = binvert_table[(__n / 2) & 0x7F]; /* 8 */ \ - if (W_TYPE_SIZE > 8) __inv = 2 * __inv - __inv * __inv * __n; \ - if (W_TYPE_SIZE > 16) __inv = 2 * __inv - __inv * __inv * __n; \ - if (W_TYPE_SIZE > 32) __inv = 2 * __inv - __inv * __inv * __n; \ - \ - if (W_TYPE_SIZE > 64) \ - { \ - int __invbits = 64; \ - do { \ - __inv = 2 * __inv - __inv * __inv * __n; \ - __invbits *= 2; \ - } while (__invbits < W_TYPE_SIZE); \ - } \ - \ - (inv) = __inv; \ - } while (0) - -/* q = u / d, assuming d|u. */ -#define divexact_21(q1, q0, u1, u0, d) \ - do { \ - uintmax_t _di, _q0; \ - binv (_di, (d)); \ - _q0 = (u0) * _di; \ - if ((u1) >= (d)) \ - { \ - uintmax_t _p1, _p0; \ - umul_ppmm (_p1, _p0, _q0, d); \ - (q1) = ((u1) - _p1) * _di; \ - (q0) = _q0; \ - } \ - else \ - { \ - (q0) = _q0; \ - (q1) = 0; \ - } \ - } while(0) - -/* x B (mod n). */ -#define redcify(r_prim, r, n) \ - do { \ - uintmax_t _redcify_q; \ - udiv_qrnnd (_redcify_q, r_prim, r, 0, n); \ - } while (0) - -/* x B^2 (mod n). Requires x > 0, n1 < B/2 */ -#define redcify2(r1, r0, x, n1, n0) \ - do { \ - uintmax_t _r1, _r0, _i; \ - if ((x) < (n1)) \ - { \ - _r1 = (x); _r0 = 0; \ - _i = W_TYPE_SIZE; \ - } \ - else \ - { \ - _r1 = 0; _r0 = (x); \ - _i = 2*W_TYPE_SIZE; \ - } \ - while (_i-- > 0) \ - { \ - lsh2 (_r1, _r0, _r1, _r0, 1); \ - if (ge2 (_r1, _r0, (n1), (n0))) \ - sub_ddmmss (_r1, _r0, _r1, _r0, (n1), (n0)); \ - } \ - (r1) = _r1; \ - (r0) = _r0; \ - } while (0) - -/* Modular two-word multiplication, r = a * b mod m, with mi = m^(-1) mod B. - Both a and b must be in redc form, the result will be in redc form too. */ -static inline uintmax_t -mulredc (uintmax_t a, uintmax_t b, uintmax_t m, uintmax_t mi) -{ - uintmax_t rh, rl, q, th, tl, xh; - - umul_ppmm (rh, rl, a, b); - q = rl * mi; - umul_ppmm (th, tl, q, m); - xh = rh - th; - if (rh < th) - xh += m; - - return xh; -} - -/* Modular two-word multiplication, r = a * b mod m, with mi = m^(-1) mod B. - Both a and b must be in redc form, the result will be in redc form too. - For performance reasons, the most significant bit of m must be clear. */ -static uintmax_t -mulredc2 (uintmax_t *r1p, - uintmax_t a1, uintmax_t a0, uintmax_t b1, uintmax_t b0, - uintmax_t m1, uintmax_t m0, uintmax_t mi) -{ - uintmax_t r1, r0, q, p1, p0, t1, t0, s1, s0; - mi = -mi; - assert ( (a1 >> (W_TYPE_SIZE - 1)) == 0); - assert ( (b1 >> (W_TYPE_SIZE - 1)) == 0); - assert ( (m1 >> (W_TYPE_SIZE - 1)) == 0); - - /* First compute a0 * <b1, b0> B^{-1} - +-----+ - |a0 b0| - +--+--+--+ - |a0 b1| - +--+--+--+ - |q0 m0| - +--+--+--+ - |q0 m1| - -+--+--+--+ - |r1|r0| 0| - +--+--+--+ - */ - umul_ppmm (t1, t0, a0, b0); - umul_ppmm (r1, r0, a0, b1); - q = mi * t0; - umul_ppmm (p1, p0, q, m0); - umul_ppmm (s1, s0, q, m1); - r0 += (t0 != 0); /* Carry */ - add_ssaaaa (r1, r0, r1, r0, 0, p1); - add_ssaaaa (r1, r0, r1, r0, 0, t1); - add_ssaaaa (r1, r0, r1, r0, s1, s0); - - /* Next, (a1 * <b1, b0> + <r1, r0> B^{-1} - +-----+ - |a1 b0| - +--+--+ - |r1|r0| - +--+--+--+ - |a1 b1| - +--+--+--+ - |q1 m0| - +--+--+--+ - |q1 m1| - -+--+--+--+ - |r1|r0| 0| - +--+--+--+ - */ - umul_ppmm (t1, t0, a1, b0); - umul_ppmm (s1, s0, a1, b1); - add_ssaaaa (t1, t0, t1, t0, 0, r0); - q = mi * t0; - add_ssaaaa (r1, r0, s1, s0, 0, r1); - umul_ppmm (p1, p0, q, m0); - umul_ppmm (s1, s0, q, m1); - r0 += (t0 != 0); /* Carry */ - add_ssaaaa (r1, r0, r1, r0, 0, p1); - add_ssaaaa (r1, r0, r1, r0, 0, t1); - add_ssaaaa (r1, r0, r1, r0, s1, s0); - - if (ge2 (r1, r0, m1, m0)) - sub_ddmmss (r1, r0, r1, r0, m1, m0); - - *r1p = r1; - return r0; -} - -static uintmax_t _GL_ATTRIBUTE_CONST -powm (uintmax_t b, uintmax_t e, uintmax_t n, uintmax_t ni, uintmax_t one) -{ - uintmax_t y = one; - - if (e & 1) - y = b; - - while (e != 0) - { - b = mulredc (b, b, n, ni); - e >>= 1; - - if (e & 1) - y = mulredc (y, b, n, ni); - } - - return y; -} - -static uintmax_t -powm2 (uintmax_t *r1m, - const uintmax_t *bp, const uintmax_t *ep, const uintmax_t *np, - uintmax_t ni, const uintmax_t *one) -{ - uintmax_t r1, r0, b1, b0, n1, n0; - unsigned int i; - uintmax_t e; - - b0 = bp[0]; - b1 = bp[1]; - n0 = np[0]; - n1 = np[1]; - - r0 = one[0]; - r1 = one[1]; - - for (e = ep[0], i = W_TYPE_SIZE; i > 0; i--, e >>= 1) - { - if (e & 1) - { - r0 = mulredc2 (r1m, r1, r0, b1, b0, n1, n0, ni); - r1 = *r1m; - } - b0 = mulredc2 (r1m, b1, b0, b1, b0, n1, n0, ni); - b1 = *r1m; - } - for (e = ep[1]; e > 0; e >>= 1) - { - if (e & 1) - { - r0 = mulredc2 (r1m, r1, r0, b1, b0, n1, n0, ni); - r1 = *r1m; - } - b0 = mulredc2 (r1m, b1, b0, b1, b0, n1, n0, ni); - b1 = *r1m; - } - *r1m = r1; - return r0; -} - -static bool _GL_ATTRIBUTE_CONST -millerrabin (uintmax_t n, uintmax_t ni, uintmax_t b, uintmax_t q, - unsigned int k, uintmax_t one) -{ - uintmax_t y = powm (b, q, n, ni, one); - - uintmax_t nm1 = n - one; /* -1, but in redc representation. */ - - if (y == one || y == nm1) - return true; - - for (unsigned int i = 1; i < k; i++) - { - y = mulredc (y, y, n, ni); - - if (y == nm1) - return true; - if (y == one) - return false; - } - return false; -} - -static bool -millerrabin2 (const uintmax_t *np, uintmax_t ni, const uintmax_t *bp, - const uintmax_t *qp, unsigned int k, const uintmax_t *one) -{ - uintmax_t y1, y0, nm1_1, nm1_0, r1m; - - y0 = powm2 (&r1m, bp, qp, np, ni, one); - y1 = r1m; - - if (y0 == one[0] && y1 == one[1]) - return true; - - sub_ddmmss (nm1_1, nm1_0, np[1], np[0], one[1], one[0]); - - if (y0 == nm1_0 && y1 == nm1_1) - return true; - - for (unsigned int i = 1; i < k; i++) - { - y0 = mulredc2 (&r1m, y1, y0, y1, y0, np[1], np[0], ni); - y1 = r1m; - - if (y0 == nm1_0 && y1 == nm1_1) - return true; - if (y0 == one[0] && y1 == one[1]) - return false; - } - return false; -} - -#if HAVE_GMP -static bool -mp_millerrabin (mpz_srcptr n, mpz_srcptr nm1, mpz_ptr x, mpz_ptr y, - mpz_srcptr q, unsigned long int k) -{ - mpz_powm (y, x, q, n); - - if (mpz_cmp_ui (y, 1) == 0 || mpz_cmp (y, nm1) == 0) - return true; - - for (unsigned long int i = 1; i < k; i++) - { - mpz_powm_ui (y, y, 2, n); - if (mpz_cmp (y, nm1) == 0) - return true; - if (mpz_cmp_ui (y, 1) == 0) - return false; - } - return false; -} -#endif - -/* Lucas' prime test. The number of iterations vary greatly, up to a few dozen - have been observed. The average seem to be about 2. */ -static bool -prime_p (uintmax_t n) -{ - int k; - bool is_prime; - uintmax_t a_prim, one, ni; - struct factors factors; - - if (n <= 1) - return false; - - /* We have already casted out small primes. */ - if (n < (uintmax_t) FIRST_OMITTED_PRIME * FIRST_OMITTED_PRIME) - return true; - - /* Precomputation for Miller-Rabin. */ - uintmax_t q = n - 1; - for (k = 0; (q & 1) == 0; k++) - q >>= 1; - - uintmax_t a = 2; - binv (ni, n); /* ni <- 1/n mod B */ - redcify (one, 1, n); - addmod (a_prim, one, one, n); /* i.e., redcify a = 2 */ - - /* Perform a Miller-Rabin test, finds most composites quickly. */ - if (!millerrabin (n, ni, a_prim, q, k, one)) - return false; - - if (flag_prove_primality) - { - /* Factor n-1 for Lucas. */ - factor (0, n - 1, &factors); - } - - /* Loop until Lucas proves our number prime, or Miller-Rabin proves our - number composite. */ - for (unsigned int r = 0; r < PRIMES_PTAB_ENTRIES; r++) - { - if (flag_prove_primality) - { - is_prime = true; - for (unsigned int i = 0; i < factors.nfactors && is_prime; i++) - { - is_prime = powm (a_prim, (n - 1) / factors.p[i], n, ni, one) != one; - } - } - else - { - /* After enough Miller-Rabin runs, be content. */ - is_prime = (r == MR_REPS - 1); - } - - if (is_prime) - return true; - - a += primes_diff[r]; /* Establish new base. */ - - /* The following is equivalent to redcify (a_prim, a, n). It runs faster - on most processors, since it avoids udiv_qrnnd. If we go down the - udiv_qrnnd_preinv path, this code should be replaced. */ - { - uintmax_t dummy, s1, s0; - umul_ppmm (s1, s0, one, a); - if (LIKELY (s1 == 0)) - a_prim = s0 % n; - else - udiv_qrnnd (dummy, a_prim, s1, s0, n); - } - - if (!millerrabin (n, ni, a_prim, q, k, one)) - return false; - } - - fprintf (stderr, "Lucas prime test failure. This should not happen\n"); - abort (); -} - -static bool -prime2_p (uintmax_t n1, uintmax_t n0) -{ - uintmax_t q[2], nm1[2]; - uintmax_t a_prim[2]; - uintmax_t one[2]; - uintmax_t na[2]; - uintmax_t ni; - unsigned int k; - struct factors factors; - - if (n1 == 0) - return prime_p (n0); - - nm1[1] = n1 - (n0 == 0); - nm1[0] = n0 - 1; - if (nm1[0] == 0) - { - count_trailing_zeros (k, nm1[1]); - - q[0] = nm1[1] >> k; - q[1] = 0; - k += W_TYPE_SIZE; - } - else - { - count_trailing_zeros (k, nm1[0]); - rsh2 (q[1], q[0], nm1[1], nm1[0], k); - } - - uintmax_t a = 2; - binv (ni, n0); - redcify2 (one[1], one[0], 1, n1, n0); - addmod2 (a_prim[1], a_prim[0], one[1], one[0], one[1], one[0], n1, n0); - - /* FIXME: Use scalars or pointers in arguments? Some consistency needed. */ - na[0] = n0; - na[1] = n1; - - if (!millerrabin2 (na, ni, a_prim, q, k, one)) - return false; - - if (flag_prove_primality) - { - /* Factor n-1 for Lucas. */ - factor (nm1[1], nm1[0], &factors); - } - - /* Loop until Lucas proves our number prime, or Miller-Rabin proves our - number composite. */ - for (unsigned int r = 0; r < PRIMES_PTAB_ENTRIES; r++) - { - bool is_prime; - uintmax_t e[2], y[2]; - - if (flag_prove_primality) - { - is_prime = true; - if (factors.plarge[1]) - { - uintmax_t pi; - binv (pi, factors.plarge[0]); - e[0] = pi * nm1[0]; - e[1] = 0; - y[0] = powm2 (&y[1], a_prim, e, na, ni, one); - is_prime = (y[0] != one[0] || y[1] != one[1]); - } - for (unsigned int i = 0; i < factors.nfactors && is_prime; i++) - { - /* FIXME: We always have the factor 2. Do we really need to handle it - here? We have done the same powering as part of millerrabin. */ - if (factors.p[i] == 2) - rsh2 (e[1], e[0], nm1[1], nm1[0], 1); - else - divexact_21 (e[1], e[0], nm1[1], nm1[0], factors.p[i]); - y[0] = powm2 (&y[1], a_prim, e, na, ni, one); - is_prime = (y[0] != one[0] || y[1] != one[1]); - } - } - else - { - /* After enough Miller-Rabin runs, be content. */ - is_prime = (r == MR_REPS - 1); - } - - if (is_prime) - return true; - - a += primes_diff[r]; /* Establish new base. */ - redcify2 (a_prim[1], a_prim[0], a, n1, n0); - - if (!millerrabin2 (na, ni, a_prim, q, k, one)) - return false; - } - - fprintf (stderr, "Lucas prime test failure. This should not happen\n"); - abort (); -} - -#if HAVE_GMP -static bool -mp_prime_p (mpz_t n) -{ - bool is_prime; - mpz_t q, a, nm1, tmp; - struct mp_factors factors; - - if (mpz_cmp_ui (n, 1) <= 0) - return false; - - /* We have already casted out small primes. */ - if (mpz_cmp_ui (n, (long) FIRST_OMITTED_PRIME * FIRST_OMITTED_PRIME) < 0) - return true; - - mpz_inits (q, a, nm1, tmp, NULL); - - /* Precomputation for Miller-Rabin. */ - mpz_sub_ui (nm1, n, 1); - - /* Find q and k, where q is odd and n = 1 + 2**k * q. */ - unsigned long int k = mpz_scan1 (nm1, 0); - mpz_tdiv_q_2exp (q, nm1, k); - - mpz_set_ui (a, 2); - - /* Perform a Miller-Rabin test, finds most composites quickly. */ - if (!mp_millerrabin (n, nm1, a, tmp, q, k)) - { - is_prime = false; - goto ret2; - } - - if (flag_prove_primality) - { - /* Factor n-1 for Lucas. */ - mpz_set (tmp, nm1); - mp_factor (tmp, &factors); - } - - /* Loop until Lucas proves our number prime, or Miller-Rabin proves our - number composite. */ - for (unsigned int r = 0; r < PRIMES_PTAB_ENTRIES; r++) - { - if (flag_prove_primality) - { - is_prime = true; - for (unsigned long int i = 0; i < factors.nfactors && is_prime; i++) - { - mpz_divexact (tmp, nm1, factors.p[i]); - mpz_powm (tmp, a, tmp, n); - is_prime = mpz_cmp_ui (tmp, 1) != 0; - } - } - else - { - /* After enough Miller-Rabin runs, be content. */ - is_prime = (r == MR_REPS - 1); - } - - if (is_prime) - goto ret1; - - mpz_add_ui (a, a, primes_diff[r]); /* Establish new base. */ - - if (!mp_millerrabin (n, nm1, a, tmp, q, k)) - { - is_prime = false; - goto ret1; - } - } - - fprintf (stderr, "Lucas prime test failure. This should not happen\n"); - abort (); - - ret1: - if (flag_prove_primality) - mp_factor_clear (&factors); - ret2: - mpz_clears (q, a, nm1, tmp, NULL); - - return is_prime; -} -#endif - -static void -factor_using_pollard_rho (uintmax_t n, unsigned long int a, - struct factors *factors) -{ - uintmax_t x, z, y, P, t, ni, g; - - unsigned long int k = 1; - unsigned long int l = 1; - - redcify (P, 1, n); - addmod (x, P, P, n); /* i.e., redcify(2) */ - y = z = x; - - while (n != 1) - { - assert (a < n); - - binv (ni, n); /* FIXME: when could we use old 'ni' value? */ - - for (;;) - { - do - { - x = mulredc (x, x, n, ni); - addmod (x, x, a, n); - - submod (t, z, x, n); - P = mulredc (P, t, n, ni); - - if (k % 32 == 1) - { - if (gcd_odd (P, n) != 1) - goto factor_found; - y = x; - } - } - while (--k != 0); - - z = x; - k = l; - l = 2 * l; - for (unsigned long int i = 0; i < k; i++) - { - x = mulredc (x, x, n, ni); - addmod (x, x, a, n); - } - y = x; - } - - factor_found: - do - { - y = mulredc (y, y, n, ni); - addmod (y, y, a, n); - - submod (t, z, y, n); - g = gcd_odd (t, n); - } - while (g == 1); - - n = n / g; - - if (!prime_p (g)) - factor_using_pollard_rho (g, a + 1, factors); - else - factor_insert (factors, g); - - if (prime_p (n)) - { - factor_insert (factors, n); - break; - } - - x = x % n; - z = z % n; - y = y % n; - } -} - -static void -factor_using_pollard_rho2 (uintmax_t n1, uintmax_t n0, unsigned long int a, - struct factors *factors) -{ - uintmax_t x1, x0, z1, z0, y1, y0, P1, P0, t1, t0, ni, g1, g0, r1m; - - unsigned long int k = 1; - unsigned long int l = 1; - - redcify2 (P1, P0, 1, n1, n0); - addmod2 (x1, x0, P1, P0, P1, P0, n1, n0); /* i.e., redcify(2) */ - y1 = z1 = x1; - y0 = z0 = x0; - - while (n1 != 0 || n0 != 1) - { - binv (ni, n0); - - for (;;) - { - do - { - x0 = mulredc2 (&r1m, x1, x0, x1, x0, n1, n0, ni); - x1 = r1m; - addmod2 (x1, x0, x1, x0, 0, a, n1, n0); - - submod2 (t1, t0, z1, z0, x1, x0, n1, n0); - P0 = mulredc2 (&r1m, P1, P0, t1, t0, n1, n0, ni); - P1 = r1m; - - if (k % 32 == 1) - { - g0 = gcd2_odd (&g1, P1, P0, n1, n0); - if (g1 != 0 || g0 != 1) - goto factor_found; - y1 = x1; y0 = x0; - } - } - while (--k != 0); - - z1 = x1; z0 = x0; - k = l; - l = 2 * l; - for (unsigned long int i = 0; i < k; i++) - { - x0 = mulredc2 (&r1m, x1, x0, x1, x0, n1, n0, ni); - x1 = r1m; - addmod2 (x1, x0, x1, x0, 0, a, n1, n0); - } - y1 = x1; y0 = x0; - } - - factor_found: - do - { - y0 = mulredc2 (&r1m, y1, y0, y1, y0, n1, n0, ni); - y1 = r1m; - addmod2 (y1, y0, y1, y0, 0, a, n1, n0); - - submod2 (t1, t0, z1, z0, y1, y0, n1, n0); - g0 = gcd2_odd (&g1, t1, t0, n1, n0); - } - while (g1 == 0 && g0 == 1); - - if (g1 == 0) - { - /* The found factor is one word. */ - divexact_21 (n1, n0, n1, n0, g0); /* n = n / g */ - - if (!prime_p (g0)) - factor_using_pollard_rho (g0, a + 1, factors); - else - factor_insert (factors, g0); - } - else - { - /* The found factor is two words. This is highly unlikely, thus hard - to trigger. Please be careful before you change this code! */ - uintmax_t ginv; - - binv (ginv, g0); /* Compute n = n / g. Since the result will */ - n0 = ginv * n0; /* fit one word, we can compute the quotient */ - n1 = 0; /* modulo B, ignoring the high divisor word. */ - - if (!prime2_p (g1, g0)) - factor_using_pollard_rho2 (g1, g0, a + 1, factors); - else - factor_insert_large (factors, g1, g0); - } - - if (n1 == 0) - { - if (prime_p (n0)) - { - factor_insert (factors, n0); - break; - } - - factor_using_pollard_rho (n0, a, factors); - return; - } - - if (prime2_p (n1, n0)) - { - factor_insert_large (factors, n1, n0); - break; - } - - x0 = mod2 (&x1, x1, x0, n1, n0); - z0 = mod2 (&z1, z1, z0, n1, n0); - y0 = mod2 (&y1, y1, y0, n1, n0); - } -} - -#if HAVE_GMP -static void -mp_factor_using_pollard_rho (mpz_t n, unsigned long int a, - struct mp_factors *factors) -{ - mpz_t x, z, y, P; - mpz_t t, t2; - unsigned long long int k, l; - - if (flag_verbose > 0) - { - printf ("[pollard-rho (%lu)] ", a); - } - - mpz_inits (t, t2, NULL); - mpz_init_set_si (y, 2); - mpz_init_set_si (x, 2); - mpz_init_set_si (z, 2); - mpz_init_set_ui (P, 1); - k = 1; - l = 1; - - while (mpz_cmp_ui (n, 1) != 0) - { - for (;;) - { - do - { - mpz_mul (t, x, x); - mpz_mod (x, t, n); - mpz_add_ui (x, x, a); - - mpz_sub (t, z, x); - mpz_mul (t2, P, t); - mpz_mod (P, t2, n); - - if (k % 32 == 1) - { - mpz_gcd (t, P, n); - if (mpz_cmp_ui (t, 1) != 0) - goto factor_found; - mpz_set (y, x); - } - } - while (--k != 0); - - mpz_set (z, x); - k = l; - l = 2 * l; - for (unsigned long long int i = 0; i < k; i++) - { - mpz_mul (t, x, x); - mpz_mod (x, t, n); - mpz_add_ui (x, x, a); - } - mpz_set (y, x); - } - - factor_found: - do - { - mpz_mul (t, y, y); - mpz_mod (y, t, n); - mpz_add_ui (y, y, a); - - mpz_sub (t, z, y); - mpz_gcd (t, t, n); - } - while (mpz_cmp_ui (t, 1) == 0); - - mpz_divexact (n, n, t); /* divide by t, before t is overwritten */ - - if (!mp_prime_p (t)) - { - if (flag_verbose > 0) - { - printf ("[composite factor--restarting pollard-rho] "); - } - mp_factor_using_pollard_rho (t, a + 1, factors); - } - else - { - mp_factor_insert (factors, t); - } - - if (mp_prime_p (n)) - { - mp_factor_insert (factors, n); - break; - } - - mpz_mod (x, x, n); - mpz_mod (z, z, n); - mpz_mod (y, y, n); - } - - mpz_clears (P, t2, t, z, x, y, NULL); -} -#endif - -/* FIXME: Maybe better to use an iteration converging to 1/sqrt(n)? If - algorithm is replaced, consider also returning the remainder. */ -static uintmax_t _GL_ATTRIBUTE_CONST -isqrt (uintmax_t n) -{ - uintmax_t x; - unsigned c; - if (n == 0) - return 0; - - count_leading_zeros (c, n); - - /* Make x > sqrt(n). This will be invariant through the loop. */ - x = (uintmax_t) 1 << ((W_TYPE_SIZE + 1 - c) / 2); - - for (;;) - { - uintmax_t y = (x + n/x) / 2; - if (y >= x) - return x; - - x = y; - } -} - -static uintmax_t _GL_ATTRIBUTE_CONST -isqrt2 (uintmax_t nh, uintmax_t nl) -{ - unsigned int shift; - uintmax_t x; - - /* Ensures the remainder fits in an uintmax_t. */ - assert (nh < ((uintmax_t) 1 << (W_TYPE_SIZE - 2))); - - if (nh == 0) - return isqrt (nl); - - count_leading_zeros (shift, nh); - shift &= ~1; - - /* Make x > sqrt(n) */ - x = isqrt ( (nh << shift) + (nl >> (W_TYPE_SIZE - shift))) + 1; - x <<= (W_TYPE_SIZE - shift) / 2; - - /* Do we need more than one iteration? */ - for (;;) - { - uintmax_t q, r, y; - udiv_qrnnd (q, r, nh, nl, x); - y = (x + q) / 2; - - if (y >= x) - { - uintmax_t hi, lo; - umul_ppmm (hi, lo, x + 1, x + 1); - assert (gt2 (hi, lo, nh, nl)); - - umul_ppmm (hi, lo, x, x); - assert (ge2 (nh, nl, hi, lo)); - sub_ddmmss (hi, lo, nh, nl, hi, lo); - assert (hi == 0); - - return x; - } - - x = y; - } -} - -/* MAGIC[N] has a bit i set iff i is a quadratic residue mod N. */ -#define MAGIC64 ((uint64_t) 0x0202021202030213ULL) -#define MAGIC63 ((uint64_t) 0x0402483012450293ULL) -#define MAGIC65 ((uint64_t) 0x218a019866014613ULL) -#define MAGIC11 0x23b - -/* Returns the square root if the input is a square, otherwise 0. */ -static uintmax_t _GL_ATTRIBUTE_CONST -is_square (uintmax_t x) -{ - /* Uses the tests suggested by Cohen. Excludes 99% of the non-squares before - computing the square root. */ - if (((MAGIC64 >> (x & 63)) & 1) - && ((MAGIC63 >> (x % 63)) & 1) - /* Both 0 and 64 are squares mod (65) */ - && ((MAGIC65 >> ((x % 65) & 63)) & 1) - && ((MAGIC11 >> (x % 11) & 1))) - { - uintmax_t r = isqrt (x); - if (r*r == x) - return r; - } - return 0; -} - -/* invtab[i] = floor(0x10000 / (0x100 + i) */ -static const unsigned short invtab[0x81] = - { - 0x200, - 0x1fc, 0x1f8, 0x1f4, 0x1f0, 0x1ec, 0x1e9, 0x1e5, 0x1e1, - 0x1de, 0x1da, 0x1d7, 0x1d4, 0x1d0, 0x1cd, 0x1ca, 0x1c7, - 0x1c3, 0x1c0, 0x1bd, 0x1ba, 0x1b7, 0x1b4, 0x1b2, 0x1af, - 0x1ac, 0x1a9, 0x1a6, 0x1a4, 0x1a1, 0x19e, 0x19c, 0x199, - 0x197, 0x194, 0x192, 0x18f, 0x18d, 0x18a, 0x188, 0x186, - 0x183, 0x181, 0x17f, 0x17d, 0x17a, 0x178, 0x176, 0x174, - 0x172, 0x170, 0x16e, 0x16c, 0x16a, 0x168, 0x166, 0x164, - 0x162, 0x160, 0x15e, 0x15c, 0x15a, 0x158, 0x157, 0x155, - 0x153, 0x151, 0x150, 0x14e, 0x14c, 0x14a, 0x149, 0x147, - 0x146, 0x144, 0x142, 0x141, 0x13f, 0x13e, 0x13c, 0x13b, - 0x139, 0x138, 0x136, 0x135, 0x133, 0x132, 0x130, 0x12f, - 0x12e, 0x12c, 0x12b, 0x129, 0x128, 0x127, 0x125, 0x124, - 0x123, 0x121, 0x120, 0x11f, 0x11e, 0x11c, 0x11b, 0x11a, - 0x119, 0x118, 0x116, 0x115, 0x114, 0x113, 0x112, 0x111, - 0x10f, 0x10e, 0x10d, 0x10c, 0x10b, 0x10a, 0x109, 0x108, - 0x107, 0x106, 0x105, 0x104, 0x103, 0x102, 0x101, 0x100, - }; - -/* Compute q = [u/d], r = u mod d. Avoids slow hardware division for the case - that q < 0x40; here it instead uses a table of (Euclidian) inverses. */ -#define div_smallq(q, r, u, d) \ - do { \ - if ((u) / 0x40 < (d)) \ - { \ - int _cnt; \ - uintmax_t _dinv, _mask, _q, _r; \ - count_leading_zeros (_cnt, (d)); \ - _r = (u); \ - if (UNLIKELY (_cnt > (W_TYPE_SIZE - 8))) \ - { \ - _dinv = invtab[((d) << (_cnt + 8 - W_TYPE_SIZE)) - 0x80]; \ - _q = _dinv * _r >> (8 + W_TYPE_SIZE - _cnt); \ - } \ - else \ - { \ - _dinv = invtab[((d) >> (W_TYPE_SIZE - 8 - _cnt)) - 0x7f]; \ - _q = _dinv * (_r >> (W_TYPE_SIZE - 3 - _cnt)) >> 11; \ - } \ - _r -= _q*(d); \ - \ - _mask = -(uintmax_t) (_r >= (d)); \ - (r) = _r - (_mask & (d)); \ - (q) = _q - _mask; \ - assert ( (q) * (d) + (r) == u); \ - } \ - else \ - { \ - uintmax_t _q = (u) / (d); \ - (r) = (u) - _q * (d); \ - (q) = _q; \ - } \ - } while (0) - -/* Notes: Example N = 22117019. After first phase we find Q1 = 6314, Q - = 3025, P = 1737, representing F_{18} = (-6314, 2* 1737, 3025), - with 3025 = 55^2. - - Constructing the square root, we get Q1 = 55, Q = 8653, P = 4652, - representing G_0 = (-55, 2*4652, 8653). - - In the notation of the paper: - - S_{-1} = 55, S_0 = 8653, R_0 = 4652 - - Put - - t_0 = floor([q_0 + R_0] / S0) = 1 - R_1 = t_0 * S_0 - R_0 = 4001 - S_1 = S_{-1} +t_0 (R_0 - R_1) = 706 -*/ - -/* Multipliers, in order of efficiency: - 0.7268 3*5*7*11 = 1155 = 3 (mod 4) - 0.7317 3*5*7 = 105 = 1 - 0.7820 3*5*11 = 165 = 1 - 0.7872 3*5 = 15 = 3 - 0.8101 3*7*11 = 231 = 3 - 0.8155 3*7 = 21 = 1 - 0.8284 5*7*11 = 385 = 1 - 0.8339 5*7 = 35 = 3 - 0.8716 3*11 = 33 = 1 - 0.8774 3 = 3 = 3 - 0.8913 5*11 = 55 = 3 - 0.8972 5 = 5 = 1 - 0.9233 7*11 = 77 = 1 - 0.9295 7 = 7 = 3 - 0.9934 11 = 11 = 3 -*/ -#define QUEUE_SIZE 50 - -#if STAT_SQUFOF -# define Q_FREQ_SIZE 50 -/* Element 0 keeps the total */ -static unsigned int q_freq[Q_FREQ_SIZE + 1]; -# define MIN(a,b) ((a) < (b) ? (a) : (b)) -#endif - -/* Return true on success. Expected to fail only for numbers - >= 2^{2*W_TYPE_SIZE - 2}, or close to that limit. */ -static bool -factor_using_squfof (uintmax_t n1, uintmax_t n0, struct factors *factors) -{ - /* Uses algorithm and notation from - - SQUARE FORM FACTORIZATION - JASON E. GOWER AND SAMUEL S. WAGSTAFF, JR. - - http://homes.cerias.purdue.edu/~ssw/squfof.pdf - */ - - static const unsigned int multipliers_1[] = - { /* = 1 (mod 4) */ - 105, 165, 21, 385, 33, 5, 77, 1, 0 - }; - static const unsigned int multipliers_3[] = - { /* = 3 (mod 4) */ - 1155, 15, 231, 35, 3, 55, 7, 11, 0 - }; - - const unsigned int *m; - - struct { uintmax_t Q; uintmax_t P; } queue[QUEUE_SIZE]; - - if (n1 >= ((uintmax_t) 1 << (W_TYPE_SIZE - 2))) - return false; - - uintmax_t sqrt_n = isqrt2 (n1, n0); - - if (n0 == sqrt_n * sqrt_n) - { - uintmax_t p1, p0; - - umul_ppmm (p1, p0, sqrt_n, sqrt_n); - assert (p0 == n0); - - if (n1 == p1) - { - if (prime_p (sqrt_n)) - factor_insert_multiplicity (factors, sqrt_n, 2); - else - { - struct factors f; - - f.nfactors = 0; - if (!factor_using_squfof (0, sqrt_n, &f)) - { - /* Try pollard rho instead */ - factor_using_pollard_rho (sqrt_n, 1, &f); - } - /* Duplicate the new factors */ - for (unsigned int i = 0; i < f.nfactors; i++) - factor_insert_multiplicity (factors, f.p[i], 2*f.e[i]); - } - return true; - } - } - - /* Select multipliers so we always get n * mu = 3 (mod 4) */ - for (m = (n0 % 4 == 1) ? multipliers_3 : multipliers_1; - *m; m++) - { - uintmax_t S, Dh, Dl, Q1, Q, P, L, L1, B; - unsigned int i; - unsigned int mu = *m; - unsigned int qpos = 0; - - assert (mu * n0 % 4 == 3); - - /* In the notation of the paper, with mu * n == 3 (mod 4), we - get \Delta = 4 mu * n, and the paper's \mu is 2 mu. As far as - I understand it, the necessary bound is 4 \mu^3 < n, or 32 - mu^3 < n. - - However, this seems insufficient: With n = 37243139 and mu = - 105, we get a trivial factor, from the square 38809 = 197^2, - without any corresponding Q earlier in the iteration. - - Requiring 64 mu^3 < n seems sufficient. */ - if (n1 == 0) - { - if ((uintmax_t) mu*mu*mu >= n0 / 64) - continue; - } - else - { - if (n1 > ((uintmax_t) 1 << (W_TYPE_SIZE - 2)) / mu) - continue; - } - umul_ppmm (Dh, Dl, n0, mu); - Dh += n1 * mu; - - assert (Dl % 4 != 1); - assert (Dh < (uintmax_t) 1 << (W_TYPE_SIZE - 2)); - - S = isqrt2 (Dh, Dl); - - Q1 = 1; - P = S; - - /* Square root remainder fits in one word, so ignore high part. */ - Q = Dl - P*P; - /* FIXME: When can this differ from floor(sqrt(2 sqrt(D)))? */ - L = isqrt (2*S); - B = 2*L; - L1 = mu * 2 * L; - - /* The form is (+/- Q1, 2P, -/+ Q), of discriminant 4 (P^2 + Q Q1) = - 4 D. */ - - for (i = 0; i <= B; i++) - { - uintmax_t q, P1, t, rem; - - div_smallq (q, rem, S+P, Q); - P1 = S - rem; /* P1 = q*Q - P */ - -#if STAT_SQUFOF - assert (q > 0); - q_freq[0]++; - q_freq[MIN(q, Q_FREQ_SIZE)]++; -#endif - - if (Q <= L1) - { - uintmax_t g = Q; - - if ( (Q & 1) == 0) - g /= 2; - - g /= gcd_odd (g, mu); - - if (g <= L) - { - if (qpos >= QUEUE_SIZE) - { - fprintf (stderr, "squfof queue overflow.\n"); - exit (EXIT_FAILURE); - } - queue[qpos].Q = g; - queue[qpos].P = P % g; - qpos++; - } - } - - /* I think the difference can be either sign, but mod - 2^W_TYPE_SIZE arithmetic should be fine. */ - t = Q1 + q * (P - P1); - Q1 = Q; - Q = t; - P = P1; - - if ( (i & 1) == 0) - { - uintmax_t r = is_square (Q); - if (r) - { - unsigned int j; - - for (j = 0; j < qpos; j++) - { - if (queue[j].Q == r) - { - if (r == 1) - /* Traversed entire cycle. */ - goto next_multiplier; - - /* Need the absolute value for divisibility test. */ - if (P >= queue[j].P) - t = P - queue[j].P; - else - t = queue[j].P - P; - if (t % r == 0) - { - /* Delete entries up to and including entry - j, which matched. */ - memmove (queue, queue + j + 1, - (qpos - j - 1) * sizeof (queue[0])); - qpos -= (j + 1); - } - goto next_i; - } - } - - /* We have found a square form, which should give a - factor. */ - Q1 = r; - assert (S >= P); /* What signs are possible? */ - P += r * ((S - P) / r); - - /* Note: Paper says (N - P*P) / Q1, that seems incorrect - for the case D = 2N. */ - /* Compute Q = (D - P*P) / Q1, but we need double - precision. */ - uintmax_t hi, lo; - umul_ppmm (hi, lo, P, P); - sub_ddmmss (hi, lo, Dh, Dl, hi, lo); - udiv_qrnnd (Q, rem, hi, lo, Q1); - assert (rem == 0); - - for (;;) - { - /* Note: There appears to by a typo in the paper, - Step 4a in the algorithm description says q <-- - floor([S+P]/\hat Q), but looking at the equations - in Sec. 3.1, it should be q <-- floor([S+P] / Q). - (In this code, \hat Q is Q1). */ - div_smallq (q, rem, S+P, Q); - P1 = S - rem; /* P1 = q*Q - P */ - -#if STAT_SQUFOF - q_freq[0]++; - q_freq[MIN (q, Q_FREQ_SIZE)]++; -#endif - if (P == P1) - break; - t = Q1 + q * (P - P1); - Q1 = Q; - Q = t; - P = P1; - } - - if ( (Q & 1) == 0) - Q /= 2; - Q /= gcd_odd (Q, mu); - - assert (Q > 1 && (n1 || Q < n0)); - - if (prime_p (Q)) - factor_insert (factors, Q); - else if (!factor_using_squfof (0, Q, factors)) - factor_using_pollard_rho (Q, 2, factors); - - divexact_21 (n1, n0, n1, n0, Q); - - if (prime2_p (n1, n0)) - factor_insert_large (factors, n1, n0); - else - { - if (!factor_using_squfof (n1, n0, factors)) - { - if (n1 == 0) - factor_using_pollard_rho (n0, 1, factors); - else - factor_using_pollard_rho2 (n1, n0, 1, factors); - } - } - - return true; - } - } - next_i: - ; - } - next_multiplier: - ; - } - return false; -} - -static void -factor (uintmax_t t1, uintmax_t t0, struct factors *factors) -{ - factors->nfactors = 0; - factors->plarge[1] = 0; - - if (t1 == 0 && t0 < 2) - return; - - t0 = factor_using_division (&t1, t1, t0, factors); - - if (t1 == 0 && t0 < 2) - return; - - if (prime2_p (t1, t0)) - factor_insert_large (factors, t1, t0); - else - { - if (alg == ALG_SQUFOF) - if (factor_using_squfof (t1, t0, factors)) - return; - - if (t1 == 0) - factor_using_pollard_rho (t0, 1, factors); - else - factor_using_pollard_rho2 (t1, t0, 1, factors); - } -} - -#if HAVE_GMP -static void -mp_factor (mpz_t t, struct mp_factors *factors) -{ - mp_factor_init (factors); - - if (mpz_sgn (t) != 0) - { - mp_factor_using_division (t, factors); - - if (mpz_cmp_ui (t, 1) != 0) - { - if (flag_verbose > 0) - { - printf ("[is number prime?] "); - } - if (mp_prime_p (t)) - mp_factor_insert (factors, t); - else - mp_factor_using_pollard_rho (t, 1, factors); - } - } -} -#endif - -static int -strto2uintmax (uintmax_t *hip, uintmax_t *lop, const char *s) -{ - int errcode; - unsigned int lo_carry; - uintmax_t hi, lo; - - errcode = -1; - - hi = lo = 0; - for (;;) - { - unsigned int c = *s++; - if (c == 0) - break; - - if (UNLIKELY (c < '0' || c > '9')) - { - errcode = -1; - break; - } - c -= '0'; - - errcode = 0; /* we've seen at least one valid digit */ - - if (UNLIKELY (hi > ~(uintmax_t)0 / 10)) - { - errcode = -1; /* overflow */ - break; - } - hi = 10 * hi; - - lo_carry = (lo >> (W_TYPE_SIZE - 3)) + (lo >> (W_TYPE_SIZE - 1)); - lo_carry += 10 * lo < 2 * lo; - - lo = 10 * lo; - lo += c; - - lo_carry += lo < c; - hi += lo_carry; - if (UNLIKELY (hi < lo_carry)) - { - errcode = -1; /* overflow */ - break; - } - } - - *hip = hi; - *lop = lo; - - return errcode; -} - -static void -print_uintmaxes (uintmax_t t1, uintmax_t t0) -{ - uintmax_t q, r; - - if (t1 == 0) - printf ("%ju", t0); - else - { - /* Use very plain code here since it seems hard to write fast code - without assuming a specific word size. */ - q = t1 / 1000000000; - r = t1 % 1000000000; - udiv_qrnnd (t0, r, r, t0, 1000000000); - print_uintmaxes (q, t0); - printf ("%09u", (int) r); - } -} - -static void -factor_one (const char *input) -{ - uintmax_t t1, t0; - int errcode; - - /* Try converting the number to one or two words. If it fails, use GMP or - print an error message. The 2nd condition checks that the most - significant bit of the two-word number is clear, in a typesize neutral - way. */ - errcode = strto2uintmax (&t1, &t0, input); - if (errcode == 0 && ((t1 << 1) >> 1) == t1) - { - struct factors factors; - - print_uintmaxes (t1, t0); - printf (":"); - - factor (t1, t0, &factors); - - for (unsigned int j = 0; j < factors.nfactors; j++) - for (unsigned int k = 0; k < factors.e[j]; k++) - printf (" %ju", factors.p[j]); - - if (factors.plarge[1]) - { - printf (" "); - print_uintmaxes (factors.plarge[1], factors.plarge[0]); - } - puts (""); - } - else - { -#if HAVE_GMP - mpz_t t; - struct mp_factors factors; - - mpz_init_set_str (t, input, 10); - - gmp_printf ("%Zd:", t); - mp_factor (t, &factors); - - for (unsigned int j = 0; j < factors.nfactors; j++) - for (unsigned int k = 0; k < factors.e[j]; k++) - gmp_printf (" %Zd", factors.p[j]); - - mp_factor_clear (&factors); - mpz_clear (t); - puts (""); -#else - fprintf (stderr, "error: number %s not parsable or too large\n", input); - exit (EXIT_FAILURE); -#endif - } -} - -struct inbuf -{ - char *buf; - size_t alloc; -}; - -/* Read white-space delimited items. Return 1 on success, 0 on EOF. - Exit on I/O errors. */ -int -read_item (struct inbuf *bufstruct) -{ - int c; - size_t i; - char *buf = bufstruct->buf; - - do - c = getchar_unlocked (); - while (isspace (c)); - - for (i = 0; !isspace(c); i++) - { - if (c < 0) - { - if (ferror (stdin)) - { - fprintf (stderr, "read error on stdin: %s\n", strerror(errno)); - exit (EXIT_FAILURE); - } - if (i == 0) - return 0; - else - break; - } - - if (UNLIKELY (bufstruct->alloc <= i + 1)) /* +1 byte for terminating NUL */ - { - bufstruct->alloc = bufstruct->alloc * 5 / 4 + 1; - bufstruct->buf = realloc (bufstruct->buf, bufstruct->alloc); - buf = bufstruct->buf; - } - - buf[i] = c; - c = getchar_unlocked (); - } - - buf[i] = '\0'; - return 1; -} - -int -main (int argc, char *argv[]) -{ - int c; - - alg = ALG_POLLARD_RHO; /* Default to Pollard rho */ - - while ( (c = getopt(argc, argv, "svw")) != -1) - switch (c) - { - case 's': - alg = ALG_SQUFOF; - break; - case 'v': - flag_verbose = 1; - break; - case 'w': - flag_prove_primality = 0; - break; - case '?': - printf ("Usage: %s [-s] number ...\n", argv[0]); - return EXIT_FAILURE; - } -#if STAT_SQUFOF - if (alg == ALG_SQUFOF) - memset (q_freq, 0, sizeof (q_freq)); -#endif - - if (optind < argc) - for (int i = optind; i < argc; i++) - factor_one (argv[i]); - else - { - struct inbuf bufstruct; - bufstruct.alloc = 50; /* enough unless HAVE_GMP */ - bufstruct.buf = malloc (bufstruct.alloc); - while (read_item (&bufstruct)) - factor_one (bufstruct.buf); - } - -#if STAT_SQUFOF - if (alg == ALG_SQUFOF && q_freq[0] > 0) - { - double acc_f; - printf ("q freq. cum. freq.(total: %d)\n", q_freq[0]); - for (unsigned int i = 1, acc_f = 0.0; i <= Q_FREQ_SIZE; i++) - { - double f = (double) q_freq[i] / q_freq[0]; - acc_f += f; - printf ("%s%d %.2f%% %.2f%%\n", i == Q_FREQ_SIZE ? ">=" : "", i, - 100.0 * f, 100.0 * acc_f); - } - } -#endif - - exit (EXIT_SUCCESS); -} diff --git a/src/factor.c b/src/factor.c index e63e0e01d..84d9721c4 100644 --- a/src/factor.c +++ b/src/factor.c @@ -14,17 +14,79 @@ You should have received a copy of the GNU General Public License along with this program. If not, see <http://www.gnu.org/licenses/>. */ -/* Written by Paul Rubin <phr@ocf.berkeley.edu>. +/* Originally written by Paul Rubin <phr@ocf.berkeley.edu>. Adapted for GNU, fixed to factor UINT_MAX by Jim Meyering. Arbitrary-precision code adapted by James Youngman from Torbjörn Granlund's factorize.c, from GNU MP version 4.2.2. + In 2012, the core was rewritten by Torbjörn Granlund and Niels Möller. + Contains code from GNU MP. */ + +/* Efficiently factor numbers that fit in one or two words (word = uintmax_t), + or, with GMP, numbers of any size. + + Code organisation: + + There are several variants of many functions, for handling one word, two + words, and GMP's mpz_t type. If the one-word variant is called foo, the + two-word variant will be foo2, and the one for mpz_t will be mp_foo. In + some cases, the plain function variants will handle both one-word and + two-word numbers, evidenced by function arguments. + + The factoring code for two words will fall into the code for one word when + progress allows that. + + Using GMP is optional. Define HAVE_GMP to make this code include GMP + factoring code. The GMP factoring code is based on GMP's demos/factorize.c + (last synched 2012-09-07). The GMP-based factoring code will stay in GMP + factoring code even if numbers get small enough for using the two-word + code. + + Algorithm: + + (1) Perform trial division using a small primes table, but without hardware + division since the primes table store inverses modulo the word base. + (The GMP variant of this code doesn't make use of the precomputed + inverses, but instead relies on GMP for fast divisibility testing.) + (2) Check the nature of any non-factored part using Miller-Rabin for + detecting composites, and Lucas for detecting primes. + (3) Factor any remaining composite part using the Pollard-Brent rho + algorithm or the SQUFOF algorithm, checking status of found factors + again using Miller-Rabin and Lucas. + + We prefer using Hensel norm in the divisions, not the more familiar + Euclidian norm, since the former leads to much faster code. In the + Pollard-Brent rho code and the prime testing code, we use Montgomery's + trick of multiplying all n-residues by the word base, allowing cheap Hensel + reductions mod n. + + Improvements: + + * Use modular inverses also for exact division in the Lucas code, and + elsewhere. A problem is to locate the inverses not from an index, but + from a prime. We might instead compute the inverse on-the-fly. + + * Tune trial division table size (not forgetting that this is a standalone + program where the table will be read from disk for each invocation). + + * Implement less naive powm, using k-ary exponentiation for k = 3 or + perhaps k = 4. + + * Try to speed trial division code for single uintmax_t numbers, i.e., the + code using DIVBLOCK. It currently runs at 2 cycles per prime (Intel SBR, + IBR), 3 cycles per prime (AMD Stars) and 5 cycles per prime (AMD BD) when + using gcc 4.6 and 4.7. Some software pipelining should help; 1, 2, and 4 + respectively cycles ought to be possible. + + * The redcify function could be vastly improved by using (plain Euclidian) + pre-inversion (such as GMP's invert_limb) and udiv_qrnnd_preinv (from + GMP's gmp-impl.h). The redcify2 function could be vastly improved using + similar methoods. These functions currently dominate run time when using + the -w option. */ #include <config.h> #include <getopt.h> -#include <stdarg.h> #include <stdio.h> -#include <sys/types.h> #if HAVE_GMP # include <gmp.h> #endif @@ -40,17 +102,540 @@ /* The official name of this program (e.g., no 'g' prefix). */ #define PROGRAM_NAME "factor" -#define AUTHORS proper_name ("Paul Rubin") +#define AUTHORS \ + proper_name ("Paul Rubin"), \ + proper_name_utf8 ("Torbjorn Granlund", "Torbj\303\266rn Granlund"), \ + proper_name_utf8 ("Niels Moller", "Niels M\303\266ller") /* Token delimiters when reading from a file. */ #define DELIM "\n\t " -static bool verbose = false; +#ifndef STAT_SQUFOF +# define STAT_SQUFOF 0 +#endif + +#ifndef USE_LONGLONG_H +# define USE_LONGLONG_H 1 +#endif + +#if USE_LONGLONG_H + +/* Make definitions for longlong.h to make it do what it can do for us */ +# define W_TYPE_SIZE 64 /* bitcount for uintmax_t */ +# define UWtype uintmax_t +# define UHWtype unsigned long int +# undef UDWtype +# if HAVE_ATTRIBUTE_MODE +typedef unsigned int UQItype __attribute__ ((mode (QI))); +typedef int SItype __attribute__ ((mode (SI))); +typedef unsigned int USItype __attribute__ ((mode (SI))); +typedef int DItype __attribute__ ((mode (DI))); +typedef unsigned int UDItype __attribute__ ((mode (DI))); +# else +typedef unsigned char UQItype; +typedef long SItype; +typedef unsigned long int USItype; +# if HAVE_LONG_LONG +typedef long long int DItype; +typedef unsigned long long int UDItype; +# else /* Assume `long' gives us a wide enough type. Needed for hppa2.0w. */ +typedef long int DItype; +typedef unsigned long int UDItype; +# endif +# endif +# define LONGLONG_STANDALONE /* Don't require GMP's longlong.h mdep files */ +# define ASSERT(x) /* FIXME make longlong.h really standalone */ +# define __clz_tab factor_clz_tab /* Rename to avoid glibc collision */ +# ifndef __GMP_GNUC_PREREQ +# define __GMP_GNUC_PREREQ(a,b) 1 +# endif +# if _ARCH_PPC +# define HAVE_HOST_CPU_FAMILY_powerpc 1 +# endif +# include "longlong.h" +# ifdef COUNT_LEADING_ZEROS_NEED_CLZ_TAB +const unsigned char factor_clz_tab[129] = +{ + 1,2,3,3,4,4,4,4,5,5,5,5,5,5,5,5,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6, + 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7, + 8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8, + 8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8, + 9 +}; +# endif + +#else /* not USE_LONGLONG_H */ + +# define W_TYPE_SIZE (8 * sizeof (uintmax_t)) +# define __ll_B ((uintmax_t) 1 << (W_TYPE_SIZE / 2)) +# define __ll_lowpart(t) ((uintmax_t) (t) & (__ll_B - 1)) +# define __ll_highpart(t) ((uintmax_t) (t) >> (W_TYPE_SIZE / 2)) + +#endif + +enum alg_type { ALG_POLLARD_RHO = 1, ALG_SQUFOF = 2 }; + +static enum alg_type alg; + +/* 2*3*5*7*11...*101 is 128 bits, and has 26 prime factors */ +#define MAX_NFACTS 26 + +enum +{ + VERBOSE_OPTION = CHAR_MAX + 1 +}; + +static struct option const long_options[] = +{ + {"verbose", no_argument, NULL, VERBOSE_OPTION}, + {GETOPT_HELP_OPTION_DECL}, + {GETOPT_VERSION_OPTION_DECL}, + {NULL, 0, NULL, 0} +}; + +struct factors +{ + uintmax_t plarge[2]; /* Can have a single large factor */ + uintmax_t p[MAX_NFACTS]; + unsigned char e[MAX_NFACTS]; + unsigned char nfactors; +}; + +#if HAVE_GMP +struct mp_factors +{ + mpz_t *p; + unsigned long int *e; + unsigned long int nfactors; +}; +#endif + +static void factor (uintmax_t, uintmax_t, struct factors *); + +#ifndef umul_ppmm +# define umul_ppmm(w1, w0, u, v) \ + do { \ + uintmax_t __x0, __x1, __x2, __x3; \ + unsigned long int __ul, __vl, __uh, __vh; \ + uintmax_t __u = (u), __v = (v); \ + \ + __ul = __ll_lowpart (__u); \ + __uh = __ll_highpart (__u); \ + __vl = __ll_lowpart (__v); \ + __vh = __ll_highpart (__v); \ + \ + __x0 = (uintmax_t) __ul * __vl; \ + __x1 = (uintmax_t) __ul * __vh; \ + __x2 = (uintmax_t) __uh * __vl; \ + __x3 = (uintmax_t) __uh * __vh; \ + \ + __x1 += __ll_highpart (__x0);/* this can't give carry */ \ + __x1 += __x2; /* but this indeed can */ \ + if (__x1 < __x2) /* did we get it? */ \ + __x3 += __ll_B; /* yes, add it in the proper pos. */ \ + \ + (w1) = __x3 + __ll_highpart (__x1); \ + (w0) = (__x1 << W_TYPE_SIZE / 2) + __ll_lowpart (__x0); \ + } while (0) +#endif + +#if !defined udiv_qrnnd || defined UDIV_NEEDS_NORMALIZATION +/* Define our own, not needing normalization. This function is + currently not performance critical, so keep it simple. Similar to + the mod macro below. */ +# undef udiv_qrnnd +# define udiv_qrnnd(q, r, n1, n0, d) \ + do { \ + uintmax_t __d1, __d0, __q, __r1, __r0; \ + \ + assert ((n1) < (d)); \ + __d1 = (d); __d0 = 0; \ + __r1 = (n1); __r0 = (n0); \ + __q = 0; \ + for (unsigned int __i = W_TYPE_SIZE; __i > 0; __i--) \ + { \ + rsh2 (__d1, __d0, __d1, __d0, 1); \ + __q <<= 1; \ + if (ge2 (__r1, __r0, __d1, __d0)) \ + { \ + __q++; \ + sub_ddmmss (__r1, __r0, __r1, __r0, __d1, __d0); \ + } \ + } \ + (r) = __r0; \ + (q) = __q; \ + } while (0) +#endif + +#if !defined add_ssaaaa +# define add_ssaaaa(sh, sl, ah, al, bh, bl) \ + do { \ + uintmax_t _add_x; \ + _add_x = (al) + (bl); \ + (sh) = (ah) + (bh) + (_add_x < (al)); \ + (sl) = _add_x; \ + } while (0) +#endif + +#define rsh2(rh, rl, ah, al, cnt) \ + do { \ + (rl) = ((ah) << (W_TYPE_SIZE - (cnt))) | ((al) >> (cnt)); \ + (rh) = (ah) >> (cnt); \ + } while (0) + +#define lsh2(rh, rl, ah, al, cnt) \ + do { \ + (rh) = ((ah) << cnt) | ((al) >> (W_TYPE_SIZE - (cnt))); \ + (rl) = (al) << (cnt); \ + } while (0) + +#define ge2(ah, al, bh, bl) \ + ((ah) > (bh) || ((ah) == (bh) && (al) >= (bl))) + +#define gt2(ah, al, bh, bl) \ + ((ah) > (bh) || ((ah) == (bh) && (al) > (bl))) + +#ifndef sub_ddmmss +# define sub_ddmmss(rh, rl, ah, al, bh, bl) \ + do { \ + uintmax_t _cy; \ + _cy = (al) < (bl); \ + (rl) = (al) - (bl); \ + (rh) = (ah) - (bh) - _cy; \ + } while (0) +#endif + +#ifndef count_leading_zeros +# define count_leading_zeros(count, x) do { \ + uintmax_t __clz_x = (x); \ + unsigned int __clz_c; \ + for (__clz_c = 0; \ + (__clz_x & ((uintmax_t) 0xff << (W_TYPE_SIZE - 8))) == 0; \ + __clz_c += 8) \ + __clz_x <<= 8; \ + for (; (intmax_t)__clz_x >= 0; __clz_c++) \ + __clz_x <<= 1; \ + (count) = __clz_c; \ + } while (0) +#endif + +#ifndef count_trailing_zeros +# define count_trailing_zeros(count, x) do { \ + uintmax_t __ctz_x = (x); \ + unsigned int __ctz_c = 0; \ + while ((__ctz_x & 1) == 0) \ + { \ + __ctz_x >>= 1; \ + __ctz_c++; \ + } \ + (count) = __ctz_c; \ + } while (0) +#endif + +/* Requires that a < n and b <= n */ +#define submod(r,a,b,n) \ + do { \ + uintmax_t _t = - (uintmax_t) (a < b); \ + (r) = ((n) & _t) + (a) - (b); \ + } while (0) + +#define addmod(r,a,b,n) \ + submod ((r), (a), ((n) - (b)), (n)) + +/* Modular two-word addition and subtraction. For performance reasons, the + most significant bit of n1 must be clear. The destination variables must be + distinct from the mod operand. */ +#define addmod2(r1, r0, a1, a0, b1, b0, n1, n0) \ + do { \ + add_ssaaaa ((r1), (r0), (a1), (a0), (b1), (b0)); \ + if (ge2 ((r1), (r0), (n1), (n0))) \ + sub_ddmmss ((r1), (r0), (r1), (r0), (n1), (n0)); \ + } while (0) +#define submod2(r1, r0, a1, a0, b1, b0, n1, n0) \ + do { \ + sub_ddmmss ((r1), (r0), (a1), (a0), (b1), (b0)); \ + if ((intmax_t) (r1) < 0) \ + add_ssaaaa ((r1), (r0), (r1), (r0), (n1), (n0)); \ + } while (0) + +#define HIGHBIT_TO_MASK(x) \ + (((intmax_t)-1 >> 1) < 0 \ + ? (uintmax_t)((intmax_t)(x) >> (W_TYPE_SIZE - 1)) \ + : ((x) & ((uintmax_t) 1 << (W_TYPE_SIZE - 1)) \ + ? UINTMAX_MAX : (uintmax_t) 0)) + +/* Compute r = a mod d, where r = <*t1,retval>, a = <a1,a0>, d = <d1,d0>. + Requires that d1 != 0. */ +static uintmax_t +mod2 (uintmax_t *r1, uintmax_t a1, uintmax_t a0, uintmax_t d1, uintmax_t d0) +{ + int cntd, cnta; + + assert (d1 != 0); + + if (a1 == 0) + { + *r1 = 0; + return a0; + } + + count_leading_zeros (cntd, d1); + count_leading_zeros (cnta, a1); + int cnt = cntd - cnta; + lsh2 (d1, d0, d1, d0, cnt); + for (int i = 0; i < cnt; i++) + { + if (ge2 (a1, a0, d1, d0)) + sub_ddmmss (a1, a0, a1, a0, d1, d0); + rsh2 (d1, d0, d1, d0, 1); + } + + *r1 = a1; + return a0; +} + +static uintmax_t _GL_ATTRIBUTE_CONST +gcd_odd (uintmax_t a, uintmax_t b) +{ + if ( (b & 1) == 0) + { + uintmax_t t = b; + b = a; + a = t; + } + if (a == 0) + return b; + + /* Take out least significant one bit, to make room for sign */ + b >>= 1; + + for (;;) + { + uintmax_t t; + uintmax_t bgta; + + while ((a & 1) == 0) + a >>= 1; + a >>= 1; + + t = a - b; + if (t == 0) + return (a << 1) + 1; + + bgta = HIGHBIT_TO_MASK (t); + + /* b <-- min (a, b) */ + b += (bgta & t); + + /* a <-- |a - b| */ + a = (t ^ bgta) - bgta; + } +} + +static uintmax_t +gcd2_odd (uintmax_t *r1, uintmax_t a1, uintmax_t a0, uintmax_t b1, uintmax_t b0) +{ + while ((a0 & 1) == 0) + rsh2 (a1, a0, a1, a0, 1); + while ((b0 & 1) == 0) + rsh2 (b1, b0, b1, b0, 1); + + for (;;) + { + if ((b1 | a1) == 0) + { + *r1 = 0; + return gcd_odd (b0, a0); + } + + if (gt2 (a1, a0, b1, b0)) + { + sub_ddmmss (a1, a0, a1, a0, b1, b0); + do + rsh2 (a1, a0, a1, a0, 1); + while ((a0 & 1) == 0); + } + else if (gt2 (b1, b0, a1, a0)) + { + sub_ddmmss (b1, b0, b1, b0, a1, a0); + do + rsh2 (b1, b0, b1, b0, 1); + while ((b0 & 1) == 0); + } + else + break; + } + + *r1 = a1; + return a0; +} + +static void +factor_insert_multiplicity (struct factors *factors, + uintmax_t prime, unsigned int m) +{ + unsigned int nfactors = factors->nfactors; + uintmax_t *p = factors->p; + unsigned char *e = factors->e; + + /* Locate position for insert new or increment e. */ + int i; + for (i = nfactors - 1; i >= 0; i--) + { + if (p[i] <= prime) + break; + } + + if (i < 0 || p[i] != prime) + { + for (int j = nfactors - 1; j > i; j--) + { + p[j + 1] = p[j]; + e[j + 1] = e[j]; + } + p[i + 1] = prime; + e[i + 1] = m; + factors->nfactors = nfactors + 1; + } + else + { + e[i] += m; + } +} + +#define factor_insert(f, p) factor_insert_multiplicity (f, p, 1) + +static void +factor_insert_large (struct factors *factors, + uintmax_t p1, uintmax_t p0) +{ + if (p1 > 0) + { + assert (factors->plarge[1] == 0); + factors->plarge[0] = p0; + factors->plarge[1] = p1; + } + else + factor_insert (factors, p0); +} #if HAVE_GMP -static mpz_t *factor = NULL; -static size_t nfactors_found = 0; -static size_t nfactors_allocated = 0; +static void mp_factor (mpz_t, struct mp_factors *); + +static void +mp_factor_init (struct mp_factors *factors) +{ + factors->p = xmalloc (1); + factors->e = xmalloc (1); + factors->nfactors = 0; +} + +static void +mp_factor_clear (struct mp_factors *factors) +{ + for (unsigned int i = 0; i < factors->nfactors; i++) + mpz_clear (factors->p[i]); + + free (factors->p); + free (factors->e); +} + +static void +mp_factor_insert (struct mp_factors *factors, mpz_t prime) +{ + unsigned long int nfactors = factors->nfactors; + mpz_t *p = factors->p; + unsigned long int *e = factors->e; + long i; + + /* Locate position for insert new or increment e. */ + for (i = nfactors - 1; i >= 0; i--) + { + if (mpz_cmp (p[i], prime) <= 0) + break; + } + + if (i < 0 || mpz_cmp (p[i], prime) != 0) + { + p = xrealloc (p, (nfactors + 1) * sizeof p[0]); + e = xrealloc (e, (nfactors + 1) * sizeof e[0]); + + mpz_init (p[nfactors]); + for (long j = nfactors - 1; j > i; j--) + { + mpz_set (p[j + 1], p[j]); + e[j + 1] = e[j]; + } + mpz_set (p[i + 1], prime); + e[i + 1] = 1; + + factors->p = p; + factors->e = e; + factors->nfactors = nfactors + 1; + } + else + { + e[i] += 1; + } +} + +static void +mp_factor_insert_ui (struct mp_factors *factors, unsigned long int prime) +{ + mpz_t pz; + + mpz_init_set_ui (pz, prime); + mp_factor_insert (factors, pz); + mpz_clear (pz); +} +#endif /* HAVE_GMP */ + + +#define P(a,b,c,d) a, +static const unsigned char primes_diff[] = { +#include "primes.h" +0,0,0,0,0,0,0 /* 7 sentinels for 8-way loop */ +}; +#undef P + +#define PRIMES_PTAB_ENTRIES \ + (sizeof (primes_diff) / sizeof (primes_diff[0]) - 8 + 1) + +#define P(a,b,c,d) b, +static const unsigned char primes_diff8[] = { +#include "primes.h" +0,0,0,0,0,0,0 /* 7 sentinels for 8-way loop */ +}; +#undef P + +struct primes_dtab +{ + uintmax_t binv, lim; +}; + +#define P(a,b,c,d) {c,d}, +static const struct primes_dtab primes_dtab[] = { +#include "primes.h" +{1,0},{1,0},{1,0},{1,0},{1,0},{1,0},{1,0} /* 7 sentinels for 8-way loop */ +}; +#undef P + +/* This flag is honored only in the GMP code. */ +static int verbose = 0; + +/* Prove primality or run probabilistic tests. */ +static bool flag_prove_primality = true; + +/* Number of Miller-Rabin tests to run when not proving primality. */ +#define MR_REPS 25 + +#ifdef __GNUC__ +# define LIKELY(cond) __builtin_expect ((cond), 1) +# define UNLIKELY(cond) __builtin_expect ((cond), 0) +#else +# define LIKELY(cond) (cond) +# define UNLIKELY(cond) (cond) +#endif static void debug (char const *fmt, ...) @@ -65,368 +650,1637 @@ debug (char const *fmt, ...) } static void -emit_factor (mpz_t n) +factor_insert_refind (struct factors *factors, uintmax_t p, unsigned int i, + unsigned int off) { - if (nfactors_found == nfactors_allocated) - factor = X2NREALLOC (factor, &nfactors_allocated); - mpz_init (factor[nfactors_found]); - mpz_set (factor[nfactors_found], n); - ++nfactors_found; + for (unsigned int j = 0; j < off; j++) + p += primes_diff[i + j]; + factor_insert (factors, p); } -static void -emit_ul_factor (unsigned long int i) +/* Trial division with odd primes uses the following trick. + + Let p be an odd prime, and B = 2^{W_TYPE_SIZE}. For simplicity, + consider the case t < B (this is the second loop below). + + From our tables we get + + binv = p^{-1} (mod B) + lim = floor ( (B-1) / p ). + + First assume that t is a multiple of p, t = q * p. Then 0 <= q <= lim + (and all quotients in this range occur for some t). + + Then t = q * p is true also (mod B), and p is invertible we get + + q = t * binv (mod B). + + Next, assume that t is *not* divisible by p. Since multiplication + by binv (mod B) is a one-to-one mapping, + + t * binv (mod B) > lim, + + because all the smaller values are already taken. + + This can be summed up by saying that the function + + q(t) = binv * t (mod B) + + is a permutation of the range 0 <= t < B, with the curious property + that it maps the multiples of p onto the range 0 <= q <= lim, in + order, and the non-multiples of p onto the range lim < q < B. + */ + +static uintmax_t +factor_using_division (uintmax_t *t1p, uintmax_t t1, uintmax_t t0, + struct factors *factors) { - mpz_t t; - mpz_init (t); - mpz_set_ui (t, i); - emit_factor (t); - mpz_clear (t); + if (t0 % 2 == 0) + { + unsigned int cnt; + + if (t0 == 0) + { + count_trailing_zeros (cnt, t1); + t0 = t1 >> cnt; + t1 = 0; + cnt += W_TYPE_SIZE; + } + else + { + count_trailing_zeros (cnt, t0); + rsh2 (t1, t0, t1, t0, cnt); + } + + factor_insert_multiplicity (factors, 2, cnt); + } + + uintmax_t p = 3; + unsigned int i; + for (i = 0; t1 > 0 && i < PRIMES_PTAB_ENTRIES; i++) + { + for (;;) + { + uintmax_t q1, q0, hi, lo; + + q0 = t0 * primes_dtab[i].binv; + umul_ppmm (hi, lo, q0, p); + if (hi > t1) + break; + hi = t1 - hi; + q1 = hi * primes_dtab[i].binv; + if (LIKELY (q1 > primes_dtab[i].lim)) + break; + t1 = q1; t0 = q0; + factor_insert (factors, p); + } + p += primes_diff[i + 1]; + } + if (t1p) + *t1p = t1; + +#define DIVBLOCK(I) \ + do { \ + for (;;) \ + { \ + q = t0 * pd[I].binv; \ + if (LIKELY (q > pd[I].lim)) \ + break; \ + t0 = q; \ + factor_insert_refind (factors, p, i + 1, I); \ + } \ + } while (0) + + for (; i < PRIMES_PTAB_ENTRIES; i += 8) + { + uintmax_t q; + const struct primes_dtab *pd = &primes_dtab[i]; + DIVBLOCK (0); + DIVBLOCK (1); + DIVBLOCK (2); + DIVBLOCK (3); + DIVBLOCK (4); + DIVBLOCK (5); + DIVBLOCK (6); + DIVBLOCK (7); + + p += primes_diff8[i]; + if (p * p > t0) + break; + } + + return t0; } +#if HAVE_GMP static void -factor_using_division (mpz_t t, unsigned int limit) +mp_factor_using_division (mpz_t t, struct mp_factors *factors) { - mpz_t q, r; - unsigned long int f; - int ai; - static unsigned int const add[] = {4, 2, 4, 2, 4, 6, 2, 6}; - unsigned int const *addv = add; - unsigned int failures; + mpz_t q; + unsigned long int p; - debug ("[trial division (%u)] ", limit); + debug ("[trial division] "); mpz_init (q); - mpz_init (r); - f = mpz_scan1 (t, 0); - mpz_div_2exp (t, t, f); - while (f) + p = mpz_scan1 (t, 0); + mpz_div_2exp (t, t, p); + while (p) { - emit_ul_factor (2); - --f; + mp_factor_insert_ui (factors, 2); + --p; } - while (true) + p = 3; + for (unsigned int i = 1; i <= PRIMES_PTAB_ENTRIES;) { - mpz_tdiv_qr_ui (q, r, t, 3); - if (mpz_cmp_ui (r, 0) != 0) - break; - mpz_set (t, q); - emit_ul_factor (3); + if (! mpz_divisible_ui_p (t, p)) + { + p += primes_diff[i++]; + if (mpz_cmp_ui (t, p * p) < 0) + break; + } + else + { + mpz_tdiv_q_ui (t, t, p); + mp_factor_insert_ui (factors, p); + } } - while (true) + mpz_clear (q); +} +#endif + +/* Entry i contains (2i+1)^(-1) mod 2^8. */ +static const unsigned char binvert_table[128] = +{ + 0x01, 0xAB, 0xCD, 0xB7, 0x39, 0xA3, 0xC5, 0xEF, + 0xF1, 0x1B, 0x3D, 0xA7, 0x29, 0x13, 0x35, 0xDF, + 0xE1, 0x8B, 0xAD, 0x97, 0x19, 0x83, 0xA5, 0xCF, + 0xD1, 0xFB, 0x1D, 0x87, 0x09, 0xF3, 0x15, 0xBF, + 0xC1, 0x6B, 0x8D, 0x77, 0xF9, 0x63, 0x85, 0xAF, + 0xB1, 0xDB, 0xFD, 0x67, 0xE9, 0xD3, 0xF5, 0x9F, + 0xA1, 0x4B, 0x6D, 0x57, 0xD9, 0x43, 0x65, 0x8F, + 0x91, 0xBB, 0xDD, 0x47, 0xC9, 0xB3, 0xD5, 0x7F, + 0x81, 0x2B, 0x4D, 0x37, 0xB9, 0x23, 0x45, 0x6F, + 0x71, 0x9B, 0xBD, 0x27, 0xA9, 0x93, 0xB5, 0x5F, + 0x61, 0x0B, 0x2D, 0x17, 0x99, 0x03, 0x25, 0x4F, + 0x51, 0x7B, 0x9D, 0x07, 0x89, 0x73, 0x95, 0x3F, + 0x41, 0xEB, 0x0D, 0xF7, 0x79, 0xE3, 0x05, 0x2F, + 0x31, 0x5B, 0x7D, 0xE7, 0x69, 0x53, 0x75, 0x1F, + 0x21, 0xCB, 0xED, 0xD7, 0x59, 0xC3, 0xE5, 0x0F, + 0x11, 0x3B, 0x5D, 0xC7, 0x49, 0x33, 0x55, 0xFF +}; + +/* Compute n^(-1) mod B, using a Newton iteration. */ +#define binv(inv,n) \ + do { \ + uintmax_t __n = (n); \ + uintmax_t __inv; \ + \ + __inv = binvert_table[(__n / 2) & 0x7F]; /* 8 */ \ + if (W_TYPE_SIZE > 8) __inv = 2 * __inv - __inv * __inv * __n; \ + if (W_TYPE_SIZE > 16) __inv = 2 * __inv - __inv * __inv * __n; \ + if (W_TYPE_SIZE > 32) __inv = 2 * __inv - __inv * __inv * __n; \ + \ + if (W_TYPE_SIZE > 64) \ + { \ + int __invbits = 64; \ + do { \ + __inv = 2 * __inv - __inv * __inv * __n; \ + __invbits *= 2; \ + } while (__invbits < W_TYPE_SIZE); \ + } \ + \ + (inv) = __inv; \ + } while (0) + +/* q = u / d, assuming d|u. */ +#define divexact_21(q1, q0, u1, u0, d) \ + do { \ + uintmax_t _di, _q0; \ + binv (_di, (d)); \ + _q0 = (u0) * _di; \ + if ((u1) >= (d)) \ + { \ + uintmax_t _p1, _p0; \ + umul_ppmm (_p1, _p0, _q0, d); \ + (q1) = ((u1) - _p1) * _di; \ + (q0) = _q0; \ + } \ + else \ + { \ + (q0) = _q0; \ + (q1) = 0; \ + } \ + } while (0) + +/* x B (mod n). */ +#define redcify(r_prim, r, n) \ + do { \ + uintmax_t _redcify_q ATTRIBUTE_UNUSED; \ + udiv_qrnnd (_redcify_q, r_prim, r, 0, n); \ + } while (0) + +/* x B^2 (mod n). Requires x > 0, n1 < B/2 */ +#define redcify2(r1, r0, x, n1, n0) \ + do { \ + uintmax_t _r1, _r0, _i; \ + if ((x) < (n1)) \ + { \ + _r1 = (x); _r0 = 0; \ + _i = W_TYPE_SIZE; \ + } \ + else \ + { \ + _r1 = 0; _r0 = (x); \ + _i = 2*W_TYPE_SIZE; \ + } \ + while (_i-- > 0) \ + { \ + lsh2 (_r1, _r0, _r1, _r0, 1); \ + if (ge2 (_r1, _r0, (n1), (n0))) \ + sub_ddmmss (_r1, _r0, _r1, _r0, (n1), (n0)); \ + } \ + (r1) = _r1; \ + (r0) = _r0; \ + } while (0) + +/* Modular two-word multiplication, r = a * b mod m, with mi = m^(-1) mod B. + Both a and b must be in redc form, the result will be in redc form too. */ +static inline uintmax_t +mulredc (uintmax_t a, uintmax_t b, uintmax_t m, uintmax_t mi) +{ + uintmax_t rh, rl, q, th, tl, xh; + + umul_ppmm (rh, rl, a, b); + q = rl * mi; + umul_ppmm (th, tl, q, m); + xh = rh - th; + if (rh < th) + xh += m; + + return xh; +} + +/* Modular two-word multiplication, r = a * b mod m, with mi = m^(-1) mod B. + Both a and b must be in redc form, the result will be in redc form too. + For performance reasons, the most significant bit of m must be clear. */ +static uintmax_t +mulredc2 (uintmax_t *r1p, + uintmax_t a1, uintmax_t a0, uintmax_t b1, uintmax_t b0, + uintmax_t m1, uintmax_t m0, uintmax_t mi) +{ + uintmax_t r1, r0, q, p1, p0, t1, t0, s1, s0; + mi = -mi; + assert ( (a1 >> (W_TYPE_SIZE - 1)) == 0); + assert ( (b1 >> (W_TYPE_SIZE - 1)) == 0); + assert ( (m1 >> (W_TYPE_SIZE - 1)) == 0); + + /* First compute a0 * <b1, b0> B^{-1} + +-----+ + |a0 b0| + +--+--+--+ + |a0 b1| + +--+--+--+ + |q0 m0| + +--+--+--+ + |q0 m1| + -+--+--+--+ + |r1|r0| 0| + +--+--+--+ + */ + umul_ppmm (t1, t0, a0, b0); + umul_ppmm (r1, r0, a0, b1); + q = mi * t0; + umul_ppmm (p1, p0, q, m0); + umul_ppmm (s1, s0, q, m1); + r0 += (t0 != 0); /* Carry */ + add_ssaaaa (r1, r0, r1, r0, 0, p1); + add_ssaaaa (r1, r0, r1, r0, 0, t1); + add_ssaaaa (r1, r0, r1, r0, s1, s0); + + /* Next, (a1 * <b1, b0> + <r1, r0> B^{-1} + +-----+ + |a1 b0| + +--+--+ + |r1|r0| + +--+--+--+ + |a1 b1| + +--+--+--+ + |q1 m0| + +--+--+--+ + |q1 m1| + -+--+--+--+ + |r1|r0| 0| + +--+--+--+ + */ + umul_ppmm (t1, t0, a1, b0); + umul_ppmm (s1, s0, a1, b1); + add_ssaaaa (t1, t0, t1, t0, 0, r0); + q = mi * t0; + add_ssaaaa (r1, r0, s1, s0, 0, r1); + umul_ppmm (p1, p0, q, m0); + umul_ppmm (s1, s0, q, m1); + r0 += (t0 != 0); /* Carry */ + add_ssaaaa (r1, r0, r1, r0, 0, p1); + add_ssaaaa (r1, r0, r1, r0, 0, t1); + add_ssaaaa (r1, r0, r1, r0, s1, s0); + + if (ge2 (r1, r0, m1, m0)) + sub_ddmmss (r1, r0, r1, r0, m1, m0); + + *r1p = r1; + return r0; +} + +static uintmax_t _GL_ATTRIBUTE_CONST +powm (uintmax_t b, uintmax_t e, uintmax_t n, uintmax_t ni, uintmax_t one) +{ + uintmax_t y = one; + + if (e & 1) + y = b; + + while (e != 0) { - mpz_tdiv_qr_ui (q, r, t, 5); - if (mpz_cmp_ui (r, 0) != 0) - break; - mpz_set (t, q); - emit_ul_factor (5); + b = mulredc (b, b, n, ni); + e >>= 1; + + if (e & 1) + y = mulredc (y, b, n, ni); } - failures = 0; - f = 7; - ai = 0; - while (mpz_cmp_ui (t, 1) != 0) + return y; +} + +static uintmax_t +powm2 (uintmax_t *r1m, + const uintmax_t *bp, const uintmax_t *ep, const uintmax_t *np, + uintmax_t ni, const uintmax_t *one) +{ + uintmax_t r1, r0, b1, b0, n1, n0; + unsigned int i; + uintmax_t e; + + b0 = bp[0]; + b1 = bp[1]; + n0 = np[0]; + n1 = np[1]; + + r0 = one[0]; + r1 = one[1]; + + for (e = ep[0], i = W_TYPE_SIZE; i > 0; i--, e >>= 1) { - mpz_tdiv_qr_ui (q, r, t, f); - if (mpz_cmp_ui (r, 0) != 0) + if (e & 1) { - f += addv[ai]; - if (mpz_cmp_ui (q, f) < 0) - break; - ai = (ai + 1) & 7; - failures++; - if (failures > limit) - break; + r0 = mulredc2 (r1m, r1, r0, b1, b0, n1, n0, ni); + r1 = *r1m; + } + b0 = mulredc2 (r1m, b1, b0, b1, b0, n1, n0, ni); + b1 = *r1m; + } + for (e = ep[1]; e > 0; e >>= 1) + { + if (e & 1) + { + r0 = mulredc2 (r1m, r1, r0, b1, b0, n1, n0, ni); + r1 = *r1m; + } + b0 = mulredc2 (r1m, b1, b0, b1, b0, n1, n0, ni); + b1 = *r1m; + } + *r1m = r1; + return r0; +} + +static bool _GL_ATTRIBUTE_CONST +millerrabin (uintmax_t n, uintmax_t ni, uintmax_t b, uintmax_t q, + unsigned int k, uintmax_t one) +{ + uintmax_t y = powm (b, q, n, ni, one); + + uintmax_t nm1 = n - one; /* -1, but in redc representation. */ + + if (y == one || y == nm1) + return true; + + for (unsigned int i = 1; i < k; i++) + { + y = mulredc (y, y, n, ni); + + if (y == nm1) + return true; + if (y == one) + return false; + } + return false; +} + +static bool +millerrabin2 (const uintmax_t *np, uintmax_t ni, const uintmax_t *bp, + const uintmax_t *qp, unsigned int k, const uintmax_t *one) +{ + uintmax_t y1, y0, nm1_1, nm1_0, r1m; + + y0 = powm2 (&r1m, bp, qp, np, ni, one); + y1 = r1m; + + if (y0 == one[0] && y1 == one[1]) + return true; + + sub_ddmmss (nm1_1, nm1_0, np[1], np[0], one[1], one[0]); + + if (y0 == nm1_0 && y1 == nm1_1) + return true; + + for (unsigned int i = 1; i < k; i++) + { + y0 = mulredc2 (&r1m, y1, y0, y1, y0, np[1], np[0], ni); + y1 = r1m; + + if (y0 == nm1_0 && y1 == nm1_1) + return true; + if (y0 == one[0] && y1 == one[1]) + return false; + } + return false; +} + +#if HAVE_GMP +static bool +mp_millerrabin (mpz_srcptr n, mpz_srcptr nm1, mpz_ptr x, mpz_ptr y, + mpz_srcptr q, unsigned long int k) +{ + mpz_powm (y, x, q, n); + + if (mpz_cmp_ui (y, 1) == 0 || mpz_cmp (y, nm1) == 0) + return true; + + for (unsigned long int i = 1; i < k; i++) + { + mpz_powm_ui (y, y, 2, n); + if (mpz_cmp (y, nm1) == 0) + return true; + if (mpz_cmp_ui (y, 1) == 0) + return false; + } + return false; +} +#endif + +/* Lucas' prime test. The number of iterations vary greatly, up to a few dozen + have been observed. The average seem to be about 2. */ +static bool +prime_p (uintmax_t n) +{ + int k; + bool is_prime; + uintmax_t a_prim, one, ni; + struct factors factors; + + if (n <= 1) + return false; + + /* We have already casted out small primes. */ + if (n < (uintmax_t) FIRST_OMITTED_PRIME * FIRST_OMITTED_PRIME) + return true; + + /* Precomputation for Miller-Rabin. */ + uintmax_t q = n - 1; + for (k = 0; (q & 1) == 0; k++) + q >>= 1; + + uintmax_t a = 2; + binv (ni, n); /* ni <- 1/n mod B */ + redcify (one, 1, n); + addmod (a_prim, one, one, n); /* i.e., redcify a = 2 */ + + /* Perform a Miller-Rabin test, finds most composites quickly. */ + if (!millerrabin (n, ni, a_prim, q, k, one)) + return false; + + if (flag_prove_primality) + { + /* Factor n-1 for Lucas. */ + factor (0, n - 1, &factors); + } + + /* Loop until Lucas proves our number prime, or Miller-Rabin proves our + number composite. */ + for (unsigned int r = 0; r < PRIMES_PTAB_ENTRIES; r++) + { + if (flag_prove_primality) + { + is_prime = true; + for (unsigned int i = 0; i < factors.nfactors && is_prime; i++) + { + is_prime + = powm (a_prim, (n - 1) / factors.p[i], n, ni, one) != one; + } } else { - mpz_swap (t, q); - emit_ul_factor (f); - failures = 0; + /* After enough Miller-Rabin runs, be content. */ + is_prime = (r == MR_REPS - 1); } + + if (is_prime) + return true; + + a += primes_diff[r]; /* Establish new base. */ + + /* The following is equivalent to redcify (a_prim, a, n). It runs faster + on most processors, since it avoids udiv_qrnnd. If we go down the + udiv_qrnnd_preinv path, this code should be replaced. */ + { + uintmax_t s1, s0; + umul_ppmm (s1, s0, one, a); + if (LIKELY (s1 == 0)) + a_prim = s0 % n; + else + { + uintmax_t dummy ATTRIBUTE_UNUSED; + udiv_qrnnd (dummy, a_prim, s1, s0, n); + } + } + + if (!millerrabin (n, ni, a_prim, q, k, one)) + return false; } - mpz_clear (q); - mpz_clear (r); + error (0, 0, _("Lucas prime test failure. This should not happen")); + abort (); } -/* The number of Miller-Rabin tests we require. */ -enum { MR_REPS = 25 }; +static bool +prime2_p (uintmax_t n1, uintmax_t n0) +{ + uintmax_t q[2], nm1[2]; + uintmax_t a_prim[2]; + uintmax_t one[2]; + uintmax_t na[2]; + uintmax_t ni; + unsigned int k; + struct factors factors; + + if (n1 == 0) + return prime_p (n0); + + nm1[1] = n1 - (n0 == 0); + nm1[0] = n0 - 1; + if (nm1[0] == 0) + { + count_trailing_zeros (k, nm1[1]); -static void -factor_using_pollard_rho (mpz_t n, int a_int) + q[0] = nm1[1] >> k; + q[1] = 0; + k += W_TYPE_SIZE; + } + else + { + count_trailing_zeros (k, nm1[0]); + rsh2 (q[1], q[0], nm1[1], nm1[0], k); + } + + uintmax_t a = 2; + binv (ni, n0); + redcify2 (one[1], one[0], 1, n1, n0); + addmod2 (a_prim[1], a_prim[0], one[1], one[0], one[1], one[0], n1, n0); + + /* FIXME: Use scalars or pointers in arguments? Some consistency needed. */ + na[0] = n0; + na[1] = n1; + + if (!millerrabin2 (na, ni, a_prim, q, k, one)) + return false; + + if (flag_prove_primality) + { + /* Factor n-1 for Lucas. */ + factor (nm1[1], nm1[0], &factors); + } + + /* Loop until Lucas proves our number prime, or Miller-Rabin proves our + number composite. */ + for (unsigned int r = 0; r < PRIMES_PTAB_ENTRIES; r++) + { + bool is_prime; + uintmax_t e[2], y[2]; + + if (flag_prove_primality) + { + is_prime = true; + if (factors.plarge[1]) + { + uintmax_t pi; + binv (pi, factors.plarge[0]); + e[0] = pi * nm1[0]; + e[1] = 0; + y[0] = powm2 (&y[1], a_prim, e, na, ni, one); + is_prime = (y[0] != one[0] || y[1] != one[1]); + } + for (unsigned int i = 0; i < factors.nfactors && is_prime; i++) + { + /* FIXME: We always have the factor 2. Do we really need to + handle it here? We have done the same powering as part + of millerrabin. */ + if (factors.p[i] == 2) + rsh2 (e[1], e[0], nm1[1], nm1[0], 1); + else + divexact_21 (e[1], e[0], nm1[1], nm1[0], factors.p[i]); + y[0] = powm2 (&y[1], a_prim, e, na, ni, one); + is_prime = (y[0] != one[0] || y[1] != one[1]); + } + } + else + { + /* After enough Miller-Rabin runs, be content. */ + is_prime = (r == MR_REPS - 1); + } + + if (is_prime) + return true; + + a += primes_diff[r]; /* Establish new base. */ + redcify2 (a_prim[1], a_prim[0], a, n1, n0); + + if (!millerrabin2 (na, ni, a_prim, q, k, one)) + return false; + } + + error (0, 0, _("Lucas prime test failure. This should not happen")); + abort (); +} + +#if HAVE_GMP +static bool +mp_prime_p (mpz_t n) { - mpz_t x, x1, y, P; - mpz_t a; - mpz_t g; - mpz_t t1, t2; - int k, l, c; + bool is_prime; + mpz_t q, a, nm1, tmp; + struct mp_factors factors; - debug ("[pollard-rho (%d)] ", a_int); + if (mpz_cmp_ui (n, 1) <= 0) + return false; - mpz_init (g); - mpz_init (t1); - mpz_init (t2); + /* We have already casted out small primes. */ + if (mpz_cmp_ui (n, (long) FIRST_OMITTED_PRIME * FIRST_OMITTED_PRIME) < 0) + return true; - mpz_init_set_si (a, a_int); - mpz_init_set_si (y, 2); - mpz_init_set_si (x, 2); - mpz_init_set_si (x1, 2); - k = 1; - l = 1; - mpz_init_set_ui (P, 1); - c = 0; + mpz_inits (q, a, nm1, tmp, NULL); - while (mpz_cmp_ui (n, 1) != 0) + /* Precomputation for Miller-Rabin. */ + mpz_sub_ui (nm1, n, 1); + + /* Find q and k, where q is odd and n = 1 + 2**k * q. */ + unsigned long int k = mpz_scan1 (nm1, 0); + mpz_tdiv_q_2exp (q, nm1, k); + + mpz_set_ui (a, 2); + + /* Perform a Miller-Rabin test, finds most composites quickly. */ + if (!mp_millerrabin (n, nm1, a, tmp, q, k)) { -S2: - mpz_mul (x, x, x); mpz_add (x, x, a); mpz_mod (x, x, n); + is_prime = false; + goto ret2; + } + + if (flag_prove_primality) + { + /* Factor n-1 for Lucas. */ + mpz_set (tmp, nm1); + mp_factor (tmp, &factors); + } - mpz_sub (t1, x1, x); mpz_mul (t2, P, t1); mpz_mod (P, t2, n); - c++; - if (c == 20) + /* Loop until Lucas proves our number prime, or Miller-Rabin proves our + number composite. */ + for (unsigned int r = 0; r < PRIMES_PTAB_ENTRIES; r++) + { + if (flag_prove_primality) { - c = 0; - mpz_gcd (g, P, n); - if (mpz_cmp_ui (g, 1) != 0) - goto S4; - mpz_set (y, x); + is_prime = true; + for (unsigned long int i = 0; i < factors.nfactors && is_prime; i++) + { + mpz_divexact (tmp, nm1, factors.p[i]); + mpz_powm (tmp, a, tmp, n); + is_prime = mpz_cmp_ui (tmp, 1) != 0; + } + } + else + { + /* After enough Miller-Rabin runs, be content. */ + is_prime = (r == MR_REPS - 1); } - k--; - if (k > 0) - goto S2; + if (is_prime) + goto ret1; - mpz_gcd (g, P, n); - if (mpz_cmp_ui (g, 1) != 0) - goto S4; + mpz_add_ui (a, a, primes_diff[r]); /* Establish new base. */ - mpz_set (x1, x); - k = l; - l = 2 * l; - unsigned int i; - for (i = 0; i < k; i++) + if (!mp_millerrabin (n, nm1, a, tmp, q, k)) { - mpz_mul (x, x, x); mpz_add (x, x, a); mpz_mod (x, x, n); + is_prime = false; + goto ret1; } - mpz_set (y, x); - c = 0; - goto S2; -S4: + } + + error (0, 0, _("Lucas prime test failure. This should not happen")); + abort (); + + ret1: + if (flag_prove_primality) + mp_factor_clear (&factors); + ret2: + mpz_clears (q, a, nm1, tmp, NULL); + + return is_prime; +} +#endif + +static void +factor_using_pollard_rho (uintmax_t n, unsigned long int a, + struct factors *factors) +{ + uintmax_t x, z, y, P, t, ni, g; + + unsigned long int k = 1; + unsigned long int l = 1; + + redcify (P, 1, n); + addmod (x, P, P, n); /* i.e., redcify(2) */ + y = z = x; + + while (n != 1) + { + assert (a < n); + + binv (ni, n); /* FIXME: when could we use old 'ni' value? */ + + for (;;) + { + do + { + x = mulredc (x, x, n, ni); + addmod (x, x, a, n); + + submod (t, z, x, n); + P = mulredc (P, t, n, ni); + + if (k % 32 == 1) + { + if (gcd_odd (P, n) != 1) + goto factor_found; + y = x; + } + } + while (--k != 0); + + z = x; + k = l; + l = 2 * l; + for (unsigned long int i = 0; i < k; i++) + { + x = mulredc (x, x, n, ni); + addmod (x, x, a, n); + } + y = x; + } + + factor_found: do { - mpz_mul (y, y, y); mpz_add (y, y, a); mpz_mod (y, y, n); - mpz_sub (t1, x1, y); mpz_gcd (g, t1, n); + y = mulredc (y, y, n, ni); + addmod (y, y, a, n); + + submod (t, z, y, n); + g = gcd_odd (t, n); + } + while (g == 1); + + n = n / g; + + if (!prime_p (g)) + factor_using_pollard_rho (g, a + 1, factors); + else + factor_insert (factors, g); + + if (prime_p (n)) + { + factor_insert (factors, n); + break; } - while (mpz_cmp_ui (g, 1) == 0); - mpz_div (n, n, g); /* divide by g, before g is overwritten */ + x = x % n; + z = z % n; + y = y % n; + } +} + +static void +factor_using_pollard_rho2 (uintmax_t n1, uintmax_t n0, unsigned long int a, + struct factors *factors) +{ + uintmax_t x1, x0, z1, z0, y1, y0, P1, P0, t1, t0, ni, g1, g0, r1m; + + unsigned long int k = 1; + unsigned long int l = 1; + + redcify2 (P1, P0, 1, n1, n0); + addmod2 (x1, x0, P1, P0, P1, P0, n1, n0); /* i.e., redcify(2) */ + y1 = z1 = x1; + y0 = z0 = x0; - if (!mpz_probab_prime_p (g, MR_REPS)) + while (n1 != 0 || n0 != 1) + { + binv (ni, n0); + + for (;;) { do { - mp_limb_t a_limb; - mpn_random (&a_limb, (mp_size_t) 1); - a_int = (int) a_limb; + x0 = mulredc2 (&r1m, x1, x0, x1, x0, n1, n0, ni); + x1 = r1m; + addmod2 (x1, x0, x1, x0, 0, a, n1, n0); + + submod2 (t1, t0, z1, z0, x1, x0, n1, n0); + P0 = mulredc2 (&r1m, P1, P0, t1, t0, n1, n0, ni); + P1 = r1m; + + if (k % 32 == 1) + { + g0 = gcd2_odd (&g1, P1, P0, n1, n0); + if (g1 != 0 || g0 != 1) + goto factor_found; + y1 = x1; y0 = x0; + } } - while (a_int == -2 || a_int == 0); + while (--k != 0); - debug ("[composite factor--restarting pollard-rho] "); - factor_using_pollard_rho (g, a_int); + z1 = x1; z0 = x0; + k = l; + l = 2 * l; + for (unsigned long int i = 0; i < k; i++) + { + x0 = mulredc2 (&r1m, x1, x0, x1, x0, n1, n0, ni); + x1 = r1m; + addmod2 (x1, x0, x1, x0, 0, a, n1, n0); + } + y1 = x1; y0 = x0; + } + + factor_found: + do + { + y0 = mulredc2 (&r1m, y1, y0, y1, y0, n1, n0, ni); + y1 = r1m; + addmod2 (y1, y0, y1, y0, 0, a, n1, n0); + + submod2 (t1, t0, z1, z0, y1, y0, n1, n0); + g0 = gcd2_odd (&g1, t1, t0, n1, n0); + } + while (g1 == 0 && g0 == 1); + + if (g1 == 0) + { + /* The found factor is one word. */ + divexact_21 (n1, n0, n1, n0, g0); /* n = n / g */ + + if (!prime_p (g0)) + factor_using_pollard_rho (g0, a + 1, factors); + else + factor_insert (factors, g0); } else { - emit_factor (g); + /* The found factor is two words. This is highly unlikely, thus hard + to trigger. Please be careful before you change this code! */ + uintmax_t ginv; + + binv (ginv, g0); /* Compute n = n / g. Since the result will */ + n0 = ginv * n0; /* fit one word, we can compute the quotient */ + n1 = 0; /* modulo B, ignoring the high divisor word. */ + + if (!prime2_p (g1, g0)) + factor_using_pollard_rho2 (g1, g0, a + 1, factors); + else + factor_insert_large (factors, g1, g0); + } + + if (n1 == 0) + { + if (prime_p (n0)) + { + factor_insert (factors, n0); + break; + } + + factor_using_pollard_rho (n0, a, factors); + return; } - mpz_mod (x, x, n); - mpz_mod (x1, x1, n); - mpz_mod (y, y, n); - if (mpz_probab_prime_p (n, MR_REPS)) + + if (prime2_p (n1, n0)) { - emit_factor (n); + factor_insert_large (factors, n1, n0); break; } - } - mpz_clear (g); - mpz_clear (P); - mpz_clear (t2); - mpz_clear (t1); - mpz_clear (a); - mpz_clear (x1); - mpz_clear (x); - mpz_clear (y); + x0 = mod2 (&x1, x1, x0, n1, n0); + z0 = mod2 (&z1, z1, z0, n1, n0); + y0 = mod2 (&y1, y1, y0, n1, n0); + } } -#else - +#if HAVE_GMP static void -debug (char const *fmt ATTRIBUTE_UNUSED, ...) +mp_factor_using_pollard_rho (mpz_t n, unsigned long int a, + struct mp_factors *factors) { -} + mpz_t x, z, y, P; + mpz_t t, t2; -#endif + debug ("[pollard-rho (%lu)] ", a); -/* The maximum number of factors, including -1, for negative numbers. */ -#define MAX_N_FACTORS (sizeof (uintmax_t) * CHAR_BIT) + mpz_inits (t, t2, NULL); + mpz_init_set_si (y, 2); + mpz_init_set_si (x, 2); + mpz_init_set_si (z, 2); + mpz_init_set_ui (P, 1); -/* The trial divisor increment wheel. Use it to skip over divisors that - are composites of 2, 3, 5, 7, or 11. The part from WHEEL_START up to - WHEEL_END is reused periodically, while the "lead in" is used to test - for those primes and to jump onto the wheel. For more information, see - http://www.utm.edu/research/primes/glossary/WheelFactorization.html */ + unsigned long long int k = 1; + unsigned long long int l = 1; -#include "wheel-size.h" /* For the definition of WHEEL_SIZE. */ -static const unsigned char wheel_tab[] = - { -#include "wheel.h" - }; + while (mpz_cmp_ui (n, 1) != 0) + { + for (;;) + { + do + { + mpz_mul (t, x, x); + mpz_mod (x, t, n); + mpz_add_ui (x, x, a); + + mpz_sub (t, z, x); + mpz_mul (t2, P, t); + mpz_mod (P, t2, n); + + if (k % 32 == 1) + { + mpz_gcd (t, P, n); + if (mpz_cmp_ui (t, 1) != 0) + goto factor_found; + mpz_set (y, x); + } + } + while (--k != 0); -#define WHEEL_START (wheel_tab + WHEEL_SIZE) -#define WHEEL_END (wheel_tab + ARRAY_CARDINALITY (wheel_tab)) + mpz_set (z, x); + k = l; + l = 2 * l; + for (unsigned long long int i = 0; i < k; i++) + { + mpz_mul (t, x, x); + mpz_mod (x, t, n); + mpz_add_ui (x, x, a); + } + mpz_set (y, x); + } -/* FIXME: comment */ + factor_found: + do + { + mpz_mul (t, y, y); + mpz_mod (y, t, n); + mpz_add_ui (y, y, a); -static size_t -factor_wheel (uintmax_t n0, size_t max_n_factors, uintmax_t *factors) -{ - uintmax_t n = n0, d, q; - size_t n_factors = 0; - unsigned char const *w = wheel_tab; + mpz_sub (t, z, y); + mpz_gcd (t, t, n); + } + while (mpz_cmp_ui (t, 1) == 0); - if (n <= 1) - return n_factors; + mpz_divexact (n, n, t); /* divide by t, before t is overwritten */ - /* The exit condition in the following loop is correct because - any time it is tested one of these 3 conditions holds: - (1) d divides n - (2) n is prime - (3) n is composite but has no factors less than d. - If (1) or (2) obviously the right thing happens. - If (3), then since n is composite it is >= d^2. */ + if (!mp_prime_p (t)) + { + debug ("[composite factor--restarting pollard-rho] "); + mp_factor_using_pollard_rho (t, a + 1, factors); + } + else + { + mp_factor_insert (factors, t); + } - d = 2; - do - { - q = n / d; - while (n == q * d) + if (mp_prime_p (n)) { - assert (n_factors < max_n_factors); - factors[n_factors++] = d; - n = q; - q = n / d; + mp_factor_insert (factors, n); + break; } - d += *(w++); - if (w == WHEEL_END) - w = WHEEL_START; + + mpz_mod (x, x, n); + mpz_mod (z, z, n); + mpz_mod (y, y, n); } - while (d <= q); - if (n != 1 || n0 == 1) + mpz_clears (P, t2, t, z, x, y, NULL); +} +#endif + +/* FIXME: Maybe better to use an iteration converging to 1/sqrt(n)? If + algorithm is replaced, consider also returning the remainder. */ +static uintmax_t _GL_ATTRIBUTE_CONST +isqrt (uintmax_t n) +{ + uintmax_t x; + unsigned c; + if (n == 0) + return 0; + + count_leading_zeros (c, n); + + /* Make x > sqrt(n). This will be invariant through the loop. */ + x = (uintmax_t) 1 << ((W_TYPE_SIZE + 1 - c) / 2); + + for (;;) { - assert (n_factors < max_n_factors); - factors[n_factors++] = n; - } + uintmax_t y = (x + n/x) / 2; + if (y >= x) + return x; - return n_factors; + x = y; + } } -/* Single-precision factoring */ -static void -print_factors_single (uintmax_t n) +static uintmax_t _GL_ATTRIBUTE_CONST +isqrt2 (uintmax_t nh, uintmax_t nl) { - uintmax_t factors[MAX_N_FACTORS]; - size_t n_factors = factor_wheel (n, MAX_N_FACTORS, factors); - size_t i; - char buf[INT_BUFSIZE_BOUND (uintmax_t)]; + unsigned int shift; + uintmax_t x; - printf ("%s:", umaxtostr (n, buf)); - for (i = 0; i < n_factors; i++) - printf (" %s", umaxtostr (factors[i], buf)); - putchar ('\n'); + /* Ensures the remainder fits in an uintmax_t. */ + assert (nh < ((uintmax_t) 1 << (W_TYPE_SIZE - 2))); + + if (nh == 0) + return isqrt (nl); + + count_leading_zeros (shift, nh); + shift &= ~1; + + /* Make x > sqrt(n) */ + x = isqrt ( (nh << shift) + (nl >> (W_TYPE_SIZE - shift))) + 1; + x <<= (W_TYPE_SIZE - shift) / 2; + + /* Do we need more than one iteration? */ + for (;;) + { + uintmax_t r ATTRIBUTE_UNUSED; + uintmax_t q, y; + udiv_qrnnd (q, r, nh, nl, x); + y = (x + q) / 2; + + if (y >= x) + { + uintmax_t hi, lo; + umul_ppmm (hi, lo, x + 1, x + 1); + assert (gt2 (hi, lo, nh, nl)); + + umul_ppmm (hi, lo, x, x); + assert (ge2 (nh, nl, hi, lo)); + sub_ddmmss (hi, lo, nh, nl, hi, lo); + assert (hi == 0); + + return x; + } + + x = y; + } } -#if HAVE_GMP -static int -mpcompare (const void *av, const void *bv) +/* MAGIC[N] has a bit i set iff i is a quadratic residue mod N. */ +#define MAGIC64 ((uint64_t) 0x0202021202030213ULL) +#define MAGIC63 ((uint64_t) 0x0402483012450293ULL) +#define MAGIC65 ((uint64_t) 0x218a019866014613ULL) +#define MAGIC11 0x23b + +/* Returns the square root if the input is a square, otherwise 0. */ +static uintmax_t _GL_ATTRIBUTE_CONST +is_square (uintmax_t x) { - mpz_t *const *a = av; - mpz_t *const *b = bv; - return mpz_cmp (**a, **b); + /* Uses the tests suggested by Cohen. Excludes 99% of the non-squares before + computing the square root. */ + if (((MAGIC64 >> (x & 63)) & 1) + && ((MAGIC63 >> (x % 63)) & 1) + /* Both 0 and 64 are squares mod (65) */ + && ((MAGIC65 >> ((x % 65) & 63)) & 1) + && ((MAGIC11 >> (x % 11) & 1))) + { + uintmax_t r = isqrt (x); + if (r*r == x) + return r; + } + return 0; } -static void -sort_and_print_factors (void) +/* invtab[i] = floor(0x10000 / (0x100 + i) */ +static const unsigned short invtab[0x81] = + { + 0x200, + 0x1fc, 0x1f8, 0x1f4, 0x1f0, 0x1ec, 0x1e9, 0x1e5, 0x1e1, + 0x1de, 0x1da, 0x1d7, 0x1d4, 0x1d0, 0x1cd, 0x1ca, 0x1c7, + 0x1c3, 0x1c0, 0x1bd, 0x1ba, 0x1b7, 0x1b4, 0x1b2, 0x1af, + 0x1ac, 0x1a9, 0x1a6, 0x1a4, 0x1a1, 0x19e, 0x19c, 0x199, + 0x197, 0x194, 0x192, 0x18f, 0x18d, 0x18a, 0x188, 0x186, + 0x183, 0x181, 0x17f, 0x17d, 0x17a, 0x178, 0x176, 0x174, + 0x172, 0x170, 0x16e, 0x16c, 0x16a, 0x168, 0x166, 0x164, + 0x162, 0x160, 0x15e, 0x15c, 0x15a, 0x158, 0x157, 0x155, + 0x153, 0x151, 0x150, 0x14e, 0x14c, 0x14a, 0x149, 0x147, + 0x146, 0x144, 0x142, 0x141, 0x13f, 0x13e, 0x13c, 0x13b, + 0x139, 0x138, 0x136, 0x135, 0x133, 0x132, 0x130, 0x12f, + 0x12e, 0x12c, 0x12b, 0x129, 0x128, 0x127, 0x125, 0x124, + 0x123, 0x121, 0x120, 0x11f, 0x11e, 0x11c, 0x11b, 0x11a, + 0x119, 0x118, 0x116, 0x115, 0x114, 0x113, 0x112, 0x111, + 0x10f, 0x10e, 0x10d, 0x10c, 0x10b, 0x10a, 0x109, 0x108, + 0x107, 0x106, 0x105, 0x104, 0x103, 0x102, 0x101, 0x100, + }; + +/* Compute q = [u/d], r = u mod d. Avoids slow hardware division for the case + that q < 0x40; here it instead uses a table of (Euclidian) inverses. */ +#define div_smallq(q, r, u, d) \ + do { \ + if ((u) / 0x40 < (d)) \ + { \ + int _cnt; \ + uintmax_t _dinv, _mask, _q, _r; \ + count_leading_zeros (_cnt, (d)); \ + _r = (u); \ + if (UNLIKELY (_cnt > (W_TYPE_SIZE - 8))) \ + { \ + _dinv = invtab[((d) << (_cnt + 8 - W_TYPE_SIZE)) - 0x80]; \ + _q = _dinv * _r >> (8 + W_TYPE_SIZE - _cnt); \ + } \ + else \ + { \ + _dinv = invtab[((d) >> (W_TYPE_SIZE - 8 - _cnt)) - 0x7f]; \ + _q = _dinv * (_r >> (W_TYPE_SIZE - 3 - _cnt)) >> 11; \ + } \ + _r -= _q*(d); \ + \ + _mask = -(uintmax_t) (_r >= (d)); \ + (r) = _r - (_mask & (d)); \ + (q) = _q - _mask; \ + assert ( (q) * (d) + (r) == u); \ + } \ + else \ + { \ + uintmax_t _q = (u) / (d); \ + (r) = (u) - _q * (d); \ + (q) = _q; \ + } \ + } while (0) + +/* Notes: Example N = 22117019. After first phase we find Q1 = 6314, Q + = 3025, P = 1737, representing F_{18} = (-6314, 2* 1737, 3025), + with 3025 = 55^2. + + Constructing the square root, we get Q1 = 55, Q = 8653, P = 4652, + representing G_0 = (-55, 2*4652, 8653). + + In the notation of the paper: + + S_{-1} = 55, S_0 = 8653, R_0 = 4652 + + Put + + t_0 = floor([q_0 + R_0] / S0) = 1 + R_1 = t_0 * S_0 - R_0 = 4001 + S_1 = S_{-1} +t_0 (R_0 - R_1) = 706 +*/ + +/* Multipliers, in order of efficiency: + 0.7268 3*5*7*11 = 1155 = 3 (mod 4) + 0.7317 3*5*7 = 105 = 1 + 0.7820 3*5*11 = 165 = 1 + 0.7872 3*5 = 15 = 3 + 0.8101 3*7*11 = 231 = 3 + 0.8155 3*7 = 21 = 1 + 0.8284 5*7*11 = 385 = 1 + 0.8339 5*7 = 35 = 3 + 0.8716 3*11 = 33 = 1 + 0.8774 3 = 3 = 3 + 0.8913 5*11 = 55 = 3 + 0.8972 5 = 5 = 1 + 0.9233 7*11 = 77 = 1 + 0.9295 7 = 7 = 3 + 0.9934 11 = 11 = 3 +*/ +#define QUEUE_SIZE 50 + +#if STAT_SQUFOF +# define Q_FREQ_SIZE 50 +/* Element 0 keeps the total */ +static unsigned int q_freq[Q_FREQ_SIZE + 1]; +# define MIN(a,b) ((a) < (b) ? (a) : (b)) +#endif + +/* Return true on success. Expected to fail only for numbers + >= 2^{2*W_TYPE_SIZE - 2}, or close to that limit. */ +static bool +factor_using_squfof (uintmax_t n1, uintmax_t n0, struct factors *factors) { - mpz_t **faclist; - size_t i; + /* Uses algorithm and notation from + + SQUARE FORM FACTORIZATION + JASON E. GOWER AND SAMUEL S. WAGSTAFF, JR. + + http://homes.cerias.purdue.edu/~ssw/squfof.pdf + */ + + static const unsigned int multipliers_1[] = + { /* = 1 (mod 4) */ + 105, 165, 21, 385, 33, 5, 77, 1, 0 + }; + static const unsigned int multipliers_3[] = + { /* = 3 (mod 4) */ + 1155, 15, 231, 35, 3, 55, 7, 11, 0 + }; + + const unsigned int *m; + + struct { uintmax_t Q; uintmax_t P; } queue[QUEUE_SIZE]; + + if (n1 >= ((uintmax_t) 1 << (W_TYPE_SIZE - 2))) + return false; - faclist = xcalloc (nfactors_found, sizeof *faclist); - for (i = 0; i < nfactors_found; ++i) + uintmax_t sqrt_n = isqrt2 (n1, n0); + + if (n0 == sqrt_n * sqrt_n) { - faclist[i] = &factor[i]; + uintmax_t p1, p0; + + umul_ppmm (p1, p0, sqrt_n, sqrt_n); + assert (p0 == n0); + + if (n1 == p1) + { + if (prime_p (sqrt_n)) + factor_insert_multiplicity (factors, sqrt_n, 2); + else + { + struct factors f; + + f.nfactors = 0; + if (!factor_using_squfof (0, sqrt_n, &f)) + { + /* Try pollard rho instead */ + factor_using_pollard_rho (sqrt_n, 1, &f); + } + /* Duplicate the new factors */ + for (unsigned int i = 0; i < f.nfactors; i++) + factor_insert_multiplicity (factors, f.p[i], 2*f.e[i]); + } + return true; + } } - qsort (faclist, nfactors_found, sizeof *faclist, mpcompare); - for (i = 0; i < nfactors_found; ++i) + /* Select multipliers so we always get n * mu = 3 (mod 4) */ + for (m = (n0 % 4 == 1) ? multipliers_3 : multipliers_1; + *m; m++) { - fputc (' ', stdout); - mpz_out_str (stdout, 10, *faclist[i]); + uintmax_t S, Dh, Dl, Q1, Q, P, L, L1, B; + unsigned int i; + unsigned int mu = *m; + unsigned int qpos = 0; + + assert (mu * n0 % 4 == 3); + + /* In the notation of the paper, with mu * n == 3 (mod 4), we + get \Delta = 4 mu * n, and the paper's \mu is 2 mu. As far as + I understand it, the necessary bound is 4 \mu^3 < n, or 32 + mu^3 < n. + + However, this seems insufficient: With n = 37243139 and mu = + 105, we get a trivial factor, from the square 38809 = 197^2, + without any corresponding Q earlier in the iteration. + + Requiring 64 mu^3 < n seems sufficient. */ + if (n1 == 0) + { + if ((uintmax_t) mu*mu*mu >= n0 / 64) + continue; + } + else + { + if (n1 > ((uintmax_t) 1 << (W_TYPE_SIZE - 2)) / mu) + continue; + } + umul_ppmm (Dh, Dl, n0, mu); + Dh += n1 * mu; + + assert (Dl % 4 != 1); + assert (Dh < (uintmax_t) 1 << (W_TYPE_SIZE - 2)); + + S = isqrt2 (Dh, Dl); + + Q1 = 1; + P = S; + + /* Square root remainder fits in one word, so ignore high part. */ + Q = Dl - P*P; + /* FIXME: When can this differ from floor(sqrt(2 sqrt(D)))? */ + L = isqrt (2*S); + B = 2*L; + L1 = mu * 2 * L; + + /* The form is (+/- Q1, 2P, -/+ Q), of discriminant 4 (P^2 + Q Q1) = + 4 D. */ + + for (i = 0; i <= B; i++) + { + uintmax_t q, P1, t, rem; + + div_smallq (q, rem, S+P, Q); + P1 = S - rem; /* P1 = q*Q - P */ + +#if STAT_SQUFOF + assert (q > 0); + q_freq[0]++; + q_freq[MIN (q, Q_FREQ_SIZE)]++; +#endif + + if (Q <= L1) + { + uintmax_t g = Q; + + if ( (Q & 1) == 0) + g /= 2; + + g /= gcd_odd (g, mu); + + if (g <= L) + { + if (qpos >= QUEUE_SIZE) + error (EXIT_FAILURE, 0, _("squfof queue overflow")); + queue[qpos].Q = g; + queue[qpos].P = P % g; + qpos++; + } + } + + /* I think the difference can be either sign, but mod + 2^W_TYPE_SIZE arithmetic should be fine. */ + t = Q1 + q * (P - P1); + Q1 = Q; + Q = t; + P = P1; + + if ( (i & 1) == 0) + { + uintmax_t r = is_square (Q); + if (r) + { + for (unsigned int j = 0; j < qpos; j++) + { + if (queue[j].Q == r) + { + if (r == 1) + /* Traversed entire cycle. */ + goto next_multiplier; + + /* Need the absolute value for divisibility test. */ + if (P >= queue[j].P) + t = P - queue[j].P; + else + t = queue[j].P - P; + if (t % r == 0) + { + /* Delete entries up to and including entry + j, which matched. */ + memmove (queue, queue + j + 1, + (qpos - j - 1) * sizeof (queue[0])); + qpos -= (j + 1); + } + goto next_i; + } + } + + /* We have found a square form, which should give a + factor. */ + Q1 = r; + assert (S >= P); /* What signs are possible? */ + P += r * ((S - P) / r); + + /* Note: Paper says (N - P*P) / Q1, that seems incorrect + for the case D = 2N. */ + /* Compute Q = (D - P*P) / Q1, but we need double + precision. */ + uintmax_t hi, lo; + umul_ppmm (hi, lo, P, P); + sub_ddmmss (hi, lo, Dh, Dl, hi, lo); + udiv_qrnnd (Q, rem, hi, lo, Q1); + assert (rem == 0); + + for (;;) + { + /* Note: There appears to by a typo in the paper, + Step 4a in the algorithm description says q <-- + floor([S+P]/\hat Q), but looking at the equations + in Sec. 3.1, it should be q <-- floor([S+P] / Q). + (In this code, \hat Q is Q1). */ + div_smallq (q, rem, S+P, Q); + P1 = S - rem; /* P1 = q*Q - P */ + +#if STAT_SQUFOF + q_freq[0]++; + q_freq[MIN (q, Q_FREQ_SIZE)]++; +#endif + if (P == P1) + break; + t = Q1 + q * (P - P1); + Q1 = Q; + Q = t; + P = P1; + } + + if ( (Q & 1) == 0) + Q /= 2; + Q /= gcd_odd (Q, mu); + + assert (Q > 1 && (n1 || Q < n0)); + + if (prime_p (Q)) + factor_insert (factors, Q); + else if (!factor_using_squfof (0, Q, factors)) + factor_using_pollard_rho (Q, 2, factors); + + divexact_21 (n1, n0, n1, n0, Q); + + if (prime2_p (n1, n0)) + factor_insert_large (factors, n1, n0); + else + { + if (!factor_using_squfof (n1, n0, factors)) + { + if (n1 == 0) + factor_using_pollard_rho (n0, 1, factors); + else + factor_using_pollard_rho2 (n1, n0, 1, factors); + } + } + + return true; + } + } + next_i:; + } + next_multiplier:; } - putchar ('\n'); - free (faclist); + return false; } static void -free_factors (void) +factor (uintmax_t t1, uintmax_t t0, struct factors *factors) { - size_t i; + factors->nfactors = 0; + factors->plarge[1] = 0; + + if (t1 == 0 && t0 < 2) + return; + + t0 = factor_using_division (&t1, t1, t0, factors); + + if (t1 == 0 && t0 < 2) + return; - for (i = 0; i < nfactors_found; ++i) + if (prime2_p (t1, t0)) + factor_insert_large (factors, t1, t0); + else { - mpz_clear (factor[i]); + if (alg == ALG_SQUFOF) + if (factor_using_squfof (t1, t0, factors)) + return; + + if (t1 == 0) + factor_using_pollard_rho (t0, 1, factors); + else + factor_using_pollard_rho2 (t1, t0, 1, factors); } - /* Don't actually free factor[] because in the case where we are - reading numbers from stdin, we may be about to use it again. */ - nfactors_found = 0; } -/* Arbitrary-precision factoring */ +#if HAVE_GMP static void -print_factors_multi (mpz_t t) +mp_factor (mpz_t t, struct mp_factors *factors) { - mpz_out_str (stdout, 10, t); - putchar (':'); + mp_factor_init (factors); if (mpz_sgn (t) != 0) { - /* Set the trial division limit according to the size of t. */ - size_t n_bits = mpz_sizeinbase (t, 2); - unsigned int division_limit = MIN (n_bits, 1000); - division_limit *= division_limit; - - factor_using_division (t, division_limit); + mp_factor_using_division (t, factors); if (mpz_cmp_ui (t, 1) != 0) { debug ("[is number prime?] "); - if (mpz_probab_prime_p (t, MR_REPS)) - emit_factor (t); + if (mp_prime_p (t)) + mp_factor_insert (factors, t); else - factor_using_pollard_rho (t, 1); + mp_factor_using_pollard_rho (t, 1, factors); } } - - mpz_clear (t); - sort_and_print_factors (); - free_factors (); } #endif +static strtol_error +strto2uintmax (uintmax_t *hip, uintmax_t *lop, const char *s) +{ + unsigned int lo_carry; + uintmax_t hi, lo; + + strtol_error err; + + hi = lo = 0; + for (;;) + { + unsigned int c = *s++; + if (c == 0) + break; + + if (UNLIKELY (!ISDIGIT (c))) + { + err = LONGINT_INVALID; + break; + } + c -= '0'; + + err = LONGINT_OK; /* we've seen at least one valid digit */ + + if (UNLIKELY (hi > ~(uintmax_t)0 / 10)) + { + err = LONGINT_OVERFLOW; + break; + } + hi = 10 * hi; + + lo_carry = (lo >> (W_TYPE_SIZE - 3)) + (lo >> (W_TYPE_SIZE - 1)); + lo_carry += 10 * lo < 2 * lo; + + lo = 10 * lo; + lo += c; + + lo_carry += lo < c; + hi += lo_carry; + if (UNLIKELY (hi < lo_carry)) + { + err = LONGINT_OVERFLOW; + break; + } + } + + *hip = hi; + *lop = lo; + + return err; +} + +static void +print_uintmaxes (uintmax_t t1, uintmax_t t0) +{ + uintmax_t q, r; + + if (t1 == 0) + printf ("%ju", t0); + else + { + /* Use very plain code here since it seems hard to write fast code + without assuming a specific word size. */ + q = t1 / 1000000000; + r = t1 % 1000000000; + udiv_qrnnd (t0, r, r, t0, 1000000000); + print_uintmaxes (q, t0); + printf ("%09u", (int) r); + } +} + +/* Single-precision factoring */ +static void +print_factors_single (uintmax_t t1, uintmax_t t0) +{ + struct factors factors; + + print_uintmaxes (t1, t0); + putchar (':'); + + factor (t1, t0, &factors); + + for (unsigned int j = 0; j < factors.nfactors; j++) + for (unsigned int k = 0; k < factors.e[j]; k++) + printf (" %ju", factors.p[j]); + + if (factors.plarge[1]) + { + putchar (' '); + print_uintmaxes (factors.plarge[1], factors.plarge[0]); + } + putchar ('\n'); +} /* Emit the factors of the indicated number. If we have the option of using either algorithm, we select on the basis of the length of the number. @@ -434,58 +2288,55 @@ print_factors_multi (mpz_t t) has enough digits, because the algorithm is better. The turnover point depends on the value. */ static bool -print_factors (char const *s) +print_factors (const char *input) { - uintmax_t n; - strtol_error err = xstrtoumax (s, NULL, 10, &n, ""); + uintmax_t t1, t0; -#if HAVE_GMP - enum { GMP_TURNOVER_POINT = 100000 }; + /* Try converting the number to one or two words. If it fails, use GMP or + print an error message. The 2nd condition checks that the most + significant bit of the two-word number is clear, in a typesize neutral + way. */ + strtol_error err = strto2uintmax (&t1, &t0, input); - if (err == LONGINT_OVERFLOW - || (err == LONGINT_OK && GMP_TURNOVER_POINT <= n)) + switch (err) { - mpz_t t; - mpz_init (t); - if (gmp_sscanf (s, "%Zd", t) == 1) + case LONGINT_OK: + if (((t1 << 1) >> 1) == t1) { - debug ("[%s]", _("using arbitrary-precision arithmetic")); - print_factors_multi (t); + debug ("[%s]", _("using single-precision arithmetic")); + print_factors_single (t1, t0); return true; } - err = LONGINT_INVALID; - } -#endif - - switch (err) - { - case LONGINT_OK: - debug ("[%s]", _("using single-precision arithmetic")); - print_factors_single (n); - return true; - - case LONGINT_OVERFLOW: - error (0, 0, _("%s is too large"), quote (s)); - return false; + break; default: - error (0, 0, _("%s is not a valid positive integer"), quote (s)); + error (0, 0, _("%s is not a valid positive integer"), quote (input)); return false; } -} -enum -{ - VERBOSE_OPTION = CHAR_MAX + 1 -}; +#if HAVE_GMP + debug ("[%s]", _("using arbitrary-precision arithmetic")); + mpz_t t; + struct mp_factors factors; -static struct option const long_options[] = -{ - {"verbose", no_argument, NULL, VERBOSE_OPTION}, - {GETOPT_HELP_OPTION_DECL}, - {GETOPT_VERSION_OPTION_DECL}, - {NULL, 0, NULL, 0} -}; + mpz_init_set_str (t, input, 10); + + gmp_printf ("%Zd:", t); + mp_factor (t, &factors); + + for (unsigned int j = 0; j < factors.nfactors; j++) + for (unsigned int k = 0; k < factors.e[j]; k++) + gmp_printf (" %Zd", factors.p[j]); + + mp_factor_clear (&factors); + mpz_clear (t); + putchar ('\n'); + return true; +#else + error (0, 0, _("%s is too large"), quote (input)); + return false; +#endif +} void usage (int status) @@ -535,9 +2386,6 @@ do_stdin (void) int main (int argc, char **argv) { - bool ok; - int c; - initialize_main (&argc, &argv); set_program_name (argv[0]); setlocale (LC_ALL, ""); @@ -546,12 +2394,23 @@ main (int argc, char **argv) atexit (close_stdout); + alg = ALG_POLLARD_RHO; /* Default to Pollard rho */ + + int c; while ((c = getopt_long (argc, argv, "", long_options, NULL)) != -1) { switch (c) { case VERBOSE_OPTION: - verbose = true; + verbose = 1; + break; + + case 's': + alg = ALG_SQUFOF; + break; + + case 'w': + flag_prove_primality = false; break; case_GETOPT_HELP_CHAR; @@ -563,18 +2422,36 @@ main (int argc, char **argv) } } +#if STAT_SQUFOF + if (alg == ALG_SQUFOF) + memset (q_freq, 0, sizeof (q_freq)); +#endif + + bool ok; if (argc <= optind) ok = do_stdin (); else { - int i; ok = true; - for (i = optind; i < argc; i++) + for (int i = optind; i < argc; i++) if (! print_factors (argv[i])) ok = false; } -#if HAVE_GMP - free (factor); + +#if STAT_SQUFOF + if (alg == ALG_SQUFOF && q_freq[0] > 0) + { + double acc_f; + printf ("q freq. cum. freq.(total: %d)\n", q_freq[0]); + for (unsigned int i = 1, acc_f = 0.0; i <= Q_FREQ_SIZE; i++) + { + double f = (double) q_freq[i] / q_freq[0]; + acc_f += f; + printf ("%s%d %.2f%% %.2f%%\n", i == Q_FREQ_SIZE ? ">=" : "", i, + 100.0 * f, 100.0 * acc_f); + } + } #endif + exit (ok ? EXIT_SUCCESS : EXIT_FAILURE); } diff --git a/src/local.mk b/src/local.mk index 98259fa80..6a01ef1c7 100644 --- a/src/local.mk +++ b/src/local.mk @@ -35,7 +35,10 @@ bin_PROGRAMS = @bin_PROGRAMS@ pkglibexec_PROGRAMS = @pkglibexec_PROGRAMS@ # Needed by the testsuite. -noinst_PROGRAMS = src/setuidgid src/getlimits +noinst_PROGRAMS = \ + src/getlimits \ + src/make-prime-list \ + src/setuidgid noinst_HEADERS = \ src/chown-core.h \ @@ -48,20 +51,18 @@ noinst_HEADERS = \ src/fs-is-local.h \ src/group-list.h \ src/ioblksize.h \ + src/longlong.h \ src/ls.h \ src/operand2sig.h \ src/prog-fprintf.h \ src/remove.h \ src/system.h \ - src/wheel-size.h \ - src/wheel.h \ src/uname.h EXTRA_DIST += \ src/dcgen \ src/dircolors.hin \ src/tac-pipe.c \ - src/wheel-gen.pl \ src/extract-magic \ src/c99-to-c89.diff @@ -134,6 +135,12 @@ src_link_LDADD = $(LDADD) src_ln_LDADD = $(LDADD) src_logname_LDADD = $(LDADD) src_ls_LDADD = $(LDADD) + +# This must *not* depend on anything in lib/, since it is used to generate +# src/primes.h. If it depended on libcoreutils.a, that would pull all lib/*.c +# into BUILT_SOURCES. +src_make_prime_list_LDADD = + src_md5sum_LDADD = $(LDADD) src_mkdir_LDADD = $(LDADD) src_mkfifo_LDADD = $(LDADD) @@ -371,19 +378,11 @@ src/dircolors.h: src/dcgen src/dircolors.hin $(AM_V_at)chmod a-w $@-t $(AM_V_at)mv $@-t $@ -wheel_size = 5 - -BUILT_SOURCES += src/wheel-size.h -src/wheel-size.h: Makefile.am - $(AM_V_GEN)rm -f $@ $@-t - $(AM_V_at)echo '#define WHEEL_SIZE $(wheel_size)' > $@-t - $(AM_V_at)chmod a-w $@-t - $(AM_V_at)mv $@-t $@ - -BUILT_SOURCES += src/wheel.h -src/wheel.h: src/wheel-gen.pl Makefile.am +BUILT_SOURCES += src/primes.h +CLEANFILES += src/primes.h +src/primes.h: src/make-prime-list $(AM_V_GEN)rm -f $@ $@-t - $(AM_V_at)$(srcdir)/src/wheel-gen.pl $(wheel_size) > $@-t + $(AM_V_at)src/make-prime-list 5000 > $@-t $(AM_V_at)chmod a-w $@-t $(AM_V_at)mv $@-t $@ diff --git a/src/longlong.h b/src/longlong.h index 84e32ee46..510f40ef2 100644 --- a/src/longlong.h +++ b/src/longlong.h @@ -1,7 +1,7 @@ /* longlong.h -- definitions for mixed size 32/64 bit arithmetic. Copyright 1991, 1992, 1993, 1994, 1996, 1997, 1999, 2000, 2001, 2002, 2003, -2004, 2005, 2007, 2008, 2009, 2011 Free Software Foundation, Inc. +2004, 2005, 2007, 2008, 2009, 2011, 2012 Free Software Foundation, Inc. This file is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free @@ -2055,6 +2055,10 @@ extern UWtype mpn_udiv_qrnnd_r (UWtype, UWtype, UWtype, UWtype *); __GMP_DECLSPEC UWtype __MPN(udiv_w_sdiv) (UWtype *, UWtype, UWtype, UWtype); #endif +#ifndef __GMP_DECLSPEC +#define __GMP_DECLSPEC /* empty */ +#endif + /* If udiv_qrnnd was not defined for this processor, use __udiv_qrnnd_c. */ #if !defined (udiv_qrnnd) #define UDIV_NEEDS_NORMALIZATION 1 diff --git a/src/wheel-gen.pl b/src/wheel-gen.pl deleted file mode 100755 index 65b60d12d..000000000 --- a/src/wheel-gen.pl +++ /dev/null @@ -1,114 +0,0 @@ -#!/usr/bin/perl -w -# Generate the spokes of a wheel, for wheel factorization. - -# Copyright (C) 2001-2012 Free Software Foundation, Inc. - -# This program is free software: you can redistribute it and/or modify -# it under the terms of the GNU General Public License as published by -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. - -# This program is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. - -# You should have received a copy of the GNU General Public License -# along with this program. If not, see <http://www.gnu.org/licenses/>. - -eval 'exec /usr/bin/perl -S $0 ${1+"$@"}' - if 0; - -use strict; -(my $ME = $0) =~ s|.*/||; - -# A global destructor to close standard output with error checking. -sub END -{ - defined fileno STDOUT - or return; - close STDOUT - and return; - warn "$ME: closing standard output: $!\n"; - $? ||= 1; -} - -sub is_prime ($) -{ - my ($n) = @_; - use integer; - - $n == 2 - and return 1; - - my $d = 2; - my $w = 1; - while (1) - { - my $q = $n / $d; - $n == $q * $d - and return 0; - $d += $w; - $q < $d - and last; - $w = 2; - } - return 1; -} - -{ - @ARGV == 1 - or die "$ME: missing argument\n"; - - my $wheel_size = $ARGV[0]; - - my @primes = (2); - my $product = $primes[0]; - my $n_primes = 1; - for (my $i = 3; ; $i += 2) - { - if (is_prime $i) - { - push @primes, $i; - $product *= $i; - ++$n_primes == $wheel_size - and last; - } - } - - my $ws_m1 = $wheel_size - 1; - print <<EOF; -/* The first $ws_m1 elements correspond to the incremental offsets of the - first $wheel_size primes (@primes). The $wheel_size(th) element is the - difference between that last prime and the next largest integer - that is not a multiple of those primes. The remaining numbers - define the wheel. For more information, see - http://www.utm.edu/research/primes/glossary/WheelFactorization.html. */ -EOF - - my @increments; - my $prev = 2; - for (my $i = 3; ; $i += 2) - { - my $rel_prime = 1; - foreach my $divisor (@primes) - { - $i != $divisor && $i % $divisor == 0 - and $rel_prime = 0; - } - - if ($rel_prime) - { - #warn $i, ' ', $i - $prev, "\n"; - push @increments, $i - $prev; - $prev = $i; - - $product + 1 < $i - and last; - } - } - - print join (",\n", @increments), "\n"; - - exit 0; -} diff --git a/tests/factor/run.sh b/tests/factor/run.sh new file mode 100755 index 000000000..6ff24c370 --- /dev/null +++ b/tests/factor/run.sh @@ -0,0 +1,30 @@ +#!/bin/sh +# Test the factor rewrite. +# Expect to be invoked via a file whose basename matches +# /^(\d+)\-(\d+)\-([\da-f]{40})\.sh$/ +# The test is to run this command +# seq $1 $2 | factor | shasum -c --status <(echo $3 -) +# I.e., to ensure that the factorizations of integers $1..$2 +# match what we expect. + +# Copyright (C) 2012 Free Software Foundation, Inc. + +. "${srcdir=.}/tests/init.sh"; path_prepend_ ./src + +# Don't run these tests by default. +very_expensive_ + +print_ver_ factor seq + +# Remove the ".sh" suffix: +t=${ME_%.sh} + +# Make IFS include "-", so that a simple "set" will separate the args: +IFS=-$IFS +set $t +echo "$3 -" > exp + +f=1 +seq $1 $2 | factor | shasum -c --status exp && f=0 + +Exit $f diff --git a/tests/local.mk b/tests/local.mk index 0b6d576f8..486bf312b 100644 --- a/tests/local.mk +++ b/tests/local.mk @@ -16,9 +16,9 @@ ## along with this program. If not, see <http://www.gnu.org/licenses/>. # Indirections required so that we'll still be able to know the -# complete list of our tests even if the user override TESTS from the -# command line (which he's allowed to do by the test harness API). -TESTS = $(all_tests) +# complete list of our tests even if the user overrides TESTS +# from the command line (as permitted by the test harness API). +TESTS = $(all_tests) $(factor_tests) root_tests = $(all_root_tests) EXTRA_DIST += $(all_tests) @@ -89,14 +89,15 @@ TESTS_ENVIRONMENT = \ VERBOSE = yes EXTRA_DIST += \ + init.cfg \ tests/Coreutils.pm \ tests/CuSkip.pm \ tests/CuTmpdir.pm \ tests/d_type-check \ tests/envvar-check \ + tests/factor/run.sh \ tests/filefrag-extent-compare \ tests/fiemap-capable \ - init.cfg \ tests/init.sh \ tests/lang-default \ tests/no-perl \ @@ -624,6 +625,68 @@ all_tests = \ tests/touch/trailing-slash.sh \ $(all_root_tests) +# prefix of 2^64 +p = 184467440737 +# prefix of 2^96 +q = 79228162514264337593543 + +# Each of these numbers has a Pollard rho factor larger than 2^64, +# and thus exercises some hard-to-reach code in factor.c. +t1 = 170141183460469225450570946617781744489 +t2 = 170141183460469229545748130981302223887 +# Factors of the above: +# t1: 9223372036854775421 18446744073709551709 +# t2: 9223372036854775643 18446744073709551709 + +# Each tests is a triple: lo, hi, sha1 of result. +# The test script, run.sh, runs seq lo hi|factor|sha1sum +# and verifies that the actual and expected checksums are the same. +tf = tests/factor +factor_tests = \ + $(tf)/0-10000000-a451244522b1b662c86cb3cbb55aee3e085a61a0.sh \ + $(tf)/10000000-20000000-c792a2e02f1c8536b5121f624b04039d20187016.sh \ + $(tf)/20000000-30000000-8115e8dff97d1674134ec054598d939a2a5f6113.sh \ + $(tf)/30000000-40000000-fe7b832c8e0ed55035152c0f9ebd59de73224a60.sh \ + $(tf)/40000000-50000000-b8786d66c432e48bc5b342ee3c6752b7f096f206.sh \ + $(tf)/50000000-60000000-a74fe518c5f79873c2b9016745b88b42c8fd3ede.sh \ + $(tf)/60000000-70000000-689bc70d681791e5d1b8ac1316a05d0c4473d6db.sh \ + $(tf)/70000000-80000000-d370808f2ab8c865f64c2ff909c5722db5b7d58d.sh \ + $(tf)/80000000-90000000-7978aa66bf2bdb446398336ea6f02605e9a77581.sh \ + $(tf)/$(t1)-$(t1)-4622287c5f040cdb7b3bbe4d19d29a71ab277827.sh \ + $(tf)/$(t2)-$(t2)-dea308253708b57afad357e8c0d2a111460ef50e.sh \ + $(tf)/$(p)08551616-$(p)08651615-66c57cd58f4fb572df7f088d17e4f4c1d4f01bb1.sh \ + $(tf)/$(p)08651616-$(p)08751615-729228e693b1a568ecc85b199927424c7d16d410.sh \ + $(tf)/$(p)08751616-$(p)08851615-5a0c985017c2d285e4698f836f5a059e0b684563.sh \ + $(tf)/$(p)08851616-$(p)08951615-0482295c514e371c98ce9fd335deed0c9c44a4f4.sh \ + $(tf)/$(p)08951616-$(p)09051615-9c0e1105ac7c45e27e7bbeb5e213f530d2ad1a71.sh \ + $(tf)/$(p)09051616-$(p)09151615-604366d2b1d75371d0679e6a68962d66336cd383.sh \ + $(tf)/$(p)09151616-$(p)09251615-9192d2bdee930135b28d7160e6d395a7027871da.sh \ + $(tf)/$(p)09251616-$(p)09351615-bcf56ae55d20d700690cff4d3327b78f83fc01bf.sh \ + $(tf)/$(p)09351616-$(p)09451615-16b106398749e5f24d278ba7c58229ae43f650ac.sh \ + $(tf)/$(p)09451616-$(p)09551615-ad2c6ed63525f8e7c83c4c416e7715fa1bebc54c.sh \ + $(tf)/$(p)09551616-$(p)09651615-2b6f9c11742d9de045515a6627c27a042c49f8ba.sh \ + $(tf)/$(p)09651616-$(p)09751615-54851acd51c4819beb666e26bc0100dc9adbc310.sh \ + $(tf)/$(p)09751616-$(p)09851615-6939c2a7afd2d81f45f818a159b7c5226f83a50b.sh \ + $(tf)/$(p)09851616-$(p)09951615-0f2c8bc011d2a45e2afa01459391e68873363c6c.sh \ + $(tf)/$(p)09951616-$(p)10051615-630dc2ad72f4c222bad1405e6c5bea590f92a98c.sh \ + $(tf)/$(q)940336-$(q)942335-63cbd6313d78247b04d63bbbac50cb8f8d33ff71.sh \ + $(tf)/$(q)942336-$(q)944335-0d03d63653767173182491b86fa18f8f680bb036.sh \ + $(tf)/$(q)944336-$(q)946335-ca43bd38cd9f97cc5bb63613cb19643578640f0b.sh \ + $(tf)/$(q)946336-$(q)948335-86d59545a0c13567fa96811821ea5cde950611b1.sh \ + $(tf)/$(q)948336-$(q)950335-c3740e702fa9c97e6cf00150860e0b936a141a6b.sh \ + $(tf)/$(q)950336-$(q)952335-551c3c4c4640d86fda311b5c3006dac45505c0ce.sh \ + $(tf)/$(q)952336-$(q)954335-b1b0b00463c2f853d70ef9c4f7a96de5cb614156.sh \ + $(tf)/$(q)954336-$(q)956335-8938a484a9ef6bb16478091d294fcde9f8ecea69.sh \ + $(tf)/$(q)956336-$(q)958335-d1ae6bc712d994f35edf55c785d71ddf31f16535.sh \ + $(tf)/$(q)958336-$(q)960335-2374919a89196e1fce93adfe779cb4664556d4b6.sh \ + $(tf)/$(q)960336-$(q)962335-569e4363e8d9e8830a187d9ab27365eef08abde1.sh + +$(factor_tests): tests/factor/run.sh + $(AM_V_GEN)$(MKDIR_P) $(tf) + $(AM_V_at)ln -f $(srcdir)/tests/factor/run.sh $@ + +CLEANFILES += $(factor_tests) + pr_data = \ tests/pr/0F \ tests/pr/0FF \ |