查找具有特定属性的整数 - 欧拉计划问题 221

Question

我最近对 ​​Project Euler 非常着迷，并且正在尝试这个 下一个！ 我已经开始对此进行一些分析，并且已经大大减少了问题的严重性。 这是我的工作：

A = pqr 且

1/A = 1/p + 1/q + 1/r 所以 pqr/A = pq + pr + qr

由于第一个方程：

pq+pr+qr = 1

因为 p、q 和 r 中正好有两个 为负数，我们可以简化 等式归结为：

abc，其中 ab = ac+bc+1

求解 c 我们得到：

ab-1 = (a+b)c

c = (ab-1)/(a+b)

<小时>

这意味着我们需要找到 a 和 b 其中：

ab = 1 (mod a+b)

然后我们的 A 值与 a 和 b 是：

A = abc = ab(ab-1)/(a+b)

抱歉，如果这涉及很多数学运算！ 但现在我们需要处理的只是一个条件和两个方程。 现在，由于我需要找到写为 ab(ab-1)/(a+b) 且 ab = 1 (mod a+b) 的第 150,000 个最小整数，理想情况下我想搜索 (a, b) ，其中 A 是尽可能小。

为了方便起见，我假设 a < b，我还注意到 gcd(a, b) = 1。

我的第一个实现非常简单，甚至足够快地找到 150,000 个解决方案。 然而，找到 150,000 个最小解决方案需要很长时间。 无论如何，这是代码：

n = 150000
seen = set()

