[libmath-prime-util-perl] 10/20: Add Chebyshev theta and psi functions
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:47:31 UTC 2015
This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.23
in repository libmath-prime-util-perl.
commit 1e807e10ac84cd86445a2b5eb7e24a5e811f47f9
Author: Dana Jacobsen <dana at acm.org>
Date: Fri Mar 1 16:09:43 2013 -0800
Add Chebyshev theta and psi functions
---
Changes | 9 ++--
TODO | 10 ++--
XS.xs | 6 +++
lib/Math/Prime/Util.pm | 141 +++++++++++++++++++++++++++++++++----------------
t/19-moebius.t | 46 +++++++++++++++-
t/21-conseq-lcm.t | 2 +-
util.c | 75 ++++++++++++++++++++++----
util.h | 2 +
8 files changed, 224 insertions(+), 67 deletions(-)
diff --git a/Changes b/Changes
index 945d23a..e52cf01 100644
--- a/Changes
+++ b/Changes
@@ -1,6 +1,6 @@
Revision history for Perl extension Math::Prime::Util.
-0.23 xx February 2013
+0.23 xx March 2013
- exp_mangoldt in XS, and speed up the PP version.
@@ -11,14 +11,15 @@ Revision history for Perl extension Math::Prime::Util.
- PP Zeta code replaced (for no-MPFR, non-bignums) with new series. The
new code is much more accurate for small values, and *much* faster.
- - Add consecutive_integer_lcm function, just like MPU::GMP's.
-
- - Add examples for first and second Chebyshev functions.
+ - Add consecutive_integer_lcm function, just like MPU::GMP's (though we
+ define ci_lcm(0) = 0, which should get propogated).
- Implement binary search on RiemannR for XS nth_prime when n > 2e11.
Runs ~2x faster for 1e12, 3x faster for 1e13. Thanks to Programming
Praxis for the idea and motivation.
+ - Add the first and second Chebyshev functions (theta and psi).
+
0.22 26 February 2013
- Move main factor loop out of xs and into factor.c.
diff --git a/TODO b/TODO
index 31fa191..cb3763b 100644
--- a/TODO
+++ b/TODO
@@ -37,9 +37,13 @@
- More efficient Mertens. The current version has poor growth.
-- Add first and second Chebyshev functions. See Planat and Solé (2011).
- MPU::GMP has primorial and consecutive_integer_lcm which are what we need.
-
- Keep speeding up factoring for very large values.
- More efficient totient segment. Do we really need primes to n/2?
+
+- Add a more complex system for allowing something like:
+ FOREACH_PRIME(low, high) {
+ ....
+ } END_FOREACH_PRIME
+ that (1) handles any low / high, and (2) segments. Take segment_primes
+ from XS.xs to do this.
diff --git a/XS.xs b/XS.xs
index 95e4733..83fc0ce 100644
--- a/XS.xs
+++ b/XS.xs
@@ -463,3 +463,9 @@ _XS_exp_mangoldt(IN UV n)
RETVAL = 1;
OUTPUT:
RETVAL
+
+double
+_XS_chebyshev_theta(IN UV n)
+
+double
+_XS_chebyshev_psi(IN UV n)
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index c87d125..cc7129e 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -12,25 +12,27 @@ BEGIN {
# parent is cleaner, and in the Perl 5.10.1 / 5.12.0 core, but not earlier.
# use parent qw( Exporter );
use base qw( Exporter );
-our @EXPORT_OK = qw(
- prime_get_config prime_set_config
- prime_precalc prime_memfree
- is_prime is_prob_prime is_provable_prime
- is_strong_pseudoprime is_strong_lucas_pseudoprime
- is_aks_prime
- miller_rabin
- primes
- next_prime prev_prime
- prime_count prime_count_lower prime_count_upper prime_count_approx
- nth_prime nth_prime_lower nth_prime_upper nth_prime_approx
- random_prime random_ndigit_prime random_nbit_prime
- random_strong_prime random_maurer_prime
- primorial pn_primorial consecutive_integer_lcm
- factor all_factors
- moebius mertens euler_phi jordan_totient exp_mangoldt
- divisor_sum
- ExponentialIntegral LogarithmicIntegral RiemannZeta RiemannR
- );
+our @EXPORT_OK =
+ qw( prime_get_config prime_set_config
+ prime_precalc prime_memfree
+ is_prime is_prob_prime is_provable_prime
+ is_strong_pseudoprime is_strong_lucas_pseudoprime
+ is_aks_prime
+ miller_rabin
+ primes
+ next_prime prev_prime
+ prime_count
+ prime_count_lower prime_count_upper prime_count_approx
+ nth_prime nth_prime_lower nth_prime_upper nth_prime_approx
+ random_prime random_ndigit_prime random_nbit_prime
+ random_strong_prime random_maurer_prime
+ primorial pn_primorial consecutive_integer_lcm
+ factor all_factors
+ moebius mertens euler_phi jordan_totient exp_mangoldt
+ chebyshev_theta chebyshev_psi
+ divisor_sum
+ ExponentialIntegral LogarithmicIntegral RiemannZeta RiemannR
+ );
our %EXPORT_TAGS = (all => [ @EXPORT_OK ]);
my %_Config;
@@ -46,19 +48,21 @@ sub import {
sub _import_nobigint {
$_Config{'nobigint'} = 1;
return unless $_Config{'xs'};
- undef *factor; *factor = \&_XS_factor;
- undef *is_prime; *is_prime = \&_XS_is_prime;
- undef *is_prob_prime; *is_prob_prime = \&_XS_is_prob_prime;
- undef *next_prime; *next_prime = \&_XS_next_prime;
- undef *prev_prime; *prev_prime = \&_XS_prev_prime;
- #undef *prime_count; *prime_count = \&_XS_prime_count;
- undef *nth_prime; *nth_prime = \&_XS_nth_prime;
+ undef *factor; *factor = \&_XS_factor;
+ undef *is_prime; *is_prime = \&_XS_is_prime;
+ undef *is_prob_prime; *is_prob_prime = \&_XS_is_prob_prime;
+ undef *next_prime; *next_prime = \&_XS_next_prime;
+ undef *prev_prime; *prev_prime = \&_XS_prev_prime;
+ #undef *prime_count; *prime_count = \&_XS_prime_count;
+ undef *nth_prime; *nth_prime = \&_XS_nth_prime;
undef *is_strong_pseudoprime; *is_strong_pseudoprime = \&_XS_miller_rabin;
- undef *miller_rabin; *miller_rabin = \&_XS_miller_rabin;
- undef *moebius; *moebius = \&_XS_moebius;
- undef *mertens; *mertens = \&_XS_mertens;
- undef *euler_phi; *euler_phi = \&_XS_totient;
- undef *exp_mangoldt; *exp_mangoldt = \&_XS_exp_mangoldt;
+ undef *miller_rabin; *miller_rabin = \&_XS_miller_rabin;
+ undef *moebius; *moebius = \&_XS_moebius;
+ undef *mertens; *mertens = \&_XS_mertens;
+ undef *euler_phi; *euler_phi = \&_XS_totient;
+ undef *exp_mangoldt; *exp_mangoldt = \&_XS_exp_mangoldt;
+ undef *chebyshev_theta; *chebyshev_theta = \&_XS_chebyshev_theta;
+ undef *chebyshev_psi; *chebyshev_psi = \&_XS_chebyshev_psi;
}
BEGIN {
@@ -1038,7 +1042,7 @@ sub pn_primorial {
sub consecutive_integer_lcm {
my($n) = @_;
_validate_positive_integer($n);
- return 1 if $n <= 1;
+ return 0 if $n < 1;
my $pn = 1;
if ($n >= (($_Config{'maxbits'} == 32) ? 22 : 46)) {
@@ -1335,6 +1339,31 @@ sub exp_mangoldt {
return $first;
}
+sub chebyshev_theta {
+ my($n) = @_;
+ _validate_positive_integer($n);
+ return _XS_chebyshev_theta($n) if $n <= $_XS_MAXVAL;
+ my $sum = 0.0;
+ foreach my $p (@{primes($n)}) {
+ $sum += log($p);
+ }
+ return $sum;
+}
+sub chebyshev_psi {
+ my($n) = @_;
+ _validate_positive_integer($n);
+ return 0 if $n <= 1;
+ return _XS_chebyshev_psi($n) if $n <= $_XS_MAXVAL;
+ my ($sum, $logn, $mults_are_one) = (0.0, log($n), 0);
+ foreach my $p (@{primes($n)}) {
+ my $logp = log($p);
+ $mults_are_one = 1 if !$mults_are_one && $p > int($n/$p);
+ $sum += ($mults_are_one) ? $logp : $logp * int($logn/$logp+1e-15);
+ }
+ return $sum;
+}
+
+
#############################################################################
# Front ends to functions.
@@ -2539,6 +2568,40 @@ and 0 otherwise. We return the exponential so all results are integers.
1 otherwise.
+=head2 chebyshev_theta
+
+ say chebyshev_theta(10000);
+
+Returns θ(n), the first Chebyshev function for a non-negative integer input.
+This is the sum of the logarithm of each prime where C<p E<lt>= n>. An
+alternate computation is as the logarithm of n primorial, hence:
+
+ use List::Util qw/sum/;
+ sub c1a { 0+sum( map { log($_) } @{primes(shift)} ) }
+ use Math::BigFloat;
+ sub c1b { Math::BigFloat->new(primorial(shift))->blog }
+
+yield similar results, albeit slower and using more memory.
+
+
+=head2 chebyshev_psi
+
+ say chebyshev_psi(10000);
+
+Returns ψ(n), the second Chebyshev function for a non-negative integer input.
+This is the sum of the logarithm of each prime where C<p^k E<lt>= n> for an
+integer k. An alternate computation is as the summatory Mangoldt function.
+Another alternate computation is as the logarithm of lcm(1,2,...,n).
+Hence these functions:
+
+ use List::Util qw/sum/;
+ sub c2a { 0+sum( map { log(exp_mangoldt($_)) } 1 .. shift ) }
+ use Math::BigFloat;
+ sub c2b { Math::BigFloat->new(consecutive_integer_lcm(shift))->blog }
+
+yield similar results, albeit slower and using more memory.
+
+
=head2 divisor_sum
say "Sum of divisors of $n:", divisor_sum( $n );
@@ -3128,20 +3191,6 @@ Print some primes above 64-bit range:
# Similar code using Pari:
# perl -MMath::Pari=:int,PARI,nextprime -E 'my $start = PARI "100000000000000000000"; my $end = $start+1000; my $p=nextprime($start); while ($p <= $end) { say $p; $p = nextprime($p+1); }'
-The first Chebyshev function θ(x):
-
- use List::Util qw/sum/;
- sub chebyshev1 { 0+sum( map { log($_) } @{primes(shift)} ) }
- # or using primorial and BigFloat:
- sub chebyshev1a { Math::BigFloat->new(primorial(shift))->blog }
-
-The second Chebyshev function ψ(x):
-
- use List::Util qw/sum/;
- sub chebyshev2 { 0+sum( map { log(exp_mangoldt($_)) } 1 .. shift ) }
- # or using cilcm and BigFloat:
- sub chebyshev2a { Math::BigFloat->new(consecutive_integer_lcm(shift))->blog }
-
Project Euler, problem 3 (Largest prime factor):
use Math::Prime::Util qw/factor/;
diff --git a/t/19-moebius.t b/t/19-moebius.t
index 8081ee0..433b54e 100644
--- a/t/19-moebius.t
+++ b/t/19-moebius.t
@@ -4,7 +4,8 @@ use warnings;
use Test::More;
use Math::Prime::Util
- qw/moebius mertens euler_phi jordan_totient divisor_sum exp_mangoldt/;
+ qw/moebius mertens euler_phi jordan_totient divisor_sum exp_mangoldt
+ chebyshev_theta chebyshev_psi/;
my $extra = defined $ENV{RELEASE_TESTING} && $ENV{RELEASE_TESTING};
my $use64 = Math::Prime::Util::prime_get_config->{'maxbits'} > 32;
@@ -133,6 +134,27 @@ my %mangoldt = (
130321 => 19,
);
+my %chebyshev1 = (
+ 0 => 0,
+ 1 => 0,
+ 2 => 0.693147180559945,
+ 3 => 1.79175946922805,
+ 4 => 1.79175946922805,
+ 5 => 3.40119738166216,
+ 243 => 226.593507136467,
+ 1234567 => 1233272.80087825,
+);
+my %chebyshev2 = (
+ 0 => 0,
+ 1 => 0,
+ 2 => 0.693147180559945,
+ 3 => 1.79175946922805,
+ 4 => 2.484906649788,
+ 5 => 4.0943445622221,
+ 243 => 245.274469978683,
+ 1234567 => 1234515.17962833,
+);
+
plan tests => 0 + 1
+ 1 # Small Moebius
@@ -145,7 +167,9 @@ plan tests => 0 + 1
+ 1 # Calculate J5 two different ways
+ 2 * $use64 # Jordan totient example
+ scalar(keys %sigmak)
- + scalar(keys %mangoldt);
+ + scalar(keys %mangoldt)
+ + scalar(keys %chebyshev1)
+ + scalar(keys %chebyshev2);
ok(!eval { moebius(0); }, "moebius(0)");
@@ -229,3 +253,21 @@ while (my($k, $sigmaref) = each (%sigmak)) {
while (my($n, $em) = each (%mangoldt)) {
is( exp_mangoldt($n), $em, "exp_mangoldt($n) == $em" );
}
+
+###### first Chebyshev function
+while (my($n, $c1) = each (%chebyshev1)) {
+ cmp_closeto( chebyshev_theta($n), $c1, 1e-9*abs($n), "chebyshev_theta($n)" );
+}
+###### second Chebyshev function
+while (my($n, $c2) = each (%chebyshev2)) {
+ cmp_closeto( chebyshev_psi($n), $c2, 1e-9*abs($n), "chebyshev_psi($n)" );
+}
+
+
+sub cmp_closeto {
+ my $got = shift;
+ my $expect = shift;
+ my $tolerance = shift;
+ my $message = shift;
+ cmp_ok( abs($got - $expect), '<=', $tolerance, $message );
+}
diff --git a/t/21-conseq-lcm.t b/t/21-conseq-lcm.t
index e262583..63d17d5 100644
--- a/t/21-conseq-lcm.t
+++ b/t/21-conseq-lcm.t
@@ -8,7 +8,7 @@ use Math::Prime::Util qw/consecutive_integer_lcm/;
plan tests => 101 + 1;
my @lcms = qw/
-1
+0
1
2
6
diff --git a/util.c b/util.c
index 3048e69..159d1b7 100644
--- a/util.c
+++ b/util.c
@@ -28,17 +28,33 @@
extern long double expl(long double);
extern long double logl(long double);
extern long double fabsl(long double);
+ extern long double floorl(long double);
#else
#define powl(x, y) (long double) pow( (double) (x), (double) (y) )
#define expl(x) (long double) exp( (double) (x) )
#define logl(x) (long double) log( (double) (x) )
#define fabsl(x) (long double) fabs( (double) (x) )
+ #define floorl(x) (long double) floor( (double) (x) )
#endif
#ifndef INFINITY
#define INFINITY (DBL_MAX + DBL_MAX)
#endif
+#define KAHAN_INIT(s) \
+ long double s ## _y, s ## _t; \
+ long double s ## _c = 0.0; \
+ long double s = 0.0;
+
+#define KAHAN_SUM(s, term) \
+ do { \
+ s ## _y = (term) - s ## _c; \
+ s ## _t = s + s ## _y; \
+ s ## _c = (s ## _t - s) - s ## _y; \
+ s = s ## _t; \
+ } while (0)
+
+
#include "ptypes.h"
#include "util.h"
#include "sieve.h"
@@ -737,6 +753,54 @@ IV _XS_mertens(UV n) {
#endif
}
+double _XS_chebyshev_theta(UV n)
+{
+ const unsigned char* sieve;
+ KAHAN_INIT(sum);
+
+ if (n >= 2) KAHAN_SUM(sum, 0.6931471805599453094172321214581765680755L);
+ if (n >= 3) KAHAN_SUM(sum, 1.0986122886681096913952452369225257046475L);
+ if (n >= 5) KAHAN_SUM(sum, 1.6094379124341003746007593332261876395256L);
+ if (n < 7) return (double) sum;
+
+ if (get_prime_cache(n, &sieve) < n) {
+ release_prime_cache(sieve);
+ croak("Could not generate sieve for %"UVuf, n);
+ }
+ START_DO_FOR_EACH_SIEVE_PRIME(sieve, 7, n) {
+ KAHAN_SUM(sum, logl(p));
+ } END_DO_FOR_EACH_SIEVE_PRIME
+ release_prime_cache(sieve);
+ return (double) sum;
+}
+double _XS_chebyshev_psi(UV n)
+{
+ const unsigned char* sieve;
+ UV prime, mults_are_one;
+ long double logn, logp;
+ KAHAN_INIT(sum);
+
+ logn = logl(n);
+ for (prime = 2; prime <= 5 && prime <= n; prime += prime-1) {
+ logp = logl(prime);
+ KAHAN_SUM(sum, logp * floorl(logn/logp + 1e-15));
+ }
+ if (n < 7) return (double) sum;
+
+ if (get_prime_cache(n, &sieve) < n) {
+ release_prime_cache(sieve);
+ croak("Could not generate sieve for %"UVuf, n);
+ }
+ mults_are_one = 0;
+ START_DO_FOR_EACH_SIEVE_PRIME(sieve, 7, n) {
+ logp = logl(p);
+ if (!mults_are_one && p > (n/p)) mults_are_one = 1;
+ KAHAN_SUM(sum, (mults_are_one) ? logp : logp * floorl(logn/logp + 1e-15));
+ } END_DO_FOR_EACH_SIEVE_PRIME
+ release_prime_cache(sieve);
+ return (double) sum;
+}
+
/*
@@ -762,17 +826,6 @@ IV _XS_mertens(UV n) {
static long double const euler_mascheroni = 0.57721566490153286060651209008240243104215933593992L;
static long double const li2 = 1.045163780117492784844588889194613136522615578151L;
-#define KAHAN_INIT(s) \
- long double s ## _y, s ## _t; \
- long double s ## _c = 0.0; \
- long double s = 0.0;
-
-#define KAHAN_SUM(s, term) \
- s ## _y = term - s ## _c; \
- s ## _t = s + s ## _y; \
- s ## _c = (s ## _t - s) - s ## _y; \
- s = s ## _t;
-
double _XS_ExponentialIntegral(double x) {
long double const tol = 1e-16;
long double val, term;
diff --git a/util.h b/util.h
index 4ab091f..3492469 100644
--- a/util.h
+++ b/util.h
@@ -16,6 +16,8 @@ extern UV _XS_nth_prime(UV x);
extern IV* _moebius_range(UV low, UV high);
extern IV _XS_mertens(UV n);
+extern double _XS_chebyshev_theta(UV n);
+extern double _XS_chebyshev_psi(UV n);
extern double _XS_ExponentialIntegral(double x);
extern double _XS_LogarithmicIntegral(double x);
--
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