锐客网

  1. 首页 > it技术 > >

数学|LOJ #138. 类欧几里得算法

数学离散数学数论计数

题目
推导方式比较牛逼，一劳永逸
推公式最详细的博客
消去公式：
F ( k 1 , k 2 , n , a , b , c ) = ∑ i = 0 n i k 1 ? a i + b c ? k 2 = ∑ i = 0 n i k 1 ( ? a c ? i + ? b c ? + ? ( a m o d ?? c ) i + ( b m o d ?? c ) c ? ) k 2 = ∑ i = 0 n i k 1 ∑ ∑ p = 1 3 j p = k 2 k ! ∏ p = 1 3 j p ! × ( ? a c ? i ) j 1 ? b c ? j 2 ? ( a m o d ?? c ) i + ( b m o d ?? c ) c ? j 3 = ∑ ∑ p = 1 3 j p = k 2 k ! ? a c ? j 1 ? b c ? j 2 ∏ p = 1 3 j p ! F ( k 1 + j 1 , j 3 , n , a % c , b % c , c ) \begin{aligned}F(k_1,k_2,n,a,b,c) &= \sum_{i=0}^n i^{k_1} \left\lfloor\frac{ai+b}c\right\rfloor^{k_2} \\ &= \sum_{i=0}^ni^{k_1}(\lfloor \frac ac \rfloor i+\lfloor \frac bc \rfloor+\lfloor \frac {(a\mod c)i+(b\mod c)}c\rfloor)^{k_2} \\ &=\sum_{i=0}^ni^{k_1} \sum_{\sum_{p=1}^3{j_p} = k_2} \frac{k!}{\prod_{p=1}^3 j_p!} \times\\& \left(\left\lfloor \frac ac\right \rfloor i\right)^{j_1} \left\lfloor \frac bc \right\rfloor ^{j_2} \left\lfloor \frac {(a\mod c)i+(b\mod c)}c \right\rfloor^{j_3}\\ &=\sum_{\sum_{p=1}^3{j_p} = k_2} \frac{k!\left\lfloor\frac ac \right \rfloor^{j_1}\left \lfloor\frac bc \right \rfloor^{j_2}}{\prod_{p=1}^3 j_p!} F(k_1+j_1,j_3,n,a\% c,b\% c ,c) \end{aligned} F(k1?,k2?,n,a,b,c)?=i=0∑n?ik1??cai+b??k2?=i=0∑n?ik1?(?ca??i+?cb??+?c(amodc)i+(bmodc)??)k2?=i=0∑n?ik1?∑p=13?jp?=k2?∑?∏p=13?jp?!k!?×(?ca??i)j1??cb??j2??c(amodc)i+(bmodc)??j3?=∑p=13?jp?=k2?∑?∏p=13?jp?!k!?ca??j1??cb??j2??F(k1?+j1?,j3?,n,a%c,b%c,c)?
反转公式：
令 m = ? a n + b c ? ? 1 m = \left \lfloor \frac {an+b}c \right\rfloor - 1 m=?can+b???1
F ( k 1 , k 2 , n , a , b , c ) = ∑ i = 0 n i k 1 ? a i + b c ? k 2 = ∑ i = 0 n i k 1 ∑ ? p ∈ [ 1 , k 2 ] , j p ∈ [ 0 , m ] ∏ p = 1 k 2 [ j p < ? a i + b c ? ] ( 枚 举 j p 的 最 大 值 ) = ∑ i = 0 n i k 1 ∑ q = 0 m [ q < ? a i + b c ? ] [ ( q + 1 ) k 2 ? q k 2 ] ( q < ? a i + b c ? 等 价 于 i > ? c q + c ? b ? 1 a ? ) = ∑ q = 0 m [ ( q + 1 ) k 2 ? q k 2 ] ( ∑ i = 0 n i k ? ∑ i = 0 ? c q + c ? b ? 1 a ? i k ) = [ ( m + 1 ) k 2 ? 0 k 2 ] ∑ i = 0 n i k ? ∑ q = 0 m ∑ r = 0 k 2 ? 1 ( k 2 r ) q r ∑ i = 0 ? c q + c ? b ? 1 a ? i k ( 将 函 数 S k ( x ) = ∑ i = 0 x i k 插 值 后 得 到 其 k + 1 次 表 达 式 写 作 S k ( x ) = ∑ i = 0 k + 1 a k , i x i ) = [ ( m + 1 ) k 2 ? 0 k 2 ] ∑ i = 0 n i k ? ∑ r = 0 k 2 ? 1 ∑ i = 0 k 1 + 1 ( k 2 i ) a k 1 , i ∑ q = 0 m q r ? c q + c ? b ? 1 a ? i = [ ( m + 1 ) k 2 ? 0 k 2 ] ∑ i = 0 n i k ? ∑ r = 0 k 2 ? 1 ∑ i = 0 k 1 + 1 ( k 2 i ) a k 1 , i F ( r , i , m , c , c ? b ? 1 , a ) \begin{aligned} F(k_1,k_2,n,a,b,c) &= \sum_{i=0}^n i^{k_1} \left\lfloor\frac{ai+b}c\right\rfloor^{k_2} \\ &=\sum_{i=0}^n i^{k_1}\sum_{\forall p\in[1,k_2],j_p\in[0,m]} \prod_{p=1}^{k_2} [j_p<\left \lfloor \frac{ai+b}c \right \rfloor ]\\ &(枚举j_p的最大值)\\ &=\sum_{i=0}^n i^{k_1} \sum_{q=0}^m \left[q < \left \lfloor \frac {ai+b}c \right \rfloor\right]\left[(q+1)^{k_2}-q^{k_2}\right]\\ &(q < \left \lfloor \frac {ai+b}c \right \rfloor等价于i>\left \lfloor \frac {cq+c-b-1}a \right \rfloor) \\ &=\sum_{q=0}^m\left[(q+1)^{k_2}-q^{k_2}\right ] (\sum_{i=0}^n i^k - \sum_{i=0}^{\left \lfloor \frac {cq+c-b-1}a \right \rfloor} i^k)\\ &=\left[(m+1)^{k_2}-0^{k_2}\right ]\sum_{i=0}^n i^k - \\&\sum_{q=0}^m\sum_{r=0}^{k_2-1}\binom{k_2}{r}q^r\sum_{i=0}^{\left \lfloor \frac {cq+c-b-1}a \right \rfloor} i^k\\ &(将函数S_k(x) = \sum_{i=0}^x i^k插值后得到其k+1次表达式写作\\ &S_k(x) = \sum_{i=0}^{k+1} a_{k,i}x^i)\\ &=\left[(m+1)^{k_2}-0^{k_2}\right ]\sum_{i=0}^n i^k-\\ &\sum_{r=0}^{k_2-1}\sum_{i=0}^{k_1+1} \binom{k_2}{i}a_{k_1,i}\sum_{q=0}^mq^{r}\left\lfloor\frac {cq+c-b-1}a\right\rfloor^i\\ &=\left[(m+1)^{k_2}-0^{k_2}\right ]\sum_{i=0}^n i^k-\\ &\sum_{r=0}^{k_2-1}\sum_{i=0}^{k_1+1} \binom{k_2}{i}a_{k_1,i}F(r,i,m,c,c-b-1,a)\\ \end{aligned} F(k1?,k2?,n,a,b,c)?=i=0∑n?ik1??cai+b??k2?=i=0∑n?ik1??p∈[1,k2?],jp?∈[0,m]∑?p=1∏k2??[jp??acq+c?b?1??)=q=0∑m?[(q+1)k2??qk2?](i=0∑n?ik?i=0∑?acq+c?b?1???ik)=[(m+1)k2??0k2?]i=0∑n?ik?q=0∑m?r=0∑k2??1?(rk2??)qri=0∑?acq+c?b?1???ik(将函数Sk?(x)=i=0∑x?ik插值后得到其k+1次表达式写作Sk?(x)=i=0∑k+1?ak,i?xi)=[(m+1)k2??0k2?]i=0∑n?ik?r=0∑k2??1?i=0∑k1?+1?(ik2??)ak1?,i?q=0∑m?qr?acq+c?b?1??i=[(m+1)k2??0k2?]i=0∑n?ik?r=0∑k2??1?i=0∑k1?+1?(ik2??)ak1?,i?F(r,i,m,c,c?b?1,a)?
这个公式真的一点都不吓人。
但是我们需要注意我们在第二行做了的这个操作。
它在前后并不是完全相等（完全恒等）的。
当 k 2 = 0 k_2=0 k2?=0时会出现我们 ∑ ? p ∈ [ 1 , k 2 ] , j p ∈ [ 0 , m ] ∏ p = 1 k 2 [ j p < ? a i + b c ? ] \sum_{\forall p\in[1,k_2],j_p\in[0,m]} \prod_{p=1}^{k_2} [j_p<\left \lfloor \frac{ai+b}c \right \rfloor ] ∑?p∈[1,k2?],jp?∈[0,m]?∏p=1k2??[jp? 但是我们经过简单的计算和验证。
我们只需要让
F ( k 1 , k 2 , n , a , b , c ) = ( m + 1 ) k 2 ∑ i = 0 n i k ? ∑ r = 0 k 2 ? 1 ∑ i = 0 k 1 + 1 ( k 2 i ) a k 1 , i F ( r , i , m , c , c ? b ? 1 , a ) F(k_1,k_2,n,a,b,c) =(m+1)^{k_2}\sum_{i=0}^n i^k-\\\sum_{r=0}^{k_2-1}\sum_{i=0}^{k_1+1} \binom{k_2}{i}a_{k_1,i}F(r,i,m,c,c-b-1,a) F(k1?,k2?,n,a,b,c)=(m+1)k2?i=0∑n?ik?r=0∑k2??1?i=0∑k1?+1?(ik2??)ak1?,i?F(r,i,m,c,c?b?1,a)
即可解决问题。
注：实际上边界好像只用判 a = 0 a=0 a=0，别的判了也没多快，记忆化是真的快。
【数学|LOJ #138. 类欧几里得算法】 A CC o d e \mathrm {AC \ Code} AC Code

#include #define maxn 15 #define mod 1000000007 #define LL long long #define rep(i,j,k) for(int i=(j),LIM=(k); i<=LIM; i++) #define per(i,j,k) for(int i=(j),LIM=(k); i>=LIM; i--) using namespace std; int inv[maxn]={1,1},fac[maxn]={1,1},invf[maxn]={1,1}; int Pow(int b,int k){ int r = 1; for(; k; k>>=1,b=1ll*b*b%mod) if(k&1) r=1ll*r*b%mod; return r; }int C(int n,int m){ static int c[maxn][maxn]={},f=0; if(!f){ c[0][0]=f=1; rep(i,1,maxn-1) rep(j,c[i][0]=1,i) c[i][j]=(c[i-1][j]+c[i-1][j-1])%mod; } return c[n][m]; }int a[maxn][maxn]; int S(int n,int m){ n %= mod; LL r = 0 , pw = 1; rep(j,0,m+1) r = (r + 1ll * a[m][j] * pw) , pw = 1ll * pw * n % mod; return (r % mod + mod) % mod; }int f[maxn][maxn][50]; int F(int k1,int k2,int n,int a,int b,int c,int dep=0){ if(~f[k1][k2][dep]) return f[k1][k2][dep]; if(a == 0) return f[k1][k2][dep] = (1ll * Pow(b/c,k2) * S(n,k1) % mod + mod) % mod; if(a>=c||b>=c){ int r = 0 , A = a % c , B = b % c; a /= c , b /= c; int pwa[maxn]={},pwb[maxn]={}; rep(j,pwa[0]=pwb[0]=1,k2) pwa[j]=1ll*pwa[j-1]*a%mod, pwb[j]=1ll*pwb[j-1]*b%mod; rep(j1,0,k2) rep(j2,0,k2-j1){ int j3 = k2 - j1 - j2; r = (r + 1ll * pwa[j1] * pwb[j2] % mod * invf[j1] % mod * invf[j2] % mod * invf[j3] % mod * F(k1+j1,j3,n,A,B,c,dep+1)) % mod; } return f[k1][k2][dep] = ((1ll * r * fac[k2] % mod)+mod)%mod; } int m = (1ll * a * n + b) / c - 1 , r = 1ll * Pow(m+1,k2) * S(n,k1) % mod; rep(i,0,k2-1) rep(j,0,k1+1) r = (r - 1ll * C(k2 , i) * ::a[k1][j] % mod * F(i,j,m,c,c-b-1,a,dep+1)) % mod; return f[k1][k2][dep]=(r+mod)%mod; }int main(){ rep(i,2,maxn-1) inv[i]=1ll*(mod-mod/i)*inv[mod%i]%mod, fac[i]=1ll*fac[i-1]*i%mod, invf[i]=1ll*invf[i-1]*inv[i]%mod; a[0][1] = a[0][0] = 1; rep(i,1,maxn-2){ rep(j,0,i+1) a[i][j] = C(i+1,j); rep(j,0,i-1) rep(k,0,j+1) a[i][k] = (a[i][k] - 1ll * C(i+1,j) * a[j][k]) % mod; rep(j,0,i+1) a[i][j] = 1ll * a[i][j] * inv[i+1] % mod; } int T,n,a,b,c,k1,k2; for(scanf("%d",&T); T--; ){ scanf("%d%d%d%d%d%d",&n,&a,&b,&c,&k1,&k2); memset(f,-1,sizeof f); printf("%d\n",(F(k1,k2,n,a,b,c)+mod)%mod); } }

    推荐阅读


    上一篇：模板|类欧几里得算法

    下一篇：Android|第4章 View的工作原理