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 e1d98f362f40f52ebef23f2fc924fca7d57ffd35 Author: Dana Jacobsen <d...@acm.org> Date: Sun May 11 18:05:19 2014 -0700 vecsum() --- Changes | 1 + XS.xs | 60 +++++++++++++++++++++++++------------ lib/Math/Prime/Util.pm | 25 ++++++++++++---- lib/Math/Prime/Util/PP.pm | 13 ++++++++ lib/Math/Prime/Util/PPFE.pm | 3 ++ lib/Math/Prime/Util/RandomPrimes.pm | 24 ++++++++------- t/19-moebius.t | 21 ++++++++++++- t/24-partitions.t | 9 ++++-- 8 files changed, 117 insertions(+), 39 deletions(-) diff --git a/Changes b/Changes index 5352c1b..d1e2367 100644 --- a/Changes +++ b/Changes @@ -7,6 +7,7 @@ Revision history for Perl module Math::Prime::Util - valuation(n,k) how many times does k divide n? - invmod(a,n) inverse of a modulo n - forpart { ... } n[,{...}] loop over partitions (Pari 2.6.x) + - vecsum(...) sum list of integers [FUNCTIONALITY AND PERFORMANCE] diff --git a/XS.xs b/XS.xs index 21c1f3a..3950a3f 100644 --- a/XS.xs +++ b/XS.xs @@ -496,28 +496,49 @@ gcd(...) PROTOTYPE: @ ALIAS: lcm = 1 + vecsum = 2 PREINIT: int i, status = 1; UV ret, nullv, n; PPCODE: - /* For each arg, while valid input, validate+gcd/lcm. Shortcut stop. */ - if (ix == 0) { ret = 0; nullv = 1; } - else { ret = (items == 0) ? 0 : 1; nullv = 0; } - for (i = 0; i < items && ret != nullv && status != 0; i++) { - status = _validate_int(aTHX_ ST(i), 2); - if (status == 0) - break; - n = status * my_svuv(ST(i)); /* n = abs(arg) */ - if (i == 0) { - ret = n; - } else { - UV gcd = gcd_ui(ret, n); - if (ix == 0) { - ret = gcd; + if (ix == 2) { + UV lo = 0; + IV hi = 0; + for (ret = i = 0; i < items; i++) { + status = _validate_int(aTHX_ ST(i), 2); + if (status == 0) break; + n = my_svuv(ST(i)); + if (status == 1) { + hi += (n > (UV_MAX - lo)); + } else { + if (UV_MAX-n == (UV)IV_MAX) { status = 0; break; } /* IV Overflow */ + hi -= ((UV_MAX-n) >= lo); + } + lo += n; + } + if (status != 0 && hi == -1 && lo > IV_MAX) XSRETURN_IV((IV)lo); + if (hi != 0) status = 0; /* Overflow */ + ret = lo; + } else { + /* For each arg, while valid input, validate+gcd/lcm. Shortcut stop. */ + if (ix == 0) { ret = 0; nullv = 1; } + else { ret = (items == 0) ? 0 : 1; nullv = 0; } + for (i = 0; i < items && ret != nullv && status != 0; i++) { + status = _validate_int(aTHX_ ST(i), 2); + if (status == 0) + break; + n = status * my_svuv(ST(i)); /* n = abs(arg) */ + if (i == 0) { + ret = n; } else { - n /= gcd; - if (n <= (UV_MAX / ret) ) ret *= n; - else status = 0; /* Overflow */ + UV gcd = gcd_ui(ret, n); + if (ix == 0) { + ret = gcd; + } else { + n /= gcd; + if (n <= (UV_MAX / ret) ) ret *= n; + else status = 0; /* Overflow */ + } } } } @@ -525,8 +546,9 @@ gcd(...) XSRETURN_UV(ret); switch (ix) { case 0: _vcallsub_with_gmp("gcd"); break; - case 1: - default:_vcallsub_with_gmp("lcm"); break; + case 1: _vcallsub_with_gmp("lcm"); break; + case 2: + default:_vcallsub_with_gmp("vecsum"); break; } return; /* skip implicit PUTBACK */ diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm index ff51a04..3b98318 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 valuation invmod + gcd lcm factor factor_exp all_factors divisors valuation invmod vecsum moebius mertens euler_phi jordan_totient exp_mangoldt liouville partitions chebyshev_theta chebyshev_psi @@ -1904,6 +1904,19 @@ 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 vecsum + + say "Totient sum 500,000: ", vecsum(euler_phi(0,500_000)); + +Returns the sum of all arguments, each of which must be an integer. This +is similar to List::Util's L<List::Util/sum0> function, but has a very +important difference. List::Util turns all inputs into doubles and returns +a double, which will mean incorrect results with large integers. C<vecsum> +sums (signed) integers and returns the untruncated result. Processing is +done on native integers while possible. + + =head2 invmod say "The inverse of 42 mod 2017 = ", invmod(42,2017); @@ -2274,6 +2287,10 @@ Examples of various ways to set your own irand function: # System rand. You probably don't want to do this. prime_set_config(irand => sub { int(rand(4294967296)) }); + # Math::Random::MTwist. Fastest RNG by quite a bit. + use Math::Random::MTwist; + prime_set_config(irand => \&Math::Random::MTwist::irand32); + # Math::Random::Secure. Uses ISAAC and strong seed methods. use Math::Random::Secure; prime_set_config(irand => \&Math::Random::Secure::irand); @@ -2286,7 +2303,7 @@ Examples of various ways to set your own irand function: } prime_set_config(irand => \&irand); - # Crypt::Random. Uses Pari and /dev/random. Very slow. + # Crypt::Random. Uses Pari and /dev/random. *VERY* slow. use Crypt::Random qw/makerandom/; prime_set_config(irand => sub { makerandom(Size=>32, Uniform=>1); }); @@ -2296,10 +2313,6 @@ Examples of various ways to set your own irand function: sub nr_irand { return $rng->get(1); } } prime_set_config(irand => \&nr_irand); - # Mersenne Twister. Very fast, decent RNG, auto seeding. - use Math::Random::MT::Auto; - prime_set_config(irand=>sub {Math::Random::MT::Auto::irand() & 0xFFFFFFFF}); - # Go back to MPU's default configuration prime_set_config(irand => undef); diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm index 8a133a0..4027c5e 100644 --- a/lib/Math/Prime/Util/PP.pm +++ b/lib/Math/Prime/Util/PP.pm @@ -1484,6 +1484,19 @@ sub lcm { $lcm = _bigint_to_int($lcm) if $lcm->bacmp(''.~0) <= 0; return $lcm; } +sub vecsum { + my $sum = 0; + my $neglim = -(~0 >> 1) - 1; + foreach my $v (@_) { + $sum += $v; + if ($sum > ~0 || $sum < $neglim) { + $sum = BZERO->copy; + $sum->badd("$_") for @_; + return $sum; + } + } + $sum; +} sub invmod { my($a,$n) = @_; return if $n == 0 || $a == 0; diff --git a/lib/Math/Prime/Util/PPFE.pm b/lib/Math/Prime/Util/PPFE.pm index 0346683..2951ad9 100644 --- a/lib/Math/Prime/Util/PPFE.pm +++ b/lib/Math/Prime/Util/PPFE.pm @@ -321,6 +321,9 @@ sub gcd { sub lcm { return Math::Prime::Util::PP::lcm(@_); } +sub vecsum { + return Math::Prime::Util::PP::vecsum(@_); +} sub invmod { my ($a, $n) = @_; my ($va, $vn) = ($a, $n); diff --git a/lib/Math/Prime/Util/RandomPrimes.pm b/lib/Math/Prime/Util/RandomPrimes.pm index 2e71c12..8fe99fc 100644 --- a/lib/Math/Prime/Util/RandomPrimes.pm +++ b/lib/Math/Prime/Util/RandomPrimes.pm @@ -242,19 +242,21 @@ sub get_randf_nbit { # possibility is to, if they do not supply a rand function, use the GMP MT # function with an appropriate seed. # -# Random timings for 10M calls: -# 1.92 system rand -# 2.62 Math::Random::MT::Auto -# 12.0 Math::Random::Secure w/ISAAC::XS -# 12.6 Bytes::Random::Secure OO w/ISAAC::XS <==== our -# 31.1 Bytes::Random::Secure OO <==== default -# 44.5 Bytes::Random::Secure function w/ISAAC::XS -# 44.8 Math::Random::Secure -# 71.5 Bytes::Random::Secure function -# 1840 Crypt::Random +# Random timings for 10M calls on i4770K: +# 0.39 Math::Random::MTwist 0.13 +# 0.89 system rand +# 1.76 Math::Random::MT::Auto +# 5.35 Bytes::Random::Secure OO w/ISAAC::XS <==== our +# 7.43 Math::Random::Secure w/ISAAC::XS +# 12.40 Math::Random::Secure +# 12.78 Bytes::Random::Secure OO <==== default +# 13.86 Bytes::Random::Secure function w/ISAAC::XS +# 21.95 Bytes::Random::Secure function +# 822.1 Crypt::Random # +# time perl -E 'use Math::Random::MTwist "irand32"; irand32() for 1..10000000;' # time perl -E 'sub irand {int(rand(4294967296));} irand() for 1..10000000;' -# time perl -E 'use Math::Random::MT::Auto qw/irand/; irand() for 1..10000000;' +# time perl -E 'use Math::Random::MT::Auto; sub irand { Math::Random::MT::Auto::irand() & 0xFFFFFFFF } irand() for 1..10000000;' # time perl -E 'use Math::Random::Secure qw/irand/; irand() for 1..10000000;' # time perl -E 'use Bytes::Random::Secure qw/random_bytes/; sub irand {return unpack("L",random_bytes(4));} irand() for 1..10000000;' # time perl -E 'use Bytes::Random::Secure; my $rng = Bytes::Random::Secure->new(); sub irand {return $rng->irand;} irand() for 1..10000000;' diff --git a/t/19-moebius.t b/t/19-moebius.t index 85a9432..9fd42f4 100644 --- a/t/19-moebius.t +++ b/t/19-moebius.t @@ -7,7 +7,7 @@ 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 valuation - invmod + invmod vecsum /; my $extra = defined $ENV{EXTENDED_TESTING} && $ENV{EXTENDED_TESTING}; @@ -384,6 +384,19 @@ if ($use64) { push @invmods, [ 14, 18446744073709551615, 17129119497016012214 ]; } +my @vecsums = ( + [ 0 ], + [ -1, -1 ], + [ 0, 1,-1 ], + [ 0, -1,1 ], + [ 0, -1,1 ], + [ 0, -2147483648,2147483648 ], + [ 0, "-4294967296","4294967296" ], + [ 0, "-9223372036854775808","9223372036854775808" ], + [ "18446744073709551615", "18446744073709551615","-18446744073709551615","18446744073709551615" ], + [ "55340232221128654848", "18446744073709551616","18446744073709551616","18446744073709551616" ], +); + # These are slow with XS, and *really* slow with PP. if (!$usexs) { %big_mertens = map { $_ => $big_mertens{$_} } @@ -418,6 +431,7 @@ plan tests => 0 + 1 + scalar(@legendre_sums) + scalar(@valuations) + 3 + scalar(@invmods) + + scalar(@vecsums) + scalar(keys %powers) + scalar(keys %primroots) + 2 + scalar(keys %jordan_totients) @@ -637,6 +651,11 @@ foreach my $r (@invmods) { my($a, $n, $exp) = @$r; is( invmod($a,$n), $exp, "invmod($a,$n) = ".((defined $exp)?$exp:"<undef>") ); } +###### vecsum +foreach my $r (@vecsums) { + my($exp, @vals) = @$r; + is( vecsum(@vals), $exp, "vecsum(@vals) = $exp" ); +} sub cmp_closeto { my $got = shift; diff --git a/t/24-partitions.t b/t/24-partitions.t index 6e2cff4..e117cf9 100644 --- a/t/24-partitions.t +++ b/t/24-partitions.t @@ -3,7 +3,7 @@ use strict; use warnings; use Test::More; -use Math::Prime::Util qw/partitions forpart/; +use Math::Prime::Util qw/partitions forpart is_prime/; my $extra = defined $ENV{EXTENDED_TESTING} && $ENV{EXTENDED_TESTING}; my @parts = qw/ @@ -76,7 +76,7 @@ if (!$extra) { foreach my $n (@ns) { delete $bparts{$n} } } -plan tests => scalar(@parts) + scalar(keys(%bparts)) + 14; +plan tests => scalar(@parts) + scalar(keys(%bparts)) + 15; foreach my $n (0..$#parts) { @@ -146,3 +146,8 @@ while (my($n, $epart) = each (%bparts)) { # is_deeply( [@p1], [@p2], "forpart 22 restricted n=4 and a=[3..6]" ); } { my @p=(); forpart { push @p, [@_] } 22, {amin=>2,amax=>8,n=>4}; is_deeply( [@p], [[8,8,4,2],[8,8,3,3],[8,7,5,2],[8,7,4,3],[8,6,6,2], [8,6,5,3], [8,6,4,4], [8,5,5,4], [7,7,6,2], [7,7,5,3], [7,7,4,4], [7,6,6,3], [7,6,5,4], [7,5,5,5], [6,6,6,4], [6,6,5,5]], "forpart 22 restricted n=4 and a=[3..6]" ); } + +{ my @p = (); + forpart { push @p, [@_] unless scalar grep {!is_prime($_)} @_ } 20,{amin=>3}; + is_deeply( [@p], [[17,3], [13,7], [11,3,3,3], [7,7,3,3], [7,5,5,3], [5,5,5,5], [5,3,3,3,3,3]], "forpart 20 restricted to odd primes" ); +} -- 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