编写程序找费马“二平方”素数——“滚学球”式所使用的素数

费马二平方定理指出，一个正整数可以表示为两个平方数之和的充分必要条件是，其所有形如 $4k+3$ 的质因子的幂次均为偶数。这里，“$k$”表示任意非负整数。

“滚学球”式是指以下形式的表达式：

$$
n=a^2+b^2=c^2+d^2
$$

其中，$a,b,c,d$ 均为正整数。这个式子可以用来表示一个正整数 $n$ 能否表示为两个平方数之和的两种不同方式。

为了找到费马“二平方”素数，我们可以按照以下步骤：

1. 枚举所有形如 $4k+1$ 的正整数 $n$。
2. 对于每个 $n$，枚举所有满足 $a^2+b^2=n$ 的正整数 $a,b$，并计算 $c^2+d^2=n$ 的正整数 $c,d$。
3. 判断 $n$ 是否为素数，如果是素数，且所有形如 $4k+3$ 的质因子的幂次均为偶数，则输出 $n$。

下面是用 Python 实现的代码：

import math

def is_prime(n):
    if n < 2:
        return False
    for i in range(2, int(math.sqrt(n)) + 1):
        if n % i == 0:
            return False
    return True

for n in range(1, 100000):
    if n % 4 == 1:
        for a in range(1, int(math.sqrt(n)) + 1):
            b = int(math.sqrt(n - a * a))
            if a * a + b * b == n:
                for c in range(1, int(math.sqrt(n)) + 1):
                    d = int(math.sqrt(n - c * c))
                    if c * c + d * d == n:
                        if is_prime(n):
                            odd_factors = 0
                            for factor in range(3, int(math.sqrt(n)) + 1, 2):
                                if n % factor == 0:
                                    exponent = 0
                                    while n % factor == 0:
                                        n //= factor
                                        exponent += 1
                                    if factor % 4 == 3 and exponent % 2 == 1:
                                        odd_factors += 1
                            if n > 2 and n % 4 == 3:
                                odd_factors += 1
                            if odd_factors == 0:
                                print(n)
                                break

这个程序的运行时间比较长，因为需要枚举很多数。在这个程序中，我们首先定义了一个函数 `is_prime`，用于判断一个数是否为素数。然后，我们使用两个嵌套的循环枚举 $a,b$ 和 $c,d$，并判断 $n$ 是否为素数。如果 $n$ 是素数，我们就检查它的所有质因子，判断它们的幂次是否为偶数。如果 $n$ 的所有形如 $4k+3$ 的质因子的幂次均为偶数，则输出 $n$。