标签：a* 快速 read cal 矩阵 == opera mat ali
感觉这东西挺冷门的......
总之先记个柿子，对于整数\(n,k\)，有
\[
[n \mid k] = \frac{1}{n}\sum\limits_{i = 0} ^ {n-1} \omega_{n}^{ik}
\]
证明：
如果\(n \mid k\)，那么\(\omega_{n}^{ik} = \omega_{n}^0 =1\)，原式乘上\(\frac{1}{n}\)后显然是\(1\)。
如果\(n \not\mid k\)，那么等比数列求和
\[
\sum\limits_{i = 0} ^ {n-1} \omega_n^{ik}
= \sum\limits_{i = 0} ^ {n-1} (\omega_{n}^{k})^i
= \frac{\omega_{n}^{kn} - \omega_{n}^{0}}{\omega_{n}^{k} - 1} = 0
\]
没了。
然后反演的柿子长这样
\[
f(x) = \sum\limits_{i=0}^n a_ix^i \Leftrightarrow \sum\limits_{i=0}^n [k \mid i] a_i = \frac{1}{k} \sum\limits_{i=0}^{k-1} f(\omega_{k}^{i})
\]
这个把第一个柿子代进去就证了，这里就不写了。
下面为了简单点并不会用反演的柿子，而用第一个柿子来拆柿子。
题面太妙了......这里在说一遍题意
给整数\(n,k,p\)，让你求
\[
\sum\limits_{i=0}^{\left\lfloor \frac{n}{k} \right\rfloor} \binom{n}{ik} f_{ik} \ mod \ p
\]
其中，\(1\le n\le 10^{18}, 1\le k \le 20000\)，保证\(p\)是质数且\(p \equiv 1 \ (mod \ k)\)，\(f_i\)表示斐波那契数列的第\(i\)项（\(f_0 = 1, f_1 = 1, f_i = f_{i-1}+f_{i-2}(i>1)\)）。
我们先把下取整去掉
\[
\sum\limits_{i = 0}^n [k \mid i] \binom{n}{i}f_i
\]
组合数很想直接二项式定理......但我们有几个问题
斐波那契好解决，直接用矩乘，即\(A = \begin{bmatrix} 1 & 1\\ 1 & 0 \end{bmatrix}\)，那么\(f_n = (A^n)_{1,2} + (A^n)_{2,2}\)，即\(A^n\)的最后一列上两个数的和。
下面不看第\(2\)个限制，就是一个裸的二项式定理（二项式定理对矩阵也使用哦）
\[
\sum\limits_{i = 0} ^ n \binom{n}{i} f_i = ((A + I)^n)_{1,2} + ((A+I)^n)_{2,2}
\]
那前面整除咋办呢？注意到\(k\)只有\(20000\)，单位根反演暴力拆
\[
\begin{aligned}&\quad \sum\limits_{i = 0}^n \frac{1}{k}\sum\limits_{j=0}^{k} \omega_{k}^{ij} \binom{n}{i}f_i \\&= \frac{1}{k} \sum_{j = 0} ^ k \sum\limits_{i = 0} ^ n \binom{n}{i} (\omega_{k}^{j})^if_i\end{aligned}
\]
然后直接把\(\omega_{k}^{j}\)乘到矩阵\(A\)里，再二项式定理
\[
\frac{1}{k}\sum\limits_{j=0}^k ((A\omega_{k}^j + I)^n)_{1,2} + ((A\omega_{k}^j + I)^n)_{2,2}
\]
对\(p\)取模的话直接用原根代替单位根，再做矩阵快速幂就好了，复杂度\(O(k \ log \ n)\)。
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
int mod,K,g;ll n;
inline int fpow(ll a,ll b){
int ret=1; for (a%=mod;b;b>>=1,a=1ll*a*a%mod)
if (b&1) ret=1ll*ret*a%mod;
return ret;
}
inline int check(int k1){
for (int i=2;i*i<=mod-1;i++)
if ((mod-1)%i==0&&(fpow(k1,i)==1||fpow(k1,(mod-1)/i)==1))return 0;
return 1;
}
inline void getG(){
for (int i=2;i<=100;i++)
if (check(i)){g=i;break;}
}
inline int add(int x,int y){return x+y>=mod?x+y-mod:x+y;}
struct matrix{
int a[3][3];
matrix(int R=2,int C=2){
for (int i=0;i<R;i++)
for (int j=0;j<C;j++) a[i][j]=0;
}
int* operator [] (int k1){return a[k1];}
matrix operator * (matrix k1)const{
matrix res;
for (int i=0;i<2;i++)
for (int j=0;j<2;j++)
for (int k=0;k<2;k++)res[i][j]=add(res[i][j],1ll*a[i][k]*k1[k][j]%mod);
return res;
}
matrix operator * (ll k1)const{
k1%=mod; matrix ret;
for (int i=0;i<2;i++)
for (int j=0;j<2;j++)
ret[i][j]=1ll*a[i][j]*k1%mod;
return ret;
}
matrix operator + (matrix k1)const{
matrix ret;
for (int i=0;i<2;i++)
for (int j=0;j<2;j++) ret[i][j]=add(a[i][j],k1[i][j]);
return ret;
}
void print(){
for (int i=0;i<2;i++,puts(""))
for (int j=0;j<2;j++)
printf("%d ",a[i][j]);
}
};
inline matrix fmpow(matrix A,ll t){
matrix ret; ret[0][0]=ret[1][1]=1;
for (;t;t>>=1,A=A*A) if (t&1) ret=ret*A;
return ret;
}
int main(){
matrix I; I[0][0]=I[1][1]=1;
matrix A; A[0][0]=A[0][1]=A[1][0]=1;
int T;scanf("%d",&T);while(T--){
scanf("%lld%d%d",&n,&K,&mod);
getG(); int ans=0,wk=fpow(g,(mod-1)/K);
for (int i=0,w=1;i<K;i++,w=1ll*w*wk%mod){
matrix B=fmpow(A*w+I,n); //B.print();
ans=add(ans,B[0][1]),ans=add(ans,B[1][1]);
}
ans=1ll*ans*fpow(K,mod-2)%mod;
printf("%d\n",ans);
}
return 0;
}题意简洁明了......
给\(n,s,a_0,a_1,a_2,a_3\)，让你求
\[
\sum\limits_{i = 0} ^ n \binom{n}{i}s_ia_{i \ mod \ 4}
\]
对\(998244353\)取模，其中\(1 \le n,s,a_0,a_1,a_2,a_3 \le 10^{18}\)。
怎么又是二项式定理......
直接枚举\(i\)除以\(4\)的余数\(k\)，\(i \equiv k \ (mod \ 4)\)等价于\(4 \mid i-k\)，然后再单位根反演拆柿子，我们令\(f(k) = \sum\limits_{i=0}^n [4\mid i-k]\binom{n}{i}s_ia_i\)。
\[
\begin{aligned}f(k) &= \sum\limits_{i=0}^n \frac{1}{4}\sum\limits_{j=0}^3 \binom{n}{i}\omega_{4}^{(i-k)j}s^ia_k \\&= \frac{1}{4}a_k \sum\limits_{j=0}^3 \sum\limits_{i=0}^n \binom{n}{i}\omega_{4}^{ij-kj}s^i \\&= \frac{1}{4}a_k \sum\limits_{j=0}^3 \frac{1}{\omega_{4}^{jk}} \sum\limits_{i=0}^n \binom{n}{i} (\omega_{4}^{j})^i s^i \\&= \frac{1}{4}a_k \sum\limits_{j=0}^3 \frac{1}{\omega_{4}^{jk}} (\omega_{4}^{j}s + 1)^n\end{aligned}
\]
答案就是\(f(0) + f(1) + f(2) + f(3)\)，用快速幂计算即可。复杂度\(O(log \ n)\)
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
inline ll read(){
char ch; ll x=0; while(!isdigit(ch=getchar()));
while(isdigit(ch)){x=x*10+ch-48;ch=getchar();}return x;
}
const int P=998244353,G=3,IG=(P+1)/G;
inline int fpow(int a,ll b){
int ret=1; for (;b;b>>=1,a=1ll*a*a%P)
if (b&1) ret=1ll*a*ret%P;
return ret;
}
const int i4=fpow(4,P-2);
inline int add(int x,int y){return x+y>=P?x+y-P:x+y;}
ll n; int a[4],s;
int calc(int k){
int ans=0,w4=fpow(G,(P-1)/4),iw4=fpow(IG,(P-1)/4);
for (int j=0;j<4;j++)
ans=add(ans,1ll*fpow(iw4,j*k)*fpow(add(1ll*s*fpow(w4,j)%P,1),n)%P);
ans=1ll*ans*a[k]%P;
return ans;
}
int main(){
for (int _=read();_;_--){
int ans=0;
n=read(),s=read()%P,a[0]=read()%P,a[1]=read()%P,a[2]=read()%P,a[3]=read()%P;
for (int k=0;k<4;k++)ans=add(ans,calc(k));
printf("%d\n",1ll*i4*ans%P);
}
return 0;
}还有题目先鸽着吧......看心情会更的......
标签：a* 快速 read cal 矩阵 == opera mat ali
原文地址：https://www.cnblogs.com/wxq1229/p/12242177.html
