min_25筛入门

最新推荐文章于 2024-07-14 21:17:23 发布

gmh77

最新推荐文章于 2024-07-14 21:17:23 发布

阅读量1.1k

点赞数 2

分类专栏：算法详解 min_25筛/洲阁筛数论 LOJ题解

本文链接：https://blog.csdn.net/gmh77/article/details/88142031

版权

算法详解同时被 3 个专栏收录

44 篇文章 0 订阅

订阅专栏

数论

10 篇文章 0 订阅

订阅专栏

min_25筛/洲阁筛

3 篇文章 0 订阅

订阅专栏

前言

~~我好菜啊现在才会~~
好像踩爆了洲鸽（划掉）阁筛？
（貌似min_25筛就是洲阁筛的另一种写法）
如果有错欢迎dalao指出因为~~我特么在学的时候就被坑过N次了~~

min_25筛

用来求积性函数的前缀和
$\sum_{i=1}^{n}F(i)$ ，并且F(i)满足一些性质：
①当i为质数时F(i^c)要存在一个较快的计算方法
②当i为质数时F(i)可以被表示成一个多项式的形式
（其实只要能快速计算出质数的前缀和就行了）
什么意思？就是F(i)=a₀x⁰+a₁x¹+a₂x²+…a_nxⁿ
（话说不是多项式也可以吧，只要能用一些积性函数来表示就可以了？求解）
只有满足这个性质的函数才能用min_25筛来求

考虑把1~n中i分成1、质数、合数三类
先撕烤一下如何计算质数的F(i)
。。。好像不能直接算

质数

于是引入一个函数g(n,j)
注：f(i)只是一个用来构造的函数，与最初所求的F(i)不同！！！
设f(i)=x^k，显然f(i)是个积性函数
于是可以求出每个质数i的 F(i)多项式中的每项的f(i)（即xⁱ），再乘以a_i后相加就可以拼出原来的F(i)（i为质数）了
所以求g实际上是个构造的过程

设P[i]表示第i个质因子，g(n,j)表示1~n中当i为质数或i的最小质因子>P[i]的f(i)之和，即
$g(n,j)=\sum_{i=2}^{n}f(i)[i是质数或i的最小质因子>P[i]（即≥P[i+1]）]$
（不把1算进去是因为1*质数=质数，而每次要筛的都是合数）
那么因为f(i)=x^k，所以可以快速算出g(n,0)
举个栗子
当k=0时，g(n,0)=n
当k=1时，g(n,0)=(1+n)n/2
当k=2时，g(n,0)=n(n+1)(2n+1)/6
（平方和公式有兴趣可以推一推，大概就是用暴力和高斯求和列方程再解，更高次幂的基本不会出现）
~~好像有个叫自然数幂和的东西~~
。。。

