BZOJ 3512: DZY Loves Math IV [杜教篩]

3512: DZY Loves Math IV

題意：求\(\sum_{i=1}^n \sum_{j=1}^m \varphi(ij)\)，\(n \le 10^5, m \le 10^9\)

---

n較小，考慮寫成前綴和的形式，計算\(S(n,m)=\sum_{i=1}^m \varphi(in)\)

一開始想出
\[ n= \prod_i p_i,\ \varphi(in) = \varphi(i) \cdot \varphi(\frac{n}{d})\cdot d,\ d=(n,i) \]
比較好想，共有的質因子應該乘\(p\)而不是\(p-1\)

而後帶進去枚舉gcd用莫比烏斯反演的套路，中間的函數很奇怪很差算前綴和...

orz了題解，發現題解使用\(\varphi * 1 =id\)來替換
\[ \varphi(in) = \varphi(i) \cdot \varphi(\frac{n}{d})\cdot \sum_{e\mid d} \varphi(e) = \varphi(i) \cdot \sum_{e\mid d}\varphi(\frac{n}{e}) \]
由於n是不一樣質因子的乘積，因此能夠把兩個\(\varphi\)乘起來

這一步替換和用\(\mu * 1 = \epsilon\)替換有殊途同歸之妙，都是將\(gcd\)等於的限制弱化了，變成了整除的關係

推倒後獲得
\[ S(n,m) = \sum_{d\mid n}\varphi(\frac{n}{d})\cdot S(d, \lfloor \frac{m}{d} \rfloor) \]
對於n不是不一樣質因子的乘積的，根據\(\varphi\)的公式，多的質因子次數直接提出來乘上就好了

而後記憶化搜索，\(n=1\)就是\(\varphi\)的前綴和

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <map>
using namespace std;
typedef long long ll;
const int N=1664512, U=1664510, mo = 1e9+7;
inline int read(){
    char c=getchar(); int x=0,f=1;
    while(c<'0' || c>'9') {if(c=='-')f=-1; c=getchar();}
    while(c>='0' && c<='9') {x=x*10+c-'0'; c=getchar();}
    return x*f;
}

int n, m;
inline void mod(int &x) {if(x>=mo) x-=mo; else if(x<0) x+=mo;}
bool notp[N]; int p[N/10], phi[N], pr[N];
void sieve(int n) {
    phi[1]=1; pr[1] = 1;
    for(int i=2; i<=n; i++) {
        if(!notp[i]) p[++p[0]] = i, phi[i] = i-1, pr[i] = i;
        for(int j=1; j <= p[0] && i*p[j] <= n; j++) {
            int t = i*p[j];
            notp[ t ] = 1;
            if(i % p[j] == 0) {
                phi[t] = phi[i] * p[j];
                pr[t] = pr[i];
                break;
            }
            phi[t] = phi[i] * (p[j] - 1); 
            pr[t] = pr[i] * p[j];
        }
        mod(phi[i] += phi[i-1]);
    }
}

namespace ha {
    const int p = 1001001;
    struct meow{int ne, val, r;} e[3000];
    int cnt=1, h[p];
    inline void insert(int x, int val) {
        int u = x % p;
        for(int i=h[u];i;i=e[i].ne) if(e[i].r == x) return;
        e[++cnt] = (meow){h[u], val, x}; h[u] = cnt;
    }
    inline int quer(int x) {
        int u = x % p;
        for(int i=h[u];i;i=e[i].ne) if(e[i].r == x) return e[i].val;
        return -1;
    }
} using ha::insert; using ha::quer;

int dj_s(int n) {
    if(n <= U) return phi[n];
    if(quer(n) != -1) return quer(n);
    int ans = (ll) n * (n+1) / 2 %mo, r;
    for(int i=2; i<=n; i=r+1) {
        r = n/(n/i);
        mod(ans -= (ll) dj_s(n/i) * (r-i+1) %mo);
    }
    insert(n, ans);
    return ans;
}

inline int Pow(int a, int b) {
    int ans=1;
    for(; b; b>>=1, a=a*a)
        if(b&1) ans=ans*a;
    return ans;
}
inline ll Phi(int n) {
    int ans = 1;
    if(n <= U) {mod(ans = phi[n] - phi[n-1]); return ans;}
    int m = sqrt(n);
    for(int i=1; p[i] <= m; i++) if(n % p[i] == 0) {
        int a = 0;
        while(n % p[i] == 0) a++, n /= p[i];
        ans *= Pow(p[i], a-1) * (p[i] - 1);
    }
    return ans;
}

map<int, int> Map[N];
int S(int n, int m) {
    if(m == 0) return 0; 
    if(n == 1) return dj_s(m);
    if(Map[n][m]) return Map[n][m];
    //printf("S %d %d\n", n, m);
    int ans = 0;
    for(int i=1; i*i <= n; i++) if(n%i == 0) {
        int j = n/i;
        mod(ans += Phi(j) * S(i, m/i) %mo);
        if(i != j) mod(ans += Phi(i) * S(j, m/j) %mo);
    }
    Map[n][m]=ans;
    return ans;
}
int main() {
    freopen("in", "r", stdin);
    sieve(U);
    n=read(); m=read();
    int ans = 0;
    for(int i=1; i<=n; i++) mod(ans += (ll) i / pr[i] * S(pr[i], m) %mo);
    //for(int i=1; i<=n; i++) printf("nnnnnnnn %d  %d\n", i, S(i, m));
    printf("%d\n", ans);
}