lua-users home
lua-l archive

[Date Prev][Date Next][Thread Prev][Thread Next] [Date Index] [Thread Index]


Disregard random() which is a bad idea (and which is an unnecessary external dependancy). Using a binary search with trailing recursive calls, and at most log2(52) trailing recursions, means that you'll need at most 6 recursions to get the best approximation for 64-bit doubles.

At each recursion, you normally need to evaluate the stop condition on the maximum relative error, which is  log_2(x)=1024/x, but you can just continue to perform the 6 recursions based on the lowerbound/upperbound test of the value of log2(x)*x, depending if it is lower than 1024 or not.

And this test can be optimized because actually we just need to know  if floor(log2(x)*x)<1024:
we just need the 10 first bits of log2(x)/x, i.e. testing if log2(log2(x)*x)<10, or just log2(log2(x))<10*x, or of log2(log2(x))*x/10 < 1

Further tricks will immediately follow to optimize this test: we know that x is between 3.0 and 395.0, so log2(log2(x)/x/10 is also bound
So we just need its integer upper bound to compare it with the integer 1:

In summary we need only 1 bit of precision for log2(log2(x)/x/10 for any x in (3.0..395.0) to perform the boolean test needed for the  recursion of a standard recursive binary search (which will run 6 times always). I bet that you can optimize the two branches of the "if" in a single one, using bit fiddling to determine which value to pass to the trailing recursion (so that you'll eliminate the if/else and just recurse on f(new value of x).

And don't forget that we know in advance the number (6) of recursions that can then be unrolled, each recursion computing 10 bits of mantissa in the approximation of x; you get the result then by just outputing the result of "f f f f f f(some value)".

The trick is to find the bit fiddling to compute if we need, to pass (value1 or value2) where value1 and value2 are the ower and upper bound.
But may be we don't even need this fiddling when value1 and value2 are just values of (x^x) and (x*x^x),
i.e. it's enough to pass to the recursion call the value x^x*(1+(log2(x)*x<1024 and 1 or 0))

I must not be very far from the solution !

Le jeu. 27 juin 2019 à 17:55, Roberto Ierusalimschy <roberto@inf.puc-rio.br> a écrit :
> - 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).

Not true if you are using 5.4 :-)