[libmath-prime-util-perl] 09/55: add valuation function
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:53:39 UTC 2015
This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.41
in repository libmath-prime-util-perl.
commit 508de1ac0aec05fa55546f937b667478ce56c294
Author: Dana Jacobsen <dana at acm.org>
Date: Mon Apr 28 18:12:51 2014 -0700
add valuation function
---
Changes | 6 ++++++
TODO | 2 --
XS.xs | 38 ++++++++++++++++++++++++++------------
lib/Math/Prime/Util.pm | 12 +++++++++++-
lib/Math/Prime/Util/PP.pm | 11 +++++++++++
lib/Math/Prime/Util/PPFE.pm | 8 ++++++++
t/19-moebius.t | 17 ++++++++++++++++-
t/81-bignum.t | 6 ++++++
t/92-release-pod-coverage.t | 2 +-
util.c | 12 ++++++++++++
util.h | 1 +
11 files changed, 98 insertions(+), 17 deletions(-)
diff --git a/Changes b/Changes
index 47db7fa..18ded9b 100644
--- a/Changes
+++ b/Changes
@@ -2,6 +2,12 @@ Revision history for Perl module Math::Prime::Util
0.41 2014-05
+ [ADDED]
+
+ - valuation(n,k) how many times does k divide n?
+
+ [FUNCTIONALITY AND PERFORMANCE]
+
- Small improvement to twin_prime_count_approx and nth_twin_prime_approx.
- Better AKS testing in xt/primality-aks.pl.
diff --git a/TODO b/TODO
index 180aa39..cc53f13 100644
--- a/TODO
+++ b/TODO
@@ -74,5 +74,3 @@
(10**13,10**5) takes 2.5x longer, albeit with 6x less memory.
- lucas_sequence with n = 0 or 1
-
-- valuation
diff --git a/XS.xs b/XS.xs
index 8515dce..e36026b 100644
--- a/XS.xs
+++ b/XS.xs
@@ -790,24 +790,38 @@ znlog(IN SV* sva, IN SV* svg, IN SV* svp)
void
kronecker(IN SV* sva, IN SV* svb)
+ ALIAS:
+ valuation = 1
PREINIT:
int astatus, bstatus, abpositive, abnegative;
PPCODE:
astatus = _validate_int(aTHX_ sva, 2);
bstatus = _validate_int(aTHX_ svb, 2);
- /* Are both a and b positive? */
- abpositive = astatus == 1 && bstatus == 1;
- /* Will both fit in IVs? We should use a bitmask return. */
- abnegative = !abpositive
- && (astatus != 0 && SvIOK(sva) && !SvIsUV(sva))
- && (bstatus != 0 && SvIOK(svb) && !SvIsUV(svb));
- if (abpositive || abnegative) {
- UV a = my_svuv(sva);
- UV b = my_svuv(svb);
- int k = (abpositive) ? kronecker_uu(a,b) : kronecker_ss(a,b);
- RETURN_NPARITY(k);
+ if (ix == 0) {
+ /* Are both a and b positive? */
+ abpositive = astatus == 1 && bstatus == 1;
+ /* Will both fit in IVs? We should use a bitmask return. */
+ abnegative = !abpositive
+ && (astatus != 0 && SvIOK(sva) && !SvIsUV(sva))
+ && (bstatus != 0 && SvIOK(svb) && !SvIsUV(svb));
+ if (abpositive || abnegative) {
+ UV a = my_svuv(sva);
+ UV b = my_svuv(svb);
+ int k = (abpositive) ? kronecker_uu(a,b) : kronecker_ss(a,b);
+ RETURN_NPARITY(k);
+ }
+ } else {
+ if (astatus != 0 && bstatus != 0) {
+ UV n = (astatus == -1) ? (UV)(-(my_sviv(sva))) : my_svuv(sva);
+ UV k = (bstatus == -1) ? (UV)(-(my_sviv(svb))) : my_svuv(svb);
+ XSRETURN_UV( valuation(n, k) );
+ }
+ }
+ switch (ix) {
+ case 0: _vcallsub_with_gmp("kronecker"); break;
+ case 1:
+ default: _vcallsub_with_gmp("valuation"); break;
}
- _vcallsub_with_gmp("kronecker");
return; /* skip implicit PUTBACK */
NV
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index 9b13735..987df4b 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -40,7 +40,7 @@ our @EXPORT_OK =
random_maurer_prime random_maurer_prime_with_cert
random_shawe_taylor_prime random_shawe_taylor_prime_with_cert
primorial pn_primorial consecutive_integer_lcm
- gcd lcm factor factor_exp all_factors divisors
+ gcd lcm factor factor_exp all_factors divisors valuation
moebius mertens euler_phi jordan_totient exp_mangoldt liouville
partitions
chebyshev_theta chebyshev_psi
@@ -1878,6 +1878,16 @@ follow the semantics of Mathematica, Pari, and Perl 6, re:
lcm(0, n) = 0 Any zero in list results in zero return
lcm(n,-m) = lcm(n, m) We use the absolute values
+=head2 valuation
+
+ say "$n is divisible by 2 ", valuation($n,2), " times.";
+
+Given integers C<n> and C<k>, returns the numbers of times C<n> is divisible
+by C<k>. This is a very limited version of the algebraic valuation meaning,
+just applied to integers.
+This corresponds to Pari's C<valuation> function.
+C<0> is returned if C<n> or C<k> is one of the values C<-1>, C<0>, or C<1>.
+
=head2 moebius
say "$n is square free" if moebius($n) != 0;
diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm
index 79e1d2c..0ebbadd 100644
--- a/lib/Math/Prime/Util/PP.pm
+++ b/lib/Math/Prime/Util/PP.pm
@@ -1521,6 +1521,17 @@ sub is_power {
0;
}
+sub valuation {
+ my($n, $k) = @_;
+ return 0 if $n < 2 || $k < 2;
+ my $v = 0;
+ while ( !($n % $k) ) {
+ $n /= $k;
+ $v++;
+ }
+ $v;
+}
+
sub is_pseudoprime {
my($n, $base) = @_;
diff --git a/lib/Math/Prime/Util/PPFE.pm b/lib/Math/Prime/Util/PPFE.pm
index ac5c509..eab46ff 100644
--- a/lib/Math/Prime/Util/PPFE.pm
+++ b/lib/Math/Prime/Util/PPFE.pm
@@ -347,6 +347,14 @@ sub is_power {
_validate_positive_integer($a) if defined $a;
return Math::Prime::Util::PP::is_power($n, $a);
}
+sub valuation {
+ my($n, $k) = @_;
+ $n = -$n if defined $n && $n < 0;
+ $k = -$k if defined $k && $k < 0;
+ _validate_positive_integer($n);
+ _validate_positive_integer($k);
+ return Math::Prime::Util::PP::valuation($n, $k);
+}
#############################################################################
diff --git a/t/19-moebius.t b/t/19-moebius.t
index b1cbd09..f9011bb 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 is_power
+ znprimroot znlog kronecker legendre_phi gcd lcm is_power valuation
/;
my $extra = defined $ENV{EXTENDED_TESTING} && $ENV{EXTENDED_TESTING};
@@ -261,6 +261,14 @@ if ($use64) {
push @kroneckers, [-5694706465843977004,9365273357682496999,-1];
}
+my @valuations = (
+ [-4,2, 2],
+ [0,0, 0],
+ [1,0, 0],
+ [96552,6, 3],
+ [1879048192,2, 28],
+);
+
my @legendre_sums = (
[ 0, 92372, 0],
[ 5, 15, 1],
@@ -390,6 +398,7 @@ plan tests => 0 + 1
+ scalar(@mult_orders)
+ scalar(@znlogs)
+ scalar(@legendre_sums)
+ + scalar(@valuations)
+ scalar(keys %powers)
+ scalar(keys %primroots) + 2
+ scalar(keys %jordan_totients)
@@ -595,6 +604,12 @@ while (my($e, $vals) = each (%powers)) {
ok( @fail == 0, "is_power returns $e for " . join(",", at fail) );
}
+###### valuation
+foreach my $r (@valuations) {
+ my($n, $k, $exp) = @$r;
+ is( valuation($n, $k), $exp, "valuation($n,$k) = $exp" );
+}
+
sub cmp_closeto {
my $got = shift;
my $expect = shift;
diff --git a/t/81-bignum.t b/t/81-bignum.t
index fbbe592..f535514 100644
--- a/t/81-bignum.t
+++ b/t/81-bignum.t
@@ -83,6 +83,7 @@ plan tests => 0
+ 2 # ispower
+ 15 # random primes
+ 2 # miller-rabin random
+ + 1 # valuation
+ 1;
# Using GMP makes these tests run about 2x faster on some machines
@@ -133,6 +134,7 @@ use Math::Prime::Util qw/
random_maurer_prime
miller_rabin_random
verify_prime
+ valuation
/;
# TODO: is_strong_lucas_pseudoprime
# ExponentialIntegral
@@ -330,6 +332,10 @@ is( miller_rabin_random( $randprime, 40 ), "0", "80-bit composite fails Miller-R
###############################################################################
+is( valuation(6**10000-1,5), 5, "valuation(6^10000,5) = 5" );
+
+###############################################################################
+
is( $_, 'this should not change', "Nobody clobbered \$_" );
diff --git a/t/92-release-pod-coverage.t b/t/92-release-pod-coverage.t
index b63ddd6..eecf464 100644
--- a/t/92-release-pod-coverage.t
+++ b/t/92-release-pod-coverage.t
@@ -65,7 +65,7 @@ sub mpu_public_regex {
random_maurer_prime random_maurer_prime_with_cert
random_shawe_taylor_prime random_shawe_taylor_prime_with_cert
primorial pn_primorial consecutive_integer_lcm
- gcd lcm factor factor_exp all_factors divisors
+ gcd lcm factor factor_exp all_factors divisors valuation
moebius mertens euler_phi jordan_totient exp_mangoldt liouville
partitions
chebyshev_theta chebyshev_psi
diff --git a/util.c b/util.c
index 47225b8..b914bad 100644
--- a/util.c
+++ b/util.c
@@ -1257,6 +1257,18 @@ int is_power(UV n, UV a)
return (ret == 1) ? 0 : ret;
}
+UV valuation(UV n, UV k)
+{
+ UV v = 0;
+ UV kpower = k;
+ if (k < 2 || n < 2) return 0;
+ if (k == 2) return ctz(n);
+ while ( !(n % kpower) ) {
+ kpower *= k;
+ v++;
+ }
+ return v;
+}
/* How many times does 2 divide n? */
#define padic2(n) ctz(n)
diff --git a/util.h b/util.h
index af8205c..52529b7 100644
--- a/util.h
+++ b/util.h
@@ -27,6 +27,7 @@ extern UV nth_twin_prime_approx(UV n);
extern int powerof(UV n);
extern int is_power(UV n, UV a);
+extern UV valuation(UV n, UV k);
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