仅使用 Int64 中间体的 Int32 模数 M 的方形金字塔数
Square pyramidal number of Int32 modulo some M using only Int64 intermediates
为 n 的值计算 square pyramidal number n (n + 1) (2 n + 1) / 6 mod M 到 10^9（和素数 M）会带来一些挑战，因为模数减少之前的中间结果可能超过 10^27 并且因此对于 64 位整数来说可能太大了。
在乘法之前以 M 为模约化因数会导致除以 6 出现问题，因为在以 M 为模约化之后执行该除法显然会产生无意义的结果。
A 我正在使用基于以下事实的解决方法：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);
}
如你所见，这真的很尴尬。是否有更多 elegant/efficient 方法来执行此计算而不求助于 BigInteger 或 Decimal，也许以某种方式使用中间缩减模 3 M？
背景：问题出现在解决 HackerEarth 的 Tic Tac Toe 练习题中。基于我笨拙的 hack 的提交被接受了，但我对这个半生不熟的解决方案不满意。这就是这些练习题的全部意义，不是吗：如果我接受任何基于机器人裁判刮擦的现有知识的半生不熟的解决方案，我就不会学到任何东西。因此，我一直致力于改进解决方案，直到它们达到简单和优雅的状态...
我对缩减的直觉 modulo 3 M 摇摇欲坠 - 在测试表明它有效后，我花了一段时间从数学上确定了它。
关键是 Chinese Remainder Theorem 有效保证互质 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 减少 modulo 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);
}
这使用无符号算术，这在这里是合适的并且也更有效，因为带符号算术通常需要编译器额外的努力（阅读：附加指令的发射）以实现带符号算术语义。
Benchmark 显示这比我的问题的原始代码快两倍以上 - 但仅在旧框架版本（最高 3.5）下。从 4.0 版开始，JIT 编译器不再将 unsigned 常量除法转换为乘法 + 移位。除法指令往往比乘法指令至少慢一个数量级，因此代码变得比具有较新编译器的系统上的原始代码慢很多。
在这样的系统上，最好顺其自然，使用效率低下但政治上正确的带符号整数：
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);
}
以下是在我老化的 Haswell 笔记本电脑上针对框架版本 2.0 调用 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 慢。
Environment.Version 框架版本高达 3.5 的时间完全相同，所以我猜他们都使用相同的 environment/compiler.
现在 Microsoft 新的和“改进的”编译器与框架 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 - 使用mod元乘法逆
另一种解决方案是使用除数的 modular multiplicative inverse。在本例中，这是可行的，因为已知最终产品可以被除数整除（即 3）；如果不是，那么结果将非常不准确。示例（333333336 是 3 modulo 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 大致相当，因为它在与逆相乘后需要额外减少 modulo 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 相同。
为 n 的值计算 square pyramidal number n (n + 1) (2 n + 1) / 6 mod M 到 10^9（和素数 M）会带来一些挑战，因为模数减少之前的中间结果可能超过 10^27 并且因此对于 64 位整数来说可能太大了。
在乘法之前以 M 为模约化因数会导致除以 6 出现问题，因为在以 M 为模约化之后执行该除法显然会产生无意义的结果。
A 我正在使用基于以下事实的解决方法：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);
}
如你所见，这真的很尴尬。是否有更多 elegant/efficient 方法来执行此计算而不求助于 BigInteger 或 Decimal，也许以某种方式使用中间缩减模 3 M？
背景：问题出现在解决 HackerEarth 的 Tic Tac Toe 练习题中。基于我笨拙的 hack 的提交被接受了，但我对这个半生不熟的解决方案不满意。这就是这些练习题的全部意义，不是吗：如果我接受任何基于机器人裁判刮擦的现有知识的半生不熟的解决方案，我就不会学到任何东西。因此，我一直致力于改进解决方案，直到它们达到简单和优雅的状态...
我对缩减的直觉 modulo 3 M 摇摇欲坠 - 在测试表明它有效后，我花了一段时间从数学上确定了它。
关键是 Chinese Remainder Theorem 有效保证互质 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 减少 modulo 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);
}
这使用无符号算术，这在这里是合适的并且也更有效，因为带符号算术通常需要编译器额外的努力（阅读：附加指令的发射）以实现带符号算术语义。
Benchmark 显示这比我的问题的原始代码快两倍以上 - 但仅在旧框架版本（最高 3.5）下。从 4.0 版开始，JIT 编译器不再将 unsigned 常量除法转换为乘法 + 移位。除法指令往往比乘法指令至少慢一个数量级，因此代码变得比具有较新编译器的系统上的原始代码慢很多。
在这样的系统上，最好顺其自然，使用效率低下但政治上正确的带符号整数：
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);
}
以下是在我老化的 Haswell 笔记本电脑上针对框架版本 2.0 调用 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 慢。
Environment.Version 框架版本高达 3.5 的时间完全相同，所以我猜他们都使用相同的 environment/compiler.
现在 Microsoft 新的和“改进的”编译器与框架 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 - 使用mod元乘法逆
另一种解决方案是使用除数的 modular multiplicative inverse。在本例中，这是可行的，因为已知最终产品可以被除数整除（即 3）；如果不是，那么结果将非常不准确。示例（333333336 是 3 modulo 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 大致相当，因为它在与逆相乘后需要额外减少 modulo 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 相同。