BZOJ 2671 Calc 数论

题目大意：给定N，求∑ni=1∑nj=1[i+j|ij]1
跪Nodgd= =
不妨设d=gcd(i,j),i=ad,j=bd,gcd(a,b)=1，那么有
(a+b)d|abd2
即
a+b|abd
∵gcd(a,b)=1
∴gcd(a+b,a)=1,gcd(a+b,b)=1
∴a+b|d
不妨设d=(a+b)t，那么我们求的就是三元组(a,b,t)的个数，其中满足a<b,gcd(a,b)=1,b(a+b)t≤N

那么如果我们的a和b确定了，满足要求的t显然有⌊Nb(a+b)⌋个
由于b≤⌊N−−√−1⌋，我们不妨枚举b
那么有
ans=∑⌊N√−1⌋b=1∑b−1a=1[gcd(a,b)=1]⌊Nb(a+b)⌋
对于一个确定的b，⌊Nb(a+b)⌋最多只有2N−−√个取值，我们可以分段处理
对于一段[lower,upper]，我们需要求出这一段中与b互质的数的个数
由容斥原理易知这样的数有∑k|bμ(k)(⌊upperk⌋−⌊lower−1k⌋)个，枚举b<script type="math/tex" id="MathJax-Element-125">b</script>的约数搞一搞就好了
复杂度是啥？
别管了暴力出奇迹就是了

#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;
}