// f(n) = a * f(n - 1) + b * f(n - 2);
// f(1) = c, f(2) = d // 可忽略
// 求f(n)mod p循环节长度
//c = a * a + 4b是模p的二次剩余时,枚举n = p - 1的因子
// 否则 枚举n=(p+1)(p-1)的因子
// 这里 p = 1000000007
#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
using namespace std;
#define LL long long
const int M = 2;
const LL mod = 1000000007;
LL fac[2][505];
int cnt, ct;
LL pri[6] = {2, 3, 7, 109, 167, 500000003};
LL num[6] = {4, 2, 1, 2, 1, 1};
struct mat {
LL a[M][M];
};
mat A, I = {1, 0, 0, 1};
mat multi(mat a, mat b) {
mat ret;
for(int i = 0; i < M; ++i)
for(int j = 0; j < M; ++j) {
ret.a[i][j] = 0;
for(int k = 0; k < M; ++k)
ret.a[i][j] = (ret.a[i][j] + a.a[i][k] * b.a[k][j]) % mod;
}
return ret;
}
mat qpow(mat a, LL k) {
mat ret = I;
while(k) {
if(k & 1)
ret = multi(ret, a);
a = multi(a, a), k >>= 1;
}
return ret;
}
LL qpow(LL a, LL k) {
LL ret = 1;
a %= mod;
while(k) {
if(k & 1)
ret = ret * a % mod;
a = a * a % mod;
k >>= 1;
}
return ret;
}
LL legendre(LL a, LL p) {
LL t = qpow(a, (p - 1) >> 1);
if(t == 1) return 1;
return -1;
}
void dfs(int dep, LL product = 1) {
if(dep == cnt) {
fac[1][ct++] = product;
return;
}
for(int i = 0; i <= num[dep]; ++i) {
dfs(dep + 1, product);
product *= pri[dep];
}
}
bool ok(mat a, LL n) {
mat ans = qpow(a, n);
return (ans.a[0][0] == 1 && ans.a[0][1] == 0 && ans.a[1][0] == 0 && ans.a[1][1] == 1);
}
int main() {
fac[0][0] = 1, fac[0][1] = 2;
fac[0][2] = 500000003;
fac[0][3] = 1000000006;
LL a, b, c, d;
while(scanf("%lld%lld%lld%lld", &a, &b, &c, &d) != EOF) {
LL t = a * a + 4 * b;
A.a[0][0] = a, A.a[0][1] = b;
A.a[1][0] = 1, A.a[1][1] = 0;
if(legendre(t, mod) == 1) {
for(int i = 0; i < 4; ++i) {
if(ok(A, fac[0][i])) {
printf("%lld\n", fac[0][i]);
break;
}
}
}
else {
ct = 0, cnt = 6;
dfs(0, 1);
sort(fac[1], fac[1] + ct);
for(int i = 0; i < ct; ++i) {
if(ok(A, fac[1][i])) {
printf("%lld\n", fac[1][i]);
break;
}
}
}
}
return 0;
}
(模板)广义fib循环节
最新推荐文章于 2021-02-16 22:11:33 发布