[libmath-prime-util-perl] 32/55: vecsum()
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:53:41 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 e1d98f362f40f52ebef23f2fc924fca7d57ffd35
Author: Dana Jacobsen <dana at 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
More information about the Pkg-perl-cvs-commits
mailing list