【XSY2731】Div 数论 杜教筛 莫比乌斯反演

2018/03/06 15:31
阅读数 17

题目大意

  定义复数$a+bi$为整数$k$的约数,当且仅当$a$和$b$为整数且存在整数$c$和$d$满足$(a+bi)(c+di)=k$。

  定义复数$a+bi$的实部为$a$,虚部为$b$。

  定义$f(n)$为整数$n$的所有实部大于$0$的约数的实部之和。

  给定正整数$n$,求出$\sum_{i=1}^nf(i)$对$1004535809$取模后得到的值。

  $n\leq {10}^{10}$

题解

  以前看到一个数论题就是反演预处理。

  现在看到一个数论题就是反演杜教筛。

  记$s(n)=\sum_{i|n}i$为$n$的因数和,$S(n)=\sum_{i=1}^ns(i)$

  当$b=0$时答案就是$S(n)$。以下仅考虑$b>0$的情况($b<0$也是一样的)

  设$n=(a+bi)(c+di)$,那么 $$ \begin{cases} ac-bd&=n\ ad+bc&=0 \end{cases}\ \frac{a}{b}=-\frac{c}{d} $$   因为这是一道数论题,设 $$ \begin{align} a&=px\ b&=qx\ c&=py\ d&=-qy\ \gcd(p,q)&=1\ \end{align} $$   这样一组$x,y,p,q$就唯一确定了一组$a,b,c,d$

  记 $$ \begin{align} g(n)&=\sum_{p^2+q^2=n}[\gcd(p,q)=1]p\ G(n)&=\sum_{i=1}^ng(i)\ f(n)&=\sum_{p^2+q^2=n}p\ F(n)&=\sum_{i=1}^nf(i) \end{align} $$   问题转化为求 $$ \begin{align} &\sum_{x,y,p,q>0,[(p,q)=1]}[xy(p^2+q^2)\leq n]px\ =&\sum_{i=1}^n(\sum_{p^2+q^2=i}[\gcd(p,q)=1]p)(\sum_{xy\leq\lfloor\frac{n}{i}\rfloor}x)\ =&\sum_{i=1}^ng(i)S(\lfloor\frac{n}{i}\rfloor) \end{align} $$   那么怎么求$F,G,S$呢? $$ \begin{align} S(n)&=\sum_{i=1}^n\sum_{j|i}j\ &=\sum_{i=1}^ni\lfloor\frac{n}{i}\rfloor\ G(n)&=\sum_{p^2+q^2\leq n}p\ &=\sum_{i=1}^\sqrt ni\lfloor\sqrt{n-i^2}\rfloor\ F(n)&=\sum_{p^2+q^2\leq n}[\gcd(i,j)=1]p\ &=\sum_{i=1}^\sqrt{n}i\mu(i)\sum_{j=1}^\frac{n}{i^2}\sum_{p^2+q^2\leq\frac{n}{i^2}}p\ &=\sum_{i=1}^\sqrt{n}i\mu(i)G(\lfloor\frac{n}{i^2}\rfloor) \end{align} $$   这些东西求一次是$O(\sqrt{n})$的,预处理一下,总的复杂度是$O(n^\frac{2}{3})$,因为每个$n$都是题目给的$n$除以某个东西。

  预处理大家都会,我就不讲了。

  zjt:在$O(n^\frac{2}{3})$内求出所有$F(n$除以某个东西$)$的一类算法都叫杜教筛。

  时间复杂度:$O(n^\frac{2}{3})$

代码

#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;
const ll p=1004535809;
int _gcd[3500][3500];
int gcd(int a,int b)
{
	int &s=_gcd[a][b];
	if(~s)
		return s;
	if(!b)
		return s=a;
	return s=gcd(b,a%b);
}
const int maxn=10000000;
ll n,m;
int vis[10000010];
ll f[10000010];
int miu[10000010];
int b[10000010];
int pri[1000010];
int cnt;
ll g1[10000010];
ll g2[10000010];
ll s[10000010];
int c[10000010];
void init()
{
	int i,j;
	for(i=1;i*i<=maxn;i++)
		for(j=1;i*i+j*j<=maxn;j++)
		{
			if(gcd(i,j)==1)
				(f[i*i+j*j]+=i)%=p;
			(g1[i*i+j*j]+=i)%=p;
		}
	for(i=1;i<=maxn;i++)
	{
		f[i]=(f[i]+f[i-1])%p;
		g1[i]=(g1[i]+g1[i-1])%p;
	}
	miu[1]=1;
	c[1]=1;
	s[1]=1;
	for(i=2;i<=maxn;i++)
	{
		if(!b[i])
		{
			pri[++cnt]=i;
			miu[i]=-1;
			c[i]=i;
			s[i]=i+1;
		}
		for(j=1;j<=cnt&&i*pri[j]<=maxn;j++)
		{
			b[i*pri[j]]=1;
			if(i%pri[j]==0)
			{
				miu[i*pri[j]]=0;
				c[i*pri[j]]=c[i]*pri[j];
				if(c[i]==i)
					s[i*pri[j]]=(s[i]*pri[j]+1)%p;
				else
					s[i*pri[j]]=s[c[i*pri[j]]]*s[i/c[i]]%p;
				break;
			}
			miu[i*pri[j]]=-miu[i];
			c[i*pri[j]]=pri[j];
			s[i*pri[j]]=s[i]*(pri[j]+1)%p;
		}
	}
	for(i=1;i<=maxn;i++)
		s[i]=(s[i]+s[i-1])%p;
}
const ll inv2=502267905;
ll S(ll n)
{
	if(n<=maxn)
		return s[n];
	ll s=0,s2;
	ll i,j;
	for(i=1;i<=n;i=j+1)
	{
		j=n/(n/i);
		s=(s+(i+j)%p*(j-i+1)%p*inv2%p*((n/i)%p))%p;
	}
	return s;
}
ll G(ll n)
{
	if(n<=maxn)
		return g1[n];
	if(vis[m/n]&2)
		return g2[m/n];
	vis[m/n]|=2;
	ll i,s=0,j;
	for(i=1;i*i<=n;i++);
	j=i-1;
	for(i=1;i*i<=n;i++)
	{
		while(j*j>n-i*i)
			j--;
		s=(s+i*j)%p;
	}
	g2[m/n]=s;
	return s;
}
ll F(ll n)
{
	if(n<=maxn)
		return f[n];
	ll i,s=0;
	ll now;
	for(i=1;i*i<=n;i++)
		s=(s+miu[i]*i*G(n/i/i))%p;
	return s;
}
int main()
{
#ifndef ONLINE_JUDGE
	freopen("c.in","r",stdin);
	freopen("c.out","w",stdout);
#endif
	memset(_gcd,-1,sizeof _gcd);
	scanf("%lld",&n);
	m=n;
	init();
	ll ans=0;
	ll i,j;
	ll last=0,now;
	for(i=1;i<=n;i=j+1)
	{
		j=n/(n/i);
		now=F(j);
		ans=(ans+(now-last)%p*S(n/i))%p;
		last=now;
	}
	ans=(ans*2+S(n))%p;
	ans=(ans+p)%p;
	printf("%lld\n",ans);
	return 0;
}
展开阅读全文
mu
打赏
0
0 收藏
分享
加载中
更多评论
打赏
0 评论
0 收藏
0
分享
返回顶部
顶部