不定方程式 X3 + Y3 = Z には整数解が存在します。たとえば、Z = 2 のときの解は 13 + 13 の一通りしかありませんが、Z の値によっては複数の解が存在します。そして、n 通りの解が存在する Z の中で、最小の値を「タクシー数 (taxicab number)」といい、n 番目のタクシー数を Ta(n) と書きます。Ta(1) と Ta(2) を示します。
Ta(1) = 2 = 13 + 13
Ta(2) = 1729 = 13 + 123
= 93 + 103
参考 URL 1 によると、Ta(n) はすべての整数 n に対して存在することが証明されていて、今までに Ta(1) から Ta(6) までのタクシー数が知られているそうです。
今回は問題を簡単にして、次式を満たす整数 a, b, c, d を列挙することにします。
a3 + b3 = c3 + d3
使用するプログラミング言語は Python3 (PyPy3) です。
最近のパソコンは高性能なので、単純な四重ループでも簡単に解くことができます。次のリストを見てください。
リスト : タクシー数
def taxi(n):
xs = []
for a in range(1, n + 1):
for b in range(a, n + 1):
for c in range(a + 1, n + 1):
for d in range(c, b):
e = a ** 3 + b ** 3
if c ** 3 + d ** 3 == e:
xs.append((e, a, b, c, d))
xs.sort(key = lambda x: x[0])
return xs
引数 n が整数の上限値を表します。重複解を削除するため、以下の条件を設定してます。
1. a <= b # a と b を交換した式を削除 2. c <= d # c と d を交換した式を削除 3. a < c # (a, b) と (c, d) を交換した式を削除 4. d < b # 不要な探索を削除
条件 3 で c は a よりも大きくなるので、d は必ず b よりも小さくなります。これを条件 4 で表しています。あとはとくに難しいところはないと思います。
それでは実行してみましょう。
>>>> for x in taxi(50): print(x) (1729, 1, 12, 9, 10) (4104, 2, 16, 9, 15) (13832, 2, 24, 18, 20) (20683, 10, 27, 19, 24) (32832, 4, 32, 18, 30) (39312, 2, 34, 15, 33) (40033, 9, 34, 16, 33) (46683, 3, 36, 27, 30) (64232, 17, 39, 26, 36) (65728, 12, 40, 31, 33) (110656, 4, 48, 36, 40) (110808, 6, 48, 27, 45) >>>> s = time.time(); print(len(taxi(500))); time.time() - s 566 5.301697492599487 実行環境 : PyPy3 (ver 7.3.1), Ubunts 20.04 (WSL1), Intel Core i5-6200U 2.30GHz
整数の範囲 50 から 500 に増やすと、実行時間は遅くなります。なお、taxi(500) の結果には以下に示す重複解が含まれていることに注意してください。
(87539319, 167, 436, 228, 423) (87539319, 167, 436, 255, 414) (87539319, 228, 423, 255, 414) (119824488, 11, 493, 90, 492) (119824488, 11, 493, 346, 428) (119824488, 90, 492, 346, 428)
これらは条件 a3 + b3 = c3 + d3 = e3 + f3 を満たしています。この条件を満たす最小の数が 87,539,319 = Ta(3) になります。この場合、数値 Z の個数は 566 - 4 = 562 になります。
次のように、Python の辞書を使って二重ループにすると速くなります。
リスト : タクシー数 (2)
def taxi2(n):
table = {}
for a in range(1, n + 1):
for b in range(a, n + 1):
x = a ** 3 + b ** 3
if x in table:
table[x].append((a, b))
else:
table[x] = [(a, b)]
xs = [(k, v) for k, v in table.items() if len(v) > 1]
xs.sort(key = lambda x: x[0])
return xs
>>>> s = time.time(); print(len(taxi2(500))); time.time() - s 562 0.082763671875
60 倍ちょっと高速化することができました。ご参考までに a3 + b3 = c3 + d3 = e3 + f3 の解を求めてみました。
>>>> s = time.time(); print(taxi3(2000)); time.time() - s [(87539319, [(167, 436), (228, 423), (255, 414)]), (119824488, [(11, 493), (90, 492), (346, 428)]), (143604279, [(111, 522), (359, 460), (408, 423)]), (175959000, [(70, 560), (198, 552), (315, 525)]), (327763000, [(300, 670), (339, 661), (510, 580)]), (700314552, [(334, 872), (456, 846), (510, 828)]), (804360375, [(15, 930), (198, 927), (295, 920)]), (958595904, [(22, 986), (180, 984), (692, 856)]), (1148834232, [(222, 1044), (718, 920), (816, 846)]), (1407672000, [(140, 1120), (396, 1104), (630, 1050)]), (1840667192, [(225, 1223), (372, 1214), (681, 1151)]), (1915865217, [(9, 1242), (484, 1217), (969, 1002)]), (2363561613, [(501, 1308), (684, 1269), (765, 1242)]), (2622104000, [(600, 1340), (678, 1322), (1020, 1160)]), (3080802816, [(81, 1455), (456, 1440), (904, 1328)]), (3235261176, [(33, 1479), (270, 1476), (1038, 1284)]), (3499524728, [(116, 1518), (350, 1512), (1169, 1239)]), (3623721192, [(348, 1530), (761, 1471), (1098, 1320)]), (3877315533, [(333, 1566), (1077, 1380), (1224, 1269)]), (4750893000, [(210, 1680), (594, 1656), (945, 1575)]), (5544709352, [(207, 1769), (842, 1704), (1076, 1626)]), (5602516416, [(668, 1744), (912, 1692), (1020, 1656)]), (6434883000, [(30, 1860), (396, 1854), (590, 1840)]), (7668767232, [(44, 1972), (360, 1968), (1384, 1712)])] 1.6516444683074951
taxi3() は内包表記の条件を len(v) > 2 に変更しただけです。整数の上限値を増やすと、もっと時間がかかるようになります。Ta(4) は大きな値になるので、このような力任せのプログラムだと M.Hiroi の実行環境では Ta(3) くらいが限界のように思います。タクシー数を求めるのは難しい問題だと実感しました。