huntian oy HDU

huntian oy (HDU - 6706)

题意:

(n,a,b),求(f(n,a,b)=sum_{i=1}^nsum_{j=1}^igcd(i^a-j^a,i^b-j^b)[gcd(i,j)=1]\%(10^9+7))。保证(gcd(a,b)=1)(n<=1e9)

题解:

  • (gcd(i^a-j^a,i^b-j^b) = i^{gcd(a,b)}-j^{gcd(a,b)})
  • (sum_{i=1}^ni[gcd(i,n)=1]=frac{nphi(n)+[n=1]}{2})

根据上面两个公式。

[f(n,a,b)=sum_{i=1}^nsum_{j=1}^igcd(i^a-j^a,i^b-j^b)[gcd(i,j)=1] ]

[=sum_{i=1}^nsum_{j=1}^i(i-j)[gcd(i,j)=1] ]

[=sum_{i=1}^n(i*phi(i)-sum_{j=1}^ij[gcd(i,j)=1]) ]

[=sum_{i=1}^n(i*phi(i)-frac{i*phi(i)}{2})-frac{1}{2}(这里要减去frac{1}{2},是特殊考虑了i=1的情况) ]

[=frac{1}{2}(sum_{i=1}^{n}(i*phi(i))-1) ]

对于(sum_{i=1}^{n}i*phi(i)),可以用杜教筛求:(杜教筛学习推荐这篇博客)

对于(f(n)=i*phi(i)),构造(g(n)=n)。根据迪利克雷卷积设(h=f*g),则

[h(n)=sum_{d|n}f(d)*g(frac{n}{d}) ]

[=sum_{d|n}(d*phi(d))*frac{n}{d} ]

[=sum_{d|n}n*phi(d) ]

[=nsum_{d|n}phi(d)=n^2 ]

所以根据杜教筛,设(S(n)=sum_{i=1}^{n}i*phi(i)),则

[S(n)=sum_{i=1}^{n}i^2-sum_{d=2}^{n}d*S(n) ]

数论分块+递归计算(S(n))就好。

代码:

#include <bits/stdc++.h>
#define fopi freopen("in.txt", "r", stdin)
#define fopo freopen("out.txt", "w", stdout)
using namespace std;
typedef long long LL;
const int MOD = 1e9+7;
const LL inv2 = 5e8+4;
const LL inv6 = 166666668;
typedef long long LL;

const int maxn = 1e6 + 10;
int check[maxn], phi[maxn], prime[maxn];
LL sum[maxn];
unordered_map<int, LL> d;

void init(int N) {
    memset(check, false, sizeof(check));
    phi[1] = 1;
    int tot = 0;
    for (int i = 2; i <= N; i++) {
        if (!check[i]) {
            prime[tot++] = i;
            phi[i] = i-1;
        }
        for (int j = 0; j < tot; j++) {
            if (i * prime[j] > N) break;
            check[i * prime[j]] = true;
            if (i % prime[j] == 0) {
                phi[i * prime[j]] = phi[i] * prime[j];
                break;
            }
            else {
                phi[i * prime[j]] = phi[i] * (prime[j]-1);
            }
        }
    }
    for (int i = 1; i <= N; i++)
        sum[i] = (sum[i-1] + 1ll*i*phi[i] % MOD) % MOD;
}

LL solve(int x) {
    return 1ll * x * (x+1) % MOD * (2*x+1) % MOD * inv6 % MOD;
}

int dfs(int x) {
    if (x <= 1000000) return sum[x];
    if (d[x]) return d[x];
    LL ans = solve(x);
    for (int l = 2, r; l <= x; l = r+1) {
        r = x/(x/l);
        ans = (ans - (1LL*(r-l+1)*(l+r)/2) % MOD * dfs(x/l) % MOD + MOD) % MOD;
    }
    return d[x] = ans;
}

int T;
int n, a, b;
int main() {
    init(1000000);
    scanf("%d", &T);
    for (int ca = 1; ca <= T; ca++) {
        scanf("%d%d%d", &n, &a, &b);
        printf("%lld
", inv2 * (dfs(n) - 1 + MOD) % MOD);
    }
}
原文地址:https://www.cnblogs.com/ruthank/p/11469173.html