狄利克雷卷积。
就和多项式的卷积一样。
多项式是
g
(
x
)
×
h
(
x
)
=
∑
a
+
b
=
x
a
n
d
a
,
b
∈
N
g
(
a
)
∗
h
(
b
)
g(x) \times h(x)= \sum_{a+b = x\ and \ a,b\in \N} g(a) * h(b)
g(x)×h(x)=a+b=x and a,b∈N∑g(a)∗h(b)
而狄利克雷卷积是
g
(
x
)
×
h
(
x
)
=
∑
a
∗
b
=
x
a
n
d
a
,
b
∈
N
+
g
(
a
)
∗
h
(
b
)
g(x) \times h(x) = \sum_{a*b = x \ and \ a,b\in \N_+ } g(a) * h(b)
g(x)×h(x)=a∗b=x and a,b∈N+∑g(a)∗h(b)
而莫比乌斯反演就是基于
I
×
μ
=
e
I \times \mu = e
I×μ=e
其中
I
(
x
)
=
1
,
e
(
x
)
=
[
x
=
=
1
]
I(x) = 1 , e(x) = [x == 1]
I(x)=1,e(x)=[x==1]
e
(
x
)
e(x)
e(x)是莫比乌斯的单位元。
而积性函数的前缀和的快速算法1则是基于这样一个事实:
∵
∑
a
∗
b
<
=
x
f
(
a
)
∗
g
(
b
)
=
∑
1
<
=
a
<
=
x
g
(
a
)
∗
s
u
m
f
(
⌊
x
/
a
⌋
)
∴
s
u
m
f
(
x
)
∗
g
(
1
)
=
∑
a
∗
b
<
=
x
f
(
a
)
∗
g
(
b
)
−
∑
1
<
=
a
<
x
g
(
a
)
∗
s
u
m
f
(
⌊
x
/
a
⌋
)
\because \sum_{a*b<=x} f(a) * g(b) = \sum_{1<=a<=x} g(a) * sumf(\lfloor x/a \rfloor) \\ \therefore sumf(x) * g(1) = \sum_{a*b<=x} f(a) * g(b) - \sum_{1<=a<x} g(a) * sumf(\lfloor x/a \rfloor)
∵a∗b<=x∑f(a)∗g(b)=1<=a<=x∑g(a)∗sumf(⌊x/a⌋)∴sumf(x)∗g(1)=a∗b<=x∑f(a)∗g(b)−1<=a<x∑g(a)∗sumf(⌊x/a⌋)
只要能快速算出右边的且
g
(
1
)
≠
0
g(1)\neq 0
g(1)=0就能求解。
先放一个大佬的水表
http://www.cnblogs.com/cjyyb/category/1124706.html
1.莫比乌斯函数之和
g
=
I
g = I
g=I
ACCode:
#include<bits/stdc++.h>
#define LL long long
using namespace std;
/*
g(x) = 1
f(x) * g(x) = [x == 1]
sum(f(x) * g(x)) = 1;
*/
map<LL,int>mp;
LL n,m;
#define ts 10000005
LL mu[ts];
int pr[ts/8],cnt_pr;
bool vis[ts];
void sieve()
{
mu[1] = 1;
for(int i=2,lim = min(ts*1ll,m+1);i<lim;i++)
{
if(!vis[i]) pr[cnt_pr++] = i , mu[i] = -1;
for(int j=0;j<cnt_pr && pr[j] * i < lim;j++)
{
vis[i * pr[j]] = 1;
if(i % pr[j] == 0) break;
mu[i * pr[j]] = -mu[i];
}
mu[i] += mu[i-1];
}
}
LL solve(LL n)
{
if(n <= ts) return mu[n];
LL tmp;
if(mp.count(n)) return mp[n];
tmp = 1;
for(LL i=2,nxt;i<=n;i=nxt+1)
{
nxt = n / (n / i);
tmp = (tmp - (nxt - i + 1) * solve(n/i));
}
mp[n] = tmp;
return tmp;
}
int main()
{
scanf("%lld%lld",&n,&m);
sieve();
printf("%lld\n",solve(m)-solve(n-1));
}
2.欧拉函数之和。
g
=
I
,
s
u
m
(
f
∗
g
)
(
n
)
=
n
∗
(
n
+
1
)
/
2
g = I , sum(f * g)(n) = n * (n+1) / 2
g=I,sum(f∗g)(n)=n∗(n+1)/2
#include<bits/stdc++.h>
#define LL long long
#define mod 1000000007
using namespace std;
/*
g(x) = 1
f(x) * g(x) = [x == 1]
sum(f(x) * g(x)) = 1;
*/
LL n,m;
#define ts 5000005
LL phi[ts],phi2[100005];
int pr[ts/8],cnt_pr;
bool vis[ts];
void sieve()
{
phi[1] = 1;
for(int i=2,lim = min(ts*1ll,n+1);i<lim;i++)
{
if(!vis[i]) pr[cnt_pr++] = i , phi[i] = i-1;
for(int j=0;j<cnt_pr && pr[j] * i < lim;j++)
{
vis[i * pr[j]] = 1;
if(i % pr[j] == 0){ phi[i*pr[j]] = phi[i] * pr[j];break;}
phi[i * pr[j]] = phi[i] * (pr[j] - 1);
}
phi[i] += phi[i-1];
}
}
LL N;
LL solve(LL n)
{
if(n <= ts) return phi[n];
int loc = N / n;
if(phi2[loc] != phi2[0]) return phi2[loc];
LL tmp = ((n&1) ? (n%mod*((n+1)/2)) : (n/2%mod * ((n+1)%mod))) % mod;
for(LL i=2,nxt;i<=n;i=nxt+1)
{
nxt = n / (n / i);
tmp = (tmp - (nxt - i + 1) % mod * solve(n/i)) % mod;
}
phi2[loc] = tmp;
return tmp;
}
int main()
{
scanf("%lld",&n);
sieve();
memset(phi2,0x7f,sizeof phi2);
N = n;
printf("%lld\n",(solve(n)+mod)%mod);
}
bzoj 3944 sum
LL 好恶心。
常数很大。
#include<bits/stdc++.h>
#define LL long long
using namespace std;
LL n;
#define ts 5000005
int mu[ts],pr[ts/3],cnt_pr;
bool vis[ts];
LL phi[ts];
void sieve()
{
mu[1] = phi[1] = 1;
for(int i=2;i<ts;i++)
{
if(!vis[i]) pr[cnt_pr++] = i, mu[i] = -1,phi[i] = i-1;
for(int j=0;j<cnt_pr && pr[j] * i < ts;j++)
{
vis[i * pr[j]] = 1;
if(i % pr[j] == 0){ phi[i*pr[j]] = phi[i] * pr[j]; break; }
phi[i*pr[j]] = phi[i] * (pr[j] - 1);
mu[i * pr[j]] = -mu[i];
}
mu[i] += mu[i-1];
phi[i] += phi[i-1];
}
}
namespace Phi
{
LL f[100005];
int vis[100005],tot;
LL solve(LL N)
{
if(N < ts) return phi[N];
int loc = n / N;
if(vis[loc] == tot) return f[loc];
f[loc] = N * (N+1) / 2;
for(int i=2,nxt=0;i<=N && nxt<N;i=nxt+1)
{
nxt = N/(N/i);
f[loc] = f[loc] - (nxt - i + 1) * solve(N/i);
}
vis[loc] = tot;
return f[loc];
}
}
namespace Mu
{
LL f[100005];
int vis[100005],tot;
LL solve(LL N)
{
if(N<ts) return mu[N];
int loc = n / N;
if(vis[loc] == tot) return f[loc];
f[loc] = 1;
for(int i=2,nxt=0;i<=N && nxt<N;i=nxt+1)
{
nxt = N / (N / i);
f[loc] = f[loc] - (nxt - i + 1) * solve(N/i);
}
vis[loc] = tot;
return f[loc];
}
}
int main()
{
sieve();
int T;
for(scanf("%d",&T);T--;)
{
scanf("%lld",&n);
Mu::tot++;
Phi::tot++;
printf("%lld %lld\n",Phi::solve(n),Mu::solve(n));
}
}
hdu 5608 function
我主要说一下，
∑
i
=
1
x
n
2
+
3
n
+
2
=
x
∗
(
x
−
1
)
∗
(
x
−
2
)
/
3
=
x
3
‾
/
3
\sum_{i=1}^xn^2+3n+2 = x * (x-1) * (x-2) / 3 = x^{\underline3}/3
∑i=1xn2+3n+2=x∗(x−1)∗(x−2)/3=x3/3
然后这种和因数有关的函数可以
O
(
n
l
o
g
n
)
O(nlogn)
O(nlogn)筛。。。。
我用的
O
(
n
n
)
O(n\sqrt{n})
O(nn)
这个是一个离散微积分的基本结论
#include<bits/stdc++.h>
#define LL long long
#define mod 1000000007
using namespace std;
#define ts 100005
LL f[100005];
map<int,int>mp;
int solve(int n)
{
if(n < ts) return f[n];
if(mp.count(n)) return mp[n];
LL ret = 1ll * (n-2) * (n-1) % mod * n % mod * 333333336ll % mod;
for(int i=2,nxt;i<=n;i=nxt+1)
{
nxt = n/(n/i);
ret =(ret - 1ll * (nxt - i + 1) * solve(n / i)) % mod;
}
ret = (ret%mod+mod)%mod;
mp[n] = ret;
return ret;
}
void sieve()
{
for(int i=2;i<ts;i++)
{
f[i] = (1ll * i * i - 3ll * i + 2) % mod;
for(int j=2;j*j<=i;j++)
if(i % j == 0)
{
f[i] -= f[j];
if(j * j != i) f[i] -= f[i/j];
}
f[i] =(f[i] % mod + mod) % mod;
}
for(int i=2;i<ts;i++)
f[i] = ((f[i] + f[i-1]) % mod + mod ) % mod;
}
int main()
{
sieve();
int T;
for(scanf("%d",&T);T--;)
{
int n;
scanf("%d",&n);
printf("%d\n",solve(n));
}
}
51nod 1237
给出一个数N，输出小于等于N的所有数，两两之间的最大公约数之和。
相当于计算这段程序（程序中的gcd(i,j)表示i与j的最大公约数）：
由于结果很大，输出Mod 1000000007的结果。
G = 0 ; f o r ( i = 1 ; i < N ; i + + ) f o r ( j = 1 ; j < = N ; j + + ) G = ( G + g c d ( i , j ) ) N < = 1 0 1 0 G=0;\\ for(i=1;i<N;i++)\\ for(j=1;j<=N;j++)\\ {\\ G = (G + gcd(i,j)) % 1000000007;\\ }\\ N<=10^10\\ G=0;for(i=1;i<N;i++)for(j=1;j<=N;j++)G=(G+gcd(i,j))N<=1010
这个题揭露了杜教筛的另一个性质：
实际上求
s
u
m
f
(
n
)
sumf(n)
sumf(n)用到的
s
u
m
f
(
x
)
sumf(x)
sumf(x)满足
x
=
⌊
n
k
⌋
,
k
∈
N
∗
x = \lfloor \frac nk \rfloor , k \in N^*
x=⌊kn⌋,k∈N∗
根据整除分块可以知道只有
O
(
n
)
O(\sqrt n)
O(n)个x。
这个题的答案是：
2
∗
∑
i
=
1
n
i
∗
∑
j
=
1
⌊
n
i
⌋
φ
(
j
)
−
∑
i
=
1
n
i
2 * \sum_{i=1}^ni *\sum_{j=1}^{\lfloor \frac ni \rfloor} \varphi(j) - \sum_{i=1}^ni
2∗i=1∑ni∗j=1∑⌊in⌋φ(j)−i=1∑ni
这个分块优化套杜教筛之后也只会用到
O
(
n
)
O(\sqrt n)
O(n)个
s
u
m
φ
(
x
)
sum\varphi(x)
sumφ(x)
所以还是
O
(
n
2
3
)
O(n^{\frac 23})
O(n32)的
（我代码中的记忆化存储方法也可以用这个结论证明是不会撞的，尽管没快多少，除法的常数。。。。。）
AC Code:
#include<bits/stdc++.h>
#define LL long long
using namespace std;
/*
f(d) = d * g(n/d , m/d)
= d * sigma((1<=i<=n/d , 1<=j<=n/d) [gcd(i,j)==1])
= d * (sigma((1<=i<=n/d , phi(i)) * 2) - d
sigma(1<=d<=n,f(d))
= 2 * sigma(1<=i<=n,phi(i) * (n / i)) - (n+1) * n / 2
= 2 * sigma(1<=i<=n,i * sigma(k<=n/i,phi(k))) - (n+1) * n / 2
*/
LL n;
#define ts 5000005
#define mod 1000000007
LL calc(LL a,LL b)
{
LL c = a + b;
LL d = b - a + 1;
if(c & 1) d/=2;
else c/=2;
c%=mod,d%=mod;
return d * c % mod ;
}
int pr[ts],cnt_pr;
bool vis[ts];
LL phi[ts],f[ts];
void sieve()
{
phi[1] = 1;
for(int i=2;i<ts;i++)
{
if(!vis[i]) pr[cnt_pr++] = i , phi[i] = i-1;
for(int j=0;j<cnt_pr && pr[j] * i < ts;j++)
{
vis[i * pr[j]] = 1;
if(i % pr[j] == 0){ phi[i * pr[j]] = phi[i] * pr[j]; break; }
phi[i * pr[j]] = phi[i] * (pr[j] - 1);
}
phi[i] = (phi[i] + phi[i-1]) % mod;
}
}
LL solve(LL m)
{
if(m < ts) return phi[m];
int loc = n / m;
if(f[loc]!=-1) return f[loc];
int tmp = calc(1,m);
for(LL i=2,nxt;i<=m;i=nxt+1)
{
nxt = m/(m/i);
tmp = (tmp - (nxt - i + 1) % mod * solve(m/i)) % mod;
}
return f[loc] = (tmp + mod) % mod;
}
int main()
{
sieve();
memset(f,-1,sizeof f);
scanf("%lld",&n);
LL ans = 0;
for(LL i=1,nxt;i<=n;i=nxt+1)
{
nxt = n/(n/i);
ans = (ans + calc(i,nxt) * solve(n/i)) % mod;
}
ans = (2ll * ans) % mod;
ans = (ans - calc(1,n)) % mod;
printf("%lld\n",(ans+mod)%mod);
}
51nod 1238
这个题整死我了。
在这个题中，我们能以更一般性的视角来看杜教筛，再也不是拿个
I
I
I就来卷就行的了。
出一个数N，输出小于等于N的所有数，两两之间的最小公倍数之和。
相当于计算这段程序（程序中的lcm(i,j)表示i与j的最小公倍数）：
由于结果很大，输出Mod 1000000007的结果。
G=0;
for(i=1;i<N;i++)
for(j=1;j<=N;j++)
{
G = (G + lcm(i,j)) % 1000000007;
}
在jzptab（和这个题只有数据范围不同的一个题）中，我们可以推出这个式子
a
n
s
=
∑
d
=
1
n
(
⌊
n
d
⌋
∗
(
⌊
n
d
⌋
+
1
)
/
2
)
2
∗
∑
p
∣
d
μ
(
p
)
∗
p
∗
p
∗
d
p
ans = \sum_{d=1}^n\left(\lfloor\frac nd\rfloor * (\lfloor\frac nd\rfloor + 1) / 2 \right )^2 * \sum_{p|d}\mu(p) * p * p* \frac dp
ans=d=1∑n(⌊dn⌋∗(⌊dn⌋+1)/2)2∗p∣d∑μ(p)∗p∗p∗pd
前面这个其实是不好动的，但是我们可以对前面的分块优化。
那么只需要求后面的前缀和就行了。（这里换一下正规的符号)
f
∘
g
=
∑
x
∑
p
∣
x
f
(
p
)
∗
g
(
x
p
)
f
⋅
g
=
∑
x
f
(
x
)
×
g
(
x
)
f \circ g = \sum_x\sum_{p|x}f(p) * g(\frac xp)\\ f\cdot g = \sum_x f(x) \times g(x)
f∘g=x∑p∣x∑f(p)∗g(px)f⋅g=x∑f(x)×g(x)
后面的函数其实就是
f
(
x
)
=
(
i
d
⋅
i
d
⋅
μ
)
∘
i
d
f(x) = (id \cdot id \cdot \mu) \circ id
f(x)=(id⋅id⋅μ)∘id
然后认真研究一下这个函数
你会发现：
(
i
d
⋅
i
d
⋅
μ
)
∘
(
i
d
⋅
i
d
)
=
∑
x
∑
p
∣
x
p
∗
p
∗
μ
(
p
)
∗
x
p
∗
x
p
=
∑
x
x
2
∗
∑
p
∣
x
μ
(
p
)
(id \cdot id \cdot \mu) \circ (id \cdot id) = \sum_x\sum_{p|x}p*p*\mu(p)*\frac xp*\frac xp = \sum_x x^2 * \sum_{p|x} \mu(p)
(id⋅id⋅μ)∘(id⋅id)=x∑p∣x∑p∗p∗μ(p)∗px∗px=x∑x2∗p∣x∑μ(p)
就是元函数嘛。。。。（这启示我们，寻找其他函数的逆元很重要，不要只知道
μ
\mu
μ和
I
I
I有，通过
∘
\circ
∘来把
⋅
\cdot
⋅提出sigma也是很重要的）
然后把
f
(
x
)
f(x)
f(x)和
i
d
⋅
i
d
id\cdot id
id⋅id卷起来杜教筛就是了。
注意杜教筛中如
g
(
1
)
≠
1
g(1)\neq 1
g(1)=1还得求逆元（此题没有）
另外求前缀和的时候如果想卡常，
s
u
m
(
n
x
t
)
−
s
u
m
(
i
−
1
)
sum(nxt) - sum(i-1)
sum(nxt)−sum(i−1)的后部分用一个变量存，初始时，杜教筛是从2开始的，初始值要小心。。。。。
网上别的解法用了
∑
i
=
1
d
[
(
i
,
d
)
=
=
1
]
∗
i
=
ϕ
(
d
)
∗
d
+
[
d
=
1
]
2
\sum_{i=1}^d [(i,d)==1]*i = \frac {\phi(d) * d + [d=1]} {2}
i=1∑d[(i,d)==1]∗i=2ϕ(d)∗d+[d=1]
也很巧
AC Code:
#include<cstdio>
#include<algorithm>
#include<cstring>
#define LL long long
#define mod 1000000007
using namespace std;
/*
sigma(i,j,i * j / gcd(i,j))
= sigma(d , d * sigma(i,j,[gcd(i,j)==1] * i * j))
= sigma(d , d * sigma(p,mu[p] * S(n/d/p) * p * p))
= sigma(D , sigma(mu[p] * S(n/D) * p * p * D/p))
= sigma(D , S(n/D) * D * sigma(mu[p] * p));
f(x) = sigma(p|x , mu[p] * p * p * (x/p)) = (id.id.mu)*id
f*(id.id)
= (id . id . mu) * id * (id . id)
= (id . id . mu) * (id . id) * id
= sigma(x , sigma(p | x , p * p * mu[p] * x / p * x / p) ) * id
= sigma(x , x * x * sigma(p | x , mu[p]) ) * id
= sigma(x , x * x * [x==1]) * id
= id
sigf(n) * (id.id)(1) = sigma(i<=n,f*(id.id)(i)) - sigma(2<=i<=n,(id.id)(i) * sigf(n/i))
id.id(x) = x*x
*/
LL n;
#define ts 5000005
LL pf[ts],f[ts],pr[ts],cnt_pr;
bool vis[ts];
void sieve()
{
pf[1] = 1;
for(LL i=2;i<ts;i++)
{
if(!vis[i]) pr[cnt_pr++] = i , pf[i] = 1 - i;
for(LL j=0;j<cnt_pr && pr[j] * i < ts;j++)
{
vis[i * pr[j]] = 1;
if(i % pr[j] == 0)
{
pf[i * pr[j]] = pf[i];
break;
}
pf[i * pr[j]] = 1ll * pf[i] * pf[pr[j]] % mod;
}
}
for(LL i=2;i<ts;i++)
pf[i] =( (1ll * pf[i] * i % mod + pf[i-1]) % mod + mod) % mod;
}
LL cal2(LL a,LL b)
{
LL c = a+b , d = b-a+1;
if(c & 1) d/=2;else c/=2;
return (d % mod) * (c % mod) % mod;
}
LL cap2(LL m)
{
int tmp;
if(m>ts)
tmp = m;
LL a[3] = {m , m+1 , 2 * m + 1};
for(LL i=0;i<3;i++) if(a[i] % 2 == 0)
{
a[i] /= 2;
break;
}
for(LL i=0;i<3;i++) if(a[i] % 3 == 0)
{
a[i] /= 3;
break;
}
return (a[0] % mod) * (a[1] % mod) % mod * (a[2] % mod) % mod;
}
inline LL Pow(LL base,LL k)
{
LL ret = 1;
for(;k;k>>=1,base=1ll*base*base % mod) if(k & 1) ret = 1ll * ret * base % mod;
return ret;
}
LL Solve(LL m)
{
if(m < ts) return pf[m];
LL loc = n / m;
if(f[loc] != -1) return f[loc];
LL tmp = cal2(1,m);
for(LL i=2,nxt,pre=1,tpre;i<=m;i=nxt+1)
{
nxt = m / (m / i);
tmp = (tmp - ((tpre = cap2(nxt)) - pre) % mod * Solve(m / i) % mod) % mod;
pre = tpre;
}
return f[loc] = (tmp + mod) % mod;
}
LL S(LL m)
{
return cal2(1,m) * cal2(1,m) % mod;
}
int main()
{
sieve();
scanf("%lld",&n);
memset(f,-1,sizeof f);
LL ans = 0;
for(LL i=1,nxt,pre=0,rpre;i<=n;i=nxt+1)
{
nxt = n/(n/i);
ans = (ans + ((rpre = Solve(nxt)) - Solve(i-1)) % mod * S(n/i) % mod) % mod;
pre = rpre;
}
printf("%lld\n",(ans+mod) % mod);
}
51 nod 1227
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 s = ∑ 1 < = i < = j < = n i / g c d ( i , j ) = ∑ 1 < = d < = n ∑ 1 < = i < = j < = n / d [ g c d ( i , j ) = = 1 ] ∗ i = ∑ 1 < = d < = n ∑ 1 < = p < = n / d μ ( p ) ∗ p ∗ S ( n / d / p ) = ∑ 1 < = D < = n S ( n / D ) ∗ ∑ p ∣ D μ ( p ) ∗ p S ( n ) = ∑ 1 < = i < = n i ∗ ( n − i + 1 ) = ( n + 1 ) ∗ ( n + 1 ) ∗ n / 2 − n ∗ ( n + 1 ) ∗ ( 2 n + 1 ) / 6 = ( n + 2 ) ∗ ( n + 1 ) ∗ n 6 f ( n ) = ∑ p ∣ n μ ( p ) ∗ p f ∘ i d ( x ) = I \begin{aligned} ans =& \sum_{1<=i<=j<=n} i / gcd(i,j)\\ =&\sum_{1<=d<=n}\sum_{1<=i<=j<=n/d} [gcd(i,j)==1] * i\\ =&\sum_{1<=d<=n}\sum_{1<=p<=n/d} \mu(p) * p * S(n/d/p)\\ =&\sum_{1<=D<=n}S(n/D) * \sum_{p|D} \mu(p) * p\\ S(n) = &\sum_{1<=i<=n}i*(n-i+1) \\ =&(n+1) * (n+1) * n / 2 - n * (n+1) * (2n+1) / 6\\ =& \frac {(n+2) * (n+1) * n}{6}\\ f(n) = &\sum_{p|n} \mu(p) * p\\ f \circ &id (x)= I \end{aligned} ans====S(n)===f(n)=f∘1<=i<=j<=n∑i/gcd(i,j)1<=d<=n∑1<=i<=j<=n/d∑[gcd(i,j)==1]∗i1<=d<=n∑1<=p<=n/d∑μ(p)∗p∗S(n/d/p)1<=D<=n∑S(n/D)∗p∣D∑μ(p)∗p1<=i<=n∑i∗(n−i+1)(n+1)∗(n+1)∗n/2−n∗(n+1)∗(2n+1)/66(n+2)∗(n+1)∗np∣n∑μ(p)∗pid(x)=I
一个括号啊，我调了40min
#include<bits/stdc++.h>
#define LL long long
#define mod 1000000007
using namespace std;
LL a,b;
map<LL,int>mp;
#define ts 5000005
int mu[ts] , pr[ts / 8] , cnt_pr;
bool vis[ts];
void sieve()
{
mu[1] = 1;
for(int i=2;i<ts;i++)
{
if(!vis[i]) pr[cnt_pr++] = i , mu[i] = 1-i;
for(int j=0;j<cnt_pr && pr[j] * i < ts;j++)
{
vis[i * pr[j]] = 1;
if(i % pr[j] == 0)
{
mu[i * pr[j]] = mu[i];
break;
}
mu[i * pr[j]] = (1ll * mu[i] * mu[pr[j]]) % mod;
}
}
for(int i=2;i<ts;i++) mu[i] = (mu[i] + mu[i-1]) % mod;
}
LL calc(LL a,LL b)
{
LL c = a+b;
b = b-a+1;
if(c & 1) b/=2; else c/=2;
return (c % mod) * (b % mod) % mod;
}
LL S(LL n)
{
LL a[3] = {n,n+1,n+2};
for(int i=0;i<3;i++) if(a[i] % 2 == 0) {a[i]/=2;break;}
for(int i=0;i<3;i++) if(a[i] % 3 == 0) {a[i]/=3;break;}
return (a[0] % mod) * (a[1] % mod) % mod * (a[2] % mod) % mod;
}
LL solve(LL n)
{
if(n<ts) return mu[n];
if(mp.count(n)) return mp[n];
int tmp = n % mod;
for(LL i=2,nxt;i<=n;i=nxt+1)
{
nxt = n/(n/i);
tmp = (tmp - calc(i,nxt) * solve(n/i)) % mod;
}
return mp[n] = (tmp + mod) % mod;
}
LL As(LL n)
{
LL ans = 0;
for(LL i=1,nxt,pre=0,rpre;i<=n;i=nxt+1)
{
nxt = n / (n / i);
ans = (ans + ((rpre = solve(nxt)) - pre) * S(n/i)) % mod;
pre = rpre;
}
return ans;
}
int main()
{
sieve();
scanf("%lld%lld",&a,&b);
printf("%lld\n",(As(b) - As(a-1) + mod) % mod);
}
51 nod 1222
A
n
s
=
∑
1
<
=
i
<
=
j
<
=
n
[
i
j
g
c
d
(
i
,
j
)
<
=
n
]
=
∑
1
<
=
d
<
=
n
∑
1
<
=
i
<
=
j
<
=
n
d
[
i
j
d
<
=
n
]
∗
[
g
c
d
(
i
,
j
)
=
=
1
]
=
∑
1
<
=
d
<
=
n
∑
1
<
=
p
<
=
n
d
μ
(
p
)
∑
1
<
=
i
<
=
j
<
=
n
d
p
[
i
j
d
<
=
n
p
2
]
=
∑
1
<
=
p
<
=
n
μ
(
p
)
∗
∑
1
<
=
d
<
=
n
,
1
<
=
i
<
=
j
<
=
n
[
i
j
d
<
=
n
p
2
]
=
n
+
∑
1
<
=
p
<
=
n
μ
(
p
)
∗
∑
1
<
=
d
,
i
,
j
<
=
n
[
i
j
d
<
=
n
p
2
]
2
\begin{aligned} Ans=&\sum_{1<=i<=j<=n}[\frac {ij}{gcd(i,j)} <= n ]\\ =&\sum_{1<=d<=n}\sum_{1<=i<=j<=\frac nd}[ijd<=n] * [gcd(i,j)==1]\\ =&\sum_{1<=d<=n}\sum_{1<=p<=\frac nd} \mu(p)\sum_{1<=i<=j<=\frac n{dp}}[ijd<=\frac n{p^2}]\\ =&\sum_{1<=p<=n}\mu(p) * \sum_{1<=d<=n,1<=i<=j<=n} [ijd<=\frac n{p^2}]\\ =&\frac{n + \sum_{1<=p<=n}\mu(p) * \sum_{1<=d,i,j<=n} [ijd<=\frac n{p^2}]}2 \end{aligned}
Ans=====1<=i<=j<=n∑[gcd(i,j)ij<=n]1<=d<=n∑1<=i<=j<=dn∑[ijd<=n]∗[gcd(i,j)==1]1<=d<=n∑1<=p<=dn∑μ(p)1<=i<=j<=dpn∑[ijd<=p2n]1<=p<=n∑μ(p)∗1<=d<=n,1<=i<=j<=n∑[ijd<=p2n]2n+∑1<=p<=nμ(p)∗∑1<=d,i,j<=n[ijd<=p2n]
可以for循环
O
(
n
2
3
)
O(n^{\frac 23})
O(n32)统计。。。
AC Code:
#include<bits/stdc++.h>
#define LL long long
using namespace std;
#define ts 1000005
int pr[ts/8],cnt_pr,mu[ts];
bool vis[ts];
void sieve()
{
mu[1] = 1;
for(int i=2;i<ts;i++)
{
if(!vis[i]) pr[cnt_pr++] = i, mu[i] = -1;
for(int j=0;pr[j] * i < ts;j++)
{
vis[i * pr[j]] = 1;
if(i % pr[j] == 0) break;
mu[i * pr[j]] = -mu[i];
}
}
}
LL solve(LL n)
{
LL ans = 0;
for(LL k=1;k*k<=n;k++)
if(mu[k])
{
LL m = n / (k * k) , s = 0;
for(LL i=1;i*i*i<=m;i++)
{
for(LL j=i+1;i*j*j<=m;j++)
s+=(m / (i * j) - j) * 6 + 3;
s+=(m/(i*i)-i)*3+1;
}
ans += mu[k] * s;
}
return (ans + n) / 2;
}
int main()
{
sieve();
LL a,b;
scanf("%lld%lld",&a,&b);
printf("%lld\n",solve(b) - solve(a-1));
}
SPOJ DIVCNT2 - Counting Divisors (square)
3个结论：
∑
p
∣
d
μ
(
p
)
2
=
2
ω
(
p
)
,
ω
(
p
)
=
p
的
质
因
子
个
数
\sum_{p|d}\mu(p)^2 = 2^{\omega\left(p\right)},\omega(p) = p的质因子个数
p∣d∑μ(p)2=2ω(p),ω(p)=p的质因子个数
∑
1
<
=
i
<
=
n
μ
(
p
)
2
=
∑
1
<
=
i
<
=
n
μ
[
i
]
∗
⌊
n
i
2
⌋
\sum_{1<=i<=n}\mu(p)^2 = \sum_{1<=i<=\sqrt n}\mu[i] * \left \lfloor \frac n {i^2} \right \rfloor
1<=i<=n∑μ(p)2=1<=i<=n∑μ[i]∗⌊i2n⌋
∑
1
<
=
i
,
j
<
=
n
[
i
j
<
=
n
]
=
∑
1
<
=
i
<
=
n
⌊
n
i
⌋
\sum_{1<=i,j<=n}[ij<=n] = \sum_{1<=i<=n}\left \lfloor \frac ni \right \rfloor
1<=i,j<=n∑[ij<=n]=1<=i<=n∑⌊in⌋
然后就完了？
AC Code:
#include<bits/stdc++.h>
#define LL long long
using namespace std;
#define ts 10000000
int mu[ts],pr[ts/8],cnt_pr,lp[ts],u[ts];
LL d[ts];
bool vis[ts];
void sieve()
{
u[1] = mu[1] = d[1] = 1;
for(int i=2;i<ts;i++)
{
if(!vis[i]) pr[cnt_pr++] = i , mu[i] = 1 , d[i] = 2 , lp[i] = 2 , u[i] = -1;
for(int j=0;pr[j] * i < ts;j++)
{
vis[i * pr[j]] = 1;
if(i % pr[j] == 0)
{
d[i * pr[j]] = d[i] / lp[i] * (lp[i * pr[j]] = lp[i] + 1);
break;
}
u[i * pr[j]] = -u[i];
mu[i * pr[j]] = mu[i];
d[i * pr[j]] = d[i] * 2;
lp[i * pr[j]] = 2;
}
mu[i] += mu[i-1];
d[i] += d[i-1];
}
}
inline LL sum(LL n)
{
if(n < ts) return mu[n];
LL ans = 0;
for(LL j=1;j*j<=n;j++)
if(u[j])
ans += u[j] * (n / (j * j));
return ans;
}
inline LL sum2(LL n)
{
if(n < ts) return d[n];
LL ans = 0;
for(LL i=1,nxt;i<=n;i=nxt+1)
ans += ((nxt=(n/(n/i)))-i+1) * (n/i);
return ans;
}
int main()
{
sieve();
int T;
for(scanf("%d",&T);T--;)
{
LL n,ans=0;
scanf("%lld",&n);
for(LL i=1,nxt,pre=0,rpre;i<=n;i=nxt+1)
{
nxt = n / (n / i);
ans += ((rpre = sum(nxt)) - pre) * sum2(n/i);
pre = rpre;
}
printf("%lld\n",ans);
}
}
51nod 1220 约数之和
d(k)表示k的所有约数的和。d(6) = 1 + 2 + 3 + 6 = 12。
定义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)的结果。
∑
1
<
=
i
,
j
<
=
n
d
(
i
j
)
=
∑
1
<
=
i
,
j
<
=
n
∑
a
∣
i
,
b
∣
j
[
g
c
d
(
a
,
b
)
=
1
]
=
∑
1
<
=
i
,
j
<
=
n
∑
p
∑
a
∣
i
,
b
∣
j
μ
[
p
]
∗
[
p
∣
g
c
d
(
a
,
b
)
]
=
∑
1
<
=
i
,
j
<
=
n
⌊
n
i
⌋
⌊
n
j
⌋
∗
∑
p
∣
g
c
d
(
i
,
j
)
μ
[
p
]
=
∑
1
<
=
p
<
=
n
μ
(
p
)
∑
1
<
=
i
<
=
n
p
⌊
n
p
i
⌋
∑
1
<
=
j
<
=
n
p
⌊
n
p
j
⌋
=
∑
1
<
=
p
<
=
n
μ
(
p
)
G
(
n
p
)
2
G
=
I
∗
I
I
∗
I
∗
μ
=
I
\begin{aligned} \sum_{1<=i,j<=n}d(ij) =& \sum_{1<=i,j<=n}\sum_{a|i,b|j}[gcd(a,b) = 1]\\ =&\sum_{1<=i,j<=n}\sum_{p}\sum_{a|i,b|j} \mu[p] * [p | gcd(a,b)]\\ =&\sum_{1<=i,j<=n}\lfloor\frac ni\rfloor\lfloor\frac nj\rfloor*\sum_{p|gcd(i,j)} \mu[p]\\ =&\sum_{1<=p<=n}\mu(p) \sum_{1<=i<=\frac np}\lfloor\frac n{pi}\rfloor\sum_{1<=j<=\frac np}\lfloor\frac n{pj}\rfloor\\ =&\sum_{1<=p<=n}\mu(p) G(\frac np)^2\\ G = &I * I\\ I*I*\mu=&I \end{aligned}
1<=i,j<=n∑d(ij)=====G=I∗I∗μ=1<=i,j<=n∑a∣i,b∣j∑[gcd(a,b)=1]1<=i,j<=n∑p∑a∣i,b∣j∑μ[p]∗[p∣gcd(a,b)]1<=i,j<=n∑⌊in⌋⌊jn⌋∗p∣gcd(i,j)∑μ[p]1<=p<=n∑μ(p)1<=i<=pn∑⌊pin⌋1<=j<=pn∑⌊pjn⌋1<=p<=n∑μ(p)G(pn)2I∗II
这是约数个数之和，BZOJ 4176 Lucas的数论。
真是个悲情的故事，还是放一下代码吧。
UPD：放个结论：
σ
k
(
i
j
)
=
∑
a
∣
i
,
b
∣
j
[
(
a
,
b
)
=
1
]
(
a
j
b
)
k
\sigma_k(ij) = \sum_{a|i,b|j} [(a,b)=1](\frac {aj}b)^k
σk(ij)=∑a∣i,b∣j[(a,b)=1](baj)k
AC Code：
#include<bits/stdc++.h>
#define LL long long
#define mod 1000000007
using namespace std;
#define ts 1000006
int mu[ts],lp[ts],pr[ts],cnt_pr;
LL d[ts];
bool vis[ts];
void sieve()
{
mu[1] = d[1] = lp[1] = 1;
for(int i=2;i<ts;i++)
{
if(!vis[i]) pr[cnt_pr++] = i , d[i] = lp[i] = 2 , mu[i] = -1;
for(int j=0;pr[j] * i < ts;j++)
{
vis[i * pr[j]] = 1;
if(i % pr[j] == 0)
{
d[i * pr[j]] = d[i] / lp[i] * (lp[i * pr[j]] = lp[i] + 1);
break;
}
d[i * pr[j]] = d[i] * 2;
mu[i * pr[j]] = -mu[i];
lp[i * pr[j]] = 2;
}
d[i] = (d[i] + d[i-1]) % mod;
mu[i] = (mu[i] + mu[i-1]) % mod;
}
}
map<LL,int>mmu,md;
LL smu(LL n)
{
if(n < ts) return mu[n];
if(mmu.count(n)) return mmu[n];
int tmp = 1;
for(LL i=2,nxt;i<=n;i=nxt+1)
{
nxt = n / (n / i);
tmp = (tmp - (nxt - i + 1) % mod * smu(n/i)) % mod;
}
return mmu[n] = (tmp + mod)%mod;
}
LL sd(LL n)
{
if(n < ts) return d[n];
if(md.count(n)) return md[n];
int tmp = n % mod;
for(LL i=2,nxt,pre=1,rpre;i<=n;i=nxt+1)
{
nxt = n / (n / i);
tmp = (tmp - ((rpre = smu(nxt)) - pre) * sd(n/i)) % mod;
pre = rpre;
}
return md[n] = tmp;
}
int main()
{
sieve();
LL n,ans=0;
scanf("%lld",&n);
for(LL i=1,nxt,pre=0,rpre;i<=n;i=nxt+1)
{
nxt = n / (n / i);
ans = (ans + ((rpre = smu(nxt)) - pre) * sd(n/i) %mod *sd(n/i)) % mod;
pre = rpre;
}
printf("%lld\n",(ans+mod)%mod);
}
其实对于约数之和也有一个类似的推法。
∑
1
<
=
i
,
j
<
=
n
d
(
i
j
)
=
∑
1
<
=
i
,
j
<
=
n
∑
a
∣
i
,
b
∣
j
a
j
[
g
c
d
(
a
,
b
)
=
=
1
]
b
=
∑
1
<
=
i
,
j
<
=
n
∑
a
∣
i
,
b
∣
j
a
j
∑
p
∣
g
c
d
(
a
,
b
)
μ
(
p
)
b
先枚举a和b
=
∑
1
<
=
a
,
b
<
=
n
a
∗
⌊
n
a
⌋
∗
(
b
+
b
⌊
n
b
⌋
)
∗
⌊
n
b
⌋
2
∑
p
∣
g
c
d
(
a
,
b
)
μ
(
p
)
b
=
∑
1
<
=
a
,
b
<
=
n
a
⌊
n
a
⌋
∗
(
1
+
⌊
n
b
⌋
)
∗
⌊
n
b
⌋
2
∗
∑
p
∣
g
c
d
(
a
,
b
)
μ
(
p
)
先枚举p再枚举gcd
S
(
n
)
=
(
1
+
n
)
∗
n
2
=
∑
1
<
=
p
<
=
n
μ
(
p
)
∑
1
<
=
a
<
=
n
p
a
p
⌊
n
a
p
⌋
∗
∑
1
<
=
b
<
=
n
p
S
(
b
)
=
∑
1
<
=
p
<
=
n
μ
(
p
)
p
(
∑
1
<
=
a
<
=
n
p
a
⌊
n
a
p
⌋
)
2
=
∑
1
<
=
p
<
=
n
μ
(
p
)
p
(
∑
1
<
=
a
<
=
n
p
σ
(
a
)
)
2
\begin{aligned} \sum_{1<=i,j<=n} d(ij) =& \sum_{1<=i,j<=n}\sum_{a|i,b|j} \frac {aj[gcd(a,b)==1]}{b}\\ =&\sum_{1<=i,j<=n}\sum_{a|i,b|j} \frac {aj \sum_{p|gcd(a,b)}\mu(p)}{b}\\ &\text{先枚举a和b}\\ =&\sum_{1<=a,b<=n} \frac{\frac {a * \lfloor \frac na \rfloor * (b + b\lfloor \frac nb \rfloor) * \lfloor \frac nb \rfloor } 2\sum_{p|gcd(a,b)}\mu(p)}{b}\\ =&\sum_{1<=a,b<=n} \frac{a\lfloor \frac na \rfloor * (1+\lfloor \frac nb \rfloor) * \lfloor \frac nb \rfloor}{2} * \sum_{p|gcd(a,b)}\mu(p)\\ &\text{先枚举p再枚举gcd}\\ &S(n) = \frac {(1 + n) * n}2\\ =&\sum_{1<=p<=n} \mu(p) \sum_{1<=a<=\frac np} ap\lfloor \frac n{ap} \rfloor * \sum_{1<=b<=\frac np} S(b)\\ =&\sum_{1<=p<=n} \mu(p) p\left(\sum_{1<=a<=\frac np} a\lfloor \frac n{ap} \rfloor \right)^2\\ =&\sum_{1<=p<=n} \mu(p) p\left(\sum_{1<=a<=\frac np} \sigma(a) \right)^2 \end{aligned}
1<=i,j<=n∑d(ij)=======1<=i,j<=n∑a∣i,b∣j∑baj[gcd(a,b)==1]1<=i,j<=n∑a∣i,b∣j∑baj∑p∣gcd(a,b)μ(p)先枚举a和b1<=a,b<=n∑b2a∗⌊an⌋∗(b+b⌊bn⌋)∗⌊bn⌋∑p∣gcd(a,b)μ(p)1<=a,b<=n∑2a⌊an⌋∗(1+⌊bn⌋)∗⌊bn⌋∗p∣gcd(a,b)∑μ(p)先枚举p再枚举gcdS(n)=2(1+n)∗n1<=p<=n∑μ(p)1<=a<=pn∑ap⌊apn⌋∗1<=b<=pn∑S(b)1<=p<=n∑μ(p)p⎝⎛1<=a<=pn∑a⌊apn⌋⎠⎞21<=p<=n∑μ(p)p⎝⎛1<=a<=pn∑σ(a)⎠⎞2
后面那个是约数和的前缀和的平方，为
i
d
∘
I
id\circ I
id∘I 卷一个
μ
\mu
μ就行。
前面那个卷一个
i
d
id
id就行
那么就可以筛了。
AC Code:
#include<bits/stdc++.h>
#define LL long long
#define mod 1000000007
using namespace std;
#define ts 5000005
int pr[ts],cnt_pr,mup[ts],sig[ts],mu[ts];
bool vis[ts];
void sieve()
{
mup[1] = 1 , sig[1] = 1 , mu[1] = 1;
for(int i=2;i<ts;i++)
{
if(!vis[i]) pr[cnt_pr++] = i , mu[i] = -1 , sig[i] = i+1;
for(int j=0;pr[j] * i < ts;j++)
{
vis[i * pr[j]] = 1;
if(i % pr[j] == 0)
{
sig[i * pr[j]] = (1ll * sig[i] * (pr[j] + 1) - 1ll * sig[i / pr[j]] * pr[j]) % mod;
break;
}
mu[i * pr[j]] = -mu[i];
sig[i * pr[j]] = 1ll * sig[i] * sig[pr[j]] % mod;
}
}
for(int i=2;i<ts;i++)
mup[i] = (1ll * mu[i] * i + mup[i-1]) % mod,
mu[i] = mu[i] + mu[i-1],
sig[i] = (sig[i] + sig[i-1]) % mod;
}
LL calc2(LL n)
{
static int inv2 = (mod+1) / 2;
return n % mod * ((n+1) % mod) % mod * inv2 % mod;
}
map<LL,LL>mmup,msig,mmu;
LL calmu(LL n)
{
if(n < ts) return mu[n];
if(mmu.count(n)) return mmu[n];
int tmp = 1;
for(LL i=2,nxt;i<=n;i=nxt+1)
{
nxt = n / (n / i);
tmp = (tmp - (nxt - i + 1) % mod * calmu(n/i)) % mod;
}
return mmu[n] = tmp;
}
LL calmup(LL n)
{
if(n < ts) return mup[n];
if(mmup.count(n)) return mmup[n];
int tmp = 1;
for(LL i=2,nxt;i<=n;i=nxt+1)
{
nxt = n / (n / i);
tmp = (tmp - (calc2(nxt) - calc2(i-1)) * calmup(n/i)) % mod;
}
return mmup[n] = tmp;
}
LL calsig(LL n)
{
if(n < ts) return sig[n];
if(msig.count(n)) return msig[n];
int tmp = calc2(n);
for(LL i=2,nxt;i<=n;i=nxt+1)
{
nxt = n / (n / i);
tmp = (tmp - (calmu(nxt) - calmu(i-1)) * calsig(n/i)) % mod;
}
return msig[n] = tmp;
}
inline LL sqr(LL a)
{
return a * a % mod;
}
int main()
{
sieve();
LL n;
scanf("%lld",&n);
int ans = 0;
for(LL i=1,nxt;i<=n;i=nxt+1)
{
nxt = n/(n/i);
ans = (ans + 1ll * (calmup(nxt) - calmup(i-1)) * sqr(calsig(n/i))) % mod;
}
printf("%d\n",(ans+mod) % mod);
}
打完之后发现自己好智障啊，大于
n
2
3
n^{\frac 23}
n32的可以直接
∑
1
<
=
i
<
=
n
σ
(
i
)
=
∑
1
<
=
i
<
=
n
⌊
n
i
⌋
\sum_{1<=i<=n}\sigma(i) = \sum_{1<=i<=n}\lfloor \frac ni \rfloor
1<=i<=n∑σ(i)=1<=i<=n∑⌊in⌋
然后分块优化做。。。。。。复杂度一样，我还多筛了一个
μ
\mu
μ
51nod 1584加权约数和
去年的 tangjz 非常喜欢做数论题，但是一年以后的 tangjz 却不那么会做了。
在整理以前的试题时，他发现了这样一道题目：“求 ∑σ(i) ，其中 1≤i≤N ， σ(i) 表示 i 的约数之和。”
现在他长大了，题目也变难了，所以麻烦你来帮他解决一道数论题吧。
他需要你求如下表达式的值：
∑ i = 1 N ∑ j = 1 N m a x ( i , j ) ⋅ σ ( i ⋅ j ) ∑^N_{i=1}∑^N_{j=1}max(i,j)⋅σ(i⋅j) i=1∑Nj=1∑Nmax(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 的值即可。
为啥我明知道它就是几个题套起来但是我就是不敢打。
A
=
∑
i
=
1
n
i
∑
j
=
1
i
σ
(
i
j
)
,
B
=
∑
i
=
1
n
i
σ
(
i
2
)
A
n
s
=
2
A
−
B
A
=
∑
i
=
1
n
i
∑
j
=
1
i
∑
a
∣
i
,
b
∣
j
a
j
[
(
a
,
b
)
=
1
]
b
=
∑
1
<
=
a
,
b
<
=
n
∑
i
=
1
n
a
i
a
∑
j
=
1
i
a
b
a
j
b
b
[
(
a
,
b
)
=
=
1
]
=
∑
1
<
=
a
,
b
<
=
n
[
(
a
,
b
)
=
1
]
a
2
∑
i
=
1
n
a
i
∑
j
=
1
i
a
b
j
发现层层嵌套太复杂，彻底打散,寻找适合化简的性质
=
∑
1
<
=
a
,
b
,
i
,
j
<
=
n
a
2
i
j
[
(
a
,
b
)
=
1
]
[
i
a
>
=
j
b
]
[
i
a
<
=
n
]
之前方法的病就在ia>=jb不好利用上，于是枚举ia和jb
=
∑
1
<
=
a
i
<
=
n
a
2
i
∑
1
<
=
j
b
<
=
a
i
j
[
(
a
,
b
)
=
1
]
实际上的变量是ia和bj，于是i和j得枚举
=
∑
1
<
=
a
i
<
=
n
a
i
∑
i
∣
a
i
i
∑
1
<
=
j
b
<
=
a
i
∑
j
∣
j
b
j
=
∑
1
<
=
a
i
<
=
n
a
i
σ
(
a
i
)
∑
1
<
=
j
b
<
=
a
i
σ
(
j
b
)
明显是积性函数，可以线性筛
\begin{aligned} &A = \sum_{i=1}^ni\sum_{j=1}^i \sigma(ij) , B = \sum_{i=1}^ni\sigma(i^2)\\ &Ans = 2A - B\\ &A=\sum_{i=1}^ni\sum_{j=1}^i\sum_{a|i,b|j}\frac {aj[(a,b)=1]}b\\ &=\sum_{1<=a,b<=n}\sum_{i=1}^{\frac na} ia\sum_{j=1}^{\frac {ia}b} \frac{ajb}b[(a,b)==1]\\ &=\sum_{1<=a,b<=n}[(a,b)=1]a^2\sum_{i=1}^{\frac na}i\sum_{j=1}^{\frac {ia}b}j\\ &\text{发现层层嵌套太复杂，彻底打散,寻找适合化简的性质}\\ &=\sum_{1<=a,b,i,j<=n}a^2ij[(a,b)=1][ia>=jb][ia<=n]\\ &\text{之前方法的病就在ia>=jb不好利用上，于是枚举ia和jb}\\ &=\sum_{1<=ai<=n}a^2i\sum_{1<=jb<=ai}j[(a,b)=1]\\ &\text{实际上的变量是ia和bj，于是i和j得枚举}\\ &=\sum_{1<=ai<=n}ai\sum_{i|ai}i\sum_{1<=jb<=ai}\sum_{j|jb}j\\ &=\sum_{1<=ai<=n}ai\sigma(ai)\sum_{1<=jb<=ai}\sigma(jb)\\ &\text{明显是积性函数，可以线性筛} \end{aligned}
A=i=1∑nij=1∑iσ(ij),B=i=1∑niσ(i2)Ans=2A−BA=i=1∑nij=1∑ia∣i,b∣j∑baj[(a,b)=1]=1<=a,b<=n∑i=1∑aniaj=1∑biabajb[(a,b)==1]=1<=a,b<=n∑[(a,b)=1]a2i=1∑anij=1∑biaj发现层层嵌套太复杂，彻底打散,寻找适合化简的性质=1<=a,b,i,j<=n∑a2ij[(a,b)=1][ia>=jb][ia<=n]之前方法的病就在ia>=jb不好利用上，于是枚举ia和jb=1<=ai<=n∑a2i1<=jb<=ai∑j[(a,b)=1]实际上的变量是ia和bj，于是i和j得枚举=1<=ai<=n∑aii∣ai∑i1<=jb<=ai∑j∣jb∑j=1<=ai<=n∑aiσ(ai)1<=jb<=ai∑σ(jb)明显是积性函数，可以线性筛
B也用差不多的方法就行，这题只需要线性筛，然后再离线一波就可以做到
O
(
n
l
o
g
n
)
O(nlogn)
O(nlogn)预处理，
O
(
1
)
O(1)
O(1)回答。
AC Code:
#include<bits/stdc++.h>
#define mod 1000000007
#define maxn 1000005
#define LL long long
using namespace std;
int mu[maxn],pr[maxn/8]
,cnt_pr,d[maxn],sumuf[maxn],sumsig[maxn],sumasig[maxn],sumasigf[maxn];
bool vis[maxn];
int lim[maxn],c[maxn];
inline bool cmp(const int &u,const int &v){ return lim[u] < lim[v]; }
int add[maxn],ret[maxn];
void sieve()
{
mu[1] = 1 , d[1] = 1;
for(int i=2;i<maxn;i++)
{
if(!vis[i]) pr[cnt_pr++] = i , mu[i] = -1 , d[i] = 1 + i;
for(int j=0;pr[j]*i<maxn;j++)
{
vis[i * pr[j]] = 1;
if(i % pr[j] == 0)
{
d[i * pr[j]] = (1ll * d[i] * (pr[j]+1) - 1ll * d[i/pr[j]] * pr[j]) % mod;
break;
}
mu[i * pr[j]] = -mu[i];
d[i * pr[j]] = 1ll * d[i] * d[pr[j]] % mod;
}
}
for(int i=1;i<maxn;i++)
sumuf[i] = (1ll * mu[i] * i * i) % mod,
sumsig[i] = (sumsig[i-1] + d[i]) % mod,
sumasig[i] = 1ll * d[i] * i % mod * (2ll * sumsig[i] - d[i]) % mod;
for(int i=1;i<maxn;i++)
for(int j=1;i*j<maxn;j++)
add[i*j]=(add[i*j]+1ll*sumuf[i]*sumasig[j]) % mod;
}
int main()
{
int T,ans=0;
scanf("%d",&T);
for(int i=0;i<T;i++) scanf("%d",&lim[i]),c[i]=i;
sort(c,c+T,cmp);
sieve();
for(int i=0,j=1;i<T;i++)
{
for(;j<=lim[c[i]];j++) ans=(ans+add[j])%mod;
ret[c[i]] = ans;
}
for(int i=0;i<T;i++) printf("Case #%d: %d\n",i+1,(ret[i]+mod)%mod);
}
（我居然都可以rank1，。。。。。。）
这个题告诉我们，杜教筛最重要的不是把式子化到有多简，而是时刻判断当前式子能否完成问题，必要时让步，用最暴力的办法写式子反而会有所发现，还有就是对于积性函数要敏感，线性筛是很强大的。。。。。。
到现在我们还可以发现一个事情：
∑
a
∣
i
,
b
∣
j
a
j
b
[
(
a
,
b
)
=
1
]
=
∑
p
∣
(
i
,
j
)
μ
(
p
)
p
∗
σ
(
i
p
)
∗
σ
(
j
p
)
\sum_{a|i,b|j}\frac {aj}b[(a,b)=1] = \sum_{p|(i,j)}\mu(p)p*\sigma(\frac ip)*\sigma(\frac jp)
a∣i,b∣j∑baj[(a,b)=1]=p∣(i,j)∑μ(p)p∗σ(pi)∗σ(pj)
ZOJ Problem Set - 3881
From ABC Conjecture
主要是变形神。
r
a
d
(
n
)
=
∏
i
=
1
k
P
i
,
w
h
e
r
e
n
=
∏
i
=
1
k
P
i
a
i
f
(
n
)
=
r
a
d
(
n
)
∗
φ
(
n
r
a
d
(
n
)
)
g
=
I
∘
f
c
a
l
c
H
(
n
)
=
∑
1
<
=
i
<
=
n
g
(
i
)
\begin{aligned} &rad(n) = \prod_{i=1}^{k}Pi,where \ n = \prod_{i=1}^kPi^{ai}\\ &f(n)=rad(n) * \varphi(\frac n{rad(n)})\\ &g = I \circ f\\ &calc \ H(n) = \sum_{1<=i<=n}g(i) \end{aligned}
rad(n)=i=1∏kPi,where n=i=1∏kPiaif(n)=rad(n)∗φ(rad(n)n)g=I∘fcalc H(n)=1<=i<=n∑g(i)
首先证明
f
(
n
)
是
积
性
函
数
f(n)是积性函数
f(n)是积性函数(显然啊)
然后就是套路了。
∑
i
∣
n
f
(
i
)
=
∑
i
∣
n
∏
j
=
1
k
f
(
P
j
a
j
)
w
h
e
r
e
∏
j
=
1
k
P
j
a
j
=
i
=
∏
j
=
1
k
∑
i
=
0
a
j
f
(
P
j
i
)
\sum_{i|n}f(i) = \sum_{i|n} \prod_{j=1}^kf(Pj^{aj}) \ where \ \prod_{j=1}^k Pj^{aj} = i\\ =\prod _{j=1}^k \sum_{i=0}^{aj}f(Pj^{i})
i∣n∑f(i)=i∣n∑j=1∏kf(Pjaj) where j=1∏kPjaj=i=j=1∏ki=0∑ajf(Pji)
这个例子告诉我们和式和积式也是可以交换枚举顺序的。
然后将
f
(
x
)
f(x)
f(x)带入就可以用高中生（小学生）的等比数列求和得到
g
(
n
)
=
∏
j
=
1
k
P
j
a
j
+
1
w
h
e
r
e
n
=
∏
j
=
1
k
P
j
a
j
g(n) = \prod_{j=1}^kPj^{aj}+1\ where \ n = \prod_{j=1}^k Pj^{aj}
g(n)=j=1∏kPjaj+1 where n=j=1∏kPjaj
这个式子我们可以高级（套路）一点：
g
(
n
)
=
∑
j
∣
n
[
(
j
,
n
j
)
=
1
]
j
g(n) = \sum_{j|n}[(j,\frac nj)=1]j
g(n)=j∣n∑[(j,jn)=1]j
出现gcd就莫比乌斯反演：
g
(
n
)
=
∑
j
∣
n
j
∑
p
∣
(
j
,
n
j
)
μ
(
p
)
=
∑
p
∣
n
μ
(
p
)
∑
j
∣
n
p
,
k
∣
n
p
[
p
j
∗
p
k
=
n
]
j
p
=
∑
p
2
∣
n
μ
(
p
)
p
∑
j
,
k
[
j
∗
k
=
n
p
2
]
j
=
∑
p
2
∣
n
μ
(
p
)
p
σ
(
n
p
2
)
g(n) = \sum_{j|n}j\sum_{p|(j,\frac nj)}\mu(p) =\sum_{p|n}\mu(p)\sum_{j|\frac np,k|\frac np}[pj * pk = n]jp\\ =\sum_{p^2|n}\mu(p)p\sum_{j,k}[j*k = \frac n{p^2}]j\\ =\sum_{p^2|n}\mu(p)p\sigma(\frac n{p^2})
g(n)=j∣n∑jp∣(j,jn)∑μ(p)=p∣n∑μ(p)j∣pn,k∣pn∑[pj∗pk=n]jp=p2∣n∑μ(p)pj,k∑[j∗k=p2n]j=p2∣n∑μ(p)pσ(p2n)
然后就随便搞了。
H
(
n
)
=
∑
i
=
1
n
∑
p
2
∣
i
μ
(
p
)
p
σ
(
i
p
2
)
=
∑
p
=
1
n
μ
(
p
)
p
∑
i
=
1
n
p
2
σ
(
i
)
=
∑
p
=
1
n
μ
(
p
)
p
S
(
n
p
2
)
H(n) = \sum_{i=1}^n\sum_{p^2|i}\mu(p)p\sigma(\frac i{p^2})\\ =\sum_{p=1}^n\mu(p)p\sum_{i=1}^{\frac n{p^2}}\sigma(i)\\ =\sum_{p=1}^{\sqrt n}\mu(p)pS(\frac n{p^2})
H(n)=i=1∑np2∣i∑μ(p)pσ(p2i)=p=1∑nμ(p)pi=1∑p2nσ(i)=p=1∑nμ(p)pS(p2n)
后面那个是约数前缀和，可以
O
(
n
)
O(\sqrt n)
O(n)
前面那个只需要枚举到
n
\sqrt n
n
复杂度可以积分估计：
∫
1
n
n
x
2
d
x
=
n
g
(
x
)
w
h
e
r
e
g
(
x
)
=
ln
x
\large\int_{1}^{\sqrt n}\sqrt \frac{n}{x^2}dx=\sqrt n g(\sqrt x) \ where \ g(x) = \ln x
∫1nx2ndx=ng(x) where g(x)=lnx
为
O
(
n
ln
n
)
O(\sqrt n \ln n)
O(nlnn)
AC Code：
#include<bits/stdc++.h>
#define maxn 1000006
#define mod 1000000009
#define LL long long
#define inv2 (mod+1)/2
using namespace std;
int mu[maxn],pr[maxn],cnt_pr;
bool vis[maxn];
void sieve()
{
mu[1] = 1;
for(int i=2;i<maxn;i++)
{
if(!vis[i]) pr[cnt_pr++] = i, mu[i] = -1;
for(int j=0;pr[j]*i<maxn;j++)
{
vis[i*pr[j]] = 1;
if(i % pr[j] == 0) break;
mu[i * pr[j]] = -mu[i];
}
}
}
int S(LL n)
{
int ans = 0;
for(LL i=1,nxt;i<=n;i=nxt+1)
{
nxt = n / (n / i);
ans = (ans + (nxt+i) % mod * ((nxt - i + 1)%mod) % mod * (n/i) % mod) % mod;
}
return (1ll * ans * inv2)%mod;
}
int main()
{
sieve();
LL N;
for(;~scanf("%lld",&N);)
{
int ans = 0;
for(int i=1;1ll*i*i<=N;i++)
ans =(ans + 1ll * mu[i] * i * S(N / (1ll * i * i))) % mod;
printf("%d\n",(ans+mod)%mod);
}
}
BZOJ 3512 DZY loves Math IV
给定n,m,
∑
i
=
1
n
∑
j
=
1
m
φ
(
i
j
)
\sum_{i=1}^n\sum_{j=1}^m\varphi (ij)
∑i=1n∑j=1mφ(ij)求 模10^9+7的值。
n
<
=
1
0
5
,
m
<
=
1
0
9
n<=10^5,m<=10^9
n<=105,m<=109
这个
充分告诉我们
φ
\varphi
φ的性质可以创造奇迹。
首先有
φ
(
i
j
)
=
φ
(
i
)
φ
(
j
)
d
φ
(
d
)
\varphi(ij) = \frac {\varphi(i)\varphi(j)d}{\varphi(d)}
φ(ij)=φ(d)φ(i)φ(j)d
然后发现做不下去。。。。。。
发现
n
<
=
1
0
5
n<=10^5
n<=105
可以考虑枚举一维x
接下来就只需要求
f
(
x
,
m
)
=
∑
i
=
1
m
φ
(
x
i
)
f(x,m) = \sum_{i=1}^m\varphi(xi)
f(x,m)=∑i=1mφ(xi)
但是发现好像还是不好求。
考虑欧拉函数的性质，如果
x
x
x是无平方因子数。
f
(
x
,
m
)
=
∑
i
=
1
m
φ
(
x
)
φ
(
i
)
(
i
,
x
)
φ
(
(
i
,
x
)
)
=
∑
i
=
1
m
φ
(
x
(
i
,
x
)
)
φ
(
i
)
(
i
,
x
)
f(x,m)=\sum_{i=1}^m\frac {\varphi(x)\varphi(i)(i,x)}{\varphi((i,x))}=\sum_{i=1}^m\varphi(\frac x{(i,x)})\varphi(i)(i,x)
f(x,m)=i=1∑mφ((i,x))φ(x)φ(i)(i,x)=i=1∑mφ((i,x)x)φ(i)(i,x)
然后就像用
∑
p
∣
g
c
d
μ
(
p
)
=
[
g
c
d
=
=
1
]
\sum_{p|gcd}\mu(p) = [gcd==1]
p∣gcd∑μ(p)=[gcd==1]一样用
∑
p
∣
g
c
d
φ
(
p
)
=
g
c
d
\sum_{p|gcd}\varphi(p) = gcd
p∣gcd∑φ(p)=gcd
来将和gcd有关的式子化为和式然后再调换求和顺序最终转化为枚举约数：
f
(
x
,
m
)
=
∑
i
=
1
m
φ
(
x
d
)
φ
(
i
)
∑
p
∣
d
φ
(
p
)
f(x,m)=\sum_{i=1}^m\varphi(\frac xd)\varphi(i)\sum_{p|d}\varphi(p)\\
f(x,m)=i=1∑mφ(dx)φ(i)p∣d∑φ(p)
貌似已经山穷水尽了？毕竟还有一个gcd呢。
发现
p
p
p和
x
d
\frac xd
dx是互质的然后就很简单了。
f
(
x
,
m
)
=
∑
i
=
1
m
φ
(
i
)
∑
p
∣
d
φ
(
x
p
d
)
=
∑
i
=
1
m
φ
(
i
)
∑
p
∣
d
φ
(
x
p
)
!
!
!
=
∑
p
∣
x
φ
(
x
p
)
∑
i
=
1
m
p
φ
(
i
p
)
=
∑
p
∣
x
φ
(
x
p
)
f
(
p
,
m
p
)
f(x,m)=\sum_{i=1}^m\varphi(i)\sum_{p|d} \varphi(\frac {xp}d)\\ =\sum_{i=1}^m\varphi(i)\sum_{p|d}\varphi(\frac xp)!!!\\ =\sum_{p|x}\varphi(\frac xp)\sum_{i=1}^{\frac mp}\varphi(ip) =\sum_{p|x}\varphi(\frac xp)f(p,\frac mp)
f(x,m)=i=1∑mφ(i)p∣d∑φ(dxp)=i=1∑mφ(i)p∣d∑φ(px)!!!=p∣x∑φ(px)i=1∑pmφ(ip)=p∣x∑φ(px)f(p,pm)
然后对于非无平方因子数，把多了的质因子提出来就行了。
这里面的
φ
\varphi
φ技巧好多啊。
1:
∑
p
∣
g
c
d
φ
(
p
)
=
g
c
d
\sum_{p|gcd}\varphi(p) = gcd
p∣gcd∑φ(p)=gcd
2:积性函数性质
3:提多余质因子，提成无平方因子数
4:
φ
(
x
y
)
=
φ
(
x
)
φ
(
y
)
(
x
,
y
)
φ
(
(
x
,
y
)
)
\varphi(xy)=\frac {\varphi(x)\varphi(y)(x,y)}{\varphi((x,y))}
φ(xy)=φ((x,y))φ(x)φ(y)(x,y)
然后。。。。。。
这差不多就是
φ
\varphi
φ的套路吧。
AC Code：
#include<bits/stdc++.h>
#define maxn 100005
#define LL long long
#define mod 1000000007
using namespace std;
vector<int>di[maxn];
#define ts 5000005
int phi[ts],sphi[ts],pr[ts/8],cnt_pr,isd[ts];
bool vis[ts];
void sieve()
{
phi[1] = 1 , isd[1] = 1 , sphi[1] = 1;
for(int i=2;i<ts;i++)
{
if(!vis[i]) pr[cnt_pr++] = i , phi[i] = i-1 , isd[i] = i;
for(int j=0;pr[j] * i < ts ; j++)
{
vis[i * pr[j]] = 1;
if(i % pr[j] == 0)
{
phi[i * pr[j]] = phi[i] * pr[j];
isd[i * pr[j]] = isd[i];
break;
}
isd[i * pr[j]] = isd[i] * pr[j];
phi[i * pr[j]] = phi[i] * pr[j]- phi[i];
}
sphi[i] = (phi[i] + sphi[i-1]) % mod;
}
for(int i=1;i<maxn;i++)
for(int j=i;j<maxn;j+=i)
di[j].push_back(i);
}
map<LL,int>mp,st[maxn];
int Phi(LL n)
{
if(n < ts) return sphi[n];
if(mp.count(n)) return mp[n];
static int inv2 = (mod+1) / 2;
int tmp = 1ll * (n % mod) * ((n+1) % mod) % mod * inv2 % mod;
for(LL i=2,nxt;i<=n;i=nxt+1)
{
nxt = n / (n / i);
tmp = (tmp - 1ll * (nxt - i + 1) % mod * Phi(n / i)) % mod;
}
return mp[n] = tmp;
}
LL n , m;
int S(LL n,LL m)
{
if(n == 1) return Phi(m);
if(st[n].count(m)) return st[n][m];
int w = isd[n] , y = n / w;
int tmp = 0;
for(int i=0,siz=di[w].size();i<siz && di[w][i]<=m;i++)
tmp =(tmp + 1ll * phi[w/di[w][i]] * S(di[w][i],m/di[w][i])) % mod;
return st[n][m] = 1ll * tmp * y % mod;
}
int main()
{
sieve();
scanf("%lld%lld",&n,&m);
int ans = 0;
for(int i=1;i<=n;i++)
ans = (ans + S(i , m)) % mod;
printf("%d\n",(ans+mod)%mod);
}
ZOJ 5340 - The Sum of Unitary Totient（常规积性函数求和，数据组数较多，只可分段打表）
SPOJ DIVCNT3 - Counting Divisors (cube)（常规积性函数求和，注意代码长度限制，不可打表）
51Nod 1575 - Gcd and Lcm（常规积性函数求和，可分段打表）
51Nod 1847 - 奇怪的数学题（非常规积性函数求和，综合题，可分段打表）
最后这4个实在是找不到杜教筛的题解。
准备转Min_25/洲阁筛了，撑不住了。。。。。。。
大家请不欢而散吧。
1180
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
