Next floating point number

[Bengt Richter]

On Sun, 18 Dec 2005 10:27:16 +1100, Steven D'Aprano <steve at REMOVETHIScyber.com.au> wrote:

>On Sat, 17 Dec 2005 09:26:39 +0000, Bengt Richter wrote:
>
>
>> I wonder if this won't work (for IEEE 754 double that is)
                                ^^^^^^^^^^^^^^^^^^^[1]
>> 
>> from math import frexp
>> def nextf(x, y):
>>     f,e = frexp(x)
>>     if (f==0.5 or f==-0.5) and x>=y: eps = 2.0**-54
>>     else: eps = 2.0**-53
>>     if x<y: return (f+eps)*2.0**e
>>     else: return (f-eps)*2.0**e
>
>Doesn't that assume floats are 64 bits? Is it possible that they might
>not be on some platform or another?
yes [1], and yes ;-)

I wonder if frexp is always available, and if the above would generalize
to something that could be self-tuning via a decorator or conditional defs.
ISTM I recall a FP format based on shifting 4 bits at a time to normalize,
and keeping an exponent as a multiple of 16, which might need a bunch
of epsilons. Obviously too, you wouldn't want to calculate things like 2.**-53
more than once. Also for a limited range of e integers, 2.**e is probably
faster looked up from a precalculated list p2[e]. The math module could
also expose an efficient multiply by a power of two using FSCALE if
the pentium FPU is there. And there's probably other stuff to take advantage
of, like doing it mostly without the FPU, a la Tim's struct thing. Specialized
C code to do the function could be pretty fast.
I'm a little curious what you are using this for.

Regards,
Bengt Richter