卡常的题见过，这么变态的卡常题……可能出题人过于文明，，，下面是卡常记录的一部分：

下面开始颓式子：

\[
\begin{aligned}
&\quad\sum_{i = 1}^n (i, n)^x[i, n]^y\\
&=n^y\sum_{i = 1}^n (i, n)^{x - y} i^y\\
&=n^y\sum_{k | n} k^{x - y}\sum_{i = 1}^{\frac{n}{k}} (ik)^y [(i, \frac{n}{k})=1]\\
&=n^y\sum_{k | n} k^x\sum_{i = 1}^{\frac{n}{k}} i^y\sum_{d | i, d | \frac{n}{k}} \mu(d)\\
&=n^y\sum_{k | n} k^x\sum_{d | \frac{n}{k}}\mu(d) \sum_{i = 1}^{\frac{n}{kd}} (id)^y\\
&=n^y\sum_{k | n} k^x\sum_{d | \frac{n}{k}}\mu(d)d^y \sum_{i = 1}^{\frac{n}{kd}} i^y\\
&=n^y\cdot(\mathrm{id}^x\ast (\mu\cdot\mathrm{id}^y)\ast h_y)(n)
\end{aligned}
\]

其中\(h_y\)表示指数为\(y\)的等幂求和函数。

接下来的感觉就和BZOJ 3601差不多了……还是设想如果后面是积性函数那么肥肠好做。可惜不是，但\(h_y\)是一个\(y + 1\)次多项式，所以拆开每一项考虑。然后我们还是对每一类质因子分开算，利用\(\mu\cdot\mathrm{id}^y\)这个函数的性质（当整个卷积的参数是质数的幂\(p^k\)时，这个函数只有输入为\(1\)和\(p\)时结果不为0），就肥肠好做了。

但是这个题非常的坑……\(n\)超级大，所以质因数分解要用Pollard's rho……这样一来我们又需要写Miller-Rabin，再一来我们发现当\(n\)很大的时候可能会有一些取膜乘法的模数太大，导致一般的乘法会发生溢出，所以又要写快速乘……\(O(\log n)\)快速乘又会被卡，所以你要写gyz的\(O(1)\)快速乘……

还有一些卡常的建议：

1. Pollard's rho和Miller-Rabin一定要写的精细……复杂度首先不要像我这种沙茶一样上来多挂一个log，然后常数也要注意优化。
2. \(h_y\)说是\(y + 1\)次，其实有大约一半的项的系数都为0。
3. \(x = y\)时直接利用伯努利数算就行。
4. 模数小的地方就别xjb用快速乘了……

代码：

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cctype>
#include <ctime>
#include <algorithm>
#include <utility>
#include <random>
#include <vector>
using ll = long long;
using ld = long double;
const ll ha = 1000000007LL;
std::mt19937 gen;
inline ll mul(ll x, ll y, ll z = ha){
  // x %= z; y %= z;
  if(z <= ha) return (x * y) % z;
  ll ret = (x * y - (ll)(x / (ld)z * y + 1e-3) * z + z) % z;
  return ret;
}
inline ll pow_mod(ll a, ll b, ll p = ha) {
  ll ans = 1LL, res = a;
  while(b) {
    if(1LL & b) ans = mul(ans, res, p);
    res = mul(res, res, p);
    b >>= 1;
  }
  return ans;
}
void exgcd(ll a, ll b, ll &x, ll &y) {
  if(!b) {
    x = 1LL; y = 0LL;
  } else {
    ll nx, ny; exgcd(b, a % b, nx, ny);
    x = ny; y = nx - (a / b) * ny;
  }
}
inline ll inv(ll v, ll p = ha) {
  ll x, y;
  exgcd(v, p, x, y);
  return (x + p) % p;
}

