1：完整的Min25筛学习方案
2：虽然不是但我也挂
3：网上捞了好久才捞到的纯递推写法Min25筛
终于看懂了。。。。。
就是充分利用合数一定有小于等于
n
\sqrt n
n的质因子和积性函数的性质来筛.
打了一个欧拉函数的板。。。。。。
需要注意的是
i
2
i^2
i2并不能用来判断质数。。。。。。
因为
m
o
d
mod
mod过了，只能用
i
i
i
调了
1
h
1h
1h。。。。
555
模板：51nod 1222 最小公倍数计数
Code：
#include<bits/stdc++.h>
#define maxn 1000006
#define LL long long
using namespace std;
LL n,s[maxn],a[maxn];
int id,sn,pr[maxn],cnt_pr;
inline int ID(LL a){ return a<=sn?a:id-n/a+1; }
LL Sum(LL n,int b){
if(n<pr[b]) return 0;
LL r=s[ID(n)]-s[pr[b-1]];
for(int i=b;i<=cnt_pr && 1ll * pr[i] * pr[i] <= n;i++)
for(LL x=pr[i],f=3;x*pr[i]<=n;x*=pr[i],f+=2)
r+=f*Sum(n/x,i+1)+f+2;
return r;
}
LL solve(LL n){
::n=n;sn=sqrt(n);id=cnt_pr=0;
for(LL i=1;i<=n;i++) a[++id]=(i=n/(n/i)),s[id]=3*i-3;
for(LL i=1;i<=sn;i++) if(s[i]^s[i-1]) for(LL j=id,LIM=(pr[++cnt_pr]=i)*i,t0=s[i-1];a[j]>=LIM;j--)
s[j]-=s[ID(a[j]/i)]-t0;
return Sum(n,1)+(n>=1);
}
int main(){
LL a,b;
scanf("%lld%lld",&a,&b);
printf("%lld",(solve(b)-solve(a-1)+b-a+1)/2);
}
51 nod Gcd and Lcm
求
∑
i
=
1
n
∑
j
=
1
i
∑
k
=
1
i
[
(
i
,
j
)
,
(
i
,
k
)
]
\sum_{i=1}^n\sum_{j=1}^i\sum_{k=1}^i[(i,j),(i,k)]
∑i=1n∑j=1i∑k=1i[(i,j),(i,k)]
对于Min_25筛就是裸题.
f
(
x
)
=
∑
j
=
1
x
∑
k
=
1
x
[
(
x
,
j
)
,
(
x
,
k
)
]
f(x) = \sum_{j=1}^x\sum_{k=1}^x[(x,j),(x,k)]
f(x)=∑j=1x∑k=1x[(x,j),(x,k)]
f
(
p
)
=
(
p
−
1
)
2
+
p
(
2
p
−
1
)
=
3
p
2
−
3
p
+
1
f(p)=(p-1)^2+p(2p-1)=3p^2-3p+1
f(p)=(p−1)2+p(2p−1)=3p2−3p+1
f
(
p
k
)
=
(
2
k
+
1
)
(
p
2
k
−
p
2
k
−
1
)
+
p
k
−
1
f(p^k)=(2k+1)(p^{2k}-p^{2k-1})+p^{k-1}
f(pk)=(2k+1)(p2k−p2k−1)+pk−1
UPD:附上凌乱不堪的草稿过程：
\sum_{j=1}^n \sum_{k=1}^n
[(i,j),(i,k)] =
for p itis (n-1)*(n-1)+(2*n-1)*n
= n^2 - 2n + 1 + 2n^2 - n
= 3n^2 - 3n + 1
for p^q itis
gcd(i,p^q) = p^k holds p^{q-k}-p^{q-k-1}
lcm(gcd(i,p^q),gcd(j,p^q)) = p^k holds
(\sum_{g=0}^{k-1} p^{q-g}-p^{q-g-1}) * 2 * (p^{q-k}-p^{q-k-1}) + (p^{q-k}-p^{q-k-1}) * (p^{q-k}-p^{q-k-1})
= ((p^q - p^{q-k}) * 2 + p^{q-k} - p^{q-k-1}) * (p^{q-k}-p^{q-k-1})
= (2*p^q-p^{q-k}-p^{q-k-1}) * (p^{q-k}-p^{q-k-1})
ans for p^k :
(2*p^q-p^{q-k}-p^{q-k-1}) * (p^{q-k}-p^{q-k-1}) * p^k
= (2*p^q - p^{q-k}-p^{q-k-1}) * (p^q - p^{q-1})
ans = (p^q - p^{q-1}) * (2 * q * p^q - \sum_{k=0}^{q-1} p^{q-k}+p^{q-k-1}) + (2*p^q-1)*p^q
= (p^q - p^{q-1}) * (2 * q * p^q - \frac {p^q + p^{q+1}}{p-1} + 1) + 2*p^{2q}-p^q
= (p^q-p^{q-1})/(p-1) * (2 * q * p^q * (p-1) - p^q(p+1) + p-1) + 2*p^{2q}-p^q
= p^{2q-1} * (2 * q * p - 2 * q + p - 1) + p^{q-1}
= p^{2q-1} * (2*q+1)(p-1) + p^{q-1}
= (p^{2q} - p^{2q-1}) * (2*q+1) + p^{q-1}
AC Code:
#include<bits/stdc++.h>
#define maxn 100005
#define LL unsigned int
using namespace std;
LL n;
LL pr[maxn],a[maxn],g[maxn],m2[maxn],m1[maxn],m0[maxn],sm2[maxn],sm1[maxn],sm0[maxn],sg[maxn];
int sn,id,cnt_pr;
inline int ID(LL a){ return a <= sn ? a : id - n/a + 1; }
inline LL calc3(LL n)
{
LL a[3] = {n,n+1,2*n+1};
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]*a[1]*a[2];
}
inline LL calc2(LL n)
{
LL a[2] = {n , n+1};
for(int i=0;i<2;i++) if(a[i]%2==0) {a[i]/=2;break;}
return a[0] * a[1];
}
void sieve()
{
sn = sqrt(n);
id = cnt_pr = 0;
for(LL i=1;i<=n;i++)
a[++id] = i = (n / (n / i)),
m2[id] = calc3(i) - 1,
m1[id] = calc2(i) - 1,
m0[id] = i - 1;
for(LL i=1;i<=sn;i++)
if(m0[i]!=m0[i-1])
{
pr[++cnt_pr] = i ,
sm2[cnt_pr] = sm2[cnt_pr-1] + i * i , sm1[cnt_pr] = sm1[cnt_pr-1] + i , sm0[cnt_pr] = sm0[cnt_pr-1] + 1;
for(LL j=id,lim=i*i;a[j]>=lim;j--)
m2[j] -= i*i*(m2[ID(a[j]/i)]-sm2[cnt_pr-1]),
m1[j] -= i * (m1[ID(a[j]/i)]-sm1[cnt_pr-1]),
m0[j] -= m0[ID(a[j]/i)] - sm0[cnt_pr-1];
}
for(int i=1;i<=id;i++)
g[i] = 3 * m2[i] - 3 * m1[i] + m0[i],
sg[i] = 3 * sm2[i] - 3 * sm1[i] + sm0[i];
}
LL S(LL a,LL b)
{
if(a < pr[b]) return 0;
LL ret = g[ID(a)] - sg[b-1];
for(int i=b;i<=cnt_pr && pr[i]*pr[i]<=a;i++)
for(LL x=pr[i],f1=(pr[i]*pr[i]-pr[i]),f2=1,lim=a/pr[i],k=3;x<=lim;x*=pr[i])
ret += S(a/x,i+1) * (f1*k + f2),
ret += (f1 = f1 * pr[i] * pr[i]) * (k=k+2) + (f2 = f2 * pr[i]);
return ret;
}
int main()
{
int T;
for(scanf("%d",&T);T--;)
{
cin>>n;
sieve();
cout<<S(n,1)+1<<endl;
}
}
可以发现Min_25筛挺套路的.
UPD2:更短的实现：
#include<bits/stdc++.h>
#define LL unsigned
#define inv3 2863311531
#define maxn 100005
using namespace std;
LL n,s[3][maxn],a[maxn];
int id,cnt_pr,sn,pr[maxn];
int ID(int a){return a<=sn?a:id-n/a+1; }
LL Sum(int n,int b){
if(n<pr[b]) return 0;
LL p=ID(n),r=(3*s[2][p]-3*s[1][p]+s[0][p])-(3*s[2][pr[b-1]]-3*s[1][pr[b-1]]+s[0][pr[b-1]]);
for(int i=b;i<=cnt_pr&&pr[i]*pr[i]<=n;i++) for(LL x=pr[i],f=x*x-x,c=3,q=1;1ll*x*pr[i]<=n;x*=pr[i],f*=pr[i]*pr[i],q*=pr[i],c+=2)
r+=(f*c+q)*Sum(n/x,i+1)+(f*pr[i]*pr[i])*(c+2)+q*pr[i];
return r;
}
int main(){
int T;
for(scanf("%d",&T);T--;){
scanf("%u",&n);sn=sqrt(n),id=cnt_pr=0;
for(int i=1;i<=n;i++) a[++id]=(i=n/(n/i)),s[0][id]=i-1,s[1][id]=i*(i+1ll)/2-1,s[2][id]=i*(i+1ll)/2*(2*i+1)*inv3-1;
for(int i=1;i<=sn;i++) if(s[0][i]^s[0][i-1]) for(int j=id,LIM=(pr[++cnt_pr]=i)*i;a[j]>=LIM;j--)
s[0][j]-=s[0][ID(a[j]/i)]-s[0][i-1],
s[1][j]-=(s[1][ID(a[j]/i)]-s[1][i-1])*i,
s[2][j]-=(s[2][ID(a[j]/i)]-s[2][i-1])*i*i;
printf("%u\n",Sum(n,1)+(n>=1));
}
}
51nod 1847 奇怪的数学题
做不来
ZOJ Problem Set - 3808 The Sum of Unitary Totient
ZOJ 3808 我的题解
怎么现在觉着这这么套路…
51nod 1847 奇怪的数学题
终于做完唐老师博客题了（什么
C
o
u
n
t
i
n
g
D
i
v
i
s
o
r
s
k
Counting \ Divisors \ k
Counting Divisors k就忽略吧）
给出 N,K ,请计算下面这个式子:
∑ i = 1 N ∑ j = 1 N s g c d ( i , j ) k \sum_{i=1}^N\sum_{j=1}^Nsgcd(i,j)^k ∑i=1N∑j=1Nsgcd(i,j)k
其中,sgcd(i, j)表示(i, j)的所有公约数中第二大的,特殊地,如果gcd(i, j) = 1, 那么sgcd(i, j) = 0。
考虑到答案太大,请输出答案对2^32取模的结果.
1 ≤ N ≤ 1 0 9 , 1 ≤ K ≤ 50 1≤N≤10^9,1≤K≤50 1≤N≤109,1≤K≤50
样例解释:
因为gcd(i, j)=1时sgcd(i,j)=0对答案没有贡献,所以我们只考虑gcd(i,j)>1的情况.
当i是2时,j是2时,sgcd(i,j)=1,它的K次方是1
当i是2时,j是4时,sgcd(i,j)=1,它的K次方是1
当i是3时,j是3时,sgcd(i,j)=1,它的K次方是1
当i是4时,j是2时,sgcd(i,j)=1,它的K次方是1
当i是4时,j是4时,sgcd(i,j)=2,它的K次方是8
当i是5时,j是5时,sgcd(i,j)=1,它的K次方是1
这个，首先可以发现
s
g
c
d
=
g
c
d
/
m
i
n
d
sgcd = gcd / mind
sgcd=gcd/mind,
m
i
n
d
为
g
c
d
的
最
小
质
因
数
mind为gcd的最小质因数
mind为gcd的最小质因数
然后就枚举gcd有：
A
N
S
=
∑
d
=
1
n
(
d
m
i
n
d
)
k
∑
i
=
1
n
d
∑
j
=
1
n
d
[
g
c
d
(
i
,
j
)
=
1
]
ANS = \sum_{d=1}^n \left( \frac d{mind} \right)^k\sum_{i=1}^{\frac nd}\sum_{j=1}^{\frac nd}[gcd(i,j)=1]
ANS=d=1∑n(mindd)ki=1∑dnj=1∑dn[gcd(i,j)=1]
然后发现后面是一个很套路的函数，化成
μ
\mu
μ或
ϕ
\phi
ϕ都行，不过
ϕ
\phi
ϕ更裸。
前面不是一个积性函数，但是我们可以用
M
i
n
2
5
Min_25
Min25筛的思想解决他。
质数的贡献为1，合数都可以被小于等于
n
\sqrt n
n的质数筛出来，这个还是很好筛的，就是k次方，我用的是第二类斯特林数。
完了。。。。。。
#include<bits/stdc++.h>
#define maxn 100005
#define LL unsigned int
using namespace std;
LL n,k,str2[51][51],pr[maxn],a[maxn],g[maxn],h[maxn],sg[maxn],sh[maxn],phi[maxn],sp[maxn];
LL sn,id,cnt_pr;
inline LL ID(LL a){ return a <= sn ? a : id - n / a + 1; }
inline LL calc(LL n)
{
LL ret = 0;
for(LL i=1;i<=k;i++)
{
LL C = str2[k][i];
for(LL j=0;j<=i;j++)
if((n-j+1) % (i+1) == 0)
C *= (n-j+1) / (i+1);
else
C *= (n-j+1);
ret += C;
}
return ret;
}
inline LL Pow(LL base,LL k)
{
LL ret = 1;
for(;k;k>>=1,base=base*base) if(k&1) ret=ret*base;
return ret;
}
LL f[70000][70];
LL F(LL a,LL b)
{
if(a < pr[b]) return 0;
int tmp = ID(a);
if(tmp < 70000 && b < 70 && f[tmp][b]!=0)
return f[tmp][b];
LL ret = (g[tmp]) - sg[b-1];
for(LL i=b;i<=cnt_pr && pr[i] * pr[i] <= a;i++)
for(LL pw=sg[i]-sg[i-1],f=pw,x=pr[i],lim=a/pr[i];x>=1 && x<=lim;x*=pr[i])
ret += f * F(a/x,i+1),
ret += (f = f * pw);
if(tmp < 70000 && b < 70)
f[tmp][b] = ret;
return ret;
}
LL G(LL a)
{
LL ret = 0;
for(LL i=1;i<=cnt_pr && pr[i] * pr[i] <= a;i++)
for(LL f=1,x=pr[i],lim=a/pr[i],pk=sg[i]-sg[i-1];x>=1 && x<=lim;x*=pr[i],f*=pk)
ret += f * F(a / x,i+1) + f * pk;
return ret;
}
LL ph[70000][70];
LL P(LL a,LL b)
{
if(a < pr[b]) return 0;
int tmp = ID(a);
if(tmp < 70000 && b < 70 && ph[tmp][b]!=0)
return ph[tmp][b];
LL ret = (phi[tmp] - h[ID(a)]) - (sp[b-1] - sh[b-1]);
for(LL i=b;i<=cnt_pr && pr[i] * pr[i] <= a;i++)
for(LL x = pr[i],f=pr[i]-1,lim = a/pr[i];x<=lim;x*=pr[i],f*=pr[i])
ret += f * P(a / x , i + 1) + f * pr[i];
if(tmp < 70000 && b < 70)
ph[tmp][b] = ret;
return ret;
}
int main()
{
cin>>n>>k;
str2[0][0] = 1;
for(LL i=1;i<=k;i++)
for(LL j=1;j<=i;j++)
str2[i][j] = j * str2[i-1][j] + str2[i-1][j-1];
for(LL i=1;i<=n;i++)
a[++id] = i = (n/(n/i)),
h[id] = i - 1,
g[id] = calc(i) - 1,
phi[id] = 1ll * i * (i+1) / 2 - 1;
sn = sqrt(n);
for(LL i=1;i<=sn;i++)
if(h[i] != h[i-1])
{
LL pw = Pow(i,k);
pr[++cnt_pr] = i , sg[cnt_pr] = sg[cnt_pr-1] + pw , sh[cnt_pr] = sh[cnt_pr-1] + 1;
sp[cnt_pr] = sp[cnt_pr-1] + i;
for(LL j=id,lim=i*i;a[j]>=lim;j--)
{
g[j] -= pw * (g[ID(a[j] / i)] - sg[cnt_pr-1]);
h[j] -= (h[ID(a[j] / i)] - sh[cnt_pr-1]);
phi[j] -= i * (phi[ID(a[j]/i)] - sp[cnt_pr-1]);
}
}
LL ans = 0;
for(LL i=1,nxt,pre=0,dpre;i<=n;i=nxt+1)
{
nxt = n / (n / i);
ans += ((dpre = ((G(nxt)) + h[ID(nxt)])) - pre) * (2 * (P(n / i,1)+1) - 1);
pre = dpre;
}
cout<<ans<<endl;
}
递推版：
#include<bits/stdc++.h>
#define maxn 100005
#define un unsigned
#define LL long long
using namespace std;
un s[3][maxn],S[55][55],g[maxn],h[maxn],phi[maxn];
int n,K,a[maxn],sn,id,prk[maxn];
int ID(int a){ return a<=sn?a:id-n/a+1; }
un Pow(un b,un k){ un r=1;for(;k;k>>=1,b=b*b) if(k&1) r=r*b;return r; }
un cal(int n){
un r = 0;
for(int i=1;i<=K;i++){
un p = n+1 , s=1;
for(;p%(i+1);p--) s*=p;
s*=p--/(i+1);
for(;n-p<i;p--) s*=p;
r += s * S[K][i];
}
return r;
}
int main(){
scanf("%d%d",&n,&K);S[0][0]=1;
for(int i=1;i<=K;i++) for(int j=1;j<=i;j++) S[i][j]=S[i-1][j]*j+S[i-1][j-1];
sn = sqrt(n);
for(int i=1;i<=n;i++) a[++id]=(i=(n/(n/i))),s[0][id]=i-1,s[1][id]=i*(i+1ll)/2-1,s[2][id]=cal(i)-1;
for(int i=1;i<=sn;i++) if(s[0][i]^s[0][i-1])
for(int j=id,LIM=(prk[i]=Pow(i,K),i)*i;a[j]>=LIM;j--)
s[0][j]-=s[0][ID(a[j]/i)]-s[0][i-1],
s[1][j]-=i*(s[1][ID(a[j]/i)]-s[1][i-1]),
s[2][j]-=prk[i]*(s[2][ID(a[j]/i)]-s[2][i-1]);
for(int i=1;i<=id;i++) g[i]=s[2][i],h[i]=s[0][i],phi[i]=s[1][i]-s[0][i];
for(int i=sn;i>1;i--) if(s[0][i]^s[0][i-1]) for(int j=id,LIM=i*i;a[j]>=LIM;j--)
for(LL x=i,f=1,c=i-1;x*i<=a[j];x*=i,f*=prk[i],c*=i)
g[j]+=(g[ID(a[j]/x)]-s[2][i]+prk[i])*f*prk[i],
h[j]+=(g[ID(a[j]/x)]-s[2][i]+prk[i])*f,
phi[j]+=(phi[ID(a[j]/x)]-s[1][i]+s[0][i]+i)*c;
un ans = 0;
for(int i=2,nxt;i<=n;i=nxt+1){
nxt=n/(n/i);
ans += (h[ID(nxt)] - h[ID(i-1)]) * (2 * phi[ID(n/i)] + 1);
}
printf("%u\n",ans);
}
唐老师博客题AC全达成！
2020/1/13UPD:
写51nod1238想用
M
i
n
25
Min25
Min25筛，然后发现递归版非常之慢，怕不是不支持对于
f
(
n
1
)
,
f
(
n
2
)
.
.
.
f
(
n
k
)
f(\frac n1),f(\frac n2)...f(\frac nk)
f(1n),f(2n)...f(kn)的多组询问，写了一个
m
a
p
map
map记忆化发现能过
1
e
8
1e8
1e8过不了
1
e
9
1e9
1e9。。。。。。怕不是状态数问题，然后疯狂上网搜资料，发现我学了个假的
M
i
n
25
Min25
Min25筛。
首先
M
i
n
25
Min_{25}
Min25筛的复杂度证明实际上是这样的：
对于
n
<
=
1
e
13
n<=1e13
n<=1e13我们可以通过一些放缩（打表）技巧严谨得到复杂度是
O
(
n
3
4
log
n
)
O(\frac {n^\frac 34}{\log n})
O(lognn43)的。
但是对于
n
>
1
e
13
n>1e13
n>1e13我们放缩的条件中有一个失效了，实际上应该是
O
(
n
1
−
ϵ
)
=
O
(
n
log
log
n
)
O(n^{1-\epsilon })=O(\frac n{\log\log n})
O(n1−ϵ)=O(loglognn)的。
所以说
M
i
n
25
Min25
Min25筛严格来说是一个亚线性算法，并不是
n
3
4
n^{\frac 34}
n43之类的优秀复杂度，但是足够优秀的常数保证了在CCF的机子换成超算之前我们都是可以用它碾压洲阁筛的。
所以我们还是可以心安理得的在我们的有生之年使用
M
i
n
25
Min25
Min25筛法，但是Min26何时普及？
所以我们的状态数得到了保证，但是 m a p map map的一个 log \log log无法接受，所以我们可以用递推的写法求出所有点的值。
参考代码：
#include<bits/stdc++.h>
#define maxn 400005
#define LL long long
#define mod 1000000007
#define inv6 166666668
using namespace std;
LL n,a[maxn];
int f[3][maxn],g[maxn],pr[maxn],cnt_pr,sn,id;
map<LL,int>mAp[maxn/8];
int ID(LL a){ return a<=sn?a:id-n/a+1; }
int sqr(int a){ return 1ll*a*a%mod; }
int S22(int n){ return sqr(n*(n+1ll)/2%mod); }
int S2(int n){ return n*(n+1ll)/2%mod; }
int S3(int n){ return n*(n+1ll)%mod*(2*n+1)%mod*inv6%mod; }
int main(){
scanf("%lld",&n);sn=sqrt(n);
for(LL i=1;i<=n;i++) a[++id]=(i=(n/(n/i))),
f[0][id]=i%mod-1,
f[1][id]=S2(i%mod)-1,
f[2][id]=S3(i%mod)-1;
for(LL i=1;i<=sn;i++) if(f[0][i]^f[0][i-1]) for(LL j=id,LIM=i*i,t0=f[0][i-1],t1=f[1][i-1],t2=f[2][i-1];a[j]>=LIM;j--)
f[0][j]=(f[0][j]-f[0][ID(a[j]/i)]+t0)%mod,
f[1][j]=(f[1][j]-(f[1][ID(a[j]/i)]-t1)*i)%mod,
f[2][j]=(f[2][j]-(f[2][ID(a[j]/i)]-t2)*i%mod*i)%mod;
for(int i=1;i<=id;i++) g[i]=(f[1][i]-f[2][i])%mod;
for(int i=sn;i>=1;i--) if(f[0][i]^f[0][i-1]) for(LL j=id,LIM=1ll*i*i,C=i-f[1][i]+f[2][i];a[j]>=LIM;j--)
for(LL x=i,f=i*(1ll-i)%mod;x*i<=a[j];x*=i,f=f*i%mod)
g[j]=(g[j]+f*(g[ID(a[j]/x)]+C))%mod;
int ans = 0;
for(LL i=1,nxt,p=0,rp;i<=n;i=nxt+1){
nxt=n/(n/i);
ans = (ans + 1ll * S22(n/i%mod) * (g[ID(nxt)]-g[ID(i-1)]+(i==1)))%mod;
}
printf("%d\n",(ans+mod)%mod);
}
20200804UPD:被教育了：“Min25筛是递归写法，递推的叫洲阁筛”。
UPD2:
powerful number更好的教程。。。。。。
51nod 2026 Gcd and Lcm
对于
f
(
x
)
=
∑
d
∣
x
x
μ
(
x
)
f(x)=\sum_{d|x}x\mu(x)
f(x)=∑d∣xxμ(x)这个函数。
容易发现它的值只和其质因数有关。
f
(
x
)
=
∏
p
r
∣
x
,
p
r
为
质
数
(
1
−
p
r
)
f(x)=\prod_{pr|x,pr为质数} (1-pr)
f(x)=∏pr∣x,pr为质数(1−pr)
所以
f
(
g
c
d
(
i
,
j
)
)
∗
f
(
l
c
m
(
i
,
j
)
)
=
f
(
i
)
∗
f
(
j
)
f(gcd(i,j)) * f(lcm(i,j)) = f(i) * f(j)
f(gcd(i,j))∗f(lcm(i,j))=f(i)∗f(j)
两者的和=两个的并+两者的交。。。。
然后就是求
f
f
f的前缀和即可。
这个用
M
i
n
25
Min_{25}
Min25光速完成。
A
C
C
o
d
e
\mathrm {AC \ Code}
AC Code
#include<bits/stdc++.h>
#define mod 1000000007
#define maxn 100005
#define LL long long
using namespace std;
int n,s[2][maxn],a[maxn];
int id,sn,pr[maxn],cnt_pr;
inline int ID(int a){ return a<=sn?a:id-n/a+1; }
int sqr(int a){ return 1ll * a * a % mod; }
int Sum(int n,int b){
if(n<pr[b]) return 0;
int r=(1ll*(s[0][ID(n)]-s[1][ID(n)])-(s[0][pr[b-1]]-s[1][pr[b-1]]))%mod;
for(int i=b;i<=cnt_pr && 1ll*pr[i]*pr[i]<=n;i++) for(LL x=pr[i],f=1-pr[i];x*pr[i]<=n;x*=pr[i])
r=(r+f*Sum(n/x,i+1)+f)%mod;
return r;
}
int main(){
scanf("%d",&n);sn=sqrt(n);
for(int i=1;i<=n;i++) a[++id]=(i=(n/(n/i))),s[0][id]=i-1,s[1][id]=i*(i+1ll)/2%mod-1;
for(int i=1;i<=sn;i++) if(s[0][i]^s[0][i-1]) for(int j=id,LIM=(pr[++cnt_pr]=i)*i,t0=s[0][i-1],t1=s[1][i-1];a[j]>=LIM;j--)
s[0][j]=(s[0][j]-(s[0][ID(a[j]/i)]-t0))%mod,
s[1][j]=(s[1][j]-i*1ll*(s[1][ID(a[j]/i)]-t1))%mod;
printf("%d\n",(sqr(Sum(n,1)+1)+mod)%mod);
}
HDU 5608 function
这个题的式子虽然明示杜教筛，但是我们还是可以用
M
i
n
25
Min25
Min25筛去做。
另外需要注意
f
f
f不是积性函数。
∑
i
=
n
3
n
n
z
\sum_{i=n} \frac {3n}{n^z}
∑i=nnz3n不是积性函数。
要想用
M
i
n
25
Min25
Min25就必须把这个函数拆成3部分并且还要处理常数。
但是预处理比杜教筛短多了。
A
C
C
o
d
e
\mathrm {AC Code}
ACCode
#include<bits/stdc++.h>
#define mod 1000000007
#define maxn 100005
#define LL long long
#define inv6 166666668
using namespace std;
int n,s[3][maxn],a[maxn];
int id,sn,pr[maxn],cnt_pr;
inline int ID(int a){ return a<=sn?a:id-n/a+1; }
int sqr(int a){ return 1ll * a * a % mod; }
int cal(int a,int t){
if(t == 1) return a;
if(t == 2) return a * 1ll * a % mod;
}
int Sum(int n,int b,int t){
if(n<pr[b]) return 0;
int r = (s[t][ID(n)] - 1ll * s[0][ID(n)] - 1ll * s[t][pr[b-1]] + s[0][pr[b-1]])%mod;
for(int i=b;i<=cnt_pr && 1ll*pr[i]*pr[i]<=n;i++) for(LL x=pr[i],f=cal(x,t)-cal(1,t);x*pr[i]<=n;f=(cal(x*pr[i],t)-cal(x,t))%mod,x*=pr[i])
r=(r+f*Sum(n/x,i+1,t)+cal(x*pr[i],t)-cal(x,t))%mod;
return r;
}
int main(){
int T;
for(scanf("%d",&T);T--;){
scanf("%d",&n);sn=sqrt(n);id=cnt_pr=0;
for(int i=1;i<=n;i++) a[++id]=(i=(n/(n/i))),s[0][id]=i-1,s[1][id]=i*(i+1ll)/2%mod-1,s[2][id]=i*(i+1ll)%mod*(2*i+1ll)%mod*inv6%mod-1;
for(int i=1;i<=sn;i++) if(s[0][i]^s[0][i-1]) for(int j=id,LIM=(pr[++cnt_pr]=i)*i,t0=s[0][i-1],t1=s[1][i-1],t2=s[2][i-1];a[j]>=LIM;j--)
s[0][j]=(s[0][j]-(1ll*s[0][ID(a[j]/i)]-t0))%mod,
s[1][j]=(s[1][j]-i*1ll*(s[1][ID(a[j]/i)]-t1))%mod,
s[2][j]=(s[2][j]-i*1ll*i%mod*(s[2][ID(a[j]/i)]-t2))%mod;
printf("%lld\n",((3-3ll*(Sum(n,1,1)+1)+Sum(n,1,2))%mod+mod)%mod);
}
}
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
