BZOJ 2671 Calc 数论

题目大意：给定$N$，求$\sum_{i=1}^n\sum_{j=1}^n[i+j|ij]1$
跪Nodgd= =
不妨设$d=\gcd(i,j),i=ad,j=bd,\gcd(a,b)=1$，那么有
$(a+b)d|abd^2$
即
$a+b|abd$
$\because\gcd(a,b)=1$
$\therefore\gcd(a+b,a)=1,\gcd(a+b,b)=1$
$\therefore a+b|d$
不妨设$d=(a+b)t$，那么我们求的就是三元组$(a,b,t)$的个数，其中满足$a<b,gcd(a,b)=1,b(a+b)t\leq N$

那么如果我们的$a$和$b$确定了，满足要求的$t$显然有$\lfloor\frac{N}{b(a+b)}\rfloor$个
由于$b\leq\lfloor\sqrt N-1\rfloor$，我们不妨枚举$b$
那么有
$ans=\sum_{b=1}^{\lfloor\sqrt N-1\rfloor}\sum_{a=1}^{b-1}[gcd(a,b)=1]\lfloor\frac{N}{b(a+b)}\rfloor$
对于一个确定的$b$，$\lfloor\frac{N}{b(a+b)}\rfloor$最多只有$2\sqrt N$个取值，我们可以分段处理
对于一段$[lower,upper]$，我们需要求出这一段中与$b$互质的数的个数
由容斥原理易知这样的数有$\sum_{k|b}\mu(k)\lfloor\frac{upper}k-\frac{lower-1}k\rfloor$个，枚举$b$的约数搞一搞就好了
复杂度是啥？
别管了暴力出奇迹就是了

#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
#define M 50500
using namespace std;

long long n;

long long stack[M],top;

int mu[M],prime[M],tot;
bool not_prime[M];
void Linear_Shaker()
{
    int i,j;
    mu[1]=1;
    for(i=2;i<=50000;i++)
    {
        if(!not_prime[i])
        {
            prime[++tot]=i;
            mu[i]=-1;
        }
        for(j=1;prime[j]*i<=50000;j++)
        {
            not_prime[prime[j]*i]=true;
            if(i%prime[j]==0)
            {
                mu[prime[j]*i]=0;
                break;
            }
            mu[prime[j]*i]=-mu[i];
        }
    }
}
void Get_Divisor(long long n)
{
    long long i;
    top=0;
    for(i=1;i*i<n;i++)
        if(n%i==0)
            stack[++top]=i,stack[++top]=n/i;
    if(i*i==n)
        stack[++top]=i;
    sort(stack+1,stack+top+1);
}
long long Calculate()
{
    long long a,b,k,last,re=0;
    for(b=1;b*(b+1)<=n;b++)
    {
        Get_Divisor(b);
        for(a=1;a<b&&b*(a+b)<=n;a=last+1)
        {
            last=min(n/(n/b/(a+b))/b-b,b-1);
            long long cnt=0;
            for(k=1;stack[k]<=last;k++)
                cnt+=mu[stack[k]]*(last/stack[k]-(a-1)/stack[k]);
            re+=n/b/(a+b)*cnt;
        }
    }
    return re;
}
int main()
{
    cin>>n;
    Linear_Shaker();
    cout<<Calculate()<<endl;
}