Taxicab number $\operatorname{Ta}(4)$

今回は「 $x^3+y^3=a,\ 1\le x\le y$ を満たす格子点 $(x,\ y)$ の個数が $4$ である」ような $a$ の最小値,および,それに対応する $4$ 個の格子点を求める問題です.

処理時間とメモリーとのバランスから,適当な上界を与え,出現をカウントしつつ,$x^3+y^3$ を生成する方針を採りました.利用した Map は,いわゆる array や dictionary にあたるオブジェクトですが,個数ではなく,格子点の $x$ 座標を係数とする多項式(リストを用いると 30sec 以上処理時間が延びます)を対応させることでワンパスを実現しています.

                                 GP/PARI CALCULATOR Version 2.13.1 (released)
                         amd64 running mingw (x86-64/GMP-6.1.2 kernel) 64-bit version
                          compiled: Jan 16 2021, gcc version 8.3-posix 20190406 (GCC)
                                           threading engine: single
                                (readline v8.0 enabled, extended help enabled)

                                    Copyright (C) 2000-2020 The PARI Group

PARI/GP is free software, covered by the GNU General Public License, and comes WITHOUT ANY WARRANTY
WHATSOEVER.

Type ? for help, \q to quit.
Type ?17 for how to get moral (and possibly technical) support.

parisize = 30000000000, primelimit = 500000
(18:56) gp > #
   timer = 1 (on)
(18:56) gp > n=7*10^12;d=4;
(18:56) gp > A=Map();e=d-1;for(i=1,floor((n/2)^(1/3)),i3=i^3;for(j=i,floor((n-i^3)^(1/3)),a=i3+j^3;if(mapisdefined(A,a,&p),mapput(A,a,p*x+i);P=mapget(A,a);if(poldegree(P)==e,print([a,Vec(P)])),mapput(A,a,i))));
[6963472309248, [2421, 5436, 10200, 13322]]
time = 11min,  57,031 ms.
(19:09) gp > #A
%3 = 161516889
(19:10) gp > 27792341+16331743+47676463+37255436+32460906
%4 = 161516889

また,領域 $\ n_1< x^3+y^3\le n_2\ $ により分割して

ta(d,n1,n2)={my(A=Map(),e=d-1);
for(i=1,(n1/2)^(1/3),for(j=floor((n1-i^3)^(1/3))+1,(n2-i^3)^(1/3),a=i^3+j^3;if(mapisdefined(A,a,&p),mapput(A,a,p*x+i);P=mapget(A,a);if(poldegree(P)==e,print([a,Vec(P)])),mapput(A,a,i))));
for(i=floor((n1/2)^(1/3))+1,(n2/2)^(1/3),for(j=i,(n2-i^3)^(1/3),a=i^3+j^3;if(mapisdefined(A,a,&p),mapput(A,a,p*x+i);P=mapget(A,a);if(poldegree(P)==e,print([a,Vec(P)])),mapput(A,a,i))));
#A;};

のように生成すれば

parisize = 10000000000

(18:34) gp > ta(4,1,5*10^11)
time = 2min, 9,063 ms.
%2 = 27792341

(18:35) gp > ta(4,5*10^11,1*10^12)
time = 1min, 20,875 ms.
%2 = 16331743

(18:35) gp > ta(4,1*10^12,3*10^12)
time = 3min, 57,625 ms.
%2 = 47676463

(18:35) gp > ta(4,3*10^12,5*10^12)
time = 3min, 12,968 ms.
%2 = 37255436

(18:35) gp > ta(4,5*10^12,7*10^12)
[6963472309248, [2421, 5436, 10200, 13322]]
time = 2min, 48,703 ms.
%2 = 32460906

といった並列化も可能です.

なお,今回の記事は jurupapa さんの記事 ( https://maxima.hatenablog.jp/entry/2021/06/10/003011 ) にインスパイアされたもので,以下は,最初に書いた maxima 版.

C:\maxima-5.45.0\bin>maxima.bat
Maxima 5.45.0 https://maxima.sourceforge.io
using Lisp SBCL 2.0.0
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.

(%i1)	(display2d:off,showtime:on)$
Evaluation took 0.0000 seconds (0.0000 elapsed) using 0 bytes.

(%i2)	(n:9*10^7,d:3)$ /* 87,539,319 */
Evaluation took 0.0000 seconds (0.0000 elapsed) using 0 bytes.

(%i3)	/* method 1: using a generating function */
	(q:subst([a=1,x=1],p:args(rat(expand(sum(a^i*x^(i^3)*sum(x^(j^3),j,i,(n-i^3)^(1/3)),i,1,(n/2)^(1/3)))))),
	    for i in q do if d=pop(q) then disp(pop(p)) else pop(p))$
(a^255+a^228+a^167)*x^87539319

Evaluation took 8.4690 seconds (8.4800 elapsed) using 919.738 MB.

(%i4)	/* method 2: using a counting array */
	(kill(C),C[i]:=[],for i:1 thru (n/2)^(1/3) do for j:i thru (n-i^3)^(1/3) do (k:i^3+j^3,if d=length(C[k]:append(C[k],[[i,j]])) then disp(k,C[k])))$
87539319

[[167,436],[228,423],[255,414]]

Evaluation took 2.6870 seconds (2.6810 elapsed) using 598.288 MB.

(%i5)	/* method 3: find adjacent occurrences in the sorted list */
		(A:sort(A0:create_list(i^3+j^3,i,1,floor((n/2)^(1/3)),j,i,floor((n-i^3)^(1/3)))),
		    a:pop(A),len:1,for i in A do if a=b:pop(A) then len:len+1 else (if d=len then (disp(a),return()),a:b,len:1),
		pos:sublist_indices(A0,lambda([x],x=a)),ij:create_list([i,j],i,1,floor((n/2)^(1/3)),j,i,floor((n-i^3)^(1/3))),disp(map(lambda([k],ij[k]),pos)))$
87539319

[[167,436],[228,423],[255,414]]

Evaluation took 0.9840 seconds (0.9860 elapsed) using 250.403 MB.

このうち,個数が $4$ の場合に適用できたのは method 3 のみでした.

C:\maxima-5.45.0\bin>maxima.bat -X "--dynamic-space-size 30000"
Lisp options: (--dynamic-space-size 30000)
Maxima 5.45.0 https://maxima.sourceforge.io
using Lisp SBCL 2.0.0
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.

(%i1) (display2d:off,showtime:on)$
Evaluation took 0.0000 seconds (0.0000 elapsed) using 0 bytes.

(%i2) (n:7*10^12,d:4)$
Evaluation took 0.0000 seconds (0.0000 elapsed) using 0 bytes.

(%i3) (A:sort(A0:create_list(i^3+j^3,i,1,floor((n/2)^(1/3)),j,i,floor((n-i^3)^(1/3)))),a:pop(A),len:1,for i in A do if a=b:pop(A) then len:len+1 else (if d=len then (disp(a),return()),a:b,len:1),pos:sublist_indices(A0,lambda([x],x=a)),ii:create_list(i,i,1,floor((n/2)^(1/3)),j,i,floor((n-i^3)^(1/3))),disp(map(lambda([k],[iik:ii[k],(a-iik^3)^(1/3)]),pos)))$

6963472309248

[[2421,19083],[5436,18948],[10200,18072],[13322,16630]]

Evaluation took 1295.1720 seconds (1297.1910 elapsed) using 320617.394 MB.

様々な言語によるコード https://rosettacode.org/wiki/Taxicab_numbers