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 <d...@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(",",@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 _______________________________________________ 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