3

我最近读到了一个更快的 Eratosthenes 分段筛的实现,用于非常大的数字。

以下是相同的实现:

function sieve(low, high) {

    var primeArray = [], ll = Math.sqrt(low), output = [];

    for (var i = 0; i < high; i++) {
        primeArray[i] = true;
    }

    for (var i = 2; i <= ll; i++) {
        if (primeArray[i]) {
            for (var j = i * i; j < high; j += i) {
                primeArray[j] = false;
            }
        }
    }

    for (var i = 2; i < ll; i++) {
        if(primeArray[i])
        {
            var segmentStart = Math.floor(low/i) * i;

            for(var j = segmentStart; j <= high; j+=i)
            {
                primeArray[j] = false;
            }
        }
    }

    for(var i = low; i <= high; i++)
    {
        if(primeArray[i])
        {
            output.push(i);
        }
    }

    return output;
};

我似乎无法弄清楚我在哪里弄错了。可能是工作太久了。

例如:

sieve(4,10)应该返回[5,7] 但它正在返回[5,7,9]

4

3 回答 3

15

尽管您从阅读中了解到,埃拉托色尼的页面分段筛是一种在大范围内查找素数的快速方法,但您的问题代码(即使已更正)并未实现页面分段 SoE,请在大范围内测试代码,也不是像 SoE 实现那样特别快。以下讨论将展示如何在大范围内使用真正的页面分段 SoE。

概要

以下是导致您的意图的越来越快的实施的阶段性进展,并附有解释每个步骤的原因和实施细节的评论。它包括 JavaScript 中的可运行片段,但这些技术不仅限于 JavaScript,其他语言也不会限制一些进一步的改进,例如结果页面的多线程(Web Workers 除外,这很难控制,因为处理顺序),极端循环展开中的一些进一步优化,最后一个与代码效率有限有关,因为必须由浏览器中的 JavaScript 引擎及时(JIT)编译为本机代码;这些限制与直接编译为非常高效的本机代码的语言相比,例如 C/C++、Nim、Rust、Free Pascal、Haskell、Julia、

第 1 章 - 设置基线

首先,让我们从您当前代码算法的工作版本开始,该算法在相当大的范围内使用时序信息来建立基线;下面的代码在剔除素数的平方处开始剔除每个素数,这避免了剔除给定素数值和一些冗余起始剔除的问题,并且没有理由为我们可以产生的大范围生成结果素数的输出数组直接来自剔除数组的素数;此外,答案的确定超出了时间范围,因为我们将开发更好的技术来找到一个“答案”,对于大范围来说,答案通常是找到的素数的计数、素数的总和、素数的第一次出现间隙等,都不需要实际查看找到的素数:

"use strict";

function soePrimesTo(limit) {
  var sz = limit - 1;
  var cmpsts = new Uint8Array(sz); // index 0 represents 2; sz-1 is limit
  // no need to zero above composites array; zeroed on creation...
  for (var p = 2; ; ++p) {
    var sqr = p * p;
    if (sqr > limit) break; // while p is the square root of limit -> cull...
    if (cmpsts[p - 2] == 0 >>> 0) // 0/1 is false/true; false means prime...
      for (var c = sqr - 2; c < sz; c += p) // set true for composites...
          cmpsts[c] = 1 >>> 0; // use asm.js for some extra efficiency...
  }
  var bi = 0
  return function () {
    while (bi < sz && cmpsts[bi] != 0 >>> 0) ++bi;
    if (bi >= sz) return null;
    return bi++ + 2;
  };
}

// show it works...
var gen = soePrimesTo(100);
var p = gen();
var output = [];
while (p != null) { output.push(p); p = gen(); }
console.log("Primes to 100 are:  " + output + ".");

var n = 1000000000; // set the range here...
var elpsd = -new Date();
gen = soePrimesTo(n);
elpsd += +new Date();
var count = 0;
while (gen() != null) { count++; }
console.log("Found " + count + " primes up to " + n + " in " + elpsd + " milliseconds.");

现在,引用我的上述代码筛选到十亿的运行时间几乎没有什么目的,因为你的时间几乎肯定会更快,因为我在英特尔 x5-Z8350 中使用非常低端的平板电脑 CPU,运行频率为 1.92 Gigahertz(CPU 时钟单个活动线程的速度)。我将仅引用一次我的执行时间作为如何计算每次剔除的平均 CPU 时钟周期的示例:我将上述代码的执行时间约 43350 毫秒为 43.35 秒乘以每秒 19.2 亿时钟除以约 2.514十亿次剔除操作筛选到十亿次(可以从公式或从SoE 维基百科页面上的无轮分解表计算),每次剔除大约33.1 个 CPU 时钟周期到十亿。从现在开始,我们将只使用每次剔除的 CPU 时钟周期来比较性能。

进一步注意,这些性能分数很大程度上取决于所使用的浏览器 JavaScript 引擎,上述分数是在 Google Chrome(版本 72)上运行的;Microsoft Edge(版本 44)比我们正在转向的 Page Segmented SoE 版本慢了大约 7 倍,而 Firefox 和 Safari 的性能可能接近 Chrome。

由于使用了 TypedArray 和更多 asm.js,此性能可能比之前的答案代码更好,Uint8Array但这种“一个巨大的数组筛子”(此处使用的 GigaByte 内存)的时间由于内存访问速度而受到瓶颈超出 CPU 缓存限制的主 RAM 内存。这就是为什么我们正在开发 Page Segmented Sieve 的原因,但首先让我们做一些事情来减少使用的内存量和所需的剔除周期数。

第 2 章 - 位打包和仅奇数轮分解

下面的代码进行了位打包,这需要在紧密的内部剔除循环中稍微复杂一些,但由于它只使用每个复合数的一位,因此它将内存使用量减少了八倍;同样,由于 2 是唯一的偶数素数,它使用基本的轮子分解来筛选赔率,只是为了将内存使用量进一步减少 2 倍,并将剔除操作的数量减少约 2.5 倍。

这种仅赔率的最小轮式分解的工作原理如下:

  1. 我们做了一个“轮子”,上面有两个位置,一半“命中”偶数,另一半“命中”奇数;因此,当我们将这个“轮子”“滚动”在所有潜在的素数上时,这个“轮子”的圆周跨度为两个数字,但它只“窗口”了它“滚动”的数字的二分之一或一半。
  2. 然后,我们将我们“滚动”的所有数字分成两个连续值的位平面,在一个位平面中,我们将所有偶数从 4 开始,在另一个位平面中,我们将所有赔率从 3 开始。
  3. 现在我们丢弃偶数平面,因为没有一个表示的数字可以是素数。
  4. 考虑的奇数基素数的起始索引p * p始终是奇数,因为奇数乘以奇数始终是奇数。
  5. 当我们将索引添加到奇数位平面中的素数时,我们实际上添加了两倍的值,因为添加一个奇数和一个奇数会产生一个偶数,这将在我们丢弃的位平面中,所以我们添加了素数再次回到奇数位平面。
  6. 奇数位平面索引位置自动考虑了这一点,因为我们已经丢弃了之前在每个奇数位位置索引之间的所有偶数值。
  7. 这就是为什么我们可以通过在下面的代码中重复向索引添加素数来进行剔除。

