【杜教筛】HUD-6706 huntian oy

传送门

题意:f(n,a,b)=\sum_{i=1}^n \sum_{j=1}^i gcd(i^a-j^a,i^b-j^b)[gcd(i,j)=1]\%(10^9+7)
给定n,a,b求f

打表可知f(n,a,b)=\sum_{i=1}^n \sum_{j=1}^i (i-j)[gcd(i,j)=1]\%(10^9+7)

即:\sum\limits_{i=1}^n\sum\limits_{j=1}^i\ i\ [gcd(i,j)=1]-\sum\limits_{i=1}^n\sum\limits_{j=1}^i\ j\ [gcd(i,j)=1]

在i j互质的时候才有值。前面一部分就是每个数与其互质的个数的和,\sum\limits_{i=1}^{n}i\varphi(i)

后面一部分是每个数,与其互质的数值的和。\sum\limits_{i=1}^{n}\frac{i\varphi(i)}{2}(n>1)

原因是gcd(n,x)=gcd(n,n-x),当出现一个与n互质的数x的时候必然有一个n-x成对出现与其互质。他们的和为\frac{n}{2},但这只是在n>1的时候成立。

由于n==1的时候只有一个数字1,所以式子完整来写应该是\sum\limits_{i=1}^{n}\frac{i\varphi(i)+1}{2}因为在i=1的时候答案是1,如果除2就成了1/2,所以在分子上加个1使其成立

完整式子是\sum\limits_{i=1}^{n}\frac{i\varphi(i)-1}{2}

由于n有1e9,线筛是存不下的。
考虑
f=i\varphi(i)=id\varphi(i)
g=id
f*g(n)=\sum\limits_{i|n}\varphi(i)i\frac{n}{i}=n\sum\limits_{i|n}\varphi(i)

然后套杜教筛的板子。

#include<bits/stdc++.h>
using namespace std;
const int maxn = 5e6 + 233;
typedef long long ll;
const ll mod = 1e9 + 7;
ll inv2, inv6;
bool st[maxn];
int primes[maxn], cnt;
ll phi[maxn];
unordered_map<int, ll> up;
ll ksm(ll n, ll k)
{
    ll res = 1;
    while(k)
    {
        if(k & 1) res = res * n % mod;
        n = n * n % mod;
        k >>= 1;
    }
    return res;
}
void get_primes(int n)
{
    // mu[1] = phi[1] = 1;
    phi[1] = 1;
    for(int i = 2; i <= n; i++) {
        if(!st[i]) {
            primes[cnt++] = i;
            phi[i] = i - 1;
            // mu[i] = -1;
        }
        for(int j = 0; j < cnt && i * primes[j] <= n; j++) {
            st[i * primes[j]] = 1;
            if(i % primes[j] == 0) {
                phi[i * primes[j]] = phi[i] * primes[j];
                // mu[i * primes[j]] = 0;
                break;
            }else {
                phi[i * primes[j]] = phi[i] * (primes[j] - 1);
                // mu[i * primes[j]] = -mu[i];
            }
        }
    }
    for(ll i = 1; i <= n; i++) phi[i] = (phi[i - 1] + phi[i] * i % mod) % mod;
    //phi[1] = 0;
}
ll Get_f_g(ll n)
{
    return 1LL * n * (n + 1) % mod * (n * 2 + 1) % mod * inv6 % mod;
}
ll Get_g(ll n)
{
    return n * (n + 1) % mod * inv2 % mod;
}
ll Getsum(ll n)
{
    if(n < 5000000) return phi[n];
    if(up[n]) return up[n];
    ll ans = Get_f_g(n);
    for(ll l = 2, r; l <= n; l = r + 1)
    {
        r = n / (n / l);
        ans -= (Get_g(r) - Get_g(l - 1)) * Getsum(n / l) % mod;
        ans %= mod;
    }
    return up[n] = ans;
}
int main()
{
    inv2 = ksm(2, mod - 2);
    inv6 = ksm(6, mod - 2);
    int T; cin >> T;
    get_primes(5e6 + 10);
    while(T--)
    {
        ll n, a, b;
        scanf("%lld%lld%lld", &n, &a, &b);
        ll ans = 0;
        ans = (Getsum(n) - 1) * inv2 % mod;
        printf("%lld\n", (ans + mod) % mod);
    }
}

发表评论

邮箱地址不会被公开。 必填项已用*标注