POJ1845 Sumdiv(求所有因数和+矩阵快速幂)
2016-02-02 20:16
531 查看
题目问$A^B$的所有因数和。
根据唯一分解定理将A进行因式分解可得:A = p1^a1 * p2^a2 * p3^a3 * pn^an.
A^B=p1^(a1*B)*p2^(a2*B)*...*pn^(an*B);
A^B的所有约数之和sum=[1+p1+p1^2+...+p1^(a1*B)]*[1+p2+p2^2+...+p2^(a2*B)]*[1+pn+pn^2+...+pn^(an*B)]
知道这个,问题就变成求出A的所有质因数pi以及个数n,然后$\prod(1+p_i+p_i^2+\cdots+p_i^{n-1}+p_i^n)$就行了。可以构造矩阵来求:
记$S_n=p_i+p_i^2+\cdots+p_i^{n-1}+p_i^n$
$$ \begin{bmatrix} p_i & 1 \\ 0 & 1 \end{bmatrix} \times \begin{bmatrix} S_n \\ p_i \end{bmatrix} = \begin{bmatrix} S_{n+1} \\ p_i \end{bmatrix} $$
$$ \begin{bmatrix} S_n \\ p_i \end{bmatrix} = \begin{bmatrix} p_i & 1 \\ 0 & 1 \end{bmatrix} ^n \times \begin{bmatrix} S_0 \\ p_i \end{bmatrix} $$
A忘了$\pmod {9901}$,爆intWA到头疼= =
根据唯一分解定理将A进行因式分解可得:A = p1^a1 * p2^a2 * p3^a3 * pn^an.
A^B=p1^(a1*B)*p2^(a2*B)*...*pn^(an*B);
A^B的所有约数之和sum=[1+p1+p1^2+...+p1^(a1*B)]*[1+p2+p2^2+...+p2^(a2*B)]*[1+pn+pn^2+...+pn^(an*B)]
知道这个,问题就变成求出A的所有质因数pi以及个数n,然后$\prod(1+p_i+p_i^2+\cdots+p_i^{n-1}+p_i^n)$就行了。可以构造矩阵来求:
记$S_n=p_i+p_i^2+\cdots+p_i^{n-1}+p_i^n$
$$ \begin{bmatrix} p_i & 1 \\ 0 & 1 \end{bmatrix} \times \begin{bmatrix} S_n \\ p_i \end{bmatrix} = \begin{bmatrix} S_{n+1} \\ p_i \end{bmatrix} $$
$$ \begin{bmatrix} S_n \\ p_i \end{bmatrix} = \begin{bmatrix} p_i & 1 \\ 0 & 1 \end{bmatrix} ^n \times \begin{bmatrix} S_0 \\ p_i \end{bmatrix} $$
A忘了$\pmod {9901}$,爆intWA到头疼= =
#include<cstdio> #include<cstring> using namespace std; struct Mat{ int m[2][2]; }; Mat operator*(const Mat &m1,const Mat &m2){ Mat m={0}; for(int i=0; i<2; ++i){ for(int j=0; j<2; ++j){ for(int k=0; k<2; ++k){ m.m[i][j]+=m1.m[i][k]*m2.m[k][j]; m.m[i][j]%=9901; } } } return m; } int calu(int a,int n){ a%=9901; Mat e={1,0,0,1},x={a,1,0,1}; while(n){ if(n&1) e=e*x; x=x*x; n>>=1; } return (e.m[0][1]*a+1)%9901; } bool isPrime(int n){ if(n<2) return 0; for(int i=2; i*i<=n; ++i){ if(n%i==0) return 0; } return 1; } int main(){ int a,b; scanf("%d%d",&a,&b); if(isPrime(a)){ printf("%d",calu(a,b)); return 0; } int res=1; for(int i=2; i*i<=a; ++i){ if(a%i) continue; if(isPrime(i)){ int cnt=0,tmp=a; while(tmp%i==0){ ++cnt; tmp/=i; } res*=calu(i,cnt*b); res%=9901; } if(i!=a/i && isPrime(a/i)){ int cnt=0,tmp=a; while(tmp%i==0){ ++cnt; tmp/=i; } res*=calu(a/i,cnt*b); res%=9901; } } printf("%d",res); return 0; }
相关文章推荐
- oracle体系结构(图)
- yaourt: a pacman frontend(pacman前端,翻译)
- Robot Framework与Web界面自动化测试学习笔记:如何判断单选框的选中状态
- 在mac 上配置AndroidStudio碰到的坑
- 二叉排序树 2
- 在这里主要讨论脚本运行方式和“命令”运行的方式。环境变量和命令关系
- 看了第二个人的最小生成树的最短路径的差分约束系统
- git 解决冲突的办法
- 第一个人的解释:最小生成树中的最短路问题:差分约束系统
- 深度优先算法
- python 多线程笔记(3)-- 线程的私有命名空间
- Zabbix实现告警分级
- 什么是javaScript闭包
- 分治、CDQ分治小结(need to be updated)
- 如何将 Debian Linux 中的默认的 Python 版本切换为替代版本
- [Leecode] Maximum Gap
- JavaScript基础笔记集合
- 河南第四届ACM省赛(走迷宫)
- 如何高效率协同工作-。工具
- CF 17B Hierarchy