CSU 1325 A very hard problem (莫比乌斯反演＋分块求和优化)

 A very hard problem

Time Limit: 3 Sec   Memory Limit: 160 MB
Submit: 266   Solved: 63
[ Submit][ Status][ Web Board]

Description

CX老湿经常被人黑，被黑得多了，自己也就麻木了。于是经常听到有人黑他，他都会深情地说一句：禽兽啊！

一天CX老湿突发奇想，给大家出了一个难题，并且声称谁能够准确地回答出问题才能继续黑他，否则他就要反击了。

这个难题就是：

给出两个数p和q，接下来q个询问，每个询问给出两个数A和B，请分别求出：

一、有多少个有序数对（x，y）满足1<=x<=A，1<=y<=B，并且gcd（x，y）为p的一个约数；

二、有多少个有序数对（x，y）满足1<=x<=A，1<=y<=B，并且gcd（x，y）为p的一个倍数。

Input

只有一组测试数据。

第一行两个数：p和q。（1<p<10^7 ，1<q<1000。）

接下来有q行，每行两个数A和B。（1<A,B<10^7）

Output

输出共q行。每行两个数。用空格隔开。

分别表示题目描述中的两个对应的答案。

(x,y)=(2,3)和(x,y)=(3,2)被视为两个不同有序数对哦！

Sample Input

6 3
8 8
15 32
13 77

Sample Output

58 1
423 10
883 24

HINT

 对于64位整型请用lld，或者cin，cout。T_T

CSU_LQ

Source

CSU Monthly 2013 Oct.

题目链接：http://acm.csu.edu.cn/OnlineJudge/problem.php?id=1325

题目分析：倍数直接(a / p) * (b / p)即可，不解释，主要是约数，开始枚举因子做，各种超时，这题有点类似 HDU 4746，考虑预处理出每个因子对答案的贡献，这里贡献指的是它的莫比乌斯函数的系数和，举个的例子如果p＝6，它的因子有1，2，3，6，设F(x) = (a / x) * (b / x)，x为最大公约数，f(x)为x为最大公约数时的组数

f(1) = u(1)F(1 * 1) + u(2)F(1 * 2) + u(3)F(1 * 3) + u(4)F(1 * 4) + u(5)F(1 * 5) + u(6)F(1 * 6) + ...

f(2) = u(1)F(2 * 1) + u(2)F(2 * 2) + u(3)F(2 * 3) + ...

f(3) = u(1)F(3 * 1) + u(2)F(3 * 2) + ...

f(6) = u(1)F(6 * 1) + ...

ans = f(1) + f(2) + f(3) + f(6)，我们需要预处理出F函数前面的系数，显然F(6)的系数为u(6) + u(3) + u(2) + u(1)，这些都是它的约数，因此预处理时通过类似筛法的nlogn复杂度即可完成对fac_msum数组的预处理，fac_msum[i]表示，最终答案里F(i)前面的系数

光这样还是不够的，依然超时，必须再加上分块求和优化才可通过本题，总复杂度约为p + plog(p) + q*sqrt(p)，吐槽一句，单组样例，1e7的数据，能不memset尽量不memset。。。

#include <cstdio>
#include <cstring>
#include <algorithm>
#include <vector>
#define ll long long
using namespace std;
int const MAX = 1e7 + 5;
int p, q, pnum;
short mob[MAX];
bool noprime[MAX];
int pr[MAX], sum[MAX], fac_msum[MAX];
int fac[1000], facnum;

void Mobius()
{
    pnum = 0;
    mob[1] = 1;
    for(int i = 2; i < MAX; i++)
    {
        if(!noprime[i])
        {
            pr[pnum ++] = i;
            mob[i] = -1;
        }
        for(int j = 0; j < pnum && i * pr[j] < MAX; j++)
        {
            noprime[i * pr[j]] = true;
            if(i % pr[j] == 0)
            {
                mob[i * pr[j]] = 0;
                break;
            }
            mob[i * pr[j]] = -mob[i];
        }
    }
}

void pre()
{
    for(int i = 1; i * i <= p; i++)
    {
        if(p % i == 0)
        {   
            fac[facnum ++] = i;
            if(i * i != p)
                fac[facnum ++] = p / i;
        }
    }
    for(int i = 0; i < facnum; i++)
        for(int j = 1; j * fac[i] < MAX; j++)
            fac_msum[j * fac[i]] += mob[j];
    for(int i = 1; i < MAX; i++)
        sum[i] = sum[i - 1] + fac_msum[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) (l / i) * (r / i) * (sum[last] - sum[i - 1]);  
    }
    return ans;
}   

int main()
{
    Mobius();
    scanf("%d %d", &p, &q);
    pre();
    while(q --)
    {
        int a, b;
        facnum = 0;
        ll ans = 0;
        scanf("%d %d", &a, &b);
        printf("%lld %lld\n", cal(a, b), (ll) (a / p) * (b / p));
    }   
}