BZOJ 3994 [SDOI2015]约数个数和 (神定理＋莫比乌斯反演)

3994: [SDOI2015]约数个数和

Time Limit: 20 Sec  Memory Limit:128 MB
Submit: 239  Solved: 176
[Submit][Status][Discuss]

Description

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

Input

输入文件包含多组测试数据。

第一行，一个整数T，表示测试数据的组数。

接下来的T行，每行两个整数N、M。

Output

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

Sample Input

2
7 4
5 6

Sample Output

110
121

HINT

 1<=N, M<=50000

1<=T<=50000

Source

Round 1 感谢yts1999上传

题目链接：http://www.lydsy.com/JudgeOnline/problem.php?id=3994

题目分析：和上一题类似，还是见公式： ∑ni=1∑mj=1d(ij)=∑ni=1∑mj=1⌊ni⌋⌊mj⌋[gcd(i,j)==1] ,

于是得到：

Ans=∑i=1n∑j=1m∑d|i,d|jμ(d)∗⌊n/i⌋⌊m/j⌋

变形得到

∑d=1min(n,m)μ(d)∑i=1n/d∑j=1m/d⌊nid⌋⌊mjd⌋  

继续变形

∑d=1min(n,m)μ(d)∑i=1n/d⌊ndi⌋∑j=1m/d⌊mdj⌋

令

∑i=1n/d⌊ndi⌋=f(⌊nd⌋)

公式变成

∑d=1min(n,m)μ(d)f(⌊nd⌋)f(⌊md⌋)

现在的问题是怎样计算f的值，f[i]表示的是 ∑ik=1⌊ik⌋ ，我们先考虑一个函数f'[i]，设f'[i]为i的约数的个数，那么显然f'[i]为积性函数，我们可以用线性筛得到，仔细观察可以发现f[i] = f'[1] + f'[2] + ... + f'[i]，比如f[6] = f'[1]+f'[2]+...+f'[6] = 1+2+2+3+2+4 = 6+3+2+1+1+1 = 14，下面简单解释一下线性筛是如何得到一个数约数的个数的，首先根据约数个数定理：对于一个大于1的正整数n可以分解质因数：n=p1^a1*p2^a2*p3^a3*…*pk^ak，则n的正约数有(a₁+1)(a₂+1)(a₃+1)…(ak+1)个，又我们知道线性筛筛质数每次筛的都是用最小的质因子去筛，因此我们可以记录这个最小质因子唯一分解后的次幂然后通过上面的公式求解，详细见程序注释。这样的话问题就解决了，同时把u和f都筛出来，剩下的分块求和即可

#include <cstdio>
#include <cstring>
#include <algorithm>
#define ll long long
using namespace std;
int const MAX = 50005;
//mob表示莫比乌斯函数，sum表示莫比乌斯函数前缀和，p表示素数
int mob[MAX], sum[MAX], p[MAX];
//facnum表示约数个数，f表示约数个数前缀和，d表示最小质因子的次幂
int facnum[MAX], f[MAX], d[MAX]; 
//prime是用来筛素数
bool noprime[MAX];     //前缀和数组其实可以直接由一个数组得到，这样申明只是为了让意思更清楚

void Mobius()
{
    int pnum = 0;
    mob[1] = 1;
    sum[1] = 1;
    f[1] = 1;
    facnum[1] = 1;
    for(int i = 2; i < MAX; i++)
    {
        if(!noprime[i])
        {
            p[pnum ++] = i;
            mob[i] = -1;    
            facnum[i] = 2;  //素数的因子只有本身和1
            d[i] = 1;       //素数的最小质因子的次幂显然为1
        }
        for(int j = 0; j < pnum && i * p[j] < MAX; j++)
        {
            noprime[i * p[j]] = true;
            if(i % p[j] == 0)
            {
                mob[i * p[j]] = 0;
                //这里当前最小质因子的次幂加了1，因为i是p[j]的倍数，又乘了p[j]
                facnum[i * p[j]] = facnum[i] / (d[i] + 1) * (d[i] + 2); 
                d[i * p[j]] = d[i] + 1; 
                break;
            }
            mob[i * p[j]] = -mob[i];
            //facnum[i * p[j]] = facnum[i] * facnum[p[j]]
            //积性函数的性质i和p[j]显然互质，又facnum[p[j]] = 2，素数只有两个因子
            facnum[i * p[j]] = facnum[i] * 2;
            //此时当前最小质因子的次数为1
            d[i * p[j]] = 1; 
        }
        sum[i] = sum[i - 1] + mob[i];
        f[i] = f[i - 1] + facnum[i];
    }
}

ll cal(int l, int r)
{
    ll ans = 0;
    if(l > r)
        swap(l, r);
    for(int i = 1, last = 0; i <= l; i = last + 1)
    {
        last = min(l / (l / i), r / (r / i));
        ans += (ll) f[l / i] * f[r / i] * (sum[last] - sum[i - 1]);
    }
    return ans;
}

int main()
{
    Mobius();
    int T;
    scanf("%d", &T);
    while(T--)
    {
        int n, m;
        scanf("%d %d", &n, &m);
        printf("%lld\n", cal(n, m));
    }
}