数论函数变换

对于各种积性函数，都可以通过两种方法进行计算：

1. 分解质因数
2. 线性筛

那么，我们就可以得到一些常用的积性函数值。

typedef pair<int, int> P;
int prime[MAXN], cnt, phi[MAXN], mu[MAXN], mnsum[MAXN], a[MAXN];
bool isp[MAXN];
P f[MAXN];//约数和
void Euler(int n) {
    mu[1] = 1;
    phi[1] = 1;
    f[1] = P(1, 1);
    for(int i = 2; i <= n; i++) {
        f[i].second = i;
        if(!isp[i]) {
            prime[++cnt] = a[i] = i;
            phi[i] = i - 1;
            mu[i] = -1;
            f[i].first = mnsum[i] = i + 1;
        }
        for(int j = 1; j <= cnt && i * prime[j] <= n; j++) {
            isp[i * prime[j]] = 1;
            if(i % prime[j] == 0) {
                phi[i * prime[j]] = phi[i] * prime[j];
                mu[i * prime[j]] = 0;
                mnsum[i * prime[j]] = mnsum[i] + a[i] * prime[j];
                f[i * prime[j]].first = f[i].first / mnsum[i] * mnsum[i *
                    prime[j]]; a[i * prime[j]] = a[i] * prime[j];
                break;
            }
            phi[i * prime[j]] = phi[i] * (prime[j] - 1);
            mu[i * prime[j]] = -mu[i];
            f[i * prime[j]].first = f[i].first * (prime[j] + 1);
            a[i * prime[j]] = prime[j];
            mnsum[i * prime[j]] = prime[j] + 1;
        }
    }
}

---

先来一道简单题：bzoj2005
题目大意：求$\Sigma_i^n\Sigma_j^mgcd(i, j) * 2 - 1$
然后推一推式子
$\Sigma_i^n\Sigma_j^mgcd(i, j) * 2 - 1$
=$- n * m + 2 * \Sigma_i^n\Sigma_j^mgcd(i, j)$
=$- n * m + 2 * \Sigma_i^n\Sigma_j^m\Sigma_{d|gcd(i, j)} \phi(d)$
=$-n * m + 2 * \Sigma_d^{max(n, m)}\phi(d) * \lfloor\frac{n}{d}\rfloor\lfloor\frac{m}{d}\rfloor$
然后枚举d这道题就可以A掉了，复杂度为O(n)。

#include <cstdio>
#include <algorithm>
#include <cmath>
using namespace std;

const int MAXN = 100003;
int phi[MAXN], prime[MAXN], cnt;
bool isp[MAXN];
void Euler(int n) {
    phi[1] = 1;
    for(int i = 2; i <= n; i++) {
        if(!isp[i]) {
            prime[++cnt] = i;
            phi[i] = i - 1;
        }
        for(int j = 1; j <= cnt && i * prime[j] <= n; j++) {
            isp[i * prime[j]] = 1;
            if(i % prime[j] == 0) {
                phi[i * prime[j]] = phi[i] * prime[j];
                break;
            }
            else phi[i * prime[j]] = phi[i] * (prime[j] - 1);
        }
    }
}
long long ans;
int main() {
    int n, m, k;
    scanf("%d%d", &n, &m);
    k = max(n, m);
    Euler(k);
    ans -= (long long)n * m;
    for(int d = 1; d <= k; d++) {
        ans += 2LL * phi[d] * (n / d) * (m / d);
    }
    printf("%lld\n", ans);
    return 0;
}

---

来看另一道题：bzoj1101
题目大意：对于给定的整数N,M和d，有多少正整数对x,y，满足x<=N，y<=M，并且gcd(x,y)=d。
等价于x<=N/d，y<=M/d，互质的x，y的对数。
那么我们另n=N/d, m = M/d
于是原题就变成了求$\Sigma_i^n\Sigma_j^me(gcd(i, j))$
其中$e(x) = (x == 1) = \Sigma_{d | x} \mu(d)$
然后根据套路我们再来推一波式子
$\Sigma_i^n\Sigma_j^me(gcd(i, j))$
=$\Sigma_i^n\Sigma_j^m\Sigma_{d|gcd(i, j)} \mu(d)$
=$\Sigma_d^{min(n, m)}\mu(d) * \lfloor\frac{n}{d}\rfloor\lfloor\frac{m}{d}\rfloor$
根据上一题的做法，每次询问复杂度O(n)，显然不能胜任。
但是我们观察式子，显然 $\lfloor\frac{n}{d}\rfloor\lfloor\frac{m}{d}\rfloor$的取值只有$\sqrt{n} + \sqrt{m}$个，那么我们可以预处理$\mu$的前缀和，并分块完成。

