luoguP4000 斐波那契数列
题目链接
题解
根据这个东西
https://www.cnblogs.com/sssy/p/9418732.html
我们可以找出%p意义下的循环节
然后就可以做了
人傻，自带，大，常数
代码
#include<bits/stdc++.h>
using namespace std;
#define LL long long
const LL maxn = 1000007;
LL dp[maxn * 10];
LL prime[maxn],s = 0;
bool vis[maxn];
void init_prime() {
for(LL i = 2;i < maxn;i ++) {
if(!vis[i]) prime[s ++] = i;
for(LL j = 0;j < s && i * prime[j] < maxn;j ++) {
vis[i * prime[j]] = 1;
if(i%prime[j] == 0) break;
}
}
}
LL pow_mod(LL a1,LL b1){
LL ret = 1;
for(;b1;b1 >>= 1,a1 *= a1)
if(b1 & 1) ret = ret * a1;
return ret;
}
LL pow_mod2(LL a,LL b,LL p) {
LL ret = 1;
for(;b;b >>= 1,a = a * a % p)
if(b & 1) ret = ret * a % p;
return ret;
}
LL gcd(LL a,LL b) { return b ? gcd(b,a % b) : a; }
bool f(LL n,LL p) {
if(pow_mod2(n,(p - 1) >> 1,p) != 1) return false;
return true;
}
struct matrix{
LL x1,x2,x3,x4;
};
matrix matrix_a,matrix_b;
matrix M2(matrix aa,matrix bb,LL mod){
matrix tmp;
tmp.x1 = (aa.x1 * bb.x1 % mod + aa.x2 * bb.x3 % mod) % mod;
tmp.x2 = (aa.x1 * bb.x2 % mod + aa.x2 * bb.x4 % mod) % mod;
tmp.x3 = (aa.x3 * bb.x1 % mod + aa.x4 * bb.x3 % mod) % mod;
tmp.x4 = (aa.x3 * bb.x2 % mod + aa.x4 * bb.x4 % mod) % mod;
return tmp;
}
matrix M(LL n,LL mod){
matrix a,b;
a=matrix_a;b=matrix_b;
for(;n;n >>= 1,a = M2(a,a,mod))
if(n & 1) b = M2(b,a,mod);
return b;
}
LL fac[100][2],l,x,fs[1000];
void dfs(LL count,LL step) {
if(step == l) {
fs[x ++] = count;
return ;
}
LL sum = 1;
for(LL i = 0;i < fac[step][1];++ i) {
sum *= fac[step][0];
dfs(count * sum,step + 1);
}
dfs(count,step + 1);
}
LL solve2(LL p){
if(p < 1e6 && dp[p]) return dp[p];
bool ok = f(5,p);//判断5是否为p的二次剩余
LL t;
if(ok) t = p - 1;
else t = 2 * p + 2;
l = 0;
for(LL i = 0;i < s;i ++){
if(prime[i] > t / prime[i]) break;
if(t % prime[i] == 0) {
LL count = 0;
fac[l][0] = prime[i];
while(t % prime[i] == 0) count ++ , t /= prime[i];
fac[l ++][1] = count;
}
}
if(t > 1) fac[l][0]=t, fac[l++][1]=1;
x = 0;
dfs(1 , 0); //求t的因子
sort(fs,fs + x);
for(LL i = 0;i < x;i ++) {
matrix m1 = M(fs[i],p);
if(m1.x1 == m1.x4 && m1.x1 == 1 && m1.x2 == m1.x3 && m1.x2 == 0) {
if(p < 1e6) dp[p] = fs[i];
return fs[i];
}
}
}
LL get_M(LL n){
LL ans = 1,cnt;
for(LL i = 0;i < s;i ++) {
if(prime[i] > n / prime[i]) break;
if(n % prime[i] == 0) {
LL count = 0;
while(n % prime[i] == 0) count ++,n /= prime[i];
cnt = pow_mod(prime[i],count - 1);
cnt *= solve2(prime[i]);
ans = (ans / gcd(ans , cnt)) * cnt;
}
}
if(n > 1){
cnt = 1;
cnt *= solve2(n);
ans = ans / gcd(ans,cnt) * cnt;
}
return ans;
}
char Ss[30000007];
struct bign {
int z[30000007],l;
void init() {
memset(z,0,sizeof(z));
scanf("%s",Ss + 1);
l = strlen(Ss + 1);
for(int i = 1;i <= l;i ++)
z[i] = Ss[l - i + 1] - '0';
}
LL operator % (const long long & a) const {
LL b = 0;
for (int i = l;i >= 1;i --)
b = (b * 10 + z[i]) % a;
return b;
}
}z;
LL m1;
struct Matrix {
LL a[3][3];
Matrix () { memset(a,0,sizeof a); }
Matrix operator * (const Matrix & p) const {
Matrix ret;
for(int i = 0;i <= 1;++ i)
for(int j = 0;j <= 1;++ j)
for(int k = 0;k <= 1;++ k)
ret.a[i][j] = (ret.a[i][j] + ( a[i][k] * p.a[k][j] ) % m1 ) % m1;
return ret;
}
};
LL solve(LL x) {
Matrix p,q;
p.a[0][0] = 1; p.a[0][1] = 1; p.a[1][0] = 1;
q.a[0][1] = 1;
for(;x;x >>= 1,p = p * p)
if(x & 1) q = q * p;
return q.a[0][0];
}
main() {
LL t,n;
init_prime();
matrix_a.x1 = matrix_a.x2 = matrix_a.x3 = 1;
matrix_a.x4 = 0;
matrix_b.x1 = matrix_b.x4 = 1;
matrix_b.x2 = matrix_b.x3 = 0;
dp[2] = 3;dp[3] = 8;dp[5] = 20;
z.init(); scanf("%lld",&n); m1 = n;
n = get_M(n);
//printf("%lld\n",n % m1);
n = z % n;
printf("%lld\n",solve(n));
return 0;
}
