寻找快速素数计数函数

Question

寻找快速素数计数函数

我需要计算小于或等于某个 N 的素数个数，这就是素数计数函数或 PI 函数。我有这个，但运行速度太慢：

function PI(x) {
    var primes = 4;
    for (var i = 3; i <= x; i += 2) {
        if (i % 3 === 0 || i % 5 === 0 || i % 7 === 0) continue;
        var r = ~~Math.sqrt(i), p = true;
        for (var j = 2; j <= r; j++) {
            if (i % j === 0) {
                p = false;
                break;
            }
        }
        if (p)
            primes++;
    }
    return primes;
}

Run Code Online (Sandbox Code Playgroud)

计算时间PI(1000000)差不多一秒，PI(10000000)我的计算时间长达20秒。那太慢了。有更快的方法吗？

Answer 1 · 2017-04-16T09:16:24.000090+08:00

Gor*_*ood 7

我参加聚会的时间甚至更晚，但被问题标题中的“质数计数功能”触发，因为这并没有引起太多关注。

太长了；改性活生物体如何运作？跳到最后查看代码和结果...

素数计数函数技术的主要进步历史如下：

大约在 1830 年，阿德里安-玛丽·勒让德 (Adrien-Marie Legendre) 开发了一个公式，该公式使用“包含-排除原理”来计算下面的素数数量，而x无需枚举所有单个素数，只需要知道的平方根的基本素数值x。
1870 年，EDF Meissel 改进了勒让德公式，稍微减少了运算次数，但要求知道素数，x^(2/3)而不是x^(1/2)，其中“知道素数”通常由 SoE 进行，但只需要存储，而O(x^(1/3)/(ln x))不是O(x^(1/2)/(ln x))，尽管它只减少了的时间复杂度到O(x/((ln x)^3))from O(x/((ln x)^2))。
最后，Lagarias、Miller 和 Odzlinko (LMO) 在 1985 年开发了他们的算法，改进了 Meissel 方法，仍然需要与基本版本相同的“素值知识” x^(2/3)，但大大减少了计算“Phi”所需的时间：将计算分成两组计算（至少，更多用于以后的改进），一组计算使用递归整数除法通过通常的 Meissel 方法确定“Phi”函数的“普通”叶子，但这项工作可以忽略不计，因为只有x^(1/3)通过使用部分筛分和计数技术，可以大大减少计算其余“特殊”叶子的时间。可以使用部分筛选，因为“Phi(x, a)”的定义是以下所有数字的计数，x这些数字不能被该ath素数以下的任何素数整除，这正是通过筛选x最多的倍数所产生的结果基素值直到ath素数。基本 LMO 方法的计算复杂度为O(x^(2/3)(log (log x)))。

为了理解如何使用这些技术，让我们在将素数数到 64 的小数域中连续考虑它们，如下所示：

