如何计算均匀分布的二项式系数的总和答案

【问题标题】：How to compute sum of evenly spaced binomial coefficients如何计算均匀分布的二项式系数的总和
【发布时间】：2014-04-10 21:57:26
【问题描述】：

如何以 M 为模求均匀分布的二项式系数的总和？
IE。 (ⁿC_a + ⁿC_a+r + ⁿC_a+2r + ⁿC_a+3r + ... + ⁿC_{a+kr子>) % M = ?

给定：0 5, r}

我的第一次尝试是：

int res = 0;
int mod=1000000009;
for (int k = 0; a + r*k <= n; k++) {
    res = (res + mod_nCr(n, a+r*k, mod)) % mod;
}

但这不是有效的。所以看完here 而这个paper我发现上面的总和相当于：
求和[ω^-ja * (1 + ω^j)ⁿ / r]，对于 0 i2π/r 是一个原始的 r^th 单位根。
在 Order(r) 中找到这个总和的代码是什么？

编辑： n 可以达到 10⁵，r 可以达到 100。

原问题出处：https://www.codechef.com/APRIL14/problems/ANUCBC
比赛问题的社论：https://discuss.codechef.com/t/anucbc-editorial/5113
6年后重温这篇文章后，我不记得我是如何将原始问题陈述转换为我的版本的，但是，我分享了原始解决方案的链接，以防有人想看看正确的解决方法.

【问题讨论】：

nCa 你的意思是C(k,n)，来自n的k个元素的组合吗？
是的，更具体地说（n 选择 k MOD 1000000009）
请说明问题的根源并提供链接。

标签： c++ math combinations binomial-coefficients

【解决方案1】：

二项式系数是多项式 (1+x)^n 的系数。 x^a、x^(a+r)等的系数之和就是x^a在多项式环mod x^r-1中的(1+x)^n的系数。多项式 mod x^r-1 可以由长度为 r 的系数数组指定。您可以通过repeated squaring 计算 (1+x)^n mod (x^r-1, M)，在每一步减少 mod x^r-1 和 mod M。这需要大约 log_2(n)r^2 步和 O(r) 空间以及天真的乘法。如果您使用快速傅里叶变换对多项式进行乘法或取幂，则速度会更快。

例如，假设 n=20 且 r=5。

(1+x)    = {1,1,0,0,0}
(1+x)^2  = {1,2,1,0,0}
(1+x)^4  = {1,4,6,4,1}
(1+x)^8  = {1,8,28,56,70,56,28,8,1} 
           {1+56,8+28,28+8,56+1,70}
           {57,36,36,57,70}
(1+x)^16 = {3249,4104,5400,9090,13380,9144,8289,7980,4900}
           {3249+9144,4104+8289,5400+7980,9090+4900,13380}
           {12393,12393,13380,13990,13380}

(1+x)^20 = (1+x)^16 (1+x)^4
         = {12393,12393,13380,13990,13380}*{1,4,6,4,1}
           {12393,61965,137310,191440,211585,203373,149620,67510,13380}
           {215766,211585,204820,204820,211585}

这会告诉您 a 的 5 个可能值的总和。例如，对于 a=1，211585 = 20c1+20c6+20c11+20c16 = 20+38760+167960+4845。

【讨论】：

【解决方案2】：

类似的东西，但你必须检查a、n 和r，因为我只是放了任何东西而没有考虑条件：

#include <complex>
#include <cmath>
#include <iostream>

using namespace std;

int main( void )
{
    const int r = 10;
    const int a = 2;
    const int n = 4;

    complex<double> i(0.,1.), res(0., 0.), w;

    for( int j(0); j<r; ++j )
    {
        w = exp( i * 2. * M_PI / (double)r );

        res += pow( w, -j * a ) * pow( 1. + pow( w, j ), n ) / (double)r;
    }

    return 0;

}

【讨论】：

这适用于较小的 n 值。此外，我对输出感到惊讶，因为结果绝对应该是一个适当的整数，但是使用这个公式，我们得到一个复数，其实部等于所需的总和和一些虚部！
我会计算和使用wj = exp( i * 2 * j * M_PI / r )。
是的，我忘了说结果是 Re(res)
用long long 试试你的计算，因为int 的上限是4E9 类似的东西
比 long long 更好的 int64_t 或 uint64_t。限制分别保证为 9223372036854775807 和 18446744073709551615

【解决方案3】：

mod 操作成本高，尽量避免使用它

uint64_t res = 0;
int mod=1000000009;
for (int k = 0; a + r*k <= n; k++) {
    res += mod_nCr(n, a+r*k, mod);
    if(res > mod)
        res %= mod;
}

我没有测试这段代码

【讨论】：

此方法效率不高。我需要帮助来实现该公式，因为它涉及复数和模数。

【解决方案4】：

我不知道你是否在这个问题中得到了什么，但实现这个公式的关键是要真正弄清楚 w^i 是独立的，因此可以形成一个环。简而言之，您应该考虑实现 (1+x)^n%(x^r-1) 或在环 Z[x]/(x^r-1) 中找出 (1+x)^n 如果感到困惑，我现在会给你一个简单的实现。

制作一个大小为 r 的向量。 O(r) 空间 + O(r) 时间
在 O(r) 空间 +O(r) 时间处用零初始化这个向量
使该向量的前两个元素为 1 O(1)
使用快速求幂法计算 (x+1)^n。每个乘法需要 O(r^2) 并且有 log n 乘法因此 O(r^2 log(n) )
返回向量的第一个元素。O(1) 复杂 O(r^2 log(n) ) 时间和 O(r) 空间。可以使用傅立叶变换将此 r^2 简化为 r log(r)。乘法是怎么做的，这是正则多项式乘法，幂为mod

向量 p1(r,0); 向量 p2(r,0); p1[0]=p1[1]=1； p2[0]=p2[1]=1；现在我们要做乘法向量 res(r,0); for(int i=0;i

【讨论】：