## Muenchhausen Numbers

### December 6, 2019

A radix-10 number is a Muenchhausen number if it is equal to the sum of its digits, each digit raised to the power of itself. For instance, 3435 is a Muenchhausen number because 3**3 + 4**4 + 3**3 + 5**5 = 3435. Strict Muenchhausen numbers do not permit 0 as a digit, because 0**0 is undefined. Variant Muenchhausen numbers permit 0 as a digit, with 0**0 defined as 0. Muenchhausen numbers also exist with radices other than 10.

Your task is to find all the Muenchhausen numbers in radix-10. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

Pages: 1 2

### 10 Responses to “Muenchhausen Numbers”

1. Zack said

Here is my take on this in Julia (v. 1.1.1): https://pastebin.com/eppLUrZb

The result is this (fairly short) array: 1, 3435, 438579088

Note that I used the strict definition of Muenchhausen numbers, as it makes for a more interesting collection of numbers. Also, I set a limit of 1 billion when searching through the integer space to save time. Besides, it’s quite unlikely that there are any Muenchhausen numbers larger than that.

Have a nice weekend!

2. matthew said

Here’s a variation on that algorithm by Chai Wah Wu on http://oeis.org/A166623: generate all digit combinations of length up to radix+2 (this covers the range 2r^r), calculate the digit sum, then see if that sum has the same digits as the original combination – if so, we have a Münchhausen number. Use 0 or 1 for “zeroval” to set the assumed value for 0^0:

```def nextcombo(a,limit):
"""
Lexicographically next combination (with replacement)
If last for length, change to first for length+1
"""
for i in range(1,len(a)+1):
if a[-i] < limit-1:
a[-i] += 1
for j in range(1,i):
a[-j] = a[-i]
break
else:
a[:] = *(len(a)+1)

def digits(n,r):
""" (Reversed) digits of n in radix r """
if n == 0: return 
s = []
while n != 0:
s.append(n%r)
n //= r
return s

def münch(r, zeroval):
""" Return list of Münchhausen numbers for radix r """
digivals = [zeroval]+[d**d for d in range(1,r)]
res = []
s = 
while len(s) <= r+2:
dsum = sum(digivals[d] for d in s)
if sorted(digits(dsum, r)) == s:
res.append(dsum)
nextcombo(s,r)
return sorted(res)

```
3. matthew said

That python program runs happily up to radix 15, but I got bored waiting for hexadecimal, so converted it into C++, which finds 1, c76712ffc311e6e and c4ef722b782c26f, along with c4ef722b7820c2f and 123b2309ff572f6b if we put 0^0 = 0 (not the usual convention). Calculations are done with uint64_t which gives just enough bits to check up to 16^16, which is apparently an adequate bound. There is a possibility that the sum might wrap around and we get a false positive, but that doesn’t seem to happen here. Results agree with https://github.com/elizarov/MunchausenNumbers/ which has an interesting meet-in-the-middle algorithm.

Runtime is about 40 seconds with g++ -O3

```#include <stdint.h>
#include <stdio.h>
#include <algorithm>
bool nextcombo(int *a, int len,int limit) {
for (int i = len-1; i >= 0; i--) {
if (a[i] < limit-1) {
a[i] += 1;
for (int j = i+1; j < len; j++) {
a[j] = a[i];
}
return true;
}
}
return false;
}

int getdigits(uint64_t n, int r, int *a) {
int i = 0;
if (n == 0) a[i++] = 0;
else while (n != 0) {
a[i++] = n%r;
n /= r;
}
return i;
}

uint64_t ipow(uint64_t n, uint64_t m) {
if (m == 0) return 1;
if (m == 1) return n;
uint64_t t = ipow(n,m/2);
if (m%2 == 1) return n*t*t;
else return t*t;
}

int main() {
int N = 16, a[N];
uint64_t digival[N];
digival = 1;
for (int i = 1; i < N; i++) digival[i] = ipow(i,i);
for (int n = 1; n <= N; n++) {
for (int i = 0; i < n; i++) a[i] = 0;
do {
uint64_t sum = 0;
for (int i = 0; i < n; i++) sum += digival[a[i]];
int digits[N];
int dlen = getdigits(sum,N,digits);
if (dlen == n) {
std::sort(digits,digits+dlen);
for (int i = 0; i < dlen; i++) {
if (digits[i] != a[i]) {
goto end;
}
}
printf("%lx\n",sum);
}
end:;
} while (nextcombo(a,n,N));
}
}
```
4. matthew said

Python just got finished radix 16, so here are the full results with 0^0 = 1, and numbers in their proper radixes:

```2: 1, 10
3: 1, 12, 22
4: 1, 131, 313
5: 1
6: 1, 22352, 23452
7: 1, 13454
8: 1
9: 1, 31, 156262, 1656547
10: 1, 3435
11: 1, 18453278, 18453487
12: 1, 3a67a54832
13: 1, 33661, 2aa834668a, 4ca92a233518, 4ca92a233538
14: 1, 23, 26036, 45a0a04513cc, a992b5d96720d
15: 1, 4b1648420dcd0, 5a99e538339a43, 5acbc41c19e333, 5acbc41c19e400, 5d0b197c25e056
16: 1, c4ef722b782c26f, c76712ffc311e6e
```
5. learn some maths said

