[libmath-prime-util-perl] 51/181: Add kronecker, znprimroot; XSify carmichael_lambda
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:51:05 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 47af0e2ab81110baf4c6a15015cdcc8a564cf55b
Author: Dana Jacobsen <dana at acm.org>
Date: Fri Dec 27 18:28:29 2013 -0800
Add kronecker, znprimroot; XSify carmichael_lambda
---
Changes | 9 +++++
lib/Math/Prime/Util.pm | 93 +++++++++++++++++++++++++++++++++++++++++------
lib/Math/Prime/Util/PP.pm | 73 +++++++++++++++++++------------------
3 files changed, 128 insertions(+), 47 deletions(-)
diff --git a/Changes b/Changes
index 9d2d570..2eff5b0 100644
--- a/Changes
+++ b/Changes
@@ -16,9 +16,14 @@ 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
[FUNCTIONALITY AND PERFORMANCE]
+ - Speedup for input number validation: speedup for all functions,
+ especially noticeable with fast functions (e.g. small n is_prime(n)).
+
- Some 5.6.2-is-broken workarounds.
- Some LMO edge cases: small numbers that only show up if a #define is
@@ -36,6 +41,10 @@ Revision history for Perl module Math::Prime::Util
- Rewrite sieve composite map. Segment sieving is faster. It's a little
faster than primegen for me, but still slower than primesieve and yafu.
+ - znorder uses Carmichael Lambda instead of Euler Phi. Faster.
+
+ - carmichael_lambda switched to XS->Perl from Perl->XS.
+
0.35 2013-12-08
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index 42b6b64..5f43c73 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -42,7 +42,7 @@ our @EXPORT_OK =
partitions
chebyshev_theta chebyshev_psi
divisor_sum
- carmichael_lambda znorder
+ carmichael_lambda kronecker znorder znprimroot
ExponentialIntegral LogarithmicIntegral RiemannZeta RiemannR
);
our %EXPORT_TAGS = (all => [ @EXPORT_OK ]);
@@ -105,6 +105,9 @@ BEGIN {
*next_prime = \&Math::Prime::Util::_generic_next_prime;
*prev_prime = \&Math::Prime::Util::_generic_prev_prime;
*exp_mangoldt = \&Math::Prime::Util::_generic_exp_mangoldt;
+ *carmichael_lambda = \&Math::Prime::Util::_generic_carmichael_lambda;
+ *kronecker = \&Math::Prime::Util::_generic_kronecker;
+ *znprimroot = \&Math::Prime::Util::_generic_znprimroot;
*forprimes = sub (&$;$) { _generic_forprimes(@_); }; ## no critic qw(ProhibitSubroutinePrototypes)
*fordivisors = sub (&$) { _generic_fordivisors(@_); }; ## no critic qw(ProhibitSubroutinePrototypes)
*forcomposites = sub (&$) { _generic_forcomposites(@_); }; ## no critic qw(ProhibitSubroutinePrototypes)
@@ -1595,8 +1598,6 @@ sub _generic_exp_mangoldt {
sub liouville {
my($n) = @_;
- # Note the special behavior for n = 1
- return 1 if defined $n && $n <= 1;
_validate_num($n) || _validate_positive_integer($n);
my $l = ($n <= $_XS_MAXVAL) ? (-1) ** scalar _XS_factor($n)
: (-1) ** scalar factor($n);
@@ -1657,7 +1658,7 @@ sub chebyshev_psi {
return $sum;
}
-sub carmichael_lambda {
+sub _generic_carmichael_lambda {
my($n) = @_;
_validate_num($n) || _validate_positive_integer($n);
# lambda(n) = phi(n) for n < 8
@@ -1725,6 +1726,30 @@ sub znorder {
return $k;
}
+sub _generic_znprimroot {
+ my($n) = @_;
+ _validate_num($n) || _validate_positive_integer($n);
+ return if $n <= 0;
+ return $n-1 if $n <= 4;
+ my $a = 1;
+ my $phi = euler_phi($n);
+ my @exp = map { int($phi/$_->[0]) } factor_exp($phi);
+ #print "phi: $phi factors: ", join(",",factor($phi)), "\n";
+ #print " exponents: ", join(",", @exp), "\n";
+ while (1) {
+ my $fail = 0;
+ do { $a++; } while kronecker($a,$n) == 0;
+ return if $a >= $n;
+ foreach my $f (@exp) {
+ # As usual, quotes for RT 71548
+ my $e = Math::BigInt->new($a)->bmodpow("$f", "$n");
+ #print " $a^$f mod $n = $e\n";
+ if ($e == 1) { $fail = 1; last; }
+ }
+ return $a if !$fail;
+ }
+}
+
@@ -1807,6 +1832,19 @@ sub _generic_prev_prime {
return Math::Prime::Util::PP::prev_prime($_[0]);
}
+sub _generic_kronecker {
+ my($a, $b) = @_;
+ croak "Parameter must be defined" if !defined $a;
+ croak "Parameter must be defined" if !defined $b;
+ croak "Parameter '$a' must be an integer" unless $a =~ /^-?\d+/;
+ croak "Parameter '$b' must be an integer" unless $b =~ /^-?\d+/;
+
+ return Math::BigInt->new(''.Math::Prime::Util::GMP::kronecker($a,$b))
+ if $_HAVE_GMP && defined &Math::Prime::Util::GMP::kronecker;
+
+ return Math::Prime::Util::PP::kronecker(@_);
+}
+
sub prime_count {
my($low,$high) = @_;
if (defined $high) {
@@ -3788,6 +3826,22 @@ or Carmichael λ(n)) of a positive integer argument. It is the smallest
positive integer m such that a^m = 1 mod n for every integer a coprime to n.
This is L<OEIS series A002322|http://oeis.org/A002322>.
+=head2 kronecker
+
+Returns the Kronecker symbol C<(a|n)> for two integers. The possible
+return values with their meanings for odd positive C<n> are:
+
+ 0 a = 0 mod n
+ 1 a is a quadratic residue modulo n (a = x^2 mod n for some x)
+ -1 a is a quadratic non-residue modulo n
+
+and the return value is congruent to C<a^((n-1)/2)>. The Kronecker
+symbol is an extension of the Jacobi symbol to all integer values of
+C<n> from its domain of positive odd values of C<n>. The Jacobi
+symbol is itself an extension of the Legendre symbol, which is
+only defined for odd prime values of C<n>. This corresponds to Pari's
+C<kronecker(a,n)> function and Mathematica's C<KroneckerSymbol[n,m]>
+function.
=head2 znorder
@@ -3800,6 +3854,14 @@ C<a> and C<n> are not coprime, since no value will result in 1 mod n.
This corresponds to Pari's C<znorder(Mod(a,n))> function and Mathematica's
C<MultiplicativeOrder[n]> function.
+=head2 znprimroot
+
+Given a positive integer C<n>, returns a primitive root of C<(Z/nZ)^*>.
+A root exists when C<euler_phi($n) == carmichael_lambda($n)>, which will
+be true for all prime C<n> and some composites.
+(L<OEIS A033948|http://oeis.org/A033948>) is the list of such values.
+If a primitive root does not exist for C<n>, this function returns undef.
+
=head2 random_prime
@@ -4780,16 +4842,17 @@ Similar to MPU's L</divisors>.
Similar to MPU's L</forprimes>, L</forcomposites>, L<fordivisors>, and
L<divisor_sum>.
-=item C<eulerphi>
+=item C<eulerphi>, C<moebius>
-Similar to MPU's L</euler_phi>. MPU is 2-20x faster for native integers.
-There is also support for a range, which can be much more efficient.
-Without L<Math::Prime::Util::GMP> installed, MPU is very slow with bigints.
-With it installed, it is about 2x slower than Math::Pari.
+Similar to MPU's L</euler_phi> and L</moebius>. MPU is 2-20x faster for
+native integers. There is also support for a range, which can be much
+more efficient. Without L<Math::Prime::Util::GMP> installed, MPU is
+very slow with bigints. With it installed, it is about 2x slower than
+Math::Pari.
-=item C<moebius>
+=item C<kronecker>, C<znorder>, C<znprimroot>
-Similar to MPU's L</moebius>. Comparisons are similar to C<eulerphi>.
+Similar to MPU's L</kronecker>, L</znorder>, and L</znprimroot>.
=item C<sigma>
@@ -5124,6 +5187,14 @@ thank Kim Walisch for the many discussions about prime counting.
=item *
+Henri Cohen, "A Course in Computational Algebraic Number Theory", Springer, 1996. Practical computational number theory from the team lead of L<Pari|http://pari.math.u-bordeaux.fr/>. Lots of explicit algorithms.
+
+=item *
+
+Hans Riesel, "Prime Numbers and Computer Methods for Factorization", Birkh?user, 2nd edition, 1994. Lots of information, some code, easy to follow.
+
+=item *
+
Pierre Dusart, "Estimates of Some Functions Over Primes without R.H.", preprint, 2010. Updates to the best non-RH bounds for prime count and nth prime. L<http://arxiv.org/abs/1002.0442/>
=item *
diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm
index 674fb22..9911658 100644
--- a/lib/Math/Prime/Util/PP.pm
+++ b/lib/Math/Prime/Util/PP.pm
@@ -854,6 +854,40 @@ sub miller_rabin {
}
1;
}
+# Calculate Kronecker symbol (a|b). Cohen Algorithm 1.4.10.
+# Extension of the Jacobi symbol, itself an extension of the Legendre symbol.
+sub kronecker {
+ my($a, $b) = @_;
+ return (abs($a) == 1) ? 1 : 0 if $b == 0;
+ my $k = 1;
+ if ($b % 2 == 0) {
+ return 0 if $a % 2 == 0;
+ my $v = 0;
+ do { $v++; $b /= 2; } while $b % 2 == 0;
+ $k = -$k if $v % 2 == 1 && ($a % 8 == 3 || $a % 8 == 5);
+ }
+ if ($b < 0) {
+ $b = -$b;
+ $k = -$k if $a < 0;
+ }
+ if ($a < 0) { $a = -$a; $k = -$k if $b % 4 == 3; }
+ # Now: b > 0, b odd, a >= 0
+ while ($a != 0) {
+ if ($a % 2 == 0) {
+ my $v = 0;
+ do { $v++; $a /= 2; } while $a % 2 == 0;
+ $k = -$k if $v % 2 == 1 && ($b % 8 == 3 || $b % 8 == 5);
+ }
+ $k = -$k if $a % 4 == 3 && $b % 4 == 3;
+ ($a, $b) = ($b % $a, $a);
+ # If a,b are bigints and now small enough, finish as native.
+ if ( ref($a) eq 'Math::BigInt' && $a <= ''.~0
+ && ref($b) eq 'Math::BigInt' && $b <= ''.~0) {
+ return $k * kronecker(int($a->bstr), int($b->bstr));
+ }
+ }
+ return ($b == 1) ? $k : 0;
+}
sub _is_perfect_square {
my($n) = @_;
@@ -875,39 +909,6 @@ sub _is_perfect_square {
0;
}
-# Calculate Jacobi symbol (M|N)
-sub _jacobi {
- my($n, $m) = @_;
- return 0 if $m <= 0 || ($m % 2) == 0;
- my $j = 1;
- if ($n < 0) {
- $n = -$n;
- $j = -$j if ($m % 4) == 3;
- }
- # Split loop so we can reduce n/m to non-bigints after first iteration.
- if ($n != 0) {
- while (($n % 2) == 0) {
- $n >>= 1;
- $j = -$j if ($m % 8) == 3 || ($m % 8) == 5;
- }
- ($n, $m) = ($m, $n);
- $j = -$j if ($n % 4) == 3 && ($m % 4) == 3;
- $n = $n % $m;
- $n = int($n->bstr) if ref($n) eq 'Math::BigInt' && $n <= ''.~0;
- $m = int($m->bstr) if ref($m) eq 'Math::BigInt' && $m <= ''.~0;
- }
- while ($n != 0) {
- while (($n % 2) == 0) {
- $n >>= 1;
- $j = -$j if ($m % 8) == 3 || ($m % 8) == 5;
- }
- ($n, $m) = ($m, $n);
- $j = -$j if ($n % 4) == 3 && ($m % 4) == 3;
- $n = $n % $m;
- }
- return ($m == 1) ? $j : 0;
-}
-
# Find first D in sequence (5,-7,9,-11,13,-15,...) where (D|N) == -1
sub _lucas_selfridge_params {
my($n) = @_;
@@ -920,7 +921,7 @@ sub _lucas_selfridge_params {
my $gcd = (ref($n) eq 'Math::BigInt') ? Math::BigInt::bgcd($d, $n)
: _gcd_ui($d, $n);
return (0,0,0) if $gcd > 1 && $gcd != $n; # Found divisor $d
- my $j = _jacobi($d * $sign, $n);
+ my $j = kronecker($d * $sign, $n);
last if $j == -1;
$d += 2;
croak "Could not find Jacobi sequence for $n" if $d > 4_000_000_000;
@@ -941,7 +942,7 @@ sub _lucas_extrastrong_params {
my $gcd = (ref($n) eq 'Math::BigInt') ? Math::BigInt::bgcd($D, $n)
: _gcd_ui($D, $n);
return (0,0,0) if $gcd > 1 && $gcd != $n; # Found divisor $d
- last if _jacobi($D, $n) == -1;
+ last if kronecker($D, $n) == -1;
$P += $increment;
croak "Could not find Jacobi sequence for $n" if $P > 65535;
$D = $P*$P - 4;
@@ -1171,7 +1172,7 @@ sub is_frobenius_underwood_pseudoprime {
my $fb = $ZERO + 2;
my ($x, $t, $np1, $na) = (0, -1, $n+1, undef);
- while ( _jacobi($t, $n) != -1 ) {
+ while ( kronecker($t, $n) != -1 ) {
$x++;
$t = $x*$x - 4;
}
--
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