1. 程式人生 > >51nod1238 最小公倍數之和 V3 莫比烏斯函式 杜教篩

51nod1238 最小公倍數之和 V3 莫比烏斯函式 杜教篩

題意:求\(\sum_{i = 1}^{n}\sum_{j = 1}^{n}lcm(i, j)\).
  題解:因為是用的莫比烏斯函式求的,所以推導比大部分題解多。。。而且我寫式子一般都比較詳細,所以可能看上去很多式子,實際上是因為每一步都寫了,幾乎沒有跳過的。所以應該都可以看懂的。
  末尾的\(e\)函式是指的\(e[1] = 1\),\(e[x] = 0(x != 1)\)這樣一個函式
  \[\sum_{i = 1}^{n}\sum_{j = 1}^{n}lcm(i, j)\]
  \[\sum_{i = 1}^{n} \sum_{i = 1}^{n} \frac{ij}{gcd(i, j)}\]


  列舉\(gcd\)
  \[\sum_{d = 1}^{n} \sum_{i = 1}^{\lfloor {\frac{n}{d}} \rfloor} \sum_{j = 1}^{\lfloor {\frac{n}{d}} \rfloor}[gcd(i, j) == d] \frac{ij}{d}\]
  因為\((\frac{ijd^2}{d} = ijd)\),所以:
  \[\sum_{d = 1}^{n} \sum_{i = 1}^{\lfloor {\frac{n}{d}} \rfloor} \sum_{j = 1}^{\lfloor {\frac{n}{d}} \rfloor}[gcd(i, j) == d] ijd\]

  \[\sum_{d = 1}^{n}d\sum_{i = 1}^{\lfloor {\frac{n}{d}} \rfloor} \sum_{j = 1}^{\lfloor {\frac{n}{d}} \rfloor}ij[gcd(i, j) == 1]\]
  \[\sum_{d = 1}^{n} d \sum_{i = 1}^{\lfloor {\frac{n}{d}} \rfloor} \sum_{j = 1}^{\lfloor {\frac{n}{d}} \rfloor} ij \sum_{k | gcd(i, j)} \mu(k)\]
  列舉k,再列舉k的倍數。
  \[\sum_{d = 1}^{n}d\sum_{k = 1}^{\lfloor {\frac{n}{d}} \rfloor}\mu(k) \sum_{i = 1}^{\lfloor {\frac{n}{dk}} \rfloor}ik \sum_{j = 1}^{\lfloor {\frac{n}{dk}} \rfloor}jk\]

  設\(S(n) = \sum_{i = 1}^{n}i\)
  \[\sum_{d = 1}^{n}d \sum_{k = 1}^{\lfloor {\frac{n}{d}} \rfloor} \mu(k) k ^ 2 S(\frac{n}{dk})\]
  列舉\(T = dk\)
  \[\sum_{T = 1}^{n} S(\frac{n}{T})^ 2 \sum_{k | T} \mu(k) k ^ 2 \frac{T}{k}\]
  \[\sum_{T = 1}^{n} S(\frac{n}{T})^ 2 \sum_{k | T} \mu(k) kT\]
  \[\sum_{T = 1}^{n} S(\frac{n}{k})^ 2 \cdot T \sum_{k | T} \mu(k)k\]
  設\(f(T) = T\sum_{k | T} \mu(k) k\),捲上\(id^2\),因為\(S(\frac{n}{k})\)可以數論分塊,所以我們只需要快速求出區間\([l, r]\)內的\(f\)之和即可,顯然求出\(f\)的字首和即可解決問題
  \[(f * id^2)(n) = \sum_{i |n}f(i) \frac{n^2}{i^2}=\sum_{i | n}i \sum_{k | i} \mu(k) k \frac{n ^ 2}{i ^ 2}\]
  \[\sum_{i | n}\sum_{k | i} \mu(k)k\frac{n ^ 2}{i} = n \sum_{i | n}\sum_{k | i} \mu(k) k \frac{n}{i}\]
  設\[h(i) = \sum_{k | i} \mu(k)k\],則原式:
  \[n \sum_{i | n} h(i) \frac{n}{i} = n (h * id)(n)\]
  \[(f * id ^ 2)(n) = n (h * id)(n)\]
  \[h(n) = \sum_{k | n}\mu(k)k = (\mu \cdot id) * 1\]
  \[f * id ^ 2 = n [(\mu \cdot id) * 1 * id] = n[(\mu \cdot id) * id * 1]\]
  其中\[(\mu \cdot id) * id = \sum_{i | n} \mu(i) i \frac{n}{i} = n \sum_{i | n}\mu(i) = e\]
  所以
  \[n[(\mu \cdot id) * id * 1] = n[e * 1] = n\]
  帶入杜教篩的式子:
  \[g(1)S(n) = \sum_{i = 1}^{n} (f * g)(i) - \sum_{i = 2}^{n}g(i)S(\frac{n}{i})\]
  \[= \sum_{i = 1}^{n}i - \sum_{i = 2}^{n}i ^ 2 S(\frac{n}{i})\]
  然後直接上杜教就可以了.
  其實還有一個問題。。。一開始預處理的字首和怎麼求?
  要知道字首和,首先要求出\(f\).
  因為\(f(T) = T\sum_{k | T} \mu(k) k\),所以如果我們可以快速求出\(\sum_{k | T}\mu(k)k\),然後就只需要再\(O(n)\)的乘上\(T\)就可以了.
  我們先預處理出\(\mu(k)\),然後對於每一個\(k\),列舉它的倍數,統計貢獻。那麼複雜度為 \(\frac{n}{1} + \frac{n}{2} + ... + \frac{n}{n} = nlogn\)(此處的\(n\)為原題面的\(\frac{2}{3}\)次方,即要預處理的\(f\)個數)

