「LibreOJ NOI Round #1」动态几何问题

题链

SOL

  SKLCD​​ 为整数等价于 ab为完全平方数。

  答案 =∑x=1min(N,M)μ2(x)⋅⌊⌊N/x⌋⌋⋅⌊⌊M/x⌋⌋=sum_{x=1}^{min(N, M)} mu^2(x) cdot iglfloorsqrt{lfloor N / x floor}ig floor cdot iglfloorsqrt{lfloor M / x floor}ig floor=∑​x=1​min(N,M)​​μ​2​​(x)⋅⌊√​⌊N/x⌋​​​⌋⋅⌊√​⌊M/x⌋​​​

  ∑​x=1​n​​μ​2​​(x)=∑​i=1​⌊√​n​​​⌋​​μ(i)⋅⌊​i​2​​​​n​​⌋

  然后搞一搞

#pragma GCC optimize("-O3")
#include<bits/stdc++.h>
#define LL long long
#define N 122474489
using namespace std;
inline LL min(LL x,LL y){return (x>y)?y:x;}
LL anw;
#define sight(c) ('0'<=c&&c<='9')
inline void read(LL &x){
    static char c;
    for (c=getchar();!sight(c);c=getchar());
    for (x=0;sight(c);c=getchar())x=x*10+c-48;
}
inline LL dos(LL x,LL y){
    anw=sqrt(x/y); return x/(anw*anw);
}
char u[N];
int tot;bool usd[N]; 
short sum1[N];
int sum2[N];
int pim[N/10];
LL n,m,Last,now,ans;
void pre() {
    u[1]=1;  sum1[1]=sum2[1]=1;
    for (int i=2;i<N;i++) {
        if (!usd[i]) pim[++tot]=i,u[i]=-1;
        for (int j=1;j<=tot&&i*pim[j]<N;j++) {
            usd[i*pim[j]]=1;
            if (!(i%pim[j])) {
                u[i*pim[j]]=0; break;
            }
            u[i*pim[j]]=-u[i]; 
        } sum1[i]=sum1[i-1]+u[i];
        sum2[i]=sum2[i-1]+!!u[i];
    }
}
LL pro,X;
LL dos(LL x) {
    if (x<N) return sum2[x];
    pro=0; X=cbrt(x);
    for (int i=X;i;i--) pro+=x/((LL)i*i)*u[i];
    for (LL i=x/((X+1)*(X+1)),l=X+1,r;i;--i,l=r+1) 
       r=sqrt(x/i),pro+=(sum1[r]-sum1[l-1])*i;
    return pro;
}
int main () {
//    freopen("a.in","r",stdin);
    pre();
    read(n); read(m);
    if (n>m) swap(n,m);
    Last=0;
    for (LL i=1,last;i<=n;i=last+1) {
        last=min(dos(n,i),dos(m,i));
        now=dos(last);
        ans+=(now-Last)*((int)sqrtl(n/i))*((int)sqrtl(m/i)); Last=now;
    }
    printf("%lld
",ans); 
    return 0;
}