[libmath-prime-util-perl] 19/181: Second round of factor(0), factor(1) changes
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:51:02 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 9a56656e0c4992dfb1817452b3dfee5b9d969387
Author: Dana Jacobsen <dana at acm.org>
Date: Wed Dec 18 00:17:32 2013 -0800
Second round of factor(0),factor(1) changes
---
Changes | 10 +++++++
TODO | 4 ---
factor.c | 9 +++---
lib/Math/Prime/Util.pm | 71 ++++++++++++++++++++++++-----------------------
lib/Math/Prime/Util/PP.pm | 6 ++--
t/02-can.t | 8 ++++--
t/19-moebius.t | 4 +--
t/50-factoring.t | 51 ++++++++++++++++++++++++++++++++--
t/80-pp.t | 2 +-
t/81-bignum.t | 4 +--
10 files changed, 112 insertions(+), 57 deletions(-)
diff --git a/Changes b/Changes
index 8781b34..0bccb3c 100644
--- a/Changes
+++ b/Changes
@@ -2,6 +2,16 @@ Revision history for Perl module Math::Prime::Util
0.36 2013-12
+ [API Changes]
+
+ - factor behavior for 0 and 1 more consistent. The results for factor,
+ factor_exp, divisors, and divisor_sum should now match Pari, and the
+ omega/Omega exception for 1 is removed.
+
+ Thanks to Hugo van der Sanden for bringing this up.
+
+ - all_factors changed to divisors. The old name still remains aliased.
+
[ADDED]
- forcomposites like forprimes, but on composites. See Pari 2.6.x.
diff --git a/TODO b/TODO
index 793d840..91a2698 100644
--- a/TODO
+++ b/TODO
@@ -65,7 +65,3 @@
- Document forcomposites. Consider XS for it.
- Add fordivisors as XS.
-
-- Make divisors() an alias for all_factors().
-
-- Document changes to factor(1)
diff --git a/factor.c b/factor.c
index d2a6ba3..20a8c09 100644
--- a/factor.c
+++ b/factor.c
@@ -323,9 +323,9 @@ UV* _divisor_list(UV n, UV *num_divisors)
int i, j, nfactors, ndivisors;
if (n <= 1) {
- New(0, divs, 1, UV);
- divs[0] = 1;
- *num_divisors = n;
+ New(0, divs, 2, UV);
+ if (n == 0) { divs[0] = 0; divs[1] = 1; *num_divisors = 2; }
+ if (n == 1) { divs[0] = 1; *num_divisors = 1; }
return divs;
}
/* Factor and convert to factor/exponent pair */
@@ -375,7 +375,8 @@ UV _XS_divisor_sum(UV n, UV k)
int nfac, i, j;
UV product = 1;
- if (n <= 1) return n;
+ if (n == 0) return (k == 0) ? 2 : 1; /* divisors are [0,1] */
+ if (n == 1) return 1; /* divisors are [1] */
nfac = factor(n, factors);
if (k == 0) {
for (i = 0; i < nfac; i++) {
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index 8956d5b..349625f 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
+ factor factor_exp all_factors divisors
moebius mertens euler_phi jordan_totient exp_mangoldt liouville
partitions
chebyshev_theta chebyshev_psi
@@ -62,6 +62,7 @@ sub _import_nobigint {
return unless $_Config{'xs'};
undef *factor; *factor = \&_XS_factor;
undef *factor_exp; *factor_exp = \&_XS_factor_exp;
+ undef *divisors; *divisors = \&_XS_divisors;
#undef *prime_count; *prime_count = \&_XS_prime_count;
undef *nth_prime; *nth_prime = \&_XS_nth_prime;
undef *is_pseudoprime; *is_pseudoprime = \&_XS_is_pseudoprime;
@@ -1219,17 +1220,17 @@ sub consecutive_integer_lcm {
return $pn;
}
-sub all_factors {
+sub divisors {
my $n = shift;
- # In scalar context, returns sigma_0(n). Very fast.
- return divisor_sum($n,0) unless wantarray;
-
_validate_num($n) || _validate_positive_integer($n);
- return () if $n == 1;
return _XS_divisors($n) if $n <= $_XS_MAXVAL;
+ # In scalar context, returns sigma_0(n). Very fast.
+ return divisor_sum($n,0) unless wantarray;
+ return ($n == 0) ? (0,1) : (1) if $n <= 1;
+
my %all_factors;
foreach my $f1 (factor($n)) {
next if $f1 >= $n;
@@ -1247,6 +1248,9 @@ sub all_factors {
return @divisors;
}
+# alias the old "all_factors" to the new name: divisors
+*all_factors = \&_XS_divisors;
+
# A008683 Moebius function mu(n)
# A030059, A013929, A030229, A002321, A005117, A013929 all relate.
@@ -1427,7 +1431,7 @@ my @_ds_overflow = # We'll use BigInt math if the input is larger than this.
: ( 50, 845404560, 52560, 1548, 252, 84);
sub divisor_sum {
my($n, $k) = @_;
- return (0,1)[$n] if defined $n && $n <= 1;
+ return 1 if defined $n && $n == 1;
_validate_num($n) || _validate_positive_integer($n);
# Call the XS routine for k=0 and k=1 immediately if possible.
@@ -1442,7 +1446,7 @@ sub divisor_sum {
if (defined $k && ref($k) eq 'CODE') {
my $sum = $n-$n;
- foreach my $f (all_factors($n)) {
+ foreach my $f (divisors($n)) {
$sum += $k->($f);
}
return $sum;
@@ -1666,8 +1670,8 @@ sub znorder {
return if Math::BigInt::bgcd($a, $n) > 1;
# Method 1: check all a^k 1 .. $n-1.
# Naive and terrible slow.
- # Method 2: check all k in (all_factors(euler_phi($n), $n).
- # Method 3: check all k in (all_factors(carmichael_lambda($n), $n).
+ # Method 2: check all k in (divisors(euler_phi($n), $n).
+ # Method 3: check all k in (divisors(carmichael_lambda($n), $n).
# Good for most cases, but slow when factor quantity is large.
# Method 4: Das algorithm 1.7, just enough multiples of p
# Fastest.
@@ -1693,7 +1697,7 @@ sub znorder {
# Method 3:
# my $cl = carmichael_lambda($n);
- # foreach my $k (all_factors($cl), $cl) {
+ # foreach my $k (divisors($cl), $cl) {
# my $t = $a->copy->bmodpow("$k", $n);
# return $k if $t->is_one;
# }
@@ -1833,13 +1837,15 @@ sub factor {
my($n) = @_;
_validate_num($n) || _validate_positive_integer($n);
- return () if $n == 1;
return _XS_factor($n) if $n <= $_XS_MAXVAL;
if ($_HAVE_GMP) {
- my @factors = Math::Prime::Util::GMP::factor($n);
- if (ref($_[0]) eq 'Math::BigInt') {
- @factors = map { ($_ > ~0) ? Math::BigInt->new(''.$_) : $_ } @factors;
+ my @factors;
+ if ($n != 1) {
+ @factors = Math::Prime::Util::GMP::factor($n);
+ if (ref($_[0]) eq 'Math::BigInt') {
+ @factors = map { ($_ > ~0) ? Math::BigInt->new(''.$_) : $_ } @factors;
+ }
}
return @factors;
}
@@ -2594,7 +2600,7 @@ Version 0.35
@pe = factor_exp( $n );
# Get all divisors other than 1 and n
- @divisors = all_factors( $n );
+ @divisors = divisors( $n );
# Euler phi (Euler's totient) on a large number
use bigint; say euler_phi( 801294088771394680000412 );
@@ -4038,10 +4044,8 @@ Allows setting of some parameters. Currently the only parameters are:
# returns (204518747,16476429743)
Produces the prime factors of a positive number input, in numerical order.
-The special cases of C<n = 0> and C<n = 1> will return C<n>, which
-guarantees multiplying the factors together will always result in the
-input value, though those are the only cases where the returned factors
-are not prime.
+The product of the returned factors will be equal to the input. C<n = 1>
+will return an empty list, and C<n = 0> will return 0. This matches Pari.
In scalar context, returns Ω(n), the total number of prime factors
(L<OEIS A001222|http://oeis.org/A001222>).
@@ -4049,8 +4053,6 @@ This corresponds to Pari's C<bigomega(n)> function and Mathematica's
C<PrimeOmega[n]> function.
This is same result that we would get if we evaluated the resulting
array in scalar context.
-Do note that the inputs of C<0> and C<1> will return C<1>, contrary
-to the standard definition of Omega.
The current algorithm for non-bigints is a sequence of small trial division,
a few rounds of Pollard's Rho, SQUFOF, Pollard's p-1, Hart's OLF, a long
@@ -4087,24 +4089,24 @@ This corresponds to Pari's C<omega(n)> function and Mathematica's
C<PrimeNu[n]> function.
This is same result that we would get if we evaluated the resulting
array in scalar context.
-Do note that the inputs of C<0> and C<1> will return C<1>, contrary
-to the standard definition of omega.
The internals are identical to L</factor>, so all comments there apply.
Just the way the factors are arranged is different.
-=head2 all_factors
+=head2 divisors
- my @divisors = all_factors(30); # returns (1, 2, 3, 5, 6, 10, 15, 30)
+ my @divisors = divisors(30); # returns (1, 2, 3, 5, 6, 10, 15, 30)
Produces all the divisors of a positive number input, including 1 and
the input number. The divisors are a power set of multiplications of
the prime factors, returned as a uniqued sorted list. The result is
-identical to that of Pari's C<divisors> function.
+identical to that of Pari's C<divisors> and Mathematica's C<Divisors[n]>
+functions.
In scalar context this returns the sigma0 function,
the sigma function (see Hardy and Wright section 16.7, or OEIS A000203).
+This is the same result as evaluating the array in scalar context.
=head2 trial_factor
@@ -4112,8 +4114,7 @@ the sigma function (see Hardy and Wright section 16.7, or OEIS A000203).
my @factors = trial_factor($n);
Produces the prime factors of a positive number input. The factors will be
-in numerical order. The special cases of C<n = 0> and C<n = 1> will return
-C<n>, while with all other inputs the factors are guaranteed to be prime.
+in numerical order.
For large inputs this will be very slow.
=head2 fermat_factor
@@ -4430,11 +4431,11 @@ Here is the best way for PE187. Under 30 milliseconds from the command line:
Produce the C<matches> result from L<Math::Factor::XS> without skipping:
- use Math::Prime::Util qw/all_factors/;
+ use Math::Prime::Util qw/divisors/;
use Algorithm::Combinatorics qw/combinations_with_repetition/;
my $n = 139650;
my @matches = grep { $_->[0] * $_->[1] == $n && $_->[0] > 1 }
- combinations_with_repetition( [all_factors($n)], 2 );
+ combinations_with_repetition( [divisors($n)], 2 );
=head1 PRIMALITY TESTING NOTES
@@ -4597,7 +4598,7 @@ for security. MPU can return a primality certificate.
What Crypt::Primes has that MPU does not is the ability to return a generator.
L<Math::Factor::XS> calculates prime factors and factors, which correspond to
-the L</factor> and L</all_factors> functions of MPU. These functions do
+the L</factor> and L</divisors> functions of MPU. These functions do
not support bigints. Both are implemented with trial division, meaning they
are very fast for really small values, but quickly become unusably slow
(factoring 19 digit semiprimes is over 700 times slower). The function
@@ -4702,8 +4703,8 @@ doesn't support segmenting.
=item C<factorint>
-Similar to MPU's L</factor> though with a different return. MPU offers
-L</factor> for a linear array of prime factors where
+Similar to MPU's L</factor_exp> though with a slightly different return.
+MPU offers L</factor> for a linear array of prime factors where
n = p1 * p2 * p3 * ... as (p1,p2,p3,...)
and L</factor_exp> for an array of factor/exponent pairs where:
n = p1^e1 * p2^e2 * ... as ([p1,e1],[p2,e2],...)
@@ -4717,7 +4718,7 @@ faster on average in MPU.
=item C<divisors>
-Similar to MPU's L</all_factors>.
+Similar to MPU's L</divisors>.
=item C<eulerphi>
diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm
index 2de4bda..674fb22 100644
--- a/lib/Math/Prime/Util/PP.pm
+++ b/lib/Math/Prime/Util/PP.pm
@@ -1328,7 +1328,7 @@ sub is_aks_prime {
sub _basic_factor {
# MODIFIES INPUT SCALAR
- return ($_[0]) if $_[0] < 4;
+ return ($_[0] == 1) ? () : ($_[0]) if $_[0] < 4;
my @factors;
if (ref($_[0]) ne 'Math::BigInt') {
@@ -1365,7 +1365,7 @@ sub trial_factor {
# Don't use _basic_factor here -- they want a trial forced.
my @factors;
if ($n < 4) {
- @factors = ($n);
+ @factors = ($n == 1) ? () : ($n);
return @factors;
}
while ( !($n % 2) ) { push @factors, 2; $n = int($n / 2); }
@@ -2996,7 +2996,7 @@ Operations that are slower include:
primes
random_prime / random_ndigit_prime
- factor / all_factors
+ factor / factor_exp / divisors
nth_prime
primecount
diff --git a/t/02-can.t b/t/02-can.t
index 160daa9..094fe71 100644
--- a/t/02-can.t
+++ b/t/02-can.t
@@ -20,7 +20,8 @@ my @functions = qw(
miller_rabin_random
lucas_sequence
primes
- forprimes prime_iterator prime_iterator_object
+ forprimes forcomposites
+ prime_iterator prime_iterator_object
next_prime prev_prime
prime_count
prime_count_lower prime_count_upper prime_count_approx
@@ -29,8 +30,9 @@ my @functions = qw(
random_proven_prime random_proven_prime_with_cert
random_maurer_prime random_maurer_prime_with_cert
primorial pn_primorial consecutive_integer_lcm
- factor all_factors
- moebius mertens euler_phi jordan_totient exp_mangoldt
+ factor factor_exp divisors
+ moebius mertens euler_phi jordan_totient exp_mangoldt liouville
+ partitions
chebyshev_theta chebyshev_psi
divisor_sum
carmichael_lambda znorder
diff --git a/t/19-moebius.t b/t/19-moebius.t
index 941b2ff..fc0efd7 100644
--- a/t/19-moebius.t
+++ b/t/19-moebius.t
@@ -289,8 +289,8 @@ while (my($k, $tref) = each (%jordan_totients)) {
}
{
- my @J5_moebius = map { divisor_sum($_, sub { my $d=shift; $d**5 * moebius($_/$d); }) } (0 .. 100);
- my @J5_jordan = map { jordan_totient(5, $_) } (0 .. 100);
+ my @J5_moebius = map { divisor_sum($_, sub { my $d=shift; $d**5 * moebius($_/$d); }) } (1 .. 100);
+ my @J5_jordan = map { jordan_totient(5, $_) } (1 .. 100);
is_deeply( \@J5_moebius, \@J5_jordan, "Calculate J5 two different ways");
}
diff --git a/t/50-factoring.t b/t/50-factoring.t
index c9cab00..a3af134 100644
--- a/t/50-factoring.t
+++ b/t/50-factoring.t
@@ -3,7 +3,7 @@ use strict;
use warnings;
use Test::More;
-use Math::Prime::Util qw/factor factor_exp all_factors is_prime/;
+use Math::Prime::Util qw/factor factor_exp all_factors divisor_sum is_prime/;
my $use64 = Math::Prime::Util::prime_get_config->{'maxbits'} > 32;
my $extra = defined $ENV{EXTENDED_TESTING} && $ENV{EXTENDED_TESTING};
@@ -79,11 +79,43 @@ my %all_factors = (
4 => [1,2,4],
3 => [1,3],
2 => [1,2],
+ 1 => [1],
+ 0 => [0,1],
+);
+
+my %prime_factors = (
+ 456789 => [3,43,3541],
+ 123456 => [2,2,2,2,2,2,3,643],
+ 115553 => [115553],
+ 30107 => [7,11,17,23],
+ 5 => [5],
+ 4 => [2,2],
+ 3 => [3],
+ 2 => [2],
+ 1 => [],
+ 0 => [0],
+);
+
+my %factor_exponents = (
+ 456789 => [[3,1],[43,1],[3541,1]],
+ 123456 => [[2,6],[3,1],[643,1]],
+ 115553 => [[115553,1]],
+ 30107 => [[7,1],[11,1],[17,1],[23,1]],
+ 5 => [[5,1]],
+ 4 => [[2,2]],
+ 3 => [[3,1]],
+ 2 => [[2,1]],
1 => [],
- 0 => [],
+ 0 => [[0,1]],
);
-plan tests => (3 * scalar @testn) + scalar(keys %all_factors) + 10*9 + 8 + 1;
+plan tests => (3 * scalar @testn)
+ + 2*scalar(keys %prime_factors)
+ + 4*scalar(keys %all_factors)
+ + 2*scalar(keys %factor_exponents)
+ + 10*9
+ + 8
+ + 1;
foreach my $n (@testn) {
my @f = factor($n);
@@ -105,10 +137,23 @@ foreach my $n (@testn) {
is_deeply( [factor_exp($n)], [linear_to_exp(@f)], " factor_exp looks right" );
}
+while (my($n, $factors) = each(%prime_factors)) {
+ is_deeply( [factor($n)], $factors, "factors($n)" );
+ is( scalar factor($n), scalar @$factors, "scalar factors($n)" );
+}
+
while (my($n, $divisors) = each(%all_factors)) {
is_deeply( [all_factors($n)], $divisors, "all_factors($n)" );
+ is( scalar all_factors($n), scalar @$divisors, "scalar all_factors($n)" );
+ is( divisor_sum($n,0), scalar @$divisors, "divisor_sum($n,0)" );
+ my $sum = 0; foreach my $f (@$divisors) { $sum += $f; }
+ is( divisor_sum($n), $sum, "divisor_sum($n)" );
}
+while (my($n, $factors) = each(%factor_exponents)) {
+ is_deeply( [factor_exp($n)], $factors, "factor_exp($n)" );
+ is( scalar factor_exp($n), scalar @$factors, "scalar factor_exp($n)" );
+}
extra_factor_test("trial_factor", sub {Math::Prime::Util::trial_factor(shift)});
extra_factor_test("fermat_factor", sub {Math::Prime::Util::fermat_factor(shift)});
diff --git a/t/80-pp.t b/t/80-pp.t
index bb25783..6e113a3 100644
--- a/t/80-pp.t
+++ b/t/80-pp.t
@@ -29,7 +29,7 @@ my @primes = qw/
/;
my @composites = qw/
-0 1 4 6 8 9 10 12 14 15 16 18 20 21 22
+0 4 6 8 9 10 12 14 15 16 18 20 21 22
9 2047 1373653 25326001 3215031751
561 1105 1729 2465 2821 6601 8911 10585 15841 29341 41041 46657 52633
62745 63973 75361 101101 340561 488881 852841 1857241 6733693
diff --git a/t/81-bignum.t b/t/81-bignum.t
index 3802bb0..de1038e 100644
--- a/t/81-bignum.t
+++ b/t/81-bignum.t
@@ -97,7 +97,7 @@ use Math::Prime::Util qw/
nth_prime_approx
factor
factor_exp
- all_factors
+ divisors
moebius
euler_phi
jordan_totient
@@ -207,7 +207,7 @@ SKIP: {
is_deeply( [factor_exp($n)], [linear_to_exp(@$factors)], "factor_exp($n)" );
}
while (my($n, $allfactors) = each(%allfactors)) {
- is_deeply( [all_factors($n)], $allfactors, "all_factors($n)" );
+ is_deeply( [divisors($n)], $allfactors, "divisors($n)" );
}
}
--
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