9

我最近对 ​​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 的 (a, b)尽可能小。

为方便起见，我假设 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 的不同值是否产生解决方案。该代码很复杂，虽然速度更快，但还不够快......

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

4

3 回答 3

4

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

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

工具：

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

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

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

mathstudio mathstudio.sourceforge.net 交互式方程编辑器和逐步求解器。

图书馆：

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

Blitz++ http://www.oonumerics.org/blitz/ Blitz++ 是一个用于科学计算的 C++ 类库

BigInteger C# http://msdn.microsoft.com/pt-br/magazine/cc163441.aspx

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

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

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

于 2008-12-22T18:24:09.967 回答
3

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

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 个亚历山大整数时停止。

于 2009-01-01T12:46:20.967 回答
0

好的。这是我的中国剩余定理解决方案的更多内容。事实证明，a+b 不能是任何素数 p 的乘积，除非p = 1 (mod 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)

但是，我看不到有太大的改善...

于 2008-12-23T12:46:07.963 回答