[libmath-prime-util-perl] 31/55: forpart better interface, add docs, PP, tests
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 7ac02a4b39ad009558cb0c0d1507ca536bbe47f8
Author: Dana Jacobsen <dana at acm.org>
Date: Sat May 10 17:18:16 2014 -0700
forpart better interface, add docs, PP, tests
---
Changes | 1 +
TODO | 10 +++---
XS.xs | 87 ++++++++++++++++++++-------------------------
lib/Math/Prime/Util.pm | 28 ++++++++++++++-
lib/Math/Prime/Util/PP.pm | 38 ++++++++++++++++++++
lib/Math/Prime/Util/PPFE.pm | 4 +++
t/24-partitions.t | 64 +++++++++++++++++++++++++++++++--
t/92-release-pod-coverage.t | 2 +-
8 files changed, 176 insertions(+), 58 deletions(-)
diff --git a/Changes b/Changes
index 8d49f17..5352c1b 100644
--- a/Changes
+++ b/Changes
@@ -6,6 +6,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)
[FUNCTIONALITY AND PERFORMANCE]
diff --git a/TODO b/TODO
index 5743d8e..2f0bb1c 100644
--- a/TODO
+++ b/TODO
@@ -73,8 +73,8 @@
to speed it up using some ideas from the Ohana 2011 SAGE branch. For example
(10**13,10**5) takes 2.5x longer, albeit with 6x less memory.
-- forpart:
- - alias with fordivisors to save space
- - consider use of const_int since so many entries will be 1 or 2
- - documentation (mention for partitions function, mention Pari 2.7)
- - tests
+- XS alias forprimes, forcomposites, forpart
+
+- More Pari: vecsum, parforprime, binomial
+
+- is_power should take a third argument that gets set to the root
diff --git a/XS.xs b/XS.xs
index d56da8f..21c1f3a 100644
--- a/XS.xs
+++ b/XS.xs
@@ -208,28 +208,6 @@ static int _vcallsubn(pTHX_ I32 flags, I32 stashflags, const char* name, int nar
else { PUSHs(sv_2mortal(newSViv(r_))); } \
} while (0)
-static void part_setminmax(UV* min, UV* max, SV* sv, const char* text) {
- if (!sv) return;
- if (SvIOK(sv)) {
- *min = 0;
- *max = my_svuv(sv);
- } else if (SvROK(sv) && SvTYPE(SvRV(sv)) == SVt_PVAV) {
- SV** svmin;
- SV** svmax;
- AV* array_a = (AV*) SvRV(sv);
- if (av_len(array_a)+1 != 2)
- croak("%s must have exactly 2 elements", text);
- svmin = av_fetch(array_a, 0, 0);
- svmax = av_fetch(array_a, 1, 0);
- if (!svmin || !svmax || !SvIOK(*svmin) || !SvIOK(*svmax))
- croak("Invalid data in %s", text);
- *min = my_svuv(*svmin);
- *max = my_svuv(*svmax);
- } else {
- /* Allow fall through for non-existant arg */
- /* croak("Invalid data in %s", text); */
- }
-}
MODULE = Math::Prime::Util PACKAGE = Math::Prime::Util
@@ -1252,11 +1230,10 @@ fordivisors (SV* block, IN SV* svn)
Safefree(divs);
void
-forpart (SV* block, IN SV* svn, IN SV* svra = 0, IN SV* svrn = 0)
- PROTOTYPE: &$;$$
+forpart (SV* block, IN SV* svn, IN SV* svh = 0)
+ PROTOTYPE: &$;$
PREINIT:
- UV i, m, h, n, mina, maxa, minn, maxn;
- UV *x;
+ UV i, n, amin, amax, nmin, nmax;
GV *gv;
HV *stash;
CV *cv;
@@ -1266,7 +1243,7 @@ forpart (SV* block, IN SV* svn, IN SV* svra = 0, IN SV* svrn = 0)
if (cv == Nullcv)
croak("Not a subroutine reference");
if (!_validate_int(aTHX_ svn, 0)) {
- _vcallsubn(aTHX_ G_VOID|G_DISCARD, VCALL_ROOT, "_generic_forpart", 2);
+ _vcallsub_with_pp("forpart");
return;
}
n = my_svuv(svn);
@@ -1277,42 +1254,54 @@ forpart (SV* block, IN SV* svn, IN SV* svra = 0, IN SV* svrn = 0)
SvREADONLY_on(svals[i]);
}
- /* ZS1 algorithm from Zoghbi and Stojmenovic 1998) */
- New(0, x, n+1, UV);
- for(i = 0; i <= n; i++) x[i] = 1;
- x[1] = n;
- m = (n > 0) ? 1 : 0; /* n=0 => one call with empty list */
- h = 1;
+ amin = 0; amax = n; nmin = 0; nmax = n;
+ if (svh != 0) {
+ HV* rhash;
+ SV** svp;
+ if (!SvROK(svh) || SvTYPE(SvRV(svh)) != SVt_PVHV)
+ croak("forpart second argument must be a hash reference");
+ rhash = (HV*) SvRV(svh);
+ if ((svp = hv_fetchs(rhash, "n", 0)) != NULL)
+ { nmin = my_svuv(*svp); nmax = nmin; }
+ if ((svp = hv_fetchs(rhash, "amin", 0)) != NULL) amin = my_svuv(*svp);
+ if ((svp = hv_fetchs(rhash, "amax", 0)) != NULL) amax = my_svuv(*svp);
+ if ((svp = hv_fetchs(rhash, "nmin", 0)) != NULL) nmin = my_svuv(*svp);
+ if ((svp = hv_fetchs(rhash, "nmax", 0)) != NULL) nmax = my_svuv(*svp);
+ if (amin < 1) amin = 1;
+ if (amax > n) amax = n;
+ if (nmax > n) nmax = n;
+ }
- mina = 0; maxa = n; minn = 0; maxn = n;
- if (svra != 0) part_setminmax( &mina, &maxa, svra, "forpart A arg" );
- if (svrn != 0) part_setminmax( &minn, &maxn, svrn, "forpart N arg" );
+ { /* ZS1 algorithm from Zoghbi and Stojmenovic 1998) */
+ UV *x, m, h;
+ New(0, x, n+1, UV);
+ for(i = 2; i <= n; i++) x[i] = 1;
+ x[1] = n;
+ m = (n > 0) ? 1 : 0; /* n=0 => one call with empty list */
+ h = 1;
- {
+ /* We could add some optimizations for amin, amax, nmin, nmax.
+ * Pari is seriously fast for some restrictions. */
while (1) {
- if (m >= minn && m <= maxn && x[m] >= mina && x[1] <= maxa)
+ if (m >= nmin && m <= nmax && x[m] >= amin && x[1] <= amax)
{ dSP; ENTER; PUSHMARK(SP);
EXTEND(SP, m); for (i=1; i <= m; i++) { PUSHs(svals[x[i]]); }
PUTBACK; call_sv((SV*)cv, G_VOID|G_DISCARD); LEAVE;
}
- if (x[1] <= 1)
- break;
+ if (x[1] <= 1) break;
if (x[h] == 2) {
- m++;
- x[h] = 1;
- h--;
+ m++; x[h--] = 1;
} else {
- UV r,t;
- r = x[h]-1;
- t = m-h+1;
+ UV r = x[h]-1;
+ UV t = m-h+1;
x[h] = r;
- while (t >= r) { h++; x[h] = r; t -= r; }
+ while (t >= r) { x[++h] = r; t -= r; }
if (t == 0) { m = h; }
- else { m = h+1; if (t > 1) { h++; x[h] = t; } }
+ else { m = h+1; if (t > 1) { x[++h] = t; } }
}
}
+ Safefree(x);
}
- Safefree(x);
for (i = 0; i <= n; i++)
SvREFCNT_dec(svals[i]);
Safefree(svals);
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index 14e7900..ff51a04 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -813,7 +813,7 @@ __END__
=encoding utf8
-=for stopwords forprimes forcomposites fordivisors Möbius Deléglise totient moebius mertens liouville znorder irand primesieve uniqued k-tuples von SoE pari yafu fonction qui compte le nombre nombres voor PhD superset sqrt(N) gcd(A^M k-th (10001st primegen libtommath kronecker znprimroot znlog gcd lcm invmod
+=for stopwords forprimes forcomposites fordivisors forpart Möbius Deléglise totient moebius mertens liouville znorder irand primesieve uniqued k-tuples von SoE pari yafu fonction qui compte le nombre nombres voor PhD superset sqrt(N) gcd(A^M k-th (10001st primegen libtommath kronecker znprimroot znlog gcd lcm invmod
=head1 NAME
@@ -916,6 +916,8 @@ Version 0.41
say chebyshev_psi(234984);
say chebyshev_theta(92384234);
say partitions(1000);
+ # Show all prime partitions of 25
+ forpart { say "@_" unless scalar grep { !is_prime($_) } @_ } 25;
# divisor sum
$sigma = divisor_sum( $n ); # sum of divisors
@@ -1177,6 +1179,30 @@ which are not prime: C<4, 6, 8, 9, 10, 12, 14, 15, ...>
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 forpart
+
+ forpart { say "@_" } 25; # unrestricted partitions
+ forpart { say "@_" } 25,{n=>5} # ... with exactly 5 values
+ forpart { say "@_" } 25,{nmax=>5} # ... with <=5 values
+
+Given a non-negative number C<n>, the block is called with C<@_> set to
+the array of additive integer partitions. The operation is very similar
+to the C<forpart> function in Pari/GP 2.6.x, though the ordering is
+different. The algorithm is ZS1 from Zoghbi and Stojmenović (1998), hence
+the ordering is identical to that of L<Integer::Partition>.
+Use L</partitions> to get just the count of unrestricted partitions.
+
+
+An optional hash reference may be given to produce restricted partitions.
+Each value must be a non-negative integer. The allowable keys are:
+
+ n restrict to exactly this many values
+ amin all elements must be at least this value
+ amax all elements must be at most this value
+ nmin the array must have at least this many values
+ nmax the array must have at most this many values
+
=head2 prime_iterator
diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm
index ccdf658..8a133a0 100644
--- a/lib/Math/Prime/Util/PP.pm
+++ b/lib/Math/Prime/Util/PP.pm
@@ -3542,6 +3542,44 @@ sub RiemannR {
return $sum;
}
+sub forpart {
+ my($sub, $n, $rhash) = @_;
+ _validate_positive_integer($n);
+ my($mina, $maxa, $minn, $maxn) = (0,$n,0,$n);
+ if (defined $rhash) {
+ croak "forpart second argument must be a hash reference"
+ unless ref($rhash) eq 'HASH';
+ $minn = $maxn = $rhash->{n} if defined $rhash->{n};
+ $mina = $rhash->{amin} if defined $rhash->{amin};
+ $maxa = $rhash->{amax} if defined $rhash->{amax};
+ $minn = $rhash->{nmin} if defined $rhash->{nmin};
+ $maxn = $rhash->{nmax} if defined $rhash->{nmax};
+ _validate_positive_integer($mina);
+ _validate_positive_integer($maxa);
+ _validate_positive_integer($minn);
+ _validate_positive_integer($maxn);
+ }
+ my @x = (1) x ($n+1);
+ my $m = ($n > 0) ? 1 : 0;
+ my $h = 1;
+ $x[1] = $n;
+ while (1) {
+ $sub->(@x[1 .. $m])
+ if $m >= $minn && $m <= $maxn && $x[$m] >= $mina && $x[1] <= $maxa;
+ last if $x[1] <= 1;
+ if ($x[$h] == 2) {
+ $m++; $x[$h] = 1; $h--; }
+ else {
+ my $r = $x[$h] - 1;
+ my $t = $m - $h + 1;
+ $x[$h] = $r;
+ while ($t >= $r) { $h++; $x[$h] = $r; $t -= $r; }
+ if ($t == 0) { $m = $h }
+ else { $m = $h+1; if ($t > 1) { $h++; $x[$h] = $t; } }
+ }
+ }
+}
+
1;
__END__
diff --git a/lib/Math/Prime/Util/PPFE.pm b/lib/Math/Prime/Util/PPFE.pm
index cac2d4b..0346683 100644
--- a/lib/Math/Prime/Util/PPFE.pm
+++ b/lib/Math/Prime/Util/PPFE.pm
@@ -416,6 +416,10 @@ sub fordivisors (&$) { ## no critic qw(ProhibitSubroutinePrototypes)
}
}
+sub forpart (&$;$) { ## no critic qw(ProhibitSubroutinePrototypes)
+ Math::Prime::Util::PP::forpart(@_);
+}
+
1;
__END__
diff --git a/t/24-partitions.t b/t/24-partitions.t
index cb51723..6e2cff4 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/;
+use Math::Prime::Util qw/partitions forpart/;
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));
+plan tests => scalar(@parts) + scalar(keys(%bparts)) + 14;
foreach my $n (0..$#parts) {
@@ -86,3 +86,63 @@ foreach my $n (0..$#parts) {
while (my($n, $epart) = each (%bparts)) {
is( partitions($n), $epart, "partitions($n)" );
}
+
+################### forpart
+
+{ my @p=(); forpart { push @p, [@_] } 0;
+ is_deeply( [@p], [[]], "forpart 0" ); }
+
+{ my @p=(); forpart { push @p, [@_] } 1;
+ is_deeply( [@p], [[1]], "forpart 1" ); }
+
+{ my @p=(); forpart { push @p, [@_] } 2;
+ is_deeply( [@p], [[2],[1,1]], "forpart 2" ); }
+
+{ my @p=(); forpart { push @p, [@_] } 3;
+ is_deeply( [@p], [[3],[2,1],[1,1,1]], "forpart 3" ); }
+
+{ my @p=(); forpart { push @p, [@_] } 3;
+ is_deeply( [@p], [[3],[2,1],[1,1,1]], "forpart 3" ); }
+
+{ my @p=(); forpart { push @p, [@_] } 6;
+ is_deeply( [@p], [[6],[5,1],[4,2],[4,1,1],[3,3],[3,2,1],[3,1,1,1],[2,2,2],[2,2,1,1],[2,1,1,1,1],[1,1,1,1,1,1]], "forpart 6" ); }
+
+{ my @p=(); forpart { push @p, [@_] } 17,{n=>2};
+ is_deeply( [@p], [[16,1],[15,2],[14,3],[13,4],[12,5],[11,6],[10,7],[9,8]], "forpart 17 restricted n=[2,2]" ); }
+
+{ my @p1 = (); my @p2 = ();
+ forpart { push @p1, [@_] if @_ <= 5 } 27;
+ forpart { push @p2, [@_] } 27, {nmax=>5};
+ is_deeply( [@p1], [@p2], "forpart 27 restricted nmax 5" ); }
+
+{ my @p1 = (); my @p2 = ();
+ forpart { push @p1, [@_] if @_ >= 20 } 27;
+ forpart { push @p2, [@_] } 27, {nmin=>20};
+ is_deeply( [@p1], [@p2], "forpart 27 restricted nmin 20" ); }
+
+{ my @p1 = (); my @p2 = ();
+ forpart { push @p1, [@_] if @_ >= 10 && @_ <= 13 } 19;
+ forpart { push @p2, [@_] } 19, {nmin=>10,nmax=>13};
+ is_deeply( [@p1], [@p2], "forpart 19 restricted n=[10..13]" ); }
+
+{ my @p1 = (); my @p2 = ();
+ forpart { push @p1, [@_] unless scalar grep { $_ > 4 } @_ } 20;
+ forpart { push @p2, [@_] } 20, {amax=>4};
+ is_deeply( [@p1], [@p2], "forpart 20 restricted amax 4" ); }
+
+{ my @p1 = (); my @p2 = ();
+ forpart { push @p1, [@_] unless scalar grep { $_ < 4 } @_ } 15;
+ forpart { push @p2, [@_] } 15, {amin=>4};
+ is_deeply( [@p1], [@p2], "forpart 15 restricted amin 4" ); }
+
+{ my @p1 = (); my @p2 = ();
+ forpart { push @p1, [@_] unless scalar grep { $_ < 3 || $_ > 6 } @_ } 21;
+ forpart { push @p2, [@_] } 21, {amin=>3,amax=>6};
+ is_deeply( [@p1], [@p2], "forpart 21 restricted a=[3..6]" ); }
+
+#{ my @p1 = (); my @p2 = ();
+# forpart { push @p1, [@_] unless @_ != 4 || scalar grep { $_ < 2 || $_ > 8 } @_ } 22;
+# forpart { push @p2, [@_] } 22, {amin=>2,amax=>8,n=>4};
+# 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]" ); }
diff --git a/t/92-release-pod-coverage.t b/t/92-release-pod-coverage.t
index 74eb36b..eb71b84 100644
--- a/t/92-release-pod-coverage.t
+++ b/t/92-release-pod-coverage.t
@@ -52,7 +52,7 @@ sub mpu_public_regex {
miller_rabin miller_rabin_random
lucas_sequence
primes
- forprimes forcomposites fordivisors
+ forprimes forcomposites fordivisors forpart
prime_iterator prime_iterator_object
next_prime prev_prime
prime_count
--
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