【数学】【P5150】 生日礼物

Description

给定 $n$，求

$$\sum_{i}~\sum_j~[lcm(i,j)~=~n]$$

input

一行一个整数代表 $n$

Output

一行一个整数代表答案

Hint

$1~\leq~n~\leq~10^{16}$

Solution

一开始看到这个形式以为是反演，然后看到数据范围就自闭了……

然后发现这是个唯一分解定理题……

吐槽一下标算 $O(\sqrt{n})$ 暴力卡常 范围出1e16也太[数据删除]了吧（大雾

然后用py写了一发和标算差不多的暴力，惨遭卡常

所以这里来提供一种 $O(\sqrt[3]{n})$ 的方法！

其实就是讨论里 @mrsrz 神仙的第一种踩标算做法

设

$$n~=~\prod_{i = 1}^{k} p_i^{c_i}$$

$$x~=~\prod_{i = 1}^k p_i^{d_i}$$

$$y~=~\prod_{i = 1}^{k}~p_{i}^{e_i}$$

其中 $p$ 为质数。

则显然有

$$lcm(x, y)~=~n~\Leftrightarrow~\forall~i~\in~[1,k],c_i~=~\max(d_i,e_i)$$

我们考虑固定 $x$ 第 $i$ 位指数即 $d_i~=~c_i$，则 $e_i$ 选 $[0,c_i]$ 都是合法的，共 $c_i~+~1$ 中选法。将 $x,y$ 反过来同样成立。但是注意固定 $d_i~=~c_i$ 时令 $e_i~=~c_i$ 的选法和反过来是一样的，于是要把这个方案扣除 $1$。

所以对于第 $i$ 个质因子的方案数为 $(2~\times c_i~+~1)$。根据乘法原理，总方案数为

$$\prod_{i = 1}^{k}~(2~\times~c_i~+~1)$$

于是 $O(\sqrt{n})$ 分解一下，发现py被卡常了。我们考虑一种更优秀的做法：

我们在分解质因数时，分解到 $\sqrt[3]{n}$，即当 $i^3~>~n$ 时停止。考虑现在 $n$ 除掉已经筛出的质因子后剩下的值共有如下几种情况：

剩下 $1$：这种情况对答案无贡献，无需理会

剩下的数是一个质数：显然这个剩下的数是 $n$ 的最后一个质因子，并且指数显然为 $1$，于是直接将答案乘 $3$ 即可

考虑除去这两种情况外，剩下的数只能是两个质数的积，而不可能是更多质数的积。

证明上，可以设剩下的最小的质数是 $p$，则有 $p~>~\sqrt[3]{n}$，假设是 $k$ 个质数的乘积，那么显然有剩下的数字 $dn~\geq~p^k$。由于 $p^3~>~(\sqrt[3]{n})^3~=~n$，$dn~\leq~n$，则在 $k~\geq~3$ 时产生矛盾，于是 $k~\leq~2$。

再分两种情况：

剩下的数是一个质数的平方：直接将答案乘 $5$ 即可

否则一定是两个质数相乘。考虑每个质数贡献 $3$，所以将答案乘上 $9$ 即可。

考虑如何快速判断剩下的数字是一个质数：直接进行米勒拉宾质数判定，时间复杂度 $O(\log n)$。

考虑不损失精度的判断一个数是一个完全平方数：直接进行二分开方，时间复杂度 $O(\log n)$

于是总时间复杂度 $O(\sqrt[3]{n})$，踩爆标算

Code

def mpow(x, y, p):
	_ret, _temp = 1, x
	while y:
		if y & 1:
			_ret = _ret * _temp % p
		_temp = _temp * _temp % p
		y >>= 1
	return _ret

def ML(x, n):
	if n == x: return 1
	sn = n - 1
	s, d = n - 1, 0
	while not (s & 1):
		s >>= 1
		d += 1
	t = mpow(x, s, n)
	if t == 1 or t == -1: return 1
	for i in range(d):
		if t == sn: return 1
		t = t * t % n
	return 0

def IsPrime(x):
	if not (x & 1):
		if x == 2: return 1
		else: return 0
	elif not ML(2, x): return 0
	elif not ML(7, x): return 0
	elif not ML(61, x): return 0
	else: return 1

def IsPow(x):
	l, r, mid, = 1, x, 0
	while l <= r:
		mid = (l + r) >> 1
		k = mid * mid
		if k < x: l = mid + 1
		elif k == x: return 1
		else: r = mid - 1
	return 0

n = int(input())

ans, i, dn = 1, 2, n

while (i * i * i) <= n:
	if (dn % i) == 0:
		cnt = 0
		while (dn % i) == 0:
			dn //= i
			cnt += 1
		ans *= (cnt << 1) + 1
	i += 1
if dn != 1:
	if IsPrime(dn):
		ans *= 3
	elif IsPow(dn):
		ans *= 5
	else:
		ans *= 9
		
print(ans)