令 n 为平方数。使用Python，我们如何有效地计算自然数y直到极限l，使得n+y^2再次成为平方数？

Question

使用Python，我想实现一个函数，它接受一个自然数n作为输入并输出一个自然数列表[y1, y2, y3, ...]，使得n + y1*y1和n + y2*y2等等n + y3*y3又是一个正方形。

到目前为止我尝试的是y使用以下函数获取一个值：

def find_square(n:int) -> tuple[int, int]:
    if n%2 == 1:
        y = (n-1)//2
        x = n+y*y
        return (y,x)
    return None

它工作正常，例如。find_square(13689)给了我一个正确的解决方案y=6844。y如果有一种算法能够产生所有可能的- 值（例如y=44或 ），那就太好了y=156。

Answer 1

最简单的慢速方法当然是对于给定的 N 只是迭代所有可能的 Y 并检查是否N + Y^2是平方的。

\n

但是有一种更快的方法，使用整数分解技术：

\n

让我们注意到，要求解方程N + Y^2 = X^2，即找到给定固定整数 N 的所有整数对 (X, Y)，我们可以重写该方程，N = X^2 - Y^2 = (X + Y) * (X - Y)该方程源自著名的平方差公式。

\n

现在让我们将两个因子重命名为 A, B 即N = (X + Y) * (X - Y) = A * B，这意味着X = (A + B) / 2和Y = (A - B) / 2。

\n

请注意，A 和 B 应该具有相同的奇数，要么都是奇数，要么都是偶数，否则在上面的最后一个公式中我们不能将其全部除以 2。

\n

我们将把 N 分解为所有可能的具有相同奇数的两个因子 (A, B) 对。为了在下面的代码中进行快速分解，我使用了实现简单但速度相当快的算法Pollard Rho，还需要两种额外的算法作为 Pollard Rho 的助手，一个是费马素性测试（它允许快速检查数字是否可能是素数），第二个是是试除因式分解（它帮助 Pollard Rho 分解出小因子，这可能会导致 Pollard Rho 失败）。

\n

复合数的 Pollard Rho 的时间复杂度O(N^(1/4))即使对于 64 位数字也非常快。如果需要搜索更大的空间，可以选择任何更快的分解算法。我的快速算法时间主要取决于因式分解的速度，算法的其余部分非常快，只需使用简单公式进行几次循环迭代即可。

\n

O(N^(1/8))如果你的 N 本身就是一个平方（因此我们很容易知道它的根），那么 Pollard Rho 可以在短时间内更快地因式分解 N。即使对于 128 位数字，这也意味着非常短的时间，即 2^16 次操作，我希望您能用少于 128 位数字来解决您的任务。

\n

如果您想处理一系列可能的 N 值，那么分解它们的最快方法是使用类似于埃拉托斯特尼筛法的技术，使用一组素数，它允许计算某个范围内所有 N 个数字的所有因子。对于 Ns 范围的情况，使用埃拉托斯特尼筛法比用 Pollard Rho 对每个 N 进行因式分解要快得多。

\n

将 N 分解为对 (A, B) 后，我们根据 (A, B) 通过上述公式计算 (X, Y)。并输出结果Y作为快速算法的解。

\n

以下代码作为示例，是用纯Python实现的。当然可以使用Numba来加速，Numba通常可以提供30-200倍的加速，对于Python来说它可以达到与优化的C++相同的速度。但我认为这里主要是实现快速算法，之后可以轻松完成 Numba 优化。

\n

我在以下代码中添加了时间测量。虽然它是纯Python，但与常规的暴力方法相比，我的快速算法在 1 000 000 的限制下实现了8500倍的加速。

\n

您可以更改limit变量来调整搜索空间的数量，或num_tests更改变量来调整不同测试的数量。

\n

以下代码实现了两种解决方案 -find_fast()上面描述的快速解决方案加上非常小的暴力解决方案find_slow()，该解决方案非常慢，因为它扫描所有可能的候选者。这种缓慢的解决方案仅用于比较测试中的正确性和比较加速。

\n

下面的代码除了几个标准 Python 库模块外什么也没使用，没有使用任何外部模块。

\n

在线尝试一下！

\n

