杜教筛

在这里插入图片描述
在这里插入图片描述

1.找积性函数与f卷积配对

2.求出f的前缀和

3.套用杜教筛模板公式

对于比较简单的题,一般f(n)都为已知的常用函数,比如欧拉函数,莫比乌斯函数等,则可直接求得f(n)的前缀和

#include<bits/stdc++.h>
#define ll long long
const int maxn=5e6+5;
using namespace std;
map<int,ll>ans;
int prime[maxn];
bool flag[maxn];
ll phi[maxn];
void init()
{
	phi[1]=1;
	for(int i=2;i<maxn;i++)
	{
		if(!flag[i])
		{
			prime[++prime[0]]=i;
			phi[i]=i-1;
		}
		for(int j=1;j<=prime[0];j++)
		{
			if(i*prime[j]>maxn)break;
			flag[i*prime[j]]=1;
			if(i%prime[j]==0)
			{
				phi[i*prime[j]]=phi[i]*prime[j];
				break;
			}
			else
			phi[i*prime[j]]=phi[i]*phi[prime[j]];
		}
	}
	for(int i=1;i<maxn;i++)
	phi[i]+=phi[i-1];
}//这一块为求f(n)的前缀和 

ll sumphi(int n)//此处最好用int 然后下面用lll转换为long long类型,省时间
{
	if(n<=maxn) return phi[n];
	if(ans[n]) return ans[n];
	ll tmp;
	tmp=1ll*n*(n+1)/2;
	for(int i=2,j;i<=n;i=j+1)//也用int 定义,然后1ll转换,省时间
	{
		j=n/(n/i);
		tmp=tmp-sumphi(n/i)*1ll*(j-i+1);
	}
	ans[n]=tmp;
	return tmp;
}//这一块为模板

int main()
{
	init();
	int n;
	scanf("%d",&n);
	printf("%lld\n",sumphi(n));
}

但有的题目所给出的f(n)不是像欧拉函数这种的

在这里插入图片描述
设 g(n)=n^2-3n+2,便可知 g = f * 1 (此处为卷积)
则由性质 h=f * g→ f=h * ug可得到 f = g * u
在这里插入图片描述
(其实也可由莫比乌斯公式直接得到,此处为了好理解套用杜教筛公式用卷积)
既然已找到了两个函数 g与1 那么就可以套用杜教筛公式了
S(n)=n
(n-1)*(n-2)/3-后面那特定的一串 //前面的多项式是通过计算得到的

g是已知的一个函数,由杜教筛的规定可知要想求S(n)则要求出f(n)的前缀和
这时在求f(n)的前缀和是,首先要把1~n项的f(i)求出来,代码如下

	for(int i=1;i<maxn;++i)
	{
		for(int j=i;j<maxn;j+=i)//满足i是j的约数进行 
		f[j]=(f[j]+mu[j/i]*g(i)%mod)%mod;//g为一已知的函数,防止数过大而取模
	}		
	for(int i=1;i<maxn;++i)//求出前缀和
	f[i]=(f[i-1]+f[i])%mod;

总代码如下

#include <cstdio>
#include <cstring>
#include <algorithm>
#include <map>
#define ll long long
using namespace std;
const int maxn = 1000001;
const int mod = 1e9+7;
ll g(ll n)// g=f*1 中的g函数表达式 
{
	return (n*n%mod-3*n%mod+2)%mod;
}
ll Prime[maxn],Cnt,mu[maxn],f[maxn];
bool flag[maxn];
void init()
{
	mu[1]=1;
	for(int i=2;i<maxn;i++)
	{
		if(!flag[i]) 
		{
			Prime[++Cnt]=i;
			mu[i]=-1;
		}
		for(int j=1;j<=Cnt;j++)
		{
			if(i*Prime[j]>maxn) 
				break;
			flag[i*Prime[j]]=1;
			if(i%Prime[j]==0)
			{
				mu[i*Prime[j]]=0;
				break;
			}
			else
			mu[i*Prime[j]]=-mu[i];
		}
	}
	for(int i=1;i<maxn;i++)//求出个f[i]的值 
	{
		for(int j=i;j<maxn;j+=i)//由于那个约数求和符号,所以应满足i是j的约数进行 
		f[j]=(f[j]+mu[j/i]*g(i)%mod)%mod;
	}		
	for(int i=1;i<maxn;i++)//求出前缀和 
	f[i]=(f[i-1]+f[i])%mod;
}
ll inv(ll a,ll b,ll p)//由于分数取模,所以求逆元 
{
  		ll ret=1;
  		while(b)
  		{
  			if(b&1) ret=ret*a%p;
  			a=a*a%p;
  			b>>=1;
  		}
  		return ret;
}
map<int,ll>ans;
ll solve(ll n)
{
	if(n<=maxn) return f[n];
	if(ans[n]) return ans[n];
	int ret = (n*(n-1)%mod*(n-2)%mod*inv(3,mod-2,mod)%mod)%mod;//由杜教筛公式需要求出g函数的前n项和 
	for(int i=2,j;i<=n;i=j+1)
	{
		j=n/(n/i);
		ret=(ret-(j-i+1)*solve(n/i)%mod)%mod; //由于是g=f*1 所以(j-i+1) 这个就相当于 1函数 ,符合杜教筛公式 
	}
	ans[n]=ret;
	return ret;
}

int main()
{
	init();
	ll T, n;
	scanf("%lld",&T);
	while(T--)
	{
		scanf("%lld",&n);
		printf("%lld\n", (solve(n)+mod)%mod);
	}
}



  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

henulmh

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值