51nod1223 分数等式的数量

1. 【清华冬令营模拟12.28】和与积 (Standard IO)
   Time Limits: 2000 ms Memory Limits: 262144 KB Detailed Limits

Description

Input

Output

Sample Input

2
2
15

Sample Output

0
4

Data Constraint

Hint

题解

题目要我们求的是 a< b< n 且 a+b|ab的合法的(a,b)的对数

考虑转化，设d=gcd(a,b) ，设$a'=a/d ,b'=b/d$则 $d(a'+b')| a'b'd^2$

即 $a'+b'|a'b'd$

又$gcd(a',b')=1$
那么 $gcd(a',a'+b')=gcd(a'+b',a')=1$
则 $a'+b'|d$

由于不知道怎么打下取整，下面的除法都当成是套了下取整的吧~
发现d的取值范围是 1—-n/b，那么使得$a'+b'|d$的取值就有 $\frac{n}{b'(a'+b')}$种

$ans=\sum_{i=2}^{n}\sum_{j=1}^{i-1} [gcd(i,j)=1]\frac{n}{i*(i+j)}$

当i=j时对gcd不为1所以我们可以补上来，同时我们考虑补上i=1，最后再减掉
下面开始我们的n就是$\sqrt{n}了（为了方便嘛）$
$ans'=\sum_{i=1}^{n}\sum_{j=1}^{i} [gcd(i,j)=1]\frac{n}{i*(i+j)}$

设$f[k]=\sum_{i=1}^{n}\sum_{j=1}^{i} [gcd(i,j)=k]\frac{n}{i*(i+j)}$
那么$ans'=f[1]$

考虑反演

设$g[d]=\sum_{k=1}^{n/d}f[kd]$
那么$g[d]=\sum_{i=1}^{n/d}\sum_{j=1}^{i} \frac{n}{i*(i+j)*d*d}$

$f[d]=\sum_{k=1}^{n/d}μ(k)*g[k*d]$
$f[1]=\sum_{k=1}^{n}μ(k)\sum_{i=1}^{n/k}\sum_{j=1}^{i} \frac{n}{(i+j)*i*k*k}$

其中i+j可以分块

时间复杂度不会算QvQ

贴代码

#include<iostream>
#include<algorithm>
#include<cstdio>
#include<cstring>
#include<cmath>
#define fo(i,a,b) for(i=a;i<=b;i++)
#define ll long long
using namespace std;

const int maxn=36130;

ll mu[maxn],st[maxn];
bool bz[maxn];
ll i,j,k,l,m,n,x,y,z,o,tp,ans,chu,zi,be,ed;

void yu(){
    mu[1]=1;
    fo(i,2,36123){
        if (bz[i]==false){
            st[++tp]=i; mu[i]=-1;
        }
        fo(j,1,tp){
            if (st[j]*i>36122) break;
            bz[st[j]*i]=true;
            if (i%st[j]==0){
                mu[st[j]*i]=0;
                break;
            }
            mu[st[j]*i]=mu[i]*(-1);
        }
    }
}
int main(){
//  freopen("t3.in","r",stdin);
    scanf("%d",&o);
    yu();
    while (o){
        o--;
        scanf("%d",&n);
        x=trunc(sqrt(n));
        ans=0;
        fo(i,1,x){
            if (mu[i]==0) continue;
            fo(j,1,x/i){
                chu=i*i*j*(j+1);
                be=j+1; zi=n/chu;
                while (be<=j+j){
                    if (zi==0) break;
                    ed=(n/i/i/j)/zi; if (ed>j+j) ed=j+j;
                    ans=ans+mu[i]*zi*(ed-be+1);
                    be=ed+1; chu=i*i*j*be; zi=n/chu;
                }
            }
/*          fo(j,1,x/i) 
                fo(k,1,j) ans+=mu[i]*(n/((j+k)*j*i*i));*/
        }
        printf("%lld\n",ans-n/2);
    }
    return 0;
}