约数个数和(数论，莫比乌斯反演)BZOJ3994

题目

3994: [SDOI2015]约数个数和

Time Limit: 20 Sec Memory Limit: 128 MB
Submit: 492 Solved: 319
[Submit][Status][Discuss]
Description

设d(x)为x的约数个数，给定N、M，求

$$\sum_{i=1}^N\sum_{j=1}^Md(ij)$$
Input

输入文件包含多组测试数据。
第一行，一个整数T，表示测试数据的组数。
接下来的T行，每行两个整数N、M。
Output

T行，每行一个整数，表示你所求的答案。
Sample Input

2
7 4
5 6
Sample Output

110
121

思路：
公式
$d(mn)=\sum_{i|n}{\sum_{j|m}[gcd(i,j)==1]}$
可以归纳证明：
显然$m=n=1$是成立。假设当$m',n'$都不含有素数$p$时，上式成立。
当$m,n$含有$p$的个数分别为$k_1,k_2$，即$p^k_1|m,p^k_2|n$时；
$d(mn)=(k1+k2+1)d(m'n')$;
同时式子右边中$i,j$在原来的基础上可以含有$p$的个数情况为:$(1,k_1),(1,k_1-1),...(1,1),...(k_2,1)$,共有$(k_1+k_2+1)$个。所以也是$(k_1+k_2+1)$倍。结论成立。

上式的前缀和就是所求。

然后反演一下就有：

$$ans=\sum_{i=1}^n\sum_{j=1}^m\sum_{d|i,d|j}{\mu(d)\lfloor\frac ni\rfloor\lfloor\frac mj\rfloor}\\=\sum_{d=1}^n\sum_{i=1}^{n|d}\sum_{j=1}^{n|d}\mu(d)\lfloor\frac n{id}\rfloor\lfloor\frac m{jd}\rfloor\\=\sum_{d=1}^n\mu(d)\left(\sum_{i=1}^{b|n}\lfloor\frac n{id}\rfloor\times\sum_{j=1}^{m|d}\lfloor\frac m{jd}\rfloor\right)\\=\sum_{d=1}^n\mu(d)f(\lfloor\frac nd\rfloor)f(\lfloor\frac md\rfloor)$$

这里的$f(x)$实际上就是$x$的约数个数的前缀和，然后约数个数和$\mu$函数可以打表。

/**************************************************************
    Problem: 3994
    User: ocgcn
    Language: C++
    Result: Accepted
    Time:2984 ms
    Memory:2248 kb
****************************************************************/

#include<iostream>
#include<cmath>
#include<string.h>
#include<algorithm>
#include<stdio.h>
using namespace std;
#define REP(i,a,b)  for(int i=a;i<b;++i)
#define LL      long long
#define mset(a,b)   memset(a,b,sizeof a)
#define scan(n) scanf("%d",&n)
const int maxn = 50020;
const int N = 1003;
const int mod = 1000000007;
const int INF = 0x3f3f3f3f;

int mu[maxn], prime[maxn], c[maxn];
LL f[maxn];
void init_mu()
{
    mset(prime, 0);
    mu[1] = f[1] = 1;
    REP(i, 2, maxn)
    {
        if (!prime[i])
        {
            prime[++prime[0]] = i;
            mu[i] = -1; f[i] = 2; c[i] = 1;
        }
        for (int j = 1; j <= prime[0] && prime[j] * i < maxn; ++j)
        {
            prime[i*prime[j]] = 1;
            if (i%prime[j] == 0)
            {
                mu[i*prime[j]] = 0;
                c[i*prime[j]] = c[i] + 1;
                f[i*prime[j]] = f[i] / (c[i] + 1)*(c[i] + 2);
                break;
            }
            mu[i*prime[j]] = -mu[i];
            f[i*prime[j]] = f[i] * 2;
            c[i*prime[j]] = 1;

        }

    }
    REP(i, 2, maxn) f[i] += f[i - 1], mu[i] += mu[i - 1];
}


int main()
{
    init_mu();
    int t, m, n;
    scan(t);
    while (t--)
    {
        scan(m);
        scan(n);
        int mn = min(m, n);
        LL ans = 0;
        for (int i = 1, j; i <= mn; i = j + 1)
        {
            j = min(n / (n / i), m / (m / i));
            ans += (mu[j] - mu[i - 1])*f[n / i] * f[m / i];
        }
        printf("%lld\n",ans);
    }
    return 0;
}