[libmath-prime-util-perl] 22/23: Fix 32-bit issue with lehmer

Partha P. Mukherjee ppm-guest at moszumanska.debian.org
Thu May 21 18:45:58 UTC 2015 

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

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

commit 44aef7dcc448dd1342d90a32b542fed53655c2a0
Author: Dana Jacobsen <dana at acm.org>
Date:   Thu Nov 29 17:16:11 2012 -0800

    Fix 32-bit issue with lehmer
---
 lehmer.c               | 45 +++++++++++++++++++++++++++------------------
 lib/Math/Prime/Util.pm | 15 +++++++--------
 t/18-functions.t       |  2 +-
 3 files changed, 35 insertions(+), 27 deletions(-)

diff --git a/lehmer.c b/lehmer.c
index 6bb5d59..aef6b37 100644
--- a/lehmer.c
+++ b/lehmer.c
@@ -190,21 +190,21 @@ static UV bs_prime_count(UV n, UV const* const primes, UV lastprime)
  * convenient and a little Perl script will spit this code out for whatever
  * limit we select.  It gets unwieldy with large a values.
  */
-static IV mapes(IV x, UV a)
+static UV mapes(UV x, UV a)
 {
   IV val;
   if (a == 0)  return x;
   if (a == 1)  return x-x/2;
   val = x-x/2-x/3+x/6;
-  if (a >= 3) val += -x/5+x/10+x/15-x/30;
-  if (a >= 4) val += -x/7+x/14+x/21-x/42+x/35-x/70-x/105+x/210;
-  if (a >= 5) val += -x/11+x/22+x/33-x/66+x/55-x/110-x/165+x/330+x/77-x/154-x/231+x/462-x/385+x/770+x/1155-x/2310;
-  if (a >= 6) val += -x/13+x/26+x/39-x/78+x/65-x/130-x/195+x/390+x/91-x/182-x/273+x/546-x/455+x/910+x/1365-x/2730+x/143-x/286-x/429+x/858-x/715+x/1430+x/2145-x/4290-x/1001+x/2002+x/3003-x/6006+x/5005-x/10010-x/15015+x/30030;
-  if (a >= 7) val += -x/17+x/34+x/51-x/102+x/85-x/170-x/255+x/510+x/119-x/238-x/357+x/714-x/595+x/1190+x/1785-x/3570+x/187-x/374-x/561+x/1122-x/935+x/1870+x/2805-x/5610-x/1309+x/2618+x/3927-x/7854+x/6545-x/13090-x/19635+x/39270+x/221-x/442-x/663+x/1326-x/1105+x/2210+x/3315-x/6630-x/1547+x/3094+x/4641-x/9282+x/7735-x/15470-x/23205+x/46410-x/2431+x/4862+x/7293-x/14586+x/12155-x/24310-x/36465+x/72930+x/17017-x/34034-x/51051+x/102102-x/85085+x/170170+x/255255-x/510510;
-  return val;
+  if (a >= 3) val += 0-x/5+x/10+x/15-x/30;
+  if (a >= 4) val += 0-x/7+x/14+x/21-x/42+x/35-x/70-x/105+x/210;
+  if (a >= 5) val += 0-x/11+x/22+x/33-x/66+x/55-x/110-x/165+x/330+x/77-x/154-x/231+x/462-x/385+x/770+x/1155-x/2310;
+  if (a >= 6) val += 0-x/13+x/26+x/39-x/78+x/65-x/130-x/195+x/390+x/91-x/182-x/273+x/546-x/455+x/910+x/1365-x/2730+x/143-x/286-x/429+x/858-x/715+x/1430+x/2145-x/4290-x/1001+x/2002+x/3003-x/6006+x/5005-x/10010-x/15015+x/30030;
+  if (a >= 7) val += 0-x/17+x/34+x/51-x/102+x/85-x/170-x/255+x/510+x/119-x/238-x/357+x/714-x/595+x/1190+x/1785-x/3570+x/187-x/374-x/561+x/1122-x/935+x/1870+x/2805-x/5610-x/1309+x/2618+x/3927-x/7854+x/6545-x/13090-x/19635+x/39270+x/221-x/442-x/663+x/1326-x/1105+x/2210+x/3315-x/6630-x/1547+x/3094+x/4641-x/9282+x/7735-x/15470-x/23205+x/46410-x/2431+x/4862+x/7293-x/14586+x/12155-x/24310-x/36465+x/72930+x/17017-x/34034-x/51051+x/102102-x/85085+x/170170+x/255255-x/510510;
+  return (UV) val;
 }
 
