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

2019-04-14 17:43发布

站内文章 / 模拟电子
15473 0 1780

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 using namespace std; typedef long long ll; ll A, C, P, Ans, test, Yn, Y; map 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 ", Ans); } }

热门文章