This is an automated email from the git hooks/post-receive script. ppm-guest pushed a commit to annotated tag v0.38 in repository libmath-prime-util-perl.
commit 3783549f9d032d48354d74620c8ca4b2a8c2ebf8 Author: Dana Jacobsen <[email protected]> Date: Thu Feb 27 13:34:36 2014 -0800 Custom lgamma so improved AKS works on all platforms --- Changes | 4 ++-- aks.c | 34 +++++++++++++++++++++++++--------- t/22-aks-prime.t | 6 +----- xt/primality-aks.pl | 16 +++++++++++----- 4 files changed, 39 insertions(+), 21 deletions(-) diff --git a/Changes b/Changes index 26238c2..78b1691 100644 --- a/Changes +++ b/Changes @@ -10,8 +10,8 @@ Revision history for Perl module Math::Prime::Util - Factoring powers is much faster. - - Add Bernstein+Voloch improvements to AKS. - Much faster, though of course still much slower than BPSW. + - Add Bernstein+Voloch improvements to AKS. Much faster than the v6 + implementation, though still terribly slow vs. BPSW or other proofs. [OTHER] diff --git a/aks.c b/aks.c index ed754b5..4154a8b 100644 --- a/aks.c +++ b/aks.c @@ -29,13 +29,8 @@ #define SQRTN_SHORTCUT 1 -/* Use improvements from Bornemann's 2002 implementation if we have lgamma */ -#if !defined(_MSC_VER) && \ - (defined(__USE_ISOC99) || defined(_ISOC99_SOURCE) || defined(_NETBSD_SOURCE) || defined(__cplusplus) || !defined(__STRICT_ANSI__) || __STDC_VERSION__ >= 199901L || _POSIX_C_SOURCE >= 200112L || _XOPEN_SOURCE >= 600) +/* Use improvements from Bornemann's 2002 implementation */ #define IMPL_BORNEMANN 1 -#else -#define IMPL_BORNEMANN 0 -#endif #include "ptypes.h" #include "aks.h" @@ -58,6 +53,28 @@ static int is_primitive_root(UV n, UV r) } return 1; } +/* We could use lgamma, but it isn't in MSVC and not in pre-C99. The only + * sure way to find if it is available is test compilation (ala autoconf). + * Instead, we'll just use our own implementation. + * See http://mrob.com/pub/ries/lanczos-gamma.html for alternates. */ +static double lanczos_coef[8+1] = +{ 0.99999999999980993, 676.5203681218851, -1259.1392167224028, + 771.32342877765313, -176.61502916214059, 12.507343278686905, + -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7 }; +static double log_sqrt_two_pi = 0.91893853320467274178; +static double log_gamma(double x) +{ + double base = x + 7 + 0.5; + double sum = 0; + int i; + for (i = 8; i >= 1; i--) + sum += lanczos_coef[i] / (x + (double)i); + sum += lanczos_coef[0]; + sum = log_sqrt_two_pi + logl(sum/x) + ( (x+0.5)*logl(base) - base ); + return sum; +} +#undef lgamma +#define lgamma(x) log_gamma(x) #else /* Naive znorder. Works well here because limit will be very small. */ static UV order(UV r, UV n, UV limit) { @@ -254,8 +271,6 @@ int _XS_is_aks_prime(UV n) return 1; s = (UV) floor(sqrt(r-1) * log2n); - - if (verbose) { printf("# aks r = %lu s = %lu\n", (unsigned long) r, (unsigned long) s); } } #else { @@ -293,9 +308,10 @@ int _XS_is_aks_prime(UV n) if (slim >= HALF_WORD || (slim*slim) >= n) return 1; } - if (verbose) { printf("# aks r = %lu s = %lu\n", (unsigned long) r, (unsigned long) s); } #endif + if (verbose) { printf("# aks r = %lu s = %lu\n", (unsigned long) r, (unsigned long) s); } + for (a = 1; a <= s; a++) { if (! test_anr(a, n, r) ) return 0; diff --git a/t/22-aks-prime.t b/t/22-aks-prime.t index 58089aa..7418599 100644 --- a/t/22-aks-prime.t +++ b/t/22-aks-prime.t @@ -37,11 +37,7 @@ SKIP: { # If we're pure Perl, then this is definitely too slow. # Arguably we should check to see if they have the GMP code. skip "Skipping PP AKS on PP -- just too slow.", 1 if $ispp; - # If we have 64-bit available in the compiler (e.g. uint64_t), this can - # still be quite fast. However for pretty much everyone else, this is - # just far too slow for running in a test suite. - skip "Skipping PP AKS on 32-bit -- just too slow.", 1 if !$use64; - # The first number that makes it past the sqrt test to actually run. + # The least number that performs the full test with either implementation. is( is_aks_prime(69197), 1, "is_aks_prime(69197) is true" ); } diff --git a/xt/primality-aks.pl b/xt/primality-aks.pl index c0cf7f5..4107acc 100755 --- a/xt/primality-aks.pl +++ b/xt/primality-aks.pl @@ -5,16 +5,22 @@ $| = 1; # fast pipes my $limit = shift || 10_000_000; -use Math::Prime::Util qw/is_aks_prime/; -use Math::Prime::FastSieve; -my $sieve = Math::Prime::FastSieve::Sieve->new($limit + 10_000); +use Math::Prime::Util qw/is_aks_prime next_prime/; + +# It doesn't really matter for this use, but try to use independent code for +# next_prime. Math::Prime::FastSieve works well. +my $sieve; +if (eval { require Math::Prime::FastSieve; Math::Prime::FastSieve->import(); require Math::Prime::FastSieve::Sieve; 1; }) { + $sieve = Math::Prime::FastSieve::Sieve->new($limit + 10_000); +} if (1) { my $n = 2; + my $i = 1; while ($n <= $limit) { - print "$n\n" if $n > 69000; # unless $i++ % 262144; + print "$n\n" unless $i++ % 1000; die "$n should be prime" unless is_aks_prime($n); - my $next = $sieve->nearest_ge( $n+1 ); + my $next = (defined $sieve) ? $sieve->nearest_ge( $n+1 ) : next_prime($n); my $diff = ($next - $n) >> 1; if ($diff > 1) { foreach my $d (1 .. $diff-1) { -- Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/pkg-perl/packages/libmath-prime-util-perl.git _______________________________________________ Pkg-perl-cvs-commits mailing list [email protected] http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/pkg-perl-cvs-commits