0^0 = 1, not undefined.

6. Zack said

Also, this video can help shed some light on the 0^0 topic: https://youtu.be/UKqoXMyp4_U

Whatever the case, it’s always useful to remember that the Mathematics we is largely based on conventions and assumptions we have made about certain matters (e.g. Euclid’s axiom about the parallel lines). So, regardless of what 0^0 is in reality (if there is a way to access that kind of truth), what we accept it to be may be a matter of debate. In my experience as an engineer / scientist, figuring out a solution to a problem using Math is much more useful than any philosophical discussions about Math properties. I hope that helps. Cheers!

7. Steve said

Klong version 20190209

```
t::0,1_{x^x}'!10
[0 1 4 27 256 3125 46656 823543 16777216 387420489]

f::{[n s]; n::\$x; s::0; x=+/{[n2]; t@n2::1:\$\$x}'n}

438579088{:[f'x; .p(x); ""]; x+1}:*0; :done

```

It figures out the first 3 (0, 1, 3435) as soon as I run the code. The last number (438579088) takes 2-3 hours.

8. Daniel said

Here’s a GMP-based C++ solution. Combinations are visited using a library I wrote, revdoor.

Example usage for some radixes up to 20 follows the code. The radix of the output numbers corresponds to the specified radix.

For the examples I ran, the runtime scaled approximately exponentially, base 4, for each increment in the radix. E.g., radix 16 took 1.5 minutes, radix 17 took 6 minutes, radix 18 took 23 minutes, radix 19 took 91 minutes, and radix 20 took 334 minutes.

```/*
* muenchhausen.c
*
* Build
*   \$ c++ -o muenchhausen muenchhausen.cpp -lgmpxx -lgmp
*/

#include <cassert>
#include <cstdlib>
#include <gmpxx.h>
#include <iostream>
#include <revdoor.hpp>
#include <string>
#include <vector>

using namespace std;

static bool is_muenchhausen(
const vector<int>& counter, mpz_class sum, int base, int t) {
int sizeinbase = mpz_sizeinbase(sum.get_mpz_t(), base);
int diff = sizeinbase - t;
// The result of mpz_sizeinbase "will be either exact or 1 too big".
if (diff < 0 || diff > 1) return false;
vector<int> counter2(counter);
while (sum != 0) {
mpz_class digit = sum % base;
sum /= base;
if (--(counter2[digit.get_si()]) < 0) return false;
}
for (int i = 0; i < (int)counter2.size(); ++i) {
if (counter2[i] != 0) return false;
}
return true;
}

int main(int argc, char* argv[]) {
assert(argc <= 2);
int base = 10;
if (argc == 2) base = stoi(argv);
assert(base > 2);    // for t <= base + 1 assumption
assert(base <= 62);  // requirement of mpz_sizeinbase
vector<mpz_class> lookup(base);
lookup = 1;
for (int i = 1; i < base; ++i) {
mpz_class x;
mpz_ui_pow_ui(x.get_mpz_t(), i, i);
lookup[i] = x;
}
for (int i = 1; i < base; ++i) {
if (i == lookup[i]) cout << mpz_class(i).get_str(base) << endl;
}
for (int t = 2; t <= base + 1; ++t) {
revdoor::combinations_with_replacement combs(base, t);
vector<int> init = combs.state();
vector<int> counter(base, 0);
mpz_class sum = 0;
for (int i = 0; i < (int)init.size(); ++i) {
++(counter[init[i]]);
sum += lookup[init[i]];
}
if (is_muenchhausen(counter, sum, base, t))
cout << sum.get_str(base) << endl;
int out1, in1, out2, in2;
while (combs.step(&out1, &in1, &out2, &in2)) {
if (out1 != in1) {
counter[out1] -= 1;
counter[in1] += 1;
sum -= lookup[out1];
sum += lookup[in1];
}
counter[out2] -= 1;
counter[in2] += 1;
sum -= lookup[out2];
sum += lookup[in2];
if (is_muenchhausen(counter, sum, base, t))
cout << sum.get_str(base) << endl;
}
}
return EXIT_SUCCESS;
}
```

Example usage:

```\$ ./muenchhausen 10
1
3435

\$ ./muenchhausen 16
1
c76712ffc311e6e
c4ef722b782c26f

\$ ./muenchhausen 17
1
33
2a422a28dc5f286
6a4f5cbg33ee2736
67a9b7eaa0485c0g

\$ ./muenchhausen 18
1
2900a7868c2
13abf45e3d8656d
2e78130c14a3g22f

\$ ./muenchhausen 19
1

\$ ./muenchhausen 20
1
6534
15g396h63be638d7c

```
9. Daniel said

Here’s a variant of my solution above, modified to conduct the search in parallel.

GMP and revdoor are the existing dependencies, with OpenMP added for multiprocessing functionality.

