[libmath-prime-util-perl] 03/72: Add carmichael_lambda; use F&T A2 for random_nbit_prime
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:49:35 UTC 2015
This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.32
in repository libmath-prime-util-perl.
commit 1ef154a158c6f35bc213a1300c0fa5abe7fa1f7a
Author: Dana Jacobsen <dana at acm.org>
Date: Tue Aug 13 16:40:56 2013 -0700
Add carmichael_lambda; use F&T A2 for random_nbit_prime
---
Changes | 3 +-
TODO | 2 ++
lib/Math/Prime/Util.pm | 92 +++++++++++++++++++++++++++++++++++++++++++++++
lib/Math/Prime/Util/PP.pm | 5 +++
t/19-moebius.t | 11 +++++-
5 files changed, 111 insertions(+), 2 deletions(-)
diff --git a/Changes b/Changes
index b467c1c..a0bd1b4 100644
--- a/Changes
+++ b/Changes
@@ -5,8 +5,9 @@ Revision history for Perl module Math::Prime::Util
[Functions Added]
- is_proven_prime
- is_proven_prime_with_cert
+ - carmichael_lambda
- - random_nbit_prime changed to Fouque and Tibouchi A1.
+ - random_nbit_prime now uses Fouque and Tibouchi A1 and A2.
0.31 2013-08-07
diff --git a/TODO b/TODO
index 146dc4c..394cda7 100644
--- a/TODO
+++ b/TODO
@@ -42,3 +42,5 @@
- Refactor where functions exist in .c. Lots of primality tests in factor.c.
- Rewrite 23-primality-proofs.t for new format (keep some of the old tests?).
+
+- Document carmichael_lambda
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index 2805e4c..621f42d 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -40,6 +40,7 @@ our @EXPORT_OK =
moebius mertens euler_phi jordan_totient exp_mangoldt
chebyshev_theta chebyshev_psi
divisor_sum
+ carmichael_lambda
ExponentialIntegral LogarithmicIntegral RiemannZeta RiemannR
);
our %EXPORT_TAGS = (all => [ @EXPORT_OK ]);
@@ -828,10 +829,72 @@ sub primes {
return $_random_prime->($low, $high);
}
+ my @_random_nbit_m;
+ my @_random_nbit_lambda;
+ my @_random_nbit_arange;
+
sub random_nbit_prime {
my($bits) = @_;
_validate_num($bits, 2) || _validate_positive_integer($bits, 2);
+ # Fouque and Tobouchi (2011) Algorithm 2
+ if (1 && $bits > 256) {
+ if (!defined $Math::BigInt::VERSION) {
+ eval { require Math::BigInt; Math::BigInt->import(try=>'GMP,Pari'); 1; }
+ or do { croak "Cannot load Math::BigInt"; };
+ }
+ my ($m, $lambda, $arange) = ($_random_nbit_m[$bits],
+ $_random_nbit_lambda[$bits],
+ $_random_nbit_arange[$bits]);
+ if (!defined $m) {
+ # Calculate beta and primorial
+ my $target = $bits - $_Config{'maxbits'};
+ my $beta = 2;
+ $m = Math::BigInt->new(2);
+ while ($m->copy->blog(2)->badd(1) <= $target) {
+ $beta = next_prime($beta);
+ $m *= $beta;
+ }
+ # Calculate Carmichael Lambda (used to create b) and arange.
+ $lambda = Math::BigInt::blcm( map { $_-1 } @{primes(3, $beta)} );
+ $arange = Math::BigInt->new(2)->bpow($bits)->bdiv($m)->bsub(1);
+ my $arange_bits = $arange->copy->blog(2)->badd(1);
+ die "Incorrect arange" if $arange_bits > $_Config{'maxbits'};
+ $arange = int($arange->bstr);
+ #print "For $bits, I got arange = $arange_bits using primorial($beta) primes\n";
+ ($_random_nbit_m[$bits],
+ $_random_nbit_lambda[$bits],
+ $_random_nbit_arange[$bits]) = ($m, $lambda, $arange);
+ }
+ my $irandf = _get_rand_func();
+ my $b = $irandf->($m-2) + 1;
+ while (1) {
+ my $u = Math::BigInt->bone->bsub($b->copy->bmodpow($lambda, $m))->bmod($m);
+ last if $u == 0;
+ my $r = $irandf->($m-2) + 1;
+ $b->badd($r * $u)->bmod($m);
+ }
+ _make_big_gcds() if $_big_gcd_use < 0;
+ my $loop_limit = 1_000_000;
+ while ($loop_limit-- > 0) {
+ my $a = $irandf->($arange);
+ my $p = $m * $a + $b;
+ if ($_HAVE_GMP) {
+ next unless Math::Prime::Util::GMP::is_prime($p);
+ } else {
+ if ($_big_gcd_use && $p > $_big_gcd_top) {
+ next unless Math::BigInt::bgcd($p, $_big_gcd[0]) == 1;
+ next unless Math::BigInt::bgcd($p, $_big_gcd[1]) == 1;
+ next unless Math::BigInt::bgcd($p, $_big_gcd[2]) == 1;
+ next unless Math::BigInt::bgcd($p, $_big_gcd[3]) == 1;
+ }
+ next unless is_prob_prime($p);
+ }
+ return $p;
+ }
+ croak "Random function broken?";
+ }
+
# Fouque and Tibouchi (2011) Algorithm 1 (basic)
#
# Example for random_nbit_prime(512) on 64-bit Perl:
@@ -1504,6 +1567,35 @@ sub chebyshev_psi {
return $sum;
}
+sub carmichael_lambda {
+ my($n) = @_;
+ _validate_num($n) || _validate_positive_integer($n);
+ # lambda(n) = phi(n) for n < 8
+ return euler_phi($n) if $n < 8;
+ # lambda(n) = phi(n)/2 for powers of two greater than 4
+ return euler_phi($n)/2 if ($n & ($n-1)) == 0;
+
+ my %factor_mult;
+ my @factors = grep { !$factor_mult{$_}++ }
+ ($n <= $_XS_MAXVAL) ? _XS_factor($n) : factor($n);
+ $factor_mult{2}-- if defined $factor_mult{2} && $factor_mult{2} > 2;
+
+ if (!defined $Math::BigInt::VERSION) {
+ eval { require Math::BigInt; Math::BigInt->import(try=>'GMP,Pari'); 1; }
+ or do { croak "Cannot load Math::BigInt"; };
+ }
+ my $lcm = Math::BigInt::blcm(
+ map { Math::BigInt->new("$_")
+ ->bpow($factor_mult{$_}-1)
+ ->bmul($_-1)
+ } @factors
+ );
+ $lcm = int($lcm->bstr) if $lcm->bacmp(''.~0) <= 0;
+ return $lcm;
+}
+
+
+
#############################################################################
diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm
index 5c27899..e2036ce 100644
--- a/lib/Math/Prime/Util/PP.pm
+++ b/lib/Math/Prime/Util/PP.pm
@@ -723,6 +723,11 @@ sub _gcd_ui {
$x;
}
+sub _lcm_ui {
+ my($x, $y) = @_;
+ return (abs($x) / _gcd_ui($x, $y)) * abs($y);
+}
+
sub _is_perfect_power {
my $n = shift;
return 0 if $n <= 3 || $n != int($n);
diff --git a/t/19-moebius.t b/t/19-moebius.t
index e1a870a..057f8ff 100644
--- a/t/19-moebius.t
+++ b/t/19-moebius.t
@@ -5,7 +5,7 @@ use warnings;
use Test::More;
use Math::Prime::Util
qw/moebius mertens euler_phi jordan_totient divisor_sum exp_mangoldt
- chebyshev_theta chebyshev_psi/;
+ chebyshev_theta chebyshev_psi carmichael_lambda/;
my $extra = defined $ENV{EXTENDED_TESTING} && $ENV{EXTENDED_TESTING};
my $use64 = Math::Prime::Util::prime_get_config->{'maxbits'} > 32;
@@ -155,6 +155,8 @@ my %chebyshev2 = (
1234567 => 1234515.17962833,
);
+my @A002322 = (0,1,1,2,2,4,2,6,2,6,4,10,2,12,6,4,4,16,6,18,4,6,10,22,2,20,12,18,6,28,4,30,8,10,16,12,6,36,18,12,4,40,6,42,10,12,22,46,4,42,20,16,12,52,18,20,6,18,28,58,4,60,30,6,16,12,10,66,16,22,12,70,6,72,36,20,18,30,12,78,4,54,40,82,6,16,42,28,10,88,12,12,22,30,46,36,8,96,42,30,20,100,16,102,12,12,52,106,18,108,20,36,12,112,18,44,28,12,58,48,4,110,60,40,30,100,6,126,32,42,12,130,10,18,66,36,16,136,22,138,12,46,70,60,12,28,72,42,36,148,20,150,18,48,30,60,12,156,78,52,8,66,54,162,40,20, [...]
+
plan tests => 0 + 1
+ 1 # Small Moebius
@@ -162,6 +164,7 @@ plan tests => 0 + 1
+ 1*scalar(keys %big_mertens)
+ 2 # Small Phi
+ 7 + scalar(keys %totients)
+ + 1 # Small Carmichael Lambda
+ scalar(keys %jordan_totients)
+ 2 # Dedekind psi calculated two ways
+ 1 # Calculate J5 two different ways
@@ -275,6 +278,12 @@ while (my($n, $c2) = each (%chebyshev2)) {
cmp_closeto( chebyshev_psi($n), $c2, 1e-9*abs($n), "chebyshev_psi($n)" );
}
+###### Carmichael Lambda
+{
+ my @lambda = map { carmichael_lambda($_) } (0 .. $#A002322);
+ is_deeply( \@lambda, \@A002322, "carmichael_lambda with range: 0, $#A000010" );
+}
+
sub cmp_closeto {
my $got = shift;
--
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