mersenneforum.org Requestion for fast chebyshev theta function
 Register FAQ Search Today's Posts Mark Forums Read

 2018-03-04, 17:08 #1 danaj   "Dana Jacobsen" Feb 2011 Bangkok, TH 32×101 Posts 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.
2018-03-04, 21:54   #2
danaj

"Dana Jacobsen"
Feb 2011
Bangkok, TH

38D16 Posts

Quote:
 Originally Posted by JM Montolio A ¿ 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.

 2018-03-05, 00:27 #3 CRGreathouse     Aug 2006 3×1,993 Posts I don't know of a clever algorithm. Batched logs (ideally with some guard bits) are the best I could come up with, too.
 2018-03-05, 17:07 #4 Dr Sardonicus     Feb 2017 Nowhere 22×1,459 Posts 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.
 2018-03-05, 18:30 #5 mart_r     Dec 2008 you know...around... 7·109 Posts 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
 2018-03-06, 14:22 #6 danaj   "Dana Jacobsen" Feb 2011 Bangkok, TH 32×101 Posts 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.
2018-03-06, 17:04   #7
axn

Jun 2003

19×283 Posts

Quote:
 Originally Posted by danaj 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

2018-03-06, 20:47   #8
R. Gerbicz

"Robert Gerbicz"
Oct 2005
Hungary

112·13 Posts

Quote:
 Originally Posted by danaj 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]

2018-03-08, 13:33   #9
Dr Sardonicus

Feb 2017
Nowhere

583610 Posts

Quote:
 Originally Posted by R. Gerbicz 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.

 2018-03-31, 14:59 #10 mart_r     Dec 2008 you know...around... 7·109 Posts 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.

 Similar Threads Thread Thread Starter Forum Replies Last Post Dr Sardonicus Number Theory Discussion Group 6 2022-01-15 12:13 jasong jasong 35 2016-12-11 00:57 R.D. Silverman Programming 27 2010-11-19 17:50 brownkenny Math 2 2009-01-22 17:21 Dougy Math 2 2009-01-05 05:09

All times are UTC. The time now is 16:23.

Tue Jul 5 16:23:17 UTC 2022 up 82 days, 14:24, 1 user, load averages: 2.25, 1.59, 1.34