【bzoj 3512】DZY loves math IV

求$\sum_{i=1}^n \sum_{j=1}^m \phi(ij)$，$n$不大。

$n$不大摆明了就是暴力对吧。
考虑$S(n,m)=\sum_{i=1}^m\phi(in)$，我们先来看一个搞笑做法。。。
首先$n$内幂指数大于1的质因子可以拎出来到只剩一次，于是我们只考虑$\mu(n)\neq 0$的$n$。接着考虑取$n$的任意一个质因子$p$，对于$p\mid i$的指标$i$，我们规定此时$p$对$\phi(in)$的贡献为$p$，$p-1$的贡献由$i$提供，于是也可以提出来。而对于$p\nmid i$的指标$i$，$p$的贡献即为$p-1$。综合两个考虑，如果我们直接把前一类的$i$当做后一类来处理了，那么我们就要补回$1$的贡献。于是我们有

$$S(n,m)=(p-1)\sum_{i=1}^m\phi(\frac{n}{p}i)+\sum_{i=1}^m[i\%p==0] \phi(i\frac{m}{p})\\=(p-1)S(n/p,m)+S(n,\lfloor m/p \rfloor)$$
递归终止条件是 $n==1$，然后上杜教筛。看起来很靠谱。。？
然而仔细想想这就是个。。。爆搜。。。对吧= =但是又因为 $n$很小所以不同的质因子数也很少。。最多只有6。。。。所以总共可能的状态大概只有2kw个左右。。
于是我就写了这个玩意。。。然后华丽的狗带了。。。冷静了一下发现。。。空间限制只有128MB？那我2kw个数挂链表都存不下啊。。？一番剪枝剪到了1.3kw。。再也卡不动空间了。。当然如果有谁能卡过去那他很niubia。。。于是重新推了推式子看看能不能搞搞别的。。。

下面是方法二。
考虑直接展开$\phi(in)$会是个什么东西。
取 $d=\gcd(i,n)$，因为 $\mu(n)\neq 0$所以也有 $\mu(d)\neq0$。。我们单独考虑一个 $d$的质因子$p$，而前面那个条件就允许我们直接令 $phi$中 $i$的$p$所贡献的乘积中包含了 $p-1$。。因此我们希望 $n$中的$p$贡献要恰好为 $p$。。。而实际上这也很好处理。。。当$phi(x)$中包含 $p$时，$p$的贡献为 $p-1$，否则为 $1$，而又由于他是个积性的玩意所以每个质因子的贡献独立。。。所以如果我们算$\sum_{k\mid d}\phi(n/k)$的话 $p$在这里面的贡献就刚好为$p$了。因此我们得到 $\phi(in)= \phi(i)\sum_{k\mid \gcd(i,n)}\phi(n/k)$。
然后一波代入式子化简就可以得到 $S(n,m)=\sum_{k\mid n} \phi(n/k)S(k,m/k)$，这样复杂度就是对的了。
你问我复杂度？大概是 $O(n^{3/2}+m^{2/3})$吧。。

代码：

arr mn;
adj vis , phi , pr;
int ptt;

int n , m , N;

void input() {
    n = rd() , m = rd();
    N = (int) (pow(m , 0.7) + 0.5);
    upmax(N , n);
    rep (i , 2 , N) {
        if (!vis[i]) {
            pr[++ ptt] = i;
            phi[i] = i - 1;
            if (i <= n)
                mn[i] = i;
        }
        for (int j = 1; j <= ptt && i * pr[j] <= N; j ++) {
            int t = i * pr[j];
            vis[t] = 1;
            if (t <= n) mn[t] = min(mn[i] , pr[j]);
            if (i % pr[j] == 0) {
                phi[t] = pr[j] * phi[i];
                break;
            }
            phi[t] = (pr[j] - 1) * phi[i];
        }
    }
    phi[1] = 1;
    rep (i , 2 , N) phi[i] = add(phi[i] , phi[i - 1]);
}

inline int get(int x) {
    int t = 1;
    while (x > 1) {
        int y = mn[x];
        x /= y;
        t *= y;
        while (x % y == 0) x /= y;
    }
    return t;
}

__gnu_pbds::gp_hash_table<int , int> ms;

int dls(int n) {
    if (n <= N)
        return phi[n];
    if (ms.find(n) != ms.end())
        return ms[n];
    int ret = 1ll * n * (n + 1) / 2 % mod;
    for (int l = 2 , r = 0; l <= n; l = r + 1) {
        r = n / (n / l);
        ret = dec(ret , mul(r - l + 1 , dls(n / l)));
    }
    ms[n] = ret;
    return ret;
}

__gnu_pbds::gp_hash_table<int , int> exi[maxn];

int gao(int n , int m) {
    if (n == 1)
        return dls(m);

    if (!m) return 0;
    if (m == 1) return phi[n] - phi[n - 1];
    if (exi[n].find(m)!=exi[n].end()) return exi[n][m];

    int t = 0;
    for (int i = 1; i * i <= n; i ++) if (n % i == 0) {
        int j = n / i;
        t = add(t , mul(phi[j] - phi[j - 1] , gao(i , m / i)));
        if (i != j) {
            t = add(t , mul(phi[i] - phi[i - 1] , gao(j , m / j)));
        }
    }

    exi[n][m]=t;

    return t;
}

void solve() {
    int ans = 0;
    rep (i , 1 , n) {
        int x = get(i);
        int tmp = gao(x , m);
        ans = add(ans , mul(i / x , tmp));
    }
    printf("%d\n" , ans);
}