LG5325 模板Min_25筛

Posted 2022-03-05 autoint

P5325 【模板】Min_25筛

题目背景

模板题，无背景。

题目描述

定义积性函数$f(x)$，且$f(p^k)=p^k(p^k-1)$（$p$是一个质数），求

$$\\sum_i=1^n f(x)$$

对$10^9+7$取模。

输入输出格式

输入格式：

一行一个整数$n$。

输出格式：

一个整数表示答案。

输入输出样例

输入样例#1： 复制

10

输出样例#1： 复制

263

输入样例#2： 复制

1000000000

输出样例#2： 复制

710164413

说明

$f(1)=1,f(2)=2,f(3)=6,f(4)=12,f(5)=20$

$f(6)=12,f(7)=42,f(8)=56,f(9)=72,f(10)=40$

对于$30\\%$的数据，保证$1\\le n\\le 10^6$。

对于$100\\%$的数据，保证$1\\le n\\le 10^10$

wucstdio的题解

题意

求\\(\\sum_i=1^nf(i)\\)，其中\\(f(i)\\)是一个积性函数，且\\(f(p^k)=p^k(p^k-1)\\)。

Min_25筛

为什么这个筛法叫做Min25筛呢？因为这个算法是Min25发明的。

假设我们要求一个\\(\\sum_i=1^nf(i)\\)，满足\\(f(x)\\)是一个积性函数，且\\(f(p^e)\\)是一个关于\\(p\\)的低阶多项式。

因为多项式可以拆成若干个单项式，所以我们只需要考虑求出\\(f(p)=p^k\\)的前缀和，然后每一项加起来就行了。

那么如何求出每一项的和呢？

• 分类

  让我们先对\\(i\\)按照质数和合数分类：

\\[\\sum_i=1^nf(i)=\\sum_1\\le p\\le nf(p)+\\sum_i=1\\&\\texti is not a prime^nf(i)\\]

然后，我们枚举后面合数的最小质因子以及最小质因子的次数。注意所有合数的最小质因子一定都小于等于\\(\\sqrt n\\)：

\\[\\sum_1\\le p\\le nf(p)+\\sum_1\\le p^e\\le n,1\\le p\\le \\sqrt nf(p^e)\\left(\\sum_1\\le i\\le n/p^e\\&minp>pf(i)\\right)\\]

其中\\(minp\\)表示\\(i\\)的最小质因子，因为公式中文太丑了所以就只好写英文了。

这样，整个式子就变成了两个部分，第一部分是所有质数的\\(f\\)之和，另一部分是枚举最小质因子后，求所有最小质因子大于这个质因子的\\(f\\)之和。

• 辅助函数

  由于\\(n\\)实在太大，所以我们没法用线性筛求出所有质数的\\(f\\)和。

我们考虑一个DP的思路（天哪这是怎么想到的）：我们不知道从哪里找来了一个DP数组\\(g(n,i)\\)，满足

\\[g(n,j)=\\sum_i=1^n[\\texti is a prime or minp$> p_j$]i^k\\]

这里的\\(k\\)就是前面我们说的低阶多项式的一项（此题中\\(k\\in\\1,2\\\\)）。注意\\(i^k\\)并不是我们要求的\\(f\\)，只是一个和\\(f\\)在质数处的取值一样的完全积性函数，这样后面计算起来比较方便。

这个式子说人话就是\\(g(n,j)\\)表示求\\(1\\)到\\(n\\)之间所有满足条件的数的\\(k\\)次方和，条件就是要么是质数要么最小质因子大于\\(p_j\\)。

我们考虑\\(g(n,j-1)\\)如何转移到\\(g(n,j)\\)。随着\\(j\\)的增大，满足条件的数变少了，所以我们需要减去一些原来满足条件而现在不满足条件的数。

这些数应该是最小质因子恰好为\\(p_j\\)的合数。

