只使用Int64中间体的Int32模取M的平方锥体数

本文关键字：方锥体 模取 Int32 Int64 中间体 | 更新日期: 2023-09-27 18:11:26

计算n到10^9(和素数M)的平方金字塔数n (n + 1) (2 n + 1) / 6 mod 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);
}

正如你所看到的，这真的很尴尬。有没有一种更优雅/更有效的方法来执行这个计算，而不诉诸于BigInteger或Decimal，也许以某种方式使用中间约简模3m ?

背景:这个问题是在HackerEarth解决井字棋练习问题时出现的。基于我笨拙的hack的提交被接受了，但我不满意这个不成熟的解决方案。这就是这些练习题的全部意义所在，不是吗?如果我接受任何基于已有知识的半生不熟的解决方案，那我就什么也学不到。因此，我一直致力于改进解决方案，直到它们达到简单和优雅的状态。

只使用Int64中间体的Int32模取M的平方锥体数

我对简化模3m的直觉被淘汰了——在测试表明它有效之后，我花了一些时间在数学上把它固定下来。

关键是中国剩余定理，它有效地保证了对素数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);
}

这使用了unsigned算术，这在这里是合适的，也更有效，因为有符号算术通常需要编译器额外的努力(读:发出额外的指令)来实现有符号算术语义。

基准测试显示，这比我的问题中的原始代码快两倍多，但仅在旧框架版本(高达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慢得多。

环境。版本和时间对于3.5以下的框架版本是完全相同的，所以我猜它们都使用相同的环境/编译器。

现在是微软新的和"改进的"编译器随框架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

环境。版本和时间与框架版本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模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相同。

相关文章：