LOJ6053 简单的函数【Min_25筛】【埃拉托斯特尼筛】

最新推荐文章于 2021-03-05 17:27:03 发布

weixin_30902675

最新推荐文章于 2021-03-05 17:27:03 发布

阅读量88

点赞数

原文链接：http://www.cnblogs.com/Menhera/p/9226649.html

版权

先定义几个符号：

[]：若方括号内为一个值，则向下取整，否则为布尔判断

集合P：素数集合。

题目分析：

题目是一个积性函数。做法之一是洲阁筛，也可以采用Min_25筛。

对于一个可以进行Min_25筛法的积性函数，它需要满足与洲阁筛相同的条件，即：

对于$f(p), p \in P$，它可以多项式表出。对于$f(p^k),p \in P$可以被快速计算出。

这道题中$f(p) = p-1$再对$2$进行修正即可。

对于1的情况我们单独考虑，现在我们对答案进行一些变换。

$$\sum_{i=2}^{n}f(i) = \sum_{2<=p^e<=n, p \in P} f(p^e)(1+\sum_{2<=i<=[\frac{n}{p^e}],i质因子大于p}f(i)) = \sum_{2<=p<=\sqrt n,p^e<=n, p \in P} f(p^e)(1+\sum_{2<=i<=[\frac{n}{p^e}],i质因子大于p}f(i))+\sum _{p \in P,\sqrt n<p<=n}f(p)$$

我们注意到中间刺眼的中文，不妨以此为根据设置dp状态。

令$$g_{n,m}=\sum_{i=2,i的质因子大于m}^{n}f(i)$$

它的递推式不难写出来。根据上面那串式子:

$$g_{n,m}=\sum_{i=m,i \in P}^{\sqrt n}[p^e \leq n]([e \neq 1]+\sum_{i=p+1,i的质因子大于p}^{[\frac{n}{p^e}]}f(i))+\sum_{i=m+1,i \in P}^{n}f(i) = \sum_{i=m,i \in P}^{\sqrt n}[p^e \leq n]([e \neq 1]+g_{[\frac{n}{p^e}],p})+\sum_{i=m+1,i \in P}^{n}f(i)$$ 之后我们可以发现后面那一部分等于$m+1$到$n$的所有素数.

朱老大的论文似乎在这个式子上有一些问题，这个式子是我加了修正的，可以对比原式找到区别。

现在我们重点考虑后面的部分，质数前缀和的求解。我们可以用DP模拟埃拉托斯特尼筛法来完成。现在对f(p)这个多项式的每个单项式分开考虑。

令$$h_{n,m}=\sum_{i=2,i为素数或i质因子大于m}f(i)$$

由埃氏筛可知$$h_{n,m}=h_{n,m-1}-p_m^s*(h_{[\frac{n}{p}],m-1}-h_{p_m,m-1})$$起始状态将所有数看做质数，做k次幂和。

由于n以下的合数至多的因子大小是$\sqrt n$,所以对于一个n只用筛前$\sqrt n$个素数。

g数组同理。

这样我们就完成整个过程了，现在我们来思考时间复杂度上的问题。这一部分我的证明存在争议，与论文上的有出入。希望有人可以帮助我确认出入之处。

先给出两个引理吧。

素数定理:$\pi(x) \approx \frac{x}{lnx} \approx li(x)$

数论分块:对于一个数n，它对1到n整除的结果只有$2\sqrt n$个，分别是1到$\sqrt n$与$n/\sqrt n$到$n/1$

首先是dp模拟埃氏筛。对于数论分块中的每一个数，每个小于它的根号的素数都被考虑，所以时间为，

$$O(\sum_{i=2}^{\sqrt n}\frac{\sqrt i}{lni})+O(\sum_{i=1}^{\sqrt n}\frac{\sqrt \frac{n}{i}}{ln{\frac{n}{i}}})$$

对于前面的大O,将$\sqrt i$放缩到$\sqrt n$,$lni$放缩到$lnn$,因为根号的增长远大于对数，所以放缩成立。前面被证明是$O(\frac{n^{0.75}}{lnn})$的。

后面的比较复杂，将它看做积分，即$$O(\int_1^{\sqrt n}\frac{\sqrt \frac{n}{x}}{ln\frac{n}{x}}dx)$$

由于对数函数增长极慢，不妨将分母用同样的放缩。当然也可以查积分表得到与我相同的结论。这样再对单项式函数做积分就简单了，得到$$O(\frac{\sqrt n}{lnn}*[2\sqrt x]_1^{\sqrt n})=O(\frac{n^{0.75}}{lnn})$$。

埃氏筛的复杂度得证，我们再来看看Min_25筛的表现。论文中通过证明得出一个与运行效率不相符的时间复杂度，我试着用其它方法得到它的复杂度。

观察Min_25的特点，它与DP埃氏筛的区别在于多枚举了质数的次幂。这样我们不妨对于质数的次幂分开考虑。

首先是质数本身，如果只枚举它，时间与埃氏筛DP相同。实际上它代表的是二次方，否则不被考虑。

然后是质数的平方，它要小于等于n。接着是质数的三次方，它也要小于等于n，然后是四次方，五次方。我们对第x次方小于等于n做一些转化使得他们统一。

$p^x \leq n \leftrightarrow p \leq n^{1/x} \leftrightarrow p^2 \leq n^{2/x}$

这时候我们的复杂度相当于对$n^{2/2},n^{2/3},n^{2/4}$做埃氏筛DP。每个的时间在上面分析了，所以下面直接转化成积分的形式。

