博客
关于我
强烈建议你试试无所不能的chatGPT,快点击我
BZOJ3328: PYXFIB(单位根反演?)
阅读量:6691 次
发布时间:2019-06-25

本文共 3114 字,大约阅读时间需要 10 分钟。

Sol

\[A=\begin{bmatrix}1 & 1 \\ 1 & 0\end{bmatrix}\]
那么要求的相当于是
\[\sum_{i=0}^{n}[k|i]\binom{n}{i}A^i\]
求出其中的 \(A_{0,0}\) 即可

引入单位根(单位根反演?)

\[[n \mid k]=\frac{1}{n}\sum_{i=0}^{n-1} \omega_{n}^{ik}\]

证明:
\(n|k\),那么根据单位根的消去引理可以得到就是 \(1\)
否则,等比数列求和,发现分子为 \(0\)

所以带入单位根

\[\sum_{i=0}^{n}[k|i]\binom{n}{i}A^i=\frac{1}{k}\sum_{i=0}^{n}\binom{n}{i}A^i\sum_{j=0}^{k-1}w_k^{ij}=\frac{1}{k}\sum_{i=0}^{k-1}\sum_{j=0}^{n}\binom{n}{j}A^jw_k^{ij}\]
由二项式定理得到
\[\frac{1}{k}\sum_{i=0}^{k-1}(Aw_k^i+I)^n\]
\(I\) 为单位矩阵
这道题 \(k=p-1\)
所以直接求出 \(p\) 的原根 \(g\),令 \(w_k^i=g^{\frac{p-1}{k}i}\) 即可
只要矩乘+快速幂就好了

# include 
using namespace std;typedef long long ll;int k, test, mod, g, pr[233333], cnt;ll n;struct Matrix { ll a[2][2]; inline Matrix() { a[0][0] = a[0][1] = a[1][0] = a[1][1] = 0; } inline Matrix operator *(Matrix b) const { register Matrix c; c.a[0][0] = (a[0][0] * b.a[0][0] + a[0][1] * b.a[1][0]) % mod; c.a[0][1] = (a[0][0] * b.a[0][1] + a[0][1] * b.a[1][1]) % mod; c.a[1][0] = (a[1][0] * b.a[0][0] + a[1][1] * b.a[1][0]) % mod; c.a[1][1] = (a[1][1] * b.a[1][1] + a[1][0] * b.a[0][1]) % mod; return c; } inline Matrix operator +(Matrix b) const { register Matrix c; c.a[0][0] = a[0][0] + b.a[0][0] >= mod ? a[0][0] + b.a[0][0] - mod : a[0][0] + b.a[0][0]; c.a[0][1] = a[0][1] + b.a[0][1] >= mod ? a[0][1] + b.a[0][1] - mod : a[0][1] + b.a[0][1]; c.a[1][0] = a[1][0] + b.a[1][0] >= mod ? a[1][0] + b.a[1][0] - mod : a[1][0] + b.a[1][0]; c.a[1][1] = a[1][1] + b.a[1][1] >= mod ? a[1][1] + b.a[1][1] - mod : a[1][1] + b.a[1][1]; return c; } inline Matrix operator *(int b) const { register Matrix c; c.a[0][0] = a[0][0] * b % mod, c.a[0][1] = a[0][1] * b % mod; c.a[1][0] = a[1][0] * b % mod, c.a[1][1] = a[1][1] * b % mod; return c; }} ANS, I, A;inline int Pow(ll x, int y) { register ll ret = 1; for (; y; y >>= 1, x = x * x % mod) if (y & 1) ret = ret * x % mod; return ret;}inline void Getroot() { register int i, phi = mod - 1, x = phi, j; for (cnt = 0, i = 2; i * i <= x; ++i) if (x % i == 0) { pr[++cnt] = i; while (x % i == 0) x /= i; } if (x > 1) pr[++cnt] = x; for (i = 2; i <= phi; ++i) { for (x = 0, j = 1; !x && j <= cnt; ++j) if (Pow(i, phi / pr[j]) == 1) x = 1; if (x) continue; g = i; return; }}inline Matrix MatrixPow(Matrix X, ll y) { register Matrix RET = I; for (; y; y >>= 1, X = X * X) if (y & 1) RET = RET * X; return RET;}int main() { register int i, w, wn; I.a[0][0] = I.a[1][1] = A.a[0][0] = A.a[0][1] = A.a[1][0] = 1; scanf("%d", &test); while (test) { scanf("%lld%d%d", &n, &k, &mod), --test; ANS = Matrix(), Getroot(), wn = 1, w = Pow(g, (mod - 1) / k); for (i = 0; i < k; ++i) ANS = ANS + MatrixPow(A * wn + I, n), wn = (ll)wn * w % mod; ANS = ANS * Pow(k, mod - 2); printf("%lld\n", ANS.a[0][0]); } return 0;}

转载于:https://www.cnblogs.com/cjoieryl/p/10187468.html

你可能感兴趣的文章
Android第三方开源FloatingActionButton(com.getbase.floatingactionbutton): FloatingActionsMenu【3】...
查看>>
Python3—— collections模块
查看>>
WebLogic口令猜解工具【Python脚本】
查看>>
消灭天价药,AI会成为新药神吗?
查看>>
[增删改查] 使用 React 做后台管理 CRUD 界面和 RESTful 交互
查看>>
springboot demo
查看>>
PostgreSQL 11 preview - 多阶段并行聚合array_agg, string_agg
查看>>
Linux下查看某个命令的参数
查看>>
CUBA Platform 7.0.4 发布,企业级应用开发平台
查看>>
KodExplorer 4.40 发布,权限机制优化
查看>>
北京软件造价评估联盟:开启软件成本度量新篇章
查看>>
Mac下安装eclipse(Mac 10.12/JDK/tomcat)
查看>>
Facebook的Aquila遭查,无人机事故频发皆因哪般
查看>>
4月17日云栖精选夜读:90后剁手党占了一半!天猫是如何成为奢侈品第一平台的?...
查看>>
后引力波之战已经打响,神秘伽马射线来自何方?
查看>>
ecshop运行超过30秒超时的限制解决办法
查看>>
家庭物联网:从全屋智能到数据服务
查看>>
大数据领域的新面孔!我国微生物大数据平台获得国家大力支持
查看>>
韩国现代汽车成功实测L4级别自动驾驶汽车
查看>>
大数据产业发展明确四大重点
查看>>