Solution -「LOJ #138」「模板」类欧几里得算法
2021/7/15 9:06:21
本文主要是介绍Solution -「LOJ #138」「模板」类欧几里得算法,对大家解决编程问题具有一定的参考价值,需要的程序猿们随着小编来一起学习吧!
\(\mathcal{Description}\)
Link.
\(T\) 组询问,每次给出 \(n,a,b,c,k_1,k_2\),求
\[\sum_{x=0}^nx^{k_1}\left\lfloor\frac{ax+b}{c}\right\rfloor^{k_2}\bmod(10^9+7) \]\(T=1000\),\(n,a,b,c\le10^9\),\(0\le k_1+k_2\le 10\)。
\(\mathcal{Solution}\)
类欧模板题的集大成者。
令 \(f(n,a,b,c,k_1,k_2)\) 表示所求答案,尝试通过分类讨论递归到不同的情况求解。
-
\(a=0\lor an+b<c\)
\[f(n,a,b,c,k_1,k_2)=\left\lfloor\frac{an+b}{c}\right\rfloor^{k_2}\sum_{x=0}^nx^{k_1} \]预处理 \(k+1\) 次多项式 \(g_k(x)=\sum_{i=0}^xi^k\) 的各项系数,这个就能直接算出来。
-
否则若 \(a\ge c\)
\[\begin{aligned} f(n,a,b,c,k_1,k_2)&=\sum_{x=0}^nx^{k_1} \left(\left\lfloor\frac{a}{c}\right\rfloor x+\left\lfloor\frac{(a\bmod c)x+b}{c}\right\rfloor \right)^{k_2}\\ &=\sum_{i=0}^{k_2} \binom{k_2}{i} \left\lfloor\frac{a}{c}\right\rfloor^i\sum_{x=0}^n x^{k_1+i}\left\lfloor\frac{(a\bmod c)x+b}{c}\right\rfloor^{k_2-i}\\ &=\sum_{i=0}^{k_2} \binom{k_2}{i} \left\lfloor\frac{a}{c}\right\rfloor^i f(n,a\bmod c,b,c,k_1+i,k_2-i) \end{aligned} \]递归处理。
-
否则若 \(b\ge c\)
类似 2.,最终得到
\[f(n,a,b,c,k_1,k_2)=\sum_{i=0}^{k_2} \binom{k_2}{i} \left\lfloor\frac{b}{c}\right\rfloor^i f(n,a,b\bmod c,c,k_1,k_2-i) \]递归处理。
-
对于其余情况
考虑把高斯函数丢到求和指标上,注意到
\[m^k=\sum_{j=0}^{m-1}[(j+1)^k-j^k] \]以此替代 \(\left\lfloor\frac{ax+b}{c}\right\rfloor^{k_2}\),并交换求和顺序,得到
\[\begin{aligned} f(n,a,b,c,k_1,k_2)&=\sum_{j=0}^{\left\lfloor\frac{an+b}{c}\right\rfloor-1}[(j+1)^{k_2}-j^{k_2}]\sum_{x=0}^nx^{k_1}\left[x>\left\lfloor\frac{cj+c-b-1}{a}\right\rfloor \right]\\ &= \sum_{j=0}^{\left\lfloor\frac{an+b}{c}\right\rfloor-1} [(j+1)^{k_2}-j^{k_2}]\sum_{x=0}^n x^{k_1}-\sum_{j=0}^{\left\lfloor\frac{an+b}{c}\right\rfloor-1} [(j+1)^{k_2}-j^{k_2}]\sum_{x=0}^{\left\lfloor\frac{cj+c-b-1}{a}\right\rfloor} x^{k_1} \end{aligned} \]前半部分直接计算,研究后半部分,发现最后一个求和符号是 \(g_{k_1}\left(\left\lfloor\frac{cj+c-b-1}{a}\right\rfloor\right)\),枚举次数将其展开:
\[\begin{aligned} &=\sum_{i=0}^{k_2-1}\binom{k_2}{i} \sum_{j=0}^{k_1+1}([x^j]g_{k_1})\sum_{x=0}^{\left\lfloor\frac{an+b}{c}\right\rfloor-1}x^i\left\lfloor\frac{cx+c-b-1}{a}\right\rfloor^j\\ &=\sum_{i=0}^{k_2-1}\binom{k_2}{i} \sum_{j=0}^{k_1+1}([x^j]g_{k_1})f\left(\left\lfloor\frac{an+b}{c}\right\rfloor-1,c,c-b-1,a,i,j\right) \end{aligned} \]其中 \(i\) 枚举 \((j+1)^{k_2}\) 的次数;\(j\) 枚举 \(g_{k_1}\) 的次数;\(x\) 枚举原来的 \(j\)。这样也能递归求解了。
递归中指标变换形如欧几里得算法,设指标值域为 \(V\),\(K=k_1+k_2\),则复杂度为 \(\mathcal O(TK^4\log V)\)。
\(\mathcal{Code}\)
实现时,可以对于 \(n,a,b,c\),一次算出所有可能的 \(f(n,a,b,c,k_1,k_2)\) 的值。
/*~Rainybunny~*/ #include <cstdio> #include <cassert> #include <iostream> #define rep( i, l, r ) for ( int i = l, rep##i = r; i <= rep##i; ++i ) #define per( i, r, l ) for ( int i = r, per##i = l; i >= per##i; --i ) typedef long long LL; const int MOD = 1e9 + 7; int comb[15][15]; inline void subeq( int& a, const int b ) { ( a -= b ) < 0 && ( a += MOD ); } inline int sub( int a, const int b ) { return ( a -= b ) < 0 ? a + MOD : a; } inline int mul( const long long a, const int b ) { return int( a * b % MOD ); } inline int add( int a, const int b ) { return ( a += b ) < MOD ? a : a - MOD; } inline void addeq( int& a, const int b ) { ( a += b ) >= MOD && ( a -= MOD ); } inline int mpow( int a, int b ) { int ret = 1; for ( ; b; a = mul( a, a ), b >>= 1 ) ret = mul( ret, b & 1 ? a : 1 ); return ret; } struct PowerPoly { int n, a[12]; PowerPoly() {} PowerPoly( const int k ): n( k ) { static int mat[15][15]; for ( int i = 0, p = 0; i <= n + 1; ++i ) { addeq( p, mpow( i, n ) ); mat[i][n + 2] = p; } rep ( i, 0, n + 1 ) { for ( int j = 0, p = 1; j <= n + 1; ++j, p = mul( p, i ) ) { mat[i][j] = p; } } rep ( i, 0, n + 1 ) { int p = -1; rep ( j, i, n + 1 ) if ( mat[j][i] ) { p = j; break; } assert( ~p ); if ( p != i ) std::swap( mat[i], mat[p] ); int iv = mpow( mat[i][i], MOD - 2 ); rep ( j, i + 1, n + 1 ) { int t = mul( iv, mat[j][i] ); rep ( l, i, n + 2 ) subeq( mat[j][l], mul( mat[i][l], t ) ); } } per ( i, n + 1, 0 ) { a[i] = mul( mat[i][n + 2], mpow( mat[i][i], MOD - 2 ) ); rep ( j, 0, i - 1 ) subeq( mat[j][n + 2], mul( a[i], mat[j][i] ) ); } } inline int operator () ( LL x ) const { int ret = 0; x %= MOD; for ( int i = 0, p = 1; i <= n + 1; ++i, p = mul( x, p ) ) { addeq( ret, mul( a[i], p ) ); } return ret; } }; const PowerPoly pwr[11] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 }; struct Atom { int res[11][11]; Atom(): res{} {} inline int* operator [] ( const int k ) { return res[k]; } }; /* * * calculate \sum_{x=0}^n x^{k_1} \lfloor \frac{ax+b}{c} \rfloor^{k_2}, * where k_1+k_2 <= 10. * */ inline Atom solve( const LL n, const LL a, const LL b, const LL c ) { Atom ret; LL m = ( a * n + b ) / c; if ( !a || 1ll * a * n + b < c ) { rep ( k1, 0, 10 ) { int t = int( m % MOD ), p = 1; rep ( k2, 0, 10 - k1 ) { ret[k1][k2] = mul( p, pwr[k1]( n ) ); p = mul( p, t ); } } return ret; } if ( a >= c ) { Atom tmp( solve( n, a % c, b, c ) ); rep ( k1, 0, 10 ) rep ( k2, 0, 10 - k1 ) { int t = int( ( a / c ) % MOD ), p = 1; rep ( i, 0, k2 ) { addeq( ret[k1][k2], mul( comb[k2][i], mul( p, tmp[k1 + i][k2 - i] ) ) ); p = mul( p, t ); } } return ret; } if ( b >= c ) { Atom tmp( solve( n, a, b % c, c ) ); rep ( k1, 0, 10 ) rep ( k2, 0, 10 ) { int t = int( ( b / c ) % MOD ), p = 1; rep ( i, 0, k2 ) { addeq( ret[k1][k2], mul( comb[k2][i], mul( p, tmp[k1][k2 - i] ) ) ); p = mul( p, t ); } } return ret; } Atom tmp( solve( m - 1, c, c - b - 1, a ) ); rep ( k1, 0, 10 ) { int t = int( m % MOD ), p = 1; rep ( k2, 0, 10 - k1 ) { ret[k1][k2] = mul( p, pwr[k1]( n ) ); rep ( i, 0, k2 - 1 ) rep ( j, 0, k1 + 1 ) { subeq( ret[k1][k2], mul( comb[k2][i], mul( pwr[k1].a[j], tmp[i][j] ) ) ); } p = mul( t, p ); } } return ret; } int main() { comb[0][0] = 1; rep ( i, 1, 10 ) { comb[i][0] = 1; rep ( j, 1, i ) comb[i][j] = add( comb[i - 1][j - 1], comb[i - 1][j] ); } int n, a, b, c, k1, k2, T; for ( scanf( "%d", &T ); T--; ) { scanf( "%d %d %d %d %d %d", &n, &a, &b, &c, &k1, &k2 ); printf( "%d\n", solve( n, a, b, c )[k1][k2] ); } return 0; }
这篇关于Solution -「LOJ #138」「模板」类欧几里得算法的文章就介绍到这儿,希望我们推荐的文章对大家有所帮助,也希望大家多多支持为之网!
- 2024-11-01后台管理开发学习:新手入门指南
- 2024-11-01后台管理系统开发学习:新手入门教程
- 2024-11-01后台开发学习:从入门到实践的简单教程
- 2024-11-01后台综合解决方案学习:从入门到初级实战教程
- 2024-11-01接口模块封装学习入门教程
- 2024-11-01请求动作封装学习:新手入门教程
- 2024-11-01登录鉴权入门:新手必读指南
- 2024-11-01动态面包屑入门:轻松掌握导航设计技巧
- 2024-11-01动态权限入门:新手必读指南
- 2024-11-01动态主题处理入门:新手必读指南