Ciao, Shane,
On Sat, Nov 6, 2021 at 1:31 PM Marco Bodrato <bodr...@mail.dm.unipi.it>
wrote:
Il 2021-11-03 19:18 Shane Neph ha scritto:
> mpz_primorial_ui( ) is a great function. Why not return the actual
> number of primes that go into that calculation?
It may make sense, but only if we add also another function: a way to
compute the number of primes in a range; or at least in the
range [1..n], i.e. pi(n) for a given unsigned long n.
Il 2021-11-06 22:29 Shane Neph ha scritto:
That would be a nice solution.
If you really need both the primorial and the prime count, my solution
is somehow a waste, because no memory of the sieved numbers is kept
between the library calls. The numbers are sieved twice...
I got inspired by a piece of code that Seth Troisi sent a couple of
years ago and I wrote a gmp_pi_ui(n) function with two variants, a
trivial one (it calls gmp_primesieve) and an implementation of the
Legendre strategy.
I attach it for comments. Compiled with -DMAIN it prints the number of
primes up to increasing powers of two. On my computer, with small ranges
it is slower than GP/Pari, for large ranges, it is faster...
Ĝis,
m
/* pi (N) -- Returns the number of primes in the range [1..N].
Copyright 2019-2021 Marco Bodrato.
*/
#include "gmp-impl.h"
#include "longlong.h"
#ifndef PI_LEGENDRE_THRESHOLD
#define PI_LEGENDRE_THRESHOLD 80000
#endif
/* n_fto_bit (n) = ((n-1)&(-CNST_LIMB(2)))/3U-1 */
static mp_limb_t
n_fto_bit (mp_limb_t n) { return ((n-5)|1)/3U; }
/* n_fto_bit(x) + 2 = gmp_legendre_pi_phi_2 (x) */
static mp_size_t
primesieve_size (mp_limb_t n) { return n_fto_bit(n) / GMP_LIMB_BITS + 1; }
static mp_limb_t
limb_apprsqrt (mp_limb_t x)
{
int s;
ASSERT (x > 2);
count_leading_zeros (s, x);
s = (GMP_LIMB_BITS - s) >> 1;
return ((CNST_LIMB(1) << (s - 1)) + (x >> 1 >> s));
}
/* phi(a, 0) = a
* phi(a, b) = phi(a, b-1) - phi(a/p_b, b-1)
* phi(a, 1) = phi(a, 0) - phi(a/2, 0) = a - a/2
* phi(a, 2) = phi(a, 1) - phi(a/3, 1) = a - a/2 - (a/3 - a/6) ~= a/3
* phi(a, b) = phi(a, 2) - sum_{1<c<b} phi(a/p_c, c)
*
* If a < p_b^2, phi(a, b) = pi(a) - b + 1. Not jet implemented.
* Only the case a < p_{b+1}, then phi(a, b) = 1, is implemented;
* slightly extended to 0 <= a-p_{b+1} is very small, phi (a,b) = 1+1.
*
* This implementation receives primes = b - 2, the number of primes
* in sieve (2 and 3 are precomputed and not counted in the sieve)
*/
static long
gmp_legendre_pi_phi_2 (unsigned long x /*, 0 (= b-2) */ )
{
return ((x + 1) | 1) / 3UL;
}
static long
gmp_legendre_pi_phi(unsigned long n, unsigned primes, mp_ptr sieve)
{
/* n -= 1 - (n & 1); /\* Reduce to odd *\/ */
long res = gmp_legendre_pi_phi_2 (n);
/* res precomputed for primes 2 and 3 */
ASSERT (primes > 0);
/* if (primes == 0) */
/* return res; */
/* if (res * 3 == n) n -= 2; */
unsigned count = 0;
unsigned long i = 4, *sp = sieve;
do { /* Loop on primes */
for (unsigned long b = i, x = ~ *(sp++); x != 0; b += 3, x >>= 1)
if (x & 1)
{
/* b = 4 + 3*k :
* if k = 2n, b = 4 (mod 6) is even, the prime is b+1,
* and b+2 = 0 (mod 6) is composite;
* if k = 2n+1, b = 1 (mod 6) is odd, b is the prime,
* and b+1 (even), b+2 = 3 (mod 6), b+2+(b%2) (even) are composite.
*/
unsigned long frac = n / (b | 1);
if (frac <= (b + 2 + (b & 1)))
return res - primes + count - (frac >= (b | 1));
if (count < 2) /* Directly compute, avoid recursion. */
res -= gmp_legendre_pi_phi_2 (frac) /* count == 0 */
- (count ? gmp_legendre_pi_phi_2 (frac / 5) : 0); /* count == 1 */
else
res -= gmp_legendre_pi_phi(frac, count, sieve);
if (primes == ++count)
return res;
}
i += GMP_LIMB_BITS * 3;
} while (1);
}
/* Takes a pre-sieved range of primes [5..maxprime].
* The range contains primes_in_sieve primes.
*
* Legendre showed that, if maxprime^2 >= n, then
* primepi(n) = primepi(maxprime) + phi(n, primepi(maxprime)) - 1
* where phi (a, b) = numbers x <= a, not divisible by the first b primes.
*
* The sieve does not count 2 and 3, so that
* primepi(maxprime) = primes_in_sieve + 2
*
* Here we pass to pi_phi the prime index of b: primes_in_sieve.
*/
static long
gmp_legendre_pi_ps (unsigned long n, unsigned primes_in_sieve, mp_ptr sieve)
{
return primes_in_sieve + 1 + gmp_legendre_pi_phi(n, primes_in_sieve, sieve);
}
static long
gmp_legendre_pi_ui(unsigned long n)
{
unsigned long sqrt_n = limb_apprsqrt (n);
unsigned sieve_size = primesieve_size (sqrt_n);
TMP_SDECL;
TMP_SMARK;
mp_ptr sieve = TMP_SALLOC_LIMBS (sieve_size);
unsigned count = gmp_primesieve (sieve, sqrt_n);
long res = gmp_legendre_pi_ps (n, count, sieve);
TMP_SFREE;
return res;
}
/* TODO: the _primesieve interface should be extended, to allow
piece-wise sieving. Storing the full sieve is wasteful here.
*/
static long
gmp_persieve_pi_ui (unsigned long n)
{
mp_size_t size = primesieve_size (n);
TMP_DECL;
TMP_MARK;
long res = gmp_primesieve (TMP_ALLOC_LIMBS (size), n);
TMP_FREE;
return res + 2;
}
long
gmp_pi_ui (unsigned long n)
{
if (n < 7)
return (n>1) + (n>2) + (n>4);
if (BELOW_THRESHOLD (n, PI_LEGENDRE_THRESHOLD))
return gmp_persieve_pi_ui (n);
else
return gmp_legendre_pi_ui (n);
}
#ifdef MAIN
int main() {
unsigned step = 2;
for (unsigned long n=1; n< ~0UL / step ; n *= step)
gmp_printf ("%lu -> %ld\n", n, gmp_pi_ui (n));
}
#endif
_______________________________________________
gmp-devel mailing list
gmp-devel@gmplib.org
https://gmplib.org/mailman/listinfo/gmp-devel