"use strict";

function soeOddPrimesTo(limit) {
  var lmti = (limit - 3) >> 1; // bit index for limit value
  var sz = (lmti >> 3) + 1; // size in bytes
  var cmpsts = new Uint8Array(sz); // index 0 represents 3
  // no need to zero above composites array; zeroed on creation...
  for (var i = 0; ; ++i) {
    var p = i + i + 3; // the square index is (p * p - 3) / 2 but we
    var sqri = (i << 1) * (i + 3) + 3; // calculate start index directly 
    if (sqri > lmti) break; // while p is < square root of limit -> cull...
    // following does bit unpacking to test the prime bit...
    // 0/1 is false/true; false means prime...
    // use asm.js with the shifts to make uint8's for some extra efficiency...
    if ((cmpsts[i >> 3] & ((1 >>> 0) << (i & 7))) == 0 >>> 0)
      for (var c = sqri; c <= lmti; c += p) // set true for composites...
          cmpsts[c >> 3] |= (1 >>> 0) << (c & 7); // masking in the bit
  }
  var bi = -1
  return function () { // return function to return successive primes per call...
    if (bi < 0) { ++bi; return 2 } // the only even prime is a special case
    while (bi <= lmti && (cmpsts[bi >> 3] & ((1 >>> 0) << (bi & 7))) != 0 >>> 0) ++bi;
    if (bi > lmti) return null; // return null following the last prime
    return (bi++ << 1) + 3; // else calculate the prime from the index
  };
}

// show it works...
var gen = soeOddPrimesTo(100);
var p = gen();
var output = [];
while (p != null) { output.push(p); p = gen(); }
console.log("Primes to 100 are:  " + output + ".");

var n = 1000000000; // set the range here...
var elpsd = -new Date();
gen = soeOddPrimesTo(n);
elpsd += +new Date();
var count = 0;
while (gen() != null) { count++; }
console.log("Found " + count + " primes up to " + n + " in " + elpsd + " milliseconds.");

现在,性能约为34.75 CPU 时钟周期,对于 10.257 亿次剔除操作,仅在 10 亿范围内(来自 Wikipedia),这意味着时间的减少几乎完全是由于剔除操作数量的减少,由于内存使用减少了 16 倍,“位旋转”位打包的额外复杂性只需要与节省的时间相同的额外时间。因此,这个版本使用了十六分之一的内存,速度提高了大约 2.5 倍。

但是我们还没有完成,正如您的消息来源告诉您的那样,页面分割可以更快地加快我们的速度。

第 3 章 - 包括更快素数计数的页面分割

那么什么是页面分割应用于 SoE 以及它对我们有什么作用?

页面分割是将筛选工作从一次筛选的一个巨大数组划分为一系列连续筛选的较小页面。然后它需要更多的复杂性,因为必须有一个单独的基本素数流可用,这可以通过使用内部筛子递归筛选来获得,该内部筛子生成要由主筛子使用的记忆基本素数列表。同样,输出结果的生成通常稍微复杂一些,因为它涉及对每个生成的筛选页面进行连续扫描和缩减。

页面分割有很多好处,如下所示:

  1. 它进一步降低了内存需求。使用以前的纯概率代码筛选到十亿需要大约 64 兆字节的 RAM,但不能使用该算法筛选超过 16 到 320 亿的范围(即使我们想等待相当长的时间)由于使用 JavaScript 可用的所有内存。通过(例如)该数量的平方根的页面大小进行筛选,我们可以筛选到该数量的平方或超出我们有时间等待的任何范围。我们还必须将基本素数存储到所需范围的平方根,但这对于我们想要考虑的任何范围都只有几十兆字节。
  2. 它提高了内存访问速度,这对每次剔除操作的 CPU 时钟周期数有直接影响。如果剔除操作主要发生在 CPU 缓存中,则内存访问会从每次访问主 RAM 内存的数十个 CPU 时钟周期(取决于 CPU 以及 RAM 的质量和类型)变为 CPU L1 缓存的大约一个时钟周期(较小大约 16 KB)到 CPU L2 缓存的大约 8 个时钟周期(更大至少大约 128 KB),我们可以制定剔除算法,以便我们充分利用所有缓存,使用小型快速 L1 缓存大多数具有小基本素数的操作,对于中等大小的基本素数使用较大的一点点慢速二级缓存,并且仅使用主 RAM 访问进行少数操作,使用大型基本素数进行非常大的范围。
  3. 它通过将筛选每个较大页面的工作分配给不同的 CPU 内核以实现相当粗粒度的并发性,开启了多线程的可能性,尽管这不适用于 JavaScript,除非通过使用 Web Workers(混乱)。

除了增加的复杂性之外,页面分段还有另一个问题需要解决:与“一个巨大的数组”筛不同,其中开始索引很容易计算一次然后用于整个数组,分段筛要求起始地址通过以下方式计算每页每个素数的模(除)运算(计算成本高),或者需要使用额外的内存来存储每页每个基本素数达到的索引,因此不必重新计算起始索引,但最后一种技术排除了多线程,因为这些数组不同步。将在 Ultimate 版本中使用的最佳解决方案是结合使用这些技术,其中几个页面段被分组以形成一个相当大的线程工作单元,以便这些地址计算占用总时间的可接受的一小部分,并且索引存储表用于每个线程的这些较大工作单元的基本素数,以便每个较大的工作单元只需进行一次复杂的计算。因此,我们既获得了多线程的可能性,又减少了开销。然而,下面的代码并没有减少这种开销,它筛分到十亿时,成本约为 10% 至 20% 。

除了页面分割之外,以下代码通过使用计数查找表 (CLUT) 人口计数算法增加了对找到的素数的有效计数,该算法一次计数 32 位,以便连续查找结果的开销发现素数只占筛分时间的一小部分。如果不这样做,枚举单个找到的素数只是为了确定有多少,至少需要在给定范围内进行筛选所需的时间。可以轻松开发类似的快速例程来执行诸如求和找到的素数、查找素数间隙等事情。

START_EDIT:

下面的代码增加了另一个加速:对于较小的素数(这种优化是有效的),代码通过识别剔除操作遵循八步模式来执行一种循环分离形式。这是因为一个字节的位数是偶数,我们通过奇数素数进行剔除,每八次剔除将返回到一个字节中的相同位位置;这意味着对于每个位位置,我们可以简化内部剔除循环以屏蔽一个常量位,从而大大简化内部循环并使剔除速度提高大约两倍,因为模式中的每个剔除循环不需要执行“位旋转”位打包操作。这一变化节省了大约 35%的执行时间,筛选到十亿。可以通过更改来禁用640. 它还为八循环的本机代码极端展开奠定了基础,因为这种模式可以在使用本机代码编译器时将剔除操作速度提高大约两倍。

进一步的小修改通过对掩码值使用查找表 (LUT) 而不是左移操作来使较大素数(大于 8192)的循环更快,从而在每次剔除操作时平均节省大约一半的 CPU 时钟周期剔除十亿的范围;随着范围从十亿增加,这种节省会略有增加,但在 JavaScript 中不是那么有效,因此已被删除。

END_EDIT

ANOTHER_EDIT:

