LOJ 572 「LibreOJ Round #11」Misaka Network 与求和——min_25筛

题目：https://loj.ac/problem/572

莫比乌斯反演得 ( ans=sumlimits_{D=1}^{n}leftlfloorfrac{n}{D} ight floor^2sumlimits_{d|D}f(d)^kmu (frac{D}{d}) )

计算 ( S(n)=sumlimits_{i=1}^{n}f×mu )

像杜教筛（https://blog.csdn.net/a1799342217/article/details/80328510）一样写一写式子，因为有 ( mu ) 所以把补的 g 函数设为 1 函数。

　　( sumlimits_{i=1}^{n}f×mu×1(i) )

　( = sumlimits_{i=1}^{n}sumlimits_{d|D}1(d)*(f×mu)(frac{i}{d}) )

　( = sumlimits_{d=1}^{n}1(d)sumlimits_{i=1}^{frac{n}{d}}f×mu (i) )

　( = sumlimits_{d=1}^{n}1(d)S(frac{n}{d}) )

用这个表示 ( S(n) )，则

　　( S(n)=sumlimits_{d=1}^{n}1(d)S(frac{n}{d}) - sumlimits_{d=2}^{n}1(d)S(frac{n}{d}) )

　　( S(n)=sumlimits_{i=1}^{n}f×mu×1(i)- sumlimits_{d=2}^{n}1(d)S(frac{n}{d}) )

　　( S(n)=sumlimits_{i=1}^{n}f(i)- sumlimits_{d=2}^{n}1(d)S(frac{n}{d}) )

用和 UOJ 188 一样的方法求 f 的前缀和即可。

 S 的可能角标一定是某个 ( leftlfloorfrac{n}{i} ight floor ) 。所以 S 可以预处理。预处理的时候别忘了记忆化。

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define int unsigned int
#define ll long long
using namespace std;
int pw(int x,int k)
{int ret=1;while(k){if(k&1)ret*=x;x*=x;k>>=1;}return ret;}
const int N=9e4+5;
int n,k,base,s[N],g[N],w[N],m,p[N],pk[N],cnt,ans[N];
bool vis[N];ll p2[N];
int Id(int x){return x<=base?m-x+1:n/x;}
int S(int n)
{
  int k=Id(n);if(ans[k]!=-1)return ans[k];
  int ret=s[k]+g[k];
  for(int i=2,j;i<=n;i=n/j+1)
    {j=n/i;ret-=(n/j-i+1)*S(j);}
  return ans[k]=ret;
}
void init(int n)
{
  base=sqrt(n);
  for(int i=2;i<=base;i++)
    {
      if(!vis[i])p[++cnt]=i,p2[cnt]=(ll)i*i,pk[cnt]=pw(i,k);
      for(int j=1,d;j<=cnt&&(d=i*p[j])<=base;j++)
    {vis[d]=1;if(i%p[j]==0)break;}
    }
  for(int i=1,j;i<=n;i=n/j+1)w[++m]=j=n/i;

  for(int i=1;i<=m;i++)g[i]=w[i]-1;
  for(int j=1,tp=0;j<=cnt;j++,tp++)
    for(int i=1;i<=m&&p2[j]<=w[i];i++)
      g[i]-=g[Id(w[i]/p[j])]-tp;

  int p0=1;
  for(int j=cnt;j;j--)
    {
      while(p0<=m&&p2[j]<=w[p0])p0++;
      for(int i=p0-1;i;i--)
    {
      int k=Id(w[i]/p[j]);
      s[i]+=s[k]+pk[j]*(g[k]-(j-1));
    }
    }
  memset(ans,-1,sizeof ans);
  S(n);
}
signed main()
{
  scanf("%u%u",&n,&k);init(n);int prn=0;
  for(int i=1,j,lst=0,nw;i<=n;i=n/j+1)
    {
      j=n/i;nw=ans[Id(n/j)];
      prn+=j*j*(nw-lst);lst=nw;
    }
  printf("%u
",prn);
  return 0;
}