[libmath-prime-util-perl] 110/181: Add gcd function. Fix forcomposites. Fix test for GMP function definition.
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:51:12 UTC 2015
This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.36
in repository libmath-prime-util-perl.
commit d3f0a0550cb3affe106cbcc355c27174809ba739
Author: Dana Jacobsen <dana at acm.org>
Date: Fri Jan 3 19:13:31 2014 -0800
Add gcd function. Fix forcomposites. Fix test for GMP function definition.
---
Changes | 6 ++++-
XS.xs | 66 +++++++++++++++++++++++++++++++++++-----------
lib/Math/Prime/Util.pm | 67 ++++++++++++++++++++++++++---------------------
lib/Math/Prime/Util/PP.pm | 5 ++++
t/19-moebius.t | 33 ++++++++++++++++++++++-
t/32-iterators.t | 13 +++++++--
t/81-bignum.t | 6 +++++
7 files changed, 147 insertions(+), 49 deletions(-)
diff --git a/Changes b/Changes
index 582abf3..4712ef8 100644
--- a/Changes
+++ b/Changes
@@ -17,7 +17,8 @@ Revision history for Perl module Math::Prime::Util
- forcomposites like forprimes, but on composites. See Pari 2.6.x.
- fordivisors calls a block for each divisor
- kronecker Kronecker symbol (extension of Jacobi symbol)
- - znprimroot Primative root modulo n
+ - znprimroot Primitive root modulo n
+ - gcd Greatest common divisor
[FUNCTIONALITY AND PERFORMANCE]
@@ -59,6 +60,9 @@ Revision history for Perl module Math::Prime::Util
- znorder uses Carmichael Lambda instead of Euler Phi. Faster.
+ - While Math::BigInt has the bgcd function, it is slow for native numbers,
+ even with the Pari or GMP back ends. The gcd here is 20-100x faster.
+
0.35 2013-12-08
diff --git a/XS.xs b/XS.xs
index 6cc7cb7..bc31566 100644
--- a/XS.xs
+++ b/XS.xs
@@ -16,6 +16,7 @@
#include "ptypes.h"
#include "cache.h"
#include "sieve.h"
+#define FUNC_gcd_ui
#include "util.h"
#include "primality.h"
#include "factor.h"
@@ -160,9 +161,14 @@ static int _vcallsubn(pTHX_ I32 flags, const char* gmp_name, const char* name, i
/* If given a GMP function, and GMP enabled, and function exists, use it. */
int use_gmp = gmp_name != 0 && _XS_get_callgmp();
if (use_gmp) {
+ CV* cv;
strncat(fullname, gmp_name, 60);
- if (get_cv(fullname, flags) == 0)
+ cv = get_cv(fullname, flags);
+ /* This isn't covering every case for arbitrary functions */
+ if (cv == 0 || (!CvROOT(cv) && !CvXSUB(cv))) {
use_gmp = 0;
+ fullname[19] = '\0';
+ }
}
if (!use_gmp)
strncat(fullname, name, 60);
@@ -220,8 +226,9 @@ prime_memfree()
_XS_get_verbose = 2
_XS_get_callgmp = 3
_get_prime_cache_size = 4
- PPCODE:
+ PREINIT:
UV ret;
+ PPCODE:
switch (ix) {
case 0: prime_memfree(); goto return_nothing;
case 1: _prime_memfreeall(); goto return_nothing;
@@ -401,14 +408,17 @@ trial_factor(IN UV n, ...)
void
is_strong_pseudoprime(IN SV* svn, ...)
PREINIT:
- int status;
+ int c, status = 1;
PPCODE:
if (items < 2)
croak("No bases given to miller_rabin");
- status = _validate_int(aTHX_ svn, 0);
+ /* Check all arguments */
+ for (c = 0; c < items && status == 1; c++)
+ if (_validate_int(aTHX_ ST(c), 0) != 1)
+ status = 0;
if (status == 1) {
UV n = my_svuv(svn);
- int b, c, ret = 1;
+ int b, ret = 1;
if (n < 4) { ret = (n >= 2); } /* 0,1 composite; 2,3 prime */
else if ((n % 2) == 0) { ret = 0; } /* evens composite */
else {
@@ -425,6 +435,27 @@ is_strong_pseudoprime(IN SV* svn, ...)
return; /* skip implicit PUTBACK */
void
+gcd(...)
+ PROTOTYPE: @
+ PREINIT:
+ int i, status = 1;
+ UV ret = 0;
+ PPCODE:
+ /* For each arg, while valid input, validate+gcd. Shortcut stop if 1. */
+ for (i = 0; i < items && status == 1 && ret != 1; i++) {
+ if (_validate_int(aTHX_ ST(i), 1) != 1) {
+ status = 0;
+ } else {
+ UV n = my_svuv(ST(i));
+ ret = (ret == 0) ? n : gcd_ui(ret, n);
+ }
+ }
+ if (status == 1)
+ XSRETURN_UV(ret);
+ _vcallsub_with_gmp("gcd");
+ return; /* skip implicit PUTBACK */
+
+void
_XS_lucas_sequence(IN UV n, IN IV P, IN IV Q, IN UV k)
PREINIT:
UV U, V, Qk;
@@ -578,10 +609,13 @@ factor(IN SV* svn)
void
divisor_sum(IN SV* svn, ...)
+ PREINIT:
+ SV* svk;
+ int nstatus, kstatus;
PPCODE:
- SV* svk = (items > 1) ? ST(1) : 0;
- int nstatus = _validate_int(aTHX_ svn, 0);
- int kstatus = (items == 1 || (SvIOK(svk) && SvIV(svk))) ? 1 : 0;
+ svk = (items > 1) ? ST(1) : 0;
+ nstatus = _validate_int(aTHX_ svn, 0);
+ kstatus = (items == 1 || (SvIOK(svk) && SvIV(svk))) ? 1 : 0;
if (nstatus == 1 && kstatus == 1) {
UV n = my_svuv(svn);
UV k = (items > 1) ? my_svuv(svk) : 1;
@@ -591,30 +625,32 @@ divisor_sum(IN SV* svn, ...)
_vcallsub("_generic_divisor_sum");
return; /* skip implicit PUTBACK */
-
void
znorder(IN SV* sva, IN SV* svn)
ALIAS:
legendre_phi = 1
+ PREINIT:
+ int astatus, nstatus;
PPCODE:
- int astatus = _validate_int(aTHX_ sva, 0);
- int nstatus = _validate_int(aTHX_ svn, 0);
+ astatus = _validate_int(aTHX_ sva, 0);
+ nstatus = _validate_int(aTHX_ svn, 0);
if (astatus == 1 && nstatus == 1) {
UV a = my_svuv(sva);
UV n = my_svuv(svn);
UV ret;
switch (ix) {
- case 0: ret = znorder(a,n);
+ case 0: ret = znorder(a, n);
if (ret == 0) XSRETURN_UNDEF; /* not defined */
break;
- default: ret = legendre_phi(a, n);
- break;
+ case 1:
+ default: ret = legendre_phi(a, n); break;
}
XSRETURN_UV(ret);
}
switch (ix) {
case 0: _vcallsub("_generic_znorder"); break;
/* TODO: Fixup public PP legendre_phi */
+ case 1:
default: _vcallsub("PP::_legendre_phi"); break;
}
return; /* skip implicit PUTBACK */
@@ -946,7 +982,7 @@ forcomposites (SV* block, IN SV* svbeg, IN SV* svend = 0)
START_DO_FOR_EACH_SIEVE_PRIME( segment, seg_low - seg_base, seg_high - seg_base ) {
cbeg = prevprime+1; if (cbeg < beg) cbeg = beg;
prevprime = seg_base + p;
- cend = prevprime-1; if (cend > end) cbeg = end;
+ cend = prevprime-1; if (cend > end) cend = end;
for (c = cbeg; c <= cend; c++) {
sv_setuv(svarg, c); MULTICALL;
}
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index 41a706b..26453f7 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -37,7 +37,7 @@ our @EXPORT_OK =
random_proven_prime random_proven_prime_with_cert
random_maurer_prime random_maurer_prime_with_cert
primorial pn_primorial consecutive_integer_lcm
- factor factor_exp all_factors divisors
+ gcd factor factor_exp all_factors divisors
moebius mertens euler_phi jordan_totient exp_mangoldt liouville
partitions
chebyshev_theta chebyshev_psi
@@ -118,12 +118,13 @@ BEGIN {
*znorder = \&Math::Prime::Util::_generic_znorder;
*znprimroot = \&Math::Prime::Util::_generic_znprimroot;
*legendre_phi = \&Math::Prime::Util::PP::_legendre_phi;
+ *gcd = \&Math::Prime::Util::PP::gcd;
*factor = \&Math::Prime::Util::_generic_factor;
*factor_exp = \&Math::Prime::Util::_generic_factor_exp;
*divisors = \&Math::Prime::Util::_generic_divisors;
*forprimes = sub (&$;$) { _generic_forprimes(@_); }; ## no critic qw(ProhibitSubroutinePrototypes)
+ *forcomposites = sub (&$;$) { _generic_forcomposites(@_); }; ## no critic qw(ProhibitSubroutinePrototypes)
*fordivisors = sub (&$) { _generic_fordivisors(@_); }; ## no critic qw(ProhibitSubroutinePrototypes)
- *forcomposites = sub (&$) { _generic_forcomposites(@_); }; ## no critic qw(ProhibitSubroutinePrototypes)
*_prime_memfreeall = \&Math::Prime::Util::PP::_prime_memfreeall;
*prime_memfree = \&Math::Prime::Util::PP::prime_memfree;
@@ -2552,8 +2553,16 @@ block for each composite in the inclusive range. Starting at 2, the
composites are the non-primes (C<0> and C<1> are neither prime nor composite).
+=head2 fordivisors
+
+ fordivisors { $prod *= $_ } $n;
+
+Given a block and a non-negative number C<n>, the block is called with
+C<$_> set to each divisor in sorted order. Also see L</divisor_sum>.
+
=head2 prime_iterator
+
my $it = prime_iterator;
$sum += $it->() for 1..100000;
@@ -3170,13 +3179,10 @@ The following conditions must hold:
- C<< n >= 2 >>
-=head2 fordivisors
-
- fordivisors { $prod *= $_ } $n;
-
-Given a block and a non-negative number C<n>, the block is called with
-C<$_> set to each divisor in sorted order. Also see L</divisor_sum>.
+=head2 gcd
+Given a list of integers, returns the greatest common divisor. This is
+often used to test for coprimality.
=head2 moebius
@@ -4027,17 +4033,17 @@ is 40 by default). Accuracy without MPFR should be 35 digits.
Print strong pseudoprimes to base 17 up to 10M:
- # Similar to A001262's isStrongPsp function, but over 4x faster
- perl -MMath::Prime::Util=:all -E 'my $n=3; while($n <= 10000000) { print "$n " if is_strong_pseudoprime($n,$base) && !is_prime($n); $n+=2; } BEGIN {$|=1; $base=17}'
-
-or, slightly faster, use forprimes and loop over the odds between primes:
-
- perl -MMath::Prime::Util=:all -E '$|=1; $base=17; my $prev = 1; forprimes { $prev += 2; while ($prev < $_) { print "$prev " if is_strong_pseudoprime($prev,$base); $prev += 2; } } 3,10000000'
+ # Similar to A001262's isStrongPsp function, but much faster
+ perl -MMath::Prime::Util=:all -E 'forcomposites { say if is_strong_pseudoprime($_,17) } 10000000;'
Print some primes above 64-bit range:
perl -MMath::Prime::Util=:all -Mbigint -E 'my $start=100000000000000000000; say join "\n", @{primes($start,$start+1000)}'
- # Similar code using Pari:
+
+ # Another way
+ perl -MMath::Prime::Util=:all -E 'forprimes { say } "100000000000000000039", "100000000000000000993"'
+
+ # Similar using Math::Pari:
# perl -MMath::Pari=:int,PARI,nextprime -E 'my $start = PARI "100000000000000000000"; my $end = $start+1000; my $p=nextprime($start); while ($p <= $end) { say $p; $p = nextprime($p+1); }'
Examining the η3(x) function of Planat and Solé (2011):
@@ -4072,7 +4078,7 @@ Project Euler, problem 3 (Largest prime factor):
use Math::Prime::Util qw/factor/;
use bigint; # Only necessary for 32-bit machines.
- say "", (factor(600851475143))[-1]
+ say 0+(factor(600851475143))[-1]
Project Euler, problem 7 (10001st prime):
@@ -4081,9 +4087,9 @@ Project Euler, problem 7 (10001st prime):
Project Euler, problem 10 (summation of primes):
- use Math::Prime::Util qw/primes/;
+ use Math::Prime::Util qw/forprimes/;
my $sum = 0;
- $sum += $_ for @{primes(2_000_000)};
+ forprimes { $sum += $_ } 2_000_000;
say $sum;
Project Euler, problem 21 (Amicable numbers):
@@ -4103,14 +4109,13 @@ Project Euler, problem 41 (Pandigital prime), brute force command line:
Project Euler, problem 47 (Distinct primes factors):
- use Math::Prime::Util qw/pn_primorial factor/;
- use List::MoreUtils qw/distinct/;
- sub nfactors { scalar distinct factor(shift); }
+ use Math::Prime::Util qw/pn_primorial factor_exp/;
my $n = pn_primorial(4); # Start with the first 4-factor number
- $n++ while (nfactors($n) != 4 || nfactors($n+1) != 4 || nfactors($n+2) != 4 || nfactors($n+3) != 4);
+ # factor_exp in scalar context returns the number of distinct prime factors
+ $n++ while (factor_exp($n) != 4 || factor_exp($n+1) != 4 || factor_exp($n+2) != 4 || factor_exp($n+3) != 4);
say $n;
-Project Euler, problem 69, stupid brute force solution (about 2 seconds):
+Project Euler, problem 69, stupid brute force solution (about 1 second):
use Math::Prime::Util qw/euler_phi/;
my ($n, $max) = (0,0);
@@ -4127,7 +4132,7 @@ Here is 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, ~4 minutes:
+Project Euler, problem 187, stupid brute force solution, ~3 minutes:
use Math::Prime::Util qw/factor/;
my $nsemis = 0;
@@ -4161,13 +4166,15 @@ Compute L<OEIS A054903|http://oeis.org/A054903> just like CRG4's Pari example:
Construct the table shown in L<OEIS A046147|http://oeis.org/A046147>:
- use Math::Prime::Util qw/znorder euler_phi/;
+ use Math::Prime::Util qw/znorder euler_phi gcd/;
foreach my $n (1..100) {
- my $phi = euler_phi($n);
- my @r = grep { Math::BigInt::bgcd($_,$n) == 1
- && znorder($_,$n) == $phi
- } 1..$n-1;
- say "$n ", join(" ", @r);
+ if (!znprimroot($n)) {
+ say "$n -";
+ } else {
+ my $phi = euler_phi($n);
+ my @r = grep { gcd($_,$n) == 1 && znorder($_,$n) == $phi } 1..$n-1;
+ say "$n ", join(" ", @r);
+ }
}
=head1 PRIMALITY TESTING NOTES
diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm
index b1b8894..98ea9bd 100644
--- a/lib/Math/Prime/Util/PP.pm
+++ b/lib/Math/Prime/Util/PP.pm
@@ -895,6 +895,11 @@ sub _powmod {
$t;
}
+sub gcd {
+ # As usual, bend over backwards to work around RT71548
+ return Math::BigInt::bgcd(map { ($_ < 2147483647 || ref($_) eq 'Math::BigInt') ? $_ : "$_" } @_);
+}
+# unsigned, no validation
sub _gcd_ui {
my($x, $y) = @_;
if ($y < $x) { ($x, $y) = ($y, $x); }
diff --git a/t/19-moebius.t b/t/19-moebius.t
index 321ca30..6b2299f 100644
--- a/t/19-moebius.t
+++ b/t/19-moebius.t
@@ -6,7 +6,7 @@ use Test::More;
use Math::Prime::Util
qw/moebius mertens euler_phi jordan_totient divisor_sum exp_mangoldt
chebyshev_theta chebyshev_psi carmichael_lambda znorder liouville
- znprimroot kronecker legendre_phi
+ znprimroot kronecker legendre_phi gcd
/;
my $extra = defined $ENV{EXTENDED_TESTING} && $ENV{EXTENDED_TESTING};
@@ -274,6 +274,30 @@ my @legendre_sums = (
[1000000, 168, 78331],
);
+my @gcds = (
+ [ [], 0],
+ [ [8], 8],
+ [ [9,9], 9],
+ [ [0,0], 0],
+ [ [1, 0, 0], 1],
+ [ [0, 0, 1], 1],
+ [ [17,19], 1 ],
+ [ [54,24], 6 ],
+ [ [42,56], 14],
+ [ [ 9,28], 1 ],
+ [ [48,180], 12],
+ [ [2705353758,2540073744,3512215098,2214052398], 18],
+ [ [2301535282,3609610580,3261189640], 106],
+ [ [694966514,510402262,195075284,609944479], 181],
+ [ [294950648,651855678,263274296,493043500,581345426], 58 ],
+ [ [-30,-90,90], 30],
+ [ [-3,-9,-18], 3],
+);
+if ($use64) {
+ push @gcds, [ [12848174105599691600,15386870946739346600,11876770906605497900], 700];
+ push @gcds, [ [9785375481451202685,17905669244643674637,11069209430356622337], 117];
+}
+
# These are slow with XS, and *really* slow with PP.
if (!$usexs) {
%big_mertens = map { $_ => $big_mertens{$_} }
@@ -301,6 +325,7 @@ plan tests => 0 + 1
+ 7 + scalar(keys %totients)
+ 1 # Small Carmichael Lambda
+ scalar(@kroneckers)
+ + scalar(@gcds)
+ scalar(@mult_orders)
+ scalar(@legendre_sums)
+ scalar(keys %primroots) + 2
@@ -450,6 +475,12 @@ foreach my $karg (@kroneckers) {
my $k = kronecker($a, $n);
is( $k, $exp, "kronecker($a, $n) = $exp" );
}
+###### gcd
+foreach my $garg (@gcds) {
+ my($aref, $exp) = @$garg;
+ my $gcd = gcd(@$aref);
+ is( $gcd, $exp, "gcd(".join(",",@$aref).") = $exp" );
+}
###### znorder
foreach my $moarg (@mult_orders) {
my ($a, $n, $exp) = @$moarg;
diff --git a/t/32-iterators.t b/t/32-iterators.t
index 9823f30..c810bb7 100644
--- a/t/32-iterators.t
+++ b/t/32-iterators.t
@@ -13,8 +13,8 @@ my $use64 = Math::Prime::Util::prime_get_config->{'maxbits'} > 32;
my $broken64 = (18446744073709550592 == ~0);
plan tests => 8 # forprimes errors
- + 12 + 5 # forprimes simple
- + 1 # forcomposites simple
+ + 12 + 6 # forprimes simple
+ + 3 # forcomposites simple
+ 2 # fordivisors simple
+ 3 # iterator errors
+ 7 # iterator simple
@@ -64,9 +64,18 @@ ok(!eval { forprimes { 1 } 5.6; }, "forprimes abc");
{ my @t; forprimes { push @t, $_ } 31398, 31468;
is_deeply( [@t], [], "forprimes 31398,31468 (empty region)" );
}
+{ my @t; forprimes { push @t, $_ } 2147483647,2147483659;
+ is_deeply( [@t], [2147483647,2147483659], "forprimes 2147483647,2147483659" );
+}
+{ my @t; forcomposites { push @t, $_ } 2147483647,2147483659;
+ is_deeply( [@t], [qw/2147483648 2147483649 2147483650 2147483651 2147483652 2147483653 2147483654 2147483655 2147483656 2147483657 2147483658/], "forcomposites 2147483647,2147483659" );
+}
{ my @t; forcomposites { push @t, $_ } 50;
is_deeply( [@t], [qw/4 6 8 9 10 12 14 15 16 18 20 21 22 24 25 26 27 28 30 32 33 34 35 36 38 39 40 42 44 45 46 48 49 50/], "forcomposites 50" );
}
+{ my @t; forcomposites { push @t, $_ } 200,410;
+ is_deeply( [@t], [qw/200 201 202 203 204 205 206 207 208 209 210 212 213 214 215 216 217 218 219 220 221 222 224 225 226 228 230 231 232 234 235 236 237 238 240 242 243 244 245 246 247 248 249 250 252 253 254 255 256 258 259 260 261 262 264 265 266 267 268 270 272 273 274 275 276 278 279 280 282 284 285 286 287 288 289 290 291 292 294 295 296 297 298 299 300 301 302 303 304 305 306 308 309 310 312 314 315 316 318 319 320 321 322 323 324 325 326 327 328 329 330 332 333 334 335 336 338 3 [...]
+}
{
my $a = 0;
fordivisors { $a += $_ + $_*$_ } 54321;
diff --git a/t/81-bignum.t b/t/81-bignum.t
index eab4985..86f1809 100644
--- a/t/81-bignum.t
+++ b/t/81-bignum.t
@@ -78,6 +78,7 @@ plan tests => 0
+ scalar(keys %allfactors)
+ 13 # moebius, euler_phi, jordan totient, divsum, znorder, etc.
+ 2 # liouville
+ + 3 # gcd
+ 15 # random primes
+ 2 # miller-rabin random
+ 1;
@@ -107,6 +108,7 @@ use Math::Prime::Util qw/
znorder
znprimroot
liouville
+ gcd
pn_primorial
ExponentialIntegral
LogarithmicIntegral
@@ -252,6 +254,10 @@ SKIP: {
###############################################################################
is( liouville( 560812147176208202656339069), -1, "liouville(a x b x c) = -1" );
is( liouville(10571644062695614514374497899), 1, "liouville(a x b x c x d) = 1" );
+###############################################################################
+is( gcd(921166566073002915606255698642,1168315374100658224561074758384,951943731056111403092536868444), 14, "gcd(a,b,c)" );
+is( gcd(1214969109355385138343690512057521757303400673155500334102084,1112036111724848964580068879654799564977409491290450115714228), 42996, "gcd(a,b)" );
+is( gcd(745845206184162095041321,61540282492897317017092677682588744425929751009997907259657808323805386381007), 1, "gcd of two primes = 1" );
###############################################################################
--
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