CE玩家

人生最大的快乐,就是做别人说你做不到的事。

[莫比乌斯反演] 51nod1584. 加权约数和

杜教筛学傻了…一道反演题…我竟然差点推成杜教筛??

A(n)=i=1nσ(in)

这个东西推一下就变成

t|nμ(t)tf(nt)g(nt)

f(n)=d|nd,g(n)=i=1nf(n)

这个可O(nlnn)

那么答案就是

2i=1niA(i)i=1niσ(i2)

预处理后 O(1) 输出就好了

最后还被卡常啊……
最后能线性筛的东西全改成线性筛才过……

#include <cstdio>
#include <iostream>
#include <algorithm>
#include <ctime>

using namespace std;

typedef long long ll;

const int N=1000010,P=1e9+7;

int p[N],mu[N],f[N],g[N],h[N],ans[N],fac[N],w[N];

inline void Pre(const int n){
  mu[1]=1; f[1]=1;
  for(int i=2;i<=n;i++){
    if(!p[i])
      p[++*p]=i,mu[i]=-1,fac[i]=i,w[i]=(1LL*i*i+i+1)%P,f[i]=i+1;
    for(int j=1;j<=*p && 1LL*i*p[j]<=n;j++){
      p[i*p[j]]=1;
      if(i%p[j]){
    mu[i*p[j]]=-mu[i];
    fac[i*p[j]]=p[j];
    w[i*p[j]]=1LL*w[i]*w[p[j]]%P;
    f[i*p[j]]=1LL*f[i]*f[p[j]]%P;
      }
      else{
    fac[i*p[j]]=fac[i]*p[j];
    if(fac[i]==i){
      fac[i*p[j]]=fac[i]*p[j];
      w[i*p[j]]=(w[i]+1LL*fac[i]*fac[i]%P*p[j]+1LL*fac[i]*fac[i]%P*p[j]%P*p[j])%P;
      f[i*p[j]]=(f[i]+fac[i]*p[j])%P;
    }
    else{
      fac[i*p[j]]=fac[i]*p[j];
      w[i*p[j]]=1LL*w[i/fac[i]]*w[p[j]*fac[i]]%P;
      f[i*p[j]]=1LL*f[i/fac[i]]*f[p[j]*fac[i]]%P;
    }
    break;
      }
    }
  }
  for(int i=1;i<=n;i++) g[i]=(g[i-1]+f[i])%P;
  for(int i=1;i<=n;i++)
    for(int j=i;j<=n;j+=i)
      h[j]=(h[j]+1LL*mu[i]*i*f[j/i]%P*g[j/i])%P;
  for(int i=1;i<=n;i++)
    ans[i]=(ans[i-1]+1LL*i*h[i])%P;
  w[1]=1;
  for(int i=1;i<=n;i++) w[i]=(1LL*i*w[i]+w[i-1])%P;
  for(int i=1;i<=n;i++)
    ans[i]=(2LL*ans[i]-w[i])%P,ans[i]=(ans[i]+P)%P;
}

int main(){
  freopen("1.in","r",stdin);
  freopen("1.out","w",stdout);
  int t,n; scanf("%d",&t); Pre(1000000);
  for(int i=1;i<=t;i++)
    scanf("%d",&n),printf("Case #%d: %d\n",i,ans[n]);
  return 0;
}
阅读更多
个人分类: 莫比乌斯反演
想对作者说点什么? 我来说一句

没有更多推荐了,返回首页

不良信息举报

[莫比乌斯反演] 51nod1584. 加权约数和

最多只允许输入30个字

加入CSDN,享受更精准的内容推荐,与500万程序员共同成长!
关闭
关闭