mersenneforum.org  

Go Back   mersenneforum.org > Math Stuff > Computer Science & Computational Number Theory

Reply
 
Thread Tools
Old 2018-03-04, 17:08   #1
danaj
 
"Dana Jacobsen"
Feb 2011
Bangkok, TH

32·101 Posts
Default Requestion for fast chebyshev theta function

I am looking for ideas of making a fast implementation of the Chebyshev theta function. Ideally a method substantially better than the straightforward sum-all-the-logs.

We can compute the Chebyshev psi function easily given a theta function so no need to worry about that.

The theta function is the sum of the logs of the primes up to n. Using log(ab) = log(a)+log(b) it can be seen this is also the log of the primorial.

The simple way to do this is to loop through primes, summing up log(p). Assuming finite precision, we can use a sparse table to accelerate it a bit. I'm looking for some ideas on significant improvements. Without tables, using with ~18 digit precision (long double on x86-64), it takes 3.7s for 2^32 and 1 minute for 2^36.

Going all the way to log(p#), even though we have fast methods for primorials and only have one log to do, seems like a bad idea. p# grows far too quickly. However the idea of batching some of these seems like it could help a little. Given p1 and p2 we can add log(p1)+log(p2) or log(p1*p2). The latter is likely faster.

Precision is a consideration. I am initially doing this using a native type (64-, 80-, or 128-bit float). We need to be careful to not lose too much precision in the sum. This also puts a limit on the batching idea, both in the input size and the impact on output precision from the log. Depending on the float type this also limits the function input to ~ 2^53 or 2^64 depending on the float type.

I'm doing this in C, and currently it's ~200x faster than Pari/GP forprime(p=2,n,s+=log(p)) even with sub-18-digit precision. But of course Pari/GP's log() is quite different than libm's, so this isn't a really fair comparison.

I was hoping there would be a great solution akin to prime counting / summing, but that works on completely multiplicative functions, which log(n) is not.
danaj is offline   Reply With Quote
Old 2018-03-04, 21:54   #2
danaj
 
"Dana Jacobsen"
Feb 2011
Bangkok, TH

32×101 Posts
Default

Quote:
Originally Posted by JM Montolio A View Post
¿ 200x faster that PARIGP ?

I believe PARIGP is unbeatable.
Trolling?

Pari is about 10x slower looping through the primes -- some of this is just the obvious C vs. GP. primesieve is probably even faster. The real kicker is logl or logq (quadmath) vs. Pari's log. If I wanted more precision I'd have to use MPFR or my own GMP/mpf version (which is slower than MPFR), and I bet most of the difference would disappear.

Regardless, I can microoptimize loops and so on, but if there is a better *algorithm* it could easily win. Similar to how even a so-so Lehmer's prime count method can run much faster than counting primes with the excellent primesieve. The latter is seriously optimized for its task, but there are better ways to count primes than generating them.
danaj is offline   Reply With Quote
Old 2018-03-05, 00:27   #3
CRGreathouse
 
CRGreathouse's Avatar
 
Aug 2006

3·1,993 Posts
Default

I don't know of a clever algorithm. Batched logs (ideally with some guard bits) are the best I could come up with, too.
CRGreathouse is offline   Reply With Quote
Old 2018-03-05, 17:07   #4
Dr Sardonicus
 
Dr Sardonicus's Avatar
 
Feb 2017
Nowhere

5×13×79 Posts
Default

Offhand, I can't think of anything faster than "batched logs and guard bits." As you say, precision is key. Alas, Chebyshev psi seems to be the easier one to work with.

There are any number of tantalizing results: It is of course well known that either Chebyshev theta (sum of log p) or Chebyshev psi (sum of log(p^e)) is asymptotic to the upper bound. There is an "explicit formula" relating psi to the zeroes of the Riemann \zeta-function.

There is also a RH-dependent error term

\psi(x) = x + O(\sqrt{x} log^2(x))

which could conceivably help. Of course, you would have to formulate the difference between Chebyshev psi and Chebyshev theta suitably.
Dr Sardonicus is offline   Reply With Quote
Old 2018-03-05, 18:30   #5
mart_r
 
mart_r's Avatar
 
Dec 2008
you know...around...

3×229 Posts
Default

I think I came up with something original on this issue in 2011, but got stuck somewhere in-between: See posts 1 - 4 in
http://www.mersenneforum.org/showthread.php?t=15107
mart_r is offline   Reply With Quote
Old 2018-03-06, 14:22   #6
danaj
 
"Dana Jacobsen"
Feb 2011
Bangkok, TH

32×101 Posts
Default

One thing I found useful with precision was using Kahan summation. While it slows it down about 20%, it adds at least 3 correct digits at 2^32. With Pari or MPFR it's cheaper and easier to just add more bits.

I tested using MPFR and it is a massive slowdown, and the times are quite close to Pari/GP. There is little time impact of adding some more digits. I did not do a timing comparison with quadmath, but it's quite noticible slower than gcc/Intel long doubles.
danaj is offline   Reply With Quote
Old 2018-03-06, 17:04   #7
axn
 
axn's Avatar
 
Jun 2003

5,197 Posts
Default

Quote:
Originally Posted by danaj View Post
IGoing all the way to log(p#), even though we have fast methods for primorials and only have one log to do, seems like a bad idea. p# grows far too quickly. However the idea of batching some of these seems like it could help a little. Given p1 and p2 we can add log(p1)+log(p2) or log(p1*p2). The latter is likely faster.
Code:
s=1.0; forprime(p=2,2^26, s*=p); log(s)
is much faster and more accurate than
Code:
s=0.0; forprime(p=2,2^26, s+=log(p)); s

Last fiddled with by axn on 2018-03-06 at 17:10 Reason: code correction
axn is offline   Reply With Quote
Old 2018-03-06, 20:47   #8
R. Gerbicz
 
R. Gerbicz's Avatar
 
"Robert Gerbicz"
Oct 2005
Hungary

27578 Posts
Default

Quote:
Originally Posted by danaj View Post
I am looking for ideas of making a fast implementation of the Chebyshev theta function. Ideally a method substantially better than the straightforward sum-all-the-logs.
One obvious method would be to use (large) precomputed table for sum(log(p)), for example with spacing=2^20 up to n=2^40.

Even without tables, there is a way to gain a speedup: use the central binomial coefficient:
binomial(2*n,n) is divisible by all primes between (n,2*n), and has got exponent<=1 for all primes in (sqrt(2*n),n], and we know that this exponent is constant in "long" intervals, say it is 1 in (n/2,2*n/3]. For p<=sqrt(2*n) compute the exponent for each prime and finally use a good Stirling approx. to get log(binomial(2*n,n))=log((2*n)!)-2*log(n!).

[Paul Erdos used exactly this idea to prove the Chebyshev's theorem]
R. Gerbicz is offline   Reply With Quote
Old 2018-03-08, 13:33   #9
Dr Sardonicus
 
Dr Sardonicus's Avatar
 
Feb 2017
Nowhere

5·13·79 Posts
Default

Quote:
Originally Posted by R. Gerbicz View Post
One obvious method would be to use (large) precomputed table for sum(log(p)), for example with spacing=2^20 up to n=2^40.

Even without tables, there is a way to gain a speedup: use the central binomial coefficient:
binomial(2*n,n) is divisible by all primes between (n,2*n), and has got exponent<=1 for all primes in (sqrt(2*n),n], and we know that this exponent is constant in "long" intervals, say it is 1 in (n/2,2*n/3]. For p<=sqrt(2*n) compute the exponent for each prime and finally use a good Stirling approx. to get log(binomial(2*n,n))=log((2*n)!)-2*log(n!).

[Paul Erdos used exactly this idea to prove the Chebyshev's theorem]
I dimly recall having seen this approach presented in an Analytic Number Theory course many years ago. It is a sweet piece of reasoning.

In general, for the prime p, the exponent of the p-power factor of the binomial coefficient binomial(N, m) is the number of carries in the base-p addition

m + (N - m) = N.

For binomial(2*n, n) this is the number of carries in the base-p addition n + n = 2*n.

This result makes it easy to confirm the statements about exponents.

It also makes the following an easy exercise: The binomial coefficients binomial(N, m), m = 0 to N, are all odd if and only if N is one less than a power of two.
Dr Sardonicus is offline   Reply With Quote
Old 2018-03-31, 14:59   #10
mart_r
 
mart_r's Avatar
 
Dec 2008
you know...around...

3·229 Posts
Default

I'm kind of intrigued by this ansatz with the binomial numbers.
I see why we get the prime numbers in intervals of (n/(2k-1), n/2k] for k>=1 in the factorization of \frac{n!}{(\frac{n}{2})!^2}, but is it possible to "fill the gaps"?
For instance, with \frac{n!}{(\frac{n}{2})!(\frac{n}{3})!(\frac{n}{6})!}, we have most of the primes except those in intervals (n/6k, n/(6k+1)], and the prime factors in intervals (n/(6k-1), n/6k] have exponent 2.
With \frac{n!}{(\frac{n}{2})!(\frac{n}{4})!(\frac{n}{8})!(\frac{n}{16})!...}, all intervals are covered, but the exponent in (n/k, n/(k+1)] is equal to the "hamming weight", i.e. the number of 1's in the binary representation of k.
Any further ideas there?

And what about x-\theta(x) ~ log(x)*[Li(x)-\pi(x)] + c - \sum_{y=2}^{x-1}\frac{1}{y}[Li(y)-\pi(y)]?
With an approximation of Li(y)-\pi(y) ~ (1+\frac{2}{\log{y}})*\frac{\sqrt{y}}{\log{y}}, an error bound of O(x^{\frac{1}{3}}) seems possible.
I still have no clue how useful this might be. Not only for this particular thread, but in general, and for large n.
I mean, for small x, say <=10^16, some already known values for \theta(x) can be used to determine the summation term, and then interpolate accordingly. Precomputed terms could go with a spacing of 10^10 or even 10^11 to get an error of ± 1 in x-\theta(x). But I'll wait for answers before I work out any details there.
mart_r is offline   Reply With Quote
Reply

Thread Tools


Similar Threads
Thread Thread Starter Forum Replies Last Post
The Fibonacci-Chebyshev connection Dr Sardonicus Number Theory Discussion Group 5 2017-07-09 13:18
Do normal adults give themselves an allowance? (...to fast or not to fast - there is no question!) jasong jasong 35 2016-12-11 00:57
Fast Approximate Ceiling Function R.D. Silverman Programming 27 2010-11-19 17:50
Chebyshev's Estimates brownkenny Math 2 2009-01-22 17:21
Birkhoff and Hall's theta function Dougy Math 2 2009-01-05 05:09

All times are UTC. The time now is 20:02.


Sun Dec 5 20:02:03 UTC 2021 up 135 days, 14:31, 1 user, load averages: 4.25, 2.94, 2.12

Powered by vBulletin® Version 3.8.11
Copyright ©2000 - 2021, Jelsoft Enterprises Ltd.

This forum has received and complied with 0 (zero) government requests for information.

Permission is granted to copy, distribute and/or modify this document under the terms of the GNU Free Documentation License, Version 1.2 or any later version published by the Free Software Foundation.
A copy of the license is included in the FAQ.