Aks*_*pal 2 javascript primes sieve sieve-of-eratosthenes
我最近阅读了有关大量数字的更快的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;
};
Run Code Online (Sandbox Code Playgroud)
我似乎无法弄清楚哪里弄错了。可能工作了太久了。
例如:
sieve(4,10)应该回来[5,7]
但它正在回来[5,7,9]
Gor*_*ood 10
这扩展了我之前的答案,即添加承诺但在每个答案限制 30,000 个字符中没有空间的内容:
上一个答案中第 3 章的 Eratosthenes 版本的非最大轮因式分页筛选被编写为质数生成器,其输出递归反馈为基质数馈送的输入;尽管这非常优雅且可扩展,但在接下来的工作中,我退回到了更命令式的代码风格,以便读者可以更容易地理解最大轮因式分解的核心概念。在以后的第 4.5b 章中,我会将以下示例中开发的概念重新组合成 Prime Generator 样式,并添加一些额外的改进,这些改进不会使其在数十亿的较小范围内变得更快,但会使该概念在没有大部分损失速度可达数万亿至数百或数千万亿;
以下示例的主要额外改进是在用于有效寻址轮模余数的各种查找表 (LUT) 中,特殊起始地址 LUT 的生成非常简单地允许人们计算每个模余数位的剔除起始地址平面给出了整个结构中第一个剔除的起始地址轮索引和第一个模余数位平面索引,没有额外的整数除法(慢),以及使用这些的筛选缓冲区复合数表示剔除函数。
此示例基于使用 2、3、5 和 7 的小素数的 210 个数字圆周轮,因为它似乎达到了阵列大小和位平面数量的效率“最佳点”,但实验表明通过将 11 的下一个素数添加到 2310 个数字的圆周上,可以再获得大约 5% 的性能;没有这样做的原因是初始化时间大大增加,并且很难为“只有”十亿的较小范围计时,因为只有大约四个段才能到达那个点,并且粒度成为一个问题。
请注意,第一个筛选的数字是 23,它是经过轮素数和预剔除素数后的第一个素数;通过使用它,我们避免了处理从“一”开始的数组的问题,以及恢复被某些算法消除并必须重新添加的轮素数的问题。
EDIT_ADD:添加对各种 WHL LUT 目的的解释 - 其中大部分与将计算时间昂贵的除法运算的需求减少到每个大页段的每个基素数仅约两个有关。如下所述,WHLSTRTS LUT 用于将起始地址转换为筛页段(在本例中为 48)的每个基素数所需的位索引,以进行非常简单的查找、乘法和加法操作如下面所描述的。 END_EDIT_ADD
基本上,对于每个页段剔除扫描,都有一个起始循环,该循环使用小于最大值平方根的每个基素数的段内第一个剔除地址的轮索引和模残数索引填充起始地址数组页段中表示的数字,然后使用此起始地址阵列依次完全筛选每个模余数位平面(其中 48 个),扫描每个位平面的所有基本素数,并根据每个基本素数的段起始地址计算出适当的偏移量通过使用 WHLSTARTS LUT 中的乘数和偏移量。这是通过将基本素数的轮索引与查找乘数相乘并添加查找偏移量以获得给定模残差位平面的轮索引起始位置来完成的。从此,每位平面剔除就像第三章奇数位平面一样。这对每个位平面进行了 48 次,但对于 16 KB 缓冲区(每个位平面),每个页段的有效范围是 210 轮跨度的 131072 倍或每个页面段的 27,525,120 个数字,并且从该值线性乘以更大的位平面的大小。
与第 3 章仅赔率筛选相比,使用此筛选将内存使用量减少为 105 的 48 倍或不到一半,但因为每个段都有全部 48位平面,完整的 Sieve Buffer 是 48 位平面的 16 KB 乘以 768 KB(兆字节的四分之三),并且是更大尺寸的倍数。然而,使用这种大小的 Sieve Buffer 最多只能有效处理 160 亿,我们下一章中的下一个示例将调整缓冲区的大小以适应大范围,以便最大范围增长到大约 100 兆字节。对于多线程语言(不是 JavaScript),这将是每个线程的要求。
额外的内存要求用于存储 32 位值的基本素数数组,这些值表示基本素数的轮索引及其模余数索引(对于如上所述的模地址计算是必需的);对于十亿的范围,大约有 3600 个基本素数乘以 4 个字节,每个大约 14,000 个字节,附加的起始地址数组的大小相同。这些数组随着要筛选的最大值的平方根增长,因此对于小于一亿的基本素数,筛选到 10^16(一万万亿)或大约每个 23 兆字节。
进一步的改进适用于以下使用“组合”筛子的示例,其中筛子缓冲区是从更大的轮子模式中预填充的,其中已经消除了 11、13、17 和 19 的素数因子;比这更大的消除范围是不切实际的,因为保存的预剔除模式从每个模位平面仅大约 64 KB 增长到大约 20 倍,或者对于 48 个模残数平面中的每个平面或大约 60 个大约 1.5 兆字节仅仅通过添加 23 的质数的额外因素来实现兆字节 - 同样,它在内存和初始化方面的成本相当大,而性能仅占百分之几。请注意,此数组可以共享以供所有线程使用。
在实施时,WHLPTRN 阵列大约是 64 千字节乘以 48 模位平面大约是 3 兆字节,这不是那么大,而且是固定大小,不会随着筛分范围的增加而改变;对于访问速度和初始化时间来说,这是一个非常可行的大小。
这些“最大轮因数分解”改进将用于筛选十亿范围的合数剔除操作循环的总数从第 3 章仅几率示例的大约十亿次操作减少到这个“组合”筛选器的大约 25 亿次操作,其中目标是尝试使每个剔除操作的 CPU 时钟周期数保持相同,从而使速度提高 4 倍。
已编辑: 以下代码段已进行调整以添加基本的 HTML 单网页应用程序用户界面,以便可以轻松调整参数以进行实验。为了最好的使用,点击“运行代码片段”按钮后应该使用右上角的“全页面”链接,完成后可以用右上角的链接关闭整个页面。要在智能手机上运行(最好在 Chrome 中),请使用设置菜单(三点菜单)中的“桌面站点”复选框。
EDIT_CORRECTION:由于可以在指定的上限内轻松更改范围限制,因此虚拟基础基础索引基础素数数组大小不再足以覆盖 362 的指定上限 LIMIT 的平方根的平方根(以前只有 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;Run Code Online (Sandbox Code Playgroud)
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;
}Run Code Online (Sandbox Code Playgroud)
<!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>Run Code Online (Sandbox Code Playgroud)
请注意,上述代码仅适用于超过 10 亿 (1e9) 的非平凡筛分范围,因为相当大的筛分页面大小的粒度使得对于低于此大小的范围计时不太准确。 该程序最适用于 1E11 及以上的筛分范围。
由于以下原因,上面的代码可能不像人们预期的那么快:
使用具有固定掩码模式的特殊简化剔除循环的加速技术不再有效,因为几乎不再有任何小的剔除基数;这将每个合数剔除操作的平均时钟周期数增加了约 20%,尽管这更适用于较慢的语言,如 JavaScript,而不是更高效的机器代码编译语言,因为它们可以使用进一步的技术,如极端循环展开来进一步将每个剔除操作循环的 CPU 时钟周期数减少到每个剔除循环约 1.25 个时钟周期。
尽管由于模位平面较少(大约减少了那个因素),计算结果素数的开销减少了大约 2 倍,但这不是所需的 4 倍;这在使用 JavaScript 时变得更糟,JavaScript 无法利用 CPU POP_COUNT 机器指令,这将比此处使用的计数 LUT (CLUT) 技术快十倍。
虽然此处使用的 LUT 技术将起始地址计算开销从所需的更复杂模计算的开销减少了大约五倍,但这些计算比所需的复杂模计算复杂了大约两倍半到三倍。第 3 章中的“仅赔率”筛选,因此不足以给我们减少比率度量;我们需要一种技术来将时间进一步减少两倍左右,以便这些计算不会对减少的比率做出贡献。相信我,我尝试过,但一直无法让这一切变得更好。也就是说,在比 JavaScript 和/或比我的极低端 Atom CPU 处理器更好的 CPU 中,这些计算可能在更高效的计算机语言中更高效,
尽管如此,几乎是四倍的加速,而代码行数仅增加了约 50% 还不错,不是吗?这段JavaScript代码是唯一的三到五倍,速度较慢(取决于CPU,在性能上更接近对高端CPU的),当节点/谷歌Chrome浏览器的新版本上运行(Chrome版本75仍是约25%,比火狐版本68更快)比Kim Walisch 的“primesieve”用“C”编写并编译为 x86_64 本机代码!
我添加了一个附录答案,表明当项目的大小增加到超过几百行时,实际上并不需要编写 JavaScript 来生成 JavaScript,如下所示。在未来,我预计像 Fable 这样的转译器可能会发出 WebAssembly 代码而不是 JavaScript,现在用 Fable 编写的优点是人们不必对代码进行更改(或至少很少更改)以获取新技术的优势,它支持更快的代码执行和(最终)JavaScript 不支持的多线程。
即将到来的第 4.5b 章将是大约两倍半的代码,但将是一个能够筛选极大范围的主要生成器,部分受限于 JavaScript 只能有效地表示高达 64 位浮点尾数 53 的数字位或大约 9e15 以及人们想要等待的时间:在更好的 CPU 上筛选到 1e14 需要一天的时间,但这不是什么大问题 - 只需打开浏览器选项卡即可!
尽管您从阅读中发现,Eratosthenes的页面分段筛是在较大范围内查找质数的快速方法,但您的问题代码(即使已更正)也未实现页面分段SoE,请在较大范围内测试代码,随着SoE实施的进行,它也不是特别快。以下讨论将显示如何在较大范围内逐步使用真正的页面分段SoE。
以下是逐步实现您的意图的阶段性进展,并通过注释解释了每个步骤的原因和实施细节。它包括JavaScript中的可运行代码段,但这些技术不仅限于JavaScript,其他语言也没有对进一步的改进施加限制,例如对结果页面进行多线程处理(Web Worker除外,Web Workers难以控制)到处理顺序),对极端循环展开进行了一些进一步的优化,而最后的优化与代码效率有限有关,这是因为必须由浏览器的JavaScript引擎将即时(JIT)编译为本地代码;这些限制与直接编译为非常高效的本地代码(例如C / C ++,Nim,Rust,Free Pascal,Haskell,Julia,
首先,让我们从在相当大范围内使用当前代码算法的工作版本开始,并结合时序信息来建立基线。下面的代码在剔除素数的平方处开始每个素数的剔除,这避免了剔除给定素数值和一些多余的起始剔除的问题,并且没有理由为我们产生的大范围生成素数的输出数组直接从剔除阵列中获得质数;同样,答案的确定不在时间范围内,因为我们将开发更好的技术来查找“答案”,对于较大范围的“答案”通常是找到的素数的数量,素数的总和,素数的首次出现缺口等,这些都不需要实际查看找到的素数:
"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.");Run Code Online (Sandbox Code Playgroud)
现在,引用上述代码过筛到十亿的运行时间几乎没有用,因为您的运行时间几乎肯定会更快,因为我在Intel x5-Z8350中使用的是低端平板电脑CPU,运行频率为1.92 GHz(CPU时钟)单个活动线程的速度)。我将仅引用一次执行时间作为计算每个剔除的平均CPU时钟周期的示例:我将上述代码的执行时间大约为43350毫秒乘以43.35秒乘以每秒19.2亿个时钟除以大约2.514 10亿次筛选操作筛选到10亿次筛选(可以从公式或从表上计算出来的SoE Wikipedia页面上的无轮子分解),每个剔除可获得约33.1 CPU时钟周期到十亿。从现在开始,我们将仅使用每个剔除的CPU时钟周期来比较性能。
还要注意,这些性能得分非常取决于所使用的浏览器JavaScript引擎,以上得分是在Google Chrome(版本72)上运行的;对于我们正在朝着的Page Segmented SoE版本而言,Microsoft Edge(版本44)要慢大约七倍,并且Firefox和Safari的性能可能接近于Chrome。
由于使用Uint8ArrayTypedArray和更多的asm.js,因此此性能可能比以前的答案代码更好,但是由于“内存分配速度”快,此类“一个巨大的数组筛子”(此处使用的千兆字节内存)的时间受到瓶颈的影响。主RAM内存超出CPU缓存限制。这就是为什么我们要努力实现“页面分段筛”的原因,但首先让我们做些关于减少内存使用量和所需剔除周期数的事情。
下面的代码进行位打包,这在紧密的内部剔除循环中需要稍微更多的复杂性,但是由于考虑到每个复合数字仅使用一位,因此将内存使用量减少了八倍;同样,由于两个是唯一的偶数素数,因此它仅使用基本的轮分解来筛选赔率,仅是将内存使用量进一步减少了两倍,而剔除操作的次数也减少了约2.5。
仅赔率的最小轮分解原理如下:
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.");Run Code Online (Sandbox Code Playgroud)
现在,对于10.257亿次剔除操作,每剔除操作的性能约为34.75个CPU时钟周期,仅十亿分之一的赔率范围(来自Wikipedia),这意味着时间的减少几乎完全是由于剔除操作次数的减少,由于“位纠结”位打包的额外复杂性,由于将内存使用减少了16倍,因此仅节省了大约与节省时间相同的额外时间。因此,此版本使用内存的十六分之一,并且快约2.5倍。
但是我们还没有完成,页面分段可以进一步加快我们的速度,就像您的消息来源告诉您的那样。
那么什么是页面细分应用于SoE,它对我们有什么作用?
页面分割将筛分工作从一次筛分的庞大数组划分为一系列连续筛分的较小页面。然后,这就需要更多的复杂性,因为必须有单独的基础素数流可用,可以通过使用内部筛子递归筛分生成可记忆的基础素数列表以供主筛子使用的方法来获得。同样,输出结果的生成通常要稍微复杂一点,因为它涉及到每个生成的筛选页面的连续扫描和缩小。
页面细分具有许多优点,如下所示:
除了增加的复杂性外,页面分段还有另一个问题要解决:与“一个巨大的数组”筛子不同,它可以轻松地一次计算开始索引,然后将其用于整个数组,分段筛子需要通过以下方法计算开始地址:每页每素数的取模(除法)操作(计算上很昂贵),或者需要使用额外的内存来存储每页每素数素数所达到的索引,因此不必重新计算起始索引,但是最后一种技术排除了多这些阵列不同步时发生线程化。旗舰版将使用的最佳解决方案是结合使用这些技术,其中将几个页面段组合在一起以形成一个相当大的线程工作单元,以便这些地址计算占用总时间的一小部分,并且索引存储表用于每个线程这些较大工作单元的基本素数,以便每个较大的工作单位只需进行一次复杂的计算。因此,我们既可以进行多线程处理,又可以减少开销。但是,以下代码不会减少此开销,筛分到十亿时,成本大约为10%到20%。
除了页面分段之外,以下代码通过使用一次计数32位的计数查找表(CLUT)填充计数算法来添加对找到的素数的有效计数,这样就可以连续查找计数结果的开销发现素数仅占筛分时间的一小部分。如果不这样做,则枚举单个找到的素数只是为了确定有多少素数,至少要花费给定范围内筛分的时间。可以很容易地开发出类似的快速例程来执行诸如求和发现素数,发现素数缺口等的工作。
START_EDIT:
以下代码增加了另一个速度:对于较小的质数(在此优化有效的情况下),该代码通过识别出剔除操作遵循八步模式来进行某种形式的循环分离。这是因为一个字节具有偶数个位数,并且我们通过奇数素数进行剔除,该奇数素数将每八个剔除返回一个字节中的相同位位置;这意味着对于每个位位置,我们都可以简化内部剔除循环以掩盖恒定的位,从而极大简化内部剔除循环,并且由于图案中的每个剔除循环都不需要进行剔除,因此剔除速度提高了大约两倍。 “位旋转”位打包操作。此更改节省了大约10%的执行时间,节省了35%的时间。它可以通过改变被禁用64,以0。由于这种模式,当使用本机代码编译器时,可以将剔除操作速度提高大约两倍,这也为八循环的本机代码极端展开奠定了基础。
进一步的细微修改通过使用查找表(LUT)作为掩码值而不是向左移位操作,使较大的质数(大于8192)的循环更快,从而平均每次清除操作可节省大约一半的CPU时钟周期。淘汰范围十亿;随着范围从十亿起,这种节省将略有增加,但是在JavaScript中效果不佳,因此已被删除。
END_EDIT
ANOTHER_EDIT:
除了上述编辑之外,我们还删除了LUT BITMASK,但现在通过从相同大小的零缓冲区进行快速字节复制来将筛分缓冲区归零,并添加了Counting LUT填充计数技术,总体速度提高了10%。
END_ANOTHER_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 ZEROSPTRN = new Uint8Array(16384);
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(131072);
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.");Run Code Online (Sandbox Code Playgroud)
如此处实现的那样,每个剔除筛分到十亿个代码需要大约13.8个CPU时钟周期。当使用以下最大轮子分解算法时,由于改进了地址计算算法而节省的约20%会自动得到改善,这是由于有效页面大小增加了105倍,因此该开销仅占百分之几,与百分之几相当用于数组填充和结果计数。
现在,我们考虑使用最大车轮分解系数的更广泛的更改(对于仅奇数的车轮,不仅2,对于覆盖210个潜在质数的车轮而不是2的车轮,分别为3、5和7) ),并预先筛选小筛分阵列的初始化,这样就不必剔除以下11、13、17和19的质数。这减少了使用分页筛网时复合数量剔除操作的数量。大约四分之一到十亿分之一的范围(如表中所示/根据Wikipedia文章中的公式-“组合轮”计算得出),并且由于运行时间的缩短,其运行速度大约是其四倍每次剔除操作的速度与上述代码大致相同。
有效地进行210跨度轮分解的方法是遵循“仅奇数”方法:上面的当前算法可以认为是将一个位填充的平面从两个平面中筛选出来,如上一章所述,而另一个平面可以消除,因为它仅包含两个以上的偶数;对于210跨度,我们可以定义48个这种大小的位数组,它们表示11或以上的可能素数,其中所有其他162平面都包含因数为2、3、5或7的数字,因此不需要被考虑。然后,仅通过基数素数的重复索引就可以剔除每个位平面,就像奇数平面是由结构自动处理乘法一样进行的,就像仅针对奇数一样。
尽管上面描述的扩展代码有几百行,并且要在此处发布很长的时间,但是在使用Google V8 JavaScript引擎的我的低端Intel 1.92 Gigahertz CPU上,它可以在大约两秒钟内将质数计算为十亿。比以本机代码运行的相同算法慢五倍。
尽管上面的代码在大约160亿范围内是非常有效的,但是其他改进可以帮助将效率保持在甚至几万亿的更大范围,例如1e14或更高。我们通过向上调整有效页面大小来实现此目的,以使它们永远不会小于被筛分的整个范围的平方根,而对于较小的质数,按16 KiloByte块递增筛选,对于中等质数的则按128 KiloByte块递增筛选,而仅对根据我们的基准实施方案,对于用于最大基本素数的极少数剔除操作,阵列非常庞大。这样,对于我们可能考虑的最大范围,我们的每只剔除时钟不会增加最多约2的小倍。
由于此答案接近于30,000个字符的有限大小,因此,在我的后续第4.5a章和(将来)第4.5b章中针对上述技术的实际实现,继续进行关于最大车轮分解的进一步讨论。
对于JavaScript和其他虚拟机语言,最小剔除循环时间约为每个剔除循环10个CPU周期,并且变化不大。这是一个关于慢三至四倍比可以很容易地与直接编译到诸如C / C ++,稔,锈,自由帕斯卡,Haskell中,朱使用相同的算法高效的机器代码语言实现的大约三个CPU时钟周期,等等
此外,还有一些极端的循环展开技术可以与这些语言中的至少一些一起使用,这些语言可以使平均剔除操作周期减少大约两倍,这是JavaScript所无法实现的。
多线程可以将执行时间减少大约所使用的有效CPU内核的倍数,但是使用JavaScript的唯一方法是通过使用Web Workers来实现,而且同步也很麻烦。在我的机器上,我有四个核心,但是由于所有核心都处于活动状态时,CPU时钟速率降低到四分之三,因此速度只能提高三倍。对于JavaScript而言,这三分之二不容易获得。
因此,这是关于使用JavaScript的最新技术,其他当前的VM语言具有相同的局限性,除了它们可以轻松地使用多线程外,上述因素的结合意味着本机代码编译器可以达到大约二十倍。比JavaScript更快(包括多线程,甚至在具有大量内核的新CPU上更快)。
但是,我相信三到五年内Web编程的未来将是Web Assembly,这有可能克服所有这些限制。现在,它非常接近支持多线程,尽管目前在Chrome上该算法仅比JavaScript快30%,但是当使用某些Web使用某些语言从某些语言编译时,它仅比某些当前浏览器中的本机代码慢一点。汇编编译器。对于Web汇编程序来说,高效的编译器和对本机代码的有效浏览器的编译仍处于发展初期,但是由于Web Assembly与大多数VM相比更接近于本机代码,因此可以很容易地对其进行改进,以产生与本机代码一样快或几乎相同的本机代码。速度与其他通知代码编译语言中的代码一样快。
但是,除了将JavaScript库和Frameworks编译为Web Assembly之外,我不认为Web的未来将是JavaScript到Web Assembly编译器,而是从其他某种语言进行编译。对于Web编程的未来,我最喜欢的选择是F#,也许将Fable实现转换为可生成Web Assembly而不是JavaScript(asm.js)或Nim。甚至有可能产生Web组件,以支持并显示极端循环展开的优势,使其非常接近最快的已知页面分段SoE速度。
我们使用JavaScript构建了Eratosthenes的页面分段筛网,该筛网适用于筛分数十亿的大范围,并具有进一步扩展这项工作的方法。生成的代码的剔除操作大约减少了十倍(完全轮分解时),剔除操作大约快了三倍,这意味着每个给定(较大)范围的代码快了大约30倍,而减少的内存使用意味着人们可以筛选大约9e15的53位尾数的范围(大约一年)(只需打开Chrome浏览器标签并备份电源)。
尽管还有其他一些细微的调整,但这是关于使用JavaScript筛选质数的最新技术:虽然由于给定的原因,它不如本机代码快,但它足以找到1e14的质数在一两天内(甚至在中端智能手机上)通过打开浏览器选项卡打开所需的计算时间;这非常令人惊讶,因为直到1985年才知道该范围内的质数,然后使用数值分析技术而不是使用筛子,因为当时的计算机使用最快的编码技术还不够快在合理且经济的时间内。尽管我们可以使用这些算法在短短几个小时内完成此操作,以获得现代台式计算机上最好的本机代码编译器,