
注意: @David_Burnham的挑战条目在截止日期之前提交,但被我的电子邮件系统隔离。我只是在原始博客文章发布后才发现它。我正在更新结果以反映基于新条目的变化。
感谢所有为挑战 2 提交代码的人。得分是一项挑战,因为有些代码在较小的值时效果更好,而有些代码在较大的值时效果更好。下表给出了挑战每个部分的结果。isPrime 运行了三次,1000 个不能被 5 整除的 5 位随机奇数、不能被 5 整除的 8 和 9 位奇数以及不能被 5 整除的 12 到 15 位奇数。结果基于最后两个条件。分数在挑战之间分配,第一名得 1.5 分,第二名得 1 分,第三名得 0.5 分。还给出了最快提交的时间(以秒为单位)。
@leelaim isPrime(Small) 的时间与 JMP 工作人员的任何时间一样快(+3 分)。 @landon isPrime(Big) 脚本使用 Parallel Assign 函数来提高速度。两个条目缓存了它们的结果,以便使用命名空间获取素数。当需要多次访问结果时,这种方法很有用(但是,运行时间基于未缓存的结果)。所有这四个条目都包含在下面。
以下是两次挑战后的排行榜:
挑战 3 将于 6 月 1 日星期三发布。
挑战 2 要求用两种可能的方法寻找素数。给定一个正整数值n ,你必须寻找:1) n 个素数或 2) 所有小于n的素数。我们从 2) 开始,因为它更容易实现。寻找素数的最简单方法是试除法。以 2 作为候选数和一个空素数集开始,检查当前素数集中是否有任何东西可以整除候选数。如果不是,则将其添加到素数集中并继续下一个候选数。由于素数集中没有任何数字,因此添加 2。检查 3,它不能整除 {2} 中的任何东西,因此将其添加到素数集中。然而,数字 4 可以整除 {2,3} 中的某个数字,因此没有添加。下面显示了一个示例实现。
getPrimes2 = Function( {n}, {primeSet,i,j},
primeSet = [];
For(i=2,i<=n,i++,
For(j=2,j<=N Row(primeSet),j++,
If(Mod(i,j)==0, Break());
);
If(j>N Row(primeSet),primeSet = primeSet |/ i);
);
primeSet;
);
这种方法有几个优点:它很容易理解,代码简洁,并且可以轻松适应素数查找挑战。不幸的是,它的速度非常慢。查找 100,000 以下的素数花了我 19.43 秒。
稍微复杂一点的方法是使用筛子。筛子不是测试数字是否为素数,而是测试数字是否为素数。最古老、最著名的筛子是埃拉托斯特尼筛法。首先假设一切都是素数。以 2 作为第一个素数,找到所有可以整除它的数字,并将它们从素数集中删除。继续下一个仍在素数集中的数字。以下是非优化版本。
getPrimes2 = Function( {n}, {primeSet,i,j},
primeSet = J( n, 1, 1 );
primeSet[1] = 0;
For( i = 2, i <= n, i++,
For( j = i, j <= n / i, j++,
primeSet[j * i] = 0
)
);
Loc( primeSet );
);
最明显的优势是速度。我找到 100,000 以下的所有素数的时间是 0.2811 秒,大约快了 70 倍,找到 1,000,000 的所有素数的时间是 3.4025 秒。代码仍然简洁易懂。一个小缺点是您需要创建一个向量来跟踪非素数(也就是说,您不能通过连接来增加列表——我们稍后会看到这不一定是一个优势)。这使得寻找n 个素数的问题更加困难,因为您不知道确切需要多大的向量。
下面是两个优化版本。第一个版本首先删除所有 2 的倍数,然后仅循环 3 以后的奇数。Sqrt( n ) 用于外循环的上限,因为所有小于 100 的非素数都已被筛选掉。例如,要找到 100 以下的所有素数,只需筛选到 10。任何大于 10 的素数的筛选数字将是该素数乘以小于 10 的数字的乘积,该数字将与其他素数一起被删除(例如,11*5 被 5 的筛选删除,11*9 被 3 的筛选删除,等等)。这些改进将运行时间缩短了一半以上,降至 1.3545。第二段展示了使用 JSL 的矩阵运算进行内循环的优势,1,000,000 的时间为 0.01625,100,000,000 的时间为 2.6694,提高了近 100 倍!
getPrimes2 = Function( {n}, {primeSet,i,j},
primeSet = J( n, 1, 1 );
primeSet[4::n::2] = 0;
primeSet[1] = 0;
For( i = 3, i <= Sqrt( n ), i+=2,
For( j = i, j <= n / i, j++,
primeSet[j * i] = 0
)
);
Loc( primeSet );
);
getPrimes2 = Function( {n}, {primeSet,i,j,upper},
primeSet = J( n, 1, 1 );
primeSet[4::n::2] = 0;
primeSet[1] = 0;
i = 3;
upper = Sqrt(n);
For(i=3,i<=upper,i+=2,
If(primeSet[i] == 1,
primeSet[(i*2)::n::i] = 0;
);
);
Loc( primeSet );
);
您能做得更好吗? Atkin 筛法在理论上更高效,并在Primes Atkins 示例.JSL和Primes Atkins(优化).JSL中实现。前者提供了非优化且更易于理解的实现(1,000,000 需要 2.236 秒)。比 Eratosthenes 的非优化版本更好。但是,优化版本无法击败优化的 Eratosthenes(100,000,000 需要 4.914 秒)。此外,它在编程和理论上都更长更复杂。
那么从数字表中读取呢?虽然这不符合编写算法的精神,但我并没有说你不能这样做。 @Craige_Hales提供了一个使用此方法的解决方案。鉴于问题规范n不超过 1e8,包含压缩值的脚本大小需要超过 400 万个字符的代码、4.6MB 的磁盘空间和大约 5.5 秒的加载时间,无论n的大小如何。比其他解决方案慢。虽然它既不是最快也不是最短的,但我还是包含了 Craige 的脚本,因为他对问题的独特处理方式、压缩的使用以及并行分配。
寻找n 个素数是一个更具挑战性的问题。假设您使用筛子,则需要创建一个包含最多n个素数的值的向量。为此,需要估计该值的上限。应避免低估,因为可能需要重新运行部分或全部算法才能找到未筛选的值。Estimate value of nth prime.jmp显示 6 个不同的估计值。Est1给出了一个在n较小时值被低估的情况。许多估计仅适用于大于特定值的n ,在这些情况下需要小心提供替代估计。Est2 适用于n ≥ 13,而Est3 则针对更大的值进行了改进。不幸的是,好的界限很少是简单的。Est4 用于Primes 快速提交.jsl 。虽然它做得不是很好,并且低估了非常大的n值,但它简单易记。Est5 是使用实际值的非线性得出的。 Est6通过对值进行轻微扰动来改进这是5 ,这样表格中的任何值都不会被低估。不能保证 100 和 1e8 之间不存在低估。其他表格列给出了每种技术的低估/高估。
确定素数可以遵循与生成素数类似的路径。给定一个数字x ,你可以让i = 2。如果Mod( n , i ) == 0 则停止,否则增加i 。继续,直到Mod( n , i ) == 0或i == n 。下面的代码的好处在于,在i = 2 的初始测试之后只检查奇数i ,并在i > = Ceiling(Sqrt( n ))时停止。
isPrime = Function({inVals},{i,j,outVec,nextVal,uBound},
outVec = J(N Row(inVals),1,-1);
For(i=1,i<=N Row(inVals),i++,
nextVal = inVals[i];
If(Mod(nextVal,2) == 0,
outVec[i] == 0;
,//ELSE
j = 3;
uBound = Ceiling(Sqrt(nextVal));
While(j < uBound & Mod(nextVal,j) != 0, j+=2);
If(j >= uBound, outVec[i] = 1,outVec[i] = 0);
)
);
outVec;
);
上述脚本大约需要 1.4 秒来检查表大数表.jmp ( 小的 ) 中的所有 1,000 个八位和九位数字。
为了进一步改进这一点,您可以生成一个素数向量(最多为Ceiling(Sqrt( n )) )并对其进行迭代。
getPrimes = Function( {n}, {primeSet,i,j,upper},
primeSet = J( n, 1, 1 );
primeSet[4::n::2] = 0;
primeSet[1] = 0;
i = 3;
upper = Sqrt(n);
For(i=3,i<=upper,i+=2,
If(primeSet[i] == 1,
primeSet[(i*2)::n::i] = 0;
);
);
Loc( primeSet );
);
isPrime = Function({inVals},{i,j,outVec,primeVec,nextVal},
outVec = J(N Row(inVals),1,-1);
primeVec = getPrimes(Ceiling(Sqrt(Max(inVals))));
For(i=1,i<=N Row(inVals),i++,
nextVal = inVals[i];
j=1;
While(j < N Row(primeVec) & Mod(nextVal,primeVec[j]) != 0, j++);
If(j >= N Row(primeVec), outVec[i]=1,outVec[i]=0);
);
outVec;
);
这将时间缩短了一半以上,降至 0.6 秒以下。
提高速度的两个一般注意事项:在需要多次使用矩阵值的情况下,将值分配给变量(使用变量而不是索引矩阵)会更快。虽然这可能会使代码更长,但可以节省时间。尝试将上述代码中的下一个值替换为其索引矩阵值 ( 入参[i] )。其次,随着矩阵大小的增加, 地方和矩阵连接会变慢。如果可以使用排序矩阵, 本地排序会执行二分搜索,速度会快得多。有关矩阵引用和连接的更多示例,请参阅速度测试矩阵操作.jsl 。
尝试在大数 ( 大的 ) 列中查找素数更加困难。查找前 10 个大约需要 1.9 秒。此时,您必须依靠专门的算法进行改进。Miller -Rabin (MR)是一种常用且相对容易实现的素数测试技术。millerRabin1.jsl中给出了确定性版本。MR 依赖于模幂运算,必须小心谨慎,以免发生溢出错误。MR 本身需要 5.28 秒才能处理所有 1000 个大的数字。您可能会惊讶地发现,MR 在处理小数字时比上面的第一个 isPrime 算法(2.34 秒)慢。这可能是因为随着数字的增加,它的计算效率会更高,这表明可以进行改进:使用上面的 2) 测试某个大小的素数。如果数字是素数,则停止。否则,切换到 MR。此方法可在millerRabin2.jsl中找到。我们进行了一些调整,以找到理想的素数检查数量。这将时间缩短了一半,降至 2.33 秒。