prime_counts.c

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "ptypes.h"
#include "sieve.h"
#include "cache.h"
#include "lmo.h"
#include "constants.h"
#include "prime_counts.h"
#include "util.h"
#include "util_bits.h"
#include "util_math.h"
#include "real.h"
#include "mathl.h"

#if defined(__GNUC__)
 #define word_unaligned(m,wordsize)  ((uintptr_t)m & (wordsize-1))
#else  /* uintptr_t is part of C99 */
 #define word_unaligned(m,wordsize)  ((unsigned int)m & (wordsize-1))
#endif

static const unsigned char byte_zeros[256] =
  {8,7,7,6,7,6,6,5,7,6,6,5,6,5,5,4,7,6,6,5,6,5,5,4,6,5,5,4,5,4,4,3,
   7,6,6,5,6,5,5,4,6,5,5,4,5,4,4,3,6,5,5,4,5,4,4,3,5,4,4,3,4,3,3,2,
   7,6,6,5,6,5,5,4,6,5,5,4,5,4,4,3,6,5,5,4,5,4,4,3,5,4,4,3,4,3,3,2,
   6,5,5,4,5,4,4,3,5,4,4,3,4,3,3,2,5,4,4,3,4,3,3,2,4,3,3,2,3,2,2,1,
   7,6,6,5,6,5,5,4,6,5,5,4,5,4,4,3,6,5,5,4,5,4,4,3,5,4,4,3,4,3,3,2,
   6,5,5,4,5,4,4,3,5,4,4,3,4,3,3,2,5,4,4,3,4,3,3,2,4,3,3,2,3,2,2,1,
   6,5,5,4,5,4,4,3,5,4,4,3,4,3,3,2,5,4,4,3,4,3,3,2,4,3,3,2,3,2,2,1,
   5,4,4,3,4,3,3,2,4,3,3,2,3,2,2,1,4,3,3,2,3,2,2,1,3,2,2,1,2,1,1,0};
static UV count_zero_bits(const unsigned char* m, UV nbytes)
{
  UV count = 0;
#if BITS_PER_WORD == 64
  if (nbytes >= 16) {
    while ( word_unaligned(m,sizeof(UV)) && nbytes--)
      count += byte_zeros[*m++];
    if (nbytes >= 8) {
      UV* wordptr = (UV*)m;
      UV nwords = nbytes / 8;
      UV nzeros = nwords * 64;
      m += nwords * 8;
      nbytes %= 8;
      while (nwords--)
        nzeros -= popcnt(*wordptr++);
      count += nzeros;
    }
  }
#endif
  while (nbytes--)
    count += byte_zeros[*m++];
  return count;
}

/* Given a sieve of size nbytes, walk it counting zeros (primes) until:
 *
 * (1) we counted them all: return the count, which will be less than maxcount.
 *
 * (2) we hit maxcount: set position to the index of the maxcount'th prime
 *     and return count (which will be equal to maxcount).
 */
static UV count_segment_maxcount(const unsigned char* sieve, UV base, UV nbytes, UV maxcount, UV* pos)
{
  UV count = 0;
  UV byte = 0;
  const unsigned char* sieveptr = sieve;
  const unsigned char* maxsieve = sieve + nbytes;

  MPUassert(sieve != 0, "count_segment_maxcount incorrect args");
  MPUassert(pos != 0, "count_segment_maxcount incorrect args");
  *pos = 0;
  if ( (nbytes == 0) || (maxcount == 0) )
    return 0;

  /* Do fixed-length word counts to start, with possible overcounting */
  while ((count+64) < maxcount && sieveptr < maxsieve) {
    UV top = base + 3*maxcount;
    UV div = (top <       8000) ? 8 :     /* 8 cannot overcount */
             (top <    1000000) ? 4 :
             (top <   10000000) ? 3 : 2;
    UV minbytes = (maxcount-count)/div;
    if (minbytes > (UV)(maxsieve-sieveptr)) minbytes = maxsieve-sieveptr;
    count += count_zero_bits(sieveptr, minbytes);
    sieveptr += minbytes;
  }
  /* Count until we reach the end or >= maxcount */
  while ( (sieveptr < maxsieve) && (count < maxcount) )
    count += byte_zeros[*sieveptr++];
  /* If we went too far, back up. */
  while (count >= maxcount)
    count -= byte_zeros[*--sieveptr];
  /* We counted this many bytes */
  byte = sieveptr - sieve;

  MPUassert(count < maxcount, "count_segment_maxcount wrong count");

  if (byte == nbytes)
    return count;

  /* The result is somewhere in the next byte */
  START_DO_FOR_EACH_SIEVE_PRIME(sieve, 0, byte*30+1, nbytes*30-1)
    if (++count == maxcount)  { *pos = p; return count; }
  END_DO_FOR_EACH_SIEVE_PRIME;

  MPUassert(0, "count_segment_maxcount failure");
  return 0;
}

