bzoj5104 Fib数列

bzoj5104 Fib数列

求在模 \(10^9+9\) 的意义下 \(n\) 在斐波那契数列中出现的最小位置，无解输出 \(-1\)

\(n<10^9+9\)

---

由于 \(\sqrt5\) 在模 \(10^9+9\) 的意义下的二次剩余为 \(383008016\) ，因此可以将 \(\sqrt5\) 带入斐波那契数列的通项公式

令 \(y=\frac{1+\sqrt5}{2}\)

即求最小的 \(x\) 满足 $$\frac{y^{x-(\frac{-1}{y})}x}{\sqrt5}\equiv n\pmod {10^9+9}$$

化一下式子，得到 $$y^x+(-1)\frac{1}{y^x}\equiv\sqrt5n\pmod {10^9+9}$$

令 \(t=y^x\) ，原式可以看作为一个二元一次方程 \(t^2-\sqrt5nt+(-1)^{n+1}=0\)

于是 分类讨论 \(n\) 的奇偶，暴力带入求根公式， \(\Delta\) 用二次剩余计算，即可求出 \(t\)

现在的问题就变成了求 \(y^x\equiv t\pmod{10^9+9}\) ，BSGS 即可

时间复杂度 \(O(\sqrt n\log n)\)

代码

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

const int P = 1e9 + 9, inv_2 = 500000005, sqrt_5 = 383008016;
int N, omg;

inline int inc(int x, int y) {
  return x + y < P ? x + y : x + y - P;
}

struct domain {
  int x, y;

  domain() {}
  domain(int _x, int _y) :
    x(_x), y(_y) {}

  inline domain operator + (const domain &o) const {
    return domain(inc(x, o.x), inc(y, o.y));
  }

  inline domain operator * (const domain &o) const {
    return domain((1ll * x * o.x + 1ll * y * o.y % P * omg) % P, (1ll * x * o.y + 1ll * y * o.x) % P);
  }
};

inline int rnd() {
  return (rand() << 15 | rand()) % P;
}

inline int qp(int a, int k) {
  int res = 1;
  for (; k; k >>= 1, a = 1ll * a * a % P) {
    if (k & 1) res = 1ll * res * a % P;
  }
  return res;
}

inline domain qp(domain a, int k) {
  domain res = domain(1, 0);
  for (; k; k >>= 1, a = a * a) {
    if (k & 1) res = res * a;
  }
  return res;
}

int Cipolla(int n) {
  int k = (P - 1) >> 1;
  if (qp(n, k) == -1) {
    return -1;
  }
  int a = rnd();
  while (qp((1ll * a * a - n + P) % P, k) == 1) {
    a = rnd();
  }
  omg = (1ll * a * a - n + P) % P;
  return qp(domain(a, 1), k + 1).x;
}

int BSGS(int a, int b) {
  if (!a && b) return -1;
  map <int, int> s;
  int sz = sqrt(P), inv_a = qp(a, P - 2);
  for (int i = 0, cur = 1; i <= sz; i++) {
    s.insert(make_pair(1ll * b * cur % P, i)), cur = 1ll * cur * inv_a % P;
  }
  int pw = qp(a, sz);
  map <int, int> :: iterator it;
  for (int i = 0, cur = 1; i <= sz; i++, cur = 1ll * cur * pw % P) {
    if ((it = s.find(cur)) != s.end()) {
      return i * sz + (it -> second);
    }
  }
  return -1;
}

signed main() {
  srand(time(0));
  scanf("%d", &N);
  int x = 1ll * N * sqrt_5 % P;
  int y = 1ll * (sqrt_5 + 1) * inv_2 % P;
  int delta, tmp, res, ans = INT_MAX;
  delta = Cipolla((1ll * x * x - 4 + P) % P);
  if (~delta) {
    tmp = 1ll * (x + delta) * inv_2 % P, res = BSGS(y, tmp);
    if (~res && res & 1) ans = min(ans, res);
    tmp = 1ll * (x - delta + P) * inv_2 % P, res = BSGS(y, tmp);
    if (~res && res & 1) ans = min(ans, res);
  }
  delta = Cipolla((1ll * x * x + 4) % P);
  if (~delta) {
    tmp = 1ll * (x + delta) * inv_2 % P, res = BSGS(y, tmp);
    if (~res && ~res & 1) ans = min(ans, res);
    tmp = 1ll * (x - delta + P) * inv_2 % P, res = BSGS(y, tmp);
    if (~res && ~res & 1) ans = min(ans, res);
  }
  printf("%d", ans < INT_MAX ? ans : -1);
  return 0;
}