51nod 1362 搬箱子——[ 推式子+組合數計算方法 ] [ 拉格朗日插值 ]
題目：http://www.51nod.com/Challenge/Problem.html#!#problemId=1362
方法一：
設 a 是向下走的步數、 b 是向右下走的步數、 c 是向下走的步數。如果是走到第 j 列的方案數的話，有：
\( a+b = n \) \( b+c = j \)
所以列舉 a 和 j ， b 和 c 的值就是確定的，可以用組合數算：
\( \sum\limits_{i=0}^{n}\sum\limits_{j=0}^{m}C_{i+j}^{i}*C_{j}^{n-i} \)
\( = \sum\limits_{i=0}^{n}\sum\limits_{j=0}^{m} \frac{(i+j)!}{i!j!} * \frac{j!}{(n-i)!(i+j-n)!} \)
\( = \sum\limits_{i=0}^{n}\sum\limits_{j=0}^{m} \frac{(i+j)!}{i!(n-i)!(i+j-n)!} \)
一看就是要把只和 i 有關的提到前面。而且分子和分母好像能湊成新的組合數，所以弄一個 \( n! \)
\( = \sum\limits_{i=0}^{n} \frac{n!}{i!(n-i)!} \sum\limits_{j=0}^{m} \frac{(i+j)!}{n!(i+j-n)!} \)
\( = \sum\limits_{i=0}^{n} C_{n}^{i} \sum\limits_{j=0}^{m} C_{i+j}^{n} \)
\( = \sum\limits_{i=0}^{n} C_{n}^{i} * C_{i+m+1}^{n+1} \)
但 i+m+1 很大，而且模數也不是質數。所以用擴充套件Lucas。
但它的 pk 可以很大，會TLE。總之弄了半天還是會T一個點。反正本來複雜度也不對……
#include<iostream> #include<cstdio> #include<cstring> #include<algorithm> #define ll long long using namespace std; constView Codeint N=805,M=1e5+5; int n,m,mod,tot,p[M],pk[M],ans; void init() { int d=mod; tot=0; for(int i=2;i*i<=d;i++) if(d%i==0) { p[++tot]=i;pk[tot]=1; while(d%i==0)d/=i,pk[tot]*=i; } if(d>1)p[++tot]=pk[tot]=d; } int pw(int x,int k,int md) {int ret=1;while(k){if(k&1)ret=(ll)ret*x%mod;x=(ll)x*x%mod;k>>=1;}return ret;} void exgcd(int a,int b,int &x,int &y) {if(!b){x=1;y=0;return;} exgcd(b,a%b,y,x); y-=a/b*x;} int inv(int a,int md){int x,y;exgcd(a,md,x,y);return (x+md)%md;} int multi(int n,int p,int pk) { if(!n)return 1; int sm=1; for(int i=2;i<pk;i++)if(i%p)sm=(ll)sm*i%pk;//if sm=pw(sm,n/pk,pk); for(int i=2,j=n%pk;i<=j;i++)if(i%p)sm=(ll)sm*i%pk;// return (ll)sm*multi(n/p,p,pk)%pk; } int exlcs(int n,int m,int p,int pk) { int sm=0; for(int i=n;i;i/=p)sm+=i/p; for(int i=m;i;i/=p)sm-=i/p; for(int i=n-m;i;i/=p)sm-=i/p; return (ll)pw(p,sm,pk)*multi(n,p,pk)%pk*inv(multi(m,p,pk),pk)%pk*inv(multi(n-m,p,pk),pk)%pk; } int C(int n,int m,int p) { if(n<m)return 0; int ret,sm=1; for(int i=2;i<=n;i++)sm=(ll)sm*i%p; ret=sm; sm=1; for(int i=2;i<=m;i++)sm=(ll)sm*i%p; ret=(ll)ret*pw(sm,p-2,p)%p; sm=1; for(int i=2;i<=n-m;i++)sm=(ll)sm*i%p; ret=(ll)ret*pw(sm,p-2,p)%p; return ret; } int lcs(int n,int m,int p) { if(n<m)return 0; if(!m||n==m)return 1; return (ll)C(n%p,m%p,p)*lcs(n/p,m/p,p)%p; } int calc(int n,int m) { if(n<m)return 0;///// int ret=0; for(int i=1;i<=tot;i++) { if(p[i]==pk[i]) { ret=(ret+(ll)lcs(n,m,p[i])*(mod/pk[i])%mod*inv(mod/pk[i],pk[i]))%mod; } else ret=(ret+(ll)exlcs(n,m,p[i],pk[i])*(mod/pk[i])%mod*inv(mod/pk[i],pk[i]))%mod; } return ret; } int main() { while(scanf("%d%d%d",&n,&m,&mod)==3) { init(); ans=0; for(int i=0;i<=n;i++) { ans=(ans+(ll)calc(n,i)*calc(i+m+1,n+1))%mod; } printf("%d\n",ans); } return 0; }
發現別人都不是這樣寫的。
因為 \( C_{i+m+1}^{n+1} \) 雖然 i+m+1 很大，但n+1 很小，所以可以列舉分子上的數（約掉分母裡很大的那一項之後只剩下 n 項了！），一邊把含的 p 拿出來之類的。
注意不要用求階乘裡 p 的個數的方法求好總的 p 的個數之後在列舉的時候跳過 %p==0 的數。因為可能那個數除掉一些 p 之後有剩下的。
#include<iostream> #include<cstdio> #include<cstring> #include<algorithm> #define ll long long using namespace std; const int N=805,M=1e5+5; int n,m,mod,tot,p[M],pk[M],ans; void init() { int d=mod; tot=0; for(int i=2;i*i<=d;i++) if(d%i==0) { p[++tot]=i;pk[tot]=1; while(d%i==0)d/=i,pk[tot]*=i; } if(d>1)p[++tot]=pk[tot]=d; } int pw(int x,int k,int md) {int ret=1;while(k){if(k&1)ret=(ll)ret*x%mod;x=(ll)x*x%mod;k>>=1;}return ret;} void exgcd(int a,int b,int &x,int &y) {if(!b){x=1;y=0;return;} exgcd(b,a%b,y,x); y-=a/b*x;} int inv(int a,int md){int x,y;exgcd(a,md,x,y);return (x%md+md)%md;} int C(int n,int m,int p,int pk) { int nm=0,ret=1; for(int i=n,j=1;j<=m;i--,j++) { int a=i,b=j; while(a%p==0)nm++,a/=p; while(b%p==0)nm--,b/=p; ret=(ll)ret*a%pk*inv(b,pk)%pk; } ret=(ll)ret*pw(p,nm,pk)%pk; return ret; } int calc(int n,int m) { if(n<m)return 0;///// int ret=0; for(int i=1;i<=tot;i++) { ret=(ret+(ll)C(n,m,p[i],pk[i])*(mod/pk[i])%mod*inv(mod/pk[i],pk[i]))%mod; } return ret; } int main() { while(scanf("%d%d%d",&n,&m,&mod)==3) { init(); ans=0; for(int i=0;i<=n;i++) { ans=(ans+(ll)calc(n,i)*calc(i+m+1,n+1))%mod; } printf("%d\n",ans); } return 0; }View Code
方法二：
設 dp[ i ][ j ] 表示走到第 i 行第 j 列的方案數，則有 dp[ i ][ j ] = dp[ i-1 ][ j ] + dp[ i ][ j-1 ] + dp[ i-1 ][ j-1 ] ；
用差分的方法判斷一下，發現 dp[ i ] 是一個 i 次的多項式。（設原數列為第0次差分後數列，若差分 k 次後數列變成常數列（即每一項值都一樣），則為 k 次多項式）
設 \( s[ i ][ j ] = \sum\limits_{k=0}^{j} dp[ i ][ k ] \) ，則 s[ i ] 是一個 i+1 次多項式。
#include<iostream> #include<cstdio> #include<cstring> #include<algorithm> #define ll long long using namespace std; const int N=805,M=35; int n,m,mod,dp[N][N]; void upd(int &x){x>=mod?x-=mod:0;} int main() { while(scanf("%d%d%d",&n,&m,&mod)==3) { memset(dp,0,sizeof dp); int d=min(n+1,m); dp[0][0]=1; for(int i=0;i<=n;i++) for(int j=0;j<=d;j++) { if(i)dp[i][j]+=dp[i-1][j],upd(dp[i][j]); if(j)dp[i][j]+=dp[i][j-1],upd(dp[i][j]); if(i&&j)dp[i][j]+=dp[i-1][j-1],upd(dp[i][j]); } for(int j=0;j<=d;j++)dp[n][j]+=dp[n][j-1],upd(dp[n][j]); int lm=1000; for(int i=1,R=d-1;i<=lm;i++,R--) { int flag=1; for(int j=0;j<=R;j++) { dp[n][j]=dp[n][j+1]-dp[n][j]; if(j&&dp[n][j]!=dp[n][j-1])flag=0; } if(flag){printf("%d\n",i);break;} } } return 0; }判斷
所以可以用拉格朗日插值。
注意 ( i - j ) 可能沒有逆元，同樣採用把 mod 質因數（只有log(mod)個質因數！）分解、最後乘上 pk 的方法。
注意要除的東西不能先累乘起來最後再除掉。
#include<iostream> #include<cstdio> #include<cstring> #include<algorithm> #define ll long long using namespace std; const int N=805,M=35; int n,m,mod,tot,dp[N][N],p[M],nm[M],phi; void upd(int &x){x>=mod?x-=mod:0;} int pw(int x,int k) {int ret=1;while(k){if(k&1)ret=(ll)ret*x%mod;x=(ll)x*x%mod;k>>=1;}return ret;} void init() { tot=0; phi=mod; int k=mod; for(int i=2;i*i<=k;i++) if(k%i==0) { p[++tot]=i; phi/=i; phi*=(i-1); while(k%i==0)k/=i; } if(k>1)p[++tot]=k,phi/=k,phi*=(k-1); } void add(int a,int fx,int &ret) { for(int i=1;i<=tot;i++) { while(a%p[i]==0)nm[i]+=fx,a/=p[i]; } fx==1? ret=(ll)ret*a%mod : ret=(ll)ret*pw(a,phi-1)%mod ; } int cz(int ret) { for(int i=1;i<=tot;i++) { ret=(ll)ret*pw(p[i],nm[i])%mod; nm[i]=0; } return ret; } int main() { while(scanf("%d%d%d",&n,&m,&mod)==3) { memset(dp,0,sizeof dp); int d=min(n+1,m); dp[0][0]=1; for(int i=0;i<=n;i++) for(int j=0;j<=d;j++) { if(i)dp[i][j]+=dp[i-1][j],upd(dp[i][j]); if(j)dp[i][j]+=dp[i][j-1],upd(dp[i][j]); if(i&&j)dp[i][j]+=dp[i-1][j-1],upd(dp[i][j]); } for(int i=1;i<=d;i++) dp[n][i]+=dp[n][i-1],upd(dp[n][i]); if(m==d){printf("%d\n",dp[n][m]);continue;} init(); int ans=0; for(int i=0;i<=d;i++) { int ret=1; for(int j=0;j<=d;j++) { if(i==j)continue; add(m-j,1,ret); add(i-j,-1,ret); } ans=(ans+(ll)dp[n][i]*cz(ret))%mod; } if(ans<0)ans+=mod; printf("%d\n",ans); } return 0; }View Code