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