Re: [FUN] CodeGolf: solve equation x^x=C

excellent ! (note that your did not include the function header "function f(x)" and trailer "end" and the code "=x" to print the result on the Lua console).

It works because the range from 0 to 395 is small enough and it assumes that math.random() will return an exact value.

This can create an infinite loop if this condition x*x=C is never reached exactly, because:

- random() will not explore all the possible bits, as random() does not return all the possible values for the 52 bits of numbers in (0,1).

- even if you have the best approximation possible of x, x^x=C may still be wrong (and there's no upper bound of the relative error on the value of x).

So let's assume that random() is able to generate at least 31 bits of precision with a resonably flat distribution, the exact value of log2(x^x) may be different from log2(C) by a relative factor of up to 15.5*x

i.e. you must stop the loop when log2(x^x)/log2(C) <= 15.5*x,

i.e. equivalently when ln(x^x)/ln(C)/15.5 <= x

or when ln(x^x)/x <= 15.5*ln(C)

You can see here why 15.5 is used in the problem, it's not "magic" but it's log2(31)

----

And why using random(), which is unnecessarily slow ?

We could just enumerate the 31 bits of precision in the open range (from 3.0 to 395.0). And this exploration can be much faster if you have a starting approximate (but lower) bound of the solution so that it will enumerate from a better start position (in which case it will require much less loops to reach the best lower bound that matches the expected precision of the result.

What is needed is only to estimate the lower bound of the base-2 integral exponent of x.

Interestingly, we have : x= 1024/log_2(x) so log_2(x)=1024/x

This allows starting the "for" loop to explore the 31 bits of precision of values of x (from 3.0 to 395.0) at a known integer position which is (1024 div x).

You just need an integer "for" loop (and your loop will not explore bits for very long because you are quite near from the expected precision at the start !).

Lets say that (exponent=1024 div x) is the starting value for your integer for-loop, the for-loop will run for at most log2(exponent), and given that exponent is below 1024, it will never run for more than 10 times (and 5 times on average)... and it is warrantied to give a correct result of x with at last 31 bits of precision.

now the 10 loops can be easily unrolled into inline code performing successive tests. But an even better solution is to perform a binary search in the small constant table of 10 constant values.

Le jeu. 27 juin 2019 à 17:11, Coda Highland <chighland@gmail.com> a écrit :

On Thu, Jun 27, 2019 at 9:57 AM Philippe Verdy <verdy_p@wanadoo.fr> wrote:
but a basic asymptotic approcimation cannot work for any number x that sastisfies x^x > 15.5 (which can be very large).
So let's see where x^x does not return positive infinite (starting at 2^1024).
We need the mathematic solution for x^x = 2^1024.

But let's start by how we can compute the base-2 exponent (which will be betwen 3 and 1024 for this problem) to get a rough estimate of the order of magnitude
log_2(x^x) =x*log_2(x) = 1024
so x= 1024/log_2(x). and then x must be below

The problem states that x^x>15.5, so x >15.5/log_2(x) (this will be used for tuning the value of the mantissa and correct the effect of the approximate base-2 exponent.

There's not a lot of values of x that can qualify, it is in a small range in the open range (3.3, 395)
If x is to be restricted to be an integer, this x is in the integer range 4..394 (because 395^395>2^1024 and gives a positive infinite value for 64-bit IEEE doubles)

But in that case a table lookup is not possible (the integer range 4..394 would contain 390 double values in a static table)

So what is to approximate if the base-2 integer exponent of x^x and then we can correct the mantissa; and this requires amultiplication from the correct mantissa, by the approximation of the base-2 exponent.

This will still produce large errors compared to the recursive method using the derivative function of (x => x^x), which is (x => (ln(x) + 1) * x^x).

The recursive function will then loop by dividing the current approximation of x^x, by (log(x)+1), to form the new approximation by taking the average (it will have a fast convergence, but now you need a way to compute ln(x) or you can use math.ln(x).

Thanks to your discovery of the upper bound, I can write a solution that (technically) works in 52 bytes.

x=0;while x*x~=C do x=math.random()*395 end return x

;)

/s/ Adam