[libmath-prime-util-perl] 04/05: Add factoring
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:44:12 UTC 2015
This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.01
in repository libmath-prime-util-perl.
commit b918b92808123e073e61a08c43b0c6a50f0c2677
Author: Dana Jacobsen <dana at acm.org>
Date: Mon Jun 4 05:16:23 2012 -0600
Add factoring
---
MANIFEST | 4 +-
Makefile.PL | 2 +-
TODO | 2 +
XS.xs | 116 ++++++++++++++
factor.c | 416 +++++++++++++++++++++++++++++++++++++++++++++++++
factor.h | 21 +++
lib/Math/Prime/Util.pm | 71 +++++++++
t/50-factoring.t | 45 ++++++
util.c | 7 -
9 files changed, 675 insertions(+), 9 deletions(-)
diff --git a/MANIFEST b/MANIFEST
index a0a6eb4..9e0266f 100644
--- a/MANIFEST
+++ b/MANIFEST
@@ -7,6 +7,8 @@ README
TODO
XS.xs
bitarray.h
+factor.h
+factor.c
sieve.h
sieve.c
util.h
@@ -20,7 +22,7 @@ t/12-nextprime.t
t/13-primecount.t
t/14-nthprime.t
t/30-relations.t
+t/50-factoring.t
t/90-release-perlcritic.t
t/91-release-pod-syntax.t
t/92-release-pod-coverage.t
-
diff --git a/Makefile.PL b/Makefile.PL
index afb943b..88587fd 100644
--- a/Makefile.PL
+++ b/Makefile.PL
@@ -10,7 +10,7 @@ WriteMakefile1(
LIBS => [''], # e.g., '-lm'
DEFINE => '', # e.g., '-DHAVE_SOMETHING'
INC => '', # e.g., '-I/usr/include/other'
- OBJECT => 'sieve.o util.o XS.o',
+ OBJECT => 'factor.o sieve.o util.o XS.o',
PREREQ_PM => {
'Test::More' => '0.45',
diff --git a/TODO b/TODO
index 4ddd74e..180084b 100644
--- a/TODO
+++ b/TODO
@@ -3,4 +3,6 @@
- GMP versions of all routines.
+- prime_count and nth_prime with very large inputs.
+
- factoring (Fermat, Pollard Rho, HOLF, etc.)
diff --git a/XS.xs b/XS.xs
index a2ecd9b..f91b9da 100644
--- a/XS.xs
+++ b/XS.xs
@@ -6,6 +6,7 @@
#include "sieve.h"
#include "util.h"
#include "bitarray.h"
+#include "factor.h"
MODULE = Math::Prime::Util PACKAGE = Math::Prime::Util
@@ -214,3 +215,118 @@ erat_simple_primes(IN UV low, IN UV high)
RETVAL = newRV_noinc( (SV*) av );
OUTPUT:
RETVAL
+
+void
+factor(IN UV n)
+ PREINIT:
+ UV const maxtrials = UV_MAX;
+ UV factors[MPU_MAX_FACTORS+1];
+ int nfactors;
+ int i;
+ PPCODE:
+#if BITS_PER_WORD == 32
+ nfactors = trial_factor(n, factors, maxtrials);
+ if (nfactors < 1)
+ croak("No factors from trial_factor");
+ for (i = 0; i < nfactors; i++) {
+ XPUSHs(sv_2mortal(newSVuv( factors[i] )));
+ }
+#else
+ /* TODO: worry about squfof overflow */
+ if ( n < UVCONST(0xFFFFFFFF) ) {
+ /* small number */
+ nfactors = trial_factor(n, factors, maxtrials);
+ } else {
+ UV sqfactors, f1, f2;
+ /* First factor out small numbers */
+ nfactors = trial_factor(n, factors, 29);
+ /* Use SQUFOF to separate the remainder */
+ n = factors[--nfactors];
+ sqfactors = squfof_factor(n, factors+nfactors, 800000);
+ assert(sqfactors <= 2);
+ if (sqfactors == 1) {
+ n = factors[nfactors];
+ nfactors += trial_factor(n, factors+nfactors, maxtrials);
+ } else {
+ UV n1 = factors[nfactors];
+ n = factors[nfactors+1];
+ nfactors += trial_factor(n1, factors+nfactors, maxtrials);
+ nfactors += trial_factor(n, factors+nfactors, maxtrials);
+ }
+ }
+ if (nfactors < 1) croak("No factors");
+ for (i = 0; i < nfactors; i++) {
+ XPUSHs(sv_2mortal(newSVuv( factors[i] )));
+ }
+#endif
+
+void
+fermat_factor(IN UV n)
+ PREINIT:
+ UV factors[MPU_MAX_FACTORS+1];
+ int nfactors;
+ int i;
+ PPCODE:
+ nfactors = fermat_factor(n, factors);
+ if (nfactors < 1)
+ croak("No factors from fermat_factor");
+ for (i = 0; i < nfactors; i++) {
+ XPUSHs(sv_2mortal(newSVuv( factors[i] )));
+ }
+
+void
+squfof_factor(IN UV n)
+ PREINIT:
+ UV factors[MPU_MAX_FACTORS+1];
+ int nfactors;
+ int i;
+ PPCODE:
+ nfactors = squfof_factor(n, factors, 800000);
+ if (nfactors < 1)
+ croak("No factors from squfof_factor");
+ for (i = 0; i < nfactors; i++) {
+ XPUSHs(sv_2mortal(newSVuv( factors[i] )));
+ }
+
+
+void
+pbrent_factor(IN UV n)
+ PREINIT:
+ UV factors[MPU_MAX_FACTORS+1];
+ int nfactors;
+ int i;
+ PPCODE:
+ nfactors = pbrent_factor(n, factors, 50000000);
+ if (nfactors < 1)
+ croak("No factors from pbrent_factor");
+ for (i = 0; i < nfactors; i++) {
+ XPUSHs(sv_2mortal(newSVuv( factors[i] )));
+ }
+
+void
+prho_factor(IN UV n)
+ PREINIT:
+ UV factors[MPU_MAX_FACTORS+1];
+ int nfactors;
+ int i;
+ PPCODE:
+ nfactors = prho_factor(n, factors, 50000000);
+ if (nfactors < 1)
+ croak("No factors from prho_factor");
+ for (i = 0; i < nfactors; i++) {
+ XPUSHs(sv_2mortal(newSVuv( factors[i] )));
+ }
+
+void
+pminus1_factor(IN UV n)
+ PREINIT:
+ UV factors[MPU_MAX_FACTORS+1];
+ int nfactors;
+ int i;
+ PPCODE:
+ nfactors = pminus1_factor(n, factors, 50000000);
+ if (nfactors < 1)
+ croak("No factors from pminus1_factor");
+ for (i = 0; i < nfactors; i++) {
+ XPUSHs(sv_2mortal(newSVuv( factors[i] )));
+ }
diff --git a/factor.c b/factor.c
new file mode 100644
index 0000000..4957642
--- /dev/null
+++ b/factor.c
@@ -0,0 +1,416 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <assert.h>
+#include <limits.h>
+#include <math.h>
+
+#include "factor.h"
+#include "util.h"
+#include "sieve.h"
+
+
+int trial_factor(UV n, UV *factors, UV maxtrial)
+{
+ UV f, d, m;
+ UV limit;
+ int nfactors = 0;
+
+ if (maxtrial == 0) maxtrial = UV_MAX;
+
+ if ( (n < 2) || (maxtrial < 2) ) {
+ factors[0] = n;
+ return 1;
+ }
+
+ while ( (n & 1) == 0 ) {
+ factors[nfactors++] = 2;
+ n /= 2;
+ }
+
+ for (f = 3; (n > 1) && (f <= 7) && (f <= maxtrial); f += 2) {
+ while ( (n % f) == 0 ) {
+ factors[nfactors++] = f;
+ n /= f;
+ }
+ }
+
+ if ( (n < (7*7)) || (maxtrial < 11) ) {
+ if (n != 1)
+ factors[nfactors++] = n;
+ return nfactors;
+ }
+
+ limit = sqrt((double) n);
+ if (limit > maxtrial)
+ limit = maxtrial;
+
+ /* wheel 30 */
+ f = 11;
+ d = 0;
+ m = 11;
+ while (f <= limit) {
+ if ( (n%f) == 0 ) {
+ UV newlimit;
+ do {
+ factors[nfactors++] = f;
+ n /= f;
+ } while ( (n%f) == 0 );
+ newlimit = sqrt(n);
+ if (newlimit < limit) limit = newlimit;
+ }
+ m = nextwheel30[m]; if (m == 1) d++;
+ f = d*30 + m;
+ }
+ if (n != 1)
+ factors[nfactors++] = n;
+ return nfactors;
+}
+
+static UV gcd_ui(UV x, UV y) {
+ UV t;
+
+ if (y < x) { t = x; x = y; y = t; }
+
+ while (y > 0) {
+ x = x % y;
+ t = x; x = y; y = t;
+ }
+ return x;
+}
+
+/* n^power + a mod m */
+static UV powmod(UV n, UV power, UV add, UV m) {
+ UV t = 1;
+ while (power) {
+ if (power & 1)
+ t = ((t % m) * (n % m)) % m;
+ n = ((n % m) * (n % m)) % m;
+ power >>= 1;
+ }
+ return (t+add) % m;
+}
+
+/* Knuth volume 2, algorithm C.
+ * Very fast for small numbers, grows rapidly.
+ * SQUFOF is better for numbers nearing the 64-bit limit.
+ */
+int fermat_factor(UV n, UV *factors)
+{
+ int nfactors = 0;
+ IV sqn, x, y, r;
+
+ if (n < 2) {
+ factors[0] = n;
+ return 1;
+ }
+
+ while ((n & 1) == 0) {
+ factors[nfactors++] = 2;
+ n /= 2;
+ }
+
+ if (n == 1)
+ return nfactors;
+
+ sqn = sqrt((double) n);
+ x = 2 * sqn + 1;
+ y = 1;
+ r = (sqn*sqn) - n;
+
+ while (r != 0) {
+ r += x;
+ x += 2;
+ do {
+ r -= y;
+ y += 2;
+ } while (r > 0);
+ }
+ r = (x-y)/2;
+ if (r != 1)
+ factors[nfactors++] = r;
+ n /= r;
+ if (n != 1)
+ factors[nfactors++] = n;
+ return nfactors;
+}
+
+/* Pollard / Brent
+ *
+ * Probabilistic. If you give this a prime number, it will loop
+ * until it runs out of rounds.
+ */
+int pbrent_factor(UV n, UV *factors, UV rounds)
+{
+ int nfactors = 0;
+ UV a, f, Xi, Xm, i;
+
+ if (n < 2) {
+ factors[0] = n;
+ return 1;
+ }
+
+ while ((n & 1) == 0) {
+ factors[nfactors++] = 2;
+ n /= 2;
+ }
+
+ if (n == 1)
+ return nfactors;
+
+ Xi = 2;
+ Xm = 2;
+ switch (n%4) {
+ case 0: a = 1; break;
+ case 1: a = 3; break;
+ case 2: a = 5; break;
+ case 3: a = 7; break;
+ default: a = 11; break;
+ }
+
+ for (i = 1; i < rounds; i++) {
+ Xi = powmod(Xi, 2, a, n);
+ f = gcd_ui(Xi - Xm, n);
+ if ( (f != 1) && (f != n) ) {
+ factors[nfactors++] = f;
+ factors[nfactors++] = n/f;
+ return nfactors;
+ }
+ if ( (i & (i-1)) == 0) /* i is a power of 2 */
+ Xm = Xi;
+ }
+ factors[nfactors++] = n;
+ return nfactors;
+}
+
+/* Pollard's Rho
+ *
+ * Probabilistic. If you give this a prime number, it will loop
+ * until it runs out of rounds.
+ */
+int prho_factor(UV n, UV *factors, UV rounds)
+{
+ int nfactors = 0;
+ UV a, f, t, U, V, i;
+
+ if (n < 2) {
+ factors[0] = n;
+ return 1;
+ }
+
+ while ((n & 1) == 0) {
+ factors[nfactors++] = 2;
+ n /= 2;
+ }
+
+ if (n == 1)
+ return nfactors;
+
+ switch (n%4) {
+ case 0: a = 5; break;
+ case 1: a = 7; break;
+ case 2: a = 11; break;
+ case 3: a = 1; break;
+ default: a = 3; break;
+ }
+
+ U = 7;
+ V = 7;
+
+ for (i = 1; i < rounds; i++) {
+ U = powmod(U, 2, a, n);
+ V = powmod(V, 2, a, n);
+ V = powmod(V, 2, a, n);
+
+ f = gcd_ui( (U > V) ? U-V : V-U, n);
+ if ( (f != 1) && (f != n) ) {
+ factors[nfactors++] = f;
+ factors[nfactors++] = n/f;
+ return nfactors;
+ }
+ }
+ factors[nfactors++] = n;
+ return nfactors;
+}
+
+/* Pollard's P-1
+ *
+ * Probabilistic. If you give this a prime number, it will loop
+ * until it runs out of rounds.
+ */
+int pminus1_factor(UV n, UV *factors, UV rounds)
+{
+ int nfactors = 0;
+ UV f, b, i;
+
+ if (n < 2) {
+ factors[0] = n;
+ return 1;
+ }
+
+ while ((n & 1) == 0) {
+ factors[nfactors++] = 2;
+ n /= 2;
+ }
+
+ if (n == 1)
+ return nfactors;
+
+ b = 13;
+
+ for (i = 1; i < rounds; i++) {
+ b = powmod(b+1, i, 0, n);
+ if (b == 0) b = n;
+ b--;
+ f = gcd_ui(b, n);
+ if (f == n) {
+ factors[nfactors++] = n;
+ return nfactors;
+ } else if (f != 1) {
+ factors[nfactors++] = f;
+ factors[nfactors++] = n/f;
+ return nfactors;
+ }
+ }
+ factors[nfactors++] = n;
+ return nfactors;
+}
+
+/* My modification of Ben Buhrow's modification of Bob Silverman's SQUFOF code.
+ * I like Jason P's code a lot -- I should put it in. */
+static long qqueue[100];
+static long qpoint;
+static void enqu(long q, long *iter) {
+ qqueue[qpoint] = q;
+ if (++qpoint >= 100) *iter = -1;
+}
+
+int squfof_factor(UV n, UV *factors, UV rounds)
+{
+ int nfactors = 0;
+ UV temp;
+ long iq,ll,l2,p,pnext,q,qlast,r,s,t,i;
+ long jter, iter;
+ int reloop;
+
+ if ( (n < 2) ) {
+ factors[0] = n;
+ return 1;
+ }
+
+ while ((n & 1) == 0) {
+ factors[nfactors++] = 2;
+ n /= 2;
+ }
+
+ if (n == 1)
+ return nfactors;
+
+ /* TODO: What value of n leads to overflow? */
+
+ qlast = 1;
+ s = sqrt(n);
+
+ p = s;
+ temp = n - (s*s); /* temp = n - floor(sqrt(n))^2 */
+ if (temp == 0) {
+ factors[nfactors++] = s;
+ factors[nfactors++] = s;
+ return nfactors;
+ }
+
+ q = temp; /* q = excess of n over next smaller square */
+ ll = 1 + 2*(long)sqrt((double)(p+p));
+ l2 = ll/2;
+ qpoint = 0;
+
+ /* In the loop below, we need to check if q is a square right before */
+ /* the end of the loop. Is there a faster way? The current way is */
+ /* EXPENSIVE! (many branches and double prec sqrt) */
+
+ for (jter=0; jter < rounds; jter++) {
+ iq = (s + p)/q;
+ pnext = iq*q - p;
+ if (q <= ll) {
+ if ((q & 1) == 0)
+ enqu(q/2,&jter);
+ else if (q <= l2)
+ enqu(q,&jter);
+ if (jter < 0) {
+ factors[nfactors++] = n;
+ return nfactors;
+ }
+ }
+
+ t = qlast + iq*(p - pnext);
+ qlast = q;
+ q = t;
+ p = pnext; /* check for square; even iter */
+ if (jter & 1) continue; /* jter is odd:omit square test */
+ r = (int)sqrt((double)q); /* r = floor(sqrt(q)) */
+ if (q != r*r) continue;
+ if (qpoint == 0) break;
+ reloop = 0;
+ for (i=0; i<qpoint-1; i+=2) { /* treat queue as list for simplicity*/
+ if (r == qqueue[i]) { reloop = 1; break; }
+ if (r == qqueue[i+1]) { reloop = 1; break; }
+ }
+ if (reloop || (r == qqueue[qpoint-1])) continue;
+ break;
+ } /* end of main loop */
+
+ qlast = r;
+ p = p + r*((s - p)/r);
+ q = (n - (p*p)) / qlast; /* q = (n - p*p)/qlast (div is exact)*/
+ for (iter=0; iter<(rounds/16); iter++) { /* unrolled second main loop */
+ iq = (s + p)/q;
+ pnext = iq*q - p;
+ if (p == pnext) break;
+ t = qlast + iq*(p - pnext);
+ qlast = q;
+ q = t;
+ p = pnext;
+ iq = (s + p)/q;
+ pnext = iq*q - p;
+ if (p == pnext) break;
+ t = qlast + iq*(p - pnext);
+ qlast = q;
+ q = t;
+ p = pnext;
+ iq = (s + p)/q;
+ pnext = iq*q - p;
+ if (p == pnext) break;
+ t = qlast + iq*(p - pnext);
+ qlast = q;
+ q = t;
+ p = pnext;
+ iq = (s + p)/q;
+ pnext = iq*q - p;
+ if (p == pnext) break;
+ t = qlast + iq*(p - pnext);
+ qlast = q;
+ q = t;
+ p = pnext;
+ }
+
+ if (iter >= (rounds/20)) {
+ factors[nfactors++] = n;
+ return nfactors;
+ }
+
+ if ((q & 1) == 0) q/=2; /* q was factor or 2*factor */
+ p = n/q;
+
+ if (p < q) {
+ factors[nfactors++] = p;
+ factors[nfactors++] = q;
+ } else {
+ factors[nfactors++] = q;
+ factors[nfactors++] = p;
+ }
+ return nfactors;
+}
+
+
+/* TODO: Add Jason Papadopoulos's racing SQUFOF */
+
diff --git a/factor.h b/factor.h
new file mode 100644
index 0000000..39d2001
--- /dev/null
+++ b/factor.h
@@ -0,0 +1,21 @@
+#ifndef MPU_FACTOR_H
+#define MPU_FACTOR_H
+
+#include "EXTERN.h"
+#include "perl.h"
+
+#define MPU_MAX_FACTORS 64
+
+extern int trial_factor(UV n, UV *factors, UV maxtrial);
+
+extern int fermat_factor(UV n, UV *factors);
+
+extern int squfof_factor(UV n, UV *factors, UV rounds);
+
+extern int pbrent_factor(UV n, UV *factors, UV maxrounds);
+
+extern int prho_factor(UV n, UV *factors, UV maxrounds);
+
+extern int pminus1_factor(UV n, UV *factors, UV maxrounds);
+
+#endif
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index f71ce1f..aa95762 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -18,6 +18,7 @@ our @EXPORT_OK = qw(
next_prime prev_prime
prime_count prime_count_lower prime_count_upper prime_count_approx
nth_prime nth_prime_lower nth_prime_upper nth_prime_approx
+ factor
);
BEGIN {
@@ -321,6 +322,66 @@ If you call another function that uses a MemFree object, the cache will stay
in place because you still have an object.
+
+=head1 FACTORING FUNCTIONS
+
+=head2 factor
+
+ my @factors = factor(3_369_738_766_071_892_021);
+
+Produces the prime factors of a positive number input. They will typically
+but necessarily be in numerical order. The special cases of C<n = 0> and
+C<n = 1> will return C<n>, which guarantees multiplying the factors together
+will always result in the input value, though those are the only cases where
+the returned factors are not prime.
+
+The current algorithm is to use trial division for 32-bit numbers, while for
+larger numbers a small set of trial divisions is performed, followed by a
+single run of SQUFOF, then trial division of the results. This results in
+faster factorization of most large numbers. More sophisticated methods could
+be used.
+
+=head2 fermat_factor
+
+ my @factors = fermat_factor($n);
+
+Produces factors, not necessarily prime, of the positive number input. The
+particular algorithm is Knuth's algorithm C. For small inputs this will be
+very fast, but it slows down quite rapidly as the number of digits increases.
+If there is a factor close to the midpoint (e.g. a semiprime p*q where p and
+q are the same number of digits), then this will be very fast.
+
+=head2 squfof_factor
+
+ my @factors = squfof_factor($n);
+
+Produces factors, not necessarily prime, of the positive number input. It
+is possible the factorization will fail, in which case you will need to use
+another method. Failure in this case means a single factor is returned that
+is not prime.
+
+=head2 prho_factor
+
+=head2 pbrent_factor
+
+=head2 pminus1_factor
+
+Attempts to find a factor using one of the probabilistic algorigthms of
+Pollard Rho, Brent's modification of Pollard Rho, or Pollard's C<p - 1>.
+These are more specialized algorithms usually used for pre-factoring very
+large inputs, or checking very large inputs for naive mistakes. If given
+a prime input, or if just unlucky, these will take a long time to return
+back the single input value. There are cases where these can result in
+finding a factor or very large inputs in remarkably short time, similar to
+how Fermat's method works very well for factors near C<sqrt(n)>. They are
+also amenable to massively parallel searching.
+
+For 64-bit input, these are unlikely to be of too much use. An optimized
+SQUFOF implementation takes under 20 milliseconds to find a factor for any
+62-bit input on modern desktop computers. Lightweight quadratic sieves
+are typically much faster for moderate in the 19+ digit range.
+
+
=head1 LIMITATIONS
The functions C<prime_count> and C<nth_prime> have not yet transitioned to
@@ -330,6 +391,10 @@ called with very large numbers (C<10^11> and up).
I have not completed testing all the functions near the word size limit
(e.g. C<2^32> for 32-bit machines). Please report any problems you find.
+The factoring algorithms are mildly interesting but really limited by not
+being big number aware. Factoring 64-bit numbers is not much of a challenge
+these days.
+
=head1 PERFORMANCE
@@ -374,6 +439,12 @@ Tomás Oliveira e Silva has released the source for a very fast segmented sieve.
The current implementation does not use these ideas, but future versions likely
will.
+The SQUFOF implementation being used is my modifications to Ben Buhrow's
+modifications to Bob Silverman's code. I may experiment with some other
+implementations (Ben Buhrows and Jason Papadopoulos both have published their
+excellent versions in the public domain).
+
+
=head1 COPYRIGHT
diff --git a/t/50-factoring.t b/t/50-factoring.t
new file mode 100644
index 0000000..15da997
--- /dev/null
+++ b/t/50-factoring.t
@@ -0,0 +1,45 @@
+#!/usr/bin/env perl
+use strict;
+use warnings;
+
+use Test::More;
+use Math::Prime::Util qw/factor is_prime/;
+
+my $use64 = Math::Prime::Util::_maxbits > 32;
+my $extra = defined $ENV{RELEASE_TESTING} && $ENV{RELEASE_TESTING};
+
+my @testn = qw/0 1 2 3 4 5 6 7 8 16 57 64 377 9592 30107 78498 664579 5761455
+ 114256942 2214143 999999929 50847534 455052511 2147483647
+ 4118054813/;
+
+my @testn64 = qw/37607912018 346065536839 600851475143
+ 3204941750802 29844570422669
+ 279238341033925 2623557157654233 24739954287740860
+ 3369738766071892021 10023859281455311421
+ 9007199254740991 9007199254740992 9007199254740993
+ /;
+
+
+push @testn, @testn64 if $use64;
+
+push @testn, qw/9999986200004761 99999989237606677 999999866000004473/
+ if $use64 && $extra;
+
+plan tests => 2 * scalar @testn;
+
+foreach my $n (@testn) {
+ my @f = factor($n);
+ my $facstring = join(' * ', at f);
+
+ # Do they multiply to the number?
+ my $product = 1; $product *= $_ for @f;
+ is( $product, $n, "$n = [ $facstring ]" );
+
+ # Are they all prime?
+ my $isprime = 1; $isprime *= is_prime($_) for @f;
+ if ($n < 2) {
+ ok( !$isprime, "All factors [ $facstring ] of $n are not prime" );
+ } else {
+ ok( $isprime, "All factors [ $facstring ] of $n are prime" );
+ }
+};
diff --git a/util.c b/util.c
index a775239..3947ddd 100644
--- a/util.c
+++ b/util.c
@@ -467,13 +467,6 @@ UV nth_prime(UV n)
if (n < NPRIMES_SMALL)
return primes_small[n];
-if (n == 100000) return 1299709;
-if (n == 1000000) return 15485863;
-if (n == 10000000) return 179424673;
-if (n == 100000000) return 2038074743;
-if (n == 1000000000) return 22801763489;
-if (n == 10000000000) return 252097800623;
-
upper_limit = nth_prime_upper(n);
if (upper_limit == 0) {
croak("nth_prime(%lu) would overflow", (unsigned long)n);
--
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