[libmath-prime-util-perl] 17/181: _divisor_list in XS, use it. Add simple 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 c6a1930985c483639261101caf4e3a8d8d62740b
Author: Dana Jacobsen <dana at acm.org>
Date: Mon Dec 16 17:14:34 2013 -0800
_divisor_list in XS, use it. Add simple forcomposites.
---
Changes | 10 ++++++++
TODO | 6 +++++
XS.xs | 13 ++++++++++
factor.c | 70 ++++++++++++++++++++++++++++++++++++++++++++++++++
factor.h | 3 ++-
lib/Math/Prime/Util.pm | 60 +++++++++++++++++++++++++------------------
6 files changed, 136 insertions(+), 26 deletions(-)
diff --git a/Changes b/Changes
index 109c4f3..8781b34 100644
--- a/Changes
+++ b/Changes
@@ -2,6 +2,12 @@ Revision history for Perl module Math::Prime::Util
0.36 2013-12
+ [ADDED]
+
+ - forcomposites like forprimes, but on composites. See Pari 2.6.x.
+
+ [FUNCTIONALITY AND PERFORMANCE]
+
- Some 5.6.2-is-broken workarounds.
- Some LMO edge cases: small numbers that only show up if a #define is
@@ -10,6 +16,10 @@ Revision history for Perl module Math::Prime::Util
- LMO much faster if -march=native is used for gcc on a platform with
asm popcnt (e.g. Nahalem+, Barcelona+, ARM Neon, SPARC, Power7, etc.).
+ - all_factors is faster for native size numbers with many factors. This
+ is mostly in preparation for a alldivisors function, and moving a few
+ more functions to XS/Perl from Perl/XS to cut overhead.
+
0.35 2013-12-08
diff --git a/TODO b/TODO
index 6c8682a..0be771d 100644
--- a/TODO
+++ b/TODO
@@ -61,3 +61,9 @@
- look at sieve.c style prime walking
- Add Inverse Li and Legendre Phi to API?
+
+- Document forcomposites. Consider XS for it.
+
+- Add fordivisors as XS.
+
+- Make divisors() an alias for all_factors().
diff --git a/XS.xs b/XS.xs
index 4832367..cb98997 100644
--- a/XS.xs
+++ b/XS.xs
@@ -449,6 +449,19 @@ _XS_factor_exp(IN UV n)
}
void
+_XS_divisors(IN UV n)
+ PPCODE:
+ if (GIMME_V == G_SCALAR) {
+ PUSHs(sv_2mortal(newSVuv( _XS_divisor_sum(n, 0) )));
+ } else {
+ UV i, ndivisors;
+ UV* divs = _divisor_list(n, &ndivisors);
+ for (i = 0; i < ndivisors; i++)
+ XPUSHs(sv_2mortal(newSVuv(divs[i])));
+ Safefree(divs);
+ }
+
+void
trial_factor(IN UV n, IN UV maxfactor = 0)
PPCODE:
SIMPLE_FACTOR(trial_factor, n, maxfactor);
diff --git a/factor.c b/factor.c
index a425a38..0af17bd 100644
--- a/factor.c
+++ b/factor.c
@@ -292,6 +292,76 @@ static int is_perfect_square(UV n, UV* sqrtn)
return 1;
}
+static int _divisors_from_factors(UV v, UV npe, UV* fp, UV* fe, UV* res) {
+ UV p, e, i;
+ if (npe == 0) return 0;
+ p = *fp++;
+ e = *fe++;
+ if (npe == 1) {
+ for (i = 0; i <= e; i++) {
+ *res++ = v;
+ v *= p;
+ }
+ return e+1;
+ } else {
+ int nret = 0;;
+ for (i = 0; i <= e; i++) {
+ int nres = _divisors_from_factors(v, npe-1, fp, fe, res);
+ v *= p;
+ res += nres;
+ nret += nres;
+ }
+ return nret;
+ }
+}
+
+UV* _divisor_list(UV n, UV *num_divisors)
+{
+ UV factors[MPU_MAX_FACTORS+1];
+ UV exponents[MPU_MAX_FACTORS+1];
+ UV* divs;
+ int i, j, nfactors, ndivisors;
+
+ if (n <= 1) {
+ New(0, divs, 1, UV);
+ divs[0] = 1;
+ *num_divisors = n;
+ return divs;
+ }
+ /* Factor and convert to factor/exponent pair */
+ nfactors = factor(n, factors);
+ exponents[0] = 1;
+ for (i = 1, j = 1; i < nfactors; i++) {
+ if (factors[i] != factors[i-1]) {
+ exponents[j] = 1;
+ factors[j++] = factors[i];
+ } else {
+ exponents[j-1]++;
+ }
+ }
+ nfactors = j;
+ /* Calculate number of divisors, allocate space, fill with divisors */
+ ndivisors = exponents[0] + 1;
+ for (i = 1; i < nfactors; i++)
+ ndivisors *= (exponents[i] + 1);
+ New(0, divs, ndivisors, UV);
+ (void) _divisors_from_factors(1, nfactors, factors, exponents, divs);
+ { /* Sort (Shell sort is easy and efficient) */
+ static UV gaps[] = {301, 132, 57, 23, 10, 4, 1, 0};
+ UV gap, gapi = 0;
+ for (gap = gaps[gapi]; gap > 0; gap = gaps[++gapi]) {
+ for (i = gap; i < ndivisors; i++) {
+ UV v = divs[i];
+ for (j = i; j >= gap && divs[j-gap] > v; j -= gap)
+ divs[j] = divs[j-gap];
+ divs[j] = v;
+ }
+ }
+ }
+ *num_divisors = ndivisors;
+ return divs;
+}
+
/*
* The usual method, on OEIS for instance, is:
diff --git a/factor.h b/factor.h
index 06ea83d..862fb65 100644
--- a/factor.h
+++ b/factor.h
@@ -18,6 +18,7 @@ extern int pplus1_factor(UV n, UV *factors, UV B);
extern int squfof_factor(UV n, UV *factors, UV rounds);
extern int racing_squfof_factor(UV n, UV *factors, UV rounds);
-extern UV _XS_divisor_sum(UV n, UV k);
+extern UV _XS_divisor_sum(UV n, UV k);
+extern UV* _divisor_list(UV n, UV *num_divisors);
#endif
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index 0bae93f..7c84a53 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -27,7 +27,8 @@ our @EXPORT_OK =
miller_rabin miller_rabin_random
lucas_sequence
primes
- forprimes prime_iterator prime_iterator_object
+ forprimes forcomposites
+ prime_iterator prime_iterator_object
next_prime prev_prime
prime_count
prime_count_lower prime_count_upper prime_count_approx
@@ -1218,41 +1219,31 @@ sub consecutive_integer_lcm {
return $pn;
}
-
sub all_factors {
my $n = shift;
# In scalar context, returns sigma_0(n). Very fast.
return divisor_sum($n,0) unless wantarray;
- my $use_bigint = !($n < $_Config{'maxparam'} || int($n) eq $_Config{'maxparam'});
- my @factors = factor($n);
- if ($n <= 0) { @factors = (); return @factors; }
+ _validate_num($n) || _validate_positive_integer($n);
+
+ return _XS_divisors($n) if $n <= $_XS_MAXVAL;
+
my %all_factors;
- if ($use_bigint) {
- foreach my $f1 (@factors) {
- next if $f1 >= $n;
- my $big_f1 = Math::BigInt->new("$f1");
- my @to_add = map { ($_ <= ~0) ? int($_->bstr) : $_ }
- grep { $_ < $n }
- map { $big_f1 * $_ }
- keys %all_factors;
- undef @all_factors{ $f1, @to_add };
- }
- } else {
- foreach my $f1 (@factors) {
- next if $f1 >= $n;
- my @to_add = grep { $_ < $n }
- map { int($f1 * $_) }
- keys %all_factors;
- undef @all_factors{ $f1, @to_add };
- }
+ foreach my $f1 (factor($n)) {
+ next if $f1 >= $n;
+ my $big_f1 = Math::BigInt->new("$f1");
+ my @to_add = map { ($_ <= ~0) ? int($_->bstr) : $_ }
+ grep { $_ < $n }
+ map { $big_f1 * $_ }
+ keys %all_factors;
+ undef @all_factors{ $f1, @to_add };
}
# Add 1 and n
undef $all_factors{1};
undef $all_factors{$n};
- @factors = sort {$a<=>$b} keys %all_factors;
- return @factors;
+ my @divisors = sort {$a<=>$b} keys %all_factors;
+ return @divisors;
}
@@ -1523,6 +1514,25 @@ sub _generic_forprimes (&$;$) { ## no critic qw(ProhibitSubroutinePrototypes)
}
}
+sub forcomposites(&$;$) { ## no critic qw(ProhibitSubroutinePrototypes)
+ my($sub, $beg, $end) = @_;
+ if (!defined $end) { $end = $beg; $beg = 4; }
+ _validate_num($beg) || _validate_positive_integer($beg);
+ _validate_num($end) || _validate_positive_integer($end);
+ $beg = 4 if $beg < 4;
+ $end = Math::BigInt->new(''.~0) if ref($end) ne 'Math::BigInt' && $end == ~0;
+ {
+ my $pp;
+ local *_ = \$pp;
+ for ( ; $beg <= $end ; $beg++ ) {
+ if (!is_prime($beg)) {
+ $pp = $beg;
+ $sub->();
+ }
+ }
+ }
+}
+
sub prime_iterator {
my($start) = @_;
$start = 0 unless defined $start;
--
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