[杜教篩 約數和字首和] 51Nod 1220 約數之和
阿新 • • 發佈:2019-02-07
套用公式後反演 然後杜教篩求和
比較有意思的是其間我算了兩個值
後來發現這兩個值竟然是相等的 都可以由約數和字首和推導過來
#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; }