Exhausting decimal Munichausen numbers in Prolog under convention 0^0 = 1

[kuniaki.mukai]

A decimal number 3435 has a special property that

3435 = 3^3+4^4+3^3+5^5,

which is called a Munichausen Number: Munichausen Number (Wikipedia) As 1=1^1 1 is a Munichausen.

Is there any other Munichausen number ? It is not difficult to see that a decimal Munichausen must be less than 10^{10}, provided that convnention 0^0=1 is exploited. Note that
SWI-Prolog exploits it.

To see that 1 and 3435 are all Munichausen, I tried to check it by naive generate-test way for all positive integers < 10^{10.} But 10^{10} seemed to be too big for
that naive method, though I did not try any constraint based technology yet.

Anyway finally building a (irreducible) BDD tree for all 10^{10} decimal digits lists and check
Munichausen for each list of digits. Fortunately it terminates showing the expected result that a decimal Munichausen are only 1 and 3435, after running the query below on m2 mac mini for more than one hour,

However, I am afraid that using BDD is too much brute-force for the problem. I would like to hear about any clever prolog method to exhaust Munichausen numbers less than specified limit (e.g. 10^{10}) . Thanks in advance.

% ?- time( (zdd { K=10, numlist(0, 9, Ns) }, bdd_list_raise(Ns, K, X),
% munichausen(0, 0, X), {fail})).

%@ munichausen_number(3435)
%@ munichausen_number(1)
%@ % 100,000,132,272 inferences, 4716.645 CPU in 4737.086 seconds (100% CPU, 21201541 Lips)
%@ false.

munichausen(X, X, 1, _):- X>0, writeln(munichausen_number(X)).
munichausen(X, Y, U, S):- U>1,
	cofact(U, t(I, L, R), S),
	munichausen(X, Y, I, L, R, S).
%
munichausen(X, Y, _, L, _, S):-	munichausen(X, Y, L, S).
munichausen(0, _, 0, _, R, S):-!, munichausen(0, 0, R, S).
munichausen(X, Y, I, _, R, S):-
	X0 is 10*X + I,
	Y0 is Y + I^I,
	munichausen(X0, Y0, R, S).

[j4n_bur53]

Looks like using 0^0 = 0 is easier, but possibly wrong?

?- digits(4, 0, 0, R).
R = 0 ;
R = 1 ;
R = 3435 ;
false.

The code reads:

digits(0, P, Q, R) :- !, P = Q, R = P.
digits(N, P, Q, R) :-
   M is N-1,
   between(0, 9, D),
   S is P*10+D,
   (D = 0 -> T = Q; T is Q+D^D),
   digits(M, S, T, R).

Edit 13.03.2023
BTW: You can speed up the search by interval estimate and check,
and bring it below 1 hour. But I found a new Münchhausen number.
Possibly some post processing filtering to 0^0 = 1 can help exclude it,

but we are not sure whether some 0^0 = 1 solution is missing:

?- time((digits2(10, 0, 0, R), write(R), nl, fail; true)).
0
1
3435
438579088
% 2,211,759,677 inferences, 311.531 CPU in 312.360 seconds
(100% CPU, 7099640 Lips)
true.

The new number is indeed Münchhausen under the 0^0 = 0 convention:

?- X is 4^4+3^3+8^8+5^5+7^7+9^9+8^8+8^8.
X = 438579088.

Thats the interval estimate and check version, basically what
CLP(FD) could do through its interval propagation, but still CLP(FD)
might not do it the way this hand crafted code does it:

digits2(0, P, Q, R) :- !, P = Q, R = P.
digits2(N, P, Q, R) :-
   M is N-1,
   H is 10^M,
   between(0, 9, D),
   S is P*10+D,
   (D = 0 -> T = Q; T is Q+D^D),
   S*H =< T+M*387420489,
   T =< S*H+H-1,
   digits2(M, S, T, R).

[kuniaki.mukai]

Initial version of codes of mine did not work for the case N = 9, where N is the length of decimal number, i.e. N=3 for 123.

It is strange for me why BDD version works even for N=10, which is enough for the exhaustive search for the Munichausen numbers. As the BDD version is also based on generate and test method, I doubt that my solution is good.

Does you codes terminate well for N=10 ?

[j4n_bur53]

@kuniaki.mukai

Does you codes terminate well for K=10 ?

Yes, it takes 312.360 seconds on my machine. So its around 5 minutes,
see my previous post. I edited it. But you might want to fix the 0^0 = 0
issue, unless you are happy with this convention as well.

According to this paper, if you search the interval [1,2*R^R],
where R is the radix, you know whether there is a Münchhausen
Number or not, in the infinite interval [1,oo), since there are

only finitely many Münchhausen Numbers and 2*R^R is an
upper bound. I don’t know, is this paper correct? Whats the argument,
I only picked it up by chance, dont understand it yet:

On a curious property of 3435
Daan van Berkel - 2009
https://arxiv.org/abs/0911.3038

So Kuniaki Mukais choice of interval [1,10^10] for R=10 is
already quite close in settling Münchhausen numbers question
once and for all, right?

[kuniaki.mukai]

@j4n_bur53

Yes, it takes 312.360 seconds on my machine. So its around 5 minutes,

Sound great. I will check your codes.

My codes is flexible between the two conventions, and for choice of base (radix). I will play around on these variants.

@j4n_bur53

So Kuniaki Mukais choice of interval [1,10^10] for R=10 is
already quite close in settling Münchhausen numbers question
once and for all, right?