#include<bits/stdc++.h>
using namespace std;
#define R register int
#define LL long long
#define RL register LL
#define AC 3000
#define ac 5000000
#define p 1000000007LL
//#define h(x) ((x <= block) ? sum[x] : S[n / x])

LL n, ans, block; 
LL mu[ac], S[AC], sum[ac], inv[AC];
int pri[ac], tot;
bool z[ac], vis[AC];

inline LL h(LL x)
{
    return ((x <= block) ? sum[x] : S[n / x]);
}

inline void up(LL & a, LL b)
{
    a += b;
    if(a >= p) a -= p;
    if(a <= -p) a += p;
}

LL count(LL l, LL r){
    return (r - l + 1) % p * ((r + l) % p) % p * inv[2] % p;
}

void pre()
{
    scanf("%lld", &n), block = pow(n, 0.66666);
    mu[1] = 1;
    for(R i = 2; i <= block; i ++)
    {
        if(!z[i]) pri[++ tot] = i, mu[i] = -1;
        for(R j = 1; j <= tot; j ++)
        {
            int now = pri[j];
            if(i * now > block) break;
            z[i * now] = true;
            if(!(i % now)) break;
            mu[i * now] = - mu[i];
        }
    }
    inv[1] = 1;
    for(R i = 2; i <= 10; i ++) inv[i] = (p - p / i) * inv[p % i] % p;
    for(R i = 1; i <= block; i ++)//列舉mu(i)
        for(R j = 1; j; j ++)//列舉i的倍數
        {
            if(j * i > block) break;
            up(sum[i * j], mu[i] * i % p);
        }
    for(R i = 1; i <= block; i ++) sum[i] = sum[i] * i % p; 
    for(R i = 1; i <= block; i ++) up(sum[i], sum[i - 1]);//算出f陣列後還要統計字首和
}

LL get(LL x)
{
    x %= p;
    return x * (x + 1) % p * (2 * x + 1) % p * inv[6] % p;
}

void cal(LL x)
{
    if(x <= block || vis[n / x]) return ;
    LL rnt = count(1, x);
    for(RL i = 2, lim, now; i <= x; i = lim + 1)
    {
        lim = x / (x / i), now = x / i, cal(now);
        up(rnt, - ((get(lim) - get(i - 1)) % p * h(now) % p));
    }
    S[n / x] = rnt, vis[n / x] = true;
}

void work()
{
    for(RL i = 1, lim, now, x; i <= n; i = lim + 1)
    {
        lim = n / (n / i), now = n / i, x = count(1, now);
        up(ans, (x % p * x % p * ((h(lim) - h(i - 1)) % p) % p));
    }
    printf("%lld\n", (ans + p) % p);
}

int main()
{
    //freopen("in.in", "r", stdin);
    pre();
    cal(n);
    work();
//  fclose(stdin);
    return 0;
}