【问题标题】:Implementing the Page Segmented Sieve of Eratosthenes in Javascript在 Javascript 中实现 Eratosthenes 的页面分段筛
【发布时间】:2017-01-11 17:35:54
【问题描述】:

我最近读到了一个更快的 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 开始时。
  • 您能否发布这些案例,包括预期结果和实际结果?
  • 其实不只是少数情况。它跳过了几个数字来找出所有情况下的素数。已更新。

标签: javascript primes sieve-of-eratosthenes sieve


【解决方案1】:

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

概要

以下是逐步加快实施的阶段性进展,导致您的意图,并通过评论解释每个步骤的原因和实施细节。它包括 JavaScript 中可运行的 sn-ps,但这些技术不仅限于 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.");

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

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

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

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

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

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

  1. 我们制作了一个“轮子”,上面有两个位置,一半“命中”偶数,另一半“命中”奇数;因此,当我们将这个“轮子”“滚动”在所有潜在的素数上时,这个“轮子”的圆周跨度有两个数字,但它只“窗口”了它“滚动”的数字的二分之一或一半。
  2. 然后,我们将我们“滚动”的所有数字分成两个连续值的位平面,在一个位平面中,我们将所有偶数从 4 开始,在另一个位平面中,我们将所有赔率从三点钟。
  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 缓存大多数操作使用较小的基本素数,较大的、稍慢一点的 L2 高速缓存用于中等大小的基本素数,并且仅使用主 RAM 访问进行少数操作,使用大型基本素数来处理非常大的范围。
  3. 它通过将筛选每个较大页面的工作分配给不同的 CPU 内核以实现相当粗粒度的并发性,从而开启了多线程的可能性,尽管这不适用于 JavaScript,除非通过使用 Web Workers(混乱)。

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

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

START_EDIT:

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

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

END_EDIT

ANOTHER_EDIT:

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

END_ANOTHER_EDIT

FINAL_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 个潜在素数而不是跨度的轮子,使用 3、5 和 7仅 2) 并且还在小筛选数组的初始化时进行预剔除,这样就不需要通过以下素数 11、13、17 和 19 进行剔除。这减少了使用页面分段筛分大约四到十亿的范围(如表格中所示/根据维基百科文章 - “组合轮”中的公式计算),并且可以编写为它的运行速度大约是由于速度的四倍减少操作,每个剔除操作的速度与上述代码大致相同。

有效地进行 210 跨度车轮分解的方法是遵循“odds-only”方法:上面的当前算法可以被认为是从两个平面中筛选出一个位压缩平面,如前一章所述其中另一个平面可以被消除,因为它只包含两个以上的偶数;对于 210 跨度,我们可以定义 48 个这种大小的位压缩数组,表示 11 及以上的可能素数,其中所有其他 162 个平面包含的数字是 2、3、5 或 7 的因数,因此不需要被考虑。每个位平面都可以通过以基素数为增量重复索引来剔除,就像奇数平面是通过结构自动处理的乘法一样完成的,就像仅针对奇数一样。以这种方式,筛选内存需求更少(与 1/2 的“仅赔率”相比,筛选 48/210)和仅赔率的效率一样高,其中一个 48 平面“页面”表示 16 千字节 = 每个平面大小 131072 位乘以 210 是每个筛选页段的 27,525,120 个数字的范围,因此只有近 40 个页段可以筛选到十亿(而不是上面的近四千),因此在每个页段的每个基本素数的起始地址计算,以进一步提高效率。

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

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

由于此答案接近 30,000 个字符的限制大小,因此在 my followup Chapter 4.5a 和(未来)第 4.5b 章对上述技术的实际实现进行了进一步讨论。

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

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

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

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

