這個式子可以遞迴計算，邊界為：

m=0:0

m=1:u(n);

n=1:segma{u(i)|1<=i<=m}

其中segma{u(i)|1<=i<=m}可以用杜教篩在O(m^(2/3))內求出

求解過程中要求的u(i)都可以O(1)求出，做法是先將n唯一分解，然後列舉子集，判一下元素個數，儲存下來以後就可以O(1)求了

注意判一下n是否有平方因子，有的話直接輸出0

#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
#include <map>
using namespace std;
typedef long long ll;
const int maxn = 5e6 + 10;
const int N = 5e6;
bool vis[maxn];
int prime_cnt, prime[maxn], mu[maxn];
map<ll, ll> dp;
void init() {
    mu[1] = 1;
    for (int i = 2; i <= N; i++) {
        if (!vis[i]) prime[++prime_cnt] = i, mu[i] = -1;
        for (int j = 1; j <= prime_cnt && i * prime[j] <= N; j++) {
            vis[i*prime[j]] = 1;
            if (i % prime[j] == 0) break;
            mu[i*prime[j]] = -mu[i];
        }
    }
    for (int i = 2; i <= N; i++) mu[i] += mu[i-1];
}
ll cal(ll x) {
    if (x <= N) return mu[x];
    if (dp[x]) return dp[x];
    ll ans = 1;
    for (ll i = 2, t; i <= x; i = t + 1) {
        t = x/(x/i);
        ans -= cal(x/i)*(t-i+1);
    }
    return dp[x] = ans;
}
int p_cnt;
ll p[20];
int cnt[5000];
ll val[5000];
ll u(ll x) { return ((cnt[x]&1) ? -1 : 1); }
ll solve(ll m, ll n) {
    if (cnt[n] == 0) return cal(m);
    if (m == 0) return 0;
    if (m == 1) u(n);
    ll ans = 0, i = n;
    for (ll i = n; ; i = (i-1)&n) {
        ans += u(i) * solve(m/val[i], i);
        if (!i) break;
    }
    return ans * u(n);
}
int main() {
    init();
    ll n, m;
    scanf("%lld%lld", &m, &n);
    ll x = n;
    p_cnt = 0;
    bool flag = true;
    for (ll i = 2; i * i <= x; i++) {
        if (x % i == 0) {
            if ((x/i) % i == 0) { flag = 0; break; }
            p[p_cnt++] = i;
            while (x % i == 0) x /= i;
        }
    }
    if (!flag) {
        printf("0\n"); return 0;
    }
    if (x > 1) p[p_cnt++] = x;
    for (int i = 0; i < (1<<p_cnt); i++) {
        cnt[i] = 0, val[i] = 1;
        for (int j = 0; j < p_cnt; j++) {
            if (i & (1<<j)) cnt[i]++, val[i] *= p[j];
        }
    }
    printf("%lld\n", solve(m, (1<<p_cnt)-1));
    return 0;
}