Autumn Challenge - Three Cubes Sum

[j4n_bur53]

This is already a post from comp.lang.prolog. But it seems it is reasonably challenging and could be fun. So I am also crossposting here.

It took Deep Thought 7½ million years to compute and check the answer, which turns out to be 42. But can we also compute via Prolog positive integer numbers x, y and z such that:

 x^3 + y^3 + 42 = z^3

When I use a nearly quadratic algorithm, I can compute k=87. Here is an example code that does the computation and the time measurement is from a slow Mac laptop:

?- time((above(0,S), H is S//2, between(0,H,X), 
         Y is S-X, C is X^3+Y^3+87, rootrem(C,3,Z,0))).
% Up 70,783 ms, GC 612 ms, Threads 69,768 ms (Current 09/17/19 21:31:34)
X = 1972,
Y = 4126,
Z = 4271 

Other practical ways to speed up a Prolog solution so that we could reach k=42, or is this out of reach concerning computing time?

Edit 18.09.2019:
SWI-Prolog timing:

?- time((between(0,inf,S), H is S//2, between(0,H,X),
    Y is S-X, C is X^3+Y^3+87, nth_integer_root_and_remainder(3,C,Z,0))).
% 37,217,891 inferences, 25.142 CPU in 25.408 seconds (99% CPU, 1480280 Lips)
X = 1972,
Y = 4126,
Z = 4271

[EricGT]

You forgot to add the link: 42 is the new 33 - Numberphile

It shows a simplified version of the concept of how to dramatically simplify the problem, but because the suspected numbers to the solution are 10*16 or greater it sill is a brute force problem.

See: Cracking the problem with 33

[PaulBrownMagic]

According to the Numberphile interview, they increased the computing power by an order of magnitude to go from 33 to 42. Rumor is they changed little about the algorithm, to the point where they’re not going to publish a paper about it. In the linked video Andrew says the changes made were to make it run at a larger scale, they’re focusing on the use of community powered parallel computing more.

[j4n_bur53]

@j4n_bur53

positive integer numbers x, y and z such that:

 x^3 + y^3 + 42 = z^3

Can this be solved via ZZD and the axiom of extensionality?
Or maybe the new CLP(Z) from Scryer Prolog can do it?

Or are we in the mist of Hilbert’s 10-th problem. Maybe Poincaré
would have objected we simply need new problem specific rules?

[kuniaki.mukai]

@j4n_bur53

Can this be solved via ZZD and the axiom of extensionality?

I don’t find any positive clue at first glance. IMHO, ZDD might be useless for the 3 cube problem. Your solutions looks smarter enough, which I checked works well.

[j4n_bur53]

Furtther question, does SWI-Prolog have modinv/3, also a function popular
in cryptography? It can be used for a Chinese Remainder algorithm.

But currently banging my head on a slow (^)/2 for smallints. In modular
algorithms the arguments for (^)/2 are small. So I now get:

/* SWI-Prolog 8.5.17 */
?- time((modular([11,7,5], X, Y, Z), fail; true)).
% 15,388,917 inferences, 3.969 CPU in 3.963 seconds (100% CPU, 3877522 Lips)
true.

/* Jekejeke Prolog 1.5.4 */
?- time((modular([11,7,5], X, Y, Z), fail; true)).
% Threads 2,078 ms, GC 10 ms, Up 2,122 ms (Current 09/20/22 11:03:10)
true.

But I only solved this k=9 equation, not yet the k=42 problem:

/* x^3 + y^3 + 9 = z^3 */
?- modular([11,7,5], X, Y, Z).
X = 216,
Y = 52,
Z = 217 ;
X = 52,
Y = 216,
Z = 217 ;
false.

[brebs]

Hopefully of interest: clpBNR solves that quickly:

?- time(([X,Y,Z]::integer(0, _), {X**3 + Y**3 + 9 == Z**3, Y >= X}, solve([X,Y,Z]))).
% 1,438,518 inferences, 0.147 CPU in 0.147 seconds (100% CPU, 9783542 Lips)
X = 52,
Y = 216,
Z = 217 ;

