[libmath-prime-util-perl] 01/13: Fix segmented euler_phi

[Partha P. Mukherjee]

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

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

commit 6390df317eea8be534707656939a4b335545d904
Author: Dana Jacobsen <dana at acm.org>
Date:   Tue Mar 5 16:00:30 2013 -0800

    Fix segmented euler_phi
---
 Changes             |  4 ++++
 MANIFEST            |  1 +
 t/13-primecount.t   |  3 +++
 t/19-moebius.t      |  3 ++-
 util.c              | 10 ++++++----
 xt/totient-range.pl | 31 +++++++++++++++++++++++++++++++
 6 files changed, 47 insertions(+), 5 deletions(-)

diff --git a/Changes b/Changes
index db49369..8b20caf 100644
--- a/Changes
+++ b/Changes
@@ -1,5 +1,9 @@
 Revision history for Perl extension Math::Prime::Util.
 
+0.24  5 March 2013
+
+    - euler_phi on a range wasn't working right with some ranges.
+
 0.23  5 March 2013
 
     - Replace XS Zeta for x > 5 with series from Cephes.  It is 1 eps more
diff --git a/MANIFEST b/MANIFEST
index 8506ef8..8df46f1 100644
--- a/MANIFEST
+++ b/MANIFEST
@@ -86,6 +86,7 @@ t/90-release-perlcritic.t
 t/91-release-pod-syntax.t
 t/92-release-pod-coverage.t
 xt/moebius-mertens.pl
+xt/totient-range.pl
 xt/primality-small.pl
 xt/primality-aks.pl
 xt/factor-holf.pl
diff --git a/t/13-primecount.t b/t/13-primecount.t
index 9448451..09fe825 100644
--- a/t/13-primecount.t
+++ b/t/13-primecount.t
@@ -62,6 +62,9 @@ my %intervals = (
   "24689 to 7973249" => 535368,
   "1e10 +2**16" => 2821,
   "17 to 13"    => 0,
+  "0 to 1"      => 0,
+  "0 to 2"      => 1,
+  "1 to 3"      => 2,
   "3 to 17"     => 6,
   "4 to 17"     => 5,
   "4 to 16"     => 4,
diff --git a/t/19-moebius.t b/t/19-moebius.t
index 9b976cd..5353ed8 100644
--- a/t/19-moebius.t
+++ b/t/19-moebius.t
@@ -161,7 +161,7 @@ plan tests => 0 + 1
                 + 3*scalar(keys %mertens)
                 + 1*scalar(keys %big_mertens)
                 + 2 # Small Phi
-                + 6 + scalar(keys %totients)
+                + 7 + scalar(keys %totients)
                 + scalar(keys %jordan_totients)
                 + 2  # Dedekind psi calculated two ways
                 + 1  # Calculate J5 two different ways
@@ -210,6 +210,7 @@ is_deeply( [euler_phi(0,1)], [0,1],   "euler_phi 0-1" );
 is_deeply( [euler_phi(1,2)], [1,1],   "euler_phi 1-2" );
 is_deeply( [euler_phi(1,3)], [1,1,2], "euler_phi 1-3" );
 is_deeply( [euler_phi(2,3)], [1,2],   "euler_phi 2-3" );
+is_deeply( [euler_phi(10,20)], [4,10,4,12,6,8,8,16,6,18,8], "euler_phi 10-20" );
 
 ###### Jordan Totient
 while (my($k, $tref) = each (%jordan_totients)) {
diff --git a/util.c b/util.c
index 1da3426..368fc24 100644
--- a/util.c
+++ b/util.c
@@ -652,6 +652,8 @@ UV _XS_nth_prime(UV n)
 /* Return an IV array with lo-hi+1 elements.  mu[k-lo] = µ(k) for k = lo .. hi.
  * It is the callers responsibility to call Safefree on the result. */
 #define PGTLO(p,lo)  ((p) >= lo) ? (p) : ((p)*(lo/(p)) + ((lo%(p))?(p):0))
+#define P2GTLO(pinit, p, lo) \
+   ((pinit) >= lo) ? (pinit) : ((p)*(lo/(p)) + ((lo%(p))?(p):0))
 char* _moebius_range(UV lo, UV hi)
 {
   char* mu;
@@ -756,16 +758,16 @@ UV* _totient_range(UV lo, UV hi) {
     croak("Could not get memory for %"UVuf" totients\n", hi);
   for (i = lo; i <= hi; i++)
     totients[i-lo] = i;
-  for (i = PGTLO(2*2, lo); i <= hi; i += 2) totients[i-lo] -= totients[i-lo]/2;
-  for (i = PGTLO(2*3, lo); i <= hi; i += 3) totients[i-lo] -= totients[i-lo]/3;
-  for (i = PGTLO(2*5, lo); i <= hi; i += 5) totients[i-lo] -= totients[i-lo]/5;
   sievehi = hi/2;
+  for (i=P2GTLO(2*2,2,lo); i <= hi; i += 2) totients[i-lo] -= totients[i-lo]/2;
+  for (i=P2GTLO(2*3,3,lo); i <= hi; i += 3) totients[i-lo] -= totients[i-lo]/3;
+  for (i=P2GTLO(2*5,5,lo); i <= hi; i += 5) totients[i-lo] -= totients[i-lo]/5;
   if (get_prime_cache(sievehi, &sieve) < sievehi) {
     release_prime_cache(sieve);
     croak("Could not generate sieve for %"UVuf, sievehi);
   } else {
     START_DO_FOR_EACH_SIEVE_PRIME( sieve, 7, sievehi ) {
-      for (i = PGTLO(2*p, lo); i <= hi; i += p)
+      for (i = P2GTLO(2*p,p,lo); i <= hi; i += p)
         totients[i-lo] -= totients[i-lo]/p;
     } END_DO_FOR_EACH_SIEVE_PRIME
     release_prime_cache(sieve);
diff --git a/xt/totient-range.pl b/xt/totient-range.pl
new file mode 100755
index 0000000..6371a6c
--- /dev/null
+++ b/xt/totient-range.pl
@@ -0,0 +1,31 @@
+#!/usr/bin/env perl
+use strict;
+use warnings;
+$| = 1;  # fast pipes
+
+use Math::Prime::Util qw/euler_phi/;
+use List::Util qw/sum/;
+
+my $limit = shift || 1_000_000;
+
+print "Calculating totients from 1 to $limit...";
+my @phi = map { euler_phi($_) } 1 .. $limit;
+print "...";
+unshift @phi, 0;
+print "...done\n";
+
+while (1) {
+  my $beg = 1 + int(rand($limit));
+  my $end = 1 + int(rand($limit));
+  ($beg,$end) = ($end,$beg) if $beg > $end;
+
+  # Does range return the same values?
+  my @phi_range = @phi[ $beg .. $end ];
+  my @totients = euler_phi($beg,$end);
+
+  my $sum1 = sum(@phi_range);
+  my $sum2 = sum(@totients);
+  warn "\nbeg $beg  end $end  sum $sum1  range sum $sum2\n"
+       unless $sum1 == $sum2;
+  print ".";
+}

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