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

先定义几个符号：

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

集合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 eq 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 eq 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整除的结果只有$2sqrt 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}}{lnfrac{n}{x}}dx)$$

由于对数函数增长极慢，不妨将分母用同样的放缩。当然也可以查积分表得到与我相同的结论。这样再对单项式函数做积分就简单了，得到$$O(frac{sqrt n}{lnn}*[2sqrt 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 aoperatorname{li}(a^{frac{1}{x}})}{4lnn}]_{2}^{log_{2}n})$$

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

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

代码：

#include<bits/stdc++.h>
using namespace std;

const int SQR = 100000;
const int mod = 1000000007;

long long n;

int prime[SQR],flag[SQR+5],num;
long long h[SQR*2+5],sqr;
long long fh[SQR*2+5];

void getprime(int N){
    for(int i=2;i<=N;i++){
    if(!flag[i]) prime[++num] = i;
    for(int j=1;j<=num&&i*prime[j] <= N;j++){
        flag[i*prime[j]] = 1;
        if(i%prime[j] == 0) break;
    }
    }
}

void geth(int ft){
    for(int i=1;i<=num;i++){int zt = (ft > 0?prime[i]:1);
    for(int j=1;j<=sqr&&(1ll*prime[i]*prime[i]<=(n/j));j++){
        long long nxt = 1ll*j*prime[i]; int plu = 2*sqr-j+1;
        if(nxt > sqr) h[plu] -= 1ll*zt*(h[n/nxt]-h[prime[i]-1])%mod;
        else h[plu] -= 1ll*zt*(h[2*sqr-nxt+1]-h[prime[i]-1])%mod;
        h[plu]<0?h[plu]+=mod:1;
    }
    for(int j=sqr;j>=2&&(1ll*prime[i]*prime[i]<=j);j--){
        h[j] -= 1ll*zt*(h[j/prime[i]]-h[prime[i]-1])%mod;
        h[j]<0?h[j]+=mod:1;
    }
    }
}

int f(int now,int zt){
    if(now <= sqr){
    if(now <= 1 || prime[zt] >= now) return 0;
    int ans = fh[now]-fh[prime[zt]]; if(ans < 0) ans += mod;
    if(1ll*prime[zt]*prime[zt] >= now) return ans;
    for(int i=zt+1;1ll*prime[i]*prime[i] <= now;i++){
        long long dt = prime[i],k=1;
        while(dt <= now){
        ans += ((prime[i]^k)*(f(now/dt,i)+1))%mod; ans %= mod;
        dt = dt*prime[i];k++;
        }
        ans-=(prime[i]^1); if(ans < 0) ans += mod;
    }
    return ans;
    }else{
    long long rd = n/(2*sqr-now+1);
    int ans = fh[now]-fh[prime[zt]]; if(ans < 0) ans += mod;
    if(1ll*prime[zt]*prime[zt] >= rd) return ans;
    for(int i=zt+1;1ll*prime[i]*prime[i] <= rd&&i<=num;i++){
        long long dt = prime[i],k=1;
        while(dt <= rd){
        if(rd/dt <= sqr) ans += ((prime[i]^k)*(f(rd/dt,i)+1))%mod;
        else ans += ((prime[i]^k)*(f(2*sqr-(2*sqr-now+1)*dt+1,i)+1))%mod;
        ans %= mod; dt = dt*prime[i];k++;
        }
        ans-=(prime[i]^1); if(ans < 0) ans += mod;
    }
    return ans;
    }
}

void work(){
    for(int i=2;i<=sqr;i++){h[i] = (h[i-1] + i)%mod;}
    for(int i=1;i<=sqr;i++){
    long long fak=(n/i)%mod;h[sqr*2-i+1]=((2+fak)*(fak-1)/2)%mod;
    }
    geth(1); for(int i=2;i<=2*sqr;i++) fh[i] += h[i];
    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;
    geth(0); for(int i=2;i<=2*sqr;i++) fh[i] = fh[i]-h[i]+2,fh[i] =(fh[i]+mod)%mod;
    printf("%d",f(2*sqr,0)+1);
}

int main(){
    scanf("%lld",&n); sqr = sqrt(n);
    if(n == 1){puts("1");return 0;}
    getprime(sqr);work();
    return 0;
}