所以这是关于使用 JavaScript 和其他当前 VM 语言的最新技术的限制,除了它们可以轻松使用多线程之外,它们具有大致相同的限制,结合上述因素 意味着本机代码编译器的速度可能比 JavaScript 快大约 20 倍(包括多线程,如果使用具有大量内核的新 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 速度。

结论

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

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

【讨论】:

  • 抱歉,我的空间用完了,无法发布代码 sn-p,该代码 最大轮式分解 速度大约快四倍,并且能够在一些很大的范围内进行筛选超过160亿;我正在制作另一个答案并在此处链接到它...
  • 这是一个很棒的教程,就像答案一样。我认为它涉及到关于筛选素数的所有知识,但不仅如此,它还展示了对类型化数组和位打包的有效使用。这些都不容易找到信息。我会好奇地将我修改后的 SoS 协议转换为这里找到的知识。这个话题应该受到保护。等待看到最大的车轮分解版本。
  • @Redu,我已经添加了一个带有 a further expansion of the Maximally Wheel Factorized technique 的答案,并可能会添加至少一个能够将这些技术扩展到 huge 范围的答案.这并没有比之前的最佳示例快四倍,但也不错,大约快了三倍半。
  • @Redu,关于您与 Sundaram 筛 (SoS) 的进一步合作,正如私下讨论的那样,SoS 相对于 SoE 的好处主要是它(稍微)更容易实现,但是总是会慢至少 ln(range) 的比率,因为它使用所有奇数而不是仅奇素数
  • 我在这里尝试了所有实现,但它们都受到相同问题的困扰——未能超过 32 位,这意味着大约 1.03 亿个素数。我无法理解如何在这里修改您的解决方案以超越 32 位。如果您可以将其包含在您的答案中,这将使其更加有用,IMO(JavaScript 可以处理 53 位整数)。干杯!
【解决方案2】:

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

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

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

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

这个例子是基于一个 210 的圆周轮,它使用了 2、3、5 和 7 的小素数,因为它似乎达到了阵列大小和位平面数量的效率的“最佳点”,但是实验表明,通过将下一个 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 千字节乘以 48 位平面或 768 千字节(四分之三兆字节)和更大尺寸的倍数。然而,使用这种大小的 Sieve Buffer 的效率最高只有 160 亿左右,我们下一章的下一个示例将针对巨大的范围调整缓冲区的大小,使其增长到最大范围的大约 100 兆字节。对于多线程语言(不是 JavaScript),这将是每个线程的要求。

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

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

在实现时,WHLPTRN 数组大约为 64 KB 乘以 48 模位平面大约为 3 MB,这不是那么大,并且它是一个固定大小,不会随着筛分范围的增加而改变;就访问速度和初始化时间而言,这是一个相当可行的大小。

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

已编辑: 已调整以下 sn-p 以添加基本的 HTML 单网页应用程序用户界面,以便可以轻松调整参数以进行实验。为方便使用,点击“Run code sn-p”按钮后应使用右上角的“Full page”链接,完成后可通过右上角的链接关闭整页。要在智能手机上运行(最好在 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 章中的“仅赔率”筛子是必需的,因此仅给我们一个比率度量缩减是不够的;我们需要一种技术来将时间再减少两倍左右,以便这些计算不会对降低的比率做出贡献。相信我,我试过了,但没能做得更好。也就是说,这些计算可能在比 JavaScript 更高效的计算机语言和/或比我的极低端 Atom CPU 处理器更好的 CPU 中更有效,但是复合数剔除速度也可能更高效!

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

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

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

【讨论】:

  • 哇! :O 我的筛分记录是每 1e7 220 毫秒。您的代码需要 15 毫秒。太棒了!
  • @GordoneBGood 我没有仔细阅读。我试图通过阅读代码来理解,但我发现这对我来说很难。我需要花更多的时间来理解这个想法。我希望我能找到一天深入研究。
  • @vitaly-t,不客气。很遗憾听到您放弃目前的筛分水平。这不是关于聪明,而是关于坚持:我在这里发布的各种代码是近 50 年编写筛子经验的结果 - 相信我,当我开始时,我不知道算法会导致这。至于以贡献者的身份加入你的 repo,我必须拒绝,因为我经常收到这样的请求,并且有其他项目我必须优先考虑;但是,如果你做一些新的事情,只需在提交中标记我,我会尝试找时间审查和评论。
  • @vitaly-t,你似乎错过了我的意思——我使用组合数值算法引用的技术并不慢,而且确实是目前已知最快的找到素数的通用算法以相当复杂为代价的大限制。 Pierre Dusart 的工作及其后来的其他扩展发现了素数计数的越来越精细的近似值,直到一个极限,但从未声称解析到一个精确的正确数字,否则我们将使用这些方法来设置记录,例如素数计数10^28 以上,尚不可知。
  • @vitaly-t,如果你真的想要最快,你需要超越 Meissel-Lehmer 到至少 Lagarias-Miller-Odlyzko (LMO),这至少是一个更快的订单,更大的订单更多范围。由于您的目标是拥有一个可以计算 53 位数字范围内的素数的网页,因此您不希望花费数天,但最多需要几分钟的进度条。该技术还可以通过查找素数计数达到近似极限,然后仅使用常规筛分来找到该近似极限与所需索引之间的差值处的素数来找到第 n 个素数...
【解决方案3】:

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

截至my Chapter 4a answer,根据问题的请求和代码,此线程中仅使用了 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,我不会再做 OOP(函数式编程 - FP - 是我的最佳选择,就像这里一样),但是 here is a Fable version of the Chapter 4a code(在移动设备上,在浏览器中查看作为桌面站点使用);对于第 4a 章的代码,应该按几次 Sieve 按钮,让 JavaScript 引擎对生成的代码进行热调优化,从而提高速度,四五次迭代后就会达到。使用 Fable,通过使用基于 Elmish React 的库,可以完全避免命令式代码,即使是用户界面 (UI),为了不混淆示例代码,我在这里没有这样做;同样,为了速度,我继续使用筛子剔除缓冲区作为可变数组——即使是最终的 FP 语言,Haskell,尽管受到“ST”monad 的保护,但也允许这种突变(可以在Fable/F#,但它们性能不佳,因为 Haskell 没有自定义优化)。

CORRECTIONS_TO_COMMENTS: Tt 事实证明,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 代码一样快!

【讨论】:

    猜你喜欢
    • 1970-01-01
    • 2016-04-06
    • 1970-01-01
    • 2022-11-28
    • 2017-01-05
    • 2014-07-20
    • 2015-09-06
    • 2021-12-19
    • 2019-04-14
    相关资源
    最近更新 更多