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 の 2 つだけが負でなければならないため、式を次のように単純化できます。

ab = ac+bc+1 の abc

c について解くと、次のようになります。

ab-1 = (a+b)c

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

これは、次の a と b を見つける必要があることを意味します。

ab = 1 (mod a+b)

そして、これらの a と b の A 値は次のとおりです。

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

それが多くの数学である場合は申し訳ありません！しかし、ここで取り扱わなければならないのは、1 つの条件と 2 つの方程式だけです。ab = 1（mod a + b）でab（ab-1）/（a + b）として記述された150,000番目に小さい整数を見つける必要があるため、理想的には、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 の異なる値が解をもたらすかどうかを確認することでした。そのコードは複雑で、高速ではありますが、十分に高速ではありません...

だから私は絶対にアイデアがありません！誰でもアイデアはありますか？

4

3 に答える 3

4

中国語の残り、迅速な実装に関するこの記事は、あなたを助けることができます: www.codeproject.com/KB/recipes/CRP.aspx

これは、ツールとライブラリへのその他のリンクです。

ツール:

Maxima http://maxima.sourceforge.net/ Maxima は、微分、積分、テイラー級数、ラプラス変換、常微分方程式、線形方程式系、多項式、集合、リストなど、記号式および数値式を操作するためのシステムです。 、ベクトル、行列、およびテンソル。Maxima は、正確な分数、任意精度の整数、および可変精度の浮動小数点数を使用して、高精度の数値結果を生成します。Maxima は関数とデータを 2 次元と 3 次元でプロットできます。

Mathomatic http://mathomatic.org/math/ Mathomatic は、無料で移植可能な汎用 CAS (Computer Algebra System) および計算ソフトウェアであり、方程式を記号的に解き、単純化し、結合し、比較し、複素数および多項式演算を実行できます。など。微積分を行い、非常に使いやすいです。

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 プロジェクトのホームページへようこそ。このプロジェクトの目的は、最も使いやすい任意精度の 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

Project Euler の問題の多くと同様に、秘訣は、ブルート フォース ソリューションをより単純なものに還元する手法を見つけることです。

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 の 2 乗に 1 を加えた値の約数である値のみを取ることができます。

各解をセットに追加すると、必要な 150000 番目のアレクサンドリア整数が見つかったときに停止できます。

于 2009-01-01T12:46:20.967 に答える
0

わかった。これが私の中国剰余定理の解決策でもう少し遊んでいます。p = 1（mod 4）でない限り、a+bは素数pの積にはなり得ないことがわかります。これにより、2、5、13、17、29、37などの素数の倍数であるa + bをチェックするだけでよいため、計算が高速になります。

したがって、可能な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 に答える