loj#572. 「LibreOJ Round #11」Misaka Network 与求和

胡策题，迫使懒得点技能树选手学算法。。。

至少杜教和min_25会了

不想写小结啊啊啊啊

那么对于这个题先上个莫反套路

然后用杜教乘个I把f*u给消掉，发现h就是f，这个可以用类似min_25的思想搞。。。

#include<cstdio>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<algorithm>
#include<cmath>
#include<map>
#include<ctime>
using namespace std;
typedef long long LL;
typedef unsigned int uint;
const int _=1e2;
const int maxs=1e5+_;

namespace MIN_25
{
    int N,sqN;
    uint quick_pow(uint A,int p)
    {
        uint ret=1;
        while(p!=0)
        {
            if(p%2==1)ret=ret*A;
            A=A*A;p/=2;
        }
        return ret;
    }
    int pr,prime[maxs];uint kp[maxs]; bool v[maxs];
    void getprime(int K)
    {
        for(int i=2;i<=sqN;i++)
        {
            if(v[i]==false)prime[++pr]=i,kp[pr]=quick_pow(i,K);
            for(int j=1;j<=pr&&i*prime[j]<=sqN;j++)
            {
                v[i*prime[j]]=true;
                if(i%prime[j]==0)break;
            }
        }
    }
    
    int ulen,u[2*maxs],id1[maxs],id2[maxs];
    int ID(int x){return (x<=sqN)?id1[x]:id2[N/x];}
    uint g[2*maxs];//质数个数前缀和 
    void getg()
    {
        ulen=0;
        for(int i=1;i<=N;i=N/(N/i)+1)
        {
            u[++ulen]=N/i;
            if(u[ulen]<=sqN)id1[u[ulen]]=ulen;
            else id2[N/u[ulen]]=ulen;
            g[ulen]=u[ulen]-1;
        }
        
        for(int j=1;j<=pr;j++)
            for(int i=1;i<=ulen&&(LL)prime[j]*prime[j]<=u[i];i++)
                g[i]-=g[ID(u[i]/prime[j])]-(j-1);
    }
    
    uint S(int n,int j)
    {
        if(n<=1||prime[j]>n)return 0;
        
        uint sum=0;
        for(int i=j;(LL)prime[i]*prime[i]<=n&&i<=pr;i++)
        {
            uint p=prime[i];
            for(int q=1;(LL)p*prime[i]<=n;q++)
            {
                sum+=S(n/p,i+1);
                sum+=kp[i]*(g[ID(n/p)]-(i-1));
                p*=prime[i];
            }
        }
        if(j==1)sum+=g[ID(n)];
        return sum;
    }
}

namespace DJ
{
    int N,sqN;
    uint mp1[maxs],mp2[maxs];
    uint MP(uint x){return (x<=sqN)?mp1[x]:mp2[N/x];}
    void upd(int x,uint d)
    {
        if(x<=sqN)mp1[x]=d;
        else mp2[N/x]=d;
    }
    uint S(int n)
    {
        if(MP(n)!=0)return MP(n);
        uint ret=MIN_25::S(n,1); int last;
        for(int i=2;i<=n;i=last+1)
        {
            last=n/(n/i);
            ret-=(last-i+1)*S(n/last);
        }
        upd(n,ret);
        return ret;
    }
}

int main()
{
    int n,K;
    scanf("%d%d",&n,&K);
    
    DJ::N=MIN_25::N=n;
    DJ::sqN=MIN_25::sqN=int(sqrt(double(n+1)));
    MIN_25::getprime(K);
    MIN_25::getg();
    
    uint ans=0,ps=0,ns; int last;
    for(int i=2;i<=n;i=last+1)
    {
        last=n/(n/i);
        ns=DJ::S(last);
        ans+=(uint)(n/i)*(n/i)*(ns-ps);
        ps=ns;
    }
    printf("%u\n",ans);
    
    return 0;
}