首页 > 代码库 > HDU 5895 Mathematician QSC(矩阵乘法+循环节降幂+除法取模小技巧+快速幂)
HDU 5895 Mathematician QSC(矩阵乘法+循环节降幂+除法取模小技巧+快速幂)
传送门:HDU 5895 Mathematician QSC
这是一篇很好的题解,我想讲的他基本都讲了http://blog.csdn.net/queuelovestack/article/details/52577212
【分析】
一开始想简单了,对于a^x mod p这种形式的直接用欧拉定理的数论定理降幂了
结果可想而知,肯定错,因为题目并没有保证gcd(x,s+1)=1,而欧拉定理的数论定理是明确规定的
所以得另谋出路
那么网上提供了一种指数循环节降幂的方法
具体证明可以自行从网上找一找
有了这种降幂的方法之后,我们要分析一下如何求g(n)
由于f(0)=0,f(1)=1,f(n)=f(n−2)+2∗f(n−1)(n≥2)
可得,g(n)=f(n)*f(n+1)/2
这个是很好发现的
如果你发现不了的话,可以直接丢到OEIS里搜一下
然后,要求出g(n*y),就需要先求出f(n*y)和f(n*y+1)
这时,我们可以考虑用矩阵乘法
构造矩阵
套一下矩阵快速幂的模板就可以求出f(n*y)和f(n*y+1)
然后要求g(n)还有个除以2的操作,显然除法取模要用逆元
但考虑到2与模数不一定互质,无法用乘法逆元,所以要采用一点小技巧转化一下
这样我们就可以得到简化好的最终的指数部分
这样我们用快速幂就可以求x的幂次对(s+1)取模了
【时间复杂度&&优化】
O(1ogn)
/************************************************************** Problem:hdu 5895 Mathematician QSC User: youmi Language: C++ Result: Accepted Time:31MS Memory:1584K****************************************************************///#pragma comment(linker, "/STACK:1024000000,1024000000")//#include<bits/stdc++.h>#include <iostream>#include <cstdio>#include <cstring>#include <algorithm>#include <map>#include <stack>#include <set>#include <sstream>#include <cmath>#include <queue>#include <deque>#include <string>#include <vector>#define zeros(a) memset(a,0,sizeof(a))#define ones(a) memset(a,-1,sizeof(a))#define sc(a) scanf("%d",&a)#define sc2(a,b) scanf("%d%d",&a,&b)#define sc3(a,b,c) scanf("%d%d%d",&a,&b,&c)#define scs(a) scanf("%s",a)#define sclld(a) scanf("%I64d",&a)#define pt(a) printf("%d\n",a)#define ptlld(a) printf("%I64d\n",a)#define rep(i,from,to) for(int i=from;i<=to;i++)#define irep(i,to,from) for(int i=to;i>=from;i--)#define Max(a,b) ((a)>(b)?(a):(b))#define Min(a,b) ((a)<(b)?(a):(b))#define lson (step<<1)#define rson (lson+1)#define eps 1e-6#define oo 0x3fffffff#define TEST cout<<"*************************"<<endlconst double pi=4*atan(1.0);using namespace std;typedef long long ll;template <class T> inline void read(T &n){ char c; int flag = 1; for (c = getchar(); !(c >= ‘0‘ && c <= ‘9‘ || c == ‘-‘); c = getchar()); if (c == ‘-‘) flag = -1, n = 0; else n = c - ‘0‘; for (c = getchar(); c >= ‘0‘ && c <= ‘9‘; c = getchar()) n = n * 10 + c - ‘0‘; n *= flag;}ll Pow(ll base, ll n, ll mo){ ll res=1; while(n) { if(n&1) res=res*base%mo; n>>=1; base=base*base%mo; } return res;}//***************************ll n,y,x,s;const ll mod=1000000007;ll modp,modq;const int maxn=2;ll euler(ll nn){ ll res=nn,a=nn; for(ll i=2;i*i<=a;i++){ if(a%i==0){ res=res/i*(i-1);//先进行除法是为了防止中间数据的溢出 while(a%i==0) a/=i; } } if(a>1) res=res/a*(a-1); return res;}struct matrix{ ll mat[maxn][maxn]; matrix operator*(const matrix & rhs)const { matrix ans; rep(i,0,maxn-1) rep(j,0,maxn-1) ans.mat[i][j]=0; rep(i,0,maxn-1) rep(j,0,maxn-1) rep(k,0,maxn-1) ans.mat[i][j]=(ans.mat[i][j]+mat[i][k]*rhs.mat[k][j])%modp; return ans; } matrix operator^(ll k)const { matrix rhs=*this; matrix res; rep(i,0,maxn-1) rep(j,0,maxn-1) res.mat[i][j]=(i==j); while(k) { if(k&1) res=res*rhs; rhs=rhs*rhs; k>>=1; } return res; }}xx;int main(){ #ifndef ONLINE_JUDGE freopen("in.txt","r",stdin); #endif int T_T; scanf("%d",&T_T); for(int kase=1;kase<=T_T;kase++) { read(n),read(y),read(x),read(s); modp=euler(s+1)*2; modq=s+1; xx.mat[0][0]=2,xx.mat[0][1]=1,xx.mat[1][0]=1,xx.mat[1][1]=0; matrix temp=xx^(n*y); ll fn1=temp.mat[0][0]; ll fn=temp.mat[1][0]; ll gn=fn*fn1%modp/2; ll ans=Pow(x,gn,modq); ptlld(ans); } return 0;}
HDU 5895 Mathematician QSC(矩阵乘法+循环节降幂+除法取模小技巧+快速幂)
声明:以上内容来自用户投稿及互联网公开渠道收集整理发布,本网站不拥有所有权,未作人工编辑处理,也不承担相关法律责任,若内容有误或涉及侵权可进行投诉: 投诉/举报 工作人员会在5个工作日内联系你,一经查实,本站将立刻删除涉嫌侵权内容。