(模板)广义fib循环节

来源:互联网 发布:若人欲了知三世一切佛 编辑:程序博客网 时间:2024/05/13 12:23
// 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 longconst 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;}

0 0