[数论][二次剩余][BSGS] CodeChef FN

Description

设Fn是斐波那契数列第n项，给定FnmodP和P求最小的n。

保证5在膜P意义下有二次剩余。

Solution

经过一系列xjb推会发现FnmodP是循环的。所以答案小于P。

首先用Cipolla算法求出5的二次剩余。

考虑Fn的通项公式。设x=5√，y=5√+１2就有Fn=1x[yn+(−y)−n]通过解一元二次方程，可以得到yn之后就BSGS求一下n，求的时候判下奇偶性就好啦。

#include <bits/stdc++.h>
using namespace std;

typedef long long ll;

ll A, C, P, Ans, test, Yn, Y;
map<ll, ll> ap;
ll w, rtf, RT, res;
struct Complex {
ll r, i;
Complex(ll _r = 0, ll _i = 0):r(_r), i(_i) {}
inline Complex operator +(ll x) {
return Complex((x + r) % P, i);
}
inline Complex operator *(Complex a) {
ll rr = (a.r * r % P + w * a.i % P * i % P + P * 3) % P,
ii = (a.i * r % P + a.r * i % P) % P;
return Complex(rr, ii);
}
};

inline ll Pow(Complex a, ll b) {
Complex c = 1;
while (b) {
if (b & 1) c = c * a;
a = a * a; b >>= 1;
}
return c.r;
}
inline ll Pow(ll a, ll b) {
ll c = 1;
while (b) {
if (b & 1) c = c * a % P;
a = a * a % P; b >>= 1;
}
return (c + P) % P;
}
inline ll Inv(ll a, ll P) {
return Pow(a, P - 2);
}
ll Log(ll a, ll b, ll p, ll par) {
ap.clear();
ll m = ceil(sqrt(p)), ia = Inv(Pow(a, m), P), res;
res = Inv(b, p); par ^= 1;
for (int i = 0; i < m; i++) {
if (!ap.count(res)) ap[res] = i;
res = res * a % P;
}
res = 1;
for (int i = 0; i <= m; i++) {
if (ap.count(res) && (((ap[res] + i * m) & 1) ^ par)) return ap[res] + i * m;
res = res * ia % P;
}
return -1;
}
inline ll Qr(ll x, ll p) {
if (Pow(x, (P - 1) / 2) == P - 1) return -1;
if (Pow(x, (P - 1) / 2) == 0) return 0;
ll a;
Complex k(0, 1);
while (true) {
a = rand();
if (Pow((a * a - x + P) % P, (P - 1) / 2) == P - 1) break;
}
w = (a * a - x + P) % P;
return Pow(k + a, (P + 1) / 2);
}
inline ll Min(ll a, ll b) {
return a < b ? a : b;
}

int main(void) {
freopen("1.in", "r", stdin);
freopen("1.out", "w", stdout);
scanf("%d", &test);
srand(19260817);
while (test--) {
scanf("%lld%lld", &C, &P);
rtf = Qr(5, P);
Ans = P;
A = C * rtf % P;
Y = (rtf + 1) * (P + 1) / 2 % P;

RT = Qr((A * A % P + 4) % P, P);
if (~RT) {
Yn = (A + RT) * (P + 1) / 2 % P;
res = Log(Y, Yn, P, 0);
if (~res) Ans = Min(Ans, res);
Yn = (A - RT + P) * (P + 1) / 2 % P;
res = Log(Y, Yn, P, 0);
if (~res) Ans = Min(Ans, res);
}

RT = Qr((A * A % P - 4 + P) % P, P);
if (~RT) {
Yn = (A + RT) * (P + 1) / 2 % P;
res = Log(Y, Yn, P, 1);
if (~res) Ans = Min(Ans, res);
Yn = (A - RT + P) * (P + 1) / 2 % P;
res = Log(Y, Yn, P, 1);
if (~res) Ans = Min(Ans, res);
}
if (Ans == P) Ans = -1;
printf("%lld\n", Ans);
}
}