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