Luogu P3768 简单的数学题

题目链接

https://www.luogu.org/problemnew/show/P3768

题解

反演一发得到
$$\sum _{T=1}^n(\sum _{i=1}^{\lfloor n/T\rfloor }i{)}^2{T}^2\phi (T)$$
设
$$f(T)=T^2\phi (T)$$
显然$f$可以杜教筛，因此答案可以整除分块求。

代码

#include <map>
#include <cstdio>
#include <algorithm>

template<typename T>
T read()
{
  T x=0;
  int f=1;
  char ch=getchar();
  while((ch<'0')||(ch>'9'))
    {
      if(ch=='-')
        {
          f=-f;
        }
      ch=getchar();
    }
  while((ch>='0')&&(ch<='9'))
    {
      x=x*10+ch-'0';
      ch=getchar();
    }
  return x*f;
}

const int maxn=4641588;

int mod,inv4,inv6;

int quickpow(int a,int b)
{
  int res=1;
  while(b)
    {
      if(b&1)
        {
          res=1ll*res*a%mod;
        }
      a=1ll*a*a%mod;
      b>>=1;
    }
  return res;
}

int p[maxn+10],prime[maxn+10],cnt,f[maxn+10];

int getprime()
{
  p[1]=f[1]=1;
  for(int i=2; i<=maxn; ++i)
    {
      if(!p[i])
        {
          prime[++cnt]=i;
          f[i]=i-1;
        }
      for(int j=1; (j<=cnt)&&(i*prime[j]<=maxn); ++j)
        {
          int x=i*prime[j];
          p[x]=1;
          if(i%prime[j]==0)
            {
              f[x]=f[i]*prime[j];
              break;
            }
          f[x]=f[i]*(prime[j]-1);
        }
    }
  for(int i=1; i<=maxn; ++i)
    {
      f[i]=(f[i-1]+1ll*i*i%mod*f[i])%mod;
    }
  return 0;
}

std::map<long long,int> mp;

int tri(long long n)
{
  int i=n%mod,ans=1ll*i*(i+1)%mod;
  return 1ll*ans*ans%mod*inv4%mod;
}

int dou(long long n)
{
  int i=n%mod;
  return 1ll*i*(i+1)%mod*(2*i+1)%mod*inv6%mod;
}

int getsum(long long n)
{
  if(n<=maxn)
    {
      return f[n];
    }
  if(mp.count(n))
    {
      return mp[n];
    }
  int ans=tri(n);
  for(long long l=2,r; l<=n; l=r+1)
    {
      r=n/(n/l);
      ans-=1ll*(dou(r)-dou(l-1)+mod)*getsum(n/l)%mod;
      if(ans<0)
        {
          ans+=mod;
        }
    }
  return mp[n]=ans;
}

long long n;

int main()
{
  mod=read<int>();
  n=read<long long>();
  getprime();
  inv4=quickpow(4,mod-2);
  inv6=quickpow(6,mod-2);
  int ans=0;
  for(long long l=1,r; l<=n; l=r+1)
    {
      r=n/(n/l);
      ans+=1ll*(getsum(r)-getsum(l-1)+mod)*tri(n/l)%mod;
      if(ans>=mod)
        {
          ans-=mod;
        }
    }
  printf("%d\n",ans);
  return 0;
}