[libmath-prime-util-perl] 07/54: Change to is_power, still needs work
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:52:07 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 f174ddec326be9d142b6db4e86d7f1ad4df73598
Author: Dana Jacobsen <dana at acm.org>
Date: Thu Jan 30 09:16:07 2014 -0800
Change to is_power, still needs work
---
TODO | 2 ++
XS.xs | 25 +++++++++----------------
lib/Math/Prime/Util.pm | 19 ++++++++++++++++---
util.c | 32 ++++++++++++++++++++------------
util.h | 3 +--
5 files changed, 48 insertions(+), 33 deletions(-)
diff --git a/TODO b/TODO
index 31f4abc..21f8915 100644
--- a/TODO
+++ b/TODO
@@ -78,3 +78,5 @@
- 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 eba9916..666e28a 100644
--- a/XS.xs
+++ b/XS.xs
@@ -13,7 +13,6 @@
#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"
@@ -542,11 +541,9 @@ is_prime(IN SV* svn, ...)
is_extra_strong_lucas_pseudoprime = 5
is_frobenius_underwood_pseudoprime = 6
is_aks_prime = 7
- is_perfect_square = 8
- is_perfect_cube = 9
- is_perfect_power = 10
- is_pseudoprime = 11
- is_almost_extra_strong_lucas_pseudoprime = 12
+ is_power = 8
+ is_pseudoprime = 9
+ is_almost_extra_strong_lucas_pseudoprime = 10
PREINIT:
int status;
PPCODE:
@@ -565,11 +562,9 @@ 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 = 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:
+ case 8: ret = (a == 0) ? is_power(n) : !(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
(n, (items == 1) ? 1 : a); break;
}
@@ -585,11 +580,9 @@ 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_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:
+ case 8: _vcallsub_with_gmp("is_power"); break;
+ case 9:_vcallsub_with_gmp("is_pseudoprime"); break;
+ case 10:
default:_vcallsub_with_gmp("is_almost_extra_strong_lucas_pseudoprime"); break;
}
return; /* skip implicit PUTBACK */
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index 661daf0..51519b7 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -23,7 +23,7 @@ our @EXPORT_OK =
is_almost_extra_strong_lucas_pseudoprime
is_frobenius_underwood_pseudoprime
is_aks_prime
- is_perfect_power
+ is_power
miller_rabin miller_rabin_random
lucas_sequence
primes
@@ -1761,9 +1761,22 @@ Brent, 2010: "AKS is not a practical algorithm. ECPP is much faster."
We have ECPP, and indeed it is much faster.
-=head2 is_perfect_power
+=head2 is_power
-Given a positive integer input n, returns 1 if C<n = p^k> for C<k E<gt> 1>.
+ say "$n is a perfect square" if is_power($n, 2);
+ say "$n is a perfect cube" if is_power($n, 3);
+ say "$n is a ", is_power($n), "-th power";
+
+Given a single positive integer input C<n>, returns k if C<n = p^k> for
+some integer C<p E<gt> 1, k E<gt> 1>, and 0 otherwise. The k returned is
+the largest possible. This can be used in a boolean statement to
+determine if C<n> is a perfect power.
+
+If given two arguments C<n> and C<k>, returns 1 if C<n> is a C<k-th> power,
+and 0 otherwise.
+
+This corresponds to Pari/GP's C<ispower> function, with the limitations of
+only integer arguments and no third argument may be given returning the root.
=head2 lucas_sequence
diff --git a/util.c b/util.c
index a31ac02..9ad1a76 100644
--- a/util.c
+++ b/util.c
@@ -1015,11 +1015,19 @@ static const UV perfect_powers[] =
};
#define NPOWERS (sizeof(perfect_powers)/sizeof(perfect_powers[0]))
-int is_perfect_power(UV n) {
+/* 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;
- 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 */
+ 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_cube(n)) {
+ next = is_power(icbrt(n));
+ return (next == 0) ? 3 : 3*next;
+ }
{
UV lo = 0;
UV hi = NPOWERS-1;
@@ -1033,17 +1041,17 @@ int is_perfect_power(UV 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;
+ int ib;
for (ib = 3; ib <= 6; ib++) { /* prime exponents from 5 to 13 */
- UV b = primes_small[ib];
+ UV k, 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;
+ 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 == pk) return 1;
}
}
#endif
diff --git a/util.h b/util.h
index 61627d7..1243fa0 100644
--- a/util.h
+++ b/util.h
@@ -21,8 +21,7 @@ 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 int is_power(UV n);
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
More information about the Pkg-perl-cvs-commits
mailing list