Standardizing IBM-Float as a Julia minipackage: any interest?

[jpjones76]

Hi everyone,

I'm the creator/maintainer of the "other" SeisIO package, and I'm trying to standardize IBM-Float in Julia (same reason you use it, really -- we all need full SEG Y support). Are you interested in turning IBMFloat support into a separate package with me? Should be short, not much work required.

At issue: the solution you use, from JuliaSeis.jl (radix shift + reinterpret), yields wrong answers outside the range 1.0f-45 ≤ abs(x) ≤ 3.4028235e38; the range of single-precision IBM-Float is 5.39761e-79 ≤ abs(x) ≤ 7.23701e75. So it works for most seismic data, but it can fail.

Proposed solution: convert IBMFloat32 ==> Float64 and IBMFloat64 ==> BigFloat to support their full ranges.

Below is a trial function that correctly parses all test bit representations in https://en.wikipedia.org/wiki/IBM_hexadecimal_floating_point. However, it's a factor of four slower than JuliaSeis.jl because of the last step. (slower as presented here the function below doesn't define a new primitive Type)

I haven't worked out the radix shifts, etc., to convert IBMFloat32 ==> Float64, but I think that could be done by modifying your code.

Is this of any interest?

function ibmfloat(x::UInt32)
  local fra = ntoh(x)
  local sgn = UInt8((fra & 0x80000000)>>31)
  fra <<= 1
  local exp = Int16(fra >> 25) - Int16(64)
  fra <<= 7
  y = (sgn == 0x00 ? 1.0 : -1.0) * 16.0^exp * signed(fra >> 8)/16777216
  return y
end