找到最小 N 使得 N! 的算法能被素数的幂整除

Question

是否有一种有效的算法来计算最小整数 N 使得 N!可被 p^k 整除，其中 p 是一个相对较小的素数，k 是一个非常大的整数。换句话说，

factorial(N) mod p^k == 0

如果给定 N 和 p，我想找出 p 能除以 N 的次数！，我会使用众所周知的公式，

k = Sum(floor(N/p^i) for i=1,2,...

我已经对 k 的小值进行了强力搜索，但这种方法非常失败随着 k 的增加，速度很快，而且似乎没有一种模式可以推断出更大的值。

编辑于 6/13/2011

根据 Fiver 和 Hammar 提出的建议，我使用准二进制搜索来解决问题，但并不完全按照他们建议的方式。使用上面第二个公式的截断版本，我将 N 的上限计算为 k 和 p 的乘积（仅使用第一项）。我使用 1 作为下限。使用经典的二分搜索算法，我计算了这两个值之间的中点，并使用该中点值作为第二个公式中的 N 来计算 k 的值，这次使用了所有项。

如果计算出的 k 太小，我会调整下限并重复。太大了，我首先测试一下在 midpoint-1 计算的 k 是否小于所需的 k。如果是这样，则返回中点作为最接近的 N。否则，我调整高点并重复。

如果计算出的 k 相等，我测试 midpoint-1 处的值是否等于 midpoint 处的值。如果是这样，我将高点调整为中点并重复。如果 midpoint-1 小于所需的 k，则返回中点作为所需的答案。

即使 k 值非常大（10 位或更多数字），这种方法也能以 O(n log(n)) 的速度运行。

Answer 1

好吧，这很有趣。

定义 f(i) = (p^i - 1) / (p - 1)

将 k 写成一种有趣的“基数”，其中位置 i 的值就是这个 f(i)。

您可以从最高有效数字到最低有效数字执行此操作。首先，找到满足 f(j) <= k 的最大 j。然后计算 k / f(j) 的商和余数。将商存储为 q_j，将余数存储为 r。现在计算 r / f(j-1) 的商和余数。将商存储为 q_{j-1}，将余数再次存储为 r。现在计算 r / f(j-2) 的商和余数。等等。

这会生成序列 q_j, q_{j-1}, q_{j-2}, ..., q_1。 （请注意，序列以 1 结尾，而不是 0。）然后计算 q_j*p^j + q_{j-1}*p^(j-1) + ... q_1*p。这就是您的 N。

示例：k = 9，p = 3。因此 f(i) = (3^i - 1) / 2。f(1) = 1，f(2) = 4，f(3) = 13因此，f(j) <= 9 时最大的 j 是 i = 2，且 f(2) = 4。取 9 / 4 的商和余数。这是 2 的商（即2 位置上的数字）和 1 的余数。

对于 1 的余数，找到 1 / f(1) 的商和余数。商为 1，余数为零，这样我们就完成了。

所以q_2 = 2，q_1 = 1。2*3^2 + 1*3^1 = 21，这是正确的N。

我在纸上有一个解释为什么这样有效，但我不知道如何用它来表达它文本...请注意，f(i) 回答了问题“(p^i) 中有 p 有多少个因子！”。一旦你找到了最大的 i,j 使得 j*f(i) 小于 k，并且意识到你真正在做的是找到小于 N 的最大的 j*p^i，剩下的就变得不那么重要了。例如，在我们的 p=3 示例中，我们得到 1-9 的乘积贡献的 4 个 p，10-18 的乘积贡献的另外 4 个 p，以及 21 的乘积贡献的另外一个。前两个只是 p^ 的倍数2； f(2) = 4 告诉我们，p^2 的每个倍数都会为乘积贡献 4 个以上的 p。

[更新]

代码总是有助于澄清。将以下 perl 脚本保存为 foo.pl 并以 foo.pl

运行它。请注意，** 是 Perl 的幂运算符，bdiv 计算 BigInts（无限精度整数）的商和余数，use bigint 告诉 Perl到处使用 BigInt。

#!/usr/bin/env perl

use warnings;
use strict;
use bigint;

@ARGV == 2
    or die "Usage: $0 <p> <k>\n";

my ($p, $k) = map { Math::BigInt->new($_) } @ARGV;

sub f {
    my $i = shift;
    return ($p ** $i - 1) / ($p - 1);
}

my $j = 0;
while (f($j) <= $k) {
    $j++;
}
$j--;

my $N = 0;

my $r = $k;
while ($r > 0) {
    my $val = f($j);
    my ($q, $new_r) = $r->bdiv($val);
    $N += $q * ($p ** $j);
    $r = $new_r;
    $j--;
}

print "Result: $N\n";

exit 0;

Answer 2

使用您提到的公式，给定固定 p 和 N = 1,2... 的 k 值序列是非递减的。这意味着您可以使用二分搜索的变体来查找给定所需 k 的 N。

• 从N = 1开始，计算k。
• 双倍 N 直到 k 大于或等于您想要的 k 以获得上限。
• 对剩余间隔进行二分搜索以找到您的 k。

Answer 3

为什么不尝试使用您提到的第二个公式进行二分搜索来寻找答案？

您只需要考虑 N 的值，其中 p 除以 N，因为如果不是，则为 N！和（N-1）！除以 p 的同次幂，因此 N 不可能是最小的。

Answer 4

考虑

I = (pⁿ)！

并忽略 p 以外的素因数。结果看起来像

I = pⁿ * p^n-1 * p^n-2 * ... * p * 1
I = p^{n + (n-1) + (n-2) + ... 2 + 1}
I = p^{(n2 +n)/2}

所以我们试图找到最小的 n 使得

(n² +n) /2 >= k

如果我还记得二次方程的话，就会得到

N = pⁿ，其中 n >= (sqrt(1+8k) -1)/2

（PS有谁知道如何在markdown中显示偏旁符号？）

编辑：

这是错误的。让我看看是否可以挽救它......