This is an automated email from the git hooks/post-receive script.

ppm-guest pushed a commit to annotated tag v0.22
in repository libmath-prime-util-perl.

commit 49f6ec70e0601889f77626a11da41813f3780111
Author: Dana Jacobsen <d...@acm.org>
Date:   Sun Feb 24 15:27:21 2013 -0800

    Totient and Mobius changes, move factor loop out of XS.xs
---
 Changes                          |   8 ++
 XS.xs                            | 263 ++++++++++++++-------------------------
 examples/test-primality-small.pl |  42 -------
 factor.c                         | 123 +++++++++++++++++-
 factor.h                         |   2 +
 lib/Math/Prime/Util.pm           | 118 ++++++++++--------
 util.c                           |   7 +-
 xt/primality-small.pl            |  44 +++++++
 8 files changed, 329 insertions(+), 278 deletions(-)

diff --git a/Changes b/Changes
index c842455..91c7e2d 100644
--- a/Changes
+++ b/Changes
@@ -1,5 +1,13 @@
 Revision history for Perl extension Math::Prime::Util.
 
+0.22 xx February 2012
+
+    - Move main factor loop out of xs and into factor.c.
+
+    - Totient and Moebius now have complete XS implementations.
+
+    - Ranged totient uses less memory when segmented.
+
 0.21 22 February 2012
 
     - Switch to using Bytes::Random::Secure for random primes.  This is a
diff --git a/XS.xs b/XS.xs
index d79b4f7..a712910 100644
--- a/XS.xs
+++ b/XS.xs
@@ -214,131 +214,6 @@ erat_primes(IN UV low, IN UV high)
     RETVAL
 
 
