540. 本原毕达哥拉斯三元组计数(Counting primitive Pythagorean triples)

毕达哥拉斯三元组包含有三个整数\(a,b,c\)且满足等式\(a^2+b^2=c^2\),当\(a,b,c\)互素时,这个三元组被称为本原的。记\(P(n)\)为满足\(a<b<c \le n\)的本原毕达哥拉斯三元组的数目。例如\(P(20)=3\),这三个三元组是:\((3,4,5),(5,12,13)\)\((8,15,17)\)。已知\(P(10^6)=159139\),求\(P(3141592653589793)\)

分析:针对这道题我虽然尝试了很多方法,最终也给出了正确答案,但是程序的运行时间达到了36分钟,显然不符合欧拉计划网站的标准,所以不算完全解决了这个问题。但是也不妨记录一下思路,等到将来有更好的方法再做改进。题目要求满足斜边长度不超过\(n\)的本原毕达哥拉斯三元组的个数,首先我们需要知道这些三元组如何产生出来。我们可以使用在以前的题目中就介绍过的欧几里德公式,即当\(gcd(x,y)=1,x>y\ge1\)\(x,y\)不同为奇数时,本原毕达哥拉斯三元组可以使用以下公式产生:

$$ a=x^2-y^2,b=2xy,c=x^2+y^2 $$

因此我们的问题可以转化成求满足以下条件的\((x,y)\)对的个数:

我们可以用图示让上述要求看得更清楚,针对第一个条件即\(x^2+y^2\le n=(\sqrt{n})^2\),即\(x,y\)要在圆心在原点,半径为\(\sqrt{n}\)的圆内及圆上,同时要求\(x>y\),则这些点还要在\(y=x\)的直线下方。在下面的图示中,我们要求三元组的斜边长度不超过100,即\(x,y\)都在半径为10原圆内或圆上及位于绿色直线下方,即位于圆\(x^2+y^2=10\),直线\(y=x\)以及正横轴三者围住的区域,而\((x,y)\)正对应这个区域内的格点(即横纵坐标均为整数的点)。

需要注意的是,橙色直线\(BC\)把这个区域分成了左右两个部分,在左边绿色直线位于圆周下方,则起作用的约束条件是\(x>y\);而在右边绿色直线位于圆周上方,这时起作用的约束条件是\(x^2+y^2\le n\),即\(y\le \sqrt{n-x^2}\)

在以上区域中,我们找到了16个格点,也即斜边长度不超过100的本原毕达哥拉斯三元组共有16个。比如当\(x=5\)时,我们首先看\((5,1)\)这个格点,因为条件\(x,y\)不能同时为奇数,所以这个格点排除,同理可以把\((5,3)\)排除;再看\((5,2)\)这个格点,满足\(gcd(5,2)=1\),则一个符合要求的格点,同理\((5,4)\)也符合。再往上就不满足\(x>y\)的要求了,所以格点\((5,5)\)以及的格点都排除点。基于以上思路,我们可以把满足要求的格点分为四种情况:

欧拉函数\(\varphi(x)\)可以比较容易求得,但考虑到计算中需要大量的欧拉函数值,所以可以用筛法预先计算出这次值,这样比每次需要时再做计算效率要更高一些,具体实现可以参见代码中的\(euler\_phi\_sieve()\)函数。而对于广义欧拉函数,使用莫比乌斯函数,我们也可以计算出来:

$$ \phi ( x,y) =\sum _{d|x,d\leq y} \mu ( d)\left\lfloor \frac{y}{d}\right\rfloor $$

为加快计算,我们可以使用\(spymp\)模块中的\(mobiusrange()\)函数预先计算出所需的莫比乌斯函数值。完整代码如下:

# time cost = 35min 37s

from sympy import sieve,divisors

def euler_phi_sieve(n):
    arr = list(range(n+1))
    for i in range(2,n+1):
        if arr[i] == i:
            for j in range(i,n+1,i):
                arr[j] -= arr[j]//i
    return arr

def phi(m,n,premob):
    res = 0
    for i in [x for x in divisors(m) if x<=n]:
        res += (premob[i-1]*(n//i))
    return res

def main(N=3141592653589793):
    res = 0
    u = int((N/2)**0.5)+1
    prephi = euler_phi_sieve(u)
    premob = list(sieve.mobiusrange(1,u))
    for m in range(2,u):
        if m % 2 == 0:
            res += prephi[m]
        else:
            res += (prephi[m]//2)
    for m in range(u,int(N**0.5)+1):
        v = int((N-m**2)**0.5)
        if m % 2 == 0:
            res += phi(m,v,premob)
        else:
            res += phi(m,v//2,premob)
    return res