除了上述编辑之外,我们还删除了 LUT BITMASK,但现在通过从相同大小的零缓冲区快速复制字节将筛缓冲区归零,并添加了计数 LUT 人口计数技术,总体速度提高了约 10%。

END_ANOTHER_EDIT

最终编辑:

使筛缓冲区大小约为 CPU L2 缓存大小而不是 L1,因为剔除循环速度永远不会受到 L2 缓存访问速度的限制。这导致速度略有提高,航程增加 64 倍,同时保持效率。

END_FINAL_EDIT

// JavaScript implementation of Page Segmented Sieve of Eratosthenes...
// This takes almost no memory as it is bit-packed and odds-only,
// and only uses memory proportional to the square root of the range;
// it is also quite fast for large ranges due to imrproved cache associativity...

"use strict";

const PGSZBYTES = 16384 * 8;
const PGSZBITS = PGSZBYTES * 8;
const ZEROSPTRN = new Uint8Array(PGSZBYTES);

function soePages(bitsz) {
  let len = bitsz >> 3;
  let bpa = [];
  let buf =  new Uint8Array(len);
  let lowi = 0;
  let gen;
  return function () {
let nxt = 3 + ((lowi + bitsz) << 1); // just beyond the current page
buf.set(ZEROSPTRN.subarray(0,buf.length)); // clear sieve array!
if (lowi <= 0 && bitsz < 131072) { // special culling for first page as no base primes yet:
  for (let i = 0, p = 3, sqr = 9; sqr < nxt; ++i, p += 2, sqr = p * p)
    if ((buf[i >> 3] & (1 << (i & 7))) === 0)
      for (let j = (sqr - 3) >> 1; j < 131072; j += p)
        buf[j >> 3] |= 1 << (j & 7);
} else { // other than the first "zeroth" page:
  if (!bpa.length) { // if this is the first page after the zero one:
    gen = basePrimes(); // initialize separate base primes stream:
    bpa.push(gen()); // keep the next prime (3 in this case)
  }
  // get enough base primes for the page range...
  for (let p = bpa[bpa.length - 1], sqr = p * p; sqr < nxt;
       p = gen(), bpa.push(p), sqr = p * p);
  for (let i = 0; i < bpa.length; ++i) { // for each base prime in the array
    let p = bpa[i] >>> 0;
    let s = (p * p - 3) >>> 1; // compute the start index of the prime squared
    if (s >= lowi) // adjust start index based on page lower limit...
      s -= lowi;
    else { // for the case where this isn't the first prime squared instance
      let r = (lowi - s) % p;
      s = (r != 0) ? p - r : 0;
    }
    if (p <= 32) {
      for (let slmt = Math.min(bitsz, s + (p << 3)); s < slmt; s += p) {
        let msk = ((1 >>> 0) << (s & 7)) >>> 0;
        for (let c = s >>> 3, clmt = bitsz >= 131072 ? len : len; c < clmt | 0; c += p)
          buf[c] |= msk;
      }
    }
    else
      // inner tight composite culling loop for given prime number across page
      for (let slmt = bitsz >= 131072 ? bitsz : bitsz; s < slmt; s += p)
        buf[s >> 3] |=  ((1 >>> 0) << (s & 7)) >>> 0;
  }
}
let olowi = lowi;
lowi += bitsz;
return [olowi, buf];
  };
}

function basePrimes() {
  var bi = 0;
  var lowi;
  var buf;
  var len;
  var gen = soePages(256);
  return function () {
while (true) {
  if (bi < 1) {
    var pg = gen();
    lowi = pg[0];
    buf = pg[1];
    len = buf.length << 3;
  }
  //find next marker still with prime status
  while (bi < len && buf[bi >> 3] & ((1 >>> 0) << (bi & 7))) bi++;
  if (bi < len) // within buffer: output computed prime
    return 3 + ((lowi + bi++) << 1);
  // beyond buffer range: advance buffer
  bi = 0;
  lowi += len; // and recursively loop to make a new page buffer
}
  };
}

