hdu5780 gcd

hdu5780 gcd

给定 (n, x) ，求 $$displaystylesum_{i=1}^nsum_{j=1}ngcd(x^{i-1, x}j-1)pmod{10^9+7}$$

共 (T) 组询问

(Tleq300, n, xleq10^7)

数论，数论分块

---

首先有一个公式： (gcd(x^a-1, x^b-1)=x^{gcd(a, b)}-1)

于是

[egin{aligned}displaystylesum_{i=1}^nsum_{j=1}^ngcd(x^i-1, x^j-1)&=displaystylesum_{i=1}^nsum_{j=1}^nx^{gcd(i, j)-1}\&=displaystylesum_{i=1}^nsum_{j=1}^nsum_{d=1}^n{x^{d[gcd(i, j)=d]}-1}\&=displaystylesum_{i=1}^nsum_{j=1}^nsum_{d=1}^n{[gcd(i, j)=d](x^d-1)}\&=displaystylesum_{d=1}^n(x^d-1)sum_{i=1}^{frac{n}{d}}sum_{j=1}^{frac{n}{d}}{[gcd(i, j)=1]}\&=displaystylesum_{d=1}^n(x^d-1) imes(2 imes(sum_{i=1}^{frac{n}{d}}varphi(i))-1)end{aligned} ]

注意到第二项的 (frac{n}{d}) ，可以数论分块，第一项直接等比数列求和

需预处理 (displaystylesum_{i=1}^n{varphi(i)}pmod{10^9+7})

时间复杂度 (O(n+Tsqrt nlog n))

代码

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

const int maxn = 1e7 + 10, P = 1e9 + 7;
int tot, p[maxn], phi[maxn];

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

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

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;
}

void sieve() {
  int N = 1e7;
  for (int i = 2; i <= N; i++) {
    if (!p[i]) p[++tot] = i, phi[i] = i - 1;
    for (int j = 1; j <= tot && i * p[j] <= N; j++) {
      p[i * p[j]] = 1;
      if (i % p[j] == 0) {
        phi[i * p[j]] = phi[i] * p[j]; break;
      }
      phi[i * p[j]] = phi[i] * (p[j] - 1);
    }
  }
  phi[1] = 1;
  for (int i = 2; i <= N; i++) {
    phi[i] = inc(phi[i], phi[i - 1]);
  }
}

inline int calc(int x, int k) {
  if (x == 1) return k + 1;
  return 1ll * dec(qp(x, k + 1), 1) * qp(x - 1, P - 2) % P;
}

inline void solve() {
  int n, x, ans = 0;
  scanf("%d %d", &x, &n);
  for (int l = 1, r; l <= n; l = r + 1) {
    r = n / (n / l);
    ans = (ans + 1ll * dec(dec(calc(x, r), calc(x, l - 1)), r - l + 1) * dec(inc(phi[n / l], phi[n / l]), 1)) % P;
  }
  printf("%d
", ans);
}

int main() {
  sieve();
  int Tests;
  scanf("%d", &Tests);
  while (Tests--) solve();
  return 0;
}