以下是问题描述:
让我们c[n]
成为加泰罗尼亚数n
并且p
是一个大素数,例如。1000000007
我需要计算范围c[n] % p
从n
{1,2,3,...,1000}
我遇到的问题是,在 32 位机器上,当你为这么大的整数计算加泰罗尼亚数时,你会溢出。我熟悉模运算。还
(a.b) % p = ((a % p)(b % p)) % p
这个公式可以帮助我分别解决分子中的溢出问题,但我不知道如何处理分母。
对于 1000000007 的模数,仅使用 32 位整数避免溢出是很麻烦的。但是任何体面的 C 实现都提供 64 位整数(任何体面的 C++ 实现也提供),所以这不是必需的。
然后处理分母,一种可能性是,正如 KerrekSB 在他的评论中所说,计算分母模数的模逆p = 1000000007
。您可以使用扩展的欧几里得算法或等效的 的连分数展开来计算模逆k/p
。然后,不是在计算中除以k
,而是乘以它的模逆。
另一种选择是对加泰罗尼亚数使用 Segner 的递推关系,它给出了没有除法的计算:
C(0) = 1
n
C(n+1) = ∑ C(i)*C(n-i)
0
由于您只需要 的加泰罗尼亚数字C(k)
,k <= 1000
您可以预先计算它们,或在程序启动时快速计算它们并将它们存储在查找表中。
如果与预期相反,没有可用的 64 位整数类型,您可以通过将因子拆分为低 16 位和高 16 位来计算模积,
a = a1 + (a2 << 16) // 0 <= a1, a2 < (1 << 16)
b = b1 + (b2 << 16) // 0 <= b1, b2 < (1 << 16)
a*b = a1*b1 + (a1*b2 << 16) + (a2*b1 << 16) + (a2*b2 << 32)
计算a*b (mod m)
,m <= (1 << 31)
减少四个产品中的每一个模m
,
p1 = (a1*b1) % m;
p2 = (a1*b2) % m;
p3 = (a2*b1) % m;
p4 = (a2*b2) % m;
合并这些变化的最简单方法是
for(i = 0; i < 16; ++i) {
p2 *= 2;
if (p2 >= m) p2 -= m;
}
p3
和 32 次迭代相同p4
。然后
s = p1+p2;
if (s >= m) s -= m;
s += p3;
if (s >= m) s -= m;
s += p4;
if (s >= m) s -= m;
return s;
这种方式不是很快,但是对于这里需要的少数乘法,它已经足够快了。应通过减少班次数来获得小幅加速;首先计算(p4 << 16) % m
,
for(i = 0; i < 16; ++i) {
p4 *= 2;
if (p4 >= m) p4 -= m;
}
那么所有的p2
,p3
的当前值p4
需要乘以 2 16模m
,
p4 += p3;
if (p4 >= m) p4 -= m;
p4 += p2;
if (p4 >= m) p4 -= m;
for(i = 0; i < 16; ++i) {
p4 *= 2;
if (p4 >= m) p4 -= m;
}
s = p4+p1;
if (s >= m) s -= m;
return s;
如果您使用动态编程存储结果并且在填充查找表时,您可以在每个步骤中使用 MODULO 除法。它将处理 1000 个加泰罗尼亚人的溢出问题,并且比 BigDecimal/BigInteger 更快。
public class Catalan {
private static long [] catalan= new long[1001];
private static final int MOD=1000000007;
public static void main(String[] args) {
precalc();
for (int i=1;i<=1000;i++){
System.out.println("Catalan number for "+i+" is: "+catalan[i]);
}
}
private static void precalc(){
for (int i=0;i<=1000;i++){
if (i==0 || i==1){
catalan[i]=1;
}
else{
long sum =0;long left, right;
for (int k=1;k<=i;k++){
left = catalan[k-1] % MOD;
right= catalan[i-k] % MOD;
sum =(sum+ (left * right)%MOD)%MOD;
}
catalan[i]=sum;
}
}
}
}
使用大整数库怎么样?尝试谷歌搜索它...
#include <stdio.h>
#include <stdlib.h>
/*
C(n) = (2n)!/(n+1)!n!
= (2n)(2n-1)(2n-2)..(n+2)/n!
*/
int p = 1000000007;
int gcd(int x, int y){
while(y!=0){
int wk = x % y;
x = y;
y = wk;
}
return x;
}
int catalanMod(n){
long long c = 1LL;
int i;
int *list,*wk;
//make array [(2n),(2n-1),(2n-2)..(n+2)]
wk = list = (int*)malloc(sizeof(int)*(n-1));
for(i=n+2;i<=2*n;++i){
*wk++ = i;
}
wk=list;
//[(2n),(2n-1),(2n-2)..(n+2)] / [1,2,3,..n]
//E.g C(10)=[13,17,19,4]
for(i=2;i<=n;++i){
int j,k,w;
for(w=i,j=0;j<n-1;++j){
while(1!=(k = gcd(wk[j], w))){
wk[j] /= k;
w /= k;
}
if(w == 1) break;
}
}
wk=list;
//Multiplication and modulo reduce
for(i=0;i<n-1;++i){
if(wk[i]==1)continue;
c = c * wk[i] % p;
}
free(list);
return c;
}
简单地说,使用属性,(a * b) % mod = (a % mod) * (b % mod)