原文链接https://www.cnblogs.com/zhouzhendong/p/Min-25.html
前置技能
埃氏筛法
整除分块(这里有提到)
本文概要
1. 问题模型
2. Min_25 筛
3. 模板题以及模板代码
问题模型
有一个积性函数 $f$ ,对于所有质数 $p$,$f(p)$ 是关于 $p$ 的多项式,$f(p^k)$ 非常容易计算(不一定是关于 p 的多项式)。
求
$$sum_{i=1}^{n} f(i)$$
$nleq 10^{10}$
${ m Time Limit} = 1s$
Min_25 筛
设集合 $P$ 表示素数集合。
设
$$g_{n,m} = sum_{2leq ileq n, forall pin P and pleq m,p mid i} f(i)$$
则假设 $pin P$。
$$g(n,m) = sum_{m<pleq sqrt n ,p^eleq n,egeq 1} f(p^e) left([e>1] + sum_{2leq x leq lfloor frac n {p^e} floor, forall p'in P and p'leq p ,p' mid x}f(x) ight)+sum_{m<pleq n} f(p)$$
设
$$h(n) = sum_{1leq pleq n} f(p)$$
则
$$g(n,m)=sum_{m<pleq sqrt n ,p^eleq n,egeq 1} f(p^e) left([e>1] + g(lfloor frac n {p^e} floor,p) ight)+h(n)-h(m)$$
(以上公式以及下图摘自 集训队论文2018 - 朱震霆 - 一些特殊的数论函数求和问题)
接下来我们考虑如何求 $h(x)$ 。
时间复杂度积分算一算就可以知道是 $O(frac {n^{frac 3 4}}{log n})$。
在求 $g(n,m)$ 的直接爆搜就好了,连记忆化都不用!(但这个我不会证明,为什么是对的自己看论文)
具体代码实现主要参见模板部分。
模板题以及模板代码
51Nod1222 最小公倍数计数
题意
给定 $a,b$, 求
$$sum_{n=a}^b sum_{i=1}^n sum_{j=1}^i [{ m lcm } (i,j) = n]$$
$$a,bleq 10^{11}$$
$${ m Time Limit } = 6s$$
题解
先差分一下,转化成求前缀和。
先把原题的统计无序数对转化成统计有序数对,最终 $ans' = (ans+n)/2$ 即可。
设集合 $P$ 表示素数集合。
设 $c(n,p)$ 表示最大的使得 $p^{c(n,p)}|n$ 的数。
若 ${ m lcm } (i,j) = n$ ,则
$$forall p in P, c(n,p)=max(c(i,p),c(j,p))$$
所以,$forall pin P$ ,$c(i,p)$ 和 $c(j,p)$ 共有 $2c(n,p) +1 $ 种取值方法。
所以,设
$$n=prod_i p_i^{k_i} (p_iin P)$$
则
$$ sum_{i=1}^n sum_{j=1}^i [{ m lcm } (i,j) = n] = prod_t (2k_t+1) $$
显然这个式子满足 Min_25 筛的条件,直接筛就好了。
关于本题,还有一些其他做法,详见https://www.cnblogs.com/zhouzhendong/p/51Nod1222.html
代码
#pragma GCC optimize("Ofast","inline") #include <bits/stdc++.h> #define clr(x) memset(x,0,sizeof (x)) using namespace std; typedef long long LL; LL read(){ LL x=0,f=0; char ch=getchar(); while (!isdigit(ch)) f|=ch=='-',ch=getchar(); while (isdigit(ch)) x=(x<<1)+(x<<3)+(ch^48),ch=getchar(); return f?-x:x; } const int Base=1000005,N=Base*2+5; LL n,cn,a,b,base; LL h[N],ps[N],cnt; LL p[N],pcnt; #define ID(i) ((i)<=base?i:cnt-cn/(i)+1) LL f(int e){ return e*2+1; } LL g(LL n,LL m){ LL ans=max(0LL,h[ID(n)]-h[ID(p[m-1])]); for (int i=m;i<=pcnt&&p[i]*p[i]<=n;i++){ LL nn=n/p[i]; for (int e=1;nn>0;e++,nn/=p[i]) ans+=f(e)*((e>1)+g(nn,i+1)); } return ans; } LL _solve(LL _n){ cn=n=_n,base=(LL)sqrt(n),cnt=pcnt=0; for (LL i=1;i<=n;i=ps[cnt]+1) ps[++cnt]=n/(n/i),h[cnt]=ps[cnt]-1; p[0]=1; for (LL i=2;i<=base;i++) if (h[i]!=h[i-1]){ p[++pcnt]=i;//顺便把质数筛出来 LL i2=i*i; for (LL j=cnt;ps[j]>=i2;j--) h[j]-=h[ID(ps[j]/i)]-(pcnt-1); } for (LL i=1;i<=cnt;i++) h[i]*=3; return g(n,1)+1; } LL solve(LL n){ return (_solve(n)+n)/2; } int main(){ a=read(),b=read(); cout<<solve(b)-solve(a-1)<<endl; return 0; }