特定のプロパティを持つ整数を求める - プロジェクト オイラー問題 221
https://stackoverflow.com/questions/386862
-
23-08-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
Russian質問
私は最近 Project Euler にとてもハマっていて、次のことに挑戦しています。 これ 次に一つ！私はそれについて分析を開始し、すでに問題を大幅に軽減しました。私の仕事は次のとおりです。
A = pqr および
1/a = 1/p + 1/q + 1/r so pqr/a = pq + pr + qr
そして、最初の方程式により、次のようになります。
pq+pr+qr = 1
正確に2つの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とbの私たちの値は次のとおりです。
A = abc = ab(ab-1)/(a+b)
計算が多すぎたらごめんなさい!しかし、今私たちが対処しなければならないのは、1 つの条件と 2 つの方程式だけです。ここで、ab = 1 (mod a+b) で ab(ab-1)/(a+b) として書かれた 150,000 番目に小さい整数を見つける必要があるので、理想的には A が何であるか (a, b) を検索したいと思います。できるだけ小さいもの。
簡単にするために、a < b と仮定しましたが、gcd(a, b) = 1 であることにも気付きました。
私の最初の実装は簡単で、150,000 個のソリューションも十分な速度で検出できました。しかし、15万件を見つけるには時間がかかりすぎます。 最小 ソリューション。とにかくコードは次のとおりです。
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 の異なる値が解を生み出すかどうかを確認することでした。そのコードは複雑で、高速ではありますが、十分な速さではありません...
だから完全にアイデアが足りない！何かアイデアがある人はいますか?
解決
プロジェクトオイラーの問題の多くと同じように、トリックはまっすぐ進むより何かに強引なソリューションを減らす技術を見つけることです。
A = pqr and
1/A = 1/p + 1/q + 1/r
だから、
pq + qr + rp = 1 or -r = (pq - 1)/(p + q)
一般性を失うことなく、0
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)
だから我々は、pを反復処理する必要がある、とAを計算し、kは唯一のpの約数の二乗プラスされている値を取ることができる場所1。
我々は必要な第十五万アレクサンドリア整数を見つけたときに、のセットに各溶液を加え、我々は停止することができます。
他のヒント
この中国の残りの部分についての記事、高速な実装、あなたを助けることができる：www.codeproject.com/KB/recipes/CRP.aspx
これは、ツールとライブラリのためのより多くのリンクです。
ツールます：
マキシマ http://maxima.sourceforge.net/する 最大値は、微分、積分、テイラー級数、ラプラス変換、常微分方程式、線形方程式のシステム、多項式、およびセット、リスト、ベクトル、行列、テンソルとを含むシンボリックと数値表現を操作するためのシステムです。最大値は、正確な画分を、任意の精度の整数、および可変精度浮動小数点数を使用して高精度の数値結果が得られます。 Maximaは、2次元および3次元での関数やデータをプロットすることができます。
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 インタラクティブな数式エディタとステップバイステップのソルバー。
ライブラリます：
アルマジロC ++ライブラリ http://arma.sourceforge.net/する アルマジロC ++ライブラリは、線形代数演算（行列及びベクトル数学）簡単で使いやすいインターフェースを持ちながらするための効率的な基盤を提供することを目的とします。
ブリッツ++ http://www.oonumerics.org/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 ++で実装有限要素法ライブラリです。これは、流体 - 構造相互作用のシミュレーションに適しており、それは数学的に高次の四面体要素に基づいています。
わかりました。ここで、中国剰余定理の解法をさらに試してみます。そうでない限り、a+b はどの素数 p の積にもなり得ないことがわかります。 p = 1 (mod 4). 。これにより、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)
しかし、これではあまり改善が見られません…。
