51nod1965. 奇怪的式子（min_25筛）

题目链接

http://www.51nod.com/Challenge/Problem.html#!#problemId=1965

题解

需要求的式子显然是个二合一形式，我们将其拆开，分别计算 (prod_limits{i = 1}^n sigma_0(i)^i) 与 (prod_limits{i = 1}^n sigma_0(i)^{mu(i)})，再将两部分乘起来得到答案。

对于第一部分 (prod_limits{i = 1}^n sigma_0(i)^i)：

由于若 (x) 的唯一分解式为 (p_1^{alpha_1}p_2^{alpha_2} cdots p_k^{alpha_k})，那么有 (sigma_0(x) = prod_{i = 1}^k (alpha_i + 1))，因此我们可以考虑求每一个质数的贡献。令 (P) 表示质数集合，不难得到：

[prod_{i = 1}^n sigma_0(i)^i = prod_{p in P} prod_{p^k leq n} (k + 1)^w ]

其中的 (w) 即为所有唯一分解式中包含 (p^k) 的数的和，做个简单容斥可以得到 (w = p^k cdot s(leftlfloorfrac{n}{p^k} ight floor) - p^{k + 1} cdot s(leftlfloorfrac{n}{p^{k + 1}} ight floor))，(s(x)) 表示 (1 sim x) 的所有数的和，即 (frac{x(x + 1)}{2})。

当 (p leq sqrt n) 时，我们可以直接暴力枚举 (p) 与 (k) 来计算式子的值。

当 (p > sqrt n) 时，对应的 (k) 只可能为 (1)，因此式子可以化为 (2^{c}) 的形式，其中 (c = sum_limits{p in P, p > sqrt n} p leftlfloorfrac{n}{p} ight floor)。考虑如何求 (c)。我们对 (leftlfloorfrac{n}{p} ight floor) 整除分块后，问题转化为了计算一段区间内所有质数的和，即求质数的前缀和。使用 min_25 筛即可。

对于第二部分 (prod_limits{i = 1}^n sigma_0(i)^{mu(i)})：

由于当 (x) 的唯一分解式中存在一个 (p^{alpha}) 满足 (alpha > 1) 时，(mu(x) = 0)，因此对答案式子 (prod_limits{i = 1}^n sigma_0(i)^{mu(i)}) 有贡献的 (i) 一定满足 (i) 是若干个互不相同的质数的乘积。显然，此时 (sigma_0(i)) 又可以写成 (2) 的次幂的形式，且若 (i) 是 (k) 个互不相同的质数的乘积，那么 (sigma_0(i) = 2^k)。因此若令 (g(x)) 表示 (x) 包含的不同的质因数的个数，那么有：

[prod_limits{i = 1}^n sigma_0(i)^{mu(i)} = 2^w ]

其中 (w = sum_limits{i = 1}^n mu(i)g(i))。

令 (P_i) 表示从小到大第 (i) 个质数，({ m minp}(x)) 表示 (x) 的最小质因子。使用 min_25 筛的思想，设 (f_1(x, k) = sum_limits{i = 1}^x [i in P 或 { m minp}(i) > P_k]mu(i)g(i))，(f_2(x, k) = sum_limits{i = 1}^x [i in P 或 { m minp}(i) > P_k]mu(i))，那么 (w) 即为 (f_1(n, 0))。

我们按 (k) 从大到小求 (f_1) 与 (f_2)。考虑加上最小质因子为 (P_k) 的数的贡献，那么可以得到 (f_2) 的转移如下（注意随着质因子的增加，(mu) 的符号会改变，因此新增的贡献前带有负号）：

[f_2(x, k - 1) = egin{cases} f_2(x, k) + (-f_2(leftlfloorfrac{x}{P_k} ight floor, k)) - (-f_2(P_k, k)), & P_k^2 leq x \ f_2(x, k), & P_k^2 > xend{cases} ]

同理可得 (f_1) 的转移如下：

[f_1(x, k - 1) = egin{cases} f_1(x, k) + (-f_1(leftlfloorfrac{x}{P_k} ight floor, k)) - (-f_1(P_k, k)) + (-f_2(leftlfloorfrac{x}{P_k} ight floor, k)) - (-f_2(P_k, k)), & P_k^2 leq x \ f_1(x, k), & P_k^2 > xend{cases} ]

总时间复杂度即为 min_25 筛的时间复杂度，为 (O(frac{n^{frac{3}{4}}}{log n}))。

