acdream 1148 GCD SUM 莫比乌斯反演 ansx,ansy

GCD SUM

Time Limit: 8000/4000MS (Java/Others)Memory Limit: 128000/64000KB (Java/Others)

SubmitStatisticNext Problem

Problem Description

给出N,M
执行如下程序：
long long  ans = 0,ansx = 0,ansy = 0;
for(int i = 1; i <= N; i ++)
   for(int j = 1; j <= M; j ++)
       if(gcd(i,j) == 1) ans ++,ansx += i,ansy += j;
cout << ans << " " << ansx << " " << ansy << endl;

Input

多组数据，每行两个数N,M(1 <= N,M <= 100000)。

Output

如题所描述，每行输出3个数，ans,ansx,ansy，空格隔开

Sample Input

5 5
1 3

Sample Output

19 55 55
3 3 6

Hint

总数据小于50000

 

对于第一个数字比较容易，用莫比乌斯反演直接就能做。

对于ansx，和ansy是同类问题。现在讨论ansx的做法.

我们设f(d) 代表1<=x<=N,1<=y<=M,gcd(x,y)=d 时，x的求和.

我们设F(d) 代表 1<=x<=N,1<=y<=M,gcd(x,y)为d的倍数时，x的求和

所以F(1) = f(1)+f(2)+f(3)......f(min(N,M))

     F(2) = f(2)+f(4)+f(6)....f(min(N/2,M/2));

     F(i) = f(i)+f(i*2)+f(i*3)....f(min(N/i,M/i));

 

由此我们可得 

由于d = 1 ；所以

对于的反演为

然后可以对mu[i]*i进行筛选，求前N项和。用分块来完成。

代码：

#include<iostream>
#include<stdio.h>
#include<cstring>
#include<cstdlib>
using namespace std;
typedef long long LL;

const int maxn = 1e5+3;
bool s[maxn];
int prime[maxn],len = 0;
int mu[maxn];
LL hxl [maxn];
int sum1[maxn];
void  init()
{
    memset(s,true,sizeof(s));
    mu[1] = 1;
    for(int i=2;i<maxn;i++)
    {
        if(s[i] == true)
        {
            prime[++len]  = i;
            mu[i] = -1;
        }
        for(int j=1;j<=len && (long long)prime[j]*i<maxn;j++)
        {
            s[i*prime[j]] = false;
            if(i%prime[j]!=0)
                mu[i*prime[j]] = -mu[i];
            else
            {
                mu[i*prime[j]] = 0;
                break;
            }
        }
    }
    for(int i=1;i<maxn;i++)
        sum1[i] = sum1[i-1]+mu[i];
    hxl[1] = mu[1];
    for(int i=2;i<maxn;i++){
        hxl[i] = i*mu[i]+hxl[i-1];
    }
}
int main()
{
    init();
    int n,m;
    while(scanf("%d%d",&n,&m)>0)
    {
        LL sum = 0;
        LL ansi = 0,ansj = 0;
        int a = n;
        int b = m;
        if(a>b) swap(a,b);
        for(int i=1,la = 0;i<=a;i++,i = la+1)
        {
            la = min(a/(a/i),b/(b/i));
            sum = sum + ((LL)(a/i))*(b/i)*(sum1[la]-sum1[i-1]);
            ansi = ansi +(hxl[la]-hxl[i-1])*(((LL)(n/i+1)*(n/i))/2)*(m/i);
            ansj = ansj +(hxl[la]-hxl[i-1])*(((LL)(m/i+1)*(m/i))/2)*(n/i);
        }
        printf("%lld %lld %lld
",sum,ansi,ansj);
    }
    return 0;
}