[Tsinsen A1300]JZPKIL

danihao123 posted @ 2018年5月29日 16:16 in 题解 with tags 莫比乌斯反演 Tsinsen 狄利克雷卷积 伯努利数 Pollard's rho Miller-Rabin , 1062 阅读
转载请注明出处:http://danihao123.is-programmer.com/

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

卡常记录(三分之一)

下面开始颓式子:

\[
\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;
}
AP 2nd Inter Economi 说:
Sep 08, 2022 10:07:02 PM

The AP Intermediate students can download the Economics question bank with solved study material with practice questions in chapter wise to every TM, EM, UM student, and the economics subject paper-1 AP 2nd Inter Economics Model Paper and paper-2 important questions with suggestions are available through AP Jr and Sr inter Economics Model Paper 2023 Pdf with guess paper. The AP Intermediate students can download the Economics question bank with solved study material with practice questions in chapter wise to every TM, EM, UM student


登录 *


loading captcha image...
(输入验证码)
or Ctrl+Enter