代码

#include<bits/stdc++.h>

using namespace std;

const int N = 1e6 + 10;
const long long mod = 1e12 + 39;

long long n, powk[N];
int tt, sq;
bool is_prime[N];
vector<int> primes;
vector<long long> values;

template<typename T>
struct my_array {
  T a[N];

  T& operator [] (long long x) {
    return x <= sq ? a[x] : a[n / x + sq];
  }
};

my_array<long long> g0, g1, f1, f2;

long long add(long long x, long long y, long long md) {
  long long t = x + y;
  if (t >= md) {
    t -= md;
  }
  return t;
}

long long sub(long long x, long long y, long long md) {
  long long t = x - y;
  if (t < 0) {
    t += md;
  }
  return t;
}

long long mul(long long x, long long y, long long md) {
  long double t = (long double) x * y;
  return (x * y - (long long) (t / md) * md) % md;
}

long long qpow(long long v, long long p) {
  long long result = 1;
  for (; p; p >>= 1, v = mul(v, v, mod)) {
    if (p & 1) {
      result = mul(result, v, mod);
    }
  }
  return result;
}

void sieve(int n) {
  memset(is_prime, true, sizeof is_prime);
  for (int i = 2; i <= n; ++i) {
    if (is_prime[i]) {
      primes.push_back(i);
    }
    for (auto v : primes) {
      long long d = (long long) v * i;
      if (d > n) {
        break;
      }
      is_prime[d] = false;
      if (i % v == 0) {
        break;
      }
    }
  }
}

long long same_diff(long long x) {
  long long y = x + 1;
  if (x & 1) {
    y >>= 1;
  } else {
    x >>= 1;
  }
  return mul(x, y, mod - 1);
}

void min_25_sieve() {
  for (long long i = 1; i <= n; i = n / (n / i) + 1) {
    values.push_back(n / i);
  }
  for (auto x : values) {
    g0[x] = x - 1;
    g1[x] = sub(same_diff(x), 1, mod - 1);
  }
  for (auto p : primes) {
    for (auto x : values) {
      if (x < (long long) p * p) {
        break;
      }
      long long y = x / p;
      long long g0_t = g0[y] - g0[p - 1];
      long long g1_t = sub(g1[y], g1[p - 1], mod - 1);
      g0[x] -= g0_t;
      g1[x] = sub(g1[x], mul(p, g1_t, mod - 1), mod - 1);
    }
  }
  reverse(primes.begin(), primes.end());
  for (auto x : values) {
    f1[x] = f2[x] = sub(0, g0[x], mod - 1);
  }
  for (auto p : primes) {
    for (auto x : values) {
      if (x < (long long) p * p) {
        break;
      }
      long long y = x / p;
      long long f2_t = sub(f2[p], f2[y], mod - 1);
      long long f1_t = sub(f1[p], f1[y], mod - 1);
      f2[x] = add(f2[x], f2_t, mod - 1);
      f1[x] = add(f1[x], add(f1_t, f2_t, mod - 1), mod - 1);
    }
  }
}

int main() {
  scanf("%d", &tt);
  while (tt--) {
    scanf("%lld", &n);
    sq = sqrt(n);
    primes.clear();
    values.clear();
    memset(powk, 0, sizeof powk);
    sieve(sq);
    min_25_sieve();
    long long answer = 1;
    int maxk = 1;
    for (auto p : primes) {
      long long t = p;
      for (int j = 1; t <= n; ++j, t *= p) {
        long long num1 = mul(t, same_diff(n / t), mod - 1);
        long long num2 = mul(t * p, same_diff(n / t / p), mod - 1);
        maxk = max(maxk, j + 1);
        powk[j + 1] = add(powk[j + 1], sub(num1, num2, mod - 1), mod - 1);
      }
    }
    for (int i = 2; i <= maxk; ++i) {
      answer = mul(answer, qpow(i, powk[i]), mod);
    }
    long long powm = 0;
    for (long long i = sq + 1; i <= n; i = n / (n / i) + 1) {
      long long j = n / (n / i);
      powm = add(powm, mul(sub(g1[j], g1[i - 1], mod - 1), same_diff(n / i), mod - 1), mod - 1);
    }
    answer = mul(answer, qpow(2, powm), mod);
    answer = mul(answer, qpow(2, f1[n]), mod);
    printf("%lld
", answer);
  }
  return 0;
}