a = 3
while len(seen) < n:
    for b in range(2, a):
        if (a*b)%(a+b) != 1: continue

        seen.add(a*b*(a*b-1)//(a+b))
        print(len(seen), (a, b), a*b*(a*b-1)//(a+b))

    a += 1

我的下一个想法是使用 Stern-Brocot 树，但这太慢了，无法找到解决方案。 我的最终算法是使用中国剩余定理来检查 a+b 的不同值是否产生解决方案。 该代码很复杂，虽然速度更快，但还不够快......

所以我完全没有想法！ 有人有什么想法吗？

Answer 1

与许多欧拉计划问题一样，诀窍是找到一种技术，将蛮力解决方案简化为更直接的解决方案：

A = pqr and 
1/A = 1/p + 1/q + 1/r

因此，

pq + qr + rp = 1  or  -r = (pq - 1)/(p + q)

不失一般性，0 << p＜ -q < -r

存在 k , 1 <= k <= p

-q = p + k
-r = (-p(p + k) – 1) / (p + -p – k)  = (p^2 + 1)/k + p

但是 r 是一个整数，所以 k 整除 p^2 + 1

pqr = p(p + q)((p^2 + 1)/k + p)

因此要计算 A 我们需要迭代 p，其中 k 只能取除数p 的平方加 1。

将每个解添加到一个集合中，当我们找到所需的第 150000 个亚历山大整数时，我们就可以停止。

Answer 2

这篇关于中文余数的文章，快速实现，可以帮助您：www.codeproject.com/KB/recipes/CRP.aspx

这是更多工具和库的链接：

工具：

Maxima
http://maxima.sourceforge.net/
Maxima 是一个用于操作符号和数值表达式的系统，包括微分、积分、泰勒级数、拉普拉斯变换、常微分方程、线性方程组、多项式以及集合、列表、向量、矩阵和张量。 Maxima 通过使用精确分数、任意精度整数和可变精度浮点数来生成高精度数值结果。 Maxima 可以绘制二维和三维的函数和数据。

数学的
http://mathomatic.org/math/
Mathomatic 是一款免费、可移植、通用的 CAS（计算机代数系统）和计算器软件，可以符号求解、化简、组合和比较方程，执行复数和多项式运算等。它可以进行一些微积分运算，非常容易使用。

科学实验室
www.scilab.org/download/index_download.php
Scilab 是一个类似于 Matlab 或 Simulink 的数值计算系统。 Scilab 包含数百个数学函数，并且可以交互地添加各种语言（例如 C 或 Fortran）的程序。

数学工作室
mathstudio.sourceforge.net
交互式方程编辑器和逐步求解器。

库：

犰狳 C++ 库
http://arma.sourceforge.net/
Armadillo C++ 库旨在为线性代数运算（矩阵和向量数学）提供有效的基础，同时拥有简单易用的界面。

闪电战++
http://www.oonumerics.org/blitz/
的 C++ 类库

Blitz++ 是一个用于科学计算BigInteger C#
http://msdn.microsoft.com/pt-br/magazine/cc163441。 aspx

libapmath
http://freshmeat.net/projects/libapmath
欢迎来到 APMath 项目的主页。 该项目的目标是实现一个任意精度的C++库，这是使用最方便的，这意味着所有操作都实现为运算符重载，命名大多与.

库马特
http://freshmeat.net/projects/libmat
MAT是一个C++数学模板类库。 使用该库进行各种矩阵运算、求多项式的根、求解方程等。该库仅包含 C++ 头文件，因此无需编译。

动画
http://www.yonsen.bz/animath/animath.html
Animath 是一个完全用 C++ 实现的有限元方法库。 它适用于流固耦合模拟，并且在数学上基于高阶四面体单元。

Answer 3

好的。 这是我的中国剩余定理解决方案的更多内容。 事实证明，a+b 不能是任何素数 p 的乘积，除非 p = 1（模 4）。 这允许更快的计算，因为我们只需要检查 a+b，它们是素数的倍数，例如 2、5、13、17、29、37...

可能的 a+b 值的示例：

[5, 8, 10, 13, 16, 17, 20, 25, 26, 29, 32, 34, 37, 40, 41, 50, 52, 53, 58, 61, 64, 65, 68, 73, 74, 80, 82, 85, 89, 97, 100]

所以这里是 使用中国余数定理的完整程序：

cachef = {}
def factors(n):
    if n in cachef: cachef[n]

    i = 2
    while i*i <= n:
        if n%i == 0:
            r = set([i])|factors(n//i)
            cachef[n] = r
            return r

        i += 1

    r = set([n])
    cachef[n] = r
    return r

cachet = {}
def table(n):
    if n == 2: return 1
    if n%4 != 1: return

    if n in cachet: return cachet[n]

    a1 = n-1
    for a in range(1, n//2+1):
        if (a*a)%n == a1:
            cachet[n] = a
            return a

cacheg = {}
def extended(a, b):
    if a%b == 0:
        return (0, 1)
    else:
        if (a, b) in cacheg: return cacheg[(a, b)]

        x, y = extended(b, a%b)
        x, y = y, x-y*(a//b)

        cacheg[(a, b)] = (x, y)
        return (x, y)

def go(n):
    f = [a for a in factors(n)]
    m = [table(a) for a in f]

    N = 1
    for a in f: N *= a

    x = 0
    for i in range(len(f)):
        if not m[i]: return 0

        s, t = extended(f[i], N//f[i])
        x += t*m[i]*N//f[i]

    x %= N
    a = x
    while a < n:
        b = n-a
        if (a*b-1)%(a+b) == 0: return a*b*(a*b-1)//(a+b)
        a += N




li = [5, 8, 10, 13, 16, 17, 20, 25, 26, 29, 32, 34, 37, 40, 41, 50, 52, 53, 58, 61, 64, 65, 68, 73, 74, 80, 82, 85, 89, 97, 100]

r = set([6])
find = 6

for a in li:
    g = go(a)
    if g:
        r.add(g)
        #print(g)
    else:
        pass#print(a)


r = list(r)
r.sort()

print(r)
print(len(r), 'with', len(li), 'iterations')

这更好，但我希望进一步改进它（例如 a+b = 2^n 似乎永远不是解决方案）。

我也开始考虑基本的替代品，例如：

a = u+1 且 b = v+1

ab = 1 (mod a+b)

uv+u+v = 0 (mod u+v+2)

但是，我看不出有多大改进......