[libmath-prime-util-perl] 21/181: Add fordivisors, forcomposites
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:51:02 UTC 2015
This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.36
in repository libmath-prime-util-perl.
commit 6e901e6534a73594bcd162bc50dea54c2caad2dd
Author: Dana Jacobsen <dana at acm.org>
Date: Wed Dec 18 21:06:09 2013 -0800
Add fordivisors, forcomposites
---
Changes | 1 +
TODO | 4 +-
XS.xs | 135 ++++++++++++++++++++++++++++++++++++++++++++++++-
lib/Math/Prime/Util.pm | 61 ++++++++++++++++++++--
t/19-moebius.t | 30 +++++++----
t/32-iterators.t | 22 +++++++-
6 files changed, 233 insertions(+), 20 deletions(-)
diff --git a/Changes b/Changes
index 0bccb3c..547926b 100644
--- a/Changes
+++ b/Changes
@@ -15,6 +15,7 @@ Revision history for Perl module Math::Prime::Util
[ADDED]
- forcomposites like forprimes, but on composites. See Pari 2.6.x.
+ - fordivisors calls a block for each divisor
[FUNCTIONALITY AND PERFORMANCE]
diff --git a/TODO b/TODO
index 91a2698..a7d6a01 100644
--- a/TODO
+++ b/TODO
@@ -62,6 +62,4 @@
- Add Inverse Li and Legendre Phi to API?
-- Document forcomposites. Consider XS for it.
-
-- Add fordivisors as XS.
+- More efficient forcomposites
diff --git a/XS.xs b/XS.xs
index 0bdcc17..67e2567 100644
--- a/XS.xs
+++ b/XS.xs
@@ -39,7 +39,7 @@
#endif
/* multicall compatibility stuff */
-#if PERL_REVISION <= 5 && PERL_VERSION < 7
+#if (PERL_REVISION <= 5 && PERL_VERSION < 7) || !defined(dMULTICALL)
# define USE_MULTICALL 0 /* Too much trouble to work around it */
#else
# define USE_MULTICALL 1
@@ -47,7 +47,7 @@
#if PERL_VERSION < 13 || (PERL_VERSION == 13 && PERL_SUBVERSION < 9)
# define FIX_MULTICALL_REFCOUNT \
- if (CvDEPTH(multicall_cv) > 1) SvREFCNT_inc(multicall_cv);
+ if (CvDEPTH(multicall_cv) > 1) SvREFCNT_inc_simple_void_NN(multicall_cv);
#else
# define FIX_MULTICALL_REFCOUNT
#endif
@@ -921,3 +921,134 @@ forprimes (SV* block, IN SV* svbeg, IN SV* svend = 0)
#endif
XSRETURN_UNDEF;
}
+
+void
+forcomposites (SV* block, IN SV* svbeg, IN SV* svend = 0)
+ PROTOTYPE: &$;$
+ CODE:
+ {
+ UV i, beg, end;
+ GV *gv;
+ HV *stash;
+ SV* svarg; /* We use svarg to prevent clobbering $_ outside the block */
+ CV *cv = sv_2cv(block, &stash, &gv, 0);
+
+ if (cv == Nullcv)
+ croak("Not a subroutine reference");
+ if (items <= 1) XSRETURN_UNDEF;
+
+ if (!_validate_int(svbeg, 0) || (items >= 3 && !_validate_int(svend,0))) {
+ dSP;
+ PUSHMARK(SP);
+ XPUSHs(block); XPUSHs(svbeg); XPUSHs(svend);
+ PUTBACK;
+ (void) call_pv("Math::Prime::Util::_generic_forcomposites", G_VOID|G_DISCARD);
+ SPAGAIN;
+ XSRETURN_UNDEF;
+ }
+
+ if (items < 3) {
+ beg = 4;
+ set_val_from_sv(end, svbeg);
+ } else {
+ set_val_from_sv(beg, svbeg);
+ set_val_from_sv(end, svend);
+ if (beg < 4) beg = 4;
+ }
+ if (beg > end)
+ XSRETURN_UNDEF;
+
+ /* TODO: segment sieve and walk composites */
+ SAVESPTR(GvSV(PL_defgv));
+ svarg = newSVuv(0);
+#if USE_MULTICALL
+ if (!CvISXSUB(cv)) {
+ dMULTICALL;
+ I32 gimme = G_VOID;
+ PUSH_MULTICALL(cv);
+ for (i = beg; i <= end; i++) {
+ if (!_XS_is_prob_prime(i)) {
+ sv_setuv(svarg, i);
+ GvSV(PL_defgv) = svarg;
+ MULTICALL;
+ }
+ }
+ FIX_MULTICALL_REFCOUNT;
+ POP_MULTICALL;
+ }
+ else
+#endif
+ {
+ for (i = beg; i <= end; i++) {
+ if (!_XS_is_prob_prime(i)) {
+ dSP;
+ sv_setuv(svarg, i);
+ GvSV(PL_defgv) = svarg;
+ PUSHMARK(SP);
+ call_sv((SV*)cv, G_VOID|G_DISCARD);
+ }
+ }
+ }
+ SvREFCNT_dec(svarg);
+ XSRETURN_UNDEF;
+ }
+
+void
+fordivisors (SV* block, IN SV* svn)
+ PROTOTYPE: &$
+ CODE:
+ {
+ UV i, n, ndivisors;
+ UV *divs;
+ GV *gv;
+ HV *stash;
+ SV* svarg; /* We use svarg to prevent clobbering $_ outside the block */
+ CV *cv = sv_2cv(block, &stash, &gv, 0);
+
+ if (cv == Nullcv)
+ croak("Not a subroutine reference");
+ if (items <= 1) XSRETURN_UNDEF;
+
+ if (!_validate_int(svn, 0)) {
+ dSP;
+ PUSHMARK(SP);
+ XPUSHs(block); XPUSHs(svn);
+ PUTBACK;
+ (void) call_pv("Math::Prime::Util::_generic_fordivisors", G_VOID|G_DISCARD);
+ SPAGAIN;
+ XSRETURN_UNDEF;
+ }
+
+ set_val_from_sv(n, svn);
+ divs = _divisor_list(n, &ndivisors);
+
+ SAVESPTR(GvSV(PL_defgv));
+ svarg = newSVuv(0);
+#if USE_MULTICALL
+ if (!CvISXSUB(cv)) {
+ dMULTICALL;
+ I32 gimme = G_VOID;
+ PUSH_MULTICALL(cv);
+ for (i = 0; i < ndivisors; i++) {
+ sv_setuv(svarg, divs[i]);
+ GvSV(PL_defgv) = svarg;
+ MULTICALL;
+ }
+ FIX_MULTICALL_REFCOUNT;
+ POP_MULTICALL;
+ }
+ else
+#endif
+ {
+ for (i = 0; i < ndivisors; i++) {
+ dSP;
+ sv_setuv(svarg, divs[i]);
+ GvSV(PL_defgv) = svarg;
+ PUSHMARK(SP);
+ call_sv((SV*)cv, G_VOID|G_DISCARD);
+ }
+ }
+ SvREFCNT_dec(svarg);
+ Safefree(divs);
+ XSRETURN_UNDEF;
+ }
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index afab153..75f17bb 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -27,7 +27,7 @@ our @EXPORT_OK =
miller_rabin miller_rabin_random
lucas_sequence
primes
- forprimes forcomposites
+ forprimes forcomposites fordivisors
prime_iterator prime_iterator_object
next_prime prev_prime
prime_count
@@ -106,6 +106,8 @@ BEGIN {
*prev_prime = \&Math::Prime::Util::_generic_prev_prime;
*exp_mangoldt = \&Math::Prime::Util::_generic_exp_mangoldt;
*forprimes = sub (&$;$) { _generic_forprimes(@_); }; ## no critic qw(ProhibitSubroutinePrototypes)
+ *fordivisors = sub (&$) { _generic_fordivisors(@_); }; ## no critic qw(ProhibitSubroutinePrototypes)
+ *forcomposites = sub (&$) { _generic_forcomposites(@_); }; ## no critic qw(ProhibitSubroutinePrototypes)
*_prime_memfreeall = \&Math::Prime::Util::PP::_prime_memfreeall;
*prime_memfree = \&Math::Prime::Util::PP::prime_memfree;
@@ -1446,8 +1448,11 @@ sub divisor_sum {
if (defined $k && ref($k) eq 'CODE') {
my $sum = $n-$n;
- foreach my $f (divisors($n)) {
- $sum += $k->($f);
+ if (ref($n) eq 'Math::BigInt') {
+ # If the original number was a bigint, make sure all divisors are.
+ fordivisors { $sum += $k->(Math::BigInt->new("$_")); } $n;
+ } else {
+ fordivisors { $sum += $k->($_); } $n;
}
return $sum;
}
@@ -1519,7 +1524,7 @@ sub _generic_forprimes (&$;$) { ## no critic qw(ProhibitSubroutinePrototypes)
}
}
-sub forcomposites(&$;$) { ## no critic qw(ProhibitSubroutinePrototypes)
+sub _generic_forcomposites(&$;$) { ## no critic qw(ProhibitSubroutinePrototypes)
my($sub, $beg, $end) = @_;
if (!defined $end) { $end = $beg; $beg = 4; }
_validate_num($beg) || _validate_positive_integer($beg);
@@ -1538,6 +1543,20 @@ sub forcomposites(&$;$) { ## no critic qw(ProhibitSubroutinePrototypes)
}
}
+sub _generic_fordivisors (&$) { ## no critic qw(ProhibitSubroutinePrototypes)
+ my($sub, $n) = @_;
+ _validate_num($n) || _validate_positive_integer($n);
+ my @divisors = ($n <= $_XS_MAXVAL) ? _XS_divisors($n) : divisors($n);
+ {
+ my $pp;
+ local *_ = \$pp;
+ foreach my $d (@divisors) {
+ $pp = $d;
+ $sub->();
+ }
+ }
+}
+
sub prime_iterator {
my($start) = @_;
$start = 0 unless defined $start;
@@ -1693,6 +1712,7 @@ sub znorder {
$k *= $pi;
}
}
+ $k = int($k->bstr) if $k->bacmp(''.~0) <= 0;
return $k;
# Method 3:
@@ -2549,6 +2569,8 @@ Version 0.35
# You can do something for every prime in a range. Twin primes to 10k:
forprimes { say if is_prime($_+2) } 10000;
+ # Or for the composites in a range
+ forcomposites { say if is_strong_pseudoprime($_,2) } 10000, 10**6;
# For non-bigints, is_prime and is_prob_prime will always be 0 or 2.
# They return 0 (composite), 2 (prime), or 1 (probably prime)
@@ -2601,6 +2623,8 @@ Version 0.35
# Get all divisors other than 1 and n
@divisors = divisors( $n );
+ # Or just apply a block for each one
+ fordivisors { $sum += $_ + $_*$_ } $n;
# Euler phi (Euler's totient) on a large number
use bigint; say euler_phi( 801294088771394680000412 );
@@ -2885,6 +2909,16 @@ exceptions. Here is a clumsy L</forprimes> exception example:
my $n = 0+$@;
+=head2 forcomposites
+
+ forcomposites { say } 1000;
+ forcomposites { say } 2000,2020;
+
+Given a block and either an end number or a start and end pair, calls the
+block for each composite in the inclusive range. Starting at 2, the
+composites are the non-primes (C<0> and C<1> are neither prime nor composite).
+
+
=head2 prime_iterator
my $it = prime_iterator;
@@ -3503,6 +3537,14 @@ The following conditions must hold:
- C<< n >= 2 >>
+=head2 fordivisors
+
+ fordivisors { $prod *= $_ } $n;
+
+Given a block and a non-negative number C<n>, the block is called with
+C<$_> set to each divisor in sorted order. Also see L</divisor_sum>.
+
+
=head2 moebius
say "$n is square free" if moebius($n) != 0;
@@ -4441,6 +4483,12 @@ Produce the C<matches> result from L<Math::Factor::XS> without skipping:
my @matches = grep { $_->[0] * $_->[1] == $n && $_->[0] > 1 }
combinations_with_repetition( [divisors($n)], 2 );
+Compute L<OEIS A054903|http://oeis.org/A054903> just like CRG4's Pari example:
+
+ use Math::Prime::Util qw/forcomposite divisor_sum/;
+ forcomposites {
+ say if divisor_sum($_)+6 == divisor_sum($_+6)
+ } 9,1e7;
=head1 PRIMALITY TESTING NOTES
@@ -4724,6 +4772,11 @@ faster on average in MPU.
Similar to MPU's L</divisors>.
+=item C<forprime>, C<forcomposite>, C<fordiv>, C<sumdiv>
+
+Similar to MPU's L</forprimes>, L</forcomposites>, L<fordivisors>, and
+L<divisor_sum>.
+
=item C<eulerphi>
Similar to MPU's L</euler_phi>. MPU is 2-20x faster for native integers.
diff --git a/t/19-moebius.t b/t/19-moebius.t
index fc0efd7..451ecc1 100644
--- a/t/19-moebius.t
+++ b/t/19-moebius.t
@@ -96,7 +96,7 @@ my %jordan_totients = (
# A059377
4 => [1, 15, 80, 240, 624, 1200, 2400, 3840, 6480, 9360, 14640, 19200, 28560, 36000, 49920, 61440, 83520, 97200, 130320, 149760, 192000, 219600, 279840, 307200, 390000, 428400, 524880, 576000, 707280, 748800, 923520, 983040, 1171200],
# A059378
- 5 => [1, 31, 242, 992, 3124, 7502, 16806, 31744, 58806, 96844, 161050, 240064, 371292, 520986, 756008, 1015808, 1419856, 1822986, 2476098, 3099008, 4067052, 4992550, 6436342, 7682048, 9762500, 11510052, 14289858, 16671552, 20511148],
+ 5 => [1, 31, 242, 992, 3124, 7502, 16806, 31744, 58806, 96844, 161050, 240064, 371292, 520986, 756008, 1015808, 1419856, 1822986, 2476098, 3099008, 4067052, 4992550, 6436342, 7682048, 9762500, 11510052, 14289858, 16671552, 20511148, 23436248, 28629150, 32505856, 38974100, 44015536, 52501944, 58335552, 69343956, 76759038, 89852664, 99168256, 115856200, 126078612, 147008442, 159761600, 183709944, 199526602, 229345006, 245825536, 282458442, 302637500, 343605152, 368321664],
# A069091
6 => [1, 63, 728, 4032, 15624, 45864, 117648, 258048, 530712, 984312, 1771560, 2935296, 4826808, 7411824, 11374272, 16515072, 24137568, 33434856, 47045880, 62995968, 85647744, 111608280, 148035888, 187858944, 244125000, 304088904, 386889048],
# A069092
@@ -114,6 +114,8 @@ my %sigmak = (
3 => [1, 9, 28, 73, 126, 252, 344, 585, 757, 1134, 1332, 2044, 2198, 3096, 3528, 4681, 4914, 6813, 6860, 9198, 9632, 11988, 12168, 16380, 15751, 19782, 20440, 25112, 24390, 31752, 29792, 37449, 37296, 44226, 43344, 55261, 50654, 61740, 61544],
);
+my @tau4 = (1,4,4,10,4,16,4,20,10,16,4,40,4,16,16,35,4,40,4,40,16,16,4,80,10,16,20,40,4,64,4,56,16,16,16,100,4,16,16,80,4,64,4,40,40,16,4,140,10,40,16,40,4,80,16,80,16,16,4,160,4,16,40,84,16,64,4,40,16,64,4,200,4,16,40,40,16);
+
my %mangoldt = (
-13 => 1,
0 => 1,
@@ -219,9 +221,9 @@ plan tests => 0 + 1
+ scalar(@mult_orders)
+ scalar(keys %jordan_totients)
+ 2 # Dedekind psi calculated two ways
- + 1 # Calculate J5 two different ways
+ + 2 # Calculate J5 two different ways
+ 2 * $use64 # Jordan totient example
- + 1 + 2*scalar(keys %sigmak) + 2
+ + 1 + 2*scalar(keys %sigmak) + 3
+ scalar(keys %mangoldt)
+ scalar(keys %chebyshev1)
+ scalar(keys %chebyshev2)
@@ -270,10 +272,7 @@ is_deeply( [euler_phi(10,20)], [4,10,4,12,6,8,8,16,6,18,8], "euler_phi 10-20" );
###### Jordan Totient
while (my($k, $tref) = each (%jordan_totients)) {
- my @tlist;
- foreach my $n (1 .. scalar @$tref) {
- push @tlist, jordan_totient($k, $n);
- }
+ my @tlist = map { jordan_totient($k, $_) } 1 .. scalar @$tref;
is_deeply( \@tlist, $tref, "Jordan's Totient J_$k" );
}
@@ -289,9 +288,11 @@ while (my($k, $tref) = each (%jordan_totients)) {
}
{
- my @J5_moebius = map { divisor_sum($_, sub { my $d=shift; $d**5 * moebius($_/$d); }) } (1 .. 100);
- my @J5_jordan = map { jordan_totient(5, $_) } (1 .. 100);
- is_deeply( \@J5_moebius, \@J5_jordan, "Calculate J5 two different ways");
+ my $J5 = $jordan_totients{5};
+ my @J5_jordan = map { jordan_totient(5, $_) } 1 .. scalar @$J5;
+ is_deeply( \@J5_jordan, $J5, "Jordan totient 5, using jordan_totient");
+ my @J5_moebius = map { my $n = $_; divisor_sum($n, sub { my $d=shift; $d**5 * moebius($n/$d); }) } 1 .. scalar @$J5;
+ is_deeply( \@J5_moebius, $J5, "Jordan totient 5, using divisor sum" );
}
if ($use64) {
@@ -330,6 +331,15 @@ while (my($k, $sigmaref) = each (%sigmak)) {
is_deeply( \@slist2, $sigmak{0}, "tau as divisor_sum(n, 0)" );
}
+{
+ # tau_4 A007426
+ my @t;
+ foreach my $n (1 .. scalar @tau4) {
+ push @t, divisor_sum($n, sub { divisor_sum($_[0],sub { divisor_sum($_[0],0) }) });
+ }
+ is_deeply( \@t, \@tau4, "Tau4 (A007426), nested divisor sums" );
+}
+
###### Exponential of von Mangoldt
while (my($n, $em) = each (%mangoldt)) {
is( exp_mangoldt($n), $em, "exp_mangoldt($n) == $em" );
diff --git a/t/32-iterators.t b/t/32-iterators.t
index 3ef3a42..9823f30 100644
--- a/t/32-iterators.t
+++ b/t/32-iterators.t
@@ -4,7 +4,8 @@ use warnings;
use Test::More;
use Math::Prime::Util qw/primes prev_prime next_prime
- forprimes prime_iterator prime_iterator_object/;
+ forprimes forcomposites fordivisors
+ prime_iterator prime_iterator_object/;
use Math::BigInt try => "GMP,Pari";
use Math::BigFloat;
@@ -13,6 +14,8 @@ my $broken64 = (18446744073709550592 == ~0);
plan tests => 8 # forprimes errors
+ 12 + 5 # forprimes simple
+ + 1 # forcomposites simple
+ + 2 # fordivisors simple
+ 3 # iterator errors
+ 7 # iterator simple
+ 2 # forprimes/iterator nesting
@@ -61,6 +64,23 @@ ok(!eval { forprimes { 1 } 5.6; }, "forprimes abc");
{ my @t; forprimes { push @t, $_ } 31398, 31468;
is_deeply( [@t], [], "forprimes 31398,31468 (empty region)" );
}
+{ my @t; forcomposites { push @t, $_ } 50;
+ is_deeply( [@t], [qw/4 6 8 9 10 12 14 15 16 18 20 21 22 24 25 26 27 28 30 32 33 34 35 36 38 39 40 42 44 45 46 48 49 50/], "forcomposites 50" );
+}
+{
+ my $a = 0;
+ fordivisors { $a += $_ + $_*$_ } 54321;
+ is($a, 3287796520, "fordivisors: d|54321: a+=d+d^2");
+ # Matches Math::Pari:
+ # my $a = PARI(0); my $j; fordiv(54321,$j,sub { $a += $j + $j**2 });
+}
+{
+ # Pari: v=List(); for(n=1, 50, fordiv(n, d, listput(v, d))); Vec(v)
+ my @A027750 = (1,1,2,1,3,1,2,4,1,5,1,2,3,6,1,7,1,2,4,8,1,3,9,1,2,5,10,1,11,1,2,3,4,6,12,1,13,1,2,7,14,1,3,5,15,1,2,4,8,16,1,17,1,2,3,6,9,18,1,19,1,2,4,5,10,20,1,3,7,21,1,2,11,22,1,23,1,2,3,4,6,8,12,24,1,5,25,1,2,13,26,1,3,9,27,1,2,4,7,14,28,1,29,1,2,3,5,6,10,15,30,1,31,1,2,4,8,16,32,1,3,11,33,1,2,17,34,1,5,7,35,1,2,3,4,6,9,12,18,36,1,37,1,2,19,38,1,3,13,39,1,2,4,5,8,10,20,40,1,41,1,2,3,6,7,14,21,42,1,43,1,2,4,11,22,44,1,3,5,9,15,45,1,2,23,46,1,47,1,2,3,4,6,8,12,16,24,48,1,7,49,1,2,5,10 [...]
+ my @a;
+ do { fordivisors { push @a, $_ } $_ } for 1..50;
+ is_deeply(\@a, \@A027750, "A027750 using fordivisors");
+}
ok(!eval { prime_iterator(-2); }, "iterator -2");
ok(!eval { prime_iterator("abc"); }, "iterator abc");
--
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