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

_______________________________________________
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

Reply via email to