0

计算 n 到 10^9（和质数 M）的值的方形金字塔数 n (n + 1) (2 n + 1) / 6 mod M会带来一些挑战，因为模减少之前的中间结果可能超过 10^27，因此对于 64 位整数来说可能太大了。

在乘法之前减少因子模 M 会产生除以 6 的问题，因为在减少模 M 之后执行该除法显然会产生无意义的结果。

我正在使用基于以下事实的解决方法：n (n + 1)对于任何 n 都必须是偶数，并且n (n + 1)必须(2 n + 1)可以被 3 整除：

const int M = 1000000007;

static int modular_square_pyramidal_number (int n)
{
    var a = (Int64)n * (n + 1) / 2;
    var b = 2 * n + 1;
    var q = a / 3;
    var p = q * 3 == a ? (q % M) * b : (a % M) * (b / 3);

    return (int)(p % M);
}

正如你所看到的，这真的很尴尬。是否有一种更优雅/更有效的方式来执行此计算而不使用 BigInteger 或 Decimal，也许以某种方式使用中间约简模 3 M？

背景：这个问题是在 HackerEarth解决井字游戏练习题时出现的。基于我笨拙的 hack 的提交被接受了，但我对这个半生不熟的解决方案不满意。这就是这些练习题的全部意义，不是吗：如果我接受任何基于预先存在的知识的半生不熟的解决方案，我不会学到任何东西，这些知识会被机器人法官刮掉。因此，我一直致力于改进解决方案，直到它们达到简单和优雅的状态......

4

1 回答 1

1

我对减少模 3 M 的直觉被淘汰了——在测试表明它有效后，我花了一段时间才用数学方法确定了这个东西。

关键是中国剩余定理，它有效地保证了互质 p 和 q

(x / q) mod p = ((x mod pq) / q) mod p

让我们对要计算的公式进行与我的问题相同的拆分：

n (n + 1) (2 n + 1) / 6 mod M = a b / 3 mod M

a = n (n + 1) / 2
b = 2 n + 1

a 或 b 必须能被 3 整除，但不知道是哪一个，并且a * b可能太大而无法放入 64 位整数（大约 90 位，给定 n ≤ 1e9 的原始约束）。

但是，对于M = 1000000007（即通常的 1e9 + 7），该项3 * M仅需要 32 位，对于约简a模 3 M 也是如此。由于b已经适合 31 位，这意味着可以使用 64 位算术计算乘积：

((a mod 3 M) * b) / 3 mod M

更改代码：

static int v1 (int i)
{
    var n = (uint)i;
    var a = ((UInt64)n * (n + 1) >> 1) % (M * 3U);
    var b = 2 * n + 1;

    return (int)((a * b / 3) % M);
}

这使用无符号算术，这在此处是合适的并且更有效，因为有符号算术通常需要编译器额外的努力（阅读：额外指令的发射）以实现有符号算术语义。

基准测试显示这比我的问题中的原始代码快两倍以上 - 但仅限于旧框架版本（最高 3.5）。从 4.0 版开始，JIT 编译器不再将无符号除以常量转换为乘法 + 移位。除法指令往往比乘法指令至少慢一个数量级，因此代码变得比使用较新编译器的系统上的原始代码慢很多。

在这样的系统上，最好顺其自然并使用低效但政治上正确的有符号整数：

static int v2 (int n)
{
    var a = ((Int64)n * (n + 1) >> 1) % (M * 3L);
    var b = 2 * n + 1;

    return (int)((a * b / 3) % M);
}

以下是针对框架版本 2.0 的老化 Haswell 笔记本电脑上 1000000 次调用的基准：

IntPtr.Size = 8, Environment.Version = 2.0.50727.8009
bench 1000000:    8,407 v0    3,413 v1    4,653 v2
bench 1000000:    8,017 v0    3,179 v1    5,038 v2
bench 1000000:    8,641 v0    3,114 v1    4,801 v2

时间以毫秒为单位，v0 代表我的问题中的原始代码。很容易看出有符号语义的开销如何使 v2 明显慢于在内部使用无符号算术的 v1。

对于最高 3.5 的框架版本，Environment.Version 和时序完全相同，所以我猜它们都使用相同的环境/编译器。

现在是微软新的和“改进的”框架 4.0 和更新版本的编译器的时间安排：

IntPtr.Size = 8, Environment.Version = 4.0.30319.42000
bench 1000000:    9,518 v0   20,479 v1    5,687 v2
bench 1000000:    9,225 v0   20,251 v1    5,540 v2
bench 1000000:    9,133 v0   20,333 v1    5,389 v2

Environment.Version 和时间对于框架版本 4.0 到 4.6.1 完全相同。

POST SCRIPTUM - 使用模乘逆

另一种解决方案是使用除数的模乘逆。在本例中，这是可行的，因为已知最终产品可以被除数（即 3）整除；如果不是，那么结果将非常不准确。示例（333333336 是 3 模 1000000007 的乘法逆元）：

7 * 333333336 % 1000000007 = 333333338  // 7 mod 3 != 0
8 * 333333336 % 1000000007 = 666666674  // 8 mod 3 != 0
9 * 333333336 % 1000000007 =         1  // 9 mod 3 == 0

该主题存在的理由是整数除法可能有损失，因为它会丢弃余数（如果有的话），因此如果错误的因子除以 3，则金字塔平方计算的结果将是错误的。

模除法 - 即与乘法逆相乘 - 不是有损的，因此哪个因子与逆相乘无关紧要。这在刚刚显示的示例中很容易看出，其中 7 和 8 的古怪残差有效地编码了小数余数，并且将它们相加 - 对应于计算 7/3 + 8/3 - 得到 1000000012 等于 5 mod 1000000007 就像预期的。

因此，问题的关键是最终产品可以被除数整除，但“除法”（与逆相乘）发生的时间和地点并不重要。生成的代码比 v1 效率稍低，但与 v2 大致相当，因为它需要在与逆相乘之后额外减少模 M。但是，无论如何我都会展示它，因为这种方法有时会派上用场：

static int v3 (int n)
{
    var a = n * (n + 1L) % M;
    var b = (2 * n + 1L) * 166666668 % M;

    return (int)(a * b % M);
}

注意：我放弃了右移并将除数 2 合并到逆中，因为单独除以 2 在这里不再有任何用途。时间与 v2 相同。

于 2016-10-04T17:02:54.463 回答