BZOJ 2219 数论之神

发表于 | 分类于

Description

给定正整数\(A, B, K\)，求方程\(x^A\equiv B\pmod{2K+1}\)在\(Z_{2K+1}\)中的解的数量。

多组询问，组数不超过1000。

\(1\leq A,B\leq 10^9, 1\leq K\leq 5\times 10^8\)。

Solution

算是一道综合性很强的题目了……

先将模数分解质因数，对每一种质因子幂单独求解，然后由中国剩余定理之类的知道最后的解数是每一部分的解数的积。

那么考虑每种质因子的幂\(m=p^c\)的解数。先把\(B\)取下模，要是\(B = 0\)的话那么\(x\)是\(p^{\lceil\frac{c}{A}\rceil}\)的倍数就行了，因此答案是\(p^{c-\lceil\frac{c}{A}\rceil}\)。

那么考虑\(B\neq 0\)的情况。你可以对两边取指标（因为\(m\)事奇质数的幂，所以一定有原根，原根当然很好找了），就得到了\(A\mathrm{ind}x\equiv\mathrm{ind}B\pmod{\varphi(m)}\)，然后根据线性同余方程那套理论，假设\(d=\gcd(A,\varphi(m))\)，如果\(d|\mathrm{ind}B\)才有解，且解的数量为\(d\)。

但是有一个问题就是可能有\(\gcd(B, m)>1\)，这样的话\(\mathrm{ind}B\)就不一定存在。那么我们在原始的式子里，把两边疯狂除\(p\)直到右边没有\(p\)了，那么假设\(B=p^su\)，那么可知新的方程为\(\frac{x^A}{p^s}\equiv u\pmod{p^{c - s}}\)。可以发现如果\(A\)不是\(s\)的因子的话那么凉了，要是是的话假设\(s = At\)，那么新的方程也可以写作\((\frac{x}{p^t})^A\equiv u\pmod{p^{c - s}}\)，这样的话就得到了一个新的方程，且\(u\)和现在的模数是互质的，因此可以直接用上面的方法做。需要注意的是我们现在只是求出了模\(p^{c-s}\)的解数，对应到\(p^c\)要乘上\(p^s\)。但是解肯定要乘上一个\(p^t\)，因此解数还要在除一下\(p^t\)。

Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cassert>
#include <cctype>
#include <cmath>
#include <algorithm>
#include <functional>
#include <utility>
#include <vector>
#include <map>
#include <tr1/unordered_map>
typedef long long ll;
ll gcd(ll a, ll b) {
  if(!b) {
    return a;
  } else {
    return gcd(b, a % b);
  }
}
void exgcd(ll a, ll b, ll &d, ll &x, ll &y) {
  if(!b) {
    d = a; x = 1; y = 0;
  } else {
    ll nx, ny; exgcd(b, a % b, d, nx, ny);
    x = ny; y = nx - ny * (a / b);
  }
}
ll inv(ll v, ll p) {
  ll d, x, y; exgcd(v, p, d, x, y);
  return (x + p) % p;
}
ll pow_mod(ll a, ll b, ll p) {
  ll ans = 1, res = a;
  while(b) {
    if(1LL & b) ans = ans * res % p;
    res = res * res % p; b >>= 1;
  }
  return ans;
}

void desc(ll x, std::vector<int> &V) {
  ll bd = floor(sqrt((double)x) + 0.5);
  for(ll i = 2; i <= bd; i ++) {
    if(x % i == 0LL) {
      V.push_back(i);
      while(x % i == 0LL) x /= i;
    }
  }
  if(x > 1LL) V.push_back(int(x));
}
ll find_g(ll p, ll phi, std::vector<int> &V) {
  if(p == 2LL) return 1;
  for(ll i = 2; i < p; i ++) {
    bool ok = true;
    for(int j = 0; j < V.size(); j ++) {
      ll c = phi / (ll)V[j];
      if(pow_mod(i, c, p) == 1LL) {
        ok = false; break;
      }
    }
    if(ok) return i;
  }
  return -1LL;
}
ll bsgs(ll a, ll b, ll p) {
  if(a == 0LL) {
    return (b == 0LL) ? 1LL : -1LL;
  }
  
  int bd = ceil(sqrt((double)p));
  std::tr1::unordered_map<ll, ll> ma;
  ll blk = 1LL;
  for(int i = 0; i < bd; i ++) {
    if(!ma.count(blk)) ma[blk] = i;
    blk = blk * a % p;
  }
  ll bs = 1LL; blk = inv(blk, p);
  for(int i = 0; (ll)i * (ll)bd < p; i ++) {
    ll th = b * bs % p;
    if(ma.count(th)) return ma[th] + i * bd;
    bs = bs * blk % p;
  }
  assert(false);
  return -1;
}

ll equation(ll a, ll b, ll p) {
  ll ret = gcd(a, p);
  if(b % ret != 0LL) return 0;
  else return ret;
}
ll solve(ll a, ll b, ll mod, ll p, ll c) {
  if(b == 0LL) {
    ll ret = ceil((double(c)) / (double(a)));
    return pow_mod(p, c - ret, 2000000007LL);
  } else {
    if(b % p == 0LL) {
      ll nc = 0;
      while(b % p == 0LL) nc ++, b /= p, mod /= p;
      if(nc % a != 0LL) return 0LL;
      ll t = nc / a;
      return solve(a, b, mod, p, c - nc) * pow_mod(p, (a - 1LL) * t, 2000000007LL);
    } else {
      ll phi = mod / p * (p - 1LL);
      std::vector<int> d_phi; desc(phi, d_phi);
      ll g = find_g(mod, phi, d_phi);
      return equation(a, bsgs(g, b, mod), phi);
    }
  }
}

int main() {
  int T; scanf("%d", &T);
  while(T --) {
    int a, b, k; scanf("%d%d%d", &a, &b, &k);
    k = 2 * k + 1;
    int bd = floor(sqrt((double)k) + 0.5);
    int ans = 1;
    for(int i = 2; i <= bd; i ++) {
      if(k % i == 0) {
        int c = 0, mod = 1;
        while(k % i == 0) k /= i, c ++, mod *= i;
        ans *= solve(a, b % mod, mod, i, c);
      }
    }
    if(k > 1) ans *= solve(a, b % k, k, k, 1);
    printf("%d\n", ans);
  }
  return 0;
}