def find_slow(N):\n    import math\n    \n    def is_square(x):\n        root = int(math.sqrt(float(x)) + 0.5)\n        return root * root == x, root\n    \n    l = []\n    for y in range(N):\n        if is_square(N + y ** 2)[0]:\n            l.append(y)\n    \n    return l\n    \ndef find_fast(N):\n    import itertools, functools\n    \n    Prod = lambda it: functools.reduce(lambda a, b: a * b, it, 1)\n    \n    fs = factor(N)\n    mfs = {}\n    for e in fs:\n        mfs[e] = mfs.get(e, 0) + 1\n    fs = sorted(mfs.items())\n    del mfs\n    \n    Ys = set()\n    for take_a in itertools.product(*[\n        (range(v + 1) if k != 2 else range(1, v)) for k, v in fs]):\n        A = Prod([p ** t for (p, _), t in zip(fs, take_a)])\n        B = N // A\n        assert A * B == N, (N, A, B, take_a)\n        if A < B:\n            continue\n        X = (A + B) // 2\n        Y = (A - B) // 2\n        assert N + Y ** 2 == X ** 2, (N, A, B, X, Y)\n        Ys.add(Y)\n    \n    return sorted(Ys)\n\ndef trial_div_factor(n, limit = None):\n    # https://en.wikipedia.org/wiki/Trial_division\n    fs = []\n    while n & 1 == 0:\n        fs.append(2)\n        n >>= 1\n    all_checked = False\n    for d in range(3, (limit or n) + 1, 2):\n        if d * d > n:\n            all_checked = True\n            break\n        while True:\n            q, r = divmod(n, d)\n            if r != 0:\n                break\n            fs.append(d)\n            n = q\n    if n > 1 and all_checked:\n        fs.append(n)\n        n = 1\n    return fs, n\n\ndef fermat_prp(n, trials = 32):\n    # https://en.wikipedia.org/wiki/Fermat_primality_test\n    import random\n    if n <= 16:\n        return n in (2, 3, 5, 7, 11, 13)\n    for i in range(trials):\n        if pow(random.randint(2, n - 2), n - 1, n) != 1:\n            return False\n    return True\n\ndef pollard_rho_factor(n):\n    # https://en.wikipedia.org/wiki/Pollard%27s_rho_algorithm\n    import math, random\n    fs, n = trial_div_factor(n, 1 << 7)\n    if n <= 1:\n        return fs\n    if fermat_prp(n):\n        return sorted(fs + [n])\n    for itry in range(8):\n        failed = False\n        x = random.randint(2, n - 2)\n        for cycle in range(1, 1 << 60):\n            y = x\n            for i in range(1 << cycle):\n                x = (x * x + 1) % n\n                d = math.gcd(x - y, n)\n                if d == 1:\n                    continue\n                if d == n:\n                    failed = True\n                    break\n                return sorted(fs + pollard_rho_factor(d) + pollard_rho_factor(n // d))\n            if failed:\n                break\n    assert False, f\'Pollard Rho failed! n = {n}\'\n    \ndef factor(N):\n    import functools\n    Prod = lambda it: functools.reduce(lambda a, b: a * b, it, 1)\n    fs = pollard_rho_factor(N)\n    assert N == Prod(fs), (N, fs)\n    return sorted(fs)\n\ndef test():\n    import random, time\n    limit = 1 << 20\n    num_tests = 20\n    \n    t0, t1 = 0, 0\n    for i in range(num_tests):\n        if (round(i / num_tests * 1000)) % 100 == 0 or i + 1 >= num_tests:\n            print(f\'test {i}, \', end = \'\', flush = True)\n        \n        N = random.randrange(limit)\n        \n        tb = time.time()\n        r0 = find_slow(N)\n        t0 += time.time() - tb\n        \n        tb = time.time()\n        r1 = find_fast(N)\n        t1 += time.time() - tb\n        \n        assert r0 == r1, (N, r0, r1, t0, t1)\n        \n    print(f\'\\nTime slow {t0:.05f} sec, fast {t1:.05f} sec, speedup {round(t0 / max(1e-6, t1))} times\')\n\nif __name__ == \'__main__\':\n    test()\n

\n

输出：

\n

test 0, test 2, test 4, test 6, test 8, test 10, test 12, test 14, test 16, test 18, test 19,\nTime slow 26.28198 sec, fast 0.00301 sec, speedup 8732 times\n

\n