#include <cstdio>
#include <algorithm>
using namespace std;

const int MAXN = 50003;
int mu[MAXN], prime[MAXN], cnt;
bool isp[MAXN];
void Euler(int n) {
    mu[1] = 1;
    for(int i = 2; i <= n; i++) {
        if(!isp[i]) {
            prime[++cnt] = i;
            mu[i] = -1;
        }
        for(int j = 1; j <= cnt && i * prime[j] <= n; j++) {
            isp[i * prime[j]] = 1;
            if(i % prime[j] == 0) {
                mu[i * prime[j]] = 0;
                break;
            }
            mu[i * prime[j]] = -mu[i];
        }
    }
}
int cal(int n, int m) {
    if(n > m) swap(n, m);
    int res = 0, last;
    for(int i = 1; i <= n; i = last + 1) {
        last = min(n / (n / i), m / (m / i));
        res += (n / i) * (m / i) * (mu[last] - mu[i - 1]);
    }
    return res;
}
int main() {
    Euler(50000);
    for(int i = 2; i <= 50000; i++) mu[i] += mu[i - 1];
    int T, a, b, d;
    scanf("%d", &T);
    while(T--) {
        scanf("%d%d%d", &a, &b, &d); a /= d; b /= d;
        printf("%d\n", cal(a, b));
    }
    return 0;
}

还有一题：bzoj3529
题目大意：另F(i)表示i的约数和，q次给定n,m,a，求$\Sigma_i^n\Sigma_{j and F(gcd(i, j)) <= a}^m F(gcd(i, j)) mod 2 ^ {31}$
有个a的限制，好像式子并不能推了啊。所以先把它去掉。
另$g(i) = \Sigma_x^n\Sigma_y^me(gcd(x, y) == i)$
那么可以得到$g(i) = \Sigma_{i | d}\mu(\frac{d}{i}) * \lfloor\frac{n}{d}\rfloor\lfloor\frac{m}{d}\rfloor$
于是$ans = \Sigma_i^{min(n, m)}F(i) * g(i)$
展开化简得到$\Sigma_d^{min(n, m)} \lfloor\frac{n}{d}\rfloor\lfloor\frac{m}{d}\rfloor\Sigma_{i | d}F(i)\mu(\frac{d}{i})$
然后对$\Sigma_{i | d}F(i)\mu(\frac{d}{i})$求个前缀和
枚举每个i更新i的倍数即可
现在考虑将a离线处理，询问按a排序
用树状数组维护前缀和即可。

#include <cstdio>
#include <algorithm>
using namespace std;

typedef pair<int, int> P;
const int MAXN = 100003;
int prime[MAXN], cnt, mu[MAXN], mnsum[MAXN], a[MAXN];
bool isp[MAXN];
P f[MAXN];
void Euler(int n) {
    mu[1] = 1;
    f[1] = P(1, 1);
    for(int i = 2; i <= n; i++) {
        f[i].second = i;
        if(!isp[i]) {
            prime[++cnt] = a[i] = i;
            mu[i] = -1;
            f[i].first = mnsum[i] = i + 1;
        }
        for(int j = 1; j <= cnt && i * prime[j] <= n; j++) {
            isp[i * prime[j]] = 1;
            if(i % prime[j] == 0) {
                mu[i * prime[j]] = 0;
                mnsum[i * prime[j]] = mnsum[i] + a[i] * prime[j];
                f[i * prime[j]].first = f[i].first / mnsum[i] * mnsum[i * prime[j]];
                a[i * prime[j]] = a[i] * prime[j];
                break;
            }
            mu[i * prime[j]] = -mu[i];
            f[i * prime[j]].first = f[i].first * (prime[j] + 1);
            a[i * prime[j]] = prime[j];
            mnsum[i * prime[j]] = prime[j] + 1;
        }
    }
}
int bit[MAXN];
void add(int x, int v) {
    for(; x <= 100000; x += x & -x) bit[x] += v;
}
int sum(int x) {
    int res = 0;
    for(; x; x -= x & -x) res += bit[x];
    return res;
}
struct Node {
    int n, m, a, id;
    bool operator < (const Node &x) const {
        return a < x.a;
    }
    inline void read(int i) {
        id = i; scanf("%d%d%d", &n, &m, &a);
    }
}q[MAXN];
int ans[MAXN];
int cal(int n, int m) {
    if(n > m) swap(n, m);
    int res = 0, last;
    for(int i = 1; i <= n; i = last + 1) {
        last = min(n / (n / i), m / (m / i));
        res += (n / i) * (m / i) * (sum(last) - sum(i - 1));
    }
    return res & 0x7fffffff;
}
int main() {
    Euler(100000);
    sort(f + 1, f + 100001);
    int T;
    scanf("%d", &T);
    for(int i = 1; i <= T; i++) q[i].read(i);
    sort(q + 1, q + T + 1);
    int cur = 1, N, M, A;
    for(int Q = 1; Q <= T; Q++) {
        N = q[Q].n, M = q[Q].m, A = q[Q].a;
        while(cur <= 100000 && f[cur].first <= A) {
            for(int i = f[cur].second; i <= 100000; i += f[cur].second)
                add(i, f[cur].first * mu[i / f[cur].second]);
            cur++;
        }
        ans[q[Q].id] = cal(N, M);
    }
    for(int i = 1; i <= T; i++) printf("%d\n", ans[i]);
    return 0;
}

