SPOJ VLATTICE 莫比乌斯反演

题目链接:https://www.spoj.com/problems/VLATTICE/en/

VLATTICE - Visible Lattice Points

Description

Consider a NNN lattice. One corner is at (0,0,0) and the opposite one is at (N,N,N). How many lattice points are visible from corner at (0,0,0) ? A point X is visible from point Y iff no other lattice point lies on the segment joining X and Y.

Input

The first line contains the number of test cases T. The next T lines contain an interger N

Output

Output T lines, one corresponding to each test case.

Sample Input

3
1
2
5

Sample Output

7
19
175

Constraints

T <= 50
1 <= N <= 1000000

Solution

一直想着学一学莫比乌斯反演，一直都没看明白，今天终于学会(套公式)了
之前在洛谷做过一道类似的题(不过那个题是在一个二维平面里),看到这个题就想到了求(gcd(i,j,k)=1)的个数，即$$sum_{i=0}^nsum_{j=0}nsum_{k=0}^n(i,j,k)=1$$
看到gcd(i,j,k)=1,想到迪利克雷卷积的单位元函数好像可以解决,(epsilonig(gcd(i,j,k)ig)),((epsilon(n))为迪利克雷卷积单位元,(epsilon(n)=[n=1]))，接下来我们再用单位元函数和莫比乌斯函数的关系(epsilon(n)=sum_{d|n}mu(d)),就可以推出如下计算公式

[sum_{i=0}^nsum_{j=0}^nsum_{k=0}^nsum_{d|(i,j,k)}mu(d)$$, 然后我们调整求和顺序，将莫比乌斯函数放到前面，那么可以得到 $$sum_{d=1}^nmu(d)sum_{i=0}^n dmid i sum_{j=0}^n dmid j sum_{k=0}^n dmid k,(mbox{注意d是从1开始的，因为除数不能为0})]

接着看后面三个和式，都是要求n里有多少个数是d的倍数，那么由整数分块儿可知n里是d的倍数的数由n/d个，那么我们可以得到如下式子

[ans=sum_{d=1}^nmu(d)(n/d)(n/d)(n/d) ]

不过显然是没有包含i,j,k=0的情况，那么这个0应该怎么处理呢，
我们先考虑ijk只有一个为0，i=0时是指在(yoz)平面，也就是说i=0时对应着(yoz)平面中的点，而上式中的倍数关系只考虑(ygt 0),(zgt 0)的点，那么我们在上式添上一个((n/d)*(n/d))，就得到了yoz平面上不在坐标轴上的答案
同理可以得知j=0，k=0的情况一样，加上((n/d)*(n/d))即可
剩下i,j,k两个为0的情况就是坐标轴上，坐标轴显然只有(0,0,1),(0,1,0),(1,0,0)三个点贡献答案，最后加上即可
那么我们就得到了最后的式子

[ans=3+sum_{d=1}^nmu(d)(n/d)^3+3(n/d)^2 ]

O(n)扫一遍即可

#include <algorithm>
#include <cctype>
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <iostream>
#include <map>
#include <queue>
#include <set>
#include <stack>
#include <unordered_map>
#include <unordered_set>
#include <vector>
#define lson rt << 1, l, mid
#define rson rt << 1 | 1, mid + 1, r
#define LONG_LONG_MAX 9223372036854775807LL
#define ll LL
using namespace std;
typedef long long ll;
typedef long double ld;
typedef unsigned long long ull;
typedef pair<int, int> P;
int n, m, k;
const int maxn = 1e6 + 10;
template <class T>
inline T read()
{
    int f = 1;
    T ret = 0;
    char ch = getchar();
    while (!isdigit(ch))
    {
        if (ch == '-')
            f = -1;
        ch = getchar();
    }
    while (isdigit(ch))
    {
        ret = (ret << 1) + (ret << 3) + ch - '0';
        ch = getchar();
    }
    ret *= f;
    return ret;
}
template <class T>
inline void write(T n)
{
    if (n < 0)
    {
        putchar('-');
        n = -n;
    }
    if (n >= 10)
    {
        write(n / 10);
    }
    putchar(n % 10 + '0');
}
int mu[maxn], prime[maxn], tot, vis[maxn];
void init()
{
    mu[1] = 1;
    for (int i = 2; i < maxn; ++i)
    {
        if (!vis[i])
            prime[++tot] = i, mu[i] = -1;
        for (int j = 1; j <= tot && i * prime[j] < maxn; ++j)
        {
            vis[i * prime[j]] = 1;
            if (i % prime[j] == 0)
            {
                mu[i * prime[j]] = 0;
                break;
            }
            mu[i * prime[j]] = -mu[i];
        }
    }
}
ll solve()
{
    ll res = 0;
    for (int i = 1; i <= n; i++)
        res += (ll)mu[i] * (n / i) * (n / i) * (n / i + 3);
    return res + 3;
}
int main(int argc, char const *argv[])
{
#ifndef ONLINE_JUDGE
    freopen("in.txt", "r", stdin);
    freopen("out.txt", "w", stdout);
#endif
    init();
    int t = read<int>();
    while (t--)
    {
        n = read<int>();
        write(solve());
        puts("");
    }
    return 0;
}