[SDOI2015]约数个数和

题目描述

设d(x)为x的约数个数，给定N、M，求 ∑_i∑_j d(i*j)  (i<=n,j<=m)

输入输出格式

输入格式：

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

输出格式：

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

输入输出样例

输入样例#1：

2
7 4
5 6

输出样例#1：

110
121

说明

1<=N, M<=50000

1<=T<=50000

题解：

已知推论：d(i*j)=∑_a∑_b[gcd(a,b)=1]  (a|i , b|j)

原式化为

∑_i∑_j∑_a∑_b[gcd(a,b)=1]   (1<=i<=n&&1<=j<=m)

=∑_i∑_j∑_a∑_b∑_dμ(d)       (d|a,d|b,d|i,d|j)

=∑_i∑_j∑_dμ(d)∑_a∑_b1     (d|i,d|j,1<=a<=n/i,1<=b<=m/j)   

=∑_i∑_j∑_dμ(d)*[n/i][m/j]

=∑_dμ(d)∑_i∑_j[n/id][m/jd]  (1<=d<=min(n,m),i<=n/d,j<=m/d)

=∑_dμ(d)*∑_i[n/id]*∑_j[m/jd]=∑_dμ(d)*∑_i[[n/d]/i]*∑_j[[m/d]/j]

令f(x)=∑[x/i] (1<=i<=x) 

接下来用分块的方法在O(N√N)时间求出f(1~n)

分块详见我写的http://www.cnblogs.com/Y-E-T-I/p/7255421.html

每一次询问的时候也用分块的方法求ans=∑_dμ(d)*[f(n/d)]*[f(m/d)]

不过第三行似乎是多余的

这个结论要记下来

延伸结论

 

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
using namespace std;
bool vis[50001];
int tot,prime[50001],mu[50001],n,m;
long long ans,d[50001];
void mobius()
{int i,j;
mu[1]=1;
    for (i=2;i<=50000;i++)
    {
        if (vis[i]==0)
        {
            tot++;
            prime[tot]=i;
            mu[i]=-1;
        }
        for (j=1;j<=tot,i*prime[j]<=50000;j++)
         {
             vis[i*prime[j]]=1;
             if (i%prime[j]==0)
             {
                 mu[i*prime[j]]=0;
                 break;
            }
                 mu[i*prime[j]]=-mu[i];
         }
    }
    for (i=1;i<=50000;i++)
    mu[i]+=mu[i-1];
}
void solve(int x)
{int i;
    int pos=1;
    for (i=1;i<=x;i=pos+1)
    {
        pos=x/(x/i);
         d[x]+=(pos-i+1)*(x/i);
    }
}
void work()
{int i;
    ans=0;
    int pos=1;
    for (i=1;i<=n;i=pos+1)
    {
        pos=min(n/(n/i),m/(m/i));
        ans+=(mu[pos]-mu[i-1])*(d[n/i])*(d[m/i]);
    }
}
int main()
{int i,T;
    mobius();
     for (i=1;i<=50000;i++)
     solve(i);
      cin>>T;
      while (T--)
      {
        scanf("%d%d",&n,&m);
        if (n>m) swap(n,m);
         work();
         printf("%lld\n",ans);
      }
}