来一道lcm的题：bzoj2154
题目大意：求$\Sigma_i^n\Sigma_j^mlcm(i, j) mod 20101009$
这里跟之前的区别在于，gcd(i,j)在分母上，那么考虑反演
另$f(x) = \frac{1}{x}$
$F(x) = f$ X $\mu$ 这里X是Dirichlet积
即$F(x) = \Sigma_{d | x}f(d) * \mu(\frac{x}{d})$
那么由莫比乌斯反演我们有$f=F$X$1$
那么$\Sigma_i^n\Sigma_j^mlcm(i, j)$
=$\Sigma_i^n\Sigma_j^mi*j*f(gcd(i, j))$
=$\Sigma_i^n\Sigma_j^mi*j*\Sigma_{d|gcd(i, j)F(d)}$
另$sum(x, y) = \Sigma_i^x\Sigma_j^yi*j$
于是可以进一步化简$ans = \Sigma_d^nF(d) * d * d * sum(\frac{n}{d}, \frac{m}{d})$

#include <iostream>
#include <cstdio>
#include <algorithm>
using namespace std;

typedef long long ll;
const int mod = 20101009, MAXN = 10000003;
int mu[MAXN], prime[MAXN], sum[MAXN];
int cnt;
bool isp[MAXN];
void getmu(int n) {
    mu[1] = 1;
    for(int i = 2; i <= n; i++) {
        if(!isp[i]) {
            mu[i] = -1;
            prime[++cnt] = i;
        }
        for(int j = 1; j <= cnt && i * prime[j] <= n; j++) {
            isp[i * prime[j]] = 1;
            if(i % prime[j] == 0) {
                mu[i * prime[j]] = 0;
                break;
            }
            mu[i * prime[j]] = -mu[i];
        }
    }
}

ll n, m, ans;

ll query(ll x, ll y) { return (x * (x + 1) / 2 % mod) * (y * (y + 1) / 2 % mod) % mod; }

ll F(ll x, ll y) {
    ll res = 0, last;
    for(ll i = 1; i <= min(x, y); i = last + 1) {
        last = min(x / (x / i), y / (y / i));
        res = (res + (sum[last] - sum[i - 1]) * query(x / i, y / i) % mod) % mod;
    }
    return res;
}

int main() {
    cin>>n>>m;
    getmu(min(n, m));
    for(ll i = 1; i <= min(n, m); i++) sum[i] = (sum[i - 1] + (i * i * mu[i]) % mod) % mod;
    ll last;
    for(ll d = 1; d <= min(n, m); d = last + 1) {
        last = min(n / (n / d), m / (m / d));
        ans = (ans + (last - d + 1) * (d + last) / 2 % mod * F(n / d, m / d) % mod) % mod;
    }
    ans = (ans + mod) % mod;
    cout<<ans<<endl;
    return 0;
}

P.S.:推式子治好了我多年的公式恐惧症啊
总结一下：当推式子陷入僵局的时候，尝试将一些东西处理出来，通过分块、反演等方法降低复杂度。对于式子本身的特性进行观察，有时也可以发现优化的地方。
在推式子过程中，将有关变量放在一块处理是化简的常用技巧。