[libmath-prime-util-perl] 08/16: Add math functions
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:44:25 UTC 2015
This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.06
in repository libmath-prime-util-perl.
commit 92269f671fa20a3e7a0329e4c24da822a858d25b
Author: Dana Jacobsen <dana at acm.org>
Date: Sun Jun 10 04:31:44 2012 -0600
Add math functions
---
Changes | 4 +-
README | 2 +-
XS.xs | 9 ++
examples/bench-pcapprox.pl | 37 +++++++
examples/test-pcapprox.pl | 72 +++++++++++++
lib/Math/Prime/Util.pm | 51 ++++++++-
t/18-functions.t | 88 +++++++++++++++
util.c | 259 ++++++++++++++++++++++++++++++++++++++++++++-
util.h | 4 +
9 files changed, 514 insertions(+), 12 deletions(-)
diff --git a/Changes b/Changes
index eaf0d16..1990528 100644
--- a/Changes
+++ b/Changes
@@ -3,8 +3,10 @@ Revision history for Perl extension Math::Prime::Util.
0.05 8 June 2012
- Use assembler for mulmod if 64-bit and gcc and x86
- nth_prime routines should now all croak on overflow in the same way.
- - Segmented prime_count, things like this are efficient:
+ - Segmented prime_count, things like this are reasonably efficient:
say prime_count( 10**16, 10**16 + 2**20 )
+ - Add Ei(x), li(x), and R(x) functions.
+ - prime_count_approx uses R(x), making it vastly more accurate.
0.04 7 June 2012
- Didn't do tests on 32-bit machine before release. Test suite caught
diff --git a/README b/README
index bffb345..36ed54e 100644
--- a/README
+++ b/README
@@ -1,4 +1,4 @@
-Math::Prime::Util version 0.04
+Math::Prime::Util version 0.05
A set of utilities related to prime numbers. These include multiple sieving
methods, is_prime, prime_count, nth_prime, approximations and bounds for
diff --git a/XS.xs b/XS.xs
index a73f3c2..e12bbe5 100644
--- a/XS.xs
+++ b/XS.xs
@@ -388,3 +388,12 @@ miller_rabin(IN UV n, ...)
int
is_prob_prime(IN UV n)
+
+double
+ExponentialIntegral(double x)
+
+double
+LogarithmicIntegral(double x)
+
+double
+RiemannR(double x)
diff --git a/examples/bench-pcapprox.pl b/examples/bench-pcapprox.pl
new file mode 100755
index 0000000..576e275
--- /dev/null
+++ b/examples/bench-pcapprox.pl
@@ -0,0 +1,37 @@
+#!/usr/bin/env perl
+use strict;
+use warnings;
+use Math::Prime::Util ":all";
+use Benchmark qw/:all/;
+use List::Util qw/min max/;
+my $maxdigits = (~0 <= 4294967295) ? 10 : 20;
+
+my $count = shift || -5;
+
+srand(29);
+test_at_digits($_) for (5 .. $maxdigits);
+
+
+sub test_at_digits {
+ my $digits = shift;
+ die "Digits must be > 0" unless $digits > 0;
+
+ my $base = ($digits == 1) ? 0 : int(10 ** ($digits-1));
+ my $max = int(10 ** $digits);
+ $max = ~0 if $max > ~0;
+ my @nums = map { $base+int(rand($max-$base)) } (1 .. 1000);
+ my $min_num = min @nums;
+ my $max_num = max @nums;
+
+ #print "miller_rabin for 1000 random $digits-digit numbers ($min_num - $max_num)\n";
+
+ my $sum;
+ cmpthese($count,{
+ 'lower' => sub { $sum += prime_count_lower($_) for @nums },
+ 'luapprox' => sub { $sum += (prime_count_lower($_)+prime_count_upper($_))/2 for @nums },
+ 'approx' => sub { $sum += prime_count_approx($_) for @nums },
+ 'li' => sub { $sum += LogarithmicIntegral($_) for @nums },
+ 'R' => sub { $sum += RiemannR($_) for @nums },
+ });
+ print "\n";
+}
diff --git a/examples/test-pcapprox.pl b/examples/test-pcapprox.pl
new file mode 100755
index 0000000..70b1041
--- /dev/null
+++ b/examples/test-pcapprox.pl
@@ -0,0 +1,72 @@
+#!/usr/bin/env perl
+use strict;
+use warnings;
+use Math::Prime::Util qw/prime_count prime_count_approx prime_count_lower prime_count_upper LogarithmicIntegral RiemannR/;
+$| = 1; # fast pipes
+my $use64 = Math::Prime::Util::_maxbits > 32;
+
+
+my %pivals = (
+ 10 => 4,
+ 100 => 25,
+ 1000 => 168,
+ 10000 => 1229,
+ 100000 => 9592,
+ 1000000 => 78498,
+ 10000000 => 664579,
+ 100000000 => 5761455,
+ 1000000000 => 50847534,
+ 10000000000 => 455052511,
+ 100000000000 => 4118054813,
+ 1000000000000 => 37607912018,
+ 10000000000000 => 346065536839,
+ 100000000000000 => 3204941750802,
+ 1000000000000000 => 29844570422669,
+ 10000000000000000 => 279238341033925,
+ 100000000000000000 => 2623557157654233,
+ 1000000000000000000 => 24739954287740860,
+10000000000000000000 => 234057667276344607,
+);
+
+printf(" N %12s %12s %12s\n", "pc_approx", "Li", "R");
+printf("----- %12s %12s %12s\n", '-'x12,'-'x12,'-'x12);
+foreach my $n (sort {$a<=>$b} keys %pivals) {
+ my $pin = $pivals{$n};
+ my $pca = prime_count_approx($n);
+
+ my $pcli = ($n < 2) ? 0 : int(LogarithmicIntegral($n)-LogarithmicIntegral(2)+0.5);
+
+ my $r = int(RiemannR($n)+0.5);
+
+ printf "10^%2d %12d %12d %12d\n", length($n)-1,
+ abs($pca-$pin), abs($pcli-$pin), abs($r-$pin);
+}
+
+# Also see http://empslocal.ex.ac.uk/people/staff/mrwatkin/zeta/encoding1.htm
+# for some ideas one how this could be made even more accurate.
+
+print "\n";
+print "Lower / Upper bounds. Percentages.\n";
+print "\n";
+
+printf(" N %12s %12s %12s %12s\n", "lower", "upper", "SchoenfeldL", "SchoenfeldU");
+printf("----- %12s %12s %12s %12s\n", '-'x12,'-'x12,'-'x12,'-'x12);
+foreach my $n (sort {$a<=>$b} keys %pivals) {
+ my $pin = $pivals{$n};
+ my $pcl = prime_count_lower($n);
+ my $pcu = prime_count_upper($n);
+ my ($scl,$scu) = schoenfeld($n);
+
+ #printf "10^%2d %12d %12d\n", length($n)-1, $pin-$pcl, $pcu-$pin;
+ printf "10^%2d %12.7f %12.7f %12.7f %12.7f\n",
+ length($n)-1, 100*($pin-$pcl)/$pin, 100*($pcu-$pin)/$pin,
+ 100*($pin-$scl)/$pin, 100*($scu-$pin)/$pin;
+}
+
+sub schoenfeld {
+ my $x = shift;
+ my $lix = LogarithmicIntegral($x);
+ my $bound = ((sqrt($x)*log($x)) / 8*3.1415926535);
+ ($lix-$bound,$lix+$bound);
+}
+
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index 217f8a4..0ca736a 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -5,7 +5,7 @@ use Carp qw/croak confess/;
BEGIN {
$Math::Prime::Util::AUTHORITY = 'cpan:DANAJ';
- $Math::Prime::Util::VERSION = '0.04';
+ $Math::Prime::Util::VERSION = '0.05';
}
# parent is cleaner, and in the Perl 5.10.1 / 5.12.0 core, but not earlier.
@@ -20,6 +20,7 @@ our @EXPORT_OK = qw(
nth_prime nth_prime_lower nth_prime_upper nth_prime_approx
random_prime random_ndigit_prime
factor
+ ExponentialIntegral LogarithmicIntegral RiemannR
);
our %EXPORT_TAGS = (all => [ @EXPORT_OK ]);
@@ -387,10 +388,9 @@ above that range.
" primes below one quintillion.\n";
Returns an approximation to the C<prime_count> function, without having to
-generate any primes. The results are very close for small numbers, but less
-so with large ranges. The current implementation is 0.00033% too small
-for the example, but takes under a microsecond and no memory to get the
-result.
+generate any primes. The current implementation uses the Riemann R function
+which is quite accurate. A slightly faster answer, but less accurate,
+can be obtained by averaging the upper and lower bounds.
=head2 nth_prime
@@ -618,6 +618,46 @@ input value. On some inputs they will take a very long time, while on
others they succeed in a remarkably short time.
+
+=head1 MATHEMATICAL FUNCTIONS
+
+=head2 ExponentialIntegral
+
+ my $Ei = ExponentialIntegral($x);
+
+Given a non-zero floating point input C<x>, this returns the exponential integral
+of C<x>, defined as the integral of C<e^t/t dt> from C<-infinity> to C<x>.
+Depending on the input, the integral
+is calculated using continued fractions (negative C<x>), a convergent series
+(small positive C<x>), or an asymptotic divergent series (large positive C<x>).
+This function only considers the real part.
+
+=head2 LogarithmicIntegral
+
+ my $li = LogarithmicIntegral($x)
+
+Given a positive floating point input, returns the floating point logarithmic
+integral of C<x>, defined as the integral of C<dt/ln t> from C<0> to C<x>.
+If given a negative input, the function will croak. The function returns
+0 at C<x = 0>, and C<-infinity> at C<x = 1>.
+
+This is often known as C<li(x)>. A related function is the offset logarithmic
+integral, sometimes known as C<Li(x)> which avoids the singularity at 1. It
+may be defined as C<Li(x) = li(x) - li(2)>.
+
+This function is implemented as C<li(x) = Ei(ln x)> after handling special
+values.
+
+=head2 RiemannR
+
+ my $r = RiemannR($x);
+
+Given a positive non-zero floating point input, returns the floating
+point value of Riemann's R function. Riemann's R function gives a very close
+approximation to the prime counting function.
+
+
+
=head1 LIMITATIONS
I have not completed testing all the functions near the word size limit
@@ -631,6 +671,7 @@ try to determine if your Perl is broken. This will show up in factoring tests.
Perl 5.6.2 32-bit works fine, as do later versions with 32-bit and 64-bit.
+
=head1 PERFORMANCE
Counting the primes to C<10^10> (10 billion), with time in seconds.
diff --git a/t/18-functions.t b/t/18-functions.t
new file mode 100644
index 0000000..1ed634a
--- /dev/null
+++ b/t/18-functions.t
@@ -0,0 +1,88 @@
+#!/usr/bin/env perl
+use strict;
+use warnings;
+
+use Test::More;
+use Math::Prime::Util qw/prime_count ExponentialIntegral LogarithmicIntegral RiemannR/;
+
+my $use64 = Math::Prime::Util::_maxbits > 32;
+my $extra = defined $ENV{RELEASE_TESTING} && $ENV{RELEASE_TESTING};
+
+#plan tests => 14*3 + ($extra ? 14 : 8) + ($use64 ? 2*18 : 0);
+
+eval { ExponentialIntegral(0); };
+like($@, qr/invalid/i, "Ei(0) is invalid");
+eval { LogarithmicIntegral(-1); };
+like($@, qr/invalid/i, "li(-1) is invalid");
+eval { RiemannR(0); };
+like($@, qr/invalid/i, "R(0) is invalid");
+eval { RiemannR(-1); };
+like($@, qr/invalid/i, "R(-1) is invalid");
+
+cmp_ok( LogarithmicIntegral(1), '<=', 0 - (~0 * ~0), "li(1) is -infinity" );
+
+
+my %eivals = (
+ -10 => -.00000415696892968532438,
+ -0.1 => -1.8229239584193906660809,
+ 0.693147180559945 => 1.0451637801174927848446, # log2
+ 1 => 1.8951178163559367554665,
+ 1.5 => 3.3012854491297978379574,
+ 2 => 4.9542343560018901633795,
+ 5 => 40.185275355803177455091,
+ 10 => 2492.2289762418777591384,
+ 12 => 14959.532666397528852292,
+ 20 => 25615652.664056588820481,
+ 40 => 6039718263611241.5783592,
+ 41 => 16006649143245041.110700,
+);
+
+while (my($n, $ein) = each (%eivals)) {
+ cmp_closeto( ExponentialIntegral($n), $ein, 0.00000001 * abs($ein), "Ei($n) ~= $ein");
+}
+
+# In pari these are: -eint1(-log($n))
+my %livals = (
+ 0 => 0,
+ 1.01 => -4.0229586739299358695031,
+ 2 => 1.0451637801174927848446,
+ 10 => 6.1655995047872979375230,
+ 24 => 11.200315795232698830550,
+ 1000 => 177.60965799015222668764,
+ 100000 => 9629.8090010507982050343,
+ 100000000 => 5762209.3754480314675691,
+ 4294967295 => 203284081.95454158906409,
+ 10000000000 => 455055614.58662307560953,
+ 100000000000 => 4118066400.6216115150394,
+);
+
+while (my($n, $lin) = each (%livals)) {
+ cmp_closeto( LogarithmicIntegral($n), $lin, 0.00000001 * abs($lin), "li($n) ~= $lin");
+}
+
+# Values from T. R. Nicely for comparison
+my %rvals = (
+ 1.01 => 1.0060697180622924796117,
+ 2 => 1.5410090161871318832885,
+ 10 => 4.5645831410050902398658,
+ 1000 => 168.35944628116734806491,
+ 1000000 => 78527.399429127704858870,
+ 10000000 => 664667.44756474776798535,
+ 4294967295 => 203280697.51326064541983,
+ 10000000000 => 455050683.30684692446315,
+18446744073709551615 => 4.25656284014012122706963685602e17,
+);
+while (my($n, $rin) = each (%rvals)) {
+ cmp_closeto( RiemannR($n), $rin, 0.00000001 * abs($rin), "R($n) ~= $rin");
+}
+
+
+done_testing();
+
+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/util.c b/util.c
index c924322..c469f05 100644
--- a/util.c
+++ b/util.c
@@ -377,13 +377,30 @@ UV prime_count_upper(UV x)
UV prime_count_approx(UV x)
{
- /* Placeholder for fancy algorithms, like Tomás Oliveira e Silva's:
+ /*
+ * A simple way:
+ * return ((prime_count_lower(x) + prime_count_upper(x)) / 2);
+ * With the current bounds, this is ~131k at 10^10 and 436B at 10^19.
+ *
+ * The logarithmic integral works quite well, with absolute errors of
+ * ~3100 at 10^10 and ~100M at 10^19.
+ *
+ * Riemann's R function works astoundingly well, with errors of ~1828
+ * at 10^10 and 24M at 10^19.
+ *
+ * Getting fancier, one try using Riemann's pi formula:
* http://trac.sagemath.org/sage_trac/ticket/8135
- * or an approximation to Li(x) plus a delta.
*/
- UV lower = prime_count_lower(x);
- UV upper = prime_count_upper(x);
- return ((lower + upper) / 2);
+ double R;
+ if (x < NPRIME_COUNT_SMALL)
+ return prime_count_small[x];
+
+ R = RiemannR(x);
+ /* We could add the additional factor:
+ * R = R - (1.0 / log(x)) + (M_1_PI * atan(M_PI/log(x)))
+ * but it's extraordinarily small, so not worth calculating here.
+ */
+ return (UV)(R+0.5);
}
@@ -762,3 +779,235 @@ UV prime_count_seg(UV low, UV high)
return count;
}
+
+/*
+ * See:
+ * "Multiple-Precision Exponential Integral and Related Functions"
+ * by David M. Smith, and
+ * "On the Evaluation of the Complex-Valued Exponential Integral"
+ * by Vincent Pegoraro and Philipp Slusallek
+ * "Numerical Recipes" 3rd edition
+ * by William H. Press et al.
+ *
+ * Any mistakes here are completely my fault. This code has not been
+ * verified for anything serious. For better reulsts, see:
+ * http://www.trnicely.net/pi/pix_0000.htm
+ * which although the author claims are demonstration programs, will
+ * produce more usable results than this code does.
+ */
+
+static double const euler_mascheroni = 0.57721566490153286060651209008240243104215933593992;
+static double const li2 = 1.045163780117492784844588889194613136522615578151;
+
+double ExponentialIntegral(double x) {
+ double const tol = 0.0000000001;
+ double val, term, x_to_the_n, fact_n;
+ double y, t;
+ double sum = 0.0;
+ double c = 0.0;
+ int n;
+
+ if (x == 0) croak("Invalid input to ExponentialIntegral: x must be != 0");
+
+ if (x < 0) {
+ /* continued fraction */
+#if 0
+ val = 0;
+ for (n = 50; n >= 1; n--) {
+ val = -n * n / (2 * n + 1 - x + val);
+ }
+ val = -exp(x) / (1 - x + val);
+#else
+ double old;
+ double lc = 0;
+ double ld = 1 / (1 - x);
+ val = ld * (-exp(x));
+ for (n = 1; n <= 2000; n++) {
+ lc = 1 / (2*n + 1 - x - n*n*lc);
+ ld = 1 / (2*n + 1 - x - n*n*ld);
+ old = val;
+ val *= ld/lc;
+ if ( fabs(val-old) <= tol*fabs(val) )
+ break;
+ }
+#endif
+ } else if (x < -log(tol)) {
+ /* Convergent series */
+ fact_n = 1;
+
+ y = euler_mascheroni-c; t = sum+y; c = (t-sum)-y; sum = t;
+ y = log(x)-c; t = sum+y; c = (t-sum)-y; sum = t;
+
+ for (n = 1; n <= 100; n++) {
+ fact_n *= x/n;
+ term = fact_n/n;
+ y = term-c; t = sum+y; c = (t-sum)-y; sum = t;
+ /* printf("C after adding %.8lf, val = %.8lf\n", term, sum); */
+ if (term < tol) break;
+ }
+ val = sum;
+ if (term > tol) printf("Tolerance not met after %d rounds for Ei(%lf)\n", n, x);
+ } else {
+ /* Asymptotic divergent series */
+#if 0
+ double last_term, sum;
+ last_term = 1.0;
+ x_to_the_n = 1;
+ fact_n = 1;
+
+ y = 1.0-c; t = sum+y; c = (t-sum)-y; sum = t;
+
+ for (n = 1; n <= 100; n++) {
+ x_to_the_n *= x;
+ fact_n *= n;
+ term = fact_n / x_to_the_n;
+ if (term > last_term) break;
+ last_term = term;
+ y = term-c; t = sum+y; c = (t-sum)-y; sum = t;
+ /* printf("A after adding %.8lf, sum = %.8lf\n", term, sum); */
+ }
+ val = (exp(x) / x) * sum - euler_mascheroni;
+#else
+ double last_term, sum;
+
+ val = exp(x) / x;
+ term = 1.0;
+ y = 1.0-c; t = sum+y; c = (t-sum)-y; sum = t;
+
+ for (n = 1; n <= 100; n++) {
+ last_term = term;
+ term *= n/x;
+ if (term < (tol/val)) break;
+ if (term < last_term) {
+ y = term-c; t = sum+y; c = (t-sum)-y; sum = t;
+ /* printf("A after adding %.8lf, sum = %.8lf\n", term, sum); */
+ } else {
+ y = (-last_term/3)-c; t = sum+y; c = (t-sum)-y; sum = t;
+ /* printf("A after adding %.8lf, sum = %.8lf\n", -last_term/3, sum); */
+ break;
+ }
+ }
+ val *= sum;
+#endif
+ }
+
+ return val;
+}
+
+double LogarithmicIntegral(double x) {
+ if (x == 0) return 0;
+ if (x == 1) return -INFINITY;
+ if (x == 2) return li2;
+ if (x <= 0) croak("Invalid input to ExponentialIntegral: x must be > 0");
+ return ExponentialIntegral(log(x));
+}
+
+/*
+ * Storing the first 10-20 Zeta values makes sense. Past that it is purely
+ * to avoid making the call to generate them ourselves. We could cache the
+ * calculated values. */
+static const double riemann_zeta_table[] = {
+ 1.6449340668482264364724151666460251892, /* zeta(2) */
+ 1.2020569031595942853997381615114499908,
+ 1.0823232337111381915160036965411679028,
+ 1.0369277551433699263313654864570341681,
+ 1.0173430619844491397145179297909205279,
+ 1.0083492773819228268397975498497967596,
+ 1.0040773561979443393786852385086524653,
+ 1.0020083928260822144178527692324120605,
+ 1.0009945751278180853371459589003190170,
+ 1.0004941886041194645587022825264699365,
+ 1.0002460865533080482986379980477396710,
+ 1.0001227133475784891467518365263573957,
+ 1.0000612481350587048292585451051353337,
+ 1.0000305882363070204935517285106450626,
+ 1.0000152822594086518717325714876367220,
+ 1.0000076371976378997622736002935630292, /* zeta(17) Past here all we're */
+ 1.0000038172932649998398564616446219397, /* zeta(18) getting is speed. */
+ 1.0000019082127165539389256569577951013,
+ 1.0000009539620338727961131520386834493,
+ 1.0000004769329867878064631167196043730,
+ 1.0000002384505027277329900036481867530,
+ 1.0000001192199259653110730677887188823,
+ 1.0000000596081890512594796124402079358,
+ 1.0000000298035035146522801860637050694,
+ 1.0000000149015548283650412346585066307,
+ 1.0000000074507117898354294919810041706,
+ 1.0000000037253340247884570548192040184,
+ 1.0000000018626597235130490064039099454,
+ 1.0000000009313274324196681828717647350,
+ 1.0000000004656629065033784072989233251,
+ 1.0000000002328311833676505492001455976,
+ 1.0000000001164155017270051977592973835,
+ 1.0000000000582077208790270088924368599,
+ 1.0000000000291038504449709968692942523,
+ 1.0000000000145519218910419842359296322,
+ 1.0000000000072759598350574810145208690,
+ 1.0000000000036379795473786511902372363,
+ 1.0000000000018189896503070659475848321,
+ 1.0000000000009094947840263889282533118,
+};
+#define NPRECALC_ZETA (sizeof(riemann_zeta_table)/sizeof(riemann_zeta_table[0]))
+
+static double evaluate_zeta(double x) {
+ double const tol = 1e-14;
+ double y, t;
+ double sum = 0.0;
+ double c = 0.0;
+ double term;
+ int k;
+
+ /* Simple method. Slow and inaccurate near x=1, but iterations and accuracy
+ * improve quickly.
+ */
+
+ term = 1.0; y = term-c; t = sum+y; c = (t-sum)-y; sum = t;
+ term = 1.0/exp2(x); y = term-c; t = sum+y; c = (t-sum)-y; sum = t;
+
+ for (k = 3; k <= 100000; k++) {
+ if (term < tol) break;
+ if (k == 4) term = 1.0 / exp2(2*x);
+ else if (k == 8) term = 1.0 / exp2(4*x);
+ else if (k == 16) term = 1.0 / exp2(8*x);
+ else term = 1.0 / pow(k, x);
+ y = term-c; t = sum+y; c = (t-sum)-y; sum = t;
+ }
+ return sum;
+}
+
+double RiemannR(double x) {
+ double const tol = 1e-14;
+ double y, t, part_term, term, flogx, zeta;
+ double sum = 0.0;
+ double c = 0.0;
+ int k;
+
+ if (x <= 0) croak("Invalid input to ReimannR: x must be > 0");
+
+ y = 1.0-c; t = sum+y; c = (t-sum)-y; sum = t;
+
+ flogx = log(x);
+ part_term = 1;
+
+ /* Do small k with zetas from table */
+ for (k = 1; k <= NPRECALC_ZETA; k++) {
+ zeta = riemann_zeta_table[k+1-2];
+ part_term *= flogx / k;
+ term = part_term / (k * zeta);
+ y = term-c; t = sum+y; c = (t-sum)-y; sum = t;
+ if (term < tol) break;
+ /* printf("R after adding %.8lf, sum = %.8lf\n", term, sum); */
+ }
+ /* Finish with function */
+ for (; k <= 10000; k++) {
+ if (term < tol) break;
+ zeta = evaluate_zeta(k+1);
+ part_term *= flogx / k;
+ term = part_term / (k * zeta);
+ y = term-c; t = sum+y; c = (t-sum)-y; sum = t;
+ /* printf("R after adding %.8lf, sum = %.8lf\n", term, sum); */
+ }
+
+
+ return sum;
+}
diff --git a/util.h b/util.h
index 6a09d98..cb7ff43 100644
--- a/util.h
+++ b/util.h
@@ -20,6 +20,10 @@ extern UV nth_prime_upper(UV x);
extern UV nth_prime_approx(UV x);
extern UV nth_prime(UV x);
+extern double ExponentialIntegral(double x);
+extern double LogarithmicIntegral(double x);
+extern double RiemannR(double x);
+
#define SEGMENT_CHUNK_SIZE UVCONST(262144)
extern unsigned char* get_prime_segment(void);
extern void free_prime_segment(void);
--
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