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 <d...@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

_______________________________________________
Pkg-perl-cvs-commits mailing list
Pkg-perl-cvs-commits@lists.alioth.debian.org
http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/pkg-perl-cvs-commits

Reply via email to