【BZOJ 4176】 Lucas的数论 - 杜教筛

　　题意：求$\sum_{i=1}^n \sum_{j=1}^n \sigma_0(ij)$
　　
　　把约数个数给展开来写，然后注意到$d\mid ij$等价于$\frac{d}{gcd(i,d)}\mid j$，于是代入，然后可以拿个求和指标替换掉$gcd(i,d)$，再然后可以用个$\mu$展开一波，最后会发现有几个相同的求和指标，等价于求一个和然后平方，于是可以得到最后的式子是
　　

$$\sum_{k=1}^n \mu(k) \left( \sum_{p=1}^{\lfloor \frac{n}{k} \rfloor} \lfloor \frac{n}{kp} \rfloor\right)^2$$
　　设 $S_r(n)=\sum_{k=1}^r \lfloor \frac{n}{k} \rfloor$，我们知道有 $\lfloor\frac{\lfloor\frac{n}{k}\rfloor}{p}\rfloor=\lfloor \frac{n}{kp} \rfloor$，那么式子可以变为
　　 $$\sum_{k=1}^n \mu(k) S_{\lfloor \frac{n}{k} \rfloor}(n)^2$$
　　记 $F_n = \sum_{i=1}^n\mu(i)$，则有 $F_n=1-\sum_{k=2}^nF(\lfloor \frac{n}{k} \rfloor)$，那么对于每一段相同的 $n/k$，两个东西都可以直接计算。所以这个的复杂度相当于是一次杜教筛。

懒癌晚期不想写公式结果还是写了那么多(?)

/*
    I will chase the meteor for you, a thousand times over.
    Please wait for me, until I fade forever.
    Just 'coz GEOTCBRL.
*/
#include <bits/stdc++.h>
using namespace std;
#define fore(i,u)  for (int i = head[u] ; i ; i = nxt[i])
#define rep(i,a,b) for (int i = a , _ = b ; i <= _ ; i ++)
#define per(i,a,b) for (int i = a , _ = b ; i >= _ ; i --)
#define For(i,a,b) for (int i = a , _ = b ; i <  _ ; i ++)
#define Dwn(i,a,b) for (int i = ((int) a) - 1 , _ = (b) ; i >= _ ; i --)
#define fors(s0,s) for (int s0 = (s) , _S = s ; s0 ; s0 = (s0 - 1) & _S)
#define foreach(it,s) for (__typeof(s.begin()) it = s.begin(); it != s.end(); it ++)

#define mp make_pair
#define pb push_back
#define pii pair<int , int>
#define fir first
#define sec second
#define MS(x,a) memset(x , (a) , sizeof (x))

#define gprintf(...) fprintf(stderr , __VA_ARGS__)
#define gout(x) std::cerr << #x << "=" << x << std::endl
#define gout1(a,i) std::cerr << #a << '[' << i << "]=" << a[i] << std::endl
#define gout2(a,i,j) std::cerr << #a << '[' << i << "][" << j << "]=" << a[i][j] << std::endl
#define garr(a,l,r,tp) rep (__it , l , r) gprintf(tp"%c" , a[__it] , " \n"[__it == _])

template <class T> inline void upmax(T&a , T b) { if (a < b) a = b ; }
template <class T> inline void upmin(T&a , T b) { if (a > b) a = b ; }

typedef long long ll;

const int maxn = 8000007;
const int maxm = 200007;
const int mod = 1000000007;
const int inf = 0x7fffffff;
const double eps = 1e-7;

typedef int arr[maxn];
typedef int adj[maxm];

inline int fcmp(double a , double b) {
    if (fabs(a - b) <= eps) return 0;
    if (a < b - eps) return -1;
    return 1;
}

inline int add(int a , int b) { return ((ll) a + b) % mod ; }
inline int mul(int a , int b) { return ((ll) a * b) % mod ; }
inline int dec(int a , int b) { return add(a , mod - b % mod) ; }
inline int Pow(int a , int b) {
    int t = 1;
    while (b) {
        if (b & 1) t = mul(t , a);
        a = mul(a , a) , b >>= 1;
    }
    return t;
}

#define gc getchar
#define idg isdigit
#define rd RD<int>
#define rdll RD<long long>
template <typename Type>
inline Type RD() {
    char c = getchar(); Type x = 0 , flag = 1;
    while (!idg(c) && c != '-') c = getchar();
    if (c == '-') flag = -1; else x = c - '0';
    while (idg(c = gc()))x = x * 10 + c - '0';
    return x * flag;
}
inline char rdch() {
    char c = gc();
    while (!isalpha(c)) c = gc();
    return c;
}
#undef idg
#undef gc

// beginning

arr mu , vis , pr;
int tot , N , n;

inline void sieve(int n) {
    mu[1] = 1;
    rep (i , 2 , n) {
        if (!vis[i])
            pr[++ tot] = i , mu[i] = -1;
        rep (j , 1 , tot) if (i * pr[j] > n) break; else {
            vis[i * pr[j]] = 1;
            if (i % pr[j] == 0) break;
            mu[i * pr[j]] = -mu[i];
        }
    }
    rep (i , 1 , n) mu[i] = (mu[i - 1] + mu[i]) % mod;
}

void input() {
    n = rd() , N = int(pow(n , 0.75) + 0.5);
    sieve(N);
}

const int cut = 10007;
map<int , int> ms[cut];

int calc(int n) {
    if (n <= N)
        return mu[n];
    if (ms[n % cut].count(n)) return ms[n % cut][n];
    int ret = 1;
    for (int l = 2 , r = 0; l <= n; l = r + 1) {
        r = n / (n / l);
        int tmp = calc(n / l);
        ret = dec(ret , mul(r - l + 1 , tmp));
    }
    ms[n % cut][n] = ret;
    return ret;
}

inline int get_sum(int n) {
    int ret = 0;
    for (int l = 1 , r = 0; l <= n; l = r + 1) {
        r = n / (n / l);
        ret = add(ret , mul(n / l , r - l + 1));
    }
    return mul(ret , ret);
}

void solve() {
    int ans = 0;
    for (int l = 1 , r = 0; l <= n; l = r + 1) {
        r = n / (n / l);
        int tmp = get_sum(n / l);
        int res = dec(calc(r) , calc(l - 1));
//      gprintf("[%d , %d]\n" , l , r);
//      gout(tmp) , gout(res);
        ans = add(ans , mul(tmp , res));
    }
    if (ans < 0) ans += mod;
    printf("%d\n" , ans);
}

int main() {
    #ifndef ONLINE_JUDGE
        freopen("data.txt" , "r" , stdin);
    #endif
    input();
    solve();
    return 0;
}