我们可以提出来一个\\(p_j\\)作为最小质因子，这样剩下的就不能有小于它的质因子了，也就是\\(g\\left(\\dfracnp_j,j\\right)-g(p_j-1,j-1)\\)，后面那个\\(g\\)是为了把所有的质数去掉。

这样我们就得到了\\(g\\)的递推式：

\\[g(n,j)=g(n,j-1)-p_j^k\\left(g\\left(\\dfracnp_j,j\\right)-g(p_j-1,j-1)\\right)\\]

完全积性函数的好处在这里就体现出来了：由于只提出了一个\\(p_i\\)，所以后面还有可能有\\(p_i\\)这个因子，如果是完全积性函数的话就可以将函数值直接相乘，而不用管是否互质。

注意到后面的\\(g(p_j-1,j-1)\\)其实就是前\\(j-1\\)个质数的\\(k\\)次方和。由于\\(p_j\\le \\sqrt n\\)，所以这一部分可以用线性筛预处理，我们设\\(sp_n=\\sum_i=1^np_i^k\\)，也就是前\\(n\\)个质数的\\(k\\)次方和。

\\(1\\)到\\(n\\)中所有质数的\\(k\\)次方和其实就是\\(g(n,x)\\)，其中\\(p_x\\)是最后一个小于等于\\(\\sqrt n\\)的质数。为了方便，我们把它记作\\(g(n)\\)。

但是因为\\(n\\)太大，我们还是没法对于每一个\\(n\\)求出\\(g(n,x)\\)，所以我们可以想到另一个重要的结论：

\\[\\left\\lfloor\\dfrac\\lfloor\\dfrac na\\rfloorb\\right\\rfloor=\\lfloor\\dfracnab\\rfloor\\]

也就是说，无论你每一次把\\(n\\)除以几，最后你能得到的数一定是某一个\\(\\lfloor\\dfrac nx\\rfloor\\)，所以我们没必要算出来所有的\\(n\\)，只需要算出可以写成\\(\\lfloor\\dfracnx\\rfloor\\)这种形式的数，这样的数一共有\\(O(\\sqrt n)\\)个。

那么我们如何存储这\\(\\sqrt n\\)个数呢？

首先我们不能直接下标访问，这样下标可以到\\(n\\)。我们需要对下标离散化。

但是离散化之后，我们还需要知道对于每一个\\(\\lfloor\\dfrac nx\\rfloor\\)，它对应的下标是什么。

如果偷懒的话可以用map，但是时间复杂度会多一个\\(\\log\\)。我们可以用\\(ind1[x]\\)表示\\(x\\)这个数对应的数组下标，\\(ind2[x]\\)表示\\(n/x\\)这个数对应的下标。这样两个\\(ind\\)数组最大都只会到\\(\\sqrt n\\)。

具体实现可以看代码。数组的记录上需要精细实现一下。

• 求解答案

  答案就是先求出所有质数的函数和，然后先枚举了一个\\(p^e\\)，再枚举最小质因子大于\\(p\\)的数。

我们还是可以考虑DP的思想。设\\(S(n,x)\\)表示求\\(1\\)到\\(n\\)中所有最小质因子大于\\(p_x\\)的函数值之和，注意这里是\\(f\\)而不是\\(k\\)次方。答案就是\\(S(n,0)\\)。

我们将满足条件的数分成两部分，第一部分是大于\\(p_x\\)的质数，也就是\\(g(n)-sp_x\\)，另一部分是最小质因子大于\\(p_x\\)的合数，枚举最小质因子：

\\[S(n,x)=g(n)-sp_x+\\sum_p_k^e\\le n\\&k>xf(p_k^e)\\left(S\\left(\\dfracnp_k^e,k\\right)+[e\\neq 1]\\right)\\]

这样问题就解决了，我们可以递归求解这个问题。根据某玄学定理，不需要记忆化。

一些细节
\\(1\\)既不是质数也不是合数，不含任何一个质因子，那么求解的过程中\\(g\\)和\\(S\\)到底是否包含\\(1\\)呢？其实是否包含都可以，但是处理上略有差别。我的\\(g\\)和\\(S\\)都没有包含\\(1\\)，只需要最后加一就行了。