考虑 64 以内的素数，我们可以很容易地确定为 2、3、5、7、11、13、17、19、23、29、31、37、41、43、47、53、57 和 61，总共 18 个素数；对于本练习，考虑到对于勒让德技术（2、3、5 和 7 的四个素数）或基于 Meissel 的技术，我们只知道 64 的平方根为 8 的值或 64 的范围^(2/3) 是 16（2、3、5、7、11 和 13 的六个素数）。
Legendre 的技术说我们只需要计算“Phi(64, 4) + 4 - 1”的值，其中 4 是平方根以下的素数数量；“Phi”函数是不超过该限制的未排除值的数量 64，这些值被包括在内，因为它们不是小于限制平方根的四个基本质数的倍数，并且通过递归技术计算如上一个答案中所使用的，在本例中表示为64 - 64/2 - 64/3 + 64/2/3 - 64/5 + 64/5/2 + 64/5/3 - 64/5/2/3 - 64/7 + 64/7/2 + 64/7/3 - 64/7/2/3 + 64/7/5 - 64/7/5/2 - 64/7/5/3 + 64/7/5/2/3 = 15，当加 4 并减 1 时，给出正确答案 18。这称为“包含排除原则”，因为素数为奇数的项被排除/相减，并且包含/添加具有偶数个质因数的项。您可以看到，这需要大量的除法运算，尽管有一些技术（例如前面答案中使用的缓存）可以减少如此小的范围内的除法运算，但不幸的是只能达到极限。
对于 Meissel 技术，计算Phi(64, 2) + 2 - 12 是最大到 4 的立方根的素数的数量，计算起来要容易得多，因为64 - 64/2 - 64/3 + 64/2/3 = 21当加到 2 并减去 1 时，结果是 22，但这不是最终答案，因为我们需要减去“P2”项，即未排除的部分，因此基于素数 5 和 7，计算为最多 64/5 - 3 + 1 的素数个数加上 64/7 - 4 + 1，即 4。所以最终答案是 22 减去这 4 对，应该是 18。
LMO 技术的工作原理与上述完全相同，只是对于这个微不足道的限制，有一个“特殊”叶子，这是64/2/3因为 2 乘以 3 等于 6 超过 64 的立方根即 4 的叶子。如果将其实现为 LMO 将这样做，为了找到这个商的“Phi”，Phi(64/6, 0)我们根本不进行筛选，因为除数的最小质因数是 2，只需计算剩余的值数量直到商即可得到“Phi”的相同答案 10 “对于叶子来说。对于这样一个微不足道的限制，“特殊”叶子的数量很小，但重要的是特殊叶子的数量只会随着的增加而增加O((limit^(2/3))/(log (limit^(1/3)))^2)。由于进行筛分的工作O(limit^(2/3)(log log limit))较大，因此这是整体计算复杂性的控制因素。

那么我们如何实现基本的LMO算法呢？有以下几个步骤：

我们制作一个素数表，直到极限的立方根（大约(limit^(1/3)) / (ln limit)大小），这花费的时间可以忽略不计。
我们可以使用这些素数通过递归除法技术来计算“Phi”的“S1”普通部分，由于计算量很少，因此所花费的时间可以忽略不计O(limit^(1/3))。
O(limit^(1/3))我们从基本素数表（大小最大的表）中制作一个从素数计数到立方根的值的表；我们还制作了每个部分筛选步骤的中间“Phi”计数直至立方根基素值的表格，以及所确定的“S2”根叶参数的表格，这些参数实际上可以在“S1”时确定上面正在计算值；所有这些花费的时间都可以忽略不计，并且最大的（例如最后一个）仅与范围的立方根成正比，这就是中使用空间的原因O(limit^(1/3))。
我们按页段进行筛选，limit^(2/3)直到每个页一次被一个基本质数部分剔除的范围（高于根“S2”叶子表中找到的最小值），累积未剔除数量的计数总和值最多为每个“S2”叶商值。
基本的页分段仅赔率位打包 SoE 需要大约 80 秒来计算素数，直到达到 53 位限制的三分之二次方，但我们需要使用修改后的版本，为每个新剔除的数更新计数器表某些剔除操作的值，因此大约需要一半的时间或大约 120 秒。
然后，在此实现中，对每个“S2”点的未剔除值进行计数，使用批计数器数组，增量的最终范围很小，小于通过快速“弹出计数”技术计数的批大小。原始的 LMO 论文描述了使用二叉索引树来实现此目的，但没有尝试减少树深度，这也是它比此实现稍慢的部分原因，因为它增加了由于此成本log x而导致的剔除时间log x剔除时更新计数器的因素，这会使剔除时间慢很多倍；对于该实现的上述最大范围，“特殊”叶子的数量的计数和累加仅花费大约几十秒的时间。Kim Walisch对此成本进行了解释并提出了解决方案。
当对于给定页段完成了达到“S2”计算所需的每个单独限制的所有基本素数值时，则将页段完全剔除到小于最大值的平方的基本素数值段中表示的“P2”计数值和该段范围内的“P2”计数值从该完成的筛分页确定，而无需再次筛分相同的范围。
当达到限制的所有页段range^(2/3)都完成时，累加器将包含“S2”和“P2”所需的值，以便可以S1 + S2 - P2 + number of base primes to the cube root of the limit - 1将输出为该的素数计数的所需答案limit。