I have to check this. Thanks.

[kuniaki.mukai]

No existence of decimal Munichausen number > 10^{10} came to me from the following simple experimental arithmetic query． Am I missing ?

% ?- between(1, 100, N), 10^N > N*(9^9), writeln(N), fail.

[kuniaki.mukai]

Thanks for links. I will go to the links soon.

I ran a script go/3 to enumerate all Munichausen numbers upto radix 10 for both conventions 0^0=0 and 0^0=1. It seems to return to correct answers. It took time about 4 hours, which should be revised.

go(R, C, Nums) enumerates Nums of all Munichausen numbers for radix R and convention flag C.

% ?- between(2, 10, I), member(J, [0,1]),
%	go(I, J, Nums), writeln(go(I,J) = Nums),
%	fail; true.
%@ go(2,0)=[1,0]
%@ go(2,1)=[2,1]
%@ go(3,0)=[8,5,1,0]
%@ go(3,1)=[8,5,1]
%@ go(4,0)=[55,29,28,1,0]
%@ go(4,1)=[55,29,1]
%@ go(5,0)=[264,28,1,0]
%@ go(5,1)=[1]
%@ go(6,0)=[3416,3164,1,0]
%@ go(6,1)=[3416,3164,1]
%@ go(7,0)=[3665,1,0]
%@ go(7,1)=[3665,1]
%@ go(8,0)=[257,256,1,0]
%@ go(8,1)=[1]
%@ go(9,0)=[16871323,923362,917139,96446,28,27,1,0]
%@ go(9,1)=[923362,96446,28,1]
%@ go(10,0)=[438579088,3435,1,0]
%@ go(10,1)=[3435,1]

[j4n_bur53]

I get a faster method with Prolog first argument indexing:

?- time((canonball, munchhausen(_), fail; true)).
% 4,633,344 inferences, 0.594 CPU in 0.589 seconds (101% CPU, 7803527 Lips)
true.

?- munchhausen(R).
R = 0 ;
R = 1 ;
R = 3435 ;
R = 438579088 ;
false.

The code is based on the observation that the fixpoint
equation for a split number m*10^k + n can be reordered
so that m is on one side and n is on the other side:

m*10^k - F(m) = F(n) - n

The rest is tabulating F(n) - n:

canonball :-
   retractall(cache(_,_)),
   between(0, 99999, N), map(N, Y), C is Y-N,
   assertz(cache(C, N)), fail; true.

munchhausen(R) :-
   between(0, 99999, M), map(M, X), B is 100000*M-X,
   cache(B, N), R is 100000*M+N.

map(0, X) :- !, X = 0.
map(N, X) :-
   M is N//10,
   map(M, Y),
   D is N mod 10,
   (D = 0 -> X=Y; X is Y+D^D).

[kuniaki.mukai]

@j4n_bur53

I get a faster method with Prolog first argument indexing:

Congratulation !

It seems worth more than to say great.

[j4n_bur53]

@kuniaki.mukai

% ?- between(2, 10, I), member(J, [0,1]),
%	go(I, J, Nums), writeln(go(I,J) = Nums),
%	fail; true.

But computing the table up to R=16 is still an open problem for me.
Need more speed-up ideas. Up to R=10 I can verify your table,
It now needs ca. 5 seconds with Dogelog Player on nodeJS:

?- time((between(2,10,R), write(R), write(': '),
(canonball(R), munchhausen(R,S), write(S), write(' '), fail; true),
nl, fail; true)).
2: 0 1
3: 0 1 5 8
4: 0 1 28 29 55
5: 0 1 28 264
6: 0 1 3164 3416
7: 0 1 3665
8: 0 1 256 257
9: 0 1 27 28 96446 917139 923362 16871323
10: 0 1 3435 438579088
% Time 5170 ms, gc 15 ms, 2877753 lips
true.

But if you have a larger radix R, the square root of the range that
needs to be first tabled and later scanned increases more than only
exponential. Its not that we have a constant C and the size is C^R,

its much worse, namely R^(R//2). Quite a challenge to go higher.

[kuniaki.mukai]

@j4n_bur53

But computing the table up to R=16 is still an open problem for me.

I have revised the estimation of upper-bound of width of possible Munichausen numbers to the radix - 1. I am fairy sure of correctness of this revision. In fact, I have proved it, which uses interval picturs [r^n, r^(n+1)) (left closed, right open) as you used in a post of of yours in this thread. I’m not sure this revision is relevant to your open problem for case R=16, but it would be good if it were the dase.

[j4n_bur53]

When R, the radix, is a parameter, I am using this code:

canonball(R) :-
   retractall(cache(_,_)),
   U is R^(R//2)-1,
   between(0, U, N), map(N, R, Y), C is Y-N,
   assertz(cache(C, N)), fail; true.

munchhausen(R, S) :-
   U is R^(R//2),
   V is R^(R-R//2)-1,
   between(0, V, M), map(M, R, X), B is U*M-X,
   cache(B, N), S is U*M+N.

map(0, _, X) :- !, X = 0.
map(N, R, X) :-
   M is N//R,
   map(M, R, Y),
   D is N mod R,
   (D = 0 -> X=Y; X is Y+D^D).

I didn’t try R=16 yet. I guess it takes quite long or it will
crash the computer. The cannonball iteration will use a
U=4294967295, that gives a quite big cache/2 table.

I think its currently eating more and more memory, and
will soon crash. I see:

?- canonball(16).

Yeah after some while I get: