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 e63256b6e19c55ab980a76a91892965905647509 Author: Dana Jacobsen <d...@acm.org> Date: Wed Jan 29 17:14:46 2014 -0800 First round: add is_perfect_power, is_perfect_square, is_perfect_cube --- XS.xs | 26 +++++++++++++---- aks.c | 50 -------------------------------- lib/Math/Prime/Util.pm | 22 +++++++------- lib/Math/Prime/Util/PP.pm | 4 +-- util.c | 73 +++++++++++++++++++++++++++++++++++++++++++++++ util.h | 3 ++ 6 files changed, 110 insertions(+), 68 deletions(-) diff --git a/XS.xs b/XS.xs index e17bf86..eba9916 100644 --- a/XS.xs +++ b/XS.xs @@ -13,6 +13,7 @@ #include "cache.h" #include "sieve.h" #define FUNC_gcd_ui 1 +#define FUNC_is_perfect_square 1 #include "util.h" #include "primality.h" #include "factor.h" @@ -541,8 +542,11 @@ is_prime(IN SV* svn, ...) is_extra_strong_lucas_pseudoprime = 5 is_frobenius_underwood_pseudoprime = 6 is_aks_prime = 7 - is_pseudoprime = 8 - is_almost_extra_strong_lucas_pseudoprime = 9 + is_perfect_square = 8 + is_perfect_cube = 9 + is_perfect_power = 10 + is_pseudoprime = 11 + is_almost_extra_strong_lucas_pseudoprime = 12 PREINIT: int status; PPCODE: @@ -561,8 +565,11 @@ is_prime(IN SV* svn, ...) case 5: ret = _XS_is_lucas_pseudoprime(n, 2); break; case 6: ret = _XS_is_frobenius_underwood_pseudoprime(n); break; case 7: ret = _XS_is_aks_prime(n); break; - case 8: ret = _XS_is_pseudoprime(n, (items == 1) ? 2 : a); break; - case 9: + case 8: ret = is_perfect_square(n); break; + case 9: ret = is_perfect_cube(n); break; + case 10: ret = is_perfect_power(n); break; + case 11: ret = _XS_is_pseudoprime(n, (items == 1) ? 2 : a); break; + case 12: default: ret = _XS_is_almost_extra_strong_lucas_pseudoprime (n, (items == 1) ? 1 : a); break; } @@ -578,8 +585,11 @@ is_prime(IN SV* svn, ...) case 5: _vcallsub_with_gmp("is_extra_strong_lucas_pseudoprime"); break; case 6: _vcallsub_with_gmp("is_frobenius_underwood_pseudoprime"); break; case 7: _vcallsub_with_gmp("is_aks_prime"); break; - case 8: _vcallsub_with_gmp("is_pseudoprime"); break; - case 9: + case 8: _vcallsub_with_gmp("is_perfect_square"); break; + case 9: _vcallsub_with_gmp("is_perfect_cube"); break; + case 10:_vcallsub_with_gmp("is_perfect_power"); break; + case 11:_vcallsub_with_gmp("is_pseudoprime"); break; + case 12: default:_vcallsub_with_gmp("is_almost_extra_strong_lucas_pseudoprime"); break; } return; /* skip implicit PUTBACK */ @@ -619,6 +629,10 @@ next_prime(IN SV* svn) } } switch (ix) { + /* + case 0: _vcallsub_with_gmp("next_prime"); break; + case 1: _vcallsub_with_gmp("prev_prime"); break; + */ case 0: _vcallsub("_generic_next_prime"); break; case 1: _vcallsub("_generic_prev_prime"); break; case 2: _vcallsub_with_pp("nth_prime"); break; diff --git a/aks.c b/aks.c index 00e8d29..85f40c7 100644 --- a/aks.c +++ b/aks.c @@ -27,60 +27,10 @@ #include "ptypes.h" #include "aks.h" #define FUNC_isqrt 1 -#define FUNC_log2floor 1 #include "util.h" #include "cache.h" #include "mulmod.h" -/* Bach and Sorenson (1993) would be better */ -static int is_perfect_power(UV n) { - UV b, last; - if ((n <= 3) || (n == UV_MAX)) return 0; - if ((n & (n-1)) == 0) return 1; /* powers of 2 */ -#if (BITS_PER_WORD == 32) || (DBL_DIG > 19) - if (1) { -#elif DBL_DIG == 10 - if (n < UVCONST(10000000000)) { -#elif DBL_DIG == 15 - if (n < UVCONST(1000000000000000)) { -#else - if ( n < (UV) pow(10, DBL_DIG) ) { -#endif - /* Simple floating point method. Fast, but need enough mantissa. */ - b = isqrt(n); - if (b*b == n) return 1; /* perfect square */ - last = log2floor(n-1) + 1; - for (b = 3; b < last; b = next_prime(b)) { - UV root = (UV) (pow(n, 1.0 / (double)b) + 0.5); - if ( ((UV)(pow(root, b)+0.5)) == n) return 1; - } - } else { - /* Dietzfelbinger, algorithm 2.3.5 (without optimized exponential) */ - last = log2floor(n-1) + 1; - for (b = 2; b <= last; b++) { - UV a = 1; - UV c = n; - while (c >= HALF_WORD) c = (1+c)>>1; - while ((c-a) >= 2) { - UV m, maxm, p, i; - m = (a+c) >> 1; - maxm = UV_MAX / m; - p = m; - for (i = 2; i <= b; i++) { - if (p > maxm) { p = n+1; break; } - p *= m; - } - if (p == n) return 1; - if (p < n) - a = m; - else - c = m; - } - } - } - return 0; -} - /* Naive znorder. Works well here because limit will be very small. */ static UV order(UV r, UV n, UV limit) { UV j; diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm index 6e96795..bd59e85 100644 --- a/lib/Math/Prime/Util.pm +++ b/lib/Math/Prime/Util.pm @@ -23,6 +23,7 @@ our @EXPORT_OK = is_almost_extra_strong_lucas_pseudoprime is_frobenius_underwood_pseudoprime is_aks_prime + is_perfect_power miller_rabin miller_rabin_random lucas_sequence primes @@ -2765,20 +2766,21 @@ Here is the right way to do PE problem 69 (under 0.03s): Project Euler, problem 187, stupid brute force solution, ~3 minutes: - use Math::Prime::Util qw/factor/; + use Math::Prime::Util qw/forcomposites factor/; my $nsemis = 0; - do { $nsemis++ if scalar factor($_) == 2; } - for 1 .. int(10**8)-1; + forcomposites { $nsemis++ if scalar factor($_) == 2; } int(10**8)-1; say $nsemis; -Here is the best way for PE187. Under 30 milliseconds from the command line: +Here is one of the best ways for PE187: under 20 milliseconds from the +command line. This is faster than Mathematica until C<10^13>. - use Math::Prime::Util qw/primes prime_count/; - use List::Util qw/sum/; - my $limit = shift || int(10**8); - my @primes = @{primes(int(sqrt($limit)))}; - say sum( map { prime_count(int(($limit-1)/$primes[$_-1])) - $_ + 1 } - 1 .. scalar @primes ); + use Math::Prime::Util qw/forprimes prime_count/; + my $limit = shift || int(10**8)-1; + my ($sum, $pc) = (0, 1); + forprimes { + $sum += prime_count(int($limit/$_)) + 1 - $pc++; + } int(sqrt($limit)); + say $sum; Produce the C<matches> result from L<Math::Factor::XS> without skipping: diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm index e1601df..c29351f 100644 --- a/lib/Math/Prime/Util/PP.pm +++ b/lib/Math/Prime/Util/PP.pm @@ -1427,7 +1427,7 @@ sub _gcd_ui { $x; } -sub _is_perfect_power { +sub is_perfect_power { my $n = shift; return 0 if $n <= 3 || $n != int($n); return 1 if ($n & ($n-1)) == 0; # Power of 2 @@ -2097,7 +2097,7 @@ sub _test_anr { sub is_aks_prime { my $n = shift; - return 0 if $n < 2 || _is_perfect_power($n); + return 0 if $n < 2 || is_perfect_power($n); my($log2n, $limit); if ($n > 2**48) { diff --git a/util.c b/util.c index d2b9363..a31ac02 100644 --- a/util.c +++ b/util.c @@ -66,9 +66,11 @@ #include "ptypes.h" #define FUNC_isqrt 1 +#define FUNC_icbrt 1 #define FUNC_lcm_ui 1 #define FUNC_ctz 1 #define FUNC_log2floor 1 +#define FUNC_is_perfect_square #define FUNC_next_prime_in_sieve 1 #define FUNC_prev_prime_in_sieve 1 #include "util.h" @@ -976,6 +978,77 @@ IV mertens(UV n) { return sum; } +int is_perfect_cube(UV n) +{ + UV r = icbrt(n); + return (r*r*r == n); +} + +/* All 32-bit perfect powers, and all for k=17,19,23,31,37 */ +static const UV perfect_powers[] = + {243,2187,3125,7776,16807,78125,100000,161051,177147,248832,279936,371293, + 537824,759375,823543,1419857,1594323,1889568,2476099,3200000,4084101,5153632, + 6436343,7962624,10000000,11881376,17210368,19487171,20511149,24300000, + 28629151,33554432,35831808,39135393,45435424,48828125,52521875,62748517, + 69343957,79235168,90224199,102400000,105413504,115856201,129140163,130691232, + 147008443,164916224,170859375,184528125,205962976,229345007,254803968, + 312500000,345025251,362797056,380204032,410338673,418195493,459165024, + 503284375,550731776,601692057,612220032,656356768,714924299,777600000, + 844596301,893871739,916132832,992436543,1160290625,1162261467,1220703125, + 1252332576,1280000000,1350125107,1453933568,1564031349,1680700000,1801088541, + 1804229351,1934917632,1977326743,2073071593, + UVCONST(2219006624), UVCONST(2373046875), UVCONST(2494357888), + UVCONST(2535525376), UVCONST(2706784157), UVCONST(2887174368), + UVCONST(3077056399), UVCONST(3276800000), UVCONST(3404825447), + UVCONST(3707398432), UVCONST(3939040643), UVCONST(4182119424) +#if BITS_PER_WORD == 64 + ,UVCONST( 94143178827), UVCONST( 762939453125), + UVCONST( 16926659444736), UVCONST( 19073486328125), + UVCONST( 68630377364883), UVCONST( 232630513987207), + UVCONST( 609359740010496), UVCONST( 617673396283947), + UVCONST( 11398895185373143), UVCONST( 11920928955078125), + UVCONST( 100000000000000000), UVCONST( 450283905890997363), + UVCONST( 505447028499293771), UVCONST( 789730223053602816), + UVCONST( 2218611106740436992), UVCONST(8650415919381337933), + UVCONST(10000000000000000000) +#endif +}; +#define NPOWERS (sizeof(perfect_powers)/sizeof(perfect_powers[0])) + +int is_perfect_power(UV n) { + if ((n <= 3) || (n == UV_MAX)) return 0; + if ((n & (n-1)) == 0) return 1; /* powers of 2 */ + if (is_perfect_square(n)) return 1; /* perfect square */ + if (is_perfect_cube(n)) return 1; /* perfect cube */ + { + UV lo = 0; + UV hi = NPOWERS-1; + while (lo < hi) { + UV mid = lo + (hi-lo)/2; + if (perfect_powers[mid] < n) lo = mid+1; + else hi = mid; + } + if (n <= UVCONST(4294967295) || perfect_powers[lo] == n) + return (perfect_powers[lo] == n); + } +#if BITS_PER_WORD == 64 + { /* n > 2**32. If n = p^k, then p in (3 .. 7131) and k in (5,7,11,13) */ + int ib, k; + for (ib = 3; ib <= 6; ib++) { /* prime exponents from 5 to 13 */ + UV b = primes_small[ib]; + UV root = (UV) ( powl(n, 1.0L / b ) + 0.01 ); + UV pk = 1; + while (b) { + if (b & 1) pk *= root; + b >>= 1; + if (b) root *= root; + } + if (n == pk) return 1; + } + } +#endif + return 0; +} /* How many times does 2 divide n? */ #define padic2(n) ctz(n) diff --git a/util.h b/util.h index 76b9fb5..61627d7 100644 --- a/util.h +++ b/util.h @@ -21,6 +21,9 @@ extern UV prime_count_lower(UV x); extern UV prime_count_upper(UV x); extern UV prime_count_approx(UV x); +extern int is_perfect_cube(UV n); +extern int is_perfect_power(UV n); + extern signed char* _moebius_range(UV low, UV high); extern UV* _totient_range(UV low, UV high); extern IV mertens(UV n); -- 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 Pkg-perl-cvs-commits@lists.alioth.debian.org http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/pkg-perl-cvs-commits