以下代码实现了基本的 LMO 算法（没有“alpha”优化，这会大大增加代码复杂性），“S2”特殊叶子的代码已松散但不完全从 Kim 的“primecount”C++ 代码翻译而来瓦利什：

"use strict";

const MAXVALUE = 9007199254740991; // 2**53 - 1

// creates a function returning a lazily memoized value from a thunk...
function lazy(thunk) {
  let value = undefined;
  return function() {
if (value === undefined) { value = thunk(); thunk = null; }
return value;
  }
}

// a page-segmented odds-only bit-packed Sieve of Eratosthenes;

const PGSZBITS = 262144; // about CPU l1 cache size in bits (power of two)

const CLUT = function () { // fast "pop count" Counting Look Up Table...
  const arr = new Uint8Array(65536);
  for (let i = 0; i < 65536; ++i) {
let nmbts = 0 | 0; let v = i;
while (v > 0) { ++nmbts; v &= (v - 1) | 0; }
arr[i] = nmbts | 0; }
  return arr;
}();

function countPageFromTo(bitstrt, bitlmt, sb) {
  const fst = bitstrt >> 5; const lst = bitlmt >> 5;
  const pg = new Uint32Array(sb.buffer);
  let v0 = (pg[fst] | ((0xFFFFFFFF >>> 0) << (bitstrt & 31))) >>> 0;
  let cnt = ((lst - fst) << 5) + CLUT[v0 & 0xFFFF]; cnt += CLUT[v0 >>> 16];
  for (let i = fst | 0; i < lst; ++i) {
let v = pg[i] >>> 0;
cnt -= CLUT[v & 0xFFFF]; cnt -= CLUT[v >>> 16];
  }
  let v1 = (pg[lst] | ((0xFFFFFFFE >>> 0) << (bitlmt & 31))) >>> 0;
  cnt -= CLUT[v1 & 0xFFFF]; cnt -= CLUT[v1 >>> 16]; return cnt | 0;
}

function partialSievePage(lwi, bp, sb) {
  const btsz = sb.length << 3;
  let s = Math.trunc((bp * bp - 3) / 2); // compute the start index...
  if (s >= lwi) s -= lwi; // adjust start index based on page lower limit...   
  else { // for the case where this isn't the first prime squared instance
let r = ((lwi - s) % bp) >>> 0;
s = (r != (0 >>> 0) ? bp - r : 0) >>> 0; }
  if (bp <= 32) {
for (let slmt = Math.min(btsz, s + (bp << 3)); s < slmt; s += bp) {
  const shft = s & 7; const msk = ((1 >>> 0) << shft) >>> 0;
  for (let c = s >> 3, clmt = sb.length; c < clmt | 0; c += bp)
    sb[c] |= msk; } }
  else
for (let slmt = sb.length << 3; s < slmt; s += bp)
  sb[s >> 3] |= ((1 >>> 0) << (s & 7)) >>> 0;
}

function partialSieveCountPage(lwi, bp, cntarr, sb) {
  const btsz = sb.length << 3; let cullcnt = 0;
  let s = Math.trunc((bp * bp - 3) / 2); // compute the start index...
  if (s >= lwi) // adjust start index based on page lower limit...
s -= lwi;
  else { // for the case where this isn't the first prime squared instance
let r = ((lwi - s) % bp) >>> 0;
s = (r != (0 >>> 0) ? bp - r : 0) >>> 0; }
  if (bp <= 32) {
for (let slmt = Math.min(btsz, s + (bp << 3)); s < slmt; s += bp) {
  const shft = s & 7; const msk = ((1 >>> 0) << shft) >>> 0;
  for (let c = s >>> 3, clmt = sb.length; c < clmt | 0; c += bp) {
    const isbit = ((sb[c] >>> shft) ^ 1) & 1;
    cntarr[c >> 6] -= isbit; cullcnt += isbit; sb[c] |= msk; }
}
  }
  else
for (let slmt = sb.length << 3; s < slmt; s += bp) {
  const sba = s >>> 3; const shft = s & 7;
  const isbit = ((sb[sba] >>> shft) ^ 1) & 1;
  cntarr[s >> 9] -= isbit; cullcnt += isbit;
  sb[sba] |= ((1 >>> 0) << shft) >>> 0; }
  return cullcnt;
}