[ridgeworks]

@brebs

Hopefully of interest: clpBNR solves that quickly:

To be fair, this isn’t an apples-to-apples comparison; it only finds one solution. And if you don’t bound the variables somehow, it could take a very long time to prove there are no other solutions. A more accurate comparison:

?- time(([X,Y,Z]::integer(0, _), {X**3 + Y**3 + 9 == Z**3, Z=<1000}, solve([X,Y,Z]))).
% 1,956,751 inferences, 0.383 CPU in 0.383 seconds (100% CPU, 5115074 Lips)
X = 52,
Y = 216,
Z = 217 ;
% 1,338,831 inferences, 0.262 CPU in 0.263 seconds (100% CPU, 5109086 Lips)
X = 216,
Y = 52,
Z = 217 ;
% 57,603,408 inferences, 11.469 CPU in 11.475 seconds (100% CPU, 5022376 Lips)
false.

BTW, clpBNR uses the builtin SWIP implementation of integer exponentiation (ar_pow) wrapped in an interval arithmetic primitive written in Prolog so hardly optimized for this particular problem. I suspect k=42 would take a very long time.

[j4n_bur53]

Using modular operations again, you can make it a tick faster by using
only two modules instead of three modules. But it still depends
heavily on the performance of ar_pow():

/* SWI-Prolog 8.5.17 */
?- time((modular([15,16], X, Y, Z), fail; true)).
% 3,816,486 inferences, 1.766 CPU in 1.762 seconds (100% CPU, 2161550 Lips)
true.

/* Jekejeke Prolog 1.5.4 */
?- modular([15,16], X, Y, Z).
X = 216, Y = 52, Z = 217;
X = 52, Y = 216, Z = 217;
fail.

?- time((modular([15,16], X, Y, Z), fail; true)).
% Threads 594 ms, GC 5 ms, Up 591 ms (Current 09/20/22 20:21:30)
true.

Thats the source code, it is ultimately intended to beat CLP(Z)
from Scryer Prolog some time in the future into dust, it also uses
ar_pow(), maybe better would be ar_modpow(), but for the small

values ar_modpow() wouldn’t give some bang:

