[杜教筛约数和前缀和] 51Nod 1220 约数之和

最新推荐文章于 2020-09-17 10:50:44 发布

里阿奴摩西

最新推荐文章于 2020-09-17 10:50:44 发布

阅读量2.4k

点赞数

CC 4.0 BY-SA版权

分类专栏：莫比乌斯反演&杜教筛

本文链接：https://blog.youkuaiyun.com/u014609452/article/details/54747741

莫比乌斯反演&杜教筛专栏收录该内容

38 篇文章

订阅专栏

本文介绍了一种利用杜教筛优化数论求和问题的方法，并通过具体实现展示了如何进行高效的数值计算。文中涉及了多项式运算、约数求和及前缀和等技巧。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

吉丽博客传送门：http://jiruyi910387714.is-programmer.com/posts/195270.html

套用公式后反演然后杜教筛求和

比较有意思的是其间我算了两个值

后来发现这两个值竟然是相等的都可以由约数和前缀和推导过来

#include<cstdio>
#include<cstdlib>
#include<algorithm>
#include <tr1/unordered_map>
typedef long long ll;
using namespace std;
using namespace std::tr1;

const int P=1000000007;
const int inv2=(P+1)/2;
const int inv3=(P+1)/3;
const int maxn=1000000;

int prime[2000000],num;
int miu[maxn+5];
ll sumd[maxn+5],t1[maxn+5],t2[maxn+5]; int d[maxn+5];

inline void Pre(){
  miu[1]=1; sumd[1]=1;
  for (int i=2;i<=maxn;i++){
    if (!d[i]) d[i]=prime[++num]=i,sumd[i]=1+i,t1[i]=1+i,t2[i]=i,miu[i]=-1;
    for (int j=1;j<=num && (ll)i*prime[j]<=maxn;j++){
      d[i*prime[j]]=prime[j]; int k=i*prime[j];
      if (i%prime[j]==0){
	miu[k]=0;
	t2[k]=t2[i]*prime[j],t1[k]=t1[i]+t2[k],sumd[k]=sumd[i]/t1[i]*t1[k];
	break;
      }
      miu[i*prime[j]]=miu[i]*miu[prime[j]];
      sumd[k]=sumd[i]*sumd[prime[j]],t1[k]=1+prime[j],t2[k]=prime[j];
    }
  }
  for (int i=1;i<=maxn;i++)
    miu[i]=((ll)i*((P+miu[i])%P)%P+miu[i-1])%P,(sumd[i]+=sumd[i-1])%=P;
}

unordered_map<ll,int> S;

inline void add(int &x,int y){
  x+=y; if (x>=P) x-=P;
}

inline int sum1(ll n){ return (n+1)%P*(n%P)%P*inv2%P;}
inline int sum1(ll l,ll r){ return (l+r)%P*((r-l+1)%P)%P*inv2%P; }

inline int Sum(ll n){
  if (n<=maxn) return miu[n];
  if (S.find(n)!=S.end()) return S[n];
  int tem=1; ll l,r;
  for (l=2;l*l<=n;l++) add(tem,P-l*Sum(n/l)%P);
  for (ll t=n/l;l<=n;l=r+1,t--)
    r=n/t,add(tem,P-(ll)sum1(l,r)*Sum(t)%P);
  return S[n]=tem;
}
/*
inline int F(ll n){
  int t1=0,t2=0; ll l,r;
  for (l=1;l*l<=n;l++) add(t1,l%P*(n/l)%P),add(t2,sum1(n/l)%P);
  for (ll t=n/l;l<=n;l=r+1,t--)
    r=n/t,add(t1,(ll)sum1(l,r)*(n/l)%P),add(t2,(r-l+1)%P*sum1(n/l)%P);
  return (ll)t1*t2%P;
}
*/

inline int F(ll n){
  if (n<=maxn) return (ll)sumd[n]*sumd[n]%P;
  int t1=0; ll l,r;
  for (l=1;l*l<=n;l++) add(t1,l%P*(n/l)%P);
  for (ll t=n/l;l<=n;l=r+1,t--)
    r=n/t,add(t1,(ll)sum1(l,r)*(n/l)%P);
  return (ll)t1*t1%P;
}

int main(){
  ll n,l,r; int Ans=0;
  freopen("t.in","r",stdin);
  freopen("t.out","w",stdout);
  Pre();
  scanf("%lld",&n);
  for (l=1;l*l<=n;l++) add(Ans,(ll)(Sum(l)+P-Sum(l-1))%P*F(n/l)%P);
  for (ll t=n/l;l<=n;l=r+1,t--)
    r=n/t,add(Ans,(ll)(Sum(r)+P-Sum(l-1))%P*F(n/l)%P);
  printf("%d\n",Ans);
  return 0;
}