// pre-culled pattern of small wheel primes...
const WHLPRMS = [ 2, 3, 5, 7, 11, 13, 17 ];
const WHLPTRNLEN = WHLPRMS.reduce((s, v) => s * v, 1) >>> 1; // odds only!
const WHLPTRN = function() { // larger than WHLPTRN by one buffer for overflow
  const len = (WHLPTRNLEN + (PGSZBITS >>> 3) + 3) & (-4); // up 2 even 32 bits!
  const arr = new Uint8Array(len);
  for (let bp of WHLPRMS.slice(1)) partialSievePage(0, bp, arr);
  arr[0] |= ~(-2 << ((WHLPRMS[WHLPRMS.length - 1] - 3) >> 1)) >>> 0; return arr;
}();

function fillPage(lwi, sb) {
  const mod = (lwi / 8) % WHLPTRNLEN;
  sb.set(new Uint8Array(WHLPTRN.buffer, mod, sb.length));
}

function cullPage(lwi, bpras, sb) {
  const btsz = sb.length << 3; let bp = 3;
  const nxti = lwi + btsz; // just beyond the current page 
  for (let bpra of bpras()) {
for (let bpri = 0; bpri < bpra.length; ++bpri) {
  const bpr = bpra[bpri]; bp += bpr + bpr;
  let s = (bp * bp - 3) / 2; // compute start index of prime squared
  if (s >= nxti) return; // enough bp's
  partialSievePage(lwi, bp, sb);
}
  }
}

function soePages(bitsz, bpras) {
  const buf =  new Uint8Array(bitsz >> 3); let lowi = 0;
  const gen = bpras === undefined ? makeBasePrimeRepArrs() : bpras;
  return function*() {
while (true) {
  fillPage(lowi, buf); cullPage(lowi, gen, buf);
  yield { lwi: lowi, sb: buf }; lowi += bitsz; }
  };
}

function makeBasePrimeRepArrs() {
  const buf = new Uint8Array(128); let gen = undefined; // avoid data race!
  fillPage(0, buf);
  for (let i = 8, bp = 19, sqr = bp * bp; sqr < 2048+3;
                                      ++i, bp += 2, sqr = bp * bp)
if (((buf[i >> 3] >>> 0) & ((1 << (i & 7)) >>> 0)) === 0 >>> 0)
  for (let c = (sqr - 3) >> 1; c < 1024; c += bp)
    buf[c >> 3] |= (1 << (c & 7)) >>> 0; // init zeroth buf
  function sb2bprs(sb) {
const btsz = sb.length << 3; let oi = 0;
const arr = new Uint8Array(countPageFromTo(0, btsz - 1, sb));
for (let i = 0, j = 0; i < btsz; ++i)
  if (((sb[i >> 3] >>> 0) & ((1 << (i & 7)) >>> 0)) === 0 >>> 0) {
    arr[j++] = (i - oi) >>> 0; oi = i; }
return { bpra: arr, lastgap: (btsz - oi) | 0 };
  }
  let { bpra, lastgap } = sb2bprs(buf);
  function next() {
const nxtpg = sb2bprs(gen.next().value.sb);
nxtpg.bpra[0] += lastgap; lastgap = nxtpg.lastgap;
return { head: nxtpg.bpra, tail: lazy(next) };
  }
  const lazylist = { head: bpra, tail: lazy(function() {
if (gen === undefined) {
  gen = soePages(1024)(); gen.next() } // past first page
return next();
  }) };
  return function*() { // return a generator of rep pages...
let ll = lazylist; while (true) {  yield ll.head; ll = ll.tail(); }
  };
}

