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 25dc64a8ab401675d451ccb747fc1abf03530520
Author: Dana Jacobsen <d...@acm.org>
Date: Thu Jan 30 14:29:46 2014 -0800
is_power updates
---
Changes | 7 +++
TODO | 6 --
XS.xs | 2 +-
aks.c | 2 +-
lib/Math/Prime/Util/PP.pm | 15 +++--
lib/Math/Prime/Util/PPFE.pm | 7 +++
t/19-moebius.t | 27 ++++++++-
t/81-bignum.t | 4 ++
util.c | 140 +++++++++++++++++++++-----------------------
util.h | 3 +-
10 files changed, 126 insertions(+), 87 deletions(-)
diff --git a/Changes b/Changes
index f3bbbc7..fedc152 100644
--- a/Changes
+++ b/Changes
@@ -1,5 +1,12 @@
Revision history for Perl module Math::Prime::Util
+0.38 2014-02
+
+ [ADDED]
+
+ - is_power Returns max k if n=p^k. See Pari 2.4.x.
+
+
0.37 2014-01-26
[FUNCTIONALITY AND PERFORMANCE]
diff --git a/TODO b/TODO
index 21f8915..904b6ca 100644
--- a/TODO
+++ b/TODO
@@ -74,9 +74,3 @@
- Ensure a fast path for Math::GMP from MPU -> MPU:GMP -> GMP, and back.
- znlog better implementation
-
-- consider following Pari and doing is_power($n) / is_power($n, $k) instead
- of the three: is_perfect_power, is_perfect_square, is_perfect_cube.
- Document power. Tests for power.
- From Sage:
ispower(3089265681159475043336839581081873360674602365963130114355701114591322241990483812812582393906477998611814245513881)
should return 14.
- Pari's ispower had some problems in 2.4
diff --git a/XS.xs b/XS.xs
index 666e28a..1523068 100644
--- a/XS.xs
+++ b/XS.xs
@@ -562,7 +562,7 @@ 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 = (a == 0) ? is_power(n) : !(is_power(n) % a); break;
+ case 8: ret = is_power(n, a); break;
case 9: ret = _XS_is_pseudoprime(n, (items == 1) ? 2 : a); break;
case 10:
default: ret = _XS_is_almost_extra_strong_lucas_pseudoprime
diff --git a/aks.c b/aks.c
index 85f40c7..d5bf5e9 100644
--- a/aks.c
+++ b/aks.c
@@ -194,7 +194,7 @@ int _XS_is_aks_prime(UV n)
if (n == 2)
return 1;
- if (is_perfect_power(n))
+ if (is_power(n, 0))
return 0;
sqrtn = isqrt(n);
diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm
index c29351f..27576e2 100644
--- a/lib/Math/Prime/Util/PP.pm
+++ b/lib/Math/Prime/Util/PP.pm
@@ -1427,16 +1427,21 @@ sub _gcd_ui {
$x;
}
-sub is_perfect_power {
- my $n = shift;
+sub is_power {
+ my ($n, $a) = @_;
return 0 if $n <= 3 || $n != int($n);
- return 1 if ($n & ($n-1)) == 0; # Power of 2
+ return !(is_power($n) % $a) if defined $a && $a != 0;
+ #return 2 if _is_perfect_square($n);
$n = Math::BigInt->new("$n") unless ref($n) eq 'Math::BigInt';
# Perl 5.6.2 chokes on this, so do it via as_bin
# my $log2n = 0; { my $num = $n; $log2n++ while $num >>= 1; }
my $log2n = length($n->as_bin) - 2;
for (my $e = 2; $e <= $log2n; $e = next_prime($e)) {
- return 1 if $n->copy()->broot($e)->bpow($e) == $n;
+ my $root = $n->copy()->broot($e);
+ if ($root->copy->bpow($e) == $n) {
+ my $next = is_power($root);
+ return ($next == 0) ? $e : $e * $next;
+ }
}
0;
}
@@ -2097,7 +2102,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_power($n);
my($log2n, $limit);
if ($n > 2**48) {
diff --git a/lib/Math/Prime/Util/PPFE.pm b/lib/Math/Prime/Util/PPFE.pm
index 5787af5..f33f78c 100644
--- a/lib/Math/Prime/Util/PPFE.pm
+++ b/lib/Math/Prime/Util/PPFE.pm
@@ -314,6 +314,13 @@ sub chebyshev_psi {
return Math::Prime::Util::PP::chebyshev_psi($n);
}
+sub is_power {
+ my($x, $a) = @_;
+ _validate_positive_integer($x);
+ _validate_positive_integer($a) if defined $a;
+ return Math::Prime::Util::PP::is_power($x, $a);
+}
+
#############################################################################
sub forprimes (&$;$) { ## no critic qw(ProhibitSubroutinePrototypes)
diff --git a/t/19-moebius.t b/t/19-moebius.t
index aee8695..eaf8975 100644
--- a/t/19-moebius.t
+++ b/t/19-moebius.t
@@ -6,7 +6,7 @@ use Test::More;
use Math::Prime::Util
qw/moebius mertens euler_phi jordan_totient divisor_sum exp_mangoldt
chebyshev_theta chebyshev_psi carmichael_lambda znorder liouville
- znprimroot znlog kronecker legendre_phi gcd lcm
+ znprimroot znlog kronecker legendre_phi gcd lcm is_power
/;
my $extra = defined $ENV{EXTENDED_TESTING} && $ENV{EXTENDED_TESTING};
@@ -339,6 +339,21 @@ if ($usexs) {
push @znlogs, [ [5678,5,10007], 8620]; # 5678 = 5^8620 mod 10007
}
+my %powers = (
+ 0 => [-2, -1, 0, 1, 2, 3, 5, 6, 7, 10, 11, 12, 13, 14, 15, 17, 18, 19],
+ 2 => [4, 9, 25, 36, 49],
+ 3 => [8, 27, 125, 343, 17576],
+ 4 => [16, 38416],
+ 9 => [19683, 1000000000],
+);
+if ($use64) {
+ push @{$powers{0}}, 9908918038843197151;
+ push @{$powers{2}}, 18446743927680663841;
+ push @{$powers{3}}, 2250923753991375;
+ push @{$powers{4}}, 1150530828529256001;
+ push @{$powers{9}}, 118587876497;
+}
+
# These are slow with XS, and *really* slow with PP.
if (!$usexs) {
%big_mertens = map { $_ => $big_mertens{$_} }
@@ -371,6 +386,7 @@ plan tests => 0 + 1
+ scalar(@mult_orders)
+ scalar(@znlogs)
+ scalar(@legendre_sums)
+ + scalar(keys %powers)
+ scalar(keys %primroots) + 2
+ scalar(keys %jordan_totients)
+ 2 # Dedekind psi calculated two ways
@@ -566,6 +582,15 @@ foreach my $r (@legendre_sums) {
is( legendre_phi($x, $a), $exp, "legendre_phi($x,$a) = $exp" );
}
+###### is_power
+while (my($e, $vals) = each (%powers)) {
+ my @fail;
+ foreach my $val (@$vals) {
+ push @fail, $val unless is_power($val) == $e;
+ }
+ ok( @fail == 0, "is_power returns $e for " . join(",",@fail) );
+}
+
sub cmp_closeto {
my $got = shift;
my $expect = shift;
diff --git a/t/81-bignum.t b/t/81-bignum.t
index 39ae952..0684139 100644
--- a/t/81-bignum.t
+++ b/t/81-bignum.t
@@ -80,6 +80,7 @@ plan tests => 0
+ 2 # liouville
+ 3 # gcd
+ 3 # lcm
+ + 1 # ispower
+ 15 # random primes
+ 2 # miller-rabin random
+ 1;
@@ -112,6 +113,7 @@ use Math::Prime::Util qw/
liouville
gcd
lcm
+ is_power
pn_primorial
ExponentialIntegral
LogarithmicIntegral
@@ -283,6 +285,8 @@ is( lcm(9999999998987,10000000001011),
99999999999979999998975857, "lcm(p1,p2)"
is( lcm(892478777297173184633,892478777297173184633), 892478777297173184633,
"lcm(p1,p1)" );
is( lcm(23498324,32497832409432,328732487324,328973248732,3487234897324),
1124956497899814324967019145509298020838481660295598696, "lcm(a,b,c,d,e)" );
+is(
is_power(3089265681159475043336839581081873360674602365963130114355701114591322241990483812812582393906477998611814245513881),
14, "ispower(150607571^14) == 14" );
+
###############################################################################
my $randprime;
diff --git a/util.c b/util.c
index 9ad1a76..c681dc6 100644
--- a/util.c
+++ b/util.c
@@ -978,85 +978,81 @@ 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]))
-
-/* TODO: This has to be redone to properly return the highest power */
-int is_power(UV n) {
- int next;
- if ((n <= 3) || (n == UV_MAX)) return 0;
+/* There are at least 4 ways to do this, plus hybrids.
+ * 1) use a table. Great for 32-bit, too big for 64-bit.
+ * 2) Use pow() to check. Relatively slow and FP is always dangerous.
+ * 3) factor or trial factor. Slow for 64-bit.
+ * 4) Dietzfelbinger algorithm 2.3.5. Quite slow.
+ * This currently uses a hybrid of 1 and 2.
+ */
+int powerof(UV n) {
+ int ib;
+ const int iblast = (n > UVCONST(4294967295)) ? 6 : 4;
+ if ((n <= 3) || (n == UV_MAX)) return 1;
if ((n & (n-1)) == 0) return ctz(n); /* powers of 2 */
- if (is_perfect_square(n)) {
- next = is_power(isqrt(n));
- return (next == 0) ? 2 : 2*next;
+ if (is_perfect_square(n)) return 2 * powerof(isqrt(n));
+ { UV cb = icbrt(n); if (cb*cb*cb==n) return 3 * powerof(cb); }
+ for (ib = 3; ib <= iblast; ib++) { /* prime exponents from 5 to 7-or-13 */
+ UV k, pk, root, b = primes_small[ib];
+ root = (UV) ( pow(n, 1.0 / b ) + 0.01 );
+ pk = root * root * root * root * root;
+ for (k = 5; k < b; k++)
+ pk *= root;
+ if (n == pk) return b * powerof(root);
}
- if (is_perfect_cube(n)) {
- next = is_power(icbrt(n));
- return (next == 0) ? 3 : 3*next;
- }
- {
- 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 > 177146) {
+ switch (n) { /* Check for powers of 11, 13, 17, 19 within 32 bits */
+ case 177147: case 48828125: case 362797056: case 1977326743: return 11;
+ case 1594323: case 1220703125: return 13;
+ case 129140163: return 17;
+ case 1162261467: return 19;
+ default: break;
}
- 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;
- for (ib = 3; ib <= 6; ib++) { /* prime exponents from 5 to 13 */
- UV k, b = primes_small[ib];
- UV root = (UV) ( powl(n, 1.0L / b ) + 0.01 );
- UV pk = root * root * root * root * root;
- for (k = 5; k < b; k++)
- pk *= root;
- if (n == pk) {
- next = is_power(root);
- return (next == 0) ? b : b*next;
- }
+ if (n > UVCONST(4294967295)) {
+ switch (n) {
+ case UVCONST(762939453125):
+ case UVCONST(16926659444736):
+ case UVCONST(232630513987207):
+ case UVCONST(100000000000000000):
+ case UVCONST(505447028499293771):
+ case UVCONST(2218611106740436992):
+ case UVCONST(8650415919381337933): return 17;
+ case UVCONST(19073486328125):
+ case UVCONST(609359740010496):
+ case UVCONST(11398895185373143):
+ case UVCONST(10000000000000000000): return 19;
+ case UVCONST(94143178827):
+ case UVCONST(11920928955078125):
+ case UVCONST(789730223053602816): return 23;
+ case UVCONST(68630377364883): return 29;
+ case UVCONST(617673396283947): return 31;
+ case UVCONST(450283905890997363): return 37;
+ default: break;
+ }
}
- }
#endif
- return 0;
+ }
+ return 1;
}
+int is_power(UV n, UV a)
+{
+ int ret;
+ if (a > 0) {
+ if (a == 1 || n <= 1) return 1;
+ if ((a % 2) == 0)
+ return is_perfect_square(n) ? is_power(isqrt(n),a>>1) : 0;
+ if ((a % 3) == 0)
+ { UV cb = icbrt(n); return (cb*cb*cb == n) ? is_power(cb, a/3) : 0; }
+ if ((a % 5) == 0)
+ { UV r5 = (UV)(pow(n,0.2) + 0.1);
+ return (r5*r5*r5*r5*r5 == n) ? is_power(r5, a/5) : 0; }
+ }
+ ret = powerof(n);
+ if (a == 0 && ret == 1) ret = 0;
+ return (a == 0) ? ret : !(ret % a);
+}
+
/* How many times does 2 divide n? */
#define padic2(n) ctz(n)
diff --git a/util.h b/util.h
index 1243fa0..dd8684e 100644
--- a/util.h
+++ b/util.h
@@ -21,7 +21,8 @@ extern UV prime_count_lower(UV x);
extern UV prime_count_upper(UV x);
extern UV prime_count_approx(UV x);
-extern int is_power(UV n);
+extern int powerof(UV n);
+extern int is_power(UV n, UV a);
extern signed char* _moebius_range(UV low, UV high);
extern UV* _totient_range(UV low, UV high);
--
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
