【XSY3532】序列计数（数论）

记录一下推导这题的时候发现的一个东西。

按照我一开始的推法，推出来是这个东西：

\[\begin{aligned} \sum_{i=1}^n\mu^2(i)\left\lfloor\sqrt{\frac{n}{i}}\right\rfloor^2 \end{aligned} \]

其中 \(\left\lfloor\sqrt{\dfrac{n}{i}}\right\rfloor^2\) 就是在 \(\left[1,\dfrac{n}{i}\right]\) 中最大的完全平方数。

这个可以 min25 筛 \(\mu^2\) 配上整除分块做，但不能过 \(n\leq 10^{14}\)，结果一直推不下去。然后换了种推法，得到：

\[2\left(\sum_{i=1}^{\sqrt n}\varphi(i)\lfloor\frac{n}{i^2}\rfloor\right)-n \]

这样就有 \(O(\sqrt n)\) 的做法了。

一个疑惑的点是：为什么上下不能互推？

因为我太菜了。 其实是可以的，首先第一个问题是这个 \(\left\lfloor\sqrt{\dfrac{n}{i}}\right\rfloor^2\) 看起来很烦。

tjy 一个非常神奇的方法，我们化简为繁：\(n^2=\sum_{a=1}^n(2a-1)\)，于是：

\[\begin{aligned} &\sum_{i=1}^n\mu^2(i)\left\lfloor\sqrt{\frac{n}{i}}\right\rfloor^2\\ =&\sum_{i=1}^n\mu^2(i)\sum_{a^2i\leq n}(2a-1)\\ =&\sum_{a=1}^{\sqrt n}(2a-1)\sum_{a^2i\leq n}\mu^2(i) \end{aligned} \]

有一个把 \(\mu^2(n)\) 转化为 \(\mu(n)\) 的式子：

\[\mu^2(n)=\sum_{d^2|n}\mu(d) \]

可以这么理解，我们考虑 \(\mu^2(n)\) 和 \(\mu(n)\) 在某个质数 \(p\) 的幂处的取值：\(\mu^2(n):\{1,1,0,\cdots\}\)，\(\mu(n):\{1,-1,0,\cdots,\}\)，容易发现 \(\mu^2(n)\) 在 \(p^k\) 处的取值可以等于 \(\mu(n)\) 在 \(p^0,\cdots,p^{\lfloor\frac{k}{2}\rfloor}\) 处取值的和。

用贝尔级数来看，\(\mu^2_p(x)=1+x=\frac{1-x^2}{1-x}=\mu_p(x^2)*1_p\)。

（这条式子还可以使得 \(\mu^2(n)\) 的前缀和在 \(O(\sqrt n)\) 的时间内求出）

然后带进去就可以互推了：

\[\begin{aligned} &\sum_{a=1}^{\sqrt n}(2a-1)\sum_{i=1}^{\lfloor\frac{n}{a^2}\rfloor}\mu^2(i)\\ =&\sum_{a=1}^{\sqrt n}(2a-1)\sum_{i=1}^{\lfloor\frac{n}{a^2}\rfloor}\sum_{d^2|i}\mu(d)\\ =&\sum_{d=1}^{\sqrt n}\mu(d)\sum_{a=1}^{\sqrt n}(2a-1)\lfloor\frac{n}{a^2d^2}\rfloor\\ =&\sum_{T=1}^{\sqrt n}\lfloor\frac{n}{T^2}\rfloor\sum_{d|T}\mu(d)(\frac{2T}{d}-1)\\ =&\sum_{T=1}^{\sqrt n}\lfloor\frac{n}{T^2}\rfloor(2\cdot \varphi (T)-[T=1])\\ =&2\left(\sum_{T=1}^{\sqrt n}\lfloor\frac{n}{T^2}\rfloor\varphi(T)\right)-n \end{aligned} \]

总的来说，感觉转化过程中比较神奇的两步：（虽然感觉有点巧）

• \(n^2\) 的转化。
• \(\mu^2(n)\) 的转化。

代码如下：

#include<bits/stdc++.h>

#define N 10000010
#define ll long long

using namespace std;

inline int read()
{
	int x=0,f=1;
	char ch=getchar();
	while(ch<'0'||ch>'9')
	{
		if(ch=='-') f=-1;
		ch=getchar();
	}
	while(ch>='0'&&ch<='9')
	{
		x=(x<<1)+(x<<3)+(ch^'0');
		ch=getchar();
	}
	return x*f;
}

ll n;
int sn;
int cnt,prime[N],phi[N];
bool notprime[N];

void init()
{
	phi[1]=1;
	for(int i=2;i<=sn;i++)
	{
		if(!notprime[i]) prime[++cnt]=i,phi[i]=i-1;
		for(int j=1,v;j<=cnt&&(v=i*prime[j])<=sn;j++)
		{
			notprime[v]=1;
			if(!(i%prime[j]))
			{
				phi[v]=phi[i]*prime[j];
				break;
			}
			phi[v]=phi[i]*phi[prime[j]];
		}
	}
}

int main()
{
	cin>>n;
	sn=sqrt(n);
	init();
	ll ans=0;
	for(int i=1;i<=sn;i++)
		ans+=1ll*phi[i]*(n/(1ll*i*i));
	printf("%lld\n",ans*2ll-n);
	return 0;
}