function *revPrimesFrom(top, bpras) {
  const topndx = (top - 3) >>> 1;
  const buf = new Uint8Array(PGSZBITS >>> 3);
  let lwi = (((topndx / PGSZBITS) >>> 0) * PGSZBITS) >>> 0;
  let si = (topndx - lwi) >>> 0;
  for (; lwi >= 0; lwi -= PGSZBITS) { // usually external limit!
const base = 3 + lwi + lwi;
fillPage(lwi, buf); cullPage(lwi, bpras, buf);
for (; si >= 0 >>> 0; --si)
  if (((buf[si >> 3] >>> 0) & ((1 << (si & 7)) >>> 0)) === (0 >>> 0))
    yield base + si + si;
si = PGSZBITS - 1;
  }
};

const TinyPrimes = [ 2, 3, 5, 7, 11, 13, 17, 19 ]; // degree eight
const TinyProduct = TinyPrimes.reduce((s, v) => s * v) >>> 0;
const TinyHalfProduct = TinyProduct >>> 1;
const TinyTotient = TinyPrimes.reduce((s, v) => s * (v - 1), 1) >>> 0;
const TinyLength = (TinyProduct + 8) >>> 2; // include zero and half point!
const TinyTotients = function() {
  const arr = new Uint32Array(TinyLength | 0);
  arr[TinyLength - 1] = 1; // mark mid point value as not prime - never is
  let spn = 3 * 5 * 7; arr[0] = 1; // mark zeroth value as not prime!
  for (let bp of [ 3, 5, 7 ]) // cull small base prime values...
for (let c = (bp + 1) >>> 1; c <= spn; c += bp) arr[c] |= 1;
  for (let bp of [ 11, 13, 17, 19 ]) {
for (let i = 1 + spn; i < TinyLength; i += spn) {
  const rng = i + spn > TinyLength ? spn >> 1 : spn;
  arr.set(new  Uint32Array(arr.buffer, 4, rng), i); }
spn *= bp;
for (let c = (bp + 1) >>> 1; // eliminate prime in pattern!
       c < (spn > TinyLength ? TinyLength : spn + 1); c += bp)
  arr[c] |= 1;
  }
  arr.reduce((s, v, i) => { // accumulate sums...
const ns = s + (v ^ 1); arr[i] = ns; return ns; }, 0);
  return arr;
}();  

function tinyPhi(m) {
  const d = Math.trunc(m / TinyProduct);
  const ti = (m - d * TinyProduct + 1) >>> 1;
  const t = ti < TinyLength
          ? TinyTotients[ti]
          : TinyTotient - TinyTotients[TinyHalfProduct - ti];
  return d * TinyTotient + t;
}