-void
-_XS_factor(IN UV n)
-  PPCODE:
-    if (n < 4) {                        /* If n is 0-3, we're done. */
-      XPUSHs(sv_2mortal(newSVuv( n )));
-    } else if (n < 10000000) {          /* For small n, just trial division */
-      int i;
-      UV facs[32];  /* maximum number of factors is log2n */
-      UV nfacs = trial_factor(n, facs, 0);
-      for (i = 1; i <= nfacs; i++) {
-        XPUSHs(sv_2mortal(newSVuv( facs[i-1] )));
-      }
-    } else {
-      int const verbose = _XS_get_verbose();
-      UV const tlim_lower = 401;  /* Trial division through this prime */
-      UV const tlim = 409;        /* This means we've checked through here */
-      UV tofac_stack[MPU_MAX_FACTORS+1];
-      UV factored_stack[MPU_MAX_FACTORS+1];
-      int ntofac = 0;
-      int nfactored = 0;
-
-      { /* Trial division, removes all factors below tlim */
-        int i;
-        UV facs[BITS_PER_WORD+1];
-        UV nfacs = trial_factor(n, facs, tlim_lower);
-        for (i = 1; i < nfacs; i++) {
-          XPUSHs(sv_2mortal(newSVuv( facs[i-1] )));
-        }
-        n = facs[nfacs-1];
-      }
-
-      do { /* loop over each remaining factor */
-        /* In theory we can try to minimize work using is_definitely_prime(n)
-         * but in practice it seems slower. */
-        while ( (n >= (tlim*tlim)) && (!_XS_is_prime(n)) ) {
-          int split_success = 0;
-          /* Adjust the number of rounds based on the number size */
-          UV br_rounds = ((n>>29) < 100000) ?  1500 :  1500;
-          UV sq_rounds = 80000; /* 20k 91%, 40k 98%, 80k 99.9%, 120k 99.99% */
-
-          /* About 94% of random inputs are factored with this pbrent call */
-          if (!split_success) {
-            split_success = pbrent_factor(n, tofac_stack+ntofac, br_rounds)-1;
-            if (verbose) { if (split_success) printf("pbrent 1:  %"UVuf" 
%"UVuf"\n", tofac_stack[ntofac], tofac_stack[ntofac+1]); else printf("pbrent 
0\n"); }
-          }
-
-          if (!split_success && n < (UV_MAX>>3)) {
-            /* SQUFOF with these parameters gets 95% of what's left. */
-            split_success = racing_squfof_factor(n, tofac_stack+ntofac, 
sq_rounds)-1;
-            if (verbose) printf("squfof %d\n", split_success);
-          }
-
-          /* Perhaps prho using different parameters will find it */
-          if (!split_success) {
-            split_success = prho_factor(n, tofac_stack+ntofac, 800)-1;
-            if (verbose) printf("prho %d\n", split_success);
-          }
-
-          /* This p-1 gets about 2/3 of what makes it through the above */
-          if (!split_success) {
-            split_success = pminus1_factor(n, tofac_stack+ntofac, 4000, 
40000)-1;
-            if (verbose) printf("pminus1 %d\n", split_success);
-          }
-
-          /* Some rounds of HOLF, good for close to perfect squares */
-          if (!split_success) {
-            split_success = holf_factor(n, tofac_stack+ntofac, 2000)-1;
-            if (verbose) printf("holf %d\n", split_success);
-          }
-
-          /* Less than 0.1% of random inputs make it here */
-          if (!split_success) {
-            split_success = prho_factor(n, tofac_stack+ntofac, 256*1024)-1;
-            if (verbose) printf("long prho %d\n", split_success);
-          }
-
-          if (split_success) {
-            MPUassert( split_success == 1, "split factor returned more than 2 
factors");
-            ntofac++; /* Leave one on the to-be-factored stack */
-            if ((tofac_stack[ntofac] == n) || (tofac_stack[ntofac] == 1))
-              croak("bad factor\n");
-            n = tofac_stack[ntofac];  /* Set n to the other one */
-          } else {
-            /* Factor via trial division.  Nothing should make it here. */
-            UV f = tlim;
-            UV m = tlim % 30;
-            UV limit = (UV) (sqrt(n)+0.1);
-            if (verbose) printf("doing trial on %"UVuf"\n", n);
-            while (f <= limit) {
-              if ( (n%f) == 0 ) {
-                do {
-                  n /= f;
-                  factored_stack[nfactored++] = f;
-                } while ( (n%f) == 0 );
-                limit = (UV) (sqrt(n)+0.1);
-              }
-              f += wheeladvance30[m];
-              m =  nextwheel30[m];
-            }
-            break;  /* We just factored n via trial division.  Exit loop. */
-          }
-        }
-        /* n is now prime (or 1), so add to already-factored stack */
-        if (n != 1)  factored_stack[nfactored++] = n;
-        /* Pop the next number off the to-factor stack */
-        if (ntofac > 0)  n = tofac_stack[ntofac-1];
-      } while (ntofac-- > 0);
-      /* Now push all the factored results in sorted order */
-      {
-        int i, j;
-        for (i = 0; i < nfactored; i++) {
-          int mini = i;
-          for (j = i+1; j < nfactored; j++)
-            if (factored_stack[j] < factored_stack[mini])
-              mini = j;
-          if (mini != i) {
-            UV t = factored_stack[mini];
-            factored_stack[mini] = factored_stack[i];
-            factored_stack[i] = t;
-          }
-          XPUSHs(sv_2mortal(newSVuv( factored_stack[i] )));
-        }
-      }
-    }
-
 #define SIMPLE_FACTOR(func, n, rounds) \
     if (n <= 1) { \
       XPUSHs(sv_2mortal(newSVuv( n ))); \
@@ -359,6 +234,18 @@ _XS_factor(IN UV n)
     }
 
 void
+_XS_factor(IN UV n)
+  PREINIT:
+    UV factors[MPU_MAX_FACTORS+1];
+    int i, nfactors;
+  PPCODE:
+    nfactors = factor(n, factors);
+    EXTEND(SP, nfactors);
+    for (i = 0; i < nfactors; i++) {
+      PUSHs(sv_2mortal(newSVuv( factors[i] )));
+    }
+
+void
 trial_factor(IN UV n, IN UV maxfactor = 0)
   PPCODE:
     SIMPLE_FACTOR(trial_factor, n, maxfactor);
@@ -463,58 +350,88 @@ double
 _XS_RiemannR(double x)
 
 
-SV*
-_XS_totient_range(IN UV lo, IN UV hi)
-  PREINIT:
-    UV* totients;
-    AV* av = newAV();
-    UV i, j;
-  CODE:
-    /* Calculate Euler's totient for all n: lo <= n <= hi. */
-    /* Return as array ref */
-    New(0, totients, hi+1, UV);
-    if (totients == 0)
-      croak("Could not get memory for %"UVuf" totients\n", hi);
-    for (i = 0; i <= hi; i++)
-      totients[i] = i;
-    if (lo <= 0 && hi >= 0) av_push(av,newSVuv(totients[0]));
-    if (lo <= 1 && hi >= 1) av_push(av,newSVuv(totients[1]));
-    if (lo <= 2 && hi >= 2) {
-      totients[2] = 1;
-      av_push(av,newSVuv(totients[2]));
-    }
-    for (j = 2; j <= hi/2; j++)
-      totients[2*j] /= 2;
-    for (i = 3; i <= hi; i++) {
-      if (totients[i] == i) {
-        totients[i] = i-1;
-        for (j = 2*i; j <= hi; j += i)
-          totients[j] = (totients[j]*(i-1))/i;
+void
+_XS_totient(IN UV lo, IN UV hi = 0)
+  PPCODE:
+    if (hi != lo && hi != 0) {
+      /* Totients in a range, returns array */
+      UV* totients;
+      UV i, p;
+
+      if (hi < lo) XSRETURN_EMPTY;
+      if (lo < 2) {
+        if (lo <= 0 && hi >= 0) XPUSHs(sv_2mortal(newSVuv(0)));
+        if (lo <= 1 && hi >= 1) XPUSHs(sv_2mortal(newSVuv(1)));
+        lo = 2;
+      }
+      New(0, totients, hi-lo+1, UV);
+      if (totients == 0)
+        croak("Could not get memory for %"UVuf" totients\n", hi);
+      for (i = lo; i <= hi; i++)
+        totients[i-lo] = i;
+      prime_precalc( hi/2 );
+      for (p = 2; p <= hi/2; p = _XS_next_prime(p)) {
+        i = 2*p;
+        if (i < lo)  i = p*(lo/p) + ( (lo%p) ? p : 0 );
+        for ( ; i <= hi; i += p)
+          totients[i-lo] -= totients[i-lo]/p;
+      }
+      /* Extend the stack to handle how many items we'll return */
+      EXTEND(SP, hi-lo+1);
+      for (i = lo; i <= hi; i++) {
+        UV t = totients[i-lo];
+        if (t == i)
+          t = i-1;
+        PUSHs(sv_2mortal(newSVuv(t)));
+      }
+      Safefree(totients);
+
+    } else {
+      UV facs[64];  /* maximum number of factors is log2n */
+      UV i, nfacs, totient, lastf;
+      UV n = lo;
+      if (n <= 1) XSRETURN_UV(n);
+      nfacs = trial_factor(n, facs, 0);
+      totient = 1;
+      lastf = 0;
+      for (i = 0; i < nfacs; i++) {
+        UV f = facs[i];
+        if (f == lastf) { totient *= f;               }
+        else            { totient *= f-1;  lastf = f; }
       }
-      if (i >= lo)
-        av_push(av,newSVuv(totients[i]));
+      PUSHs(sv_2mortal(newSVuv(totient)));
     }
-    Safefree(totients);
-    RETVAL = newRV_noinc( (SV*) av );
-  OUTPUT:
-    RETVAL
 
-SV*
-_XS_moebius_range(IN UV lo, IN UV hi)
-  PREINIT:
-    IV* mu;
-    AV* av = newAV();
+void
+_XS_moebius(IN UV lo, IN UV hi = 0)
+  PPCODE:
     UV i;
-  CODE:
-    /* Return array ref of Moebius function for all n: lo <= n <= hi. */
-    mu = _moebius_range(lo, hi);
-    MPUassert( mu != 0, "_moebius_range returned 0" );
-    for (i = lo; i <= hi; i++)
-      av_push(av,newSViv(mu[i-lo]));
-    Safefree(mu);
-    RETVAL = newRV_noinc( (SV*) av );
-  OUTPUT:
-    RETVAL
+    if (hi != lo && hi != 0) {   /* mobius in a range */
+      IV* mu = _moebius_range(lo, hi);
+      MPUassert( mu != 0, "_moebius_range returned 0" );
+      EXTEND(SP, hi-lo+1);
+      for (i = lo; i <= hi; i++)
+        PUSHs(sv_2mortal(newSViv(mu[i-lo])));
+      Safefree(mu);
+    } else {
+      UV factors[MPU_MAX_FACTORS+1];
+      UV i, nfactors, lastf;
+      UV n = lo;
+
+      if (n <= 1)
+        XSRETURN_IV(n);
+      if ( (n >= 25) && ( !(n%4) || !(n%9) || !(n%25) ) )
+        XSRETURN_IV(0);
+
+      nfactors = factor(n, factors);
+      lastf = 0;
+      for (i = 0; i < nfactors; i++) {
+        if (factors[i] == lastf)
+          XSRETURN_IV(0);
+        lastf = factors[i];
+      }
+      XSRETURN_IV( (nfactors % 2) ? -1 : 1 );
+    }
 
 IV
 _XS_mertens(IN UV n)
diff --git a/examples/test-primality-small.pl b/examples/test-primality-small.pl
deleted file mode 100755
index a3f145c..0000000
--- a/examples/test-primality-small.pl
+++ /dev/null
@@ -1,42 +0,0 @@
-#!/usr/bin/env perl
-use strict;
-use warnings;
-$| = 1;  # fast pipes
-
-# Make sure the is_prob_prime functionality is working for small inputs.
-# Good for making sure the first few M-R bases are set up correctly.
-
-use Math::Prime::Util qw/is_prob_prime/;
-use Math::Prime::Util::PrimeArray;
-
-my @primes;  tie @primes, 'Math::Prime::Util::PrimeArray';
-
-# Test just primes
-if (0) {
-  foreach my $i (1 .. 10000000) {
-    my $n = shift @primes;
-    die unless is_prob_prime($n);
-    #print "." unless $i % 16384;
-    print "$i $n\n" unless $i % 262144;
-  }
-}
-
-# Test every number up to the 100Mth prime (about 2000M)
-if (1) {
-  die "2 should be prime" unless is_prob_prime(2);
-  shift @primes;
-  my $n = shift @primes;
-  foreach my $i (2 .. 100_000_000) {
-    die "$n should be prime" unless is_prob_prime($n);
-    print "$i $n\n" unless $i % 262144;
-    my $next = shift @primes;
-    my $diff = ($next - $n) >> 1;
-    if ($diff > 1) {
-      foreach my $d (1 .. $diff-1) {
-        my $cn = $n + 2*$d;
-        die "$cn should be composite" if is_prob_prime($cn);
-      }
-    }
-    $n = $next;
-  }
-}
diff --git a/factor.c b/factor.c
index 96c7e64..d4dddaa 100644
--- a/factor.c
+++ b/factor.c
@@ -20,6 +20,117 @@
  * match the native integer type used inside our Perl, so just use those.
  */
 
+/* The main factoring loop */
+/* Puts factors in factors[] and returns the number found. */
+int factor(UV n, UV *factors)
+{
+  int nfactors = 0;           /* Number of factored in factors result */
+
+  int const verbose = _XS_get_verbose();
+  UV const tlim_lower = 401;  /* Trial division through this prime */
+  UV const tlim = 409;        /* This means we've checked through here */
+  UV tofac_stack[MPU_MAX_FACTORS+1];
+  UV fac_stack[MPU_MAX_FACTORS+1];
+  int ntofac = 0;             /* Number of items on tofac_stack */
+  int nfac = 0;               /* Number of items on fac_stack */
+
+  if (n < 10000000)
+    return trial_factor(n, factors, 0);
+
+  /* Trial division for all factors below tlim */
+  nfactors = trial_factor(n, factors, tlim_lower);
+  n = factors[--nfactors];
+
+  /* loop over each remaining factor, until ntofac == 0 */
+  do {
+    while ( (n >= (tlim*tlim)) && (!_XS_is_prime(n)) ) {
+      int split_success = 0;
+      /* Adjust the number of rounds based on the number size */
+      UV br_rounds = ((n>>29) < 100000) ?  1500 :  1500;
+      UV sq_rounds = 80000; /* 20k 91%, 40k 98%, 80k 99.9%, 120k 99.99% */
+
+      /* About 94% of random inputs are factored with this pbrent call */
+      if (!split_success) {
+        split_success = pbrent_factor(n, tofac_stack+ntofac, br_rounds)-1;
+        if (verbose) { if (split_success) printf("pbrent 1:  %"UVuf" 
%"UVuf"\n", tofac_stack[ntofac], tofac_stack[ntofac+1]); else printf("pbrent 
0\n"); }
+      }
+      /* SQUFOF with these parameters gets 95% of what's left. */
+      if (!split_success && n < (UV_MAX>>3)) {
+        split_success = racing_squfof_factor(n,tofac_stack+ntofac, 
sq_rounds)-1;
+        if (verbose) printf("squfof %d\n", split_success);
+      }
+      /* Perhaps prho using different parameters will find it */
+      if (!split_success) {
+        split_success = prho_factor(n, tofac_stack+ntofac, 800)-1;
+        if (verbose) printf("prho %d\n", split_success);
+      }
+      /* This p-1 gets about 2/3 of what makes it through the above */
+      if (!split_success) {
+        split_success = pminus1_factor(n, tofac_stack+ntofac, 4000, 40000)-1;
+        if (verbose) printf("pminus1 %d\n", split_success);
+      }
+      /* Some rounds of HOLF, good for close to perfect squares */
+      if (!split_success) {
+        split_success = holf_factor(n, tofac_stack+ntofac, 2000)-1;
+        if (verbose) printf("holf %d\n", split_success);
+      }
+      /* Less than 0.1% of random inputs make it here */
+      if (!split_success) {
+        split_success = prho_factor(n, tofac_stack+ntofac, 256*1024)-1;
+        if (verbose) printf("long prho %d\n", split_success);
+      }
+
+      if (split_success) {
+        MPUassert( split_success == 1, "split factor returned more than 2 
factors");
+        ntofac++; /* Leave one on the to-be-factored stack */
+        if ((tofac_stack[ntofac] == n) || (tofac_stack[ntofac] == 1))
+          croak("bad factor\n");
+        n = tofac_stack[ntofac];  /* Set n to the other one */
+      } else {
+        /* Factor via trial division.  Nothing should make it here. */
+        UV f = tlim;
+        UV m = tlim % 30;
+        UV limit = (UV) (sqrt(n)+0.1);
+        if (verbose) printf("doing trial on %"UVuf"\n", n);
+        while (f <= limit) {
+          if ( (n%f) == 0 ) {
+            do {
+              n /= f;
+              fac_stack[nfac++] = f;
+            } while ( (n%f) == 0 );
+            limit = (UV) (sqrt(n)+0.1);
+          }
+          f += wheeladvance30[m];
+          m =  nextwheel30[m];
+        }
+        break;  /* We just factored n via trial division.  Exit loop. */
+      }
+    }
+    /* n is now prime (or 1), so add to already-factored stack */
+    if (n != 1)  fac_stack[nfac++] = n;
+    /* Pop the next number off the to-factor stack */
+    if (ntofac > 0)  n = tofac_stack[ntofac-1];
+  } while (ntofac-- > 0);
+
+  /* Sort all the results from fac_stack and put into factors result */
+  {
+    int i, j;
+    for (i = 0; i < nfac; i++) {
+      int mini = i;
+      for (j = i+1; j < nfac; j++)
+        if (fac_stack[j] < fac_stack[mini])
+          mini = j;
+      if (mini != i) {
+        UV t = fac_stack[mini];
+        fac_stack[mini] = fac_stack[i];
+        fac_stack[i] = t;
+      }
+      factors[nfactors++] = fac_stack[i];
+    }
+  }
+  return nfactors;
+}
+
 static const unsigned short primes_small[] =
   {0,2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,
    101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191,
@@ -249,18 +360,18 @@ int _XS_is_prob_prime(UV n)
   }
 #else
 #if 1
-  /* Better bases from http://miller-rabin.appspot.com/, 8 Feb 2013 */
+  /* Better bases from http://miller-rabin.appspot.com/, 23 Feb 2013 */
   if (n < UVCONST(291831)) {
     bases[0] = UVCONST(126401071349994536);
     nbases = 1;
   } else if (n < UVCONST(520924141)) {
-    bases[0] = UVCONST( 15   );
+    bases[0] = 15;
     bases[1] = UVCONST( 750068417525532 );
     nbases = 2;
-  } else if (n < UVCONST(109134866497)) {
-    bases[0] = 2;
-    bases[1] = UVCONST( 45650740   );
-    bases[2] = UVCONST( 3722628058 );
+  } else if (n < UVCONST(154639673381)) {
+    bases[0] = 15;
+    bases[1] = UVCONST(  176006322 );
+    bases[2] = UVCONST( 4221622697 );
     nbases = 3;
   } else if (n < UVCONST(47636622961201)) {
     bases[0] = 2;
diff --git a/factor.h b/factor.h
index d2267fd..b5cf735 100644
--- a/factor.h
+++ b/factor.h
@@ -5,6 +5,8 @@
 
 #define MPU_MAX_FACTORS 64
 
+extern int factor(UV n, UV *factors);
+
 extern int trial_factor(UV n, UV *factors, UV maxtrial);
 
 extern int fermat_factor(UV n, UV *factors, UV rounds);
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index 5778fe4..22589e0 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -55,7 +55,9 @@ sub _import_nobigint {
   undef *nth_prime;     *nth_prime       = \&_XS_nth_prime;
   undef *is_strong_pseudoprime;  *is_strong_pseudoprime = \&_XS_miller_rabin;
   undef *miller_rabin;  *miller_rabin    = \&_XS_miller_rabin;
+  undef *moebius;       *moebius         = \&_XS_moebius;
   undef *mertens;       *mertens         = \&_XS_mertens;
+  undef *euler_phi;     *euler_phi       = \&_XS_totient;
 }
 
 BEGIN {
@@ -1065,37 +1067,52 @@ sub all_factors {
 # A030059, A013929, A030229, A002321, A005117, A013929 all relate.
 sub moebius {
   my($n, $nend) = @_;
-  _validate_positive_integer($n, 1);
+  _validate_positive_integer($n);
 
-  # Moebius over a range.
   if (defined $nend) {
     _validate_positive_integer($nend);
-    return () if $nend < $n;
-    return @{ _XS_moebius_range($n, $nend) } if $nend <= $_XS_MAXVAL;
-    my @mu = map { 1 } 0 .. $nend;
-    foreach my $j (2 .. $nend) {
-      next unless $mu[$j] == 1;
-      for (my $i = $j; $i <= $nend; $i += $j) {
-        $mu[$i] = ($mu[$i] == 1) ? -$j : -$mu[$i];
+    return if $nend < $n;
+  } else {
+    $nend = $n;
+  }
+  return _XS_moebius($n, $nend) if $nend <= $_XS_MAXVAL;
+
+  # Moebius over a range.
+  if ($nend != $n) {
+    my ($lo,$hi) = ($n,$nend);
+    my @mu = map { 1 } $lo .. $hi;
+    $mu[0] = 0 if $lo == 0;
+    my $sqrtn = int(sqrt($hi)+0.5);
+    foreach my $p ( @{ primes($sqrtn) } ) {
+      my $i = $p * $p;
+      $i = $i * int($lo/$i) + (($lo % $i)  ? $i : 0)  if $i < $lo;
+      while ($i <= $hi) {
+        $mu[$i-$lo] = 0;
+        $i += $p * $p;
       }
-    }
-    for (my $j = 2; $j*$j <= $nend; $j++) {
-      next unless $mu[$j] == -$j;
-      for (my $i = $j*$j; $i <= $nend; $i += $j*$j) {
-        $mu[$i] = 0;
+      $i = $p;
+      $i = $i * int($lo/$i) + (($lo % $i)  ? $i : 0)  if $i < $lo;
+      while ($i <= $hi) {
+        $mu[$i-$lo] *= -$p;
+        $i += $p;
       }
     }
-    return map { ($_>0) - ($_<0) } @mu[ $n .. $nend ];
+    foreach my $i ($lo .. $hi) {
+      my $m = $mu[$i-$lo];
+      $m *= -1 if abs($m) != $i;
+      $mu[$i-$lo] = ($m>0) - ($m<0);
+    }
+    return @mu;
   }
 
-  return 1 if $n == 1;
+  return $n if $n <= 1;
   # Quick check for small replicated factors
   return 0 if ($n >= 25) && (!($n % 4) || !($n % 9) || !($n % 25));
-
-  my @factors = ($n <= $_XS_MAXVAL) ? _XS_factor($n) : factor($n);
-  my %all_factors;
+  my @factors = factor($n);
+  my $lastf = 0;
   foreach my $factor (@factors) {
-    return 0 if $all_factors{$factor}++;
+    return 0 if $factor == $lastf;
+    $lastf = $factor;
   }
   return (((scalar @factors) % 2) == 0) ? 1 : -1;
 }
@@ -1139,14 +1156,18 @@ sub euler_phi {
   # I am following SAGE's decision for n <= 0.
   return 0 if defined $n && $n < 0;
   _validate_positive_integer($n);
+  if (defined $nend) {
+    _validate_positive_integer($nend);
+    return if $nend < $n;
+  } else {
+    $nend = $n;
+  }
+  return _XS_totient($n, $nend) if $nend <= $_XS_MAXVAL;
 
   # Totient over a range.  Could be improved, as this can use huge memory.
-  if (defined $nend && $nend != $n) {
-    _validate_positive_integer($nend);
+  if ($nend != $n) {
     return () if $nend < $n;
-    # Use XS code if at all possible.  Better memory use.
-    return @{ _XS_totient_range($n, $nend) } if $nend <= $_XS_MAXVAL;
-    my @totients = map { $_ } 0 .. $nend;
+    my @totients = (0 .. $nend);
     foreach my $i (2 .. $nend) {
       next unless $totients[$i] == $i;
       $totients[$i] = $i-1;
@@ -1158,39 +1179,29 @@ sub euler_phi {
     return @totients;
   }
 
-  return (0,1)[$n] if $n <= 1;
-
-  if ($n <= $_XS_MAXVAL) {
-    my $last_f = 0;
-    my $totient = 1;
-    foreach my $f ( _XS_factor($n) ) {
-      if ($f == $last_f) {  $totient *= $f;                   }
-      else               {  $totient *= $f-1;  $last_f = $f;  }
-    }
-    return $totient;
-  }
-
+  return $n if $n <= 1;
   my %factor_mult;
   my @factors = grep { !$factor_mult{$_}++ } factor($n);
-  my $totient = 1;
 
   if (ref($n) ne 'Math::BigInt') {
+    my $totient = 1;
     foreach my $factor (@factors) {
       $totient *= ($factor - 1);
       $totient *= $factor for (2 .. $factor_mult{$factor});
     }
-  } else {
-    # Some real wackiness to solve issues with Math::BigInt::GMP (not seen with
-    # Pari or Calc).  Results of the multiply will go negative if we don't do
-    # this.  To see if you hit the standalone bug:
-    #      perl -E 'my $a = 2931542417; use bigint lib=>'GMP'; my $n = 
49754396241690624; my $x = $n*$a; say $x;'
-    # This may be related to RT 71548 of Math::BigInt::GMP.
-    $totient = $n->copy->bone;
-    foreach my $factor (@factors) {
-      my $f = $n->copy->bzero->badd("$factor");
-      $totient->bmul($f->copy->bsub(1));
-      $totient->bmul($f)  for (2 .. $factor_mult{$factor});
-    }
+    return $totient;
+  }
+
+  # Some real wackiness to solve issues with Math::BigInt::GMP (not seen with
+  # Pari or Calc).  Results of the multiply will go negative if we don't do
+  # this.  To see if you hit the standalone bug:
+  #      perl -E 'my $a = 2931542417; use bigint lib=>'GMP'; my $n = 
49754396241690624; my $x = $n*$a; say $x;'
+  # This may be related to RT 71548 of Math::BigInt::GMP.
+  my $totient = $n->copy->bone;
+  foreach my $factor (@factors) {
+    my $f = $n->copy->bzero->badd("$factor");
+    $totient->bmul($f->copy->bsub(1));
+    $totient->bmul($f)  for (2 .. $factor_mult{$factor});
   }
   return $totient;
 }
@@ -2401,10 +2412,11 @@ primality tests.
   $sum += moebius($_) for (1..200); say "Mertens(200) = $sum";
 
 Returns the Möbius function (also called the Moebius, Mobius, or MoebiusMu
-function) for a positive non-zero integer input.  This function is 1 if
+function) for a non-negative integer input.  This function is 1 if
 C<n = 1>, 0 if C<n> is not square free (i.e. C<n> has a repeated factor),
 and C<-1^t> if C<n> is a product of C<t> distinct primes.  This is an
-important function in prime number theory.
+important function in prime number theory.  Like SAGE, we define
+C<moebius(0) = 0> for convenience.
 
 If called with two arguments, they define a range C<low> to C<high>, and the
 function returns an array with the value of the Möbius function for every n
@@ -3113,7 +3125,7 @@ Here's 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:
+Project Euler, problem 187, stupid brute force solution, ~3 minutes:
 
   use Math::Prime::Util qw/factor -nobigint/;
   my $nsemis = 0;
diff --git a/util.c b/util.c
index 4bf6e4a..37ed751 100644
--- a/util.c
+++ b/util.c
@@ -635,17 +635,16 @@ UV _XS_nth_prime(UV n)
 IV* _moebius_range(UV lo, UV hi)
 {
   IV* mu;
-  UV i, p, sqrtn, range;
+  UV i, p, sqrtn;
 
   /* This implementation follows that of Deléglise & Rivat (1996), which is
    * a segmented version of Lioen & van de Lune (1994).
    */
-  range = hi-lo+1;
   sqrtn = (UV) (sqrt(hi) + 0.5);
 
-  New(0, mu, range, IV);
+  New(0, mu, hi-lo+1, IV);
   if (mu == 0)
-    croak("Could not get memory for %"UVuf" moebius results\n", range);
+    croak("Could not get memory for %"UVuf" moebius results\n", hi-lo+1);
   for (i = lo; i <= hi; i++)
     mu[i-lo] = 1;
   if (lo == 0)  mu[0] = 0;
diff --git a/xt/primality-small.pl b/xt/primality-small.pl
new file mode 100755
index 0000000..dbed6ff
--- /dev/null
+++ b/xt/primality-small.pl
@@ -0,0 +1,44 @@
+#!/usr/bin/env perl
+use strict;
+use warnings;
+$| = 1;  # fast pipes
+
+# Make sure the is_prob_prime functionality is working for small inputs.
+# Good for making sure the first few M-R bases are set up correctly.
+my $limit = 600_000_000;
+
+use Math::Prime::Util qw/is_prob_prime/;
+# Use another code base for comparison.
+# Math::Prime::FastSieve is very fast -- far faster than Math::Primality
+use Math::Prime::FastSieve;
+my $sieve = Math::Prime::FastSieve::Sieve->new($limit + 10_000);
+
+if (0) {  # just primes using Math::Prime::FastSieve
+  my $n = 2;
+  my $i = 1;
+  while ($n < $limit) {
+    die "$n" unless is_prob_prime($n);
+    $n = $sieve->nearest_ge( $n+1 );
+    print "$i $n\n" unless $i++ % 16384;
+  }
+}
+
+# Test every number up to the 100Mth prime (about 2000M)
+if (1) {
+  my $n = 2;
+  my $i = 1;
+  while ($n <= $limit) {
+    die "$n should be prime" unless is_prob_prime($n);
+    print "$i $n\n" unless $i++ % 262144;
+    my $next = $sieve->nearest_ge( $n+1 );
+    my $diff = ($next - $n) >> 1;
+    if ($diff > 1) {
+      foreach my $d (1 .. $diff-1) {
+        my $cn = $n + 2*$d;
+        die "$cn should be composite" if is_prob_prime($cn);
+      }
+    }
+    $n = $next;
+  }
+  print "Success to $limit!\n";
+}

-- 
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