hdu5490 Simple Matrix 组合数与数列多次求和

原题目：hdu 5490 Simple Matrix

题意：

给定一个表格，一个等比数列作为第0行，一个等差数列作为第0列，剩余格子的值F(n,m)=F(n-1,m)+F(n,m-1)，求第n行第m列的值。

每个case均是不同表格，每个不同表格只求一个F(n,m)，结果mod 1000000007

分析：

O(n+m)

O(n+m)是不够快的，只是当时只能想到这步。下一段有O(n)的解法

画出表格示意图 (3*3)

0	b1	b2	b3
a1	a1+b1	a1+b1+b2	a1+b1+b2+b3
a2	a1+a2+b1	2a1+a2+2b1+b2	3a1+a2+3b1+2b2+b3
a3	a1+a2+a3+b1	3a1+2a2+a3+3b1+b2	6a1+3a2+a3+6b1+3b2+b3

一看很复杂，其实可以将等差和等比两个数列分开来算，再求和

表a

0	0	0	0
a1	a1	a1	a1
a2	a1+a2	2a1+a2	3a1+a2
a3	a1+a2+a3+b1	3a1+2a2+a3	6a1+3a2+a3

表b

0	b1	b2	b3
0	b1	b1+b2	b1+b2+b3
0	b1	2b1+b2	3b1+2b2+b3
0	b1	3b1+b2	6b1+3b2+b3

表a加上表b就是我们要的答案了。

表a怎么求？

观察表a，a1的系数，下表只是把a1的系数列出来。会发现是一个起点在(1,1)的杨辉三角形。

0	0	0	0
1	1	1	1
0	1	2	3
0	1	3	6

所以F(n,m)中a1的系数是C(n-1+m-1,n-1) 。 注：C(n,k)是组合数,C(n,k)=n!/(k!(n-k)!)

同理，对a2取系数矩阵，是起点在(2,1)的杨辉三角形，系数是C*n-2+m-1,n-2)

对所有an，bn做同样的工作，不易发现

F(n,m)等于 sum of i=1->n : C(n-i+m-1,n-i)*ai 加上 sum of j=1->m : C(n-1+m-j,m-j)*bj

这样理论的时间复杂度是O(n+m)，但是当时太天真，没发现m的范围是1-2^23，只有n才是1-10000，所以这个思路怎么写都是超时的。而且组合数太大了。1000

O(n)

表a的优化，把a2到an拆写成a1+(n-1)d的形式：

0	0	0	0
a1	a1	a1	a1
a1+d	2a1+d	3a1+d	4a1+d
a1+2d	3a1+3d	6a1+4d	10a1+5d

再对a1跟d分别取系数矩阵，得两个杨辉三角形

所以答案中关于等差数列的部分是 C(n+m-1,n-1)*a1+C(n+m-1,n-2)*d

能够很快得出等差数列部分还是不行，等比数列一共有m项，而m<=2^31

表b的优化

0	b1	b2	b3
0	b1	b1+b2	b1+b2+b3
0	b1	2b1+b2	3b1+2b2+b3
0	b1	3b1+b2	6b1+3b2+b3

重新观察表b

当n=0时，F(n,m)=bm

当n>=1时，F(n,m)= sum of i=1->m: F(n-1,i)

用一句话表达就是，每行最右边的值等于上一行所有数的求和。

于是进入了今天的重点，等比数列的求和，再求和，再求和。。一共求了n次。于是F(n,m)就变成了数列bm的前m项和的前m项和的前m项和。。一共来了n次。

当n=0时，F(0,m)=bm

当n=1时，F(1,m)=b1(1-q^m)/(1-q)=
(q*bm-b1)/(q-1)= (q*F(0,m)-b1*1)/(q-1)

当n=2时，F(2,m)=F(1,1)+F(1,2)+F(1,3)+...F(1,m)= (q*F(1,m)-b1*m)/(q-1)

当n=3时，F(3,m)=F(2,1)+F(2,2)+F(2,3)+...F(2,m)= (q*F(2,m)-b1*m*(m+1)/2)/(q-1)

后面不写了，已经的到规律了。红色公式是相同部分，第一个不同部分分别代入自身的前一项。

第二个不同部分就是：一个长度为m，初始值都为1的数列的求和，再求和，再求和。。给一个容易理解点的表格

列号	1	2	3	...	m
初始数列	1	1	1	...	1
第1次求和	1	2	3	...	m
第2次求和	1	3	6	...	m*(m+1)/2

