[libmath-prime-util-perl] 17/55: Another couple percent speedup for very small factor inputs

[Partha P. Mukherjee]

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 94f19015c476f0da117472af98040b9b70f5aec9
Author: Dana Jacobsen <dana at acm.org>
Date:   Fri May 2 13:40:25 2014 -0700

    Another couple percent speedup for very small factor inputs
---
 factor.c               | 26 +++++++++++++++++---------
 lib/Math/Prime/Util.pm |  7 ++++---
 2 files changed, 21 insertions(+), 12 deletions(-)

diff --git a/factor.c b/factor.c
index 5e4c35a..774ad31 100644
--- a/factor.c
+++ b/factor.c
@@ -61,24 +61,32 @@ int factor(UV n, UV *factors)
     while ( (n & 1) == 0 ) { factors[nfactors++] = 2; n /= 2; }
     while ( (n % 3) == 0 ) { factors[nfactors++] = 3; n /= 3; }
     while ( (n % 5) == 0 ) { factors[nfactors++] = 5; n /= 5; }
+  }
 
-    if (f*f <= n) {
-      UV sp = 3;
-      UV const lastsp = (n < 20000000)  ?  NPRIMES_SMALL-1  :  83;
-      while (++sp < lastsp) {
-        f = primes_small[sp];
+  if (f*f <= n) {
+    UV sp = 4, lastsp = 83;
+    while (sp < lastsp) {            /* Trial division from 7 to 421 */
+      if (f*f > n) break;
+      while ( (n%f) == 0 ) {
+        factors[nfactors++] = f;
+        n /= f;
+      }
+      f = primes_small[++sp];
+    }
+    /* If n is small and still composite, finish it here */
+    if (n < 2011*2011 && f*f <= n) {  /* Trial division from 431 to 2003 */
+      while (sp < NPRIMES_SMALL) {
         if (f*f > n) break;
         while ( (n%f) == 0 ) {
           factors[nfactors++] = f;
           n /= f;
         }
+        f = primes_small[++sp];
       }
-      f = primes_small[sp];
     }
   }
-  if (n < f*f) {
-    if (n != 1)
-      factors[nfactors++] = n;
+  if (f*f > n) {
+    if (n != 1) factors[nfactors++] = n;
     return nfactors;
   }
   /* Perfect squares and cubes.  Factor root only once. */
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index 793fd32..533e36a 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -2899,7 +2899,7 @@ Here is the right way to do PE problem 69 (under 0.03s):
   $n++ while pn_primorial($n+1) < 1000000;
   say pn_primorial($n);
 
-Project Euler, problem 187, stupid brute force solution, ~3 minutes:
+Project Euler, problem 187, stupid brute force solution, 2 to 3 minutes:
 
   use Math::Prime::Util qw/forcomposites factor/;
   my $nsemis = 0;
@@ -2907,10 +2907,11 @@ Project Euler, problem 187, stupid brute force solution, ~3 minutes:
   say $nsemis;
 
 Here is one of the best ways for PE187:  under 20 milliseconds from the
-command line.  This is faster than Mathematica until C<10^13>.
+command line.  Much faster than Pari, and competitive with Mathematica.
 
   use Math::Prime::Util qw/forprimes prime_count/;
-  my $limit = shift || int(10**8)-1;
+  my $limit = shift || int(10**8);
+  $limit--;
   my ($sum, $pc) = (0, 1);
   forprimes {
     $sum += prime_count(int($limit/$_)) + 1 - $pc++;

-- 
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/pkg-perl/packages/libmath-prime-util-perl.git