HDU 6706 杜教筛

传送门： http://acm.hdu.edu.cn/showproblem.php?pid=6706

题意

$$求解\sum _{i=1}^n\sum _{j=1}^{i}gcd(i^a-j^a,i^b-j^b)[gcd(i,j)=1]$$

思路

$$gcd(i^a-j^a,i^b-j^b)=i^{gcd(a,b)}-j^{gcd(a,b)}$$

$题目说gcd(a,b)=1，所以那个复杂的式子直接变成i-j.$

$$\sum _{i=1}^n\sum _{j=1}^i(i-j)[gcd(i,j)=1]$$

$$\sum _{i=1}^n\sum _{j=1}^{i}i[gcd(i,j)=1]-\sum _{i=1}^n\sum _{j=1}^{i}j[gcd(i,j)=1]$$

$左边为与i互质的个数，右边为比i小并且互质的数字和，即为：$

$$\sum _{i=1}^{n}i\phi (i)-\frac{{\sum }_{i=1}^{n}i\phi (i)+1}{2}$$

$$\frac{{\sum }_{i=1}^{n}i\phi (i)-1}{2}$$

$因为n有1e9，所以需要杜教筛优化。$

$$f\ast g(n)=\sum _{d|n}f(d)g(\frac{n}{d})=\sum _{d|n}id(d)\ast \phi (d)\ast g(\frac{n}{d})$$

$很显然，设g=id，则：$

$$f\ast g(n)=\sum _{d|n}\phi (d)id(n)=n\sum _{d|n}\phi (d)=id(n)(\phi \ast I)=i{d}^2(n)=n^2$$

$$所以得：S(n)=\sum _{i=1}^n(f\ast g)i-\sum _{i=2}^{n}g(i)S(\frac{n}{i})=\frac{n(n+1)(2n+1)}{6}-\sum _{i=2}^{n}iS(\frac{n}{i})$$

Code

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

typedef long long ll;
typedef pair<int, int> pii;
typedef pair<double, double> pdd;
typedef pair<ll, ll> pll;

#define endl "\n"
#define eb emplace_back
#define mem(a, b) memset(a , b , sizeof(a))

const ll INF = 0x3f3f3f3f;
// const ll mod = 998244353;
const ll mod = 1e9 + 7;
const double eps = 1e-6;
const double PI = acos(-1);
const double R = 0.57721566490153286060651209;

template<typename T>
inline void read(T &a) {char c = getchar();T x = 0, f = 1;
    while (!isdigit(c)) {if (c == '-')f = -1;c = getchar();}
    while (isdigit(c)) {x = (x << 1) + (x << 3) + c - '0';c = getchar();}a = f * x;
}

const int N = 1e6 + 10;

bool is_prime[N];
int prime[N], cnt;
int phi[N];
ll sum_i_phi[N];

void init() {
    phi[1] = 1;
    for(int i = 2;i < N; i++) {
        if(!is_prime[i]) prime[++cnt] = i, phi[i] = i - 1;
        for(int j = 1;j <= cnt && i * prime[j] < N; j++) {
            is_prime[i * prime[j]] = 1;
            if(i % prime[j] == 0) {
                phi[i * prime[j]] = phi[i] * prime[j];
                break;
            }
            else phi[i * prime[j]] = phi[i] * phi[prime[j]];
        }
    }
    for(int i = 1;i < N; i++) {
        sum_i_phi[i] = (sum_i_phi[i - 1] + 1ll * i * phi[i] % mod) % mod;
    }
}

ll quick_pow(ll a, ll b) {
    ll ans = 1;
    while(b) {
        if(b & 1) ans = ans * a % mod;
         a = a * a % mod;
         b >>= 1;
    }
    return ans % mod;
}

ll inv2;
ll inv6;

map<ll, ll> mp;

ll S(ll x) {
    if(x < N) return sum_i_phi[x];
    else if(mp[x]) return mp[x];
    ll ans = x * (x + 1) % mod * (2 * x + 1) % mod * inv6 % mod;
    for(int l = 2, r;l <= x; l = r + 1) {
        r = min(x, x / (x / l));
        ans = (ans - 1ll * (r + l) % mod * (r - l + 1) % mod * inv2 % mod * S(x / l)) % mod;
    }
    return mp[x] = ans;
}

void solve() {
    init();
    inv2 = quick_pow(2, mod - 2);
    inv6 = quick_pow(6, mod - 2);
    int _; scanf("%d",&_);
    while(_--) {
        ll n, a, b; scanf("%lld%lld%lld",&n,&a,&b);
        ll ans = (S(n) - 1) % mod;
        ans = ans * inv2 % mod;
        printf("%lld\n",(ans % mod + mod) % mod);
    }
}

signed main() {
    ios_base::sync_with_stdio(false);
    // cin.tie(nullptr);
    // cout.tie(nullptr);
#ifdef FZT_ACM_LOCAL
    freopen("in.txt", "r", stdin);
    freopen("out.txt", "w", stdout);
    signed test_index_for_debug = 1;
    char acm_local_for_debug = 0;
    do {
        if (acm_local_for_debug == '$') exit(0);
        if (test_index_for_debug > 20)
            throw runtime_error("Check the stdin!!!");
        auto start_clock_for_debug = clock();
        solve();
        auto end_clock_for_debug = clock();
        cout << "Test " << test_index_for_debug << " successful" << endl;
        cerr << "Test " << test_index_for_debug++ << " Run Time: "
             << double(end_clock_for_debug - start_clock_for_debug) / CLOCKS_PER_SEC << "s" << endl;
        cout << "--------------------------------------------------" << endl;
    } while (cin >> acm_local_for_debug && cin.putback(acm_local_for_debug));
#else
    solve();
#endif
    return 0;
}