108. 丢番图倒数一(Diophantine reciprocals I)

在以下方程中，\(x,y,n\)都是正整数：

$$ \dfrac{1}{x} + \dfrac{1}{y} = \dfrac{1}{n} $$

对于\(n=4\)的情况，恰好有三种不同的解：

$$ \dfrac{1}{5} + \dfrac{1}{20} = \dfrac{1}{4}， \dfrac{1}{6} + \dfrac{1}{12} = \dfrac{1}{4}， \dfrac{1}{8} + \dfrac{1}{8} = \dfrac{1}{4} $$

注：这个问题是\(problem\ 110\)的简化版本，强烈建议先解决这个问题。

分析：既然题目已经说了这个问题是第一百一十题的简化版本，所以我们把这两个问题放在一起考虑。这两个问题的区别是数据规模不同，这道题要求的是正整数解个数超过一千时\(n\)的最小值，而第一百一十题是求解的个数超过四百万时\(n\)的最小值。问题规模的提升使得第一百一十题不可能使用任何暴力方法在规定的时间内解决，必须要使用某些巧妙的方法。显然，这种巧妙的方法应该也是可以用来解决108题的。

我们首先分析下这个题目，解题的一个核心是如何求这种丢番图方程的正整数解的个数。首先把原方程左右同时乘以\(xyn\)，则有：

两边同时加上\(n^2\)，则有：

观察最右边的式子，要想\(x,y\)是正整数，则\((x-n),(y-n)\)必须是\(n^2\)的因子，这里令\(n^2\)的因子数量为\(d(n^2)\)。题目要求以上方程的不重复的正整数解的个数，因为上面方程中\(x,y\)是可以互换的，同时我们知道完全平方数的因子个数是一个奇数，所以上面方程的不重复的正整数解个数为\((d(n^2)+1)/2\)。如当\(n=4\)时，\((d(16)+1)/2=(5+1)/2=3\)。因此题目的要求可以表述为\((d(n^2)+1)/2>1000\Rightarrow d(n^2)>1999\)，即我们要求满足这个条件的最小的\(n\)。

根据初等数论的知识，我们可以从数\(n\)的质因数分解形式计算其因子个数：

则对于\(n^2\)，我们有：

观察上式可以发现，一个数的因子个数与其分解成那些质数无关，而只与这些质数的指数有关。现在假设我们已知\(d(n^2)=Q\)，则我们可以对\(Q\)进行质因数分解：

根据以上等式，我们可以推出已知因子数为\(Q\)的情况下数\(n\)的最小值，因为\(Q\)的质因数分解已经表明原数\(n\)的质因数分解中指数为\((p_j-1)/2\)的质因数有\(e_j\)个，如果我们把最大的指数安排给最小的质因数，就可以得到因子数为\(Q\)的\(n\)的最小值。

上面的推理确实有些复杂，我们看一个例子，假设我们知道\(d(n^2)=2025\)，要求\(n\)的最小值，则我们可以对\(2025\)进行质因数分解，得到\(2025=3^4\cdot5^2\)，也即：

显然，我们可以做如下对应：

因此可以知道\(a_1=a_2=a_3=a_4=1,a_5=a_6=2\)，也就是说指数为\((3-1)/2=1\)的有\(4\)个，而指数为\((5-1)/2=2\)的有\(2\)个。为了得到数\(n\)的最小值，我们必须把大的指数尽量分配给小的质因数，因此我们把指数\(2\)分配给质因数\(2,3\)，把指数\(1\)分配给质因数\(5,7,11,13\)，则：

具体到题目的要求来看，要求\(d(n^2)>1999\)时的\(n_{min}\)，我们一是需要使得\(d(n^2)\)足够小，同时其质因数分解应只包含\(3,5,7\)三个素数，这样才能使得原数\(n\)质因数分解的指数足够小，从而使得\(n\)足够小。具体的实现可见代码，这里不再赘述：

# time cost = 151 µs ± 2.16 µs

from sympy import factorint
from functools import reduce
from operator import mul,add

def prime_product(factors):
    primes = [2,3,5,7,11,13,17,19,23,29,31,37,41,43,47]
    powers = sorted(reduce(add,[[(k-1)//2]*v for k,v in factors.items()]),reverse=True)
    bases = primes[:len(powers)]
    res = reduce(mul,[b**p for b,p in zip(bases,powers)])
    return res

def main(N=1000):
    n = 2*N - 1
    primes = [2,3,5,7,11,13,17,19,23,29,31,37,41,43,47]
    while True:
        factors = factorint(n)
        if all(x in [3,5,7] for x in factors.keys()):
            return prime_product(factors)
        n += 2