在 Javascript 中实现 Eratosthenes 的页面分段筛法
Implementing the Page Segmented Sieve of Eratosthenes in Javascript
我最近读到有关埃拉托色尼分段筛法对于非常大的数字的更快实现。
以下是相同的实现:
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)
应该 return [5,7]
但它是 returning [5,7,9]
尽管您从阅读中了解到 Eratosthenes 的分页筛法 是一种在大范围内快速查找素数的方法,但您的问题代码(即使已更正)没有实现页面分段 SoE,在大范围内测试代码,也没有像 SoE 实现那样特别快。以下讨论将展示如何在大范围内使用真正的页面分段 SoE。
剧情简介
以下是逐步实现您的意图的越来越快的实施过程,并附有解释每一步的原因和实施细节的评论。它包括 JavaScript 中的 运行 可用片段,但这些技术并不仅限于 JavaScript 并且其他语言不会限制某些进一步的改进,例如 multi-threading结果页面(Web Workers 除外,它们很难控制处理顺序),在极端循环展开中进行了一些进一步的优化,最后与代码效率有限有关,因为必须及时( JIT) 由浏览器中的 JavaScript 引擎编译为本机代码;这些限制与直接编译为非常高效的本机代码的语言相比,例如 C/C++、Nim、Rust、Free Pascal、Haskell、Julia 等
第 1 章 - 设定基准
首先,让我们从在相当大范围内使用的当前代码算法的工作版本开始,并使用时序信息来建立基线;以下代码在剔除素数的平方处开始剔除每个素数,这避免了剔除给定素数和一些冗余起始剔除的问题,并且没有理由为大范围生成结果素数的输出数组,因为我们可以生成直接来自剔除数组的素数;此外,答案的确定是在时间之外,因为我们将开发更好的技术来找到一个“答案”,对于大范围通常是找到的素数的计数,素数的总和,素数的第一次出现差距等,none其中需要实际查看找到的素数:
"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.");
现在,对于筛选到十亿的上述代码,引用我的 运行 次没有什么意义,因为你的时间几乎肯定会更快,因为我使用的是非常低端的 tablet CPU 在 Intel x5-Z8350 中 运行 以 1.92 GHz 的速度运行(CPU 单个活动线程的时钟速度)。我将仅引用一次我的执行时间作为如何计算每次剔除的平均 CPU 时钟周期的示例:我将上述代码的执行时间约为 43350 毫秒,即 43.35 秒乘以每秒 19.2 亿个时钟除以约 25.14 亿次剔除操作筛选到十亿(可以从公式计算或从 table 无轮分解 on the SoE Wikipedia page)得到约 33.1 CPU 每次剔除的时钟周期 到十亿。从现在开始,我们每次剔除只使用 CPU 个时钟周期来比较性能。
进一步注意,这些性能分数非常依赖于所使用的浏览器 JavaScript 引擎,以上分数为 运行 Google Chrome(版本 72) ; Microsoft Edge(版本 44)对于我们正在转向的页面分段 SoE 版本来说大约慢七倍,而 Firefox 和 Safari 的性能可能接近 Chrome。
由于使用了 Uint8Array
TypedArray 和更多 asm.js,此性能可能比之前的答案代码更好,但是这种“一个巨大的数组筛子”的时间(使用了 1 GB 的内存此处)由于超出 CPU 缓存限制的主 RAM 内存的内存访问速度而受到瓶颈。这就是我们致力于页面分段筛选的原因,但首先让我们做一些事情来减少使用的内存量和所需的剔除周期数。
第 2 章 - 位打包和仅赔率轮分解
下面的代码进行了位打包,这在紧密的内部剔除循环中需要稍微复杂一些,但由于它只对每个合数使用一位,因此它减少了八分之一的内存使用;同样,由于 2 是唯一的偶数素数,它使用基本的轮因式分解来筛选概率,只是为了将内存使用进一步减少两倍,并将剔除操作的数量减少约 2.5 倍。
odds-only 的最小轮子分解如下:
- 我们做一个“轮子”,上面有两个位置,一半“命中”偶数,另一半命中奇数;因此这个“轮子”有两个数字的圆周跨度,因为我们将它“滚动”到所有潜在的素数上,但它只是“windows”二分之一或二分之一它“滚动”的数字。
- 然后我们将我们“滚动”的所有数字分成两个连续值的位平面,在一个位平面中我们放置所有从 4 开始的偶数,在另一个位平面中我们放置所有赔率从三开始。
- 现在我们丢弃偶数平面,因为 none 表示的数字可能是素数。
p * p
对于奇数基素数的起始索引始终是奇数,因为奇数乘以奇数始终是奇数。
- 当我们将索引增加到奇数位平面时,我们实际上是在增加两倍的值,因为添加一个奇数和一个奇数产生一个偶数,这将在我们丢弃的位平面中,所以我们再次添加质数以再次回到奇数位平面。
- 奇数位平面索引位置自动说明了这一点,因为我们已经丢弃了之前在每个奇数位位置索引之间的所有偶数值。
- 这就是为什么我们可以通过在下面的代码中重复向索引添加素数来剔除。
"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.");
对于 10.257 亿次剔除操作,每个剔除操作的性能现在约为 34.75 CPU 个时钟周期(仅针对赔率)(来自维基百科),这意味着时间的减少几乎完全是由于剔除操作数量的减少,“bit-twiddling”位打包的额外复杂性仅花费与减少所节省的时间相同的额外时间内存使用量增加了 16 倍。因此,此版本使用内存的十六分之一,速度大约快 2.5 倍。
但是我们还没有完成,正如您的消息来源告诉您的那样,页面分段可以加快我们的速度。
第 3 章 - 页面分割,包括更快的素数计数
那么什么是页面分割应用于 SoE,它为我们做了什么?
页面分割将筛选工作从一次筛选的一个巨大数组划分为一系列依次筛选的较小页面。然后它需要稍微复杂一点,因为必须有一个单独的可用基素数流,这可以通过使用内筛递归筛选生成主筛使用的记忆基素数列表来获得。同样,输出结果的生成通常稍微复杂一些,因为它涉及对每个生成的筛选页面进行连续扫描和缩减。
页面分割有很多好处,如下所示:
- 它进一步降低了内存需求。使用之前的 odds-only 代码筛选到 10 亿需要大约 64 兆字节的 RAM,但无法使用该算法筛选超过 16 到 320 亿的范围(即使我们想等待很长时间) 由于使用了 JavaScript 可用的所有内存。通过按页面大小(比如)该数量的平方根进行筛选,我们可以筛选到该数量平方根的范围或超出我们有时间等待的任何范围。我们还必须将基数存储到所需范围的平方根,但这对于我们要考虑的任何范围来说都只有几十兆字节。
- 它提高了内存访问速度,这对每个剔除操作的 CPU 敲击周期数有直接影响。如果剔除操作主要发生在 CPU 缓存中,则内存访问从每次访问主 RAM 内存的数十 CPU 个时钟周期开始变化(取决于 CPU 以及 RAM 的质量和类型) CPU L1 缓存的大约一个时钟周期(较小,大约 16 KB)到 CPU L2 缓存的大约八个时钟周期(至少大约 128 KiloByte),我们可以计算出剔除算法,以便我们充分利用所有缓存,使用小型快速 L1 缓存进行大多数具有小基素值的操作,较大的 little-bit-slower L2 缓存用于中等大小的基素数,并且仅使用 main RAM 访问使用大基素数进行非常大范围的少数操作。
- 通过将筛选每个较大页面的工作分配给不同的 CPU 核心以实现相当粗粒度的并发,它开启了 multi-threading 的可能性,尽管这不适用于 JavaScript 除了通过使用 Web Workers(混乱)之外。
除了增加的复杂性之外,页面分段还有另一个问题需要解决:与“一个巨大的数组”筛选器不同,在“一个巨大的数组”筛选器中,起始索引很容易计算一次,然后用于整个数组,分段筛选器需要起始地址通过每页每个素数的模(除)运算计算(计算成本高),或者需要使用额外的内存来存储每页每个基本素数达到的索引,这样起始索引就不会必须重新计算,但最后一项技术排除了 multi-threading,因为这些数组不同步。 Ultimate 版本中将使用的最佳解决方案是结合使用这些技术,其中将几个页面段分组以形成一个相当大的线程工作单元,以便这些地址计算只占用总时间的一小部分,并且索引存储 table 用于每个线程跨这些较大工作单元的基本素数,以便每个较大工作单元只需完成一次复杂计算。因此我们得到了 multi-threading 的可能性和减少的开销。然而,下面的代码并没有减少这种开销,在筛选到十亿时 的开销大约为 10% 到 20%。
除了页面分割之外,以下代码通过使用一次计数 32 位的计数查找 Table (CLUT) 人口计数算法添加了对找到的素数的有效计数,以便开销在一小部分筛选时间中连续查找已找到素数的结果。如果不这样做,枚举找到的单个素数只是为了确定有多少素数,所花的时间至少与筛选给定范围所需的时间一样长。可以很容易地开发出类似的快速例程来执行诸如对找到的素数求和、查找素数间隙等操作。
START_EDIT:
以下代码添加了另一个 speed-up:对于较小的素数(此优化有效),代码通过识别剔除操作遵循八步模式来执行一种循环分离形式。这是因为一个字节的位数是偶数,我们通过奇素数剔除,每八次剔除将 return 到一个字节中的相同位位置;这意味着对于每个位位置,我们可以简化内部剔除循环以屏蔽一个常量位,从而大大简化内部循环并使剔除速度提高两倍左右,因为模式中的每个剔除循环都不需要执行"bit-twiddling" 位打包操作。此更改 节省了大约 35% 的执行时间筛选到十亿。可以通过将 64
更改为 0
来禁用它。它还为八循环的本机代码极端展开奠定了基础,因为这种模式可以在使用本机代码编译器时将剔除操作速度提高大约两倍。
通过对掩码值使用查找 Table (LUT) 而不是左移操作来节省大约半个时间,进一步的小修改使循环对于较大素数(大于 8192)更快CPU 剔除十亿范围时平均每个剔除操作的时钟周期;随着范围从 10 亿增加,这种节省会略有增加,但在 JavaScript 中效果不佳,已被删除。
END_EDIT
ANOTHER_EDIT:
除了上述编辑之外,我们还删除了 LUT 位掩码,但现在通过从相同大小的零缓冲区进行快速字节复制将筛选缓冲区归零,并添加了一个计数 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 章 - 添加最大轮分解的进一步高级工作
现在我们考虑使用最大轮分解的更广泛的变化(不仅“odds-only”为 2,而且涵盖 210 个潜在素数跨度的轮为 3、5 和 7而不是仅 2) 的跨度,以及 pre-culling 对小筛分阵列的初始化,这样就没有必要通过以下素数 11、13、17 和 19 进行剔除。这减少了复合材料的数量使用页面分段筛时的数字剔除操作大约四到十亿范围(如维基百科文章 - “组合轮”中公式中的 tables/calculated 所示)并且可以这样写它 运行 由于每次剔除操作的操作减少,速度大约是 abo 的四倍e码.
有效地进行 210 跨度轮分解的方法是遵循“odds-only”方法:上面的当前算法可以被认为是从两个平面中筛选出一个 bit-packed 平面正如上一章所解释的那样,另一个平面可以被消除,因为它只包含大于 2 的偶数;对于 210 跨度,我们可以定义 48 个 bit-packed 这种大小的数组,代表 11 及以上的可能素数,其中所有其他 162 个平面包含的数字是二、三、五或七的因数,因此不需要考虑。每个位平面都可以通过以基本质数为增量的重复索引来剔除,就像奇数平面完成时一样,乘法由结构自动处理,就像只有奇数一样。以这种方式,它以更少的内存需求进行筛选(与“odds-only”的 1/2 相比,筛选效率为 48/210)并且效率与仅赔率一样高,其中一个 48 平面“页”表示 16 KB = 131072 位/平面大小乘以 210 是每个筛选页段的 27,525,120 个数字的范围,因此只有将近 40 个页段筛选到十亿(而不是上面的近四千),因此更少每个页段每个基数的起始地址计算开销,以进一步提高效率。
虽然上面描述的扩展代码是几百行并且长到 post,但它可以在我的低端 Intel 1.92 Gigahertz CPU 使用 Google V8 JavaScript 引擎,这比本机代码中的相同算法 运行 慢大约五倍。
虽然上面的代码在大约 160 亿的范围内非常有效,但其他改进可以帮助将效率保持到更大的范围,例如 1e14 或更大的数万亿。我们通过向上调整有效页面大小来实现这一点,使它们永远不会小于被筛选的整个范围的平方根,但是通过 16 KB 的小质数块、128 KB 的中质数块逐步筛选它们,并且只筛选根据我们的基线实现,一个巨大的数组用于用于最大基本素数的极少数剔除操作。通过这种方式,对于我们可能考虑的最大范围,我们每次剔除的时钟不会增加超过两倍的小因子。
由于此答案接近 30,000 个字符的限制大小,因此在 和(未来)第 4.5b 章的答案中继续进一步讨论最大轮分解,以了解上述技术的实际实现。
第 5 章 - JavaScript(和其他 VM 语言)不能做什么
对于 JavaScript 和其他虚拟机语言,最小剔除循环时间约为每个剔除循环十个 CPU 周期,并且不太可能改变太多。这大约 比大约三个 CPU 时钟周期慢 三到四倍,这可以通过使用相同算法直接编译为高效机器代码的语言轻松实现,例如 C/C++、Nim、Rust、Free Pascal、Haskell、Julia 等
此外,还有一些极端的循环展开技术,至少可以与这些语言中的一些一起使用,这些语言可以将平均剔除操作周期减少 大约两倍,这被拒绝了至 JavaScript.
Multi-threading 可以将执行时间减少大约 CPU 使用的有效核心数,但是 JavaScript 唯一的方法是通过使用 Web Workers 和同步很麻烦。在我的机器上,我有四个内核,但由于所有内核都处于活动状态,CPU 时钟速率降低到四分之三,速度只增加了三倍; 对于 JavaScript.
来说,这个三倍的系数并不容易获得
所以这是关于 state-of-the-art 使用 JavaScript 和其他当前的 VM 语言除了可以轻松使用 multi-threading 之外具有大致相同的限制,结合上述因素意味着本机代码编译器可以比 JavaScript 快大约二十倍(包括 multi-threading,新的 CPU 和 巨大个核心数)。
不过,我相信未来三到五年的 Web 编程将是 Web Assembly,它有可能克服所有这些限制。它现在非常接近支持 multi-threading,虽然目前在 Chrome 上这个算法只比 JavaScript 快大约 30%,但在某些方面它只比本机代码慢一点使用某些 Web Assembly 编译器从某些语言编译时的当前浏览器。 WebAssembly 的高效编译器和本机代码的高效浏览器编译仍处于开发初期,但由于 WebAssembly 比大多数 VM 更接近本机代码, 它可以很容易地进行改进,以生成与其他 notive-code 编译语言的代码一样快或几乎一样快的本机代码。
然而,除了将 JavaScript 库和框架编译成 Web Assembly 之外,我认为 Web 的未来不会 JavaScript 到 Web Assembly 编译器,而是从一些其他语言。对于 Web 编程的未来,我最喜欢的选择是 F#,也许 Fable 实现转换为生成 Web Assembly,而不是 JavaScript (asm.js) 或 Nim。甚至有可能生成 Web Assembly,它支持并显示极端循环展开的优势,非常接近已知最快的页面分段 SoE 速度。
结论
我们在 JavaScript 中构建了埃拉托色尼的页分段筛法,它适用于 table 以筛分数十亿的大范围,并且有进一步扩展这项工作的方法。生成的代码的剔除操作减少了大约十倍(当完全轮子分解时)并且剔除操作快了大约三倍,这意味着每个给定(大)范围的代码快了大约 30 倍,而减少的内存使用意味着可以筛选到大约 9e15 的 53 位尾数的范围在一年左右(只需打开 Chrome 浏览器选项卡并备份电源)。
虽然还有一些其他可能的小调整,但这是关于使用 JavaScript 筛选素数的最新技术:虽然由于给定原因它不如本机代码快,但它很快足以在一两天内找到从 1e14 到 1e14 的素数数量(即使是在中等范围的 smartphone 上),只需在所需的计算时间内打开浏览器选项卡即可;这是非常惊人的,因为直到 1985 年才知道这个范围内的素数数量,然后通过使用数值分析技术而不是使用筛子,因为当时的计算机使用最快的编码技术还不够快在合理和经济的时间内。虽然我们可以在短短几个小时内使用这些算法在现代台式计算机上使用最好的本机代码编译器来完成此操作,但现在我们可以在 acceptable 时间内使用 JavaScript 在智能计算机上完成此操作phone!
这扩展了 以添加承诺但每个答案 30,000 个字符限制中没有空间的内容:
第 4.5a 章 - 全轮分解的基础知识
上一个答案中第 3 章的 non-Maximally-Wheel-Factorized Page-Segmented 埃拉托色尼筛法版本被编写为素数生成器,输出递归反馈作为基本素数馈送的输入;虽然这非常优雅且可扩展,但在接下来的工作中,我退后一步,使用更命令式的代码风格,以便读者可以更容易地理解最大轮分解的核心概念。在未来的第 4.5b 章中,我会将以下示例中开发的概念重新组合到质数生成器样式中,并添加一些额外的改进,这些改进不会使其在几十亿的较小范围内更快,但将使该概念在没有大部分损失速度可达数万亿至数百或数千万亿;随着范围变大,素数生成器格式在使程序适应方面更有用。
以下示例的主要额外改进是在用于有效寻址车轮模数残差的各种查找表 (LUT) 中,特殊起始地址 LUT 的生成非常简单地允许人们计算剔除起始地址对于每个模残差位平面,给定起始地址轮索引和整个结构中第一个剔除的第一个模残差位平面索引,没有额外的整数除法(慢),以及使用这些的 Sieve Buffer 复合数表示剔除函数。
此示例基于使用二、三、五和七小素数的 210 数字圆周轮,因为它似乎达到了数组大小和位数效率的“最佳点”平面,但实验表明,通过将下一个素数 11 添加到圆周为 2310 个数字的轮子中,可以将性能再提高 5%;没有这样做的原因是初始化时间大大增加,并且很难为“只有”十亿的较小范围计时,因为只有大约四个段才能达到这一点,粒度成为一个问题。
请注意,第一个筛选的数字是 23,这是轮素数和 pre-cull 个素数之后的第一个素数;通过使用它,我们避免了处理从“一”开始的数组的问题以及恢复被消除并且必须通过某些算法重新添加的轮素数的问题。
EDIT_ADD: 添加一个词来解释各种 WHL LUT 的用途 - 其中大部分与减少对 [=116= 的需求有关] 昂贵的除法操作到每个大页面段的每个基本质数只有大约两个。如下所述,WHLSTRTS LUT 用于将起始地址转换为筛选页段(在本例中为 48)的每个基本质数每个剩余模数索引所需的位索引,以按原样进行非常简单的查找、乘法和加法运算如下面所描述的。 END_EDIT_ADD
基本上,对于每个页面段剔除扫描,都有一个起始循环,该循环使用轮索引和段内第一个剔除地址的模残基索引填充起始地址数组,用于每个小于页段中表示的最大数的平方根,然后使用此起始地址数组依次完全筛选每个模余数位平面(其中 48 个),并为每个位平面扫描所有基本素数,并从该段计算出适当的偏移量通过使用 WHLSTARTS LUT 中的乘数和偏移量来计算每个基本质数的起始地址。这是通过将基素数的轮索引与查找的乘数相乘并加上查找的偏移量以获得给定模余数位平面的轮索引起始位置来完成的。从这里开始,每位平面的剔除就像第三章奇数位平面的剔除一样。这对每个位平面完成 48 次,但是 16 KB 缓冲区(每个位平面)的每个页面段的有效范围是 210 轮跨度的 131072 倍或每个页面段 27,525,120 个数字,并从更大的范围线性乘以位平面的大小。
与第 3 章 odds-only 筛法相比,使用此筛法将内存使用量减少为 48 倍于 105 倍或不到一半,但是因为每个段都有 48 位平面,所以完整的 Sieve Buffer 是 48 位平面的 16 KB 或 768 KB(四分之三兆字节)和更大尺寸的倍数。然而,使用这种大小的 Sieve Buffer 的效率最高只能达到 160 亿左右,我们下一章的下一个示例将调整缓冲区的大小以适应巨大的范围,使其增长到大约 100 兆字节 fr 最大的范围。对于 multi-threaded 种语言(不是 JavaScript),这是每个线程的要求。
额外的内存要求用于存储 32 位值的基本素数数组,这些值表示基本素数的轮索引及其模余数索引(如上所述,模地址计算是必需的);对于 10 亿的范围,大约有 3600 个基本素数乘以四个字节,每个大约 14,000 个字节,附加的起始地址数组大小相同。这些数组随着要筛选的最大值的平方根而增长,因此对于筛选到 10^16(万亿)或大约每个 23 兆字节。
进一步细化适用于以下示例,使用“组合”筛子,其中筛子缓冲区 pre-filled 来自更大的轮子模式,从中可以得到 11、13、17、11、13、17、十九人,被淘汰;比这更大的消除范围是不切实际的,因为保存的 pre-culled 模式从每个模数位平面仅约 64 KB 增长到 48 个模数残基数平面中每个平面的大约 20 倍或大约 1.5 兆字节或仅通过添加素数 23 的额外因子约 60 兆字节 - 同样,这在内存和初始化方面的成本相当大,而性能仅占百分之几。请注意,此数组可以共享以供所有线程使用。
实施时,WHLPTRN 阵列约为 64 KB 乘以 48 模位平面约为 3 MB,这不是那么大,并且是固定大小,不会随着筛分范围的增加而改变;就访问速度和初始化时间而言,这是一个非常可行的大小。
这些“最大轮分解”改进减少了用于筛选十亿范围的合数剔除操作循环的总数,从第 3 章 odds-only 示例的大约十亿操作减少到这个“组合”筛子,目的是尽量使每个剔除操作的 CPU 时钟周期数保持相同,从而使速度提高四倍。
已编辑: 已调整以下代码段以添加基本的 HTML 单网页应用程序用户界面,以便可以轻松调整参数以进行实验。为了最好的易用性,点击“运行 代码段”按钮后,应使用右上角的“整页”link,并且可以使用右上角的 link 关闭整页等结束了。要在智能手机上 运行(最好在 Chrome 中),请使用设置菜单(三点菜单)中的“桌面站点”复选框。
EDIT_CORRECTION: 范围限制可以在指定的上限内轻松更改,虚拟基基索引的基素数数组大小不再足够覆盖 LIMIT 的指定上限 362(以前只有 229)的平方根的平方根,因此已增加到两个轮距或 439。
FURTHER_EDIT_CORRECTION: fillSieveBuffer函数在填充大于16384字节的SieveBuffer剩余位平面缓冲区时出现轻微错误,已更正
SPEED_OMISSION_CORRECTION: 从与其他语言的合作中,我们意识到由于没有有效地使用“loop unpeeling”,该版本比应有的速度慢了大约 20% “至于在应该应用的地方不计算适当的限制。已添加并应用“bplmt”来更正此问题。第一次 运行 编译代码时,应多次按下 Sieve 按钮,让 JavaScript 引擎对生成的代码进行热调以进行优化,从而提高速度,它将在四五次迭代后达到.
EDIT_POLISHING: 看来把sieve buffer size设置成CPU L1 cache size确实没什么优势如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>
请注意,上面的代码仅适用于 non-trivial 超过十亿 (1e9) 的筛分范围,因为相当大的筛分页面尺寸的粒度使得时间不那么准确范围低于此尺寸。 该程序最适用于 1E11 及以上的筛分范围。
由于以下原因,上述代码可能没有预期的那么快:
使用具有固定掩码模式的特殊简化剔除循环的加速技术不再有效,因为几乎不再有任何小的剔除基素数;这将每个合数剔除操作的平均时钟周期数增加了大约 20%,尽管这更适用于 JavaScript 等较慢的语言,而不是更高效的机器代码编译语言,因为它们可以使用更进一步的技术,例如 extreme循环展开以进一步将每个剔除操作循环的 CPU 时钟周期数减少到每个剔除循环大约 1.25 个时钟周期。
尽管由于较少的模位平面(大约那个因子),计算结果素数的开销减少了大约两倍,但这不是所需的四倍;这在使用 JavaScript 时变得更糟,它无法利用 CPU POP_COUNT 机器指令,这比使用的计数 LUT (CLUT) 技术快大约十倍这里。
虽然此处使用的 LUT 技术将起始地址计算开销减少了五分之一左右,而对于所需的更复杂的模计算,这些计算大约是两倍半比第 3 章中的“odds-only”筛法要求的复杂三倍,因此不足以减少 ratio-metric;我们需要一种技术将时间减少两倍左右,以便这些计算不会对减少的比率做出贡献。相信我,我试过了,但没能做得更好。也就是说,这些计算在更高效的计算机语言中可能比 JavaScript and/or 更好 CPU 比我非常低端的 Atom CPU 处理器更有效,但是合数剔除速度也可能更高效!
不过,速度几乎提高了四倍,代码行数只增加了大约 50%,这还算不错,不是吗?当 运行 在较新版本的 node/Google Chrome 浏览器上(Chrome 版本 75 仍然比 Firefox 版本 68 快 25%) 比 Kim Walisch 写的“primesieve” “C”并编译为 x86_64 本机代码!
我添加了一个 ,表明当项目规模增加到超过 2 个时,实际上并不需要编写 JavaScript 来生成 JavaScript一百行,就像这里一样。在未来,我希望像 Fable 这样的转译器可以发出 WebAssembly 代码而不是 JavaScript,现在用 Fable 编写的好处是不需要对代码进行更改(或至少很少更改)为了利用更新的技术,支持更快的代码执行和(最终)JavaScript 不支持的 multi-threading。
即将到来的是第 4.5b 章,这将是大约两倍半的代码,但将是一个能够筛选到极大范围的素数生成器,部分受限于 JavaScript 只能有效地表示高达 53 位的 64 位浮点位尾数或大约 9e15 的数字,以及人们想要等待的时间:在更好的 CPU 上筛选到 1e14 将占用一天的顺序,但这不是什么大问题 - 只需打开浏览器选项卡即可!
附录:使用其他JavaScript生成语言
最多 ,根据问题的要求和代码,此线程中仅使用了 JavaScript。但是,在我看来,编写超过几百行的代码时,编写 JavaScript 不再是正确的方法。这有几个原因,如下:
- JavaScript代码是很久以前设计的,所以它的编程模型已经过时了。这在很大程度上体现在它对遗留代码的支持,尽管它已经根据 ECMA2015 和更新版本进行了(大量)更新。然而,更新导致了太多的做事方式,并不是所有的方式都有效,让程序员对最好的方式感到困惑。
- JavaScript 使用原型模型支持面向对象编程 (OOP) 版本的方式非常糟糕!
- JavaScript 是动态类型的,这可能会导致意外的代码错误,并使维护和重构代码变得更加困难。
- 为了速度而手动编码(通常使用 asm.js 和更新的功能)并不总能产生最佳代码,除非人们非常了解 JavaScript 引擎使用的这些优化,这些优化实际上会执行代码。
有几个主流选项可以使用另一种语言通过转译生成 JavaScript,其中两个最常见(也是最好)如下:
- Microsoft 的 Typescript 是一种(可选)静态类型的 OOP 语言,由于上述原因而变得非常流行。
- 我最喜欢的是 Fable,它是微软主要基于 ML 的函数式语言 F# 的一个分支,开发用于生成高效的 JavaScript 代码,并且是一种更好的静态类型检查语言,但它提供类型推断以及各种简洁的功能。
如果我能再避免它,我不会做 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# 中实现 monod 的想法,但它们的性能不佳,因为没有 Haskell 具有的自定义优化)。
CORRECTIONS_TO_COMMENTS: Tt 结果表明 Chrome V8 JavaScript 引擎似乎已经优化了不可变剔除突变和差异第一次更改的速度是由于使用了命令式 JavaScript for/while 循环,而不是 Fable 使用循环模拟尾递归函数,后者速度稍慢。最好的数组复制直接调用 JavaScript 的最大效果。使用不同的位长进行计数也是一个相当小的改进。总改进约为 25%,但复制改进与其他两个组合的改进大致相同。
在按上述link打开的页面中,可以查看生成的JavaScript代码,会看到生成的纯asm.js代码更好(更一致) ) 比手写的要多,但是 Fable 代码需要一点性能帮助,迫使它在三个地方发出 JavaScript 代码,如下所示:
- Fable/F# 代码默认是不可变的,为了忠实于函数精神,我使用不可变的尾调用递归函数循环,但这些比使用命令式 for/while 循环稍微慢一些。
- 数组 copy/blit 的寓言库没有(至少)使用最快的 JavaScript 方法(设置 subarray/slice),所以我强制它发出 JavaScript 这样做。
- 在计算剔除数组中未设置的位时,Fable 不提供 TypedArray 视图的其他位宽(我不认为),所以我通过发射 JavaScript 添加它作为它对 Uint32 的使用查看比使用四个连续的 Uint8 读取和通过位移手动组装 Uint32 更快。
您会发现生成的代码与第 4a 章中手写的 JavaScript 代码一样快!
我最近读到有关埃拉托色尼分段筛法对于非常大的数字的更快实现。
以下是相同的实现:
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)
应该 return [5,7]
但它是 returning [5,7,9]
尽管您从阅读中了解到 Eratosthenes 的分页筛法 是一种在大范围内快速查找素数的方法,但您的问题代码(即使已更正)没有实现页面分段 SoE,在大范围内测试代码,也没有像 SoE 实现那样特别快。以下讨论将展示如何在大范围内使用真正的页面分段 SoE。
剧情简介
以下是逐步实现您的意图的越来越快的实施过程,并附有解释每一步的原因和实施细节的评论。它包括 JavaScript 中的 运行 可用片段,但这些技术并不仅限于 JavaScript 并且其他语言不会限制某些进一步的改进,例如 multi-threading结果页面(Web Workers 除外,它们很难控制处理顺序),在极端循环展开中进行了一些进一步的优化,最后与代码效率有限有关,因为必须及时( JIT) 由浏览器中的 JavaScript 引擎编译为本机代码;这些限制与直接编译为非常高效的本机代码的语言相比,例如 C/C++、Nim、Rust、Free Pascal、Haskell、Julia 等
第 1 章 - 设定基准
首先,让我们从在相当大范围内使用的当前代码算法的工作版本开始,并使用时序信息来建立基线;以下代码在剔除素数的平方处开始剔除每个素数,这避免了剔除给定素数和一些冗余起始剔除的问题,并且没有理由为大范围生成结果素数的输出数组,因为我们可以生成直接来自剔除数组的素数;此外,答案的确定是在时间之外,因为我们将开发更好的技术来找到一个“答案”,对于大范围通常是找到的素数的计数,素数的总和,素数的第一次出现差距等,none其中需要实际查看找到的素数:
"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.");
现在,对于筛选到十亿的上述代码,引用我的 运行 次没有什么意义,因为你的时间几乎肯定会更快,因为我使用的是非常低端的 tablet CPU 在 Intel x5-Z8350 中 运行 以 1.92 GHz 的速度运行(CPU 单个活动线程的时钟速度)。我将仅引用一次我的执行时间作为如何计算每次剔除的平均 CPU 时钟周期的示例:我将上述代码的执行时间约为 43350 毫秒,即 43.35 秒乘以每秒 19.2 亿个时钟除以约 25.14 亿次剔除操作筛选到十亿(可以从公式计算或从 table 无轮分解 on the SoE Wikipedia page)得到约 33.1 CPU 每次剔除的时钟周期 到十亿。从现在开始,我们每次剔除只使用 CPU 个时钟周期来比较性能。
进一步注意,这些性能分数非常依赖于所使用的浏览器 JavaScript 引擎,以上分数为 运行 Google Chrome(版本 72) ; Microsoft Edge(版本 44)对于我们正在转向的页面分段 SoE 版本来说大约慢七倍,而 Firefox 和 Safari 的性能可能接近 Chrome。
由于使用了 Uint8Array
TypedArray 和更多 asm.js,此性能可能比之前的答案代码更好,但是这种“一个巨大的数组筛子”的时间(使用了 1 GB 的内存此处)由于超出 CPU 缓存限制的主 RAM 内存的内存访问速度而受到瓶颈。这就是我们致力于页面分段筛选的原因,但首先让我们做一些事情来减少使用的内存量和所需的剔除周期数。
第 2 章 - 位打包和仅赔率轮分解
下面的代码进行了位打包,这在紧密的内部剔除循环中需要稍微复杂一些,但由于它只对每个合数使用一位,因此它减少了八分之一的内存使用;同样,由于 2 是唯一的偶数素数,它使用基本的轮因式分解来筛选概率,只是为了将内存使用进一步减少两倍,并将剔除操作的数量减少约 2.5 倍。
odds-only 的最小轮子分解如下:
- 我们做一个“轮子”,上面有两个位置,一半“命中”偶数,另一半命中奇数;因此这个“轮子”有两个数字的圆周跨度,因为我们将它“滚动”到所有潜在的素数上,但它只是“windows”二分之一或二分之一它“滚动”的数字。
- 然后我们将我们“滚动”的所有数字分成两个连续值的位平面,在一个位平面中我们放置所有从 4 开始的偶数,在另一个位平面中我们放置所有赔率从三开始。
- 现在我们丢弃偶数平面,因为 none 表示的数字可能是素数。
p * p
对于奇数基素数的起始索引始终是奇数,因为奇数乘以奇数始终是奇数。- 当我们将索引增加到奇数位平面时,我们实际上是在增加两倍的值,因为添加一个奇数和一个奇数产生一个偶数,这将在我们丢弃的位平面中,所以我们再次添加质数以再次回到奇数位平面。
- 奇数位平面索引位置自动说明了这一点,因为我们已经丢弃了之前在每个奇数位位置索引之间的所有偶数值。
- 这就是为什么我们可以通过在下面的代码中重复向索引添加素数来剔除。
"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.");
对于 10.257 亿次剔除操作,每个剔除操作的性能现在约为 34.75 CPU 个时钟周期(仅针对赔率)(来自维基百科),这意味着时间的减少几乎完全是由于剔除操作数量的减少,“bit-twiddling”位打包的额外复杂性仅花费与减少所节省的时间相同的额外时间内存使用量增加了 16 倍。因此,此版本使用内存的十六分之一,速度大约快 2.5 倍。
但是我们还没有完成,正如您的消息来源告诉您的那样,页面分段可以加快我们的速度。
第 3 章 - 页面分割,包括更快的素数计数
那么什么是页面分割应用于 SoE,它为我们做了什么?
页面分割将筛选工作从一次筛选的一个巨大数组划分为一系列依次筛选的较小页面。然后它需要稍微复杂一点,因为必须有一个单独的可用基素数流,这可以通过使用内筛递归筛选生成主筛使用的记忆基素数列表来获得。同样,输出结果的生成通常稍微复杂一些,因为它涉及对每个生成的筛选页面进行连续扫描和缩减。
页面分割有很多好处,如下所示:
- 它进一步降低了内存需求。使用之前的 odds-only 代码筛选到 10 亿需要大约 64 兆字节的 RAM,但无法使用该算法筛选超过 16 到 320 亿的范围(即使我们想等待很长时间) 由于使用了 JavaScript 可用的所有内存。通过按页面大小(比如)该数量的平方根进行筛选,我们可以筛选到该数量平方根的范围或超出我们有时间等待的任何范围。我们还必须将基数存储到所需范围的平方根,但这对于我们要考虑的任何范围来说都只有几十兆字节。
- 它提高了内存访问速度,这对每个剔除操作的 CPU 敲击周期数有直接影响。如果剔除操作主要发生在 CPU 缓存中,则内存访问从每次访问主 RAM 内存的数十 CPU 个时钟周期开始变化(取决于 CPU 以及 RAM 的质量和类型) CPU L1 缓存的大约一个时钟周期(较小,大约 16 KB)到 CPU L2 缓存的大约八个时钟周期(至少大约 128 KiloByte),我们可以计算出剔除算法,以便我们充分利用所有缓存,使用小型快速 L1 缓存进行大多数具有小基素值的操作,较大的 little-bit-slower L2 缓存用于中等大小的基素数,并且仅使用 main RAM 访问使用大基素数进行非常大范围的少数操作。
- 通过将筛选每个较大页面的工作分配给不同的 CPU 核心以实现相当粗粒度的并发,它开启了 multi-threading 的可能性,尽管这不适用于 JavaScript 除了通过使用 Web Workers(混乱)之外。
除了增加的复杂性之外,页面分段还有另一个问题需要解决:与“一个巨大的数组”筛选器不同,在“一个巨大的数组”筛选器中,起始索引很容易计算一次,然后用于整个数组,分段筛选器需要起始地址通过每页每个素数的模(除)运算计算(计算成本高),或者需要使用额外的内存来存储每页每个基本素数达到的索引,这样起始索引就不会必须重新计算,但最后一项技术排除了 multi-threading,因为这些数组不同步。 Ultimate 版本中将使用的最佳解决方案是结合使用这些技术,其中将几个页面段分组以形成一个相当大的线程工作单元,以便这些地址计算只占用总时间的一小部分,并且索引存储 table 用于每个线程跨这些较大工作单元的基本素数,以便每个较大工作单元只需完成一次复杂计算。因此我们得到了 multi-threading 的可能性和减少的开销。然而,下面的代码并没有减少这种开销,在筛选到十亿时 的开销大约为 10% 到 20%。
除了页面分割之外,以下代码通过使用一次计数 32 位的计数查找 Table (CLUT) 人口计数算法添加了对找到的素数的有效计数,以便开销在一小部分筛选时间中连续查找已找到素数的结果。如果不这样做,枚举找到的单个素数只是为了确定有多少素数,所花的时间至少与筛选给定范围所需的时间一样长。可以很容易地开发出类似的快速例程来执行诸如对找到的素数求和、查找素数间隙等操作。
START_EDIT:
以下代码添加了另一个 speed-up:对于较小的素数(此优化有效),代码通过识别剔除操作遵循八步模式来执行一种循环分离形式。这是因为一个字节的位数是偶数,我们通过奇素数剔除,每八次剔除将 return 到一个字节中的相同位位置;这意味着对于每个位位置,我们可以简化内部剔除循环以屏蔽一个常量位,从而大大简化内部循环并使剔除速度提高两倍左右,因为模式中的每个剔除循环都不需要执行"bit-twiddling" 位打包操作。此更改 节省了大约 35% 的执行时间筛选到十亿。可以通过将 64
更改为 0
来禁用它。它还为八循环的本机代码极端展开奠定了基础,因为这种模式可以在使用本机代码编译器时将剔除操作速度提高大约两倍。
通过对掩码值使用查找 Table (LUT) 而不是左移操作来节省大约半个时间,进一步的小修改使循环对于较大素数(大于 8192)更快CPU 剔除十亿范围时平均每个剔除操作的时钟周期;随着范围从 10 亿增加,这种节省会略有增加,但在 JavaScript 中效果不佳,已被删除。
END_EDIT
ANOTHER_EDIT:
除了上述编辑之外,我们还删除了 LUT 位掩码,但现在通过从相同大小的零缓冲区进行快速字节复制将筛选缓冲区归零,并添加了一个计数 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 章 - 添加最大轮分解的进一步高级工作
现在我们考虑使用最大轮分解的更广泛的变化(不仅“odds-only”为 2,而且涵盖 210 个潜在素数跨度的轮为 3、5 和 7而不是仅 2) 的跨度,以及 pre-culling 对小筛分阵列的初始化,这样就没有必要通过以下素数 11、13、17 和 19 进行剔除。这减少了复合材料的数量使用页面分段筛时的数字剔除操作大约四到十亿范围(如维基百科文章 - “组合轮”中公式中的 tables/calculated 所示)并且可以这样写它 运行 由于每次剔除操作的操作减少,速度大约是 abo 的四倍e码.
有效地进行 210 跨度轮分解的方法是遵循“odds-only”方法:上面的当前算法可以被认为是从两个平面中筛选出一个 bit-packed 平面正如上一章所解释的那样,另一个平面可以被消除,因为它只包含大于 2 的偶数;对于 210 跨度,我们可以定义 48 个 bit-packed 这种大小的数组,代表 11 及以上的可能素数,其中所有其他 162 个平面包含的数字是二、三、五或七的因数,因此不需要考虑。每个位平面都可以通过以基本质数为增量的重复索引来剔除,就像奇数平面完成时一样,乘法由结构自动处理,就像只有奇数一样。以这种方式,它以更少的内存需求进行筛选(与“odds-only”的 1/2 相比,筛选效率为 48/210)并且效率与仅赔率一样高,其中一个 48 平面“页”表示 16 KB = 131072 位/平面大小乘以 210 是每个筛选页段的 27,525,120 个数字的范围,因此只有将近 40 个页段筛选到十亿(而不是上面的近四千),因此更少每个页段每个基数的起始地址计算开销,以进一步提高效率。
虽然上面描述的扩展代码是几百行并且长到 post,但它可以在我的低端 Intel 1.92 Gigahertz CPU 使用 Google V8 JavaScript 引擎,这比本机代码中的相同算法 运行 慢大约五倍。
虽然上面的代码在大约 160 亿的范围内非常有效,但其他改进可以帮助将效率保持到更大的范围,例如 1e14 或更大的数万亿。我们通过向上调整有效页面大小来实现这一点,使它们永远不会小于被筛选的整个范围的平方根,但是通过 16 KB 的小质数块、128 KB 的中质数块逐步筛选它们,并且只筛选根据我们的基线实现,一个巨大的数组用于用于最大基本素数的极少数剔除操作。通过这种方式,对于我们可能考虑的最大范围,我们每次剔除的时钟不会增加超过两倍的小因子。
由于此答案接近 30,000 个字符的限制大小,因此在
第 5 章 - JavaScript(和其他 VM 语言)不能做什么
对于 JavaScript 和其他虚拟机语言,最小剔除循环时间约为每个剔除循环十个 CPU 周期,并且不太可能改变太多。这大约 比大约三个 CPU 时钟周期慢 三到四倍,这可以通过使用相同算法直接编译为高效机器代码的语言轻松实现,例如 C/C++、Nim、Rust、Free Pascal、Haskell、Julia 等
此外,还有一些极端的循环展开技术,至少可以与这些语言中的一些一起使用,这些语言可以将平均剔除操作周期减少 大约两倍,这被拒绝了至 JavaScript.
Multi-threading 可以将执行时间减少大约 CPU 使用的有效核心数,但是 JavaScript 唯一的方法是通过使用 Web Workers 和同步很麻烦。在我的机器上,我有四个内核,但由于所有内核都处于活动状态,CPU 时钟速率降低到四分之三,速度只增加了三倍; 对于 JavaScript.
来说,这个三倍的系数并不容易获得所以这是关于 state-of-the-art 使用 JavaScript 和其他当前的 VM 语言除了可以轻松使用 multi-threading 之外具有大致相同的限制,结合上述因素意味着本机代码编译器可以比 JavaScript 快大约二十倍(包括 multi-threading,新的 CPU 和 巨大个核心数)。
不过,我相信未来三到五年的 Web 编程将是 Web Assembly,它有可能克服所有这些限制。它现在非常接近支持 multi-threading,虽然目前在 Chrome 上这个算法只比 JavaScript 快大约 30%,但在某些方面它只比本机代码慢一点使用某些 Web Assembly 编译器从某些语言编译时的当前浏览器。 WebAssembly 的高效编译器和本机代码的高效浏览器编译仍处于开发初期,但由于 WebAssembly 比大多数 VM 更接近本机代码, 它可以很容易地进行改进,以生成与其他 notive-code 编译语言的代码一样快或几乎一样快的本机代码。
然而,除了将 JavaScript 库和框架编译成 Web Assembly 之外,我认为 Web 的未来不会 JavaScript 到 Web Assembly 编译器,而是从一些其他语言。对于 Web 编程的未来,我最喜欢的选择是 F#,也许 Fable 实现转换为生成 Web Assembly,而不是 JavaScript (asm.js) 或 Nim。甚至有可能生成 Web Assembly,它支持并显示极端循环展开的优势,非常接近已知最快的页面分段 SoE 速度。
结论
我们在 JavaScript 中构建了埃拉托色尼的页分段筛法,它适用于 table 以筛分数十亿的大范围,并且有进一步扩展这项工作的方法。生成的代码的剔除操作减少了大约十倍(当完全轮子分解时)并且剔除操作快了大约三倍,这意味着每个给定(大)范围的代码快了大约 30 倍,而减少的内存使用意味着可以筛选到大约 9e15 的 53 位尾数的范围在一年左右(只需打开 Chrome 浏览器选项卡并备份电源)。
虽然还有一些其他可能的小调整,但这是关于使用 JavaScript 筛选素数的最新技术:虽然由于给定原因它不如本机代码快,但它很快足以在一两天内找到从 1e14 到 1e14 的素数数量(即使是在中等范围的 smartphone 上),只需在所需的计算时间内打开浏览器选项卡即可;这是非常惊人的,因为直到 1985 年才知道这个范围内的素数数量,然后通过使用数值分析技术而不是使用筛子,因为当时的计算机使用最快的编码技术还不够快在合理和经济的时间内。虽然我们可以在短短几个小时内使用这些算法在现代台式计算机上使用最好的本机代码编译器来完成此操作,但现在我们可以在 acceptable 时间内使用 JavaScript 在智能计算机上完成此操作phone!
这扩展了
第 4.5a 章 - 全轮分解的基础知识
上一个答案中第 3 章的 non-Maximally-Wheel-Factorized Page-Segmented 埃拉托色尼筛法版本被编写为素数生成器,输出递归反馈作为基本素数馈送的输入;虽然这非常优雅且可扩展,但在接下来的工作中,我退后一步,使用更命令式的代码风格,以便读者可以更容易地理解最大轮分解的核心概念。在未来的第 4.5b 章中,我会将以下示例中开发的概念重新组合到质数生成器样式中,并添加一些额外的改进,这些改进不会使其在几十亿的较小范围内更快,但将使该概念在没有大部分损失速度可达数万亿至数百或数千万亿;随着范围变大,素数生成器格式在使程序适应方面更有用。
以下示例的主要额外改进是在用于有效寻址车轮模数残差的各种查找表 (LUT) 中,特殊起始地址 LUT 的生成非常简单地允许人们计算剔除起始地址对于每个模残差位平面,给定起始地址轮索引和整个结构中第一个剔除的第一个模残差位平面索引,没有额外的整数除法(慢),以及使用这些的 Sieve Buffer 复合数表示剔除函数。
此示例基于使用二、三、五和七小素数的 210 数字圆周轮,因为它似乎达到了数组大小和位数效率的“最佳点”平面,但实验表明,通过将下一个素数 11 添加到圆周为 2310 个数字的轮子中,可以将性能再提高 5%;没有这样做的原因是初始化时间大大增加,并且很难为“只有”十亿的较小范围计时,因为只有大约四个段才能达到这一点,粒度成为一个问题。
请注意,第一个筛选的数字是 23,这是轮素数和 pre-cull 个素数之后的第一个素数;通过使用它,我们避免了处理从“一”开始的数组的问题以及恢复被消除并且必须通过某些算法重新添加的轮素数的问题。
EDIT_ADD: 添加一个词来解释各种 WHL LUT 的用途 - 其中大部分与减少对 [=116= 的需求有关] 昂贵的除法操作到每个大页面段的每个基本质数只有大约两个。如下所述,WHLSTRTS LUT 用于将起始地址转换为筛选页段(在本例中为 48)的每个基本质数每个剩余模数索引所需的位索引,以按原样进行非常简单的查找、乘法和加法运算如下面所描述的。 END_EDIT_ADD
基本上,对于每个页面段剔除扫描,都有一个起始循环,该循环使用轮索引和段内第一个剔除地址的模残基索引填充起始地址数组,用于每个小于页段中表示的最大数的平方根,然后使用此起始地址数组依次完全筛选每个模余数位平面(其中 48 个),并为每个位平面扫描所有基本素数,并从该段计算出适当的偏移量通过使用 WHLSTARTS LUT 中的乘数和偏移量来计算每个基本质数的起始地址。这是通过将基素数的轮索引与查找的乘数相乘并加上查找的偏移量以获得给定模余数位平面的轮索引起始位置来完成的。从这里开始,每位平面的剔除就像第三章奇数位平面的剔除一样。这对每个位平面完成 48 次,但是 16 KB 缓冲区(每个位平面)的每个页面段的有效范围是 210 轮跨度的 131072 倍或每个页面段 27,525,120 个数字,并从更大的范围线性乘以位平面的大小。
与第 3 章 odds-only 筛法相比,使用此筛法将内存使用量减少为 48 倍于 105 倍或不到一半,但是因为每个段都有 48 位平面,所以完整的 Sieve Buffer 是 48 位平面的 16 KB 或 768 KB(四分之三兆字节)和更大尺寸的倍数。然而,使用这种大小的 Sieve Buffer 的效率最高只能达到 160 亿左右,我们下一章的下一个示例将调整缓冲区的大小以适应巨大的范围,使其增长到大约 100 兆字节 fr 最大的范围。对于 multi-threaded 种语言(不是 JavaScript),这是每个线程的要求。
额外的内存要求用于存储 32 位值的基本素数数组,这些值表示基本素数的轮索引及其模余数索引(如上所述,模地址计算是必需的);对于 10 亿的范围,大约有 3600 个基本素数乘以四个字节,每个大约 14,000 个字节,附加的起始地址数组大小相同。这些数组随着要筛选的最大值的平方根而增长,因此对于筛选到 10^16(万亿)或大约每个 23 兆字节。
进一步细化适用于以下示例,使用“组合”筛子,其中筛子缓冲区 pre-filled 来自更大的轮子模式,从中可以得到 11、13、17、11、13、17、十九人,被淘汰;比这更大的消除范围是不切实际的,因为保存的 pre-culled 模式从每个模数位平面仅约 64 KB 增长到 48 个模数残基数平面中每个平面的大约 20 倍或大约 1.5 兆字节或仅通过添加素数 23 的额外因子约 60 兆字节 - 同样,这在内存和初始化方面的成本相当大,而性能仅占百分之几。请注意,此数组可以共享以供所有线程使用。
实施时,WHLPTRN 阵列约为 64 KB 乘以 48 模位平面约为 3 MB,这不是那么大,并且是固定大小,不会随着筛分范围的增加而改变;就访问速度和初始化时间而言,这是一个非常可行的大小。
这些“最大轮分解”改进减少了用于筛选十亿范围的合数剔除操作循环的总数,从第 3 章 odds-only 示例的大约十亿操作减少到这个“组合”筛子,目的是尽量使每个剔除操作的 CPU 时钟周期数保持相同,从而使速度提高四倍。
已编辑: 已调整以下代码段以添加基本的 HTML 单网页应用程序用户界面,以便可以轻松调整参数以进行实验。为了最好的易用性,点击“运行 代码段”按钮后,应使用右上角的“整页”link,并且可以使用右上角的 link 关闭整页等结束了。要在智能手机上 运行(最好在 Chrome 中),请使用设置菜单(三点菜单)中的“桌面站点”复选框。
EDIT_CORRECTION: 范围限制可以在指定的上限内轻松更改,虚拟基基索引的基素数数组大小不再足够覆盖 LIMIT 的指定上限 362(以前只有 229)的平方根的平方根,因此已增加到两个轮距或 439。
FURTHER_EDIT_CORRECTION: fillSieveBuffer函数在填充大于16384字节的SieveBuffer剩余位平面缓冲区时出现轻微错误,已更正
SPEED_OMISSION_CORRECTION: 从与其他语言的合作中,我们意识到由于没有有效地使用“loop unpeeling”,该版本比应有的速度慢了大约 20% “至于在应该应用的地方不计算适当的限制。已添加并应用“bplmt”来更正此问题。第一次 运行 编译代码时,应多次按下 Sieve 按钮,让 JavaScript 引擎对生成的代码进行热调以进行优化,从而提高速度,它将在四五次迭代后达到.
EDIT_POLISHING: 看来把sieve buffer size设置成CPU L1 cache size确实没什么优势如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>
请注意,上面的代码仅适用于 non-trivial 超过十亿 (1e9) 的筛分范围,因为相当大的筛分页面尺寸的粒度使得时间不那么准确范围低于此尺寸。 该程序最适用于 1E11 及以上的筛分范围。
由于以下原因,上述代码可能没有预期的那么快:
使用具有固定掩码模式的特殊简化剔除循环的加速技术不再有效,因为几乎不再有任何小的剔除基素数;这将每个合数剔除操作的平均时钟周期数增加了大约 20%,尽管这更适用于 JavaScript 等较慢的语言,而不是更高效的机器代码编译语言,因为它们可以使用更进一步的技术,例如 extreme循环展开以进一步将每个剔除操作循环的 CPU 时钟周期数减少到每个剔除循环大约 1.25 个时钟周期。
尽管由于较少的模位平面(大约那个因子),计算结果素数的开销减少了大约两倍,但这不是所需的四倍;这在使用 JavaScript 时变得更糟,它无法利用 CPU POP_COUNT 机器指令,这比使用的计数 LUT (CLUT) 技术快大约十倍这里。
虽然此处使用的 LUT 技术将起始地址计算开销减少了五分之一左右,而对于所需的更复杂的模计算,这些计算大约是两倍半比第 3 章中的“odds-only”筛法要求的复杂三倍,因此不足以减少 ratio-metric;我们需要一种技术将时间减少两倍左右,以便这些计算不会对减少的比率做出贡献。相信我,我试过了,但没能做得更好。也就是说,这些计算在更高效的计算机语言中可能比 JavaScript and/or 更好 CPU 比我非常低端的 Atom CPU 处理器更有效,但是合数剔除速度也可能更高效!
不过,速度几乎提高了四倍,代码行数只增加了大约 50%,这还算不错,不是吗?当 运行 在较新版本的 node/Google Chrome 浏览器上(Chrome 版本 75 仍然比 Firefox 版本 68 快 25%) 比 Kim Walisch 写的“primesieve” “C”并编译为 x86_64 本机代码!
我添加了一个
即将到来的是第 4.5b 章,这将是大约两倍半的代码,但将是一个能够筛选到极大范围的素数生成器,部分受限于 JavaScript 只能有效地表示高达 53 位的 64 位浮点位尾数或大约 9e15 的数字,以及人们想要等待的时间:在更好的 CPU 上筛选到 1e14 将占用一天的顺序,但这不是什么大问题 - 只需打开浏览器选项卡即可!
附录:使用其他JavaScript生成语言
最多
- JavaScript代码是很久以前设计的,所以它的编程模型已经过时了。这在很大程度上体现在它对遗留代码的支持,尽管它已经根据 ECMA2015 和更新版本进行了(大量)更新。然而,更新导致了太多的做事方式,并不是所有的方式都有效,让程序员对最好的方式感到困惑。
- JavaScript 使用原型模型支持面向对象编程 (OOP) 版本的方式非常糟糕!
- JavaScript 是动态类型的,这可能会导致意外的代码错误,并使维护和重构代码变得更加困难。
- 为了速度而手动编码(通常使用 asm.js 和更新的功能)并不总能产生最佳代码,除非人们非常了解 JavaScript 引擎使用的这些优化,这些优化实际上会执行代码。
有几个主流选项可以使用另一种语言通过转译生成 JavaScript,其中两个最常见(也是最好)如下:
- Microsoft 的 Typescript 是一种(可选)静态类型的 OOP 语言,由于上述原因而变得非常流行。
- 我最喜欢的是 Fable,它是微软主要基于 ML 的函数式语言 F# 的一个分支,开发用于生成高效的 JavaScript 代码,并且是一种更好的静态类型检查语言,但它提供类型推断以及各种简洁的功能。
如果我能再避免它,我不会做 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# 中实现 monod 的想法,但它们的性能不佳,因为没有 Haskell 具有的自定义优化)。
CORRECTIONS_TO_COMMENTS: Tt 结果表明 Chrome V8 JavaScript 引擎似乎已经优化了不可变剔除突变和差异第一次更改的速度是由于使用了命令式 JavaScript for/while 循环,而不是 Fable 使用循环模拟尾递归函数,后者速度稍慢。最好的数组复制直接调用 JavaScript 的最大效果。使用不同的位长进行计数也是一个相当小的改进。总改进约为 25%,但复制改进与其他两个组合的改进大致相同。
在按上述link打开的页面中,可以查看生成的JavaScript代码,会看到生成的纯asm.js代码更好(更一致) ) 比手写的要多,但是 Fable 代码需要一点性能帮助,迫使它在三个地方发出 JavaScript 代码,如下所示:
- Fable/F# 代码默认是不可变的,为了忠实于函数精神,我使用不可变的尾调用递归函数循环,但这些比使用命令式 for/while 循环稍微慢一些。
- 数组 copy/blit 的寓言库没有(至少)使用最快的 JavaScript 方法(设置 subarray/slice),所以我强制它发出 JavaScript 这样做。
- 在计算剔除数组中未设置的位时,Fable 不提供 TypedArray 视图的其他位宽(我不认为),所以我通过发射 JavaScript 添加它作为它对 Uint32 的使用查看比使用四个连续的 Uint8 读取和通过位移手动组装 Uint32 更快。
您会发现生成的代码与第 4a 章中手写的 JavaScript 代码一样快!