Ciao Shane, Il 2021-11-11 23:58 Shane Neph ha scritto:
First, thank you for your work. I will be using it. I'd like to see it or similar become part of the GMP library.
I'm not sure you can use it! :-D I forgot a licence, so the default is "all rights reserved".I send the code again, with a GPLv3+ licence, and a correction such that the function works also for the value 2^32-1 in a 32-bits environment.
Moreover, remember that this implementation uses functions that are internal, undocumented, and because of this subject to unexpected changes in the future.
Before seeing something like this in GMP, we have to at least decide the interface :-) To write the function, I got inspiration from a code by Seth Troisi, who used more than a call to the prime counting function to get the n-th prime number. Should we allow recycling parts of the computations between calls? How?
I also would like to see the primordial function return that same value because it knows the value.
True, but this is not the only information that the function "knows" and might be useful. Why not returning, e.g., the larger prime?
In my work, I use primordial in ways similar to factorials, and I use it over a factorial because it simply grows faster.
Not the implementation in GMP :-) mpz_primorial_ui is the "square free" version of mpz_fac_ui :-D Strictly smaller, when the argument is larger than 3. Ĝis, m -- http://bodrato.it/papers/
/* pi (N) -- Returns the number of primes in the range [1..N]. Contributed to the GNU project by Marco Bodrato. THE FUNCTION IN THIS FILE IS INTERNAL WITH A MUTABLE INTERFACE. IT IS ONLY SAFE TO REACH IT THROUGH DOCUMENTED INTERFACES. IN FACT, IT IS ALMOST GUARANTEED THAT IT WILL CHANGE OR DISAPPEAR IN A FUTURE GNU MP RELEASE. Copyright 2019-2021 Free Software Foundation, Inc. This file was written for the GNU MP Library, and it is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 3 of the License, or (at your option) any later version. The GNU MP Library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received copies of the GNU General Public License and the GNU Lesser General Public License along with the GNU MP Library. If not, see https://www.gnu.org/licenses/. */ #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) = gmp_legendre_pi_phi_2 (x) - 2 */ 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) */ ) { /* The following condition should never be true. */ #if (ULONG_MAX % 3 != 0) && (ULONG_MAX % 7 != 0) && (ULONG_MAX % 31 != 0) && (ULONG_MAX % 127 != 0) && (ULONG_MAX % 23 != 0) && (ULONG_MAX % 8191 != 0) && (ULONG_MAX % 47 != 0) && (ULONG_MAX % 233 != 0) if (UNLIKEY (x == ULONG_MAX)) return (x - 2 + (x & 1)) / 3UL + 1; else #endif 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) { #if (ULONG_MAX % 3 == 0) || (ULONG_MAX % 7 == 0) || (ULONG_MAX % 31 == 0) || (ULONG_MAX % 127 == 0) || (ULONG_MAX % 23 == 0) || (ULONG_MAX % 8191 == 0) || (ULONG_MAX % 47 == 0) || (ULONG_MAX % 233 == 0) n -= (n == ULONG_MAX); /* Avoid overflow in gmp_legendre_pi_phi_2 */ #endif 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 long n = 1, odd = 1, step = 2; do { n *= step; n |= odd; gmp_printf ("%lu -> %ld\n", n, gmp_pi_ui (n)); } while (n <= ~0UL / step); } #endif
_______________________________________________ gmp-devel mailing list gmp-devel@gmplib.org https://gmplib.org/mailman/listinfo/gmp-devel