$$O(\int_{2}^{log_{2}n}\frac{n^{\frac{3}{2x}}}{ln{n^{\frac{2}{x}}}}dx)$$

$$=O(\frac{1}{lnn}\int_{2}^{log_{2}n}\frac{n^{\frac{3}{2x}}}{\frac{2}{x}}dx)$$

$$=O(\frac{1}{2lnn}\int_{2}^{log_{2}n}{x*n^{\frac{3}{2x}}}dx)$$

$$=O(\frac{1}{2lnn}\int_{2}^{log_{2}n}{x*a^{\frac{1}{x}}}dx),a = n^{\frac{3}{2}}$$

接下俩数学好的可以换元，令$y=\frac{1}{x}$,再令$z = a^y$，求$$\int \frac{1}{ln^3{z}dz}$$即可。

现在跳过这些步骤，得到结论是

$$O([\frac{x(x+\ln a)a^\frac{1}{x}-\ln^2 a\operatorname{li}(a^{\frac{1}{x}})}{4lnn}]_{2}^{log_{2}n})$$

这个答案是负的，证明我可能某个步骤写错了。但是代入2可以得到时间是$O(n^{\frac{3}{4}})$，递归的过程不会重复计算。

继续写下去我也写不动了。到此打止。

代码：

 1 #include<bits/stdc++.h>
 2 using namespace std;
 3 
 4 const int SQR = 100000;
 5 const int mod = 1000000007;
 6 
 7 long long n;
 8 
 9 int prime[SQR],flag[SQR+5],num;
10 long long h[SQR*2+5],sqr;
11 long long fh[SQR*2+5];
12 
13 void getprime(int N){
14     for(int i=2;i<=N;i++){
15     if(!flag[i]) prime[++num] = i;
16     for(int j=1;j<=num&&i*prime[j] <= N;j++){
17         flag[i*prime[j]] = 1;
18         if(i%prime[j] == 0) break;
19     }
20     }
21 }
22 
23 void geth(int ft){
24     for(int i=1;i<=num;i++){int zt = (ft > 0?prime[i]:1);
25     for(int j=1;j<=sqr&&(1ll*prime[i]*prime[i]<=(n/j));j++){
26         long long nxt = 1ll*j*prime[i]; int plu = 2*sqr-j+1;
27         if(nxt > sqr) h[plu] -= 1ll*zt*(h[n/nxt]-h[prime[i]-1])%mod;
28         else h[plu] -= 1ll*zt*(h[2*sqr-nxt+1]-h[prime[i]-1])%mod;
29         h[plu]<0?h[plu]+=mod:1;
30     }
31     for(int j=sqr;j>=2&&(1ll*prime[i]*prime[i]<=j);j--){
32         h[j] -= 1ll*zt*(h[j/prime[i]]-h[prime[i]-1])%mod;
33         h[j]<0?h[j]+=mod:1;
34     }
35     }
36 }
37 
38 int f(int now,int zt){
39     if(now <= sqr){
40     if(now <= 1 || prime[zt] >= now) return 0;
41     int ans = fh[now]-fh[prime[zt]]; if(ans < 0) ans += mod;
42     if(1ll*prime[zt]*prime[zt] >= now) return ans;
43     for(int i=zt+1;1ll*prime[i]*prime[i] <= now;i++){
44         long long dt = prime[i],k=1;
45         while(dt <= now){
46         ans += ((prime[i]^k)*(f(now/dt,i)+1))%mod; ans %= mod;
47         dt = dt*prime[i];k++;
48         }
49         ans-=(prime[i]^1); if(ans < 0) ans += mod;
50     }
51     return ans;
52     }else{
53     long long rd = n/(2*sqr-now+1);
54     int ans = fh[now]-fh[prime[zt]]; if(ans < 0) ans += mod;
55     if(1ll*prime[zt]*prime[zt] >= rd) return ans;
56     for(int i=zt+1;1ll*prime[i]*prime[i] <= rd&&i<=num;i++){
57         long long dt = prime[i],k=1;
58         while(dt <= rd){
59         if(rd/dt <= sqr) ans += ((prime[i]^k)*(f(rd/dt,i)+1))%mod;
60         else ans += ((prime[i]^k)*(f(2*sqr-(2*sqr-now+1)*dt+1,i)+1))%mod;
61         ans %= mod; dt = dt*prime[i];k++;
62         }
63         ans-=(prime[i]^1); if(ans < 0) ans += mod;
64     }
65     return ans;
66     }
67 }
68 
69 void work(){
70     for(int i=2;i<=sqr;i++){h[i] = (h[i-1] + i)%mod;}
71     for(int i=1;i<=sqr;i++){
72     long long fak=(n/i)%mod;h[sqr*2-i+1]=((2+fak)*(fak-1)/2)%mod;
73     }
74     geth(1); for(int i=2;i<=2*sqr;i++) fh[i] += h[i];
75     for(int i=2;i<=sqr;i++) h[i] = i-1; for(int i=1;i<=sqr;i++) h[2*sqr-i+1]=(n/i-1)%mod;
76     geth(0); for(int i=2;i<=2*sqr;i++) fh[i] = fh[i]-h[i]+2,fh[i] =(fh[i]+mod)%mod;
77     printf("%d",f(2*sqr,0)+1);
78 }
79 
80 int main(){
81     scanf("%lld",&n); sqr = sqrt(n);
82     if(n == 1){puts("1");return 0;}
83     getprime(sqr);work();
84     return 0;
85 }