-static IV mapes7(IV x) {    /* A tiny bit faster setup for a=7 */
+static UV mapes7(UV x) {    /* A tiny bit faster setup for a=7 */
   IV val = x-x/2-x/3-x/5+x/6-x/7+x/10-x/11-x/13+x/14+x/15-x/17+x/21+x/22+x/26
           -x/30+x/33+x/34+x/35+x/39-x/42+x/51+x/55+x/65-x/66-x/70+x/77-x/78
           +x/85+x/91-x/102-x/105-x/110+x/119-x/130+x/143-x/154-x/165-x/170
@@ -212,17 +212,17 @@ static IV mapes7(IV x) {    /* A tiny bit faster setup for a=7 */
           -x/357-x/374-x/385+x/390-x/429-x/442-x/455+x/462+x/510+x/546-x/561
           -x/595-x/663+x/714;
   if (x >= 715) {
-    val += -x/715+x/770+x/858+x/910-x/935-x/1001-x/1105+x/1122+x/1155+x/1190
-           -x/1309+x/1326+x/1365+x/1430-x/1547+x/1785+x/1870+x/2002+x/2145
-           +x/2210-x/2310-x/2431+x/2618-x/2730+x/2805+x/3003+x/3094+x/3315
-           -x/3570+x/3927-x/4290+x/4641+x/4862+x/5005-x/5610-x/6006+x/6545
-           -x/6630+x/7293+x/7735-x/7854;
+    val += 0-x/715+x/770+x/858+x/910-x/935-x/1001-x/1105+x/1122+x/1155+x/1190
+            -x/1309+x/1326+x/1365+x/1430-x/1547+x/1785+x/1870+x/2002+x/2145
+            +x/2210-x/2310-x/2431+x/2618-x/2730+x/2805+x/3003+x/3094+x/3315
+            -x/3570+x/3927-x/4290+x/4641+x/4862+x/5005-x/5610-x/6006+x/6545
+            -x/6630+x/7293+x/7735-x/7854;
     if (x >= 9282)
-      val += -x/9282-x/10010+x/12155-x/13090-x/14586-x/15015-x/15470+x/17017
-             -x/19635-x/23205-x/24310+x/30030-x/34034-x/36465+x/39270+x/46410
-             -x/51051+x/72930-x/85085+x/102102+x/170170+x/255255-x/510510;
+      val += 0-x/9282-x/10010+x/12155-x/13090-x/14586-x/15015-x/15470+x/17017
+              -x/19635-x/23205-x/24310+x/30030-x/34034-x/36465+x/39270+x/46410
+              -x/51051+x/72930-x/85085+x/102102+x/170170+x/255255-x/510510;
   }
-  return val;
+  return (UV) val;
 }
 
 /******************************************************************************/
@@ -393,7 +393,8 @@ static UV phi(UV x, UV a)
 {
   heap_t h1, h2;
   UV val, smallv;
-  IV count, sum = 0;
+  UV sum = 0;
+  IV count;
   const UV* primes;
 
   if (a == 1)  return ((x+1)/2);
@@ -521,6 +522,14 @@ UV _XS_lehmer_pi(UV n)
   if (n < SIEVE_LIMIT)
     return _XS_prime_count(2, n);
 
+  /* Protect against overflow.  2^32-1 and 2^64-1 are both divisible by 3. */
+  if (n == UV_MAX) {
+    if ( (n%3) == 0 || (n%5) == 0 || (n%7) == 0 || (n%31) == 0 )
+      n--;
+    else
+      return _XS_prime_count(2,n);
+  }
+
   if (verbose > 0) printf("lehmer %lu stage 1: calculate a,b,c \n", n);
   TIMING_START;
   z = (UV) sqrt((double)n+0.5);
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index d701af4..0293972 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -1448,17 +1448,16 @@ sub LogarithmicIntegral {
   return 0              if $n == 0;
   return $_Neg_Infinity if $n == 1;
   return $_Infinity     if $n == $_Infinity;
+  if ($n == 2) {
+    return (defined $bignum::VERSION || ref($n) eq 'Math::BigFloat')
+           ? Math::BigFloat->new('1.045163780117492784844588889194613136522615578151201575832909144075013205210359530172717405626383356306')
+           : 1.045163780117492784844588889194613136522615578151;
+  }
+
   croak("Invalid input to LogarithmicIntegral:  x must be >= 0") if $n <= 0;
 
-  if ( defined $bignum::VERSION || ref($n) eq 'Math::BigFloat' ) {
-    return Math::BigFloat->binf('-') if $n == 1;
-    return Math::BigFloat->new('1.045163780117492784844588889194613136522615578151201575832909144075013205210359530172717405626383356306') if $n == 2;
-  } else {
-    return $_Neg_Infinity if $n == 1;
-    return 1.045163780117492784844588889194613136522615578151 if $n == 2;
-  }
   return Math::Prime::Util::PP::LogarithmicIntegral($n)
-    if defined $bignum::VERSION || ref($n) eq 'Math::BigFloat' || !$_Config{'xs'};
+   if defined $bignum::VERSION || ref($n) eq 'Math::BigFloat' || !$_Config{'xs'};
   return _XS_LogarithmicIntegral($n);
 }
 
diff --git a/t/18-functions.t b/t/18-functions.t
index 7324f6b..5cbb4f2 100644
--- a/t/18-functions.t
+++ b/t/18-functions.t
@@ -24,7 +24,7 @@ cmp_ok( ExponentialIntegral($infinity), '>=', $infinity, "Ei(inf) is infinity");
 
 cmp_ok( LogarithmicIntegral(0),         '==', 0,         "li(0) is 0");
 cmp_ok( LogarithmicIntegral(1),         '<=',-$infinity, "li(1) is -infinity");
-cmp_ok( LogarithmicIntegral($infinity)),'>=', $infinity, "li(inf) is infinity");
+cmp_ok( LogarithmicIntegral($infinity), '>=', $infinity, "li(inf) is infinity");
 
 # Example used in Math::Cephes
 cmp_closeto( ExponentialIntegral(2.2), 5.732614700, 1e-06, "Ei(2.2)");

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