(Credits for the acronym cra for Chinese Remainder Algorithm in the
context of Prolog go to Amzi! Prolog, I found the acronym here:
https://amzi.com/manuals/amzi/pro/ref_math.htm )

modular([P|L], X, Y, Z) :-
   find(P, A, B, C),
   many(L, A-P, B-P, C-P, X-_, Y-_, Z-_),
   X^3+Y^3+9-Z^3 =:= 0.

many([], A, B, C, A, B, C).
many([P|L], U, V, W, A, B, C) :-
   find(P, X, Y, Z),
   cra_combine(X-P, U, S),
   cra_combine(Y-P, V, T),
   cra_combine(Z-P, W, R),
   many(L, S, T, R, A, B, C).

find(M, X, Y, Z) :-
   N is M-1,
   between(0, N, X), A is X^3,
   between(0, N, Y), B is Y^3,
   between(0, N, Z),
   (A+B+9-Z^3) mod M =:= 0.

/**
 * Chinese Remainder Algorithm
 * https://en.wikipedia.org/wiki/Chinese_remainder_theorem#Existence_%28constructive_proof%29
 */
% cra_combine(+Pair, +Pair, -Pair)
cra_combine(A-N, B-M, C-P) :-
   cra_egcd(N, M, X, Y),
   P is N*M,
   C is (A*M*Y+B*N*X) mod P.

% cra_egcd(+Integer, +Integer, -Integer, -Integer)
cra_egcd(_, 0, 1, 0) :- !.
cra_egcd(A, B, X, Y) :-
    divmod(A, B, Q, R),
    cra_egcd(B, R, S, X),
    Y is S - Q*X.

[j4n_bur53]

@ridgeworks

clpBNR uses the builtin SWIP implementation of integer exponentiation

It easily beats CLP(FD):

?- time(([X,Y,Z]::integer(0, 239), {X**3 + Y**3 + 9 == Z**3}, solve([X,Y,Z]), fail; true)).
% 3,494,431 inferences, 0.969 CPU in 0.968 seconds (100% CPU, 3607155 Lips)
true.

?- time(([X,Y,Z] ins 0..239, X^3+Y^3+9 #= Z^3, label([X,Y,Z]), fail; true)).
% 36,572,338 inferences, 5.109 CPU in 5.117 seconds (100% CPU, 7157889 Lips)
true. 

Cool!

[CapelliC]

@j4n_bur53: running your code on SWI-Prolog 8.5.17 (Windows 10), I get this problem (posted here mainly in hope that could be useful for @ridgeworks)

ERROR: Arithmetic: evaluation error: `float_overflow'
ERROR: In:
ERROR:   [22] clpBNR:lower_bound_val_(real,-1.0Inf,_865120)
ERROR:   [21] clpBNR:int_decl_(real,(-1.0Inf,1.0Inf),_865154) at c:/users/carlo/appdata/local/swi-prolog/pack/clpbnr/prolog/clpbnr.pl:430
ERROR:   [18] clpBNR:build_args_([_865202,...|...],[_865214|_865216],[real,real|...],_865232/_865234,_865198) at c:/users/carlo/appdata/local/swi
...

[j4n_bur53]

@CapelliC

running your code

What do you mean by “your code”? I was using a slight modification (note the
integer constraint (0,239) ) of the query of ridgeworks. What was the query
that provoked your error? What version of clpBNR you were using?

Here you see, everything works fine on my side:

[CapelliC]

I wanted the number of solutions, then I factored out the query in a module…

:- module(three_cubes,
          [p/3
          ,q/3
          ]).

:- use_module(library(clpBNR)).
:- use_module(library(clpfd)).

p(X,Y,Z) :-
    [X,Y,Z]::integer(0, 239), {X**3 + Y**3 + 9 == Z**3}, solve([X,Y,Z]).

q(X,Y,Z) :-
    [X,Y,Z] ins 0..239, X^3+Y^3+9 #= Z^3, label([X,Y,Z]).

but you’re right, running the code in REPL is working…

?- use_module(library(clpBNR)).
% *** clpBNR v0.10.1 ***.
true.

?- time(([X,Y,Z]::integer(0, 239), {X**3 + Y**3 + 9 == Z**3}, solve([X,Y,Z]), fail; true)).
% 3,494,604 inferences, 1.250 CPU in 1.260 seconds (99% CPU, 2795683 Lips)
true.

(side note: I got the error from three_cubes module also before I added the clpfd related part)

[CapelliC]

follow up…

?- use_module(three_cubes).
true.

?- p(X,Y,Z).
X = 52,
Y = 216,
Z = 217 .

so, in the end, seems to be my fault, I’m not new to stumble in UX corner cases…

[j4n_bur53]

I used integer(0, 239) as a constraint, since the modular
approach, used the two modulos 15 and 16, and we have:

?- X is 15*16.
X = 240.

But I am still bugged with a slow ar_pow(), also switching
to (**)/2 doesn’t help, same slowness in SWI-Prolog:

?- time((between(1,1000000,_), _ is 33^10, fail; true)).
% 2,000,000 inferences, 0.500 CPU in 0.504 seconds (99% CPU, 4000000 Lips)
true.

?- time((between(1,1000000,_), _ is 33**10, fail; true)).
% 2,000,000 inferences, 0.500 CPU in 0.502 seconds (100% CPU, 4000000 Lips)
true.

Switching to (**)/2 would be feasible logic wise, since
it returns integer, unlike the ISO cores standard which would
rather prefer a float:

?- X is 33**10.
X = 1531578985264449.

But I guess the point of ISO core standard introduction of
a separate (^)/2 is that it would not use Math.pow() for floats
under the hood, but some fast path.

What does (**)/2 even do in SWI-Prolog?

[jan]

@j4n_bur53

What does (**)/2 even do in SWI-Prolog?

See the manual