找到最小 N 的算法，使得 N!能被提升到幂的素数整除答案

【问题标题】：Algorithm to find smallest N such that N! is divisible by a prime raised to a power找到最小 N 的算法，使得 N!能被提升到幂的素数整除
【发布时间】：2011-06-13 01:12:07
【问题描述】：

有没有一种有效的算法来计算最小整数 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 的增加，这种方法很快就会失效，而且似乎没有可以推断出更大值的模式。

于 2011 年 6 月 13 日编辑

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

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

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

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

【问题讨论】：

虽然 Nemo 的回答不是很清楚，但我相信它比二分查找要好。毕竟，它是 O(1)！或者，更准确地说，因为您必须处理数字，所以它是 O(log k)。这个问题是直接可以解决的，所以不需要做一些迭代计算。
最好将答案放入您自己问题的答案中，而不是作为对问题的编辑。
“10 位或更多位”不是“非常大”:-)。我已经编辑了答案以添加 Perl 实现。即使对于几十位数字的 k 似乎也可以正常工作，尽管我不知道答案是否正确。如果你能找到给出错误答案的案例，我想看看。
你应该发布你的编辑作为答案，因为它是一个。

标签： algorithm math numbers

【解决方案1】：

好的，这很有趣。

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

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

您从最高有效位到最低有效位执行此操作。所以首先，找到最大的 j 使得 f(j)

这会生成一个序列 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)

对于 1 的余数，求 1 / f(1) 的商和余数。商为 1，余数为零，所以我们完成了。

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

我在纸上解释了为什么这样有效，但我不确定如何用文本进行交流......请注意，f(i) 回答了这个问题，“p 中有多少个因数在 (p^一世）！”。一旦你找到最大的 i,j 使得 j*f(i) 小于 k，并意识到你真正在做的是找到小于 N 的最大 j*p^i，其余的就掉出来了.例如，在我们的 p=3 示例中，1-9 的乘积贡献了 4 个 p，10-18 的乘积贡献了 4 个，21 贡献了一个。前两个只是 p 的倍数^ 2； f(2) = 4 告诉我们 p^2 的每个倍数对乘积贡献了 4 个 p。

[更新]

代码总是有助于澄清。将以下 perl 脚本另存为 foo.pl 并以 foo.pl <p> <k> 运行。请注意，** 是 Perl 的求幂运算符，bdiv 计算 BigInts（无限精度整数）的商和余数，use bigint 告诉 Perl 在任何地方都使用 BigInts。

#!/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;

【讨论】：

这感觉像是 p-adic 范数可以帮助解释的东西

【解决方案2】：

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

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

【讨论】：

【解决方案3】：

您为什么不尝试使用您提到的第二个公式来二分查找答案？

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

【讨论】：

【解决方案4】：

考虑

我 = (pⁿ)！

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

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

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

(n² +n)/2 >= k

如果我记得正确的二次方程给我们的话

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

（P.S. 有人知道如何在 markdown 中显示激进符号吗？）

编辑：

这是错误的。让我看看我能不能挽救它......

【讨论】：