static UV count_in_one_sieve_byte(const unsigned char* sieve, UV lowp, UV highp)
{
  UV count = 0, d = lowp/30;
  unsigned char bit, s;

  /* Count unmarked candidates in [lowp,highp].
   * Both bounds must lie in the same wheel-30 sieve byte. */
  MPUassert(d == highp/30, "count_in_one_sieve_byte must not span bytes");
  s = sieve[d] | clearprev30[lowp-d*30];
  while (s != 0xFF) {
    UV p;
    bit = nextzero30[s];
    s |= 1 << bit;
    p = d*30 + wheel30[bit];
    if (p > highp) break;
    count++;
  }
  return count;
}

/* Given a sieve of size nbytes, counting zeros (primes) but excluding the
 * areas outside lowp and highp.
 */
static UV count_segment_ranged(const unsigned char* sieve, UV nbytes, UV lowp, UV highp)
{
  UV count, hi_d, lo_d, lo_m;

  MPUassert( sieve != 0, "count_segment_ranged incorrect args");
  if (nbytes == 0) return 0;

  count = 0;
  hi_d = highp/30;

  if (hi_d >= nbytes) {
    hi_d = nbytes-1;
    highp = hi_d*30+29;
  }

  if (highp < lowp)
    return 0;

#if 0
  /* Dead simple way */
  START_DO_FOR_EACH_SIEVE_PRIME(sieve, 0, lowp, highp)
    count++;
  END_DO_FOR_EACH_SIEVE_PRIME;
  return count;
#endif

  lo_d = lowp/30;
  lo_m = lowp - lo_d*30;
  /* Count first fragment */
  if (lo_m > 1) {
    UV upper = (highp <= (lo_d*30+29)) ? highp : (lo_d*30+29);
    count += count_in_one_sieve_byte(sieve, lowp, upper);
    lowp = upper+2;
    lo_d = lowp/30;
  }
  if (highp < lowp)
    return count;

  /* Count bytes in the middle */
  {
    UV hi_m = highp - hi_d*30;
    UV count_bytes = hi_d - lo_d + (hi_m == 29);
    if (count_bytes > 0) {
      count += count_zero_bits(sieve+lo_d, count_bytes);
      lowp += 30*count_bytes;
    }
  }
  if (highp < lowp)
    return count;

  /* Count the final partial byte without reading a complete UV past it. */
  count += count_in_one_sieve_byte(sieve, lowp, highp);

  return count;
}


/*
 * The pi(x) prime count functions.  prime_count(x) gives an exact number,
 * but requires determining all the primes up to x, so will be much slower.
 *
 * prime_count_lower(x) and prime_count_upper(x) give lower and upper limits,
 * which will bound the exact value.  These bounds should be fairly tight.
 *
 *    pi_upper(x) - pi(x)                      pi_lower(x) - pi(x)
 *    <    10   for x <         5_371          <    10   for x <        9_437
 *    <    50   for x <       295_816          <    50   for x <      136_993
 *    <   100   for x <     1_761_655          <   100   for x <      909_911
 *    <   200   for x <     9_987_821          <   200   for x <    8_787_901
 *    <   400   for x <    34_762_891          <   400   for x <   30_332_723
 *    <  1000   for x <   372_748_528          <  1000   for x <  233_000_533
 *    <  5000   for x < 1_882_595_905          <  5000   for x <  over 4300M
 *
 * The average of the upper and lower bounds is within 9 for all x < 15809, and
 * within 50 for all x < 1_763_367.
 *
 * It is common to use the following Chebyshev inequality for x >= 17:
 *    1*x/logx <-> 1.25506*x/logx
 * but this gives terribly loose bounds.
 *
 * Rosser and Schoenfeld's bound for x >= 67 of
 *    x/(logx-1/2) <-> x/(logx-3/2)
 * is much tighter.  These bounds can be tightened even more.
 *
 * The formulas of Dusart for higher x are better yet.  I recommend the paper
 * by Burde for further information.  Dusart's thesis is also a good resource.
 *
 * I have tweaked the bounds formulas for small (under 70_000M) numbers so they
 * are tighter.  These bounds are verified via trial.  The Dusart bounds
 * (1.8 and 2.51) are used for larger numbers since those are proven.
 *
 */

#include "prime_count_tables.h"

UV segment_prime_count(UV low, UV high)
{
  const unsigned char* cache_sieve;
  unsigned char* segment;
  UV segment_size, low_d, high_d;
  UV count = 0;

  if ((low <= 2) && (high >= 2)) count++;
  if ((low <= 3) && (high >= 3)) count++;
  if ((low <= 5) && (high >= 5)) count++;
  if (low < 7)  low = 7;

  if (low > high)  return count;

#if !defined(BENCH_SEGCOUNT)
  if (low == 7 && high <= 30*NPRIME_SIEVE30) {
    count += count_segment_ranged(prime_sieve30, NPRIME_SIEVE30, low, high);
    return count;
  }

  /* If we have sparse prime count tables, use them here.  These will adjust
   * 'low' and 'count' appropriately for a value slightly less than ours.
   * This should leave just a small amount of sieving left.  They stop at
   * some point, e.g. 3000M, so we'll get the answer to that point then have
   * to sieve all the rest.  We should be using LMO or Lehmer much earlier. */
#ifdef APPLY_TABLES
  APPLY_TABLES
#endif
#endif

  low_d = low/30;
  high_d = high/30;

  /* Count full bytes only -- no fragments from primary cache */
  segment_size = get_prime_cache(0, &cache_sieve) / 30;
  if (segment_size < high_d) {
    /* Expand sieve to sqrt(n) */
    UV endp = (high_d >= (UV_MAX/30))  ?  UV_MAX-2  :  30*high_d+29;
    UV newsize = (UV)isqrt(endp)+1;
    if (newsize > 2642245) newsize = 2642245;  /* Limit to icbrt(2^64) */
    release_prime_cache(cache_sieve);
    segment_size = get_prime_cache( newsize, &cache_sieve) / 30;
  }

  if ( (segment_size > 0) && (low_d <= segment_size) ) {
    /* Count all the primes in the primary cache in our range */
    count += count_segment_ranged(cache_sieve, segment_size, low, high);

    if (high_d < segment_size) {
      release_prime_cache(cache_sieve);
      return count;
    }

    low_d = segment_size;
    if (30*low_d > low)  low = 30*low_d;
  }
  release_prime_cache(cache_sieve);

  /* More primes needed.  Repeatedly segment sieve. */
  {
    void* ctx = start_segment_primes(low, high, &segment);
    UV seg_base, seg_low, seg_high;
    while (next_segment_primes(ctx, &seg_base, &seg_low, &seg_high)) {
      segment_size = seg_high/30 - seg_low/30 + 1;
      count += count_segment_ranged(segment, segment_size, seg_low-seg_base, seg_high-seg_base);
    }
    end_segment_primes(ctx);
  }

  return count;
}

UV prime_count_range(UV lo, UV hi)
{
  if (lo > hi || hi < 2)
    return 0;

#if defined(BENCH_SEGCOUNT)
  return segment_prime_count(lo, hi);
#endif

  /* We use table acceleration so this is preferable for small inputs */
  if (hi < _MPU_LMO_CROSSOVER)  return segment_prime_count(lo, hi);

  { /* Rough empirical threshold for when segment faster than LMO */
    UV range_threshold = hi / (isqrt(hi)/200);
    if ((hi-lo) < range_threshold)  /* avoid the +1 so we don't overflow */
      return segment_prime_count(lo, hi);
  }
  return LMO_prime_count(hi) - ((lo < 2) ? 0 : LMO_prime_count(lo-1));
}

