Loj#6052. 「雅礼集训 2017 Day11」DIV

scb tql！！！！！！！！！！！

不复读一遍scb博客了。简略写一些。

假如要算\(\lfloor \frac{n}{i} \rfloor\)处的函数值，如果能在\(O(\sqrt n)\)的时间内算出单点函数值，那么复杂度为

\[\sum\limits_{x=1}^{\sqrt n}\sqrt{\frac{n}{x}}+\sqrt x=\sqrt n\int_0^{\sqrt n}x^{-\frac{1}{2}}dx+\int_0^{\sqrt n}x^{\frac{1}{2}}dx=O(n^{\frac{3}{4}}) \]

基于这个，我们来考虑转化之后的式子

\[\sum\limits_{d=1}^{\sqrt n}\mu(d)d\sum\limits_{g=1}^{\lfloor \frac{n}{d^2} \rfloor}g\sum\limits_{a=1}^{\sqrt{\lfloor\frac{n}{gd^2}\rfloor}}a\sum\limits_{b=1}^{\sqrt{\lfloor\frac{n}{gd^2}\rfloor}}\lfloor\frac{n}{gd^2(a^2+b^2)}\rfloor \]

考虑设\(f(m)=\sum\limits_{a=1}^{\sqrt a}a\sum\limits_{b=1}^{\sqrt b}\lfloor\frac{m}{a^2+b^2}\rfloor\)这玩意不好求，再来一个\(g(m)=\sum\limits_{a=1}^{\sqrt m}a\sum\limits_{b=1}^{\sqrt m}[a^2+b^2\leq m]\)。这样就可以通过枚举\(\lfloor\frac{m}{a^2+b^2}\rfloor\)，用上文所说的时间复杂度求出\(f\)。

再来一个\(h(m)=\sum\limits_{i=1}^mf(\lfloor\frac{m}{i}\rfloor)i\)。这个就再来一个上文所说时间复杂度求出。用\(h\)即可算出最终答案。

其实这本质上就是杜教筛。

\(b=0\)就是约数和，随便求一求。总时间复杂度\(O(n^{\frac{3}{4}})\)。

然后再丧心病狂卡卡常就能卡过去了。

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#include<queue>
#include<vector>
#include<ctime>
using namespace std;
typedef unsigned long long ll;
#define N 200102
const int p=1004535809,iv=(p+1)>>1;
ll n,b[N][2],sqr;
int tot,c[N],f[N],id1[N],id2[N],h[N],prim[N],cnt,mu[N],fjd;
ll tr[N];
bool isp[N];
inline void upd(int &x,int y){x+=x+y>=p?y-p:y;}
inline int S(ll x){if(x>=p)x%=p;return (x*(x+1)/2)%p;}
inline int calf(int i)
{
	ll td=n/b[i][0];int la=0;
		for(register ll ii=1,j;ii<=td;ii=j+1)
		{ll dd=td/ii;
			j=td/dd;
			int pos=j>=fjd?id1[n/j]:id2[tr[j]];
			upd(f[i],dd*(c[pos]-c[la]+p)%p);la=pos;
		}
	return f[i];
}
inline int cal(ll x)
{
	int pos=x<=sqr?id1[x]:id2[n/x];
	if(h[pos]>=0)return h[pos];
	int ret=0;int la=0;
	for(register ll i=1,j;i<=x;i=j+1)
	{ll td=x/i;
		j=x/td;
		int pos=td<=sqr?id1[td]:id2[n/td];
		int ts=S(j);
		int sum=(ts-la+p)%p;la=ts;
		upd(ret,1ll*f[pos]*sum%p);
	}return h[pos]=ret;
}
int main()
{memset(h,-1,sizeof(h));
	for(int i=2;i<=200000;i++)isp[i]=1;mu[1]=1;
	for(int i=2;i<=200000;i++)
	{
		if(isp[i])prim[++cnt]=i,mu[i]=-1;
		for(int j=1;j<=cnt&&prim[j]*i<=200000;j++)
		{
			isp[prim[j]*i]=0;
			if(i%prim[j]==0){mu[prim[j]*i]=0;break;}
			mu[prim[j]*i]=-mu[i];
		}
	}
	scanf("%llu",&n);sqr=sqrt(n);
	for(ll i=1,j;i<=n;i=j+1)
	{ll td=n/i;
		j=n/td;
		b[++tot][0]=i,b[tot][1]=j;
		if(td<=sqr)id1[td]=tot;
		else id2[n/td]=tot;
		for(int ii=1;1ll*ii*ii<=j;ii++)
		{int te=sqrt(j-1ll*ii*ii);
			upd(c[tot],1ll*ii*te%p);
		}
	}
	ll l=1,r=n,mid;
	while(l<=r)
	{
		mid=(l+r)>>1;
		if(n/mid<=sqr)fjd=mid,r=mid-1;
		else l=mid+1;
	}
	for(int i=1;i<=min(100000*1ull,n);i++)tr[i]=n/(n/i);
	for(int i=1;i<=tot;i++)f[i]=calf(i);
	int ans=0;
	for(int i=1;1ll*i*i<=n;i++)
	{
		int te=1ll*(mu[i]+p)%p*i%p*cal(n/i/i)%p;
		upd(ans,te);
	}
	ans=2ll*ans%p;
	for(ll i=1,j;i<=n;i=j+1)
	{
		j=n/(n/i);
		int te=1ll*(n/i)%p*((i+j)%p)%p*((j-i+1)%p)%p*iv%p;
		upd(ans,te);
	}
	printf("%d\n",ans);
}