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