当n依次取1,2,3时，F(n,m)中的b1的系数依次为表格中第n行第m列的数的值，而这个表格第n行第m列的值就是第n-1行前m格数的求和。

幸运的是，观察左边的数字，这又是一个杨辉三角形，所以这个表格第n行第m列的数等于C(m+n-2,n-1)

所以得出一个足够快的递推式： 当n=0时，F(0,m)=bm。 当n>=1时，F(n,m)=(q*F(n-1,m)-b1*C(m+n-2,n-1))/(q-1)

有了这个递推式，就能在O(n)内求出F(n,m)关于等比部分的答案。

最后一点注意的地方

for i=1->n 的过程中，每次计算F(i,m)=(q*F(i-1,m)-b1*C(m+i-2,i-1))/(q-1)时都重新算了一次组合数，会造成复杂度退化，变成O(n^2)

解决方法是预先求出一个数组BC[i]代替C(m+i-2,i-1),
递推式变成F(i,m)=(q*F(i-1,m)-b1*BC[i])/(q-1)

求BC[i]的方法，

BC[1]=C(m-1,0)=1

for i=2->n: BC[i]=BC[i-1]*(m+i-2)/(i-1)

贴上代码，代码一些数组名用了拼音，一些模板出自ACdreamer的博客，代码后附上来源链接。

----

#include <iostream>
#include <cstdio>
using namespace std;
typedef long long int LL;
const LL MOD=1000000007;

int b1,q,a1,d,n,m;
int nJieChen[10001],nNi[10001],BC[10001];

inline LL multi(LL a,LL b)
{
return (a%MOD)*(b%MOD)%MOD;
}

LL quick_mod(LL a, LL b)
{
LL ans=1; a%=MOD;
while(b)
{
if(b&1)
{
ans=multi(ans,a);
b--;
}
b>>=1;
a=multi(a,a);
}
return ans;
}

int C(int n, int m)
{
if(m>n) return 0;
LL ans=1;
for(int i=1;i<=m;i++)
ans=multi(ans,n-i+1);
ans=multi(ans,quick_mod(nJieChen[m],MOD-2));
return ans;
}

LL Lucas(LL n, LL m)
{
if(m==0) return 1;
return C(n%MOD,m%MOD)*Lucas(n/MOD,m/MOD);
}

LL A()
{
LL ans=0;
if(n>0) ans=(ans+multi(Lucas(n+m-1,n-1),a1))%MOD;
if(n>1) ans=(ans+multi(Lucas(n+m-1,n-2),d))%MOD;
return ans;
}

LL b(int n)
{
if(n==0) return 0;
return multi(b1,quick_mod(q,n-1));
}

void setBC()
{
BC[1]=Lucas(m-1,0);
for(int i=2;i<=n;i++)
BC[i]=multi( multi(BC[i-1],(m+i-2)), nNi[i-1]);
}

LL B()
{
if(n==0) return b(m);
setBC();
LL FenMu=quick_mod(q-1,MOD-2);
LL anspre=b(m),ans;
for(int i=1;i<=n;i++)
{
//LL temp=multi(q,anspre)-multi(Lucas(m+i-2,i-1),b1);
LL temp=multi(q,anspre)-multi(BC[i],b1);
temp=(temp%MOD+MOD)%MOD;
ans=multi(FenMu,temp);
//cout<<"I"<<i<<"temp"<<temp<<"ans"<<ans<<endl;
anspre=ans;
}
return ans;
}

int main()
{
nJieChen[0]=1;
for(int i=1;i<10001;i++)
nJieChen[i]=multi(nJieChen[i-1],i);

nNi[1]=1;
for(int i=2;i<10001;i++)
nNi[i]=<span style="font-family: Arial, Helvetica, sans-serif;">multi</span><span style="font-family: Arial, Helvetica, sans-serif;">(MOD-MOD/i,nNi[MOD%i]);</span>

int t;
scanf("%d",&t);
for(int cases=1;cases<=t;cases++)
{
scanf("%d%d%d%d%d%d",&b1,&q,&a1,&d,&n,&m);
printf("Case #%d: %I64d\n",cases,(A()%MOD+B()%MOD+(MOD<<1))%MOD);
//printf("Case #%d: A%I64d B%I64d\n",cases,A(n,m),B(n,m));
}
return 0;
}

----

组合数取模 http://blog.csdn.net/acdreamers/article/details/8037918
逆远详解 http://blog.csdn.net/acdreamers/article/details/8220787
end