Why does it converges to 1/4 and 3/4

Question

I've stumbled upon an "anomaly" in my code (Python).

Which one is bigger 2^3 or 3^2. Of course 3^2 beacause it's 9 and the other one is 8. What about 2^4 and 4^2? They're equal. Basically this was my idea: How many times is x^y different or equal to y^x, when x and y ranges from 1 to n (x and y are natural numbers). I wrote a code to it and for the first 60 there was nothing unusual. Then I was only interested in ratios ((times x^y = y^x) / (n^2 (all combinations))). I ran the code for bigger numbers and plotted it. I found a weird curve.

For the first 60 numbers the "diff ratio" was increasing, then it started to decrease with a "pattern" (I couldn't figure out what could it be exactly). I also noticed that the ratios not just decreasing and increasing, they are converging, and not to 0.5-0.5, but 0.75 and 0.25.

I can't figure out why does this happen, I tried to run the code with bigger numbers, but even after n=35000 it still hasn't crossed the previously mentioned limit. (I'm guessing it has something to do with the powers of 2)

What do you think, did I make a mistake somewhere in my code or I just didn't noticed a pattern? Does it have a mathematical background? Did I stumble upon a known mathematical problem?

First code with the scatter plot:

import numpy as np
import matplotlib.pyplot as plt

eqratios = [] 
diffratios = []

settings = np.arange(1, 1001, 10)



for setting in settings:
    diff = 0
    eq = 0
    for x in np.arange(1, setting, 1):
        for y in np.arange(1, setting, 1):
            if (x**y != y**x):
                diff += 1
            else:
                eq += 1

    eqratio = eq / ((setting-1)**2)
    diffratio = diff / ((setting -1) ** 2)
    eqratios.append(eqratio)
    diffratios.append(diffratio)

plt.figure(figsize=(10, 6))
plt.scatter(settings, diffratios)
plt.scatter(settings, eqratios)
plt.xlabel("n")
plt.ylabel("Blue: diff, Orange: eq")
plt.title("Scatter plot")
plt.grid(True)
plt.show()

And here is an improved, but complicated code to calculate with bigger numbers:

import numpy as np

settings = 35000
# 35k
# Diff Ratio: 0.7502916438079795
# Eq Ratio: 0.24970835619202048

# 21 k
# Diff Ratio: 0.7504859584121709
# Eq Ratio: 0.24951404158782908
# 19k
# Google Collab: - Diff Ratio: 0.7511174618072887 Eq Ratio: 0.24888253819271133
# PyCharm: Diff Ratio: 0.7505370849271196 Eq Ratio: 0.24946291507288043

n_values = settings - 1

total_combinations = n_values ** 2

x_range = np.arange(1, settings)
y_range = np.arange(1, settings)

x_matrix, y_matrix = np.meshgrid(x_range, y_range)

comparison_result = np.where(x_matrix**y_matrix != y_matrix**x_matrix, 1, 0)

diff_ratio = np.count_nonzero(comparison_result) / total_combinations
eq_ratio = 1 - diff_ratio  # Since eq_ratio + diff_ratio = 1
print("Diff Ratio:", diff_ratio)
print("Eq Ratio:", eq_ratio)

(If you have suggestions about how can I improve my code to run with bigger numbers, I'd be thankful)

(PS.: Sorry for my bad english, I hope you can understand it)

Answer 1

It's a programming error.

The only cases where $x^y = y^x$ for integers $x$ and $y$ are when $x = y$ or when $x$ and $y$ are $2$ and $4$ (in either order). So if you try all pairs $(x,y)$ where $x, y \in \{1, 2, \ldots, n\}$, you should have exactly $n + 2$ cases of "equal" out of $n^2$ total cases, so the "eq ratio" should be $\dfrac1n + \dfrac2{n^2}$. As $n$ (your setting variable) gets larger, the "eq ratio" should rapidly converge to zero.

But because you passed integer arguments to numpy.arange, the variables x and y had integer values. Not the type of arbitrary-size integer value that you get if you enter 583 to an interactive python prompt, but signed integers stored in a finite number of bits. (When I tried it, I found numpy.int32 as the type, but I suppose it is possible that your environment used numpy.int64.)

When you raise an even integer to an integer power, you get some multiple of a power of $2$. If the power is large enough, you will get a multiple of $2^{32}$; double the power, and you will get a multiple of $2^{64}$. When you do your exponentiation on signed integers of $32$ or $64$ bits, you actually get back a result equal to the correct value modulo $2^{32}$ or $2^{64}$, respectively, within the range of numbers representable by that integer type. Any multiple of $2^{32}$ (using $32$-bit integers) or $2^{64}$ (using $64$-bit integers) is congruent to $0$ modulo $2^{32}$ or $2^{64}$, respectively, so $0$ is the answer you get. Any power of an odd number will still be odd, however, so you will never get $0$.

As a result, when you try values of $x$ and $y$ over a large range, x**y is zero for almost all values of y for almost all even values of x, and y**x is zero for almost all values of x for almost all even values of y. In other words, most of the time when x and y are both positive, you are comparing zero to zero and your program says they're equal. (There are also a few unusual cases, for example, x**y and y**x are equal to the same non-zero value modulo $2^{32}$ when $x=14$ and $y=28$.) And x and y are both positive in $25\%$ of all cases.

If you pass floating-point numbers to numpy.arange, for example, np.arange(1.0, 1001.0, 10.0), then x and y will be floating-point numbers and you won't have the phenomenon where x**y comes out to zero. Instead, you will find that when the powers get large enough, the result is inf, and whenever both x**y and y**x evaluate to inf they are considered equal. In that case the "eq ratio" decreases toward zero at first, but once settings gets to about $150$ you start getting inf results for some of the numbers, and these eventually overwhelm the correct results.

Try plain old range instead of numpy.arange. Then you will get integers without limitation on the number of bits (or at least, the limit is so large you probably will not reach it) and you will get results that are in line with the mathematical results.

But essentially any "eq ratio" different from $\dfrac1n + \dfrac2{n^2}$ will be the result of a programming error or some limitation of the computer and programming language.

Answer 2

$x^y=y^x$ if and only if $x^{1/x} = y^{1/y}$. As Robert Israel's link mentions, $x^{1/x}$ increases from $0$ to $e^{1/e} \approx 1.44$ on $[0, e]$, and decreases from $e^{1/e}$ to $1$ on $[e, \infty)$. So the only way $x^{1/x} = y^{1/y}$ is possible is if $x < e < y$ (assuming $x<y$). When $x$ is restricted to integers, we only have $x=1$ and $x=2$.

When $x=1$, there is no $y > e$ such that $y^{1/y} = 1^{1/1} = 1$ since $y^{1/y}$ approaches $1$ from above as $y \to \infty$. So $x$ must be $2$.

The only $y > e$ satisfying $y^{1/y} = 2^{1/2}$ is $y=4$.

So the only $(x,y)$ pairs satisfying $x^y=y^x$ in your simulation should be

• $x=y$,
• $(x,y) = (2, 4)$, and
• $(x,y) = (4,2)$.

---

Your code is flawed because the large powers are not being computed correctly. When you use np.arange, the x and y have type np.int64, which for some reason returns zero when trying to compute large powers. This leads to a spurious match when both $x^y$ and $y^x$ are large.

np.int64(8) ** np.int64(2)
> 64

np.int64(8) ** np.int64(32)
> 0