DoubleFloats.jl icon indicating copy to clipboard operation
DoubleFloats.jl copied to clipboard

`exponent` and `significand` functions

Open mcabbott opened this issue 3 years ago • 4 comments

The exponent function gives an error on Double64, and does not seem to be covered by tests: https://app.codecov.io/gh/JuliaMath/DoubleFloats.jl/blob/master/src/math/prearith/prearith.jl

julia> exponent(123.4)
6

julia> exponent(Double64(123.4))
ERROR: DomainError with 0.0:
Cannot be ±0.0.
Stacktrace:
 [1] (::Base.Math.var"#throw2#6")(x::Float64)
   @ Base.Math ./math.jl:846
 [2] exponent
   @ ./math.jl:851 [inlined]
 [3] exponent(x::Double64)
   @ DoubleFloats ~/.julia/packages/DoubleFloats/h3HrU/src/math/prearith/prearith.jl:64

julia> exponent(big(123.4))
6

I think it's calling exponent(0.0), to make a tuple, like that made by significand:

julia> significand(123.4)
1.928125

julia> significand(Double64(123.4))
(1.928125, 0.0)

julia> significand(big(123.4))
1.928125000000000088817841970012523233890533447265625

More generally, is this the right thing to return? I would expect the results to obey x == significand(x) * 2^exponent(x), but this isn't possible if they return tuples. Maybe significand should return another Double64? If these tuples are needed, they could be internal functions, not overloads of Base.

The context is that it would be nice if https://github.com/JuliaStats/LogExpFunctions.jl/pull/48 could work on AbstractFloat.

mcabbott avatar May 09 '22 00:05 mcabbott

Turns out there's a better function for this. Which is defined & runs without error. But returns a tuple, where it's documented to return a number, I think:

julia> frexp(123.4)
(0.9640625, 7)

julia> frexp(Double64(123.4))
((0.9640625, 7), (0.0, 0))

julia> @which frexp(Double64(123.4))
frexp(x::DoubleFloat{T}) where T<:Union{Float16, Float32, Float64} in DoubleFloats at /Users/me/.julia/packages/DoubleFloats/h3HrU/src/math/prearith/prearith.jl:44

help?> frexp
search: frexp firstindex

  frexp(val)

  Return (x,exp) such that x has a magnitude in the interval [1/2, 1) or 0, and val is equal to x
  \times 2^{exp}.

mcabbott avatar May 10 '22 17:05 mcabbott

frexp and ldexp are dual functions in this sense ldexp(frexp(x)...) == x. That relationship is behind how these functions are used to examine and modify values. As DoubleFloats are pairs, high-part and low-part (most significant part and the lesser significant part), by defining frexp to operate on each part separately, and having ldexp do what it does normally just twice over, consistency is kept and, really, doing some numerical manipulation is easier.

The reason that exponent exploded is seen by trying exponent(0.0) in the REPL. The low-part of Double64(143.5) is zero -- and that brought about the issue. Thank you for raising it. Reasonably, that should be trapped. I guess returning an exponent of 1 is as reasonable as any integer (0.0^0 == 1, and not 0 in the REPL, it is undefined according to wolfram alpha, and 1 according to Maple).

JeffreySarnoff avatar May 12 '22 20:05 JeffreySarnoff

I was surprised by exponent(0.0) because I expected it to read out the bits -- while the mathematical answer is ambiguous as you say, there any solution will do, and the bitstring must have one. Maybe I'm not sure what the function is for, and in what context you'd be sure not to hit the error, even for Float64. (The docstring does not describe it as a bit-readout.)

But frexp is more what I was looking for. Its docstring makes clear that at 0.0 or NaN, any power would be fine. Here it seems more surprising to return something other than Float, Int for frexp. If you want to do something with its result before going back through ldexp, then generic code accepting ::AbstractFloat could just work, and would I think preserve the accuracy, but won't work with tuples. The linked PR seems like a perfect example of this (although I did not check accuracy). The returned float would (I think) want to have exactly the same low-part as the original.

I wonder if the version returning two tuples, for people who are sure they have a DoubleFloat, should have another name, rather than extending Base.frexp?

mcabbott avatar May 12 '22 21:05 mcabbott

I see what you mean. It could work as you explain with a helper function available to map the 2-tuple form into two 2-tuples for more exacting work. ...

JeffreySarnoff avatar May 12 '22 21:05 JeffreySarnoff