Each process explores a different subset of combinations. A combination unranking function is used to map combination indices to the corresponding combinations, which are used to determine where each process should start its search.

```/*
* muenchhausen.cpp
*
* Build
*   \$ c++ -fopenmp -o muenchhausen muenchhausen.cpp -lgmpxx -lgmp
*
* Usage
*         [--base INT]                   // defaults to 10     \
*         [--out-base INT]               // defaults to `base`
*/

#include <algorithm>
#include <cassert>
#include <cstdlib>
#include <cstring>
#include <iostream>
#include <string>
#include <vector>

#include <gmpxx.h>
#include <omp.h>
#include <revdoor.hpp>

using namespace std;

static bool is_muenchhausen(
const vector<int>& counter, mpz_class sum, int base, int t) {
int sizeinbase = mpz_sizeinbase(sum.get_mpz_t(), base);
int diff = sizeinbase - t;
// The result of mpz_sizeinbase "will be either exact or 1 too big".
if (diff < 0 || diff > 1) return false;
vector<int> counter2(counter);
while (sum != 0) {
mpz_class digit = sum % base;
sum /= base;
if (--(counter2[digit.get_si()]) < 0) return false;
}
for (int i = 0; i < (int)counter2.size(); ++i) {
if (counter2[i] != 0) return false;
}
return true;
}

static void unrank_revdoor_with_replacement(
const mpz_class& rank, int n, int t, vector<int>* output) {
mpz_class tmp;
mpz_class r(rank);
int x = n + t - 1;
for (int i = t; i > 0; --i) {
while (true) {
mpz_bin_uiui(tmp.get_mpz_t(), x, i);
if (tmp <= r) break;
--x;
}
output->push_back(x + 1 - i);
mpz_bin_uiui(tmp.get_mpz_t(), x + 1, i);
r = tmp - r - 1;
}
reverse(output->begin(), output->end());
}

static void search(int base, int out_base, int worker, int num_workers) {
vector<mpz_class> lookup(base);
lookup = 1;
for (int i = 1; i < base; ++i) {
mpz_class x;
mpz_ui_pow_ui(x.get_mpz_t(), i, i);
lookup[i] = x;
}
for (int i = 1 + worker; i < base; i += num_workers) {
if (i == lookup[i]) {
#pragma omp critical
cout << mpz_class(i).get_str(out_base) << endl;
}
}
for (int t = 2; t <= base + 1; ++t) {
mpz_class total_combs;  // across all workers
mpz_bin_uiui(total_combs.get_mpz_t(), base + t - 1, t);
mpz_class min_checks_per_worker = total_combs / num_workers;
mpz_class additional_checks = total_combs % num_workers;
mpz_class num_checks = min_checks_per_worker + (worker < additional_checks);
mpz_class init_rank = worker * min_checks_per_worker;
if (num_checks == 0) continue;
vector<int> init_state;
unrank_revdoor_with_replacement(init_rank, base, t, &init_state);
revdoor::combinations_with_replacement combs(base, t);
combs.set_state(init_state);
vector<int> counter(base, 0);
mpz_class sum = 0;
for (int i = 0; i < (int)init_state.size(); ++i) {
++(counter[init_state[i]]);
sum += lookup[init_state[i]];
}
if (is_muenchhausen(counter, sum, base, t)) {
#pragma omp critical
cout << sum.get_str(out_base) << endl;
}
mpz_class num_checked(1);
int out1, in1, out2, in2;
while (num_checked++ < num_checks) {
assert(combs.step(&out1, &in1, &out2, &in2));
if (out1 != in1) {
counter[out1] -= 1;
counter[in1] += 1;
sum -= lookup[out1];
sum += lookup[in1];
}
counter[out2] -= 1;
counter[in2] += 1;
sum -= lookup[out2];
sum += lookup[in2];
if (is_muenchhausen(counter, sum, base, t)) {
#pragma omp critical
cout << sum.get_str(out_base) << endl;
}
}
}
}

int main(int argc, char* argv[]) {
int base = 10;
int out_base = 0;
for (int i = 1; i < argc; ++i) {
if (strcmp("--base", argv[i]) == 0) {
assert(++i < argc);
base = stoi(argv[i]);
} else if (strcmp("--out-base", argv[i]) == 0) {
assert(++i < argc);
out_base = stoi(argv[i]);
}
}
assert(base > 2);    // for t <= base + 1 assumption
assert(base <= 62);  // requirement of mpz_sizeinbase
if (out_base <= 0)
out_base = base;
#pragma omp parallel
return EXIT_SUCCESS;
}
```

```\$ time OMP_NUM_THREADS=16 ./muenchhausen --base 18
1
2900a7868c2
13abf45e3d8656d
2e78130c14a3g22f

real    3m46.127s
user    47m56.826s
sys     0m1.805s

\$ time OMP_NUM_THREADS=1 ./muenchhausen --base 18
1
2900a7868c2
13abf45e3d8656d
2e78130c14a3g22f

real    31m58.072s
user    31m56.842s
sys     0m0.084s
```