uoj221. 【NOI2016】循环之美

考完后非常令人伤心的题。。。
打表&&找规律&&数学推导……..
我们得到了目标式子：

$$\Sigma_{i=1}^n\Sigma_{j=1}^m[(i,j)=1][(j,k)=1]$$
这里的（i，j）表示gcd（i，j）
然后：

方向一：

$$\Sigma_{i=1}^n\Sigma_{j=1}^m\textbf{[(i,j)=1]}[(j,k)=1]$$
$$\textbf{[(i,j)=1]}$$
思维难度较大：
可参考 here

方向二：

$$\Sigma_{i=1}^n\Sigma_{j=1}^m[(i,j)=1]\textbf{[(j,k)=1]}$$
$$\textbf{[(j,k)=1]}$$
思维难度很小的SB方法QAQ
记f（n,m,k）表示上式
展开原式得：
$\Sigma_{i=1}^n\Sigma_{j=1}^m[(i,j)=1]\Sigma_{t|(j,k)}\mu(t)$
$\Sigma_{t|(j,k)}\mu(t)\Sigma_{i=1}^n\Sigma_{j=1}^{[m/t]}[(i,j*t)=1]$
$\Sigma_{t|(j,k)}\mu(t)\Sigma_{i=1}^n\Sigma_{j=1}^{[m/t]}[(i,j)=1][(i,t)=1]$
$\Sigma_{t|(j,k)}\mu(t)f([m/t],n,t)$
注意边界：
当K=1时是广为人知的处理方法
及min（n，m）=0时

易分析复杂度$O（lognlogm\sqrt{n}+n^\frac{2}3）$

#include<bits/stdc++.h>
#define rep(i,k,n) for(int i=k;i<=n;i++)
using namespace std;
typedef long long ll;
const int N=1e6+7;
struct data{
    int n,m,k;
    data(int n=0,int m=0,int k=0):n(n),m(m),k(k){}
};
bool operator<(const data& a,const data& b){
    if(a.n==b.n)return a.m==b.m ? a.k<b.k : a.m<b.m;
    else return a.n<b.n; 
}
int mu[N],pri[N],top=0,vis[N],sum_mu[N],lim,shu[N],cnt=0;
map<int,int> sm;
void init(){
    lim=1e6;
    mu[1]=1;
    rep(i,2,lim){
        if(!vis[i]){
            pri[++top]=i;
            mu[i]=-1;
        }rep(j,1,top){
            int p=pri[j];
            if(1ll*p*i>lim)break;
            vis[p*i]=1;
            if(i%p)mu[i*p]=-mu[i];
            else{mu[i*p]=0;break;}
        }
    }rep(i,1,lim)sum_mu[i]=sum_mu[i-1] + mu[i];
}
int Du(int n){
    if(n<=lim)return sum_mu[n];
    else if(sm.count(n))return sm[n];
    else{
        int i=2,tmp=0;
        while(i<=n){
            int j=(n/(n/i));
            tmp+=(j-i+1)*Du(n/i);
            i=j+1;
        }
        return (sm[n]=1-tmp);
    }
}
map<data,ll> mp;
ll solve(int n,int m,int K){
    if(mp.count(data(n,m,K)))return mp[data(n,m,K)];
    else{
        if(min(n,m)==0)return 0;
        else if(K==1){
            if(n>m)swap(n,m);
            int i=1;
            ll tmp=0;
            while(i<=n){
                int j=min(n/(n/i),m/(m/i));
                ll res=1ll*(n/i)*(m/i);
                ll now=Du(j)-Du(i-1);
                tmp+=now*res;
                i=j+1;
            }
            return (mp[data(n,m,K)]=mp[data(m,n,K)]=tmp);
        }else{
            ll tmp=0;
            rep(i,1,cnt){
                if(shu[i]>K)break;
                if(K%shu[i]==0){
                    if(mu[shu[i]]){
                        tmp+=1ll*mu[shu[i]]*solve(m/shu[i],n,shu[i]);
                    }
                }
            }return (mp[data(n,m,K)]=tmp);
        }
    }
}
int n,m,K;
int main(){
    freopen("out.out","w",stdout);
    init();
    scanf("%d%d%d",&n,&m,&K);
    rep(i,1,K)if(K%i==0)shu[++cnt]=i;
    ll ans=solve(n,m,K);
    printf("%I64d\n",ans);
}

启示：
1.方向&&方法
2.递归结构出现的越早越好