在Javascript中实现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;
};
我似乎无法弄清楚我哪里弄错了.
可能一直在工作太久了.
例如：
筛(4,10)应该返回[5,7]
但它正在回归[5,7,9]
解决方法:
虽然您已经从阅读中得知,页面分段的Eratosthenes筛选是一种快速查找大范围质数的方法,但您的问题代码(即使已更正)也未实现页面分段SoE,请在大范围内测试代码,SoE实施也没有特别快.以下讨论将展示如何在大范围内使用真正的页面分段SoE.
概要
以下是导致您意图的越来越快的实施的阶段性进展,评论解释了每个步骤的原因和实施细节.它包含JavaScript中的可运行片段,但这些技术不仅限于JavaScript,其他语言也不限制某些进一步的细化,例如生成页面的多线程(除了Web Workers之外,它们很难控制为为了处理的顺序),在极端循环展开中进行了一些进一步的优化,并且由于必须由浏览器中的JavaScript引擎编译为本机代码的Just In Time(JIT),最后与代码的有限效率有关;这些限制与直接编译为非常有效的本机代码的语言相比,例如C/C++,Nim,Rust,Free Pascal,Haskell,Julia等.
第1章 – 设置基线
首先,让我们从当前代码算法的工作版本开始,在相当大的范围内使用定时信息来建立基线;下面的代码开始在剔除素数的平方上剔除每个素数,这避免了剔除给定素数值和一些冗余起始剔除的问题,并且没有理由为大范围生成结果素数的输出数组,因为我们可以生成直接来自剔除阵列的素数;此外,答案的确定是在时间之外,因为我们将开发更好的技术来找到“答案”,对于大范围通常是找到的素数的数量,素数的总和,第一次出现的素数差距等,其中没有一个需要实际查看找到的素数：
"use strict";
function soePrimesTo(limit) {
var sz = limit - 1;
var cmpsts = new Uint8Array(sz); // index 0 represents 2; sz-1 is limit
// no need to zero above composites array; zeroed on creation...
for (var p = 2; ; ++p) {
var sqr = p * p;
if (sqr > limit) break; // while p is the square root of limit -> cull...
if (cmpsts[p - 2] == 0 >>> 0) // 0/1 is false/true; false means prime...
for (var c = sqr - 2; c < sz; c += p) // set true for composites...
cmpsts[c] = 1 >>> 0; // use asm.js for some extra efficiency...
}
var bi = 0
return function () {
while (bi < sz && cmpsts[bi] != 0 >>> 0) ++bi;
if (bi >= sz) return null;
return bi++ + 2;
};
}
// show it works...
var gen = soePrimesTo(100);
var p = gen();
var output = [];
while (p != null) { output.push(p); p = gen(); }
console.log("Primes to 100 are: " + output + ".");
var n = 1000000000; // set the range here...
var elpsd = -new Date();
gen = soePrimesTo(n);
elpsd += +new Date();
var count = 0;
while (gen() != null) { count++; }
console.log("Found " + count + " primes up to " + n + " in " + elpsd + " milliseconds.");现在,将上述代码的运行时间引用到十亿分之一的目的几乎没有用,因为我在英特尔x5-Z8350中运行1.92千兆赫(CPU时钟)时使用非常低端的平板电脑CPU几乎肯定会更快单个活动线程的速度).我将引用我的执行时间一次作为如何计算每个剔除的平均CPU时钟周期的示例：我将上述代码的执行时间约为43350毫秒,是43.35秒,每秒19.2亿个时钟除以约2.514十亿次剔除操作筛选到十亿(可以从公式计算或从表中计算出无车轮分解on the SoE Wikipedia page),每次剔除大约33.1个CPU时钟周期到十亿.从现在开始,我们将仅使用每个剔除的CPU时钟周期来比较性能.
进一步注意,这些性能分数非常依赖于所使用的浏览器JavaScript引擎,上述分数在Google Chrome上运行(版本72);对于我们正在移动的页面分段SoE版本,Microsoft Edge(版本44)大约慢七倍,而Firefox和Safari可能在性能上接近Chrome.
由于使用Uint8Array TypedArray和更多asm.js,这种性能可能比之前的答案代码更好,但是这种“一个巨大的阵列筛”(这里使用的GigaByte内存)的时间由于内存访问速度而受到瓶颈对于CPU缓存限制之外的主RAM内存.这就是为什么我们正在努力实现Page Segmented Sieve,但首先让我们做一些关于减少使用的内存量和所需剔除周期数的事情.
第2章 – 比特打包和赔率轮分解
下面的代码进行了位打包,这需要在紧密的内部剔除循环中稍微复杂一些,但由于每个复合数字只使用一位,因此它将内存使用量减少了8倍;同样,由于两个是唯一的偶数素数,它使用基本的轮分解来筛选几率,仅将内存使用量进一步减少两倍,并将剔除操作次数减少约2.5倍.
这种奇偶校验的最小轮分解工作原理如下：
>我们制作了一个“轮子”,上面有两个位置,一半是在偶数上“击中”,另一半是在奇数上击中;因此,当我们在所有潜在的素数上“滚动”它时,这个“轮”具有两个数字的圆周跨度,但它只“窗口”中的两个或一半数字“滚动”.
>然后我们将我们“滚动”的所有数字分成连续值的两个位平面,在一个位平面中,我们将所有的平均值从4开始,而在另一个位平面中,我们将所有的赔率从3开始.
>现在我们丢弃偶数平面,因为没有代表的数字可以是素数.
>对于所考虑的奇基本素数,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.");现在,每次剔除操作的性能约为34.75个CPU时钟周期,1025亿次剔除操作仅为10亿美元(来自维基百科),这意味着时间的减少几乎完全是由于剔除操作次数的减少,由于存储器使用减少16倍,因此“bit-twiddling”位填充的额外复杂性仅占用与节省相同的额外时间.因此,这个版本使用了十六分之一的内存,速度提高了约2.5倍.
但是我们还没有完成,页面细分可以加快我们的速度,就像你的消息来源告诉你的那样.
第3章 – 页面分割包括更快的Primes计数
那么什么是页面分割应用于SoE以及它对我们有什么作用？
页面分割将筛选工作从一次筛分的一个巨大阵列划分为一系列连续筛分的小页面.然后需要更复杂一点,因为必须有一个单独的基本素数流可用,这可以通过用内筛递归筛选来获得,该内筛产生主筛用于的记忆基础素数列表.同样,输出结果的生成通常稍微复杂一些,因为它涉及对每个生成的筛分页面的连续扫描和缩减.
页面分割具有以下许多好处：
>它进一步降低了内存要求.使用之前的赔率筛选到十亿美元的代码需要大约64兆字节的RAM,但是不能使用该算法来筛选超过16到320亿的范围(即使我们想要等待很长时间)到期使用JavaScript可用的所有内存.通过筛选(例如)该数量的平方根的页面大小,我们可以筛选到该数量的范围或超出我们将有时间等待的任何数量.我们还必须将基本素数存储到所需范围的平方根,但这只是我们需要考虑的任何范围的几十兆字节.
>它增加了内存访问速度,这直接影响了每次剔除操作的CPU克拉周期数.如果剔除操作主要发生在CPU缓存中,则内存访问从主RAM内存的每次访问的数十个CPU时钟周期(取决于CPU和RAM的质量和类型)变为CPU L1缓存的大约一个时钟周期(更小)从大约16 KiloByte到大约8个时钟周期来自CPU L2缓存(大约至少大约128 KiloByte)并且我们可以计算出剔除算法,以便我们使用所有缓存达到最佳效果,使用小型快速L1缓存大多数具有小基本素数值的操作,对于中等大小的基本素数而言,较大的小位慢二级高速缓存,并且仅针对极大范围使用大基本素数的少数操作使用主RAM访问.
>它通过将每个较大页面的筛选工作分配给不同的CPU核心来实现相当粗粒度的并发性,从而开启了多线程的可能性,尽管除了通过使用Web Workers(凌乱)之外,这不适用于JavaScript.
除了增加的复杂性之外,页面分割还有另外一个问题需要解决：与“一个巨大的阵列”筛子不同,其中起始索引很容易计算一次,然后用于整个阵列,分段筛需要通过以下方式计算起始地址：每页每个素数的模数(除法)运算(计算成本很高),或者需要使用额外的内存来存储每页每个基本素数达到的索引,因此不必重新计算起始索引,但是这最后一种技术排除了多个这些数组不同步时进行线程化.在Ultimate版本中使用的最佳解决方案是使用这些技术的组合,其中多个页面段被分组以形成用于线程的相当大的工作单元,使得这些地址计算占总时间的可接受的一小部分,并且索引存储表用于每个线程的这些较大工作单元的基本素数,因此每个较大的工作单元只需要进行一次复杂的计算.因此,我们既可以获得多线程,也可以减少开销.但是,以下代码不会减少这种开销,当筛选到十亿时,这会花费大约10％到20％.
除了页面分割之外,以下代码通过使用计数查找表(CLUT)填充计数算法添加对找到的素数的有效计数,该算法一次计数32位,以便连续查找数字的结果的开销发现素数占筛选时间的一小部分.如果没有这样做,枚举个别找到的素数只是为了确定在给定范围内筛选至少需要多少时间.可以很容易地开发类似的快速例程来执行诸如对找到的素数求和,找到素数间隙等的事情.
START_EDIT：
下面的代码增加了另一个加速：对于较小的素数(此优化有效),代码通过识别剔除操作遵循八步模式来执行循环分离的形式.这是因为一个字节具有偶数个位,我们正在通过奇数素数进行剔除,这些素数将每八个剔除返回一个字节中的相同位位置;这意味着对于每个位位置,我们可以简化内部剔除循环以屏蔽常量位,从而大大简化内部循环并使剔除速度提高约两倍,因为模式中的每个剔除循环都不需要执行“bit-twiddling”位打包操作.这一变化节省了大约35％的执行时间,筛选到10亿.它可以通过将64更改为0来禁用.它还设置了本机代码极端展开八个循环的阶段,因为这种模式可以在使用本机代码编译器时将剔除操作速度提高大约两倍.
通过使用查找表(LUT)作为掩码值而不是左移操作来进行进一步的微小修改,使得更大的素数(大于8192)的循环更快,以便平均每次剔除操作节省大约半个CPU时钟周期剔除十亿的范围;随着范围从十亿增加,这种节省将略有增加,但在JavaScript中并没有那么有效,并且已被删除.
END_EDIT
ANOTHER_EDIT：
除了上述编辑之外,我们已经删除了LUT BITMASK,但现在通过从相同大小的零缓冲区快速字节复制将Sieve Buffer归零,并添加了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.");正如此处所实现的那样,该代码每次筛选大约需要13.8个CPU时钟周期才能达到10亿次.当使用以下最大轮分解算法时,由于有效页面大小增加了105倍,因此通过改进地址计算算法节省了大约20％的额外费用,因此这种开销只有几个百分点,可与几个百分比相媲美用于数组填充和结果计数.
第4章 – 进一步高级工作添加最大轮分解
现在我们考虑使用最大车轮分解的更广泛的变化(不仅仅是“仅赔率”2,而且对于覆盖210个潜在素数的跨度而不是仅仅2的跨度的3,5和7并且还预先剔除小筛分阵列的初始化,这样就不必剔除11,13,17和19的以下质数.这减少了使用页面分段筛时的复合数剔除操作的数量通过大约4到10亿的范围(如表中所示/根据维基百科文章中的公式计算 – “组合轮”)并且可以编写,因为它减少了大约四倍的速度每次剔除操作的操作大约与上述代码相同.
有效地进行210跨度车轮分解的方法是遵循“仅有赔率”的方法：上面的当前算法可以被认为是在两个中包含一个比特填充的平面,如前一章中所解释的那样.飞机可以被淘汰,因为它只包含两个以上的偶数;对于210-span,我们可以定义48个这种大小的位填充数组,表示11和11以上的可能素数,其中所有其他162个平面包含的数字是2,3,5或7的因子,因此不需要被考虑.只能通过基本素数的增量重复索引来剔除每个位平面,就像奇数平面完成一样,乘法自动由结构处理,就像仅有几率一样.通过这种方式,对于具有较少内存要求的筛分效率(与1/2的“仅赔率”相比为48/210)和效率与仅为几率的效率相同,其中一个48平面“页面”代表每平面尺寸的16千字节= 131072比特,210每个筛网页面的数量范围为27,525,120,因此只有近40页的分段可以筛选到十亿(而不是上面的近四千),因此减少了开销.每页基本素数的起始地址计算,以进一步提高效率.
虽然上面描述的扩展代码是几百行而且很长,但是在我的低端英特尔1.92 Gigahertz CPU上使用Google V8 JavaScript引擎可以计算大约两秒内的素数到10亿.比在本机代码中运行的相同算法慢五倍.
虽然上述代码在大约160亿的范围内非常有效,但是其他改进可以帮助将效率保持在数万亿甚至更大的范围内,例如1e14或更多.我们通过向上调整有效页面大小来实现这一点,这样它们永远不会小于被筛选的整个范围的平方根,而是逐步筛选16个KiloByte块用于小素数,128个KiloByte块用于中等素数,并且仅筛选根据我们的基线实现,对于用于最大基本质量大小的极少数剔除操作,这是一个巨大的数组.通过这种方式,对于我们可能考虑的最大范围,我们的每次剔除时钟不会增加超过约两倍的小因子.
由于这个答案接近30,000个字符的有限大小,因此在my followup Chapter 4.5a和(未来)第4.5b章的上述技术的实际实现的答案中继续进一步讨论最大轮分解.
第5章 – JavaScript(和其他VM语言)不能做什么
对于JavaScript和其他虚拟机语言,最小剔除循环时间大约为每个剔除循环10个CPU周期,并且不太可能发生太大变化.这比使用相同算法(如C/C++,Nim,Rust,Free Pascal,Haskell,Julia)直接编译为高效机器代码的语言可以轻松实现的大约三个CPU时钟周期慢大约三到四倍.等等
此外,还有极端的循环展开技术,这些技术可以与至少某些语言一起使用,这些技术可以使平均剔除操作周期减少大约两倍,而JavaScript则被拒绝.
多线程可以减少大约使用有效CPU核心因素的执行时间,但使用JavaScript,获得此功能的唯一方法是使用Web Workers并且同步混乱.在我的机器上,我有四个内核,但由于CPU时钟速率降低到四分之三而所有内核都处于活动状态,因此速度只增加了三倍;这个因素不容易获得JavaScript.
所以这是关于使用JavaScript和其他当前VM语言的最新技术,除了它们可以轻松使用多线程之外,它们具有相同的限制,结合上述因素意味着本机代码编译器可以大约二十次比JavaScript更快(包括多线程,甚至更多的新CPU拥有大量内核).
但是,我相信三到五年内Web编程的未来将是Web Assembly,并且有可能克服所有这些限制.它现在非常接近支持多线程,虽然目前在Chrome上这个算法的速度比JavaScript快了约30％,但是当使用某些Web从某些语言编译时,它在某些当前浏览器上的速度仅比本机代码慢一点.装配编译器.现在仍然处于Web Assembly的高效编译器和高效的浏览器编译到本机代码的早期开发,但由于Web Assembly比大多数VM更接近本机代码,因此可以轻松地改进它以生成快速或接近的本机代码与其他通知代码编译语言的代码一样快.
但是,除了将JavaScript库和框架编译到Web程序集之外,我不相信Web的未来将是JavaScript到Web程序集编译器,而是从其他语言编译.我最喜欢的Web编程未来的选择是F#,可能将Fable实现转换为生成Web Assembly而不是JavaScript(asm.js)或Nim.甚至有可能生成Web Assembly,支持并显示极端循环展开的优势,非常接近最快的已知页面分段SoE速度.
结论
我们在JavaScript中构建了一个页面分段的Eratosthenes Sieve,适用于筛选数十亿的大范围,并且有进一步扩展这项工作的方法.结果代码的剔除操作(完全车轮分解)大约减少十倍,剔除操作快约三倍意味着每个给定(大)范围的代码快约30倍,而减少内存使用意味着可以筛选到大约9e15的53位尾数的范围大约为一年(只需打开Chrome浏览器标签并备份电源).
虽然还有其他一些小的调整可能,但这是关于使用JavaScript筛选素数的最新技术：虽然由于给定的原因它不是一个快速的本机代码,它足够快到找到1e14的素数数量在一天或两天内(即使在中档智能手机上),打开浏览器标签以获得所需的计算时间;这是非常令人惊讶的,因为这个范围内的素数的数量直到1985年才知道,然后通过使用数值分析技术而不是使用筛子,因为那天的计算机使用最快的编码技术来做得不够快在合理和经济的时间内.虽然我们可以在短短几个小时内使用这些算法为现代台式计算机上最好的本机代码编译器做到这一点,但现在我们可以在智能手机上使用JavaScript在可接受的时间内完成！
标签：sieve-of-eratosthenes,sieve,javascript,primes 来源： https://codeday.me/bug/20191001/1839647.html