[libmath-prime-util-perl] 04/54: First round: add is_perfect_power, is_perfect_square, is_perfect_cube
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:52:06 UTC 2015
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 <dana at 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
More information about the Pkg-perl-cvs-commits
mailing list