因为目标是求出所有i为质数上的f(i)，再将其组合成F(i)，所以考虑从g(1~n,0)推到g(n,|P|)
那么从g(n,j-1)转移到g(n,j)时，可以发现除去了最小质因子为P[j]的合数
因为最小质因子为P[j]的合数最小为P[j]²，如果P[j]²>n，那么显然不用删了
否则就减去 $f(P[j])g(\left \lfloor \frac{n}{P[j]} \right \rfloor,j-1)$ ，因为最小质因子为P[j]，所以删掉P[j]后的最小质因子≥P[j]，即 $g(\left \lfloor \frac{n}{P[j]} \right \rfloor,j-1)$
然而可以发现这样减多了，把1~P[j]-1中质数*P[j]构成的合数多减了一次，因为这些数在之前就已经减过了，所以要加回 $f (P [j]) g (P [j] - 1, j - 1)$
（为什么不是1~n？因为P[j]~n的质数*P[j]是要减去的合数）
综上所述，g(n,j)的转移如下
$g(n,j)=\left\{\begin{matrix} g(n,j-1) & P[j]^2>n\\ g(n,j-1)-f(P[j])(g(\left \lfloor \frac{n}{P[j]} \right \rfloor,j-1)-g(P[j]-1,j-1)) & P[j]^2 \leqslant n \end{matrix}\right.$
于是一波时间玄学的操♂作后可以求出g(n,|P|)，即所有i为质数的f(i)之和，接着把各种g(n,|P|)加减可以得出 $\sum_{i=1}^{n}F(i)[i为质数]$
可以发现求g的过程类似埃筛，能在筛质数时求出，不必另求一遍
（实测可以快很多）

这里有一个问题，那就是如果把i从1~n都枚举一遍会很虚，因为min_25筛一般解决的n都是10¹⁰级别的
所以怎么解决~~留到后面再讲~~
（注：因为P[j]≤√n，所以P[j]-1<√n，g(P[j]-1,j-1)属于可以直接算的哪一类）

合数

现在已经可以求出质数的贡献了，考虑求合数的贡献
设s(n,j)表示
$s(n,j)=\sum_{i=2}^{n}{F(i)[i的最小质因子≥P[j]]}$ （不考虑为1的特殊情况，所以最后要加上F(1)）
（注意与g的条件不同，而且这里直接求的是F(i)）
对于某个s(n,j)，首先先计算出质数的贡献，就是各种（因为F(i)是由至少一个函数拼成的，所以可能不止一种）g(n,|P|)减去对应的f(P[1~j-1])再求和，这里的f(P[1~j-1])可以前缀和预处理
~~欸预处理不是O(n)的吗~~
则接着加上合数的贡献，每次枚举最小质因子P[k]（k≥j）以及其的幂e，加上的便是 $F(P[k]^e)s(\left \lfloor \frac{n}{P[k]^e} \right \rfloor,k+1)$ ，即删去P[k]^e*剩下质因子≥P[k+1]的数（最小质因子为P[k]的合数）
然而可以发现这样并没有算i=P[k]^e的情况，因为算的要么是单个的质数要么是多个质数的若干次幂且其中最大的一定是单个质数，没有算单个质数的幂以及其他的一些情况
其实把 $F(P[k]^e)[e \geqslant 2]$ 加上就完事了

然而这样要枚举到O(n)级别的质数
~~实际并不用~~
其实跟g类似，如果P[k]²>n，那么e≤1，首先不存在最小质因子为P[k]的合数，所以 $F(P[k]^e)s(\left \lfloor \frac{n}{P[k]^e} \right \rfloor,k+1)=0$
而且因为e=1，F(P[k])已经在g中算过了，所以当P[k]²>n时就可以直接退了
综上所述，f(n,j)的转移如下：
$s(n,j)=g(n,|P|)-\sum _{i=1}^{j-1}{f(P[i])}+\sum _{k\geqslant j \; and \; P[k]^2\leqslant n}{\sum _{e\geqslant 1}{ {F(P[k]^e)s(\left \lfloor \frac{n}{P[k]^e} \right \rfloor,k+1)+F(P[k]^e)[e\geqslant 2]}}}$

那么答案就是F(1)+s(n,1)，可以发现这是一个递归的过程，只可能从n走到n div i
~~因为某些玄学的原因不用记忆化~~
空间复杂度 $O(\sqrt{n})$ ，时间复杂度 $O(\frac{n^{\frac{3}{4}}}{log\;n})$
~~我只会证空间~~

F保证是积性函数但不是完全积性函数，所以要一个个质数分别考虑
而前面所求的f是完全积性函数，所以可以一个一个加进去
（还有个类似求个非递归版本，可以慢速求出所有s(n div i,|P|)~~但不知意义何在~~）

递推版本

~~好吧还是很有意义的~~
如果要考虑所有质数的情况（如jzoj6027）就要用到递推
一般来说这样的题>√n的质数都可以用前缀和处理，剩下的≤√n的就依次递推
来自https://blog.csdn.net/XianHaoMing/article/details/80397777

各种证明

引理1

显然n div i最多只有2√n个取值，则s和g最多只有2√n个取值
对于≤√n的数可以直接存，处理也很方便，对于>√n的可以用n div x来存

引理2

首先有一个引理：n div a div b=n div (a*b)
证明（来自https://blog.csdn.net/semiwaker/article/details/73822107）：
在这里插入图片描述
这样就可以把n div x div P变为n div (x*P)，可以表示为一开始处理的形式（n div i）

引理3

对于>√n的数显然可以被表示成n div x（x≤√n）的形式，每次枚举x，从n div (x*P)转移到n div x
如果n div (x*P)>√n，那就要求出这个数对应的标号x
这里有个问题，对于≤√n的x对应的n div x怎么还原成x？
~~显然就是n div (n div x)~~
倒推法证明：
设 $n = a x + b （ 0 \leq b < x ）$
$\left \lfloor \frac{n}{\left \lfloor \frac{n}{x} \right \rfloor} \right \rfloor=x$
$\left \lfloor \frac{ax+b}{\left \lfloor \frac{ax+b}{x} \right \rfloor} \right \rfloor=x$
$\left \lfloor \frac{ax+b}{a} \right \rfloor=x$
$\left \lfloor x+\frac{b}{a} \right \rfloor=x$
如果 $\frac{b}{a}<1$ 那么结论就可以成立
即 $a > b$
因为 $n = a x + b$
所以 $a=\left \lfloor \frac{n}{x} \right \rfloor$ ， $b=n\;mod\;x$
因为 $\leqslant sqrt(n)$ ，所以 $\left \lfloor \frac{n}{x} \right \rfloor \geqslant sqrt(n)$ ，即 $a\geqslant sqrt(n)$
因为 $b=n\;mod\;x$ ，所以 $b < x$ ，即 $b < s q r t (n)$
所以 $a\geqslant sqrt(n)>b$ ，即当 $\leqslant sqrt(n)$ 时原式成立
（从下往上反推一遍就行了）~~有♂趣~~

例题1

题目描述

https://loj.ac/problem/6053
LOJ#6053. 简单的函数

在这里插入图片描述

*一样的截屏

题解

显然f并不好求，但f(x)=f(x¹)=x xor 1（x为质数）
当x>2时x xor 1=x-1，可以拆成x¹-x⁰
于是min_25分别求x¹和x⁰然后合并求出f在质数时的前缀和
（注意f(2)=3要特判）

code

#include <iostream>
#include <cstdlib>
#include <cstdio>
#include <cmath>
#define fo(a,b,c) for (a=b; a<=c; a++)
#define fd(a,b,c) for (a=b; a>=c; a--)
#define mod 1000000007
#define f(p,c) ((p)^(c)%mod)
using namespace std;

long long p[100001];
long long P[100001];
bool bz[100001];
long long b[100001];
long long g0[100001];
long long G0[100001];
long long g1[100001];
long long G1[100001];
long long ans[100001];
long long Ans[100001];
int sum[100001];
long long n;
int i,j,k,l,len,sq,Sq;

void init()
{
	long long B,S,s=floor(sqrt(n));
	int _i,i,j;
	
	fo(i,1,Sq) g0[i]=i-1,g1[i]=(long long)i*(i+1)/2%mod-1;
	fo(i,1,sq) B=b[i]%mod,G0[i]=B-1,G1[i]=B*(B+1)%mod*500000004%mod-1;
	
	fo(_i,2,s)
	{
		if (!bz[_i])
		{
			p[++len]=_i;
			P[len]=p[len]*p[len];
			
			fo(i,1,sq)
			if (P[len]<=b[i])
			{
				S=b[i]/p[len];
				
				if (S<=Sq)
				{
					G0[i]=(G0[i]-(g0[S]-g0[p[len]-1]))%mod;
					G1[i]=(G1[i]-p[len]*(g1[S]-g1[p[len]-1]))%mod;
				}
				else
				{
					S=n/S;
					G0[i]=(G0[i]-(G0[S]-g0[p[len]-1]))%mod;
					G1[i]=(G1[i]-p[len]*(G1[S]-g1[p[len]-1]))%mod;
				}
			}
			else
			break;
			
			fd(i,Sq,1)
			if (P[len]<=i)
			{
				g0[i]=(g0[i]-(g0[i/p[len]]-g0[p[len]-1]))%mod;
				g1[i]=(g1[i]-p[len]*(g1[i/p[len]]-g1[p[len]-1]))%mod;
			}
			else
			break;
		}
		
		fo(j,1,len)
		if (p[j]*_i<=s)
		{
			bz[p[j]*_i]=1;
			
			if (!(_i%p[j]))
			break;
		}
		else
		break;
	}
	
	fo(i,1,len)
	sum[i]=(sum[i-1]+p[i]-1)%mod;
}

long long dfs(long long s,int j)
{
	long long S,s2;
	int k,l,e;
	
	S=((s<=Sq?ans[s]:Ans[n/s])-sum[j-1])+((j==1)?2:0);//特判f(2)==3
	
	for (k=j; k<=len && s>=P[k]; ++k)
	{
		s2=s/p[k];l=1;
		
		while (s2)
		{
			if (!(k==len && s2<=p[len] || k<len && s2<p[k+1]))//边界dfs(s2,k+1)=0
			S=(S+f(p[k],l)*dfs(s2,k+1))%mod;
			
			if (l>1)
			S=(S+f(p[k],l))%mod;
			
			s2/=p[k];++l;
		}
	}
	
	return S;
}

int main()
{
	scanf("%lld",&n); sq=floor(sqrt(n));Sq=floor(n/sq)-1;
	
	if (n==1)
	{
		printf("1\n");
		return 0;
	}
	
	fo(i,1,sq)
	b[i]=n/i;
	
	init();
	
	fo(i,1,Sq) ans[i]=g1[i]-g0[i];
	fo(i,1,sq) Ans[i]=G1[i]-G0[i];
	
	printf("%lld\n",(dfs(n,1)+1+mod)%mod); //f(1)=1
}