找到大数的组合值

Question

我想找到(n选择r)大整数,我也必须找出那个数字的mod.

long long int choose(int a,int b)
{
    if (b > a)
        return (-1);
    if(b==0 || a==1 || b==a)
        return(1);
    else
    {
        long long int r = ((choose(a-1,b))%10000007+(choose(a-1,b-  1))%10000007)%10000007;
        return r;
    }
}

我正在使用这段代码,但我得到了TLE.如果有其他方法可以做,请告诉我.

Answer 1

我还没有评论的声誉,但我想指出rock321987的答案效果很好:

它是快速和正确的,包括C(62,31)

但无法处理具有适合uint64_t的输出的所有输入.作为证据,请尝试:

C(67,33)= 14,226,520,737,620,288,370 (验证正确性和尺寸)

不幸的是,另一个实现吐出了8,829,174,638,479,413,这是不正确的.还有其他计算nCr的方法不会像这样破坏,但这里的真正问题是没有尝试利用模数.

请注意,p = 10000007是素数,这允许我们利用所有整数具有逆模p的事实,并且该逆是唯一的.此外,我们可以很快发现逆.另一个问题有关于如何在这里做到这一点的答案,我在下面复制了这个问题.

这很方便,因为:

1. x/y mod p == x*(y inverse)mod p; 和
2. xy mod p ==(x mod p)(y mod p)

稍微修改其他代码,并概括问题我们有以下几点:

#include <iostream>
#include <assert.h>

// p MUST be prime and less than 2^63
uint64_t inverseModp(uint64_t a, uint64_t p) {
    assert(p < (1ull << 63));
    assert(a < p);
    assert(a != 0);
    uint64_t ex = p-2, result = 1;
    while (ex > 0) {
        if (ex % 2 == 1) {
            result = (result*a) % p;
        }
        a = (a*a) % p;
        ex /= 2;
    }
    return result;
}

// p MUST be prime
uint32_t nCrModp(uint32_t n, uint32_t r, uint32_t p)
{
    assert(r <= n);
    if (r > n-r) r = n-r;
    if (r == 0) return 1;
    if(n/p - (n-r)/p > r/p) return 0;

    uint64_t result = 1; //intermediary results may overflow 32 bits

    for (uint32_t i = n, x = 1; i > r; --i, ++x) {
        if( i % p != 0) {
            result *= i % p;
            result %= p;
        }
        if( x % p != 0) {
            result *= inverseModp(x % p, p);
            result %= p;
        }
    }
    return result;
}

int main() {
    uint32_t smallPrime = 17;
    uint32_t medNum = 3001;
    uint32_t halfMedNum = medNum >> 1;
    std::cout << nCrModp(medNum, halfMedNum, smallPrime) << std::endl;

    uint32_t bigPrime = 4294967291ul; // 2^32-5 is largest prime < 2^32
    uint32_t bigNum = 1ul << 24;
    uint32_t halfBigNum = bigNum >> 1;
    std::cout << nCrModp(bigNum, halfBigNum, bigPrime) << std::endl;
}

如果您愿意等待,哪个应该产生任何32位输入的结果.为了证明这一点,我已经包括了24位n和最大32位素数的计算.我的微软电脑需要大约13秒才能计算出来.检查wolfram alpha的答案,但要注意它可能超过那里的"标准计算时间".

如果p远小于(nr),其中r <= nr,则仍有改进的余地.例如,我们可以预先计算所有逆变量,而不是按需多次执行.