加权约数和
Description
去年的 tangjz 非常喜欢做数论题,但是一年以后的 tangjz 却不那么会做了。
在整理以前的试题时,他发现了这样一道题目:“求 ∑σ(i) ,其中 1≤i≤N , σ(i) 表示 i 的约数之和。”
现在他长大了,题目也变难了,所以麻烦你来帮他解决一道数论题吧。
他需要你求如下表达式的值:
∑Ni=1∑Nj=1max(i,j)⋅σ(i⋅j)
其中 max(i,j) 表示 i 和 j 里的最大值, σ(i⋅j) 表示 i⋅j 的约数之和。
例如当 N=2 的时候,由 σ(1)=1,σ(2)=1+2=3,σ(4)=1+2+4=7 可知,答案应为 1⋅σ(1⋅1)+2⋅σ(1⋅2)+2⋅σ(2⋅1)+2⋅σ(2⋅2)=27 。
他发现答案有点大,所以你只需要告诉他答案模 1000000007 的值即可。
Input
每个测试点含有多组测试数据。
第一行是一个正整数T,表示接下来有T组测试数据。(1≤T≤50000)
接下来的T行,每组测试数据占一行。
每行有一个正整数N,含义如描述所示。(1≤N≤1000000)
Output
共有T行。对于每组测试数据,输出一行信息”Case #x: y”。
其中x表示对应的是第几组测试数据,y表示相应的答案模1000000007的值。
Input示例
5
1
2
3
4
5
Output示例
Case #1: 1
Case #2: 27
Case #3: 162
Case #4: 686
Case #5: 1741
看杜教筛结果被文末习题骗了过来……
推式子共用了10min多,实现方面调了2h+…….
这真的是一道推式子的题吗……
思路:
首先,考虑如何处理max的问题。
考虑想象一个矩阵,那么可以发现其max值沿对角线对称。
那么很自然的使用两倍的上三角矩阵减去对角线来代替原式:
可以发现对角线部分可以线性筛预处理,那么只考虑前半部分,即
考虑拆开
σ(i∗j)
:
反演:
将 d 提前:
将后面的两个求和符号化简:
至此完成这一部分的化简!
于是带回原式:
将
d
提前:
带回总式子:
考虑 n 的范围不大,于是线性筛出如下内容:
d[n]=σ(n)
sumd[n]=∑i=1nσ(i)
sum[n]=∑i=1niσ(i)∑j=1iσ(j)
(代码中未出现)
ds[n]=σ(n2)
sumds[n]=∑i=1ni∗σ(i2)
此处为全文精华部分,由于篇幅原因各种实现详见代码
于是现在就可以整除分块做到单词询问
O(n√)
了——
然而这并没有什么卯月。询问数巨大,因此依旧会T。
考虑 O(nln(n)) 预处理出 ans′ 的答案数组 f[n] ,每次 O(1) 回答询问,这样就能过了~
这真的不是在考线性筛的各种操作吗……
#include<iostream>
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<algorithm>
using namespace std;
typedef long long ll;
const int N=1e6+7;
const ll md=1e9+7;
ll n;
ll mu[N],pri[N>>1],ptop;
ll sumd[N],sumds[N];
ll mi1[N],mi2[N],d[N];
ll ds[N],mi1s[N],mi2s[N],f[N];
bool npri[N];
inline void init()
{
mu[1]=ds[1]=d[1]=1;
sumds[1]=sumd[1]=1;
for(ll i=2;i<N;i++)
{
if(!npri[i])
{
pri[++ptop]=i;
mu[i]=-1;
d[i]=i+1;
mi1[i]=i+1;
mi2[i]=i;
ds[i]=(ll)i*i+i+1ll;
mi1s[i]=(ll)i*i+i+1ll;
mi2s[i]=(ll)i*i;
}
for(int j=1;j<=ptop && i*pri[j]<N;j++)
{
npri[i*pri[j]]=1;
if(i%pri[j])
{
mu[i*pri[j]]=-mu[i];
d[i*pri[j]]=d[i]*d[pri[j]];
mi1[i*pri[j]]=1+pri[j];
mi2[i*pri[j]]=pri[j];
ds[i*pri[j]]=ds[i]*ds[pri[j]];
mi1s[i*pri[j]]=(ll)pri[j]*pri[j]+(ll)pri[j]+1ll;
mi2s[i*pri[j]]=(ll)pri[j]*pri[j];
}
else
{
mu[i*pri[j]]=0;
mi2[i*pri[j]]=mi2[i]*pri[j];
mi1[i*pri[j]]=mi1[i]+mi2[i*pri[j]];
d[i*pri[j]]=d[i]/mi1[i]*mi1[i*pri[j]];
mi2s[i*pri[j]]=mi2s[i]*pri[j]*pri[j];
mi1s[i*pri[j]]=mi1s[i]+mi2s[i]*((ll)pri[j]*pri[j]+pri[j]);
ds[i*pri[j]]=ds[i]/mi1s[i]*mi1s[i*pri[j]];
break;
}
}
sumd[i]=(sumd[i-1]+d[i])%md;
sumds[i]=(sumds[i-1]+ds[i]*i%md)%md;
}
for(int i=1;i<N;i++)
if(mu[i])
{
for(int j=i;j<N;j+=i)
{
if(mu[i]==1)
f[j]+=(ll)i*j%md*sumd[j/i]%md*d[j/i]%md;
else
f[j]+=md-(ll)i*j%md*sumd[j/i]%md*d[j/i]%md;
if(f[j]>=md)
f[j]-=md;
}
}
for(int i=1;i<N;i++)
{
f[i]+=f[i-1];
if(f[i]>=md)
f[i]-=md;
}
}
int main()
{
init();
int t;
scanf("%d",&t);
for(int tim=1;tim<=t;tim++)
{
scanf("%lld",&n);
printf("Case #%d: %lld\n",tim,((2ll*f[n]%md-sumds[n])%md+md)%md);
}
return 0;
}