[libmath-prime-util-perl] 01/08: Totient and Mobius changes, move factor loop out of XS.xs
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:47:21 UTC 2015
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 <dana at 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
More information about the Pkg-perl-cvs-commits
mailing list