function *countPrimesTo(limit) {
  if (limit <= WHLPRMS[WHLPRMS.length - 1]) {
let cnt = 0; for (let p of WHLPRMS) { if (p > limit) break; else ++cnt; }
return cnt; }

  const bpras = makeBasePrimeRepArrs();
  if (limit < 1024**2 + 3) { // for limit < about a million, just sieve...
let p = 3; let cnt = WHLPRMS.length;
for (let bpra of bpras())
  for (let bpr of bpra) { // just count base prime values to limit
    p += bpr + bpr; if (p > limit) return cnt; ++cnt; }
  }

  if (limit <= 32 * 2 * PGSZBITS + 3) { // count sieve to about 32 million...
const lmti = (limit - 3) / 2;
let cnt = WHLPRMS.length; // just use page counting to limit as per usual...
for (let pg of soePages(PGSZBITS, bpras)()) {
  const nxti = pg.lwi + (pg.sb.length << 3);
  if (nxti > lmti) { cnt += countPageFromTo(0, lmti - pg.lwi, pg.sb); break; }
  cnt += countPageFromTo(0, PGSZBITS - 1, pg.sb);
}
return cnt;
  }

  // Actual LMO prime counting code starts here...
  const sqrt = Math.trunc(Math.sqrt(limit)) >>> 0;
  const cbrt = Math.trunc(Math.cbrt(limit)) >>> 0;
  const sqrtcbrt = Math.trunc(Math.sqrt(cbrt)) >>> 0;
  const top = cbrt * cbrt - 1; // sized for maximum required!
  const bsprms = function() {
let bp = 3; let cnt = WHLPRMS.length + 1; for (let bpra of bpras())
  for (let bpr of bpra) {
    bp += bpr + bpr; if (bp > cbrt) return new Uint32Array(cnt); ++cnt; }
  }();
  bsprms.set(WHLPRMS, 1); // index zero not used == 0!
  const pisqrtcbrt = function() {
let cnt = WHLPRMS.length; let i = cnt + 1; let bp = 3;
stop: for (let bpra of bpras())
  for (let bpr of bpra) {
    bp += bpr + bpr; if (bp > cbrt) break stop;
    if (bp <= sqrtcbrt) ++cnt; bsprms[i++] = bp >>> 0; }
return cnt;
  }();
  const pis = function() { // starts with index 0!
const arr = new Uint32Array(cbrt + 2); let j = 0;
for (let i = 1; i < bsprms.length; ++i)
  for (; j < bsprms[i]; ) arr[j++] = (i - 1) >>> 0;
for (; j < arr.length; ) arr[j++] = (bsprms.length - 1) >>> 0;
return arr;
  }();
  const phis = function() { // index below TinyPhi degree never used...
const arr = (new Array(bsprms.length)).fill(1);
arr[0] = 0; arr[1] = 3; arr[2] = 2; // unused
for (let i = WHLPRMS.length + 2; i < arr.length; ++i) {
  arr[i] -= i - WHLPRMS.length - 1; } // account for non phi primes!
return arr;
  }();
  // indexed by `m`, contains greatest prime factor index and
  // Moebius value bit; for Moebius, a one is negative...
  const specialroots = new Uint16Array(cbrt + 1); // filled in with S1 below...
  const S1 = function() { // it is very easy to calculate S1 recursively...
let s1acc = tinyPhi(limit);
function level(lmtlpfni, mfv, m) {
  for (let lpfni = 9; lpfni < lmtlpfni; ++lpfni) {
    const pn = bsprms[lpfni]; const nm = m * pn;
    if (nm > cbrt) { // don't split, found S2 root leaf...
      specialroots[m] = (lmtlpfni << 1) | (mfv < 0 ? 1 : 0); return; }
    else { // recurse for S1; never more than 11 levels deep...
      s1acc += mfv * tinyPhi(Math.trunc(limit / nm)); // down level...
      level(lpfni, -mfv, nm); // Moebius sign change on down level!
    } // up prime value, same level!
  }
}
level(bsprms.length, -1, 1); return s1acc;
  }();

  // at last, calculate the more complex parts of the final answer:
  function *complex() {
let s2acc = 0; let p2acc = 0; let p2cnt = 0; // for "P2" calculation
const buf = new Uint8Array(PGSZBITS >>> 3); let ttlcnt = 0;
const cnts = new Uint8Array(PGSZBITS >>> 9);
const cntaccs = new Uint32Array(cnts.length);
const revgen = revPrimesFrom(sqrt, bpras);
let p2v = Math.trunc(limit / revgen.next().value);
const lwilmt = Math.trunc((top - 3) / 2);

const updtmsk = ((PGSZBITS << 3) - 1) >>> 0;
let nxttm = Date.now() + 1000;
for (let lwi = 0; lwi <= lwilmt; lwi += PGSZBITS) { // for all pages
  if ((lwi & updtmsk) == updtmsk && Date.now() >= nxttm) {
    nxttm = Date.now() + 1000; yield lwi / lwilmt * 100.0 };
  let pgcnt = 0; const low = 3 + lwi + lwi;
  const high = Math.min(low + (PGSZBITS << 1) - 1, top);
  let cntstrti = 0 >>> 0;
  function countTo(stop) {
    const cntwrd = stop >>> 9; const bsndx = stop & (-512);
    const xtr = countPageFromTo(bsndx, stop, buf);
    while (cntstrti < cntwrd) {
      const ncnt = cntaccs[cntstrti] + cnts[cntstrti];
      cntaccs[++cntstrti] = ncnt; }
    return cntaccs[cntwrd] + xtr;
  }
  const bpilmt = pis[Math.trunc(Math.sqrt(high)) >>> 0] >>> 0;
  const maxbpi = pis[Math.min(cbrt, Math.trunc(Math.sqrt(limit/low)))]>>>0;
  const tminbpi = pis[Math.min(Math.trunc(top / (high + 1)),
                               bsprms[maxbpi]) >>> 0];
  const minbpi = Math.max(TinyPrimes.length, tminbpi) + 1;
  fillPage(lwi, buf); let bpi = (WHLPRMS.length + 1) >>> 0;
 
  if (minbpi <= maxbpi) { // jump to doing P2 if not range

    // for bpi < minbpi there are no special leaves...
    for (; bpi < minbpi; ++bpi) { // eliminate all Tiny Phi primes...
      const bp = bsprms[bpi]; const i = (bp - 3) >>> 1; // cull base primes!
      phis[bpi] += countPageFromTo(0, PGSZBITS - 1, buf);
      partialSievePage(lwi, bp, buf); }
    for (let i = 0; i < cnts.length; ++i) { // init cnts arr...
      const s = i << 9; const c = countPageFromTo(s, s + 511, buf);
      cnts[i] = c; pgcnt += c; }

    // for all base prime values up to limit**(1/6) in the page,
    // add all special leaves composed of this base prime value and
    // any number of other higher base primes, all different,
    // that qualify as special leaves...
    let brkchkr = false;
    for (; bpi <= Math.min(pisqrtcbrt, maxbpi) >>> 0; ++bpi) {
      const bp = bsprms[bpi];
      const minm = Math.max(Math.trunc(limit / (bp * (high + 1))),
                            Math.trunc(cbrt / bp)) >>> 0;
      const maxm = Math.min(Math.trunc(limit / (bp * low)), cbrt) >>> 0;
      if (bp >= maxm) { brkchkr = true; break; }
      for (let m = maxm; m > minm; --m) {
        const rt = specialroots[m];
        if (rt != 0 && bpi < rt >>> 1) {
          const stop = Math.trunc(limit / (bp * m) - low) >>> 1;
          const mu = ((rt & 1) << 1) - 1; // one bit means negative!
          s2acc -= mu * (phis[bpi] + countTo(stop));
        } }
      phis[bpi] += pgcnt; // update intermediate base prime counters
      pgcnt -= partialSieveCountPage(lwi, bp, cnts, buf);
      cntstrti = 0; cntaccs[0] = 0;
    }
    // for all base prime values > limit**(1/6) in the page,
    // add results of all special levels composed using only two primes...
    if (!brkchkr)
    for (; bpi <= maxbpi; ++bpi) {
      const bp = bsprms[bpi];
      let l = pis[Math.min(Math.trunc(limit / (bp * low)), cbrt)>>>0]>>>0;
      if (bp >= bsprms[l]) break;
      const piminm = pis[Math.max(Math.trunc(limit / (bp * (high + 1))),
                                  bp) >>> 0] >>> 0;
      for (; l > piminm; --l) {          
        const stop = Math.trunc(limit / (bp * bsprms[l]) - low) >>> 1;
        s2acc += phis[bpi] + countTo(stop);
      }

归档时间：	9 年，3 月前
查看次数：	1028 次
最近记录：	3 年，7 月前