min25筛的时间复杂度据说是\\(O\\left(\\dfracn^3/4\\log n\\right)\\)，也有人说是\\(O(n^1-\\epsilon)\\)，在这道题上大致是1e10跑1s左右的样子。但是这个算法常数很小，具体表现参加WC2019课件里面的一张图：（灰色的是min25）

#include<bits/stdc++.h>
#define il inline
#define co const
template<class T>T read()
    T data=0,w=1;char ch=getchar();
    for(;!isdigit(ch);ch=getchar())if(ch=='-') w=-w;
    for(;isdigit(ch);ch=getchar()) data=data*10+ch-'0';
    return data*w;

template<class T>il T read(T&x) return x=read<T>();
typedef long long LL;
using namespace std;

co int mod=1e9+7,i2=500000004,i6=166666668;
il int add(int a,int b)
    return (a+=b)>=mod?a-mod:a;

il int mul(int a,int b)
    return (LL)a*b%mod;

//int fpow(int a,int b) // calc inv
//  int ans=1;
//  for(;b;b>>=1,a=mul(a,a))
//      if(b&1) ans=mul(ans,a);
//  return ans;
//
co int N=2e5+1; // edit 1: size
LL n;
int sqr,pri[N],num,sp1[N],sp2[N];
LL w[N];int tot,g1[N],g2[N],ind1[N],ind2[N];
void init()
    pri[1]=1;
    for(int i=2;i<=sqr;++i)
        if(!pri[i])
            pri[++num]=i;
            sp1[num]=add(sp1[num-1],i);
            sp2[num]=add(sp2[num-1],mul(i,i));
        
        for(int j=1;j<=num&&i*pri[j]<=sqr;++j)
            pri[i*pri[j]]=1;
            if(i%pri[j]==0) break;
        
    

int S(LL x,int y)
    if(pri[y]>=x) return 0;
    int k=x<=sqr?ind1[x]:ind2[n/x];
    int ans=add(add(g2[k],mod-g1[k]),mod-add(sp2[y],mod-sp1[y]));
    for(int i=y+1;i<=num&&(LL)pri[i]*pri[i]<=x;++i)
        LL p=pri[i];
        for(int e=1;p<=x;++e,p=p*pri[i])
            int t=p%mod;
            ans=add(ans,mul(t,mul(t-1,add(S(x/p,i),e!=1))));
        
    
    return ans;

int main()
    sqr=ceil(sqrt(read(n))),init();
    // part 1
    for(LL l=1,r;l<=n;l=r+1) // init g(w,0)
        r=n/(n/l);
        w[++tot]=n/l; // big->small
        g1[tot]=g2[tot]=w[tot]%mod;
        g1[tot]=mul(i2,mul(g1[tot],g1[tot]+1)),g1[tot]=add(g1[tot],mod-1);
        g2[tot]=mul(i6,mul(g2[tot],mul(g2[tot]+1,2*g2[tot]+1))),g2[tot]=add(g2[tot],mod-1);
        if(n/l<=sqr) ind1[n/l]=tot;
        else ind2[n/(n/l)]=tot;
    
    for(int i=1;i<=num;++i) // DP g(w,i)
        for(int j=1;j<=tot&&(LL)pri[i]*pri[i]<=w[j];++j)
            int k=w[j]/pri[i]<=sqr?ind1[w[j]/pri[i]]:ind2[n/(w[j]/pri[i])];
            g1[j]=add(g1[j],mod-mul(pri[i],add(g1[k],mod-sp1[i-1])));
            g2[j]=add(g2[j],mod-mul(pri[i],mul(pri[i],add(g2[k],mod-sp2[i-1]))));
        
    // part 2
    printf("%d\\n",add(S(n,0),1));
    return 0;

这题数组空间，我开刚好\\(\\sqrtn\\)，也就是1e5，就会炸掉。试了一下开到2e5才过。大概是数据有问题。