UV prime_count(UV n)
{
  if (n < 2) return 0;

  /* We use table acceleration so this is preferable for small inputs */
  if (n < _MPU_LMO_CROSSOVER)  return segment_prime_count(0, n);

  return LMO_prime_count(n);
}

UV prime_count_approx(UV n)
{
  if (n < 3000000) return segment_prime_count(2, n);
  return (UV) (RiemannR((long double) n, 1e-6) + 0.5);
}


UV prime_count_lower(UV n)
{
  long double fn, fl1, fl2, lower, a;

  if (n < 33000) return segment_prime_count(2, n);

  fn  = (long double) n;
  fl1 = logl(n);
  fl2 = fl1 * fl1;

  /* Axler 2014: https://arxiv.org/abs/1409.1780  (v7 2016), Cor 3.6
   * show variations of this. */

  if (n <= 300070) { /* Quite accurate and avoids calling Li for speed. */
    /* Based on Axler 2022, page 9, Corollary 5.1 */
    a = (n <  69720) ? 905 :
        (n <  70120) ? 961 :
        (n <  88800) ? 918.2 :
        (n < 176000) ? 887.7 :
        (n < 299270) ? 839.46 :
                       846.66;  /* Good to 300071 */
    lower = fn / (fl1 - 1 - 1/fl1 - 2.975666/fl2 - 13.024334/(fl1*fl2) + a/(fl2*fl2));
  } else if (n < UVCONST(4000000000)) {
    /* Loose enough that FP differences in Li(n) should be ok. */
    a = (n <     88783) ?   4.0L
      : (n <    300000) ?  -3.0L
      : (n <    303000) ?   5.0L
      : (n <   1100000) ?  -7.0L
      : (n <   4500000) ? -37.0L
      : (n <  10200000) ? -70.0L
      : (n <  36900000) ? -53.0L
      : (n <  38100000) ? -29.0L
      :                   -84.0L;
    lower = Li(fn) - (sqrtl(fn)/fl1) * (1.94L + 2.50L/fl1 + a/fl2);
  } else {                         /* Dusart 2018 Lemma 2.2 */
    lower = Li(fn) - 2*sqrtl(fn)/fl1;
  }
  return (UV) ceill(lower);
}

typedef struct {
  UV thresh;
  float aval;
} thresh_t;

static const thresh_t _upper_thresh[] = {
  {     59000, 2.48f },
  {    355991, 2.54f },
  {   3550000, 2.51f },
  {   3560000, 2.49f },
  {   5000000, 2.48f },
  {   8000000, 2.47f },
  {  13000000, 2.46f },
  {  18000000, 2.45f },
  {  31000000, 2.44f },
  {  41000000, 2.43f },
  {  48000000, 2.42f },
  { 119000000, 2.41f },
  { 182000000, 2.40f },
  { 192000000, 2.395f },
  { 213000000, 2.390f },
  { 271000000, 2.385f },
  { 322000000, 2.380f },
  { 400000000, 2.375f },
  { 510000000, 2.370f },
  { 682000000, 2.367f },
  { UVCONST(2953652287), 2.362f }
};
#define NUPPER_THRESH (sizeof(_upper_thresh)/sizeof(_upper_thresh[0]))

UV prime_count_upper(UV n)
{
  int i;
  long double fn, fl1, fl2, upper, a;

  if (n < 33000) return segment_prime_count(2, n);

  fn  = (long double) n;
  fl1 = logl(n);
  fl2 = fl1 * fl1;

  /* Axler 2014: https://arxiv.org/abs/1409.1780  (v7 2016), Cor 3.5
   *
   * upper = fn/(fl1-1.0L-1.0L/fl1-3.35L/fl2-12.65L/fl3-89.6L/fl4);
   * return (UV) floorl(upper);
   *
   * Axler 2022: https://arxiv.org/pdf/2203.05917.pdf (v4 2022): Cor 4.4
   *  Cor 4.4: fl2 3.024334L fl3 12.975666 fl4 71.048668 fl5 533.594
   *  Thm 1.3: fl2,fl3,fl4 as above   fl5 461.364417856444 fl6 4331.1
   *
   * These do not improve on the bound of Li(x) which covers all 64-bit n.
   */

  if (BITS_PER_WORD == 32 || fn <= 821800000.0) {  /* Dusart 2010, page 2 */
    for (i = 0; i < (int)NUPPER_THRESH; i++)
      if (n < _upper_thresh[i].thresh)
        break;
    a = (i < (int)NUPPER_THRESH)  ?  _upper_thresh[i].aval  :  2.334L;
    upper = fn/fl1 * (1.0L + 1.0L/fl1 + a/fl2);
  } else if (fn < 101260000000.0) {
    /* For some smaller inputs we further tighten */
    a = (fn <   1100000000.0) ? 0.032    /* Empirical */
      : (fn <  10010000000.0) ? 0.027    /* Empirical */
      : (fn < 101260000000.0) ? 0.021    /* Empirical */
                              : 0.0;
    upper = Li(fn) - a * fl1*sqrtl(fn)/25.132741228718345907701147L;
  } else {
    /* Büthe  2015 1.10      : Pi(x) <= Li(x) for x <= 10^19 */
    /* Dusart 2018 Lemma 2.2 : Pi(x) <= Li(x) for x <= 10^20 */
    upper = Li(fn);
  }
  return (UV) floorl(upper);
}

static void simple_nth_limits(UV *lo, UV *hi, long double n, long double logn, long double loglogn) {
  const long double a = (n < 228) ? .6483 : (n < 948) ? .8032 : (n < 2195) ? .8800 : (n < 39017) ? .9019 : .9484;
  *lo = n * (logn + loglogn - 1.0 + ((loglogn-2.10)/logn));
  *hi = n * (logn + loglogn - a);
  if (*hi < *lo) *hi = MPU_MAX_PRIME;
}

/* The nth prime will be less or equal to this number */
UV nth_prime_upper(UV n)
{
  long double fn, flogn, flog2n, upper, c, d;

  if (n < NPRIMES_SMALL)
    return primes_small[n];
  if (n >= MPU_MAX_PRIME_IDX)
    return n == MPU_MAX_PRIME_IDX  ?  MPU_MAX_PRIME  :  0;

  fn     = (long double) n;
  flogn  = logl(n);
  flog2n = logl(flogn);    /* Note distinction between log_2(n) and log^2(n) */

  /* Binary search on prime count lower.  Good but quite slow. */
  if (n < 15360) {
    UV lo,hi;
    simple_nth_limits(&lo, &hi, fn, flogn, flog2n);
    while (lo < hi) {
      UV mid = lo + (hi-lo)/2;
      if (prime_count_lower(mid) < n) lo = mid+1;
      else                            hi = mid;
    }
    return lo;
  }

  /* Under RH, we could use Dusart 2018, Theorem 3.1.
   * It is better than the below starting at 230779190522.
   *
   * long double klogk = fn*flogn, logklogk = logl(klogk);
   * long double E = sqrt(klogk)*logklogk*logklogk/25.132741228718345907701147L;
   * return inverse_li(n) + floorl(E);
   *
   * Empirically, x > 46254833, return inverse_li(n) + floorl(E * 2.0/flogn);
   * e.g. 1e16 we get 1375558510 vs 1507803850527  3 orders better
   * about 10x slower though.
   * verified for x <= 1000000000 (1e9)
   */
#if 0
  if (n > 46254833) {
    long double klogk = fn*flogn, logklogk = logl(klogk);
    long double E = sqrt(klogk)*logklogk*logklogk/25.132741228718345907701147L;
    return inverse_li(n) + floorl( E * 2.0/flogn );
  }
#endif

  /* See: Axler 2013, Dusart 2010 */
  /*      Axler 2017: http://arxiv.org/pdf/1706.03651.pdf */

  if      (n >= 46254381) { c = 2.00; d =  10.667;  } /* Axler 2017 Cor 1.2 */
  else if (n >=  8009824) { c = 2.00; d =  10.273;  } /* Axler 2013 Kor G */
  /* This is about 3x better than Dusart (2010) for 688382-8009823:
   *
   * else if (n >=   688382) { c = 2.30; d =   0.5730; }
   *
   * but we can split the range and get another 2x improvement in MSE.
   */
  else if (n >=  5450000) { c = 2.00; d =  10.1335; }  /*5450-8009 */
  else if (n >=  3906280) { c = 1.67; d =  20.2675; }  /*3906-5450 */
  else if (n >=  2110840) { c = 2.51; d =  -5.5714; }  /*2110-3906 */
  else if (n >=   876700) { c = 2.49; d =  -4.5129; }  /* 877-2110 */
  else if (n >=   688382) { c = 3.31; d = -26.3858; }  /* 688-877 */
  /* Use the Axler framework to get good bounds for smaller inputs. */
  else if (n >=   575750) { c =-0.79; d =  83.5215; }  /* 580-688 */
  else if (n >=   467650) { c = 0.93; d =  37.1597; }  /* 467-580 */
  else if (n >=   382440) { c = 2.92; d = -15.4768; }  /* 382-467 */
  else if (n >=   301130) { c = 5.92; d = -91.3415; }  /* 301-382 */
  else if (n >=   138630) { c = 2.01; d =   7.2842; }  /* 138-301 */
  else if (n >=    85820) { c = 2.07; d =   5.2103; }  /*  86-138 */
  else if (n >=    39016) { c = 2.76; d = -11.5918; }  /*  39- 86 */
  else if (n >=    31490) { c = 1.49; d =  15.1821; }  /*  31- 39 */
  else if (n >=    25070) { c =11.89; d =-197.8951; }  /*  25- 31 */
  else if (n >=    15359) { c = 4.80; d = -51.5928; }  /*  15- 25 */
  else                    { c = 3.92; d = -33.3994; }  /*   0- 15 */

  upper = fn * ( flogn + flog2n - 1.0 + ((flog2n-c)/flogn)
                 - (flog2n*flog2n-6*flog2n+d)/(2*flogn*flogn) );

  if (upper >= (long double)UV_MAX) {
    if (n <= MPU_MAX_PRIME_IDX) return MPU_MAX_PRIME;
    croak("nth_prime_upper(%"UVuf") overflow", n);
  }

  return (UV) floorl(upper);
}

/* The nth prime will be greater than or equal to this number */
UV nth_prime_lower(UV n)
{
  double fn, flogn, flog2n;
  UV plower;

  if (n < NPRIMES_SMALL)
    return primes_small[n];
  if (n >= MPU_MAX_PRIME_IDX)
    return n == MPU_MAX_PRIME_IDX  ?  MPU_MAX_PRIME  :  0;

  fn     = (double) n;
  flogn  = log(n);
  flog2n = log(flogn);

  /* For small values, do a binary search on the inverse prime count */
  if (n < 2000000) {
    UV lo,hi;
    simple_nth_limits(&lo, &hi, fn, flogn, flog2n);
    while (lo < hi) {
      UV mid = lo + (hi-lo)/2;
      if (prime_count_upper(mid) < n) lo = mid+1;
      else                            hi = mid;
    }
    return lo;
  }

  { /* Axler 2017 http://arxiv.org/pdf/1706.03651.pdf Corollary 1.4 */
    double b1 = (n < 56000000)  ?  11.200  :  11.50800000002;
    double lower = fn * (flogn + flog2n-1.0 + ((flog2n-2.00)/flogn) - ((flog2n*flog2n-6*flog2n+b1)/(2*flogn*flogn)));
    plower = (UV) ceill(lower);
  }
  return plower < MPU_MAX_PRIME ? plower : MPU_MAX_PRIME;
}

UV nth_prime_approx(UV n)
{
  return (n < NPRIMES_SMALL)  ?  primes_small[n]  :  inverse_R(n);
}


UV nth_prime(UV n)
{
  const unsigned char* cache_sieve;
  unsigned char* segment;
  UV upper_limit, segbase, segment_size, p, count, target;

  /* If very small, return the table entry */
  if (n < NPRIMES_SMALL)
    return primes_small[n];
  if (n >= MPU_MAX_PRIME_IDX)
    return n == MPU_MAX_PRIME_IDX  ?  MPU_MAX_PRIME  :  0;

  /* Determine a bound on the nth prime.  We know it comes before this. */
  upper_limit = nth_prime_upper(n);
  MPUassert(upper_limit > 0, "nth_prime got an upper limit of 0");
  p = count = 0;
  target = n-3;

  /* For relatively small values, generate a sieve and count the results.
   *
   * For larger values, compute an approximate low estimate, use our fast
   * prime count, then segment sieve forwards or backwards for the rest.
   */
  if (upper_limit <= get_prime_cache(0, 0) || upper_limit <= 32*1024*30) {
    /* Generate a sieve and count. */
    segment_size = get_prime_cache(upper_limit, &cache_sieve) / 30;
    /* Count up everything in the cached sieve. */
    if (segment_size > 0)
      count += count_segment_maxcount(cache_sieve, 0, segment_size, target, &p);
    release_prime_cache(cache_sieve);
  } else {
    /* A binary search on RiemannR is nice, but ends up either often being
     * being higher (requiring going backwards) or biased and then far too
     * low.  Using the inverse Li is easier and more consistent. */
    UV lower_limit = inverse_li(n);
    UV correction = inverse_li(isqrt(n))/4;
    /* For even better performance, add in half the usual correction, which
     * will get us even closer, so even less sieving required.  However, it
     * is now possible to get a result higher than the value, so we'll need
     * to handle that case.  It still ends up being a better deal than R,
     * given that we don't have a fast backward sieve. */
    if (lower_limit >= upper_limit || correction > upper_limit-lower_limit)
      lower_limit = upper_limit;
    else
      lower_limit += correction;
    segment_size = lower_limit / 30;
    lower_limit = 30 * segment_size - 1;
    count = prime_count(lower_limit);

    /* printf("We've estimated %lu too %s.\n", (count>n)?count-n:n-count, (count>n)?"FAR":"little"); */
    /* printf("Our limit %lu %s a prime\n", lower_limit, is_prime(lower_limit) ? "is" : "is not"); */

    if (count >= n) { /* Too far.  Walk backwards */
      if (is_prime(lower_limit)) {
        if (count == n) return lower_limit;
        count--;
      }
      for (p = 0; p <= (count-n); p++)
        lower_limit = prev_prime(lower_limit);
      return lower_limit;
    }
    count -= 3;

    /* Make sure the segment siever won't have to keep resieving. */
    prime_precalc(isqrt(upper_limit));
  }

  if (count == target)
    return p;

  /* Start segment sieving.  Get memory to sieve into. */
  segbase = segment_size;
  segment = get_prime_segment(&segment_size);

  while (count < target) {
    /* Limit the segment size if we know the answer comes earlier */
    UV upperd = upper_limit / 30;
    if (segbase > upperd) break;  /* Something is very wrong.  Assert later. */
    if (segment_size > (upperd - segbase + 1))
      segment_size = upperd - segbase + 1;

    /* Do the actual sieving in the range */
    sieve_segment(segment, segbase, segbase + segment_size-1);

    /* Count up everything in this segment */
    count += count_segment_maxcount(segment, 30*segbase, segment_size, target-count, &p);

    if (count < target)
      segbase += segment_size;
  }
  release_prime_segment(segment);
  MPUassert(count == target, "nth_prime got incorrect count");
  return ( (segbase*30) + p );
}



/******************************************************************************/
/*                                   MISC                                     */
/******************************************************************************/


double ramanujan_axler(long double n, long double c, long double d) {
  long double res, U, c1, c2, log2 = logl(2), logn = logl(n), loglogn = logl(logn);

  c1 = 2*log2*log2 + log2 + c;
  c2 = log2*log2*log2 + 2*log2*log2 + d;

  U = (log2 * logn*loglogn*loglogn - c1*logn*loglogn + c2*logn - log2*log2*loglogn + log2*log2*log2 + log2*log2)
    / (logn*logn*logn*logn + logn*logn*logn*loglogn - logn*logn*logn*log2 - logn*logn*log2);

  res = 2*n * (1.0L + log2/logn - (log2*loglogn - log2*log2 - log2) / (logn*logn) + U);
  return res;
}