Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

exponent and significand functions #148

Open
mcabbott opened this issue May 9, 2022 · 4 comments
Open

exponent and significand functions #148

mcabbott opened this issue May 9, 2022 · 4 comments

Comments

@mcabbott
Copy link

mcabbott commented May 9, 2022

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 JuliaStats/LogExpFunctions.jl#48 could work on AbstractFloat.

@mcabbott
Copy link
Author

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}.

@JeffreySarnoff
Copy link
Member

JeffreySarnoff commented May 12, 2022

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).

@mcabbott
Copy link
Author

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?

@JeffreySarnoff
Copy link
Member

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. ...

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants