Skip to content

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

@jpjones76

Description

@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

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions