【数学】【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~ imes c_i~+~1))。根据乘法原理，总方案数为

[prod_{i = 1}^{k}~(2~ imes~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)