inline bool witness(ll a, ll u, ll t, ll p) {
  ll x = pow_mod(a, u, p);
  if(x == 1LL) return true;
  while(t --) {
    if(x != 1LL && x != (p - 1LL) && mul(x, x, p) == 1LL) return false;
    x = mul(x, x, p);
  }
  return (x == 1LL);
}
inline bool miller_rabin(ll x) {
  static const int JP[9] = {2, 3, 5, 7, 11, 13, 17, 19, 23};
  if(x <= 23LL) {
    if(x == 2LL || x == 3LL || x == 5LL || x == 7LL || x == 11LL ||
       x == 13LL || x == 17LL || x == 19LL || x == 23LL) {
      return true;
    } else {
      return false;
    }
  }
  if(x % 2LL == 0LL) return false;
  ll u = x - 1LL, t = 0LL;
  while(!(u & 1LL)) {
    u >>= 1; t ++;
  }
  for(int i = 0; i < 9; i ++) {
    if(!witness(JP[i], u, t, x)) return false;
  }
  return true;
}
inline ll gcd(ll a, ll b) {
  while(b != 0LL) {
    ll na = b, nb = a % b;
    a = na, b = nb;
  }
  return a;
}
inline ll abs(ll x) {
  if(x < 0LL) {
    return -x;
  } else {
    return x;
  }
}
inline ll pollard_rho(ll p) {
  if(p % 2LL == 0LL) return 2LL;
  while(true) {
    ll a = (gen() % (p - 1LL) + 1LL), c = (gen() % (p - 1LL) + 1LL);
    ll b = a; if(c == 2LL) c = 3LL;
    do {
      a = (mul(a, a, p) + c) % p;
      b = (mul(b, b, p) + c) % p;
      b = (mul(b, b, p) + c) % p;
      if(a == b) break;
      ll d = gcd(p, abs(b - a));
      if(d > 1LL && d < p) return d;
    } while(true);
  }
}
void find_fac(ll x, std::vector<ll> &A) {
  if(miller_rabin(x)) {
    A.push_back(x); return;
  }
  ll p = pollard_rho(x);
  // while(p == x) p = pollard_rho(x);
  find_fac(p, A); find_fac(x / p, A);
}
using pii = std::pair<ll, ll>;
inline void fact(ll x, std::vector<pii> &V) {
  std::vector<ll> A; find_fac(x, A);
  std::sort(A.begin(), A.end());
  ll las = 1LL, t = 0LL;
  for(int i = 0; i < A.size(); i ++) {
    ll v = A[i];
    if(v != las) {
      if(las > 1LL) {
        V.push_back(std::make_pair(las, t));
      }
      las = v; t = 1LL;
    } else {
      t ++;
    }
  }
  if(las > 1LL) V.push_back(std::make_pair(las, t));
#ifdef LOCAL
  printf("Brk %lld...\n", x);
  for(auto p : V) {
    printf("Pair (%lld, %lld)\n", p.first, p.second);
  }
#endif
}

const int maxn = 3505;
ll C[maxn][maxn];
inline void process_C() {
  C[0][0] = 1LL;
  for(int i = 1; i <= 3005; i ++) {
    C[i][0] = C[i][i] = 1LL;
    for(int j = 1; j < i; j ++) {
      C[i][j] = (C[i - 1][j] + C[i - 1][j - 1]) % ha;
    }
  }
}
ll B[maxn];
inline void process_B() {
  B[0] = 1LL;
  for(int i = 1; i <= 3001; i ++) {
    B[i] = 1LL;
    for(int j = 0; j < i; j ++) {
      B[i] = (B[i] - (C[i][j] * ((B[j] * inv(i - j + 1)) % ha)) % ha + ha) % ha;
    }
  }
}

std::vector<pii> V;
ll n; int x, y;
inline ll calc(ll k, int z) {
  if(k == 0LL) return 0LL;
  ll ans = 1LL;
  for(auto pr : V) {
    ll p = pr.first, t = pr.second;
    p %= ha;
    ll v1 = 0LL;
    for(ll i = 0; i <= t; i ++) {
      v1 += (pow_mod(p, i * (ll)x) * pow_mod(p, (t - i) * (ll)z)) % ha;
      v1 %= ha;
    }
    ll v2 = 0LL;
    for(ll i = 0; i <= t - 1LL; i ++) {
      v2 += (pow_mod(p, i * (ll)x) * pow_mod(p, (t - i - 1LL) * (ll)z)) % ha;
      v2 %= ha;
    }
    v2 = (v2 * (ha - pow_mod(p, y))) % ha;
    v1 = (v1 + v2) % ha;
    ans = (ans * v1) % ha;
  }
  return (ans * k) % ha;
}
inline ll calc2() {
  n %= ha;
  ll ans = 0LL;
  for(int i = 0; i <= y; i ++) {
    ll delta = (C[y + 1][i] * B[i]) % ha;
    delta = (delta * pow_mod(n, y + 1 - i)) % ha;
    ans = (ans + delta) % ha;
  }
  ans = (ans * inv(y + 1)) % ha;
  ans = (ans * pow_mod(n % ha, y)) % ha;
  return ans;
}

int main() {
  gen.seed(time(0));
  double d = clock();
  process_C(); process_B();
#ifdef LOCAL
  double dd = clock();
  printf("Time used : %.3lf\n", (dd - d) / (double(CLOCKS_PER_SEC)));
#endif
  int T; scanf("%d", &T);
  while(T --) {
    scanf("%I64d%d%d", &n, &x, &y);
    if(x == y) {
      printf("%I64d\n", calc2());
      continue;
    }
    V.clear(); fact(n, V);
    ll ans = 0LL;
    for(int i = 0; i <= y; i ++) {
      ans = (ans + calc((((C[y + 1][i] * inv(y + 1)) % ha) * B[i]) % ha, y + 1 - i)) % ha;
    }
    ans = (ans * pow_mod(n % ha, y)) % ha;
    printf("%I64d\n", ans);
  }
  return 0;
}

• 相关文章
• 全局相关文章