const CLUT = function () {
  let 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 countPage(bitlmt, sb) {
  let lst = bitlmt >> 5;
  let pg = new Uint32Array(sb.buffer);
  let cnt = (lst << 5) + 32;
  for (let i = 0 | 0; i < lst; ++i) {
let v = pg[i]; cnt -= CLUT[v & 0xFFFF]; cnt -= CLUT[v >>> 16];
  }
  var n = pg[lst] | (0xFFFFFFFE << (bitlmt & 31));
  cnt -= CLUT[n & 0xFFFF]; cnt -= CLUT[n >>> 16];
  return cnt;
}

function countSoEPrimesTo(limit) {
  if (limit < 3) {
if (limit < 2) return 0;
return 1;
  }
  var cnt = 1;
  var lmti = (limit - 3) >>> 1;
  var lowi;
  var buf;
  var len;
  var nxti;
  var gen = soePages(PGSZBITS);
  while (true) {
var pg = gen();
lowi = pg[0];
buf = pg[1];
len = buf.length << 3;
nxti = lowi + len;
if (nxti > lmti) {
  cnt += countPage(lmti - lowi, buf);
  break;
}
cnt += countPage(len - 1, buf);
  }
  return cnt;
}

var limit = 1000000000; // sieve to this limit...
var start = +new Date();
var answr = countSoEPrimesTo(limit);
var elpsd = +new Date() - start;
console.log("Found " + answr + " primes up to " + limit + " in " + elpsd + " milliseconds.");

正如这里实现的那样,该代码每次筛选到十亿需要大约13.0 个 CPU 时钟周期。由于有效页面大小增加了 105 倍,因此在使用以下最大轮分解算法时,通过改进地址计算算法节省的额外约 20% 会自动得到改进,因此该开销仅变为百分之几,与百分之几相当用于数组填充和结果计数。

第 4 章 - 添加最大车轮分解的进一步高级工作

现在我们考虑使用最大轮分解的更广泛的变化(不仅对于“仅赔率”使用 2,而且对于覆盖 210 个潜在素数的轮,而不是仅仅 2 ) 并在初始化小筛选数组时进行预剔除,这样就不必通过以下素数 11、13、17 和 19 进行剔除。这减少了使用页面分段筛时的合数剔除操作的数量大约四倍到十亿的范围(如表格中所示/根据维基百科文章中的公式计算 - “组合轮”)并且可以编写成它的运行速度大约是四倍由于减少每个剔除操作的操作与上述代码的速度大致相同。

有效地进行 210 跨度车轮分解的方法是遵循“仅奇数”方法:上面的当前算法可以被认为是从两个中筛选出一个位压缩平面,如前一章所述,其中另一个平面可以被消除,因为它只包含两个以上的偶数;对于 210 跨度,我们可以定义 48 个这种大小的位压缩数组,表示 11 及以上的可能素数,其中所有其他 162 个平面包含的数字是 2、3、5 或 7 的因数,因此不需要被考虑。每个位平面可以通过以基素数为增量重复索引来剔除,就像奇数平面是通过结构自动处理的乘法一样完成的,就像仅对奇数一样。

虽然上面描述的扩展代码有几百行而且很长,但它可以在我的低端 Intel 1.92 Gigahertz CPU 上使用 Google V8 JavaScript 引擎在大约两秒内将素数的数量计算到十亿,大约比在本机代码中运行的相同算法慢五倍。

尽管上面的代码在大约 160 亿的范围内效率很高,但其他改进可以帮助将效率保持在更大的范围内,例如 1e14 或更多的数万亿。我们通过向上调整有效页面大小来实现这一点,以使它们永远不会小于被筛选的整个范围的平方根,但是对于小素数,以 16 KiloByte 块递增地筛选它们,对于中等素数,以 128 KiloByte 块递增筛选,并且只筛选根据我们的基线实现,一个巨大的数组,用于最大基本素数大小的极少数剔除操作。通过这种方式,对于我们可能考虑的最大范围,我们每次剔除的时钟不会增加超过大约两倍的小倍数。

由于这个答案接近 30,000 个字符的限制大小,因此在我的后续第 4.5a 章和(未来)第 4.5b 章对上述技术的实际实现的回答中继续对最大轮分解进行进一步讨论。

第 5 章 - JavaScript(和其他 VM 语言)不能做什么

对于 JavaScript 和其他虚拟机语言,最小剔除循环时间约为每个剔除循环 10 个 CPU 周期,并且变化不大。这使用 C/C++、Nim、Rust、Free Pascal、Haskell、Julia、等等

此外,还有一些极端的循环展开技术可以与至少其中一些语言一起使用,这些语言可以将平均剔除操作周期减少大约两倍,这是 JavaScript 所拒绝的

多线程可以通过使用有效 CPU 内核的因子来减少执行时间,但是使用 JavaScript 实现这一点的唯一方法是使用 Web Workers,而且同步起来很麻烦。在我的机器上,我有四个内核,但由于 CPU 时钟频率降低到四分之三且所有内核都处于活动状态,因此速度仅提高了三倍;对于 JavaScript 来说,这三倍是不容易获得的

所以这是关于使用 JavaScript 和其他当前 VM 语言的最新技术的限制,除了它们可以轻松使用多线程之外,它们具有大致相同的限制,以上因素的组合意味着本机代码编译器可以大约 20 倍比 JavaScript 更快(包括多线程,甚至在具有大量内核的新 CPU 中更快)。

但是,我相信三到五年内 Web 编程的未来将是 Web Assembly,它有可能克服所有这些限制。它现在非常接近支持多线程,虽然目前在 Chrome 上这个算法只比 JavaScript 快约 30%,但当使用某些 Web 从某些语言编译时,它只比当前某些浏览器上的本机代码慢一点汇编编译器。对于 Web Assembly 的高效编译器和对本机代码的高效浏览器编译仍处于早期开发阶段,但由于 Web Assembly 比大多数 VM 更接近本机代码,因此可以很容易地改进它以生成与本机代码一样快或接近的本机代码与其他新代码编译语言的代码一样快。

但是,除了将 JavaScript 库和框架编译成 Web Assembly 之外,我不相信 Web 的未来将是 JavaScript 到 Web Assembly 编译器,而是从其他一些语言编译。对于 Web 编程的未来,我最喜欢的选择是 F#,它可能将 Fable 实现转换为生成 Web 程序集,而不是 JavaScript (asm.js) 或 Nim。甚至有可能生成 Web Assembly,它支持并显示极端循环展开的优势,非常接近已知的最快的 Page Segmented SoE 速度。

结论

我们在 JavaScript 中构建了一个页面分段的 Eratosthenes 筛选器,适用于筛选数十亿的大范围,并且有进一步扩展这项工作的方法。生成的代码的剔除操作减少了大约十倍(当完全分解车轮时),剔除操作快了大约三倍,这意味着每个给定(大)范围的代码大约快 30 倍,而减少的内存使用意味着可以筛选到大约 9e15 的 53 位尾数范围在大约一年的时间内(只需打开 Chrome 浏览器选项卡并备份电源)。

尽管还有一些其他可能的小调整,但这是关于使用 JavaScript 筛选素数的最新技术:虽然由于给定的原因它不像本机代码那样快,但它足够快找到 1e14 的素数数量在一两天内(即使在中端智能手机上),只需打开浏览器选项卡即可获得所需的计算时间;这是相当惊人的,因为直到 1985 年才知道这个范围内的素数数量,然后通过使用数值分析技术而不是使用筛子,因为当时的计算机使用最快的编码技术还不够快在合理和经济的时间内。尽管我们可以在短短几个小时内使用这些算法来实现现代台式计算机上最好的本机代码编译器,

于 2019-04-19T11:17:22.677 回答
11

这扩展了我之前的答案,以添加承诺但每个答案限制 30,000 个字符中没有空间的内容:

第 4.5a 章 - 全轮因式分解的基础

上一个答案中第 3 章的 Eratosthenes 版本的非最大轮因子页分段筛被编写为素数生成器,其输出递归地反馈作为基本素数馈送的输入;虽然这非常优雅和可扩展,但在接下来的工作中,我已经退回到更命令式的代码风格,以便读者更容易理解最大轮分解的核心概念。在未来的第 4.5b 章中,我会将以下示例中开发的概念重新组合到 Prime 生成器样式中,并添加一些额外的改进,这些改进不会使其在数十亿的较小范围内变得更快,但会使该概念在没有大部分速度损失高达数万亿到数百或数千万亿;

以下示例的主要额外改进在于用于有效寻址轮模残差的各种查找表 (LUT),特殊起始地址 LUT 的生成非常简单地允许计算每个模残差位的剔除起始地址平面给出了整个结构中第一个剔除的起始地址轮索引和第一个模剩余位平面索引,没有额外的整数除法(慢),以及使用这些的 Sieve Buffer 复合数表示剔除函数。

此示例基于使用 2、3、5 和 7 的小素数的 210 号圆周轮,因为它似乎达到了阵列大小和位平面数量的效率“最佳点”,但实验表明通过将下一个 11 的质数添加到轮子上,对于 2310 个数字的圆周,可以再获得大约 5% 的性能;没有这样做的原因是初始化时间大大增加,并且很难为“仅”十亿的较小范围计时,因为那时只有大约四个段可以达到该点并且粒度成为问题。

请注意,第一个筛出的数字是 23,它是轮素数和预剔除素数之后的第一个素数;通过使用它,我们避免了处理从“一”开始的数组的问题以及恢复已消除且必须通过某些算法重新添加的轮素数的问题。

EDIT_ADD:添加对各种 WHL LUT 目的的解释 - 其中大部分是为了将计算时间昂贵的除法运算的需求减少到每个大页面段的每个基本素数大约两个。如下所述,WHLSTRTS LUT 用于将起始地址转换为筛页段的每个残差模索引(在本例中为 48)所需的位索引,以进行非常简单的查找、乘法和加法操作如下面所描述的。 END_EDIT_ADD

基本上,对于每个页段剔除扫描,都有一个起始循环,它用轮索引和段内第一个剔除地址的模剩余索引填充起始地址数组,每个基素数小于最大值的平方根页段中表示的数字,然后使用此起始地址数组完全依次筛选每个模剩余位平面(其中 48 个),并为每个位平面扫描所有基本素数,并从每个基本素数的段起始地址计算适当的偏移量通过使用 WHLSTARTS LUT 中的乘数和偏移量。这是通过将基本素数的轮索引与查找乘数相乘并添加查找的偏移量来获得给定模余数位平面的轮索引起始位置来完成的。从此,每个位平面的剔除就像第三章奇数位平面的剔除一样。这对每个位平面执行 48 次,但是对于 16 KB 缓冲区(每个位平面),每个页面段的有效范围是 210 轮跨度的 131072 倍或每个页面段的 27,525,120 个数字,并且从更大的值线性乘以位平面的大小。

与第 3 章仅赔率的筛子相比,使用此筛子可将内存使用量作为片段有效范围的比率减少 48 比 105 或不到一半,但因为每个片段都有所有 48位平面,完整的 Sieve 缓冲区是 16 KB 乘以 48 位平面或 768 KB(四分之三兆字节)和更大尺寸的倍数。然而,使用这种大小的 Sieve Buffer 的效率最高只有 160 亿左右,我们下一章的下一个示例将针对巨大的范围调整缓冲区的大小,使其增长到最大范围的大约 100 兆字节。对于多线程语言(不是 JavaScript),这将是每个线程的要求。

额外的内存要求用于存储 32 位值的基本素数数组,表示基本素数的轮索引及其模余数索引(如上所述的模地址计算所必需的);对于十亿的范围,大约有 3600 个基本素数乘以四个字节,每个大约 14,000 个字节,附加的起始地址数组大小相同。这些数组随着要筛选的最大值的平方根增长,因此对于筛选到 10^16(一万万亿)或大约所需的小于一亿的基本素数,增长到大约 5,761,455 个值(每个乘以四个字节)每个 23 兆字节。

进一步的改进适用于使用“组合”筛的以下示例,其中筛缓冲器从较大的轮型中预填充,其中已消除了质数为 11、13、17 和 19 的因子;比这更大的消除范围是不切实际的,因为保存的预剔除模式从每个模位平面仅大约 64 KB 增长到大约 20 倍大或大约 1.5 兆字节,对于 48 个模残数平面或大约 60 个仅通过添加 23 的质数的额外因子就可以达到兆字节 - 同样,它在内存和初始化方面的成本相当高,而性能只有百分之几。请注意,此数组可以共享以供所有线程使用。

在实现时,WHLPTRN 阵列约为 64 千字节乘以 48 模位平面约为 3 兆字节,它不是那么大,它是一个固定大小,不会随着筛分范围的增加而变化;对于访问速度和初始化时间,这是一个相当可行的大小。

这些“最大轮因式分解”改进减少了用于筛选十亿范围的合数剔除操作循环的总数,从第 3 章仅赔率示例的约十亿操作减少到此“组合”筛的约 25 亿,与目标是尝试保持每次剔除操作的 CPU 时钟周期数相同,从而将速度提高四倍。

已编辑: 已调整以下代码段以添加基本的 HTML 单网页应用程序用户界面,以便可以轻松调整参数以进行实验。为方便使用,点击“运行代码片段”按钮后应使用右上角的“整页”链接,完成后可通过右上角的链接关闭整页。要在智能手机上运行(最好在 Chrome 中),请使用设置菜单(三点菜单)中的“桌面站点”复选框。

EDIT_CORRECTION:由于可以在指定的上限内轻松更改范围限制,虚拟基基索引基素数数组大小不再足以覆盖 LIMIT 的指定上限 362 的平方根的平方根(以前仅 229)因此已增加到两个轮距或 439。

FURTHER_EDIT_CORRECTION: 填充大于 16384 字节的 SieveBuffer 剩余位平面缓冲区时,fillSieveBuffer 函数出现轻微错误,已更正

SPEED_OMISSION_CORRECTION: 通过使用其他语言,我们意识到该版本比应有的速度慢 20% 左右,因为没有有效地使用“循环剥离”来计算应该应用的适当限制。已添加并应用“bplmt”来纠正此问题。在第一次运行代码时,应该按几次 Sieve 按钮,让 JavaScript 引擎对生成的代码进行热调以进行优化,从而提高速度,在四到五次迭代后就会达到。

EDIT_POLISHING: 似乎将筛子缓冲区大小设置为 CPU L1 缓存大小确实没有任何优势,因为 JavaScript 的速度不足以利用其速度,并且 CPU L2 缓存大小或更少更好。唯一的问题是“筛子增加的粒度,因为筛子缓冲区大小现在分别代表每个筛子缓冲区选择大小的大约 2 亿、4 亿和 8 亿个范围。这使得对于诸如十亿这样的琐碎范围的计时测量对于实时来说是不准确的,因为整个额外缓冲区的计算可能会出现巨大的溢出以覆盖该范围。

此外,由于筛分范围功能已大大增加,因此添加了进度指示和取消功能。

然而,虽然筛分范围能力已增加,但需要添加“桶筛”的额外改进以保持大约 1e12 左右的效率,因此不建议筛分超过这一点。

上述 JavaScript 示例实现如下:

"use strict";

const WHLPRMS = new Uint32Array([2,3,5,7,11,13,17,19]);
const FRSTSVPRM = 23;
const WHLODDCRC = 105 | 0;
const WHLHITS = 48 | 0;
const WHLODDGAPS = new Uint8Array([
  3, 1, 3, 2, 1, 2, 3, 3, 1, 3, 2, 1, 3, 2, 3, 4,
  2, 1, 2, 1, 2, 4, 3, 2, 3, 1, 2, 3, 1, 3, 3, 2,
  1, 2, 3, 1, 3, 2, 1, 2, 1, 5, 1, 5, 1, 2, 1, 2 ]);
const RESIDUES = new Uint32Array([
  23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71,
  73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 121, 127,
  131, 137, 139, 143, 149, 151, 157, 163, 167, 169, 173, 179,
  181, 187, 191, 193, 197, 199, 209, 211, 221, 223, 227, 229, 233 ]);
const WHLNDXS = new Uint8Array([
  0, 0, 0, 1, 2, 2, 2, 3, 3, 4, 5, 5, 6, 6, 6,
  7, 7, 7, 8, 9, 9, 9, 10, 10, 11, 12, 12, 12, 13, 13,
  14, 14, 14, 15, 15, 15, 15, 16, 16, 17, 18, 18, 19, 20, 20,
  21, 21, 21, 21, 22, 22, 22, 23, 23, 24, 24, 24, 25, 26, 26,
  27, 27, 27, 28, 29, 29, 29, 30, 30, 30, 31, 31, 32, 33, 33,
  34, 34, 34, 35, 36, 36, 36, 37, 37, 38, 39, 39, 40, 41, 41,
  41, 41, 41, 42, 43, 43, 43, 43, 43, 44, 45, 45, 46, 47, 47, 48 ]);
const WHLRNDUPS = new Uint8Array( // two rounds to avoid overflow, used in start address calcs...
  [ 0, 3, 3, 3, 4, 7, 7, 7, 9, 9, 10, 12, 12, 15, 15,
    15, 18, 18, 18, 19, 22, 22, 22, 24, 24, 25, 28, 28, 28, 30,
    30, 33, 33, 33, 37, 37, 37, 37, 39, 39, 40, 42, 42, 43, 45,
    45, 49, 49, 49, 49, 52, 52, 52, 54, 54, 57, 57, 57, 58, 60,
    60, 63, 63, 63, 64, 67, 67, 67, 70, 70, 70, 72, 72, 73, 75,
    75, 78, 78, 78, 79, 82, 82, 82, 84, 84, 85, 87, 87, 88, 93,
    93, 93, 93, 93, 94, 99, 99, 99, 99, 99, 100, 102, 102, 103, 105,
    105, 108, 108, 108, 109, 112, 112, 112, 114, 114, 115, 117, 117, 120, 120,
    120, 123, 123, 123, 124, 127, 127, 127, 129, 129, 130, 133, 133, 133, 135,
    135, 138, 138, 138, 142, 142, 142, 142, 144, 144, 145, 147, 147, 148, 150,
    150, 154, 154, 154, 154, 157, 157, 157, 159, 159, 162, 162, 162, 163, 165,
    165, 168, 168, 168, 169, 172, 172, 172, 175, 175, 175, 177, 177, 178, 180,
    180, 183, 183, 183, 184, 187, 187, 187, 189, 189, 190, 192, 192, 193, 198,
    198, 198, 198, 198, 199, 204, 204, 204, 204, 204, 205, 207, 207, 208, 210, 210 ]);
const WHLSTARTS = function () {
  let arr = new Array(WHLHITS);
  for (let i = 0; i < WHLHITS; ++i) arr[i] = new Uint16Array(WHLHITS * WHLHITS);
  for (let pi = 0; pi < WHLHITS; ++pi) {
    let mltsarr = new Uint16Array(WHLHITS);
    let p = RESIDUES[pi]; let i = (p - FRSTSVPRM) >> 1;
    let s = ((i << 1) * (i + FRSTSVPRM) + (FRSTSVPRM * ((FRSTSVPRM - 1) >> 1))) | 0;
    // build array of relative mults and offsets to `s`...
    for (let ci = 0; ci < WHLHITS; ++ci) {
      let rmlt = (RESIDUES[((pi + ci) % WHLHITS) | 0] - RESIDUES[pi | 0]) >> 1;
      rmlt += rmlt < 0 ? WHLODDCRC : 0; let sn = s + p * rmlt;
      let snd = (sn / WHLODDCRC) | 0; let snm = (sn - snd * WHLODDCRC) | 0;
      mltsarr[WHLNDXS[snm]] = rmlt | 0; // new rmlts 0..209!
    }
    let ondx = (pi * WHLHITS) | 0
    for (let si = 0; si < WHLHITS; ++si) {
      let s0 = (RESIDUES[si] - FRSTSVPRM) >> 1; let sm0 = mltsarr[si];
      for (let ci = 0; ci < WHLHITS; ++ci) {
        let smr = mltsarr[ci];
        let rmlt = smr < sm0 ? smr + WHLODDCRC - sm0 : smr - sm0;
        let sn = s0 + p * rmlt; let rofs = (sn / WHLODDCRC) | 0;
        // we take the multiplier times 2 so it multiplies by the odd wheel index...
        arr[ci][ondx + si] = ((rmlt << 9) | (rofs | 0)) >>> 0;
      }
    }
  }
  return arr;
}();
const PTRNLEN = (11 * 13 * 17 * 19) | 0;
const PTRNNDXDPRMS = new Int32Array([ // the wheel index plus the modulo index
  (-1 << 6) + 44, (-1 << 6) + 45, (-1 << 6) + 46, (-1 << 6) + 47 ]);

function makeSieveBuffer(szbits) { // round up to 32 bit boundary!
  let arr = new Array(WHLHITS); let sz = ((szbits + 31) >> 5) << 2;
  for (let ri = 0; ri < WHLHITS; ++ri) arr[ri] = new Uint8Array(sz);
  return arr;
}

function cullSieveBuffer(lwi, bps, prmstrts, sb) {
  let len = sb[0].length; let szbits = len << 3; let bplmt = len >> 1;
  let lowndx = lwi * WHLODDCRC; let nxti = (lwi + szbits) * WHLODDCRC;
  // set up prmstrts for use by each modulo residue bit plane...
  for (let pi = 0, bpslmt = bps.length; pi < bpslmt; ++pi) {
    let ndxdprm = bps[pi] | 0;
    let prmndx = ndxdprm & 0x3F; let pd = ndxdprm >> 6;
    let rsd = RESIDUES[prmndx] | 0; let bp = (pd * (WHLODDCRC << 1) + rsd) | 0;
    let i = (bp - FRSTSVPRM) / 2;
    let s = (i + i) * (i + FRSTSVPRM) + (FRSTSVPRM * ((FRSTSVPRM - 1) / 2));
    if (s >= nxti) { prmstrts[pi] = 0xFFFFFFFF >>> 0; break; } // enough base primes!
    if (s >= lowndx) s = (s - lowndx) | 0;
    else {
      let wp = (rsd - FRSTSVPRM) >>> 1; let r = ((lowndx - s) % (WHLODDCRC * bp)) >>> 0;
      s = r == 0
        ? 0 | 0
        : (bp * (WHLRNDUPS[(wp + ((r + bp - 1) / bp) | 0) | 0] - wp) - r) | 0;
    }
    let sd = (s / WHLODDCRC) | 0; let sn = WHLNDXS[(s - sd * WHLODDCRC) | 0];
    prmstrts[pi | 0] = ((sn << 26) | sd) >>> 0;
  }
//  if (szbits == 131072) return;
  for (let ri = 0; ri < WHLHITS; ++ri) {
    let pln = sb[ri]; let plnstrts = WHLSTARTS[ri];
    for (let pi = 0, bpslmt = bps.length; pi < bpslmt; ++pi) {
      let prmstrt = prmstrts[pi | 0] >>> 0; if (prmstrt == 0xFFFFFFFF) break;
      let ndxdprm = bps[pi | 0] | 0;
      let prmndx = ndxdprm & 0x3F; let pd = ndxdprm >> 6;
      let bp = (((pd * (WHLODDCRC << 1)) | 0) + RESIDUES[prmndx]) | 0;
      let sd = prmstrt & 0x3FFFFFF; let sn = prmstrt >>> 26;
      let adji = (prmndx * WHLHITS + sn) | 0; let adj = plnstrts[adji];
      sd += ((((adj >> 8) * pd) | 0) + (adj & 0xFF)) | 0;
      if (bp < bplmt) {
        for (let slmt = Math.min(szbits, sd + (bp << 3)) | 0; sd < slmt; sd += bp) {
          let msk = (1 << (sd & 7)) >>> 0;
//          for (let c = sd >> 3, clmt = len == 16384 ? 0 : len; c < clmt; c += bp) pln[c] |= msk;
          for (let c = sd >> 3; c < len; c += bp) pln[c] |= msk;
        }
      }
//      else for (let sdlmt = szbits == 131072 ? 0 : szbits; sd < sdlmt; sd += bp) pln[sd >> 3] |= (1 << (sd & 7)) >>> 0;
      else for (; sd < szbits; sd += bp) pln[sd >> 3] |= (1 << (sd & 7)) >>> 0;
    }
  }
}

const WHLPTRN = function () {
  let sb = makeSieveBuffer((PTRNLEN + 16384) << 3); // avoid overflow when filling!
  cullSieveBuffer(0, PTRNNDXDPRMS, new Uint32Array(PTRNNDXDPRMS.length), sb);
  return sb;
}();

const CLUT = function () {
  let 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 countSieveBuffer(bitlmt, sb) {
  let lstwi = (bitlmt / WHLODDCRC) | 0;
  let lstri = WHLNDXS[(bitlmt - lstwi * WHLODDCRC) | 0];
  let lst = lstwi >> 5; let lstm = lstwi & 31;
  let count = (lst * 32 + 32) * WHLHITS;
  for (let ri = 0; ri < WHLHITS; ++ri) {
    let pln = new Uint32Array(sb[ri].buffer);
    for (let i = 0; i < lst; ++i) {
      let v = pln[i]; count -= CLUT[v & 0xFFFF]; count -= CLUT[v >>> 16];
    }
    let msk = 0xFFFFFFFF << lstm; if (ri <= lstri) msk <<= 1;
    let v = pln[lst] | msk; count -= CLUT[v & 0xFFFF]; count -= CLUT[v >>> 16];
  }
  return count;
}

function fillSieveBuffer(lwi, sb) {
  let len = sb[0].length; let cpysz = len > 16384 ? 16384 : len;
  let mod0 = lwi / 8;
  for (let ri = 0; ri < WHLHITS; ++ri) {
    for (let i = 0; i < len; i += 16384) {
      let mod = ((mod0 + i) % PTRNLEN) | 0;
      sb[ri].set(WHLPTRN[ri].subarray(mod, (mod + cpysz) | 0), i);
    }
  }
}

// a mutable cancelled flag...
let cancelled = false;

function doit() {
  const LIMIT =  Math.floor(parseFloat(document.getElementById('limit').value));
  if (!Number.isInteger(LIMIT) || (LIMIT < 0) || (LIMIT > 1e15)) {
    document.getElementById('output').innerText = "Top limit must be an integer between 0 and 9e15!";
    return;
  }
  const SIEVEBUFFERSZ = parseInt(document.getElementById('L1').value, 10);
  let startx = +Date.now();
  let count = 0;
  for (let i = 0; i < WHLPRMS.length; ++i) {
    if (WHLPRMS[i] > LIMIT) break;
    ++count;
  }
  if (LIMIT >= FRSTSVPRM) {
    const cmpsts = makeSieveBuffer(SIEVEBUFFERSZ);
    const bparr = function () {
      let szbits = (((((((Math.sqrt(LIMIT) | 0) - 23) >> 1) + WHLODDCRC - 1) / WHLODDCRC)
                                                                         + 31) >> 5) << 5;
      let cmpsts = makeSieveBuffer(szbits); fillSieveBuffer(0, cmpsts);
      let ndxdrsds = new Int32Array(2 * WHLHITS);
      for (let i = 0; i < ndxdrsds.length; ++i)
        ndxdrsds[i] = ((i < WHLHITS ? 0 : 64) + (i % WHLHITS)) >>> 0;
      cullSieveBuffer(0, ndxdrsds, new Uint32Array(ndxdrsds.length), cmpsts);
      let len = countSieveBuffer(szbits * WHLODDCRC - 1, cmpsts);
      let ndxdprms = new Uint32Array(len); let j = 0;
      for (let i = 0; i < szbits; ++i)
        for (let ri = 0; ri < WHLHITS; ++ri)
          if ((cmpsts[ri][i >> 3] & (1 << (i & 7))) == 0) {
            ndxdprms[j++] = ((i << 6) + ri) >>> 0;
          }
      return ndxdprms;
    }();
    let lwilmt = (LIMIT - FRSTSVPRM) / 2 / WHLODDCRC;
    let strts = new Uint32Array(bparr.length);
    let lwi = 0;
    const pgfnc = function () {
      if (cancelled) {
        document.getElementById('output').innerText = "Cancelled!!!";
        document.getElementById('go').value = "Start Sieve...";
        document.getElementById('go').disabled = false;
        cancelled = false;
        return;
      }
      const smlllmt = lwi + 4194304;
      const lmt = (smlllmt < lwilmt) ? smlllmt : lwilmt;
      for (; lwi <= lmt; lwi += SIEVEBUFFERSZ) {
        const nxti = lwi + SIEVEBUFFERSZ;
        fillSieveBuffer(lwi, cmpsts);
        cullSieveBuffer(lwi, bparr, strts, cmpsts);
        if (nxti <= lwilmt) count += countSieveBuffer(SIEVEBUFFERSZ * WHLODDCRC - 1, cmpsts);
        else count += countSieveBuffer((LIMIT - FRSTSVPRM) / 2 - lwi * WHLODDCRC, cmpsts);
      }
      if (lwi <= lwilmt) {
        document.getElementById('output').innerText = "Sieved " + ((lwi / lwilmt * 100).toFixed(3)) + "%";
        setTimeout(pgfnc, 7);
      }
      else {
       const elpsdx = +Date.now() - startx;
       document.getElementById('go').onclick = strtclick;
       document.getElementById('output').innerText = "Found " + count 
+ " primes up to " + LIMIT + " in " + elpsdx + " milliseconds.";
       document.getElementById('go').value = "Start Sieve...";
       document.getElementById('go').disabled = false;
      }
    };
    pgfnc();
  }
}

const cancelclick = function () {
  cancelled = true;
  document.getElementById('go').disabled = true;
  document.getElementById('go').value = "Cancelled!!!";
  document.getElementById('go').onclick = strtclick;
}

const strtclick = function () {
  document.getElementById('output').innerText = "Sieved 0%";
  document.getElementById('go').onclick = cancelclick;
  document.getElementById('go').value = "Running, click to cancel...";
  cancelled = false;
  setTimeout(doit, 7);
};

document.getElementById('go').onclick = strtclick;
html,
body {
  justify-content: center;
  align-items: center;
  text-align: center;
  font-weight: bold;
  font-size: 120%;
  margin-bottom: 10px;
}

.title {
  font-size: 200%;
}

.input {
  font-size: 100%;
  padding:5px 15px 5px 15px;
}

.output {
  padding:7px 15px 7px 15px;
}

.doit {
  font-weight: bold;
  font-size: 110%;
  border:3px solid black;
  background:#F0E5D1;
  padding:7px 15px 7px 15px;
}
<!doctype html>
<html>
<head>
  <meta http-equiv='Content-Type' content='text/html; charset=utf-8'>
  <meta name="viewport" content="width=device-width, initial-scale=1">
  <title>Page-Segmented Sieve of Eratosthenes...</title>
</head>
<body>
  <div class="title">
    <text>
      Page-Segmented Sieve of Eratosthenes
    </text>
  </div>
  <div>
    <text>
      Top sieve limit value:
    </text>
    <input class="input" type="textbox" id="limit" value="1e9" />
  </div>
  <div class="output">
    <text>The enforced limit is zero to 9e15, but values greater than about 1e12 can take a very long time!</text>
  </div>
  <div>
    <text>Sieve buffer size (CPU L2 cache?)</text>
    <select class="input" id="L1">
      <option value="1048576">128 Kilobytes</option>
      <option value="2097152">256 Kilobytes</option>
      <option value="4194304">512 Kilobytes</option>
    </select>
  </div>
  <div class="output">
    <text>Refresh the page to reset to default values or stop processing!</text>
  </div>
  <div class="output">
    <text id="output">
      Waiting to start...
    </text>
  </div>
  <div>
    <input class="doit" type="button" id="go" value="Start Sieve..." />
  </div>
</body>
</html>

请注意,上面的代码仅适用于超过 10 亿 (1e9) 的非平凡筛选范围,因为相当大的筛选页面大小的粒度使得对于低于此大小的范围的时间不太准确。 该程序最适用于 1E11 及以上的筛分范围。

由于以下原因,上述代码可能没有预期的那么快:

  1. 使用具有固定掩码图案的特殊简化剔除循环的加速技术不再有效,因为几乎没有任何小的剔除基本素数;这将每个合数剔除操作的平均时钟周期数增加了约 20%,尽管这更适用于 JavaScript 等较慢的语言而不是更高效的机器代码编译语言,因为它们可以使用更多技术,如极端循环展开来进一步将每个剔除操作循环的 CPU 时钟周期数减少到每个剔除循环大约 1.25 个时钟周期。

  2. 尽管由于较少的模位平面(大约那个因子),计算得到的素数的开销减少了大约两倍,但这不是所需的四倍;这在使用 JavaScript 时会变得更糟,因为 JavaScript 无法利用 CPU POP_COUNT 机器指令,这将比此处使用的计数 LUT (CLUT) 技术快大约十倍。

  3. 虽然此处使用的 LUT 技术将起始地址计算开销从所需的更复杂的模计算减少了五倍左右,但这些计算的复杂性大约是所需的两倍半到三倍。第 3 章中的“odds-only”筛子,因此给我们一个比率度量缩减是不够的;我们需要一种技术来将时间再减少两倍左右,以便这些计算不会对降低的比率做出贡献。相信我,我试过了,但没能做得更好。也就是说,这些计算可能在比 JavaScript 更高效的计算机语言和/或比我的极低端 Atom CPU 处理器更好的 CPU 中更有效,

尽管如此,几乎是四倍的加速,代码行数只增加了大约 50%,这还不算太糟糕,是吗?此 JavaScript 代码在较新版本的节点/Google Chrome 浏览器(Chrome 版本 75 仍比Firefox 版本 68 快 25%)上运行时仅慢三到五倍(取决于 CPU,更接近高端 CPU 的性能)Kim Walisch 的“primesieve”用“C”编写并编译为 x86_64 本机代码

我添加了一个附录答案,表明当项目规模增加到超过几百行时,实际上并不需要编写 JavaScript 来生成 JavaScript,就像这里一样。将来,我希望像 Fable 这样的转译器可能会发出 WebAssembly 代码而不是 JavaScript,现在用 Fable 编写的好处是,人们不必对代码进行更改(或至少进行少量更改)即可新技术的优势,它支持更快的代码执行和(最终)JavaScript不支持的多线程。

即将到来的是第 4.5b 章,它将增加大约两倍半的代码,但将是一个 Prime 生成器,能够筛选到非常大的范围,部分受限于 JavaScript 只能有效地表示高达 64 位浮点尾数为 53 的数字比特或大约 9e15 以及人们想要等待的时间:在更好的 CPU 上筛选到 1e14 需要一天的时间,但这不是什么大问题 - 只需打开浏览器选项卡即可!

于 2019-07-19T08:04:35.633 回答
3

附录:使用其他 JavaScript 生成语言

直到我的第 4a 章答案,根据问题的请求和代码,此线程中仅使用了 JavaScript。然而,在我看来,当编写超过几百行代码时,编写 JavaScript 不再是正确的方法。这有几个原因,如下所示:

  1. JavaScript 代码是很久以前设计的,所以它的编程模型已经过时了。这在很大程度上体现在它对遗留代码的支持上,尽管它已经(很多)按照 ECMA2015 和更新版本进行了更新。然而,这些更新导致了太多的做事方式,并不是所有的方式都高效,并且让程序员对最佳方式感到困惑。
  2. JavaScript 在使用原型模型来支持面向对象编程 (OOP) 版本的方式上是 FUGLY!
  3. JavaScript 是动态类型的,这可能会导致意外的代码错误,并使代码的维护和重构变得更加困难。
  4. 手动编码速度(通常使用 asm.js 和更新的功能)并不总能产生最佳代码,除非人们对 JavaScript 引擎使用的这些优化有很多了解,这些优化将实际执行代码。

有几个主流选项可以使用另一种语言通过转译来生成 JavaScript,其中两个最常见(也是最好的)如下:

  1. Microsoft 的 Typescript 是一种(可选)静态类型的 OOP 语言,由于上述原因而变得非常流行。
  2. 我最喜欢的是 Fable,它是 Microsoft 的主要功能性基于 ML 的语言 F# 的一个分支,它被开发用于生成高效的 JavaScript 代码,并且是一种更好的静态类型检查语言,但它提供了类型推断和各种简洁的功能。

如果我可以再避免它,我就不做 OOP(函数式编程 - FP - 是我的方法,只要它有意义,就像这里一样),但这里是第 4a 章代码的寓言版本(在移动设备上,在浏览器中作为桌面站点查看);对于第 4a 章的代码,应该按几次 Sieve 按钮,让 JavaScript 引擎对生成的代码进行热调优化,从而提高速度,四五次迭代后就会达到。使用 Fable,通过使用基于 Elmish React 的库,可以完全避免命令式代码,即使是用户界面 (UI),为了不混淆示例代码,我在这里没有这样做;同样,为了速度,我继续使用筛子剔除缓冲区作为可变数组——即使是最终的 FP 语言,Haskell,尽管受到“ST”monad 的保护,但也允许这种突变(可以在Fable/F#,但它们性能不高,因为 Haskell 没有自定义优化)。

CORRECTIONS_TO_COMMENTS: 事实证明,Chrome V8 JavaScript 引擎似乎已经优化了不可变的剔除突变,并且第一次更改的速度差异是由于使用命令式 JavaScript for/while 循环而不是 Fable 使用循环模拟尾递归函数,稍微慢一些。最大的效果是直接调用 JavaScript 的更好的数组拷贝。使用不同的位长进行计数也是一个很小的改进。总体改进大约为 25%,但副本改进与其他两者的总和大致相同。

在根据上述链接打开的页面中,可以查看生成的 JavaScript 代码,并会看到生成的纯 asm.js 代码比手动编写的代码更好(更一致),但 Fable 代码需要通过强制它在以下三个位置发出 JavaScript 代码对性能有一点帮助:

  1. Fable/F# 代码默认是不可变的,为了忠实于功能精神,我使用不可变尾调用递归函数循环,但这些比使用命令式 for/while 循环稍慢。
  2. 用于数组复制/blit 的 Fable 库没有(至少)使用最快的 JavaScript 方法(设置子数组/切片),所以我强制它发出 JavaScript 来做到这一点。
  3. 在计算剔除数组中未设置的位时,Fable 不提供 TypedArray 视图的其他位宽(我不认为),所以我通过发出 JavaScript 来添加它,因为它使用 Uint32 视图比使用四个连续的 Uint8 读取更快并通过位移手动组装 Uint32。

您会发现生成的代码与第 4a 章中手写的 JavaScript 代码一样快!

于 2020-04-02T03:43:32.590 回答