莫比乌斯函数前缀和
51nod - 1244
令
S(n)=∑ni=1μ(i)
,求
S(n),1≤n≤1010
做法
线性筛出S(n)的前 n23 项,然后利用分块优化求解上式,总时间复杂度 O(n23)
代码
#include<cstdio>
#include<algorithm>
#include<cstring>
using namespace std;
typedef long long LL;
void Read(LL &x){
char c;
while(c=getchar(),c!=EOF)
if(c>='0'&&c<='9'){
x=c-'0';
while(c=getchar(),c>='0'&&c<='9')
x=x*10+c-'0';
ungetc(c,stdin);
return;
}
}
int p[5000010],pcnt,mu[5000000+10],sum[5000000+10];
bool f[5000010],vis[5000010];
LL g[5000000+10],l,r;
void prepare(){
int i,j;
mu[1]=sum[1]=1;
for(i=2;i<=5000000;i++){
if(!f[i]){
p[++pcnt]=i;
mu[i]=-1;
}
for(j=1;i*p[j]<=5000000;j++){
f[i*p[j]]=1;
if(i%p[j]==0){
mu[i*p[j]]=0;
break;
}
mu[i*p[j]]=-mu[i];
}
sum[i]=sum[i-1]+mu[i];
}
}
void read(){
Read(l),Read(r);
}
LL solve(LL n,LL now){
if(n<=5000000)
return sum[n];
if(vis[now])
return g[now];
vis[now]=1;
LL i,j,&ret=g[now]=1;
for(i=2;i<=n;i=j+1){
j=n/(n/i);
ret-=(j-i+1)*solve(n/i,now*i);
}
return ret;
}
int main()
{
prepare();
read();
LL t1=solve(r,1);
memset(vis,0,sizeof vis);
LL t2=solve(l-1,1);
printf("%lld\n",t1-t2);
}
这种利用狄利克雷卷积筛前缀和的思路被称为杜教筛。
欧拉函数前缀和
51nod - 1239
令
S(n)=∑ni=1φ(i)
,求
S(n),1≤n≤1010
做法
线性筛出S(n)的前 n23 项,然后利用分块优化求解上式,总时间复杂度 O(n23)
代码
#include<cstdio>
#define MOD 1000000007
#include<algorithm>
using namespace std;
typedef long long LL;
void Read(LL &x){
char c;
while(c=getchar(),c!=EOF)
if(c>='0'&&c<='9'){
x=c-'0';
while(c=getchar(),c>='0'&&c<='9')
x=x*10+c-'0';
ungetc(c,stdin);
return;
}
}
int p[5000010],pcnt,phi[5000010];
bool f[5000010],vis[5000010];
LL g[5000010],sum[5000010],n;
void prepare(){
int i,j;
sum[1]=phi[1]=1;
for(i=2;i<=5000010;i++){
if(!f[i]){
p[++pcnt]=i;
phi[i]=i-1;
}
for(j=1;i*p[j]<=5000010;j++){
f[i*p[j]]=1;
if(i%p[j]==0){
phi[i*p[j]]=phi[i]*p[j];
break;
}
phi[i*p[j]]=phi[i]*phi[p[j]];
}
sum[i]=(sum[i-1]+phi[i])%MOD;
}
}
LL solve(LL n,LL now){
if(n<=5000010)
return sum[n];
if(vis[now])
return g[now];
vis[now]=1;
LL i,j,&ret=g[now]=0;
for(i=2;i<=n;i=j+1){
j=n/(n/i);
(ret-=(j-i+1)*solve(n/i,now*i)%MOD)%=MOD;
}
ret+=(n%MOD*((n+1)%MOD)%MOD*500000004%MOD)%MOD;
return ret%MOD;
}
int main()
{
prepare();
Read(n);
printf("%lld\n",(solve(n,1)+MOD)%MOD);
}
BZOJ3944 - Sum
上面两题二合一
代码
#include<cstdio>
#include<algorithm>
#include<cstring>
using namespace std;
typedef long long LL;
template<class T>
void Read(T &x){
char c;
while(c=getchar(),c!=EOF)
if(c>='0'&&c<='9'){
x=c-'0';
while(c=getchar(),c>='0'&&c<='9')
x=x*10+c-'0';
ungetc(c,stdin);
return;
}
}
int p[1000010],pcnt,mu[1000000+10],sum[1000000+10],phi[1000000+10],T;
bool f[1000010],vis[1000010];
LL n,sum2[1000000+10];
typedef pair<LL,LL>pll;
pll g[1000000+10];
void prepare(){
int i,j;
mu[1]=sum[1]=1;
phi[1]=sum2[1]=1;
for(i=2;i<=1000000;i++){
if(!f[i]){
p[++pcnt]=i;
mu[i]=-1;
phi[i]=i-1;
}
for(j=1;i*p[j]<=1000000;j++){
f[i*p[j]]=1;
if(i%p[j]==0){
mu[i*p[j]]=0;
phi[p[j]*i]=p[j]*phi[i];
break;
}
mu[i*p[j]]=-mu[i];
phi[i*p[j]]=phi[i]*phi[p[j]];
}
sum[i]=sum[i-1]+mu[i];
sum2[i]=sum2[i-1]+phi[i];
}
}
void read(){
Read(n);
}
pll solve(LL n,LL now){
if(n<=1000000)
return pll(sum[n],sum2[n]);
if(vis[now])
return g[now];
vis[now]=1;
LL i,j;
pll &ret=g[now]=pll(1,(n+1)*n/2);
for(i=2;i<=n;i=j+1){
j=n/(n/i);
pll t=solve(n/i,now*i);
ret=pll(ret.first-(j-i+1)*t.first,ret.second-(j-i+1)*t.second);
}
return ret;
}
int main()
{
prepare();
Read(T);
while(T--){
memset(vis,0,sizeof vis);
read();
pll t=solve(n,1);
printf("%lld %lld\n",t.second,t.first);
}
}
HDU5608 - function
已知
N2−3N+2=∑d|Nf(d)
,求
∑Ni=1f(i) mod 109+7
多组数据,
T≤500,N≤109
,其中只有5组
N>106
做法
然而f(j)的前缀和并不好线性筛,所以时间复杂度为 O(n34)
有没有其他方法呢?
我们莫比乌斯反演一下。
我们就可以像求莫比乌斯函数的前缀和那样做,时间复杂度 O(n23)
代码
#include<cstdio>
#include<algorithm>
#include<cstring>
#define MAXL 1000000
#define MOD 1000000007
using namespace std;
template<class T>
void Read(T &x){
char c;
while(c=getchar(),c!=EOF)
if(c>='0'&&c<='9'){
x=c-'0';
while(c=getchar(),c>='0'&&c<='9')
x=x*10+c-'0';
ungetc(c,stdin);
return;
}
}
int n,p[MAXL+10],pcnt,mu[MAXL+10],sum[MAXL+10],g[MAXL+10],ans,T;
bool f[MAXL+10];
bool vis[MAXL+10];
typedef long long LL;
void prepare(){
int i,j;
mu[1]=sum[1]=1;
for(i=2;i<=MAXL;i++){
if(!f[i]){
p[++pcnt]=i;
mu[i]=-1;
}
for(j=1;i*p[j]<=MAXL;j++){
f[i*p[j]]=1;
if(i%p[j]==0){
mu[i*p[j]]=0;
break;
}
mu[i*p[j]]=-mu[i];
}
sum[i]=mu[i]+sum[i-1];
}
}
inline LL Sum(LL x){
return x*(x-1)%MOD*(x-2)%MOD*((MOD+1)/3)%MOD;
}
int solve(int n,int now){
if(n<=MAXL)
return sum[n];
if(vis[now])
return g[now];
vis[now]=1;
int i,j,&ret=g[now]=1;
for(i=2;i<=n;i=j+1){
j=n/(n/i);
ret=(ret-1ll*(j-i+1)*solve(n/i,now*i))%MOD;
}
return ret;
}
int main()
{
prepare();
Read(T);
while(T--){
Read(n);
if(n>=MAXL)
memset(vis,0,sizeof vis);
int i,j;
ans=0;
for(i=1;i<=n;i=j+1){
j=n/(n/i);
ans=(ans+((Sum(j)-Sum(i-1))%MOD*solve(n/i,i))%MOD)%MOD;
}
printf("%d\n",(ans+MOD)%MOD);
}
}
51nod1237 - 最大公约数之和 V3
G=0;
for(i=1;i<N;i++)
for(j=1;j<=N;j++)
{
G = (G + gcd(i,j)) % 1000000007;
}
求G
做法
令
A(n)=∑ni=1gcd(i,n),F(n)=∑ni=1A(i)
计算欧拉函数的前缀和即可,时间复杂度 O(n23)
代码
#include<cstdio>
#define MOD 1000000007
#include<algorithm>
using namespace std;
typedef long long LL;
void Read(LL &x){
char c;
while(c=getchar(),c!=EOF)
if(c>='0'&&c<='9'){
x=c-'0';
while(c=getchar(),c>='0'&&c<='9')
x=x*10+c-'0';
ungetc(c,stdin);
return;
}
}
int p[5000010],pcnt,phi[5000010];
bool f[5000010],vis[5000010];
LL g[5000010],sum[5000010],n,ans;
void prepare(){
int i,j;
sum[1]=phi[1]=1;
for(i=2;i<=5000010;i++){
if(!f[i]){
p[++pcnt]=i;
phi[i]=i-1;
}
for(j=1;i*p[j]<=5000010;j++){
f[i*p[j]]=1;
if(i%p[j]==0){
phi[i*p[j]]=phi[i]*p[j];
break;
}
phi[i*p[j]]=phi[i]*phi[p[j]];
}
sum[i]=(sum[i-1]+phi[i])%MOD;
}
}
LL solve(LL n,LL now){
if(n<=5000010)
return sum[n];
if(vis[now])
return g[now];
vis[now]=1;
LL i,j,&ret=g[now]=0;
for(i=2;i<=n;i=j+1){
j=n/(n/i);
(ret-=(j-i+1)*solve(n/i,now*i)%MOD)%=MOD;
}
ret+=(n%MOD*((n+1)%MOD)%MOD*500000004%MOD)%MOD;
return ret%MOD;
}
inline LL Sum(LL x){
return x%MOD*((x+1)%MOD)%MOD*500000004%MOD;
}
int main()
{
prepare();
Read(n);
LL i,j;
for(i=1;i<=n;i=j+1){
j=n/(n/i);
ans=(ans+(Sum(j)-Sum(i-1))*solve(n/i,i))%MOD;
}
printf("%lld\n",((((ans+MOD)%MOD)*2-Sum(n))%MOD+MOD)%MOD);
}
51nod1227 - 平均最小公倍数
Lcm(a,b)表示a和b的最小公倍数,A(n)表示Lcm(n,i)的平均数(1 <= i <= n),
例如:A(4) = (Lcm(1,4) + Lcm(2,4) + Lcm(3,4) + Lcm(4,4)) / 4 = (4 + 4 + 12 + 4) / 4 = 6。
F(a, b) = A(a) + A(a + 1) + …… A(b)。(F(a,b) = ∑A(k), a <= k <= b)
例如:F(2, 4) = A(2) + A(3) + A(4) = 2 + 4 + 6 = 12。
给出a,b,计算F(a, b),由于结果可能很大,输出F(a, b) % 1000000007的结果即可。
做法
令
A(n)=∑ni=1lcm(i,n)n,F(n)=∑ni=1A(n)
令 G(n)=2×F(n)−n
怎样快速求 ∑ni=1i×φ(i) 呢?
我们构造一下,发现 ∑d|nd×φ(d)×nd=n2
线性筛 ∑ni=1i×φ(i) 的qian n23 项,总时间复杂度 O(n23)
代码
#include<cstdio>
#include<algorithm>
#include<cstring>
using namespace std;
#define MOD 1000000007
#define MAXL 1000000
typedef long long LL;
void Read(LL &x){
char c;
while(c=getchar(),c!=EOF)
if(c>='0'&&c<='9'){
x=c-'0';
while(c=getchar(),c>='0'&&c<='9')
x=x*10+c-'0';
ungetc(c,stdin);
return;
}
}
LL l,r;
int p[MAXL+10],sum[MAXL+10],pcnt,phi[MAXL+10],ans;
LL g[MAXL+10];
bool f[MAXL+10],vis[MAXL+10];
void prepare(){
int i,j;
phi[1]=sum[1]=1;
for(i=2;i<=MAXL;i++){
if(!f[i]){
p[++pcnt]=i;
phi[i]=i-1;
}
for(j=1;i*p[j]<=MAXL;j++){
f[i*p[j]]=1;
if(i%p[j]==0){
phi[i*p[j]]=phi[i]*p[j];
break;
}
phi[i*p[j]]=phi[i]*phi[p[j]];
}
sum[i]=(sum[i-1]+1ll*i*phi[i])%MOD;
}
}
LL Sum2(LL n){
n%=MOD;
return n%MOD*((n+1)%MOD)%MOD*(2*n+1)%MOD*((MOD+1)/6)%MOD;
}
LL Sum(LL n){
n%=MOD;
return n*(n+1)%MOD*((MOD+1)/2)%MOD;
}
LL Sum3(LL n){
return 1ll*Sum(n)*Sum(n)%MOD;
}
LL solve(LL n,LL now){
if(n<=MAXL)
return sum[n];
if(vis[now])
return g[now];
vis[now]=1;
LL i,j,&ret=g[now]=Sum2(n);
for(i=2;i<=n;i=j+1){
j=n/(n/i);
ret=(ret-(Sum(j)-Sum(i-1))*solve(n/i,now*i))%MOD;
}
return ret;
}
LL cal(LL n){
memset(vis,0,sizeof vis);
LL i,j,ret(0);
for(i=1;i<=n;i=j+1){
j=n/(n/i);
ret=(ret+(j-i+1)*solve(n/i,i))%MOD;
}
return (ret+n)*((MOD+1)>>1)%MOD;
}
int main()
{
prepare();
Read(l),Read(r);
printf("%lld\n",((cal(r)-cal(l-1))%MOD+MOD)%MOD);
}
51nod1238 - 最小公倍数之和 V3
G=0;
for(i=1;i<N;i++)
for(j=1;j<=N;j++)
{
G = (G + lcm(i,j)) % 1000000007;
}
求G
做法
令 A(n)=∑ni=1lcm(i,n),F(n)=∑ni=1A(i)
令 G(n)=2×F(n)−∑ni=1i
怎样快速求 ∑ni=1i2×φ(i) 呢?
类似上一题,我们构造一下
线性筛 ∑ni=1i2×φ(i) 的前 n23 项,总的时间复杂度 O(n23) ,答案就是 G(n)
代码
#include<cstdio>
#include<algorithm>
using namespace std;
#define MOD 1000000007
#define MAXL 10000000
typedef long long LL;
void Read(LL &x){
char c;
while(c=getchar(),c!=EOF)
if(c>='0'&&c<='9'){
x=c-'0';
while(c=getchar(),c>='0'&&c<='9')
x=x*10+c-'0';
ungetc(c,stdin);
return;
}
}
LL n;
int p[MAXL+10],sum[MAXL+10],pcnt,phi[MAXL+10],ans;
LL g[MAXL+10];
bool f[MAXL+10],vis[MAXL+10];
void prepare(){
int i,j;
phi[1]=sum[1]=1;
for(i=2;i<=MAXL;i++){
if(!f[i]){
p[++pcnt]=i;
phi[i]=i-1;
}
for(j=1;i*p[j]<=MAXL;j++){
f[i*p[j]]=1;
if(i%p[j]==0){
phi[i*p[j]]=phi[i]*p[j];
break;
}
phi[i*p[j]]=phi[i]*phi[p[j]];
}
sum[i]=(sum[i-1]+1ll*i*i%MOD*phi[i])%MOD;
}
}
LL Sum2(LL n){
n%=MOD;
return n%MOD*((n+1)%MOD)%MOD*(2*n+1)%MOD*((MOD+1)/6)%MOD;
}
LL Sum(LL n){
n%=MOD;
return n*(n+1)%MOD*((MOD+1)/2)%MOD;
}
LL Sum3(LL n){
return 1ll*Sum(n)*Sum(n)%MOD;
}
LL solve(LL n,LL now){
if(n<=MAXL)
return sum[n];
if(vis[now])
return g[now];
vis[now]=1;
LL i,j,&ret=g[now]=Sum3(n);
for(i=2;i<=n;i=j+1){
j=n/(n/i);
ret=(ret-(Sum2(j)-Sum2(i-1))*solve(n/i,now*i))%MOD;
}
return ret;
}
int main()
{
prepare();
Read(n);
LL i,j;
for(i=1;i<=n;i=j+1){
j=n/(n/i);
ans=(ans+(Sum(j)-Sum(i-1))*solve(n/i,i))%MOD;
}
printf("%d\n",(ans+MOD)%MOD);
}
SPOJ DIVCNT2 - Counting Divisors(square)
求 S2(N) .
做法
ω(n) 表示n的不同质因子的数量,可以看做在 d2 中减少 d 的一个质因数集合。
时间复杂度据说 O(n23)
代码
#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
using namespace std;
template<class T>
void Read(T &x){
char c;
while(c=getchar(),c!=EOF)
if(c>='0'&&c<='9'){
x=c-'0';
while(c=getchar(),c>='0'&&c<='9')
x=x*10+c-'0';
ungetc(c,stdin);
return;
}
}
typedef long long LL;
int p[10000010],mu[100000010],cf[100000010],L=1e7,n,pcnt,T,smu[100000010];
bool f[100000010];
char cs[100000010];
LL a[10000+10],mxa;
void read(){
Read(n);
}
void prepare(int n){
int i,j;
smu[1]=mu[1]=1;
cf[1]=1;
for(i=2;i<=n;i++){
if(!f[i]){
mu[i]=-1;
p[++pcnt]=i;
cs[i]=1,cf[i]=2;
}
for(j=1;i*p[j]<=n;j++){
f[i*p[j]]=1;
if(i%p[j]==0){
mu[i*p[j]]=0;
cs[i*p[j]]=cs[i]+1;
cf[i*p[j]]=cf[i]/(cs[i]+1)*(cs[i]+2);
break;
}
mu[i*p[j]]=-mu[i];
cs[i*p[j]]=1;
cf[i*p[j]]=cf[i]<<1;
}
cf[i]+=cf[i-1];
smu[i]=abs(mu[i])+smu[i-1];
}
}
LL cal(LL x){
if(x<=L)
return cf[x];
LL ret=0,i,j;
for(i=1;i<=x;i=j+1){
j=x/(x/i);
ret+=(j-i+1)*(x/i);
}
return ret;
}
LL Sum(LL x){
if(x<=L)
return smu[x];
LL ret(0);
for(LL i=1;i*i<=x;i++)
if(mu[i])
ret+=mu[i]*(x/(i*i));
return ret;
}
LL solve(LL n){
int t=sqrt(n+0.5);
LL i,j;
LL ret(0),pre=smu[t];
for(i=1;i<=t;i++)
if(mu[i])
ret+=cal(n/i);
for(i=t+1;i<=n;i=j+1){
j=n/(n/i);
LL tt=Sum(j);
ret+=(tt-pre)*cal(n/i);
pre=tt;
}
return ret;
}
int main()
{
Read(T);
for(int i=1;i<=T;i++)
Read(a[i]),mxa=max(mxa,a[i]);
L=max(1.0*L,pow(mxa,2.0/3));
prepare(L);
for(int i=1;i<=T;i++)
printf("%lld\n",solve(a[i]));
}
51nod1222 - 最小公倍数计数
定义F(n)表示最小公倍数为n的二元组的数量。
即:如果存在两个数(二元组)X,Y(X <= Y),它们的最小公倍数为N,则F(n)的计数加1。
例如:F(6) = 5,因为[2,3] [1,6] [2,6] [3,6] [6,6]的最小公倍数等于6。
给出一个区间[a,b],求最小公倍数在这个区间的不同二元组的数量。
例如:a = 4,b = 6。符合条件的二元组包括:
[1,4] [2,4] [4,4] [1,5] [5,5] [2,3] [1,6] [2,6] [3,6] [6,6],共10组不同的组合。
Input
输入数据包括2个数:a, b,中间用空格分隔(1 <= a <= b <= 10^11)。
Output
输出最小公倍数在这个区间的不同二元组的数量。
做法
令
S(n)=∑ni=1F(i)
,
ans[a,b]=S(b)−S(a−1)
莫比乌斯反演一下,
我们发现,其实 ∑n√r=1μ(r)∑nd=1∑ni=1∑nj=1[abd<=nr2] 还是很好算的,我们枚d,i就可以了,然后分别统计一个相同,两个相同,三个相同的数量,然后就可以轻易算出答案了,据说复杂度 O(n23)
代码
#include<cstdio>
#include<algorithm>
#include<cmath>
#define MAXL 1000000
using namespace std;
typedef long long LL;
void Read(LL &x){
char c;
while(c=getchar(),c!=EOF)
if(c>='0'&&c<='9'){
x=c-'0';
while(c=getchar(),c>='0'&&c<='9')
x=x*10+c-'0';
ungetc(c,stdin);
return;
}
}
LL l,r;
int p[MAXL+10],pcnt,mu[MAXL+10],sum[MAXL+10];
bool f[MAXL+10];
void read(){
Read(l),Read(r);
}
void prepare(){
mu[1]=sum[1]=1;
int i,j;
for(i=2;i<=MAXL;i++){
if(!f[i]){
p[++pcnt]=i;
mu[i]=-1;
}
for(j=1;i*p[j]<=MAXL;j++){
f[i*p[j]]=1;
if(i%p[j]==0){
mu[i*p[j]]=0;
break;
}
mu[i*p[j]]=-mu[i];
}
sum[i]=sum[i-1]+mu[i];
}
}
LL cal(LL n){
LL i,cnt1=0,cnt2=0,cnt3=0,j,k;
for(i=1;i*i*i<=n;i++){
cnt3++;
for(j=i+1;i*j*j<=n;j++){
cnt2++;
// for(k=j+1;i*j*k<=n;k++)
// cnt1++;
cnt1+=max(0ll,n/(i*j)-j);
}
}
for(i=1;i*i*i<=n;i++)
cnt2+=max(n/(i*i)-i,0ll);
return cnt1*6+cnt2*3+cnt3;
}
LL solve(LL n){
LL i,j,ans=0;
for(i=1;i*i<=n;i=j+1){
j=n/(n/(i*i));
j=sqrt(j+0.5);
ans+=(sum[j]-sum[i-1])*cal(n/(i*i));
}
return (ans+n)/2;
}
int main()
{
prepare();
read();
printf("%lld\n",solve(r)-solve(l-1));
}
BZOJ4176 - Lucas的数论
去年的Lucas非常喜欢数论题,但是一年以后的Lucas却不那么喜欢了。
在整理以前的试题时,发现了这样一道题目“求Sigma(f(i)),其中1<=i<=N”,其中 表示i的约数个数。他现在长大了,题目也变难了。
求如下表达式的值:
其中 表示ij的约数个数。
他发现答案有点大,只需要输出模1000000007的值。
做法
d|ij
等价于
dgcd(i,d)|j
据说时间复杂度 O(n34)
代码
#include<cstdio>
#include<algorithm>
#define MOD 1000000007
#define MAXL 1000000
using namespace std;
typedef long long LL;
LL n;
int p[MAXL+10],pcnt,mu[MAXL+10],sum[MAXL+10];
bool f[MAXL+10],vis[MAXL+10];
LL g[MAXL+10],ans;
void prepare(){
mu[1]=sum[1]=1;
int i,j;
for(i=2;i<=MAXL;i++){
if(!f[i]){
p[++pcnt]=i;
mu[i]=-1;
}
for(j=1;i*p[j]<=MAXL;j++){
f[i*p[j]]=1;
if(i%p[j]==0){
mu[i*p[j]]=0;
break;
}
mu[i*p[j]]=-mu[i];
}
sum[i]=sum[i-1]+mu[i];
}
}
void Read(LL &x){
char c;
while(c=getchar(),c!=EOF)
if(c>='0'&&c<='9'){
x=c-'0';
while(c=getchar(),c>='0'&&c<='9')
x=x*10+c-'0';
ungetc(c,stdin);
return;
}
}
inline LL cal(LL n){
LL i,j,ret(0);
for(i=1;i<=n;i=j+1){
j=n/(n/i);
ret=(ret+(j-i+1)*(n/i))%MOD;
}
return ret;
}
LL solve(LL n,LL now){
if(n<=MAXL)
return sum[n];
if(vis[now])
return g[now];
vis[now]=1;
LL i,j,&ret=g[now]=1;
for(i=2;i<=n;i=j+1){
j=n/(n/i);
ret=(ret-(j-i+1)*solve(n/i,now*i))%MOD;
}
return ret;
}
int main()
{
prepare();
Read(n);
LL i,j;
for(i=1;i<=n;i=j+1){
j=n/(n/i);
LL t=cal(n/i);
ans=(ans+(solve(j,n/i)-solve(i-1,n/((i-1)?(i-1):1)))*t%MOD*t)%MOD;
}
printf("%lld\n",(ans+MOD)%MOD);
}
51nod1220 - 约数之和
d(k)表示k的所有约数的和。d(6) = 1 + 2 + 3 + 6 = 10。
定义S(N) = ∑1<=i<=N ∑1<=j<=N d(i*j)。
例如:S(3) = d(1) + d(2) + d(3) + d(2) + d(4) + d(6) + d(3) + d(6) + d(9) = 59,S(1000) = 563576517282。
给出正整数N,求S(N),由于结果可能会很大,输出Mod 1000000007(10^9 + 7)的结果。
做法
和上一题类似,做法还是看http://www.51nod.com/question/index.html#!questionId=1027
代码
#include<cstdio>
#include<algorithm>
#define MOD 1000000007
#define MAXL 1000000
using namespace std;
typedef long long LL;
LL n;
int p[MAXL+10],pcnt,mu[MAXL+10],sum[MAXL+10];
bool f[MAXL+10],vis[MAXL+10];
LL g[MAXL+10],ans;
void prepare(){
mu[1]=sum[1]=1;
int i,j;
for(i=2;i<=MAXL;i++){
if(!f[i]){
p[++pcnt]=i;
mu[i]=-1;
}
for(j=1;i*p[j]<=MAXL;j++){
f[i*p[j]]=1;
if(i%p[j]==0){
mu[i*p[j]]=0;
break;
}
mu[i*p[j]]=-mu[i];
}
sum[i]=(sum[i-1]+mu[i]*i)%MOD;
}
}
void Read(LL &x){
char c;
while(c=getchar(),c!=EOF)
if(c>='0'&&c<='9'){
x=c-'0';
while(c=getchar(),c>='0'&&c<='9')
x=x*10+c-'0';
ungetc(c,stdin);
return;
}
}
inline LL Sum(LL n){
return n*(n+1)/2%MOD;
}
inline LL cal(LL n){
LL i,j,ret(0);
for(i=1;i<=n;i=j+1){
j=n/(n/i);
ret=(ret+(Sum(j)-Sum(i-1))*(n/i))%MOD;
}
return ret;
}
inline LL cal2(LL n){
LL i,j,ret(0);
for(i=1;i<=n;i=j+1){
j=n/(n/i);
ret=(ret+(j-i+1)*Sum(n/i))%MOD;
}
return ret;
}
LL solve(LL n,LL now){
if(n<=MAXL)
return sum[n];
if(vis[now])
return g[now];
vis[now]=1;
LL i,j,&ret=g[now]=1;
for(i=2;i<=n;i=j+1){
j=n/(n/i);
ret=(ret-(Sum(j)-Sum(i-1))*solve(n/i,now*i))%MOD;
}
return ret;
}
int main()
{
prepare();
Read(n);
LL i,j;
for(i=1;i<=n;i=j+1){
j=n/(n/i);
ans=(ans+(solve(j,n/i)-solve(i-1,n/((i-1)?(i-1):1)))*cal2(n/i)%MOD*cal(n/i))%MOD;
}
printf("%lld\n",(ans+MOD)%MOD);
}
BZOJ3512 - DZY Loves Math IV
给定n,m,求 模10^9+7的值。
做法
令
S(n,m)=∑mi=1φ(i×n),ans=∑ni=1S(i,m)
然后记忆化搜索+杜教筛求解即可。
代码
#include<cstdio>
#include<algorithm>
#include<map>
#define MAXL 1000000
#define MOD 1000000007
using namespace std;
using namespace __gnu_cxx;
typedef pair<int,int>pii;
map<pii,int>h;
void Read(int &x){
char c;
while(c=getchar(),c!=EOF)
if(c>='0'&&c<='9'){
x=c-'0';
while(c=getchar(),c>='0'&&c<='9')
x=x*10+c-'0';
ungetc(c,stdin);
return;
}
}
#define MAXN 100000
typedef long long LL;
int n,m,ans[MAXN+10],sf[MAXL+10],p[MAXL+10],pcnt,sum[MAXL+10],phi[MAXL+10];
bool f[MAXL+10],vis[MAXL+10];
LL g[MAXL+10];
void read(){
Read(n),Read(m);
}
void prepare(){
int i,j;
phi[1]=sum[1]=1;
for(i=2;i<=MAXL;i++){
if(!f[i]){
sf[i]=i;
p[++pcnt]=i;
phi[i]=i-1;
}
for(j=1;i*p[j]<=MAXL;j++){
f[i*p[j]]=1;
sf[i*p[j]]=sf[i];
if(i%p[j]==0){
phi[i*p[j]]=phi[i]*p[j];
break;
}
phi[i*p[j]]=phi[i]*phi[p[j]];
}
sum[i]=(sum[i-1]+phi[i])%MOD;
}
}
LL Get_fac(int n){
int i,ret(1);
for(i=1;p[i]*p[i]<=n;i++){
if(n%p[i]==0){
n/=p[i];
while(n%p[i]==0)
ret*=p[i],n/=p[i];
}
}
return ret;
}
LL Sum(LL n){
return n*(n+1)/2%MOD;
}
LL cal(LL n){
if(n<=MAXL)
return sum[n];
if(vis[m/n])
return g[m/n];
vis[m/n]=1;
LL i,j,&ret=g[m/n]=Sum(n);
for(i=2;i<=n;i=j+1){
j=n/(n/i);
ret=(ret-(j-i+1)*cal(n/i))%MOD;
}
return ret;
}
LL solve(LL n,LL m){
if(h.count(pii(n,m)))
return h[pii(n,m)];
if(!m)
return 0;
if(n==1)
return cal(m);
if(m==1)
return phi[n];
return h[pii(n,m)]=(phi[sf[n]]*solve(n/sf[n],m)+solve(n,m/sf[n]))%MOD;
}
LL solve(){
int i,t,ret=0;
for(i=1;i<=n;i++){
// printf("%d\n",i);
t=Get_fac(i);
if(t!=1)
ret=(ret+1ll*t*ans[i/t])%MOD;
else
ret=(ret+(ans[i]=solve(i,m)))%MOD;
}
return ret;
}
int main()
{
prepare();
read();
printf("%lld\n",(solve()+MOD)%MOD);
}