If you want to keep the bit-level approach (as opposed to doing a final computation in Float64), then you either use leading_zeros or modify your loop. You also need to adjust the biases. The function now looks similar the previous bit-level code:
function juliohm2(b::UInt32)
s = b & 0x80000000
e = b & 0x7f000000
f = b & 0x00ffffff
# if fraction is zero, value is zero
iszero(f) && return zero(Float64)
# shift sign bit to the rightmost position
sieee = Int64(s >> 31)
# shift exponent bits to the right (3 bytes = 24 bits)
# and multiply by 4 to convert from base 16 to base 2
# hence, the resulting shift is 24 - 2 = 22 bits
eieee = Int64(e >> 22)
eieee += 1023 - 4*64 # adjust bias
# normalize fraction until the leading digit is equal to 1,
# and adjust the exponent accordingly (opposite shift)
while iszero(f & 0x00800000)
f <<= 1
eieee -= 1
end
f = (f << 1) & 0x00ffffff # drop leading bit
eieee -= 1 # account for dropped leading bit
fieee = UInt64(f) << (52 - 24) # pad with zeros to the right
reinterpret(Float64, sieee << 63 | eieee << 52 | fieee)
end
Thank you @matthias314 , I am learning a lot about these bit manipulations from your code snippets. I will try to simplify the code further based on the last version you shared.
@matthias314 regarding this line, shouldnât it just be 1023 - 64? Where is the 4 coming from?
My understanding from your previous comments:
# exponent of 64 means 0 in the IBM standard, and
# exponent of 1023 means 0 in the IEEE standard
# hence, the bias adjustment is -64 + 1023 = 959
eieee += 959
function Base.convert(::Type{Float64}, x::IBMFloat32)
# extract sign, exponent, and fraction bits
b = reinterpret(UInt32, x)
s = b & 0x80000000
e = b & 0x7f000000
f = b & 0x00ffffff
# if fraction is zero, value is zero
iszero(f) && return zero(Float64)
# position sign bit at the leftmost bit of 8 bytes
sieee = Int64(s) << 32
# shift exponent bits to the right (3 bytes = 24 bits)
# and multiply by 4 to convert from base 16 to base 2
# hence, the resulting shift is 24 - 2 = 22 bits
eieee = Int64(e >> 22)
# exponent of 64 in base 16 means 0 in the IBM standard,
# and exponent of 1023 in base 2 means 0 in the IEEE standard
# hence, the bias adjustment is 1023 - 4*64 = 767
eieee += 767
# shift fraction bits to the left until the leading digit is 1
# and then drop the leading 1 as it is implicit in IEEE
shift = leading_zeros(f) - 8 # ignore sign and exponent bits
fieee = UInt64((f << (shift + 1)) & 0x00ffffff)
eieee -= (shift + 1)
# pad fraction with zeros to the right, and since IEEE
# uses 52 bits for the fraction and IBM uses 24 bits,
# we need to shift 52 - 24 = 28 bits to the left
fieee <<= 28
# now that the exponent bits have been adjusted, shift
# them to their final position in the IEEE representation
eieee <<= 52
# reinterpret bits as Float64
reinterpret(Float64, sieee | eieee | fieee)
end
Hi @Palli, thank you for sharing the links, I am aware of them already. I am writing a new package that supports all SEGY revisions from revision 0.0 (before 1975) to revision 2.1 (2023). It is in very good shape already, except for the IBM float that is almost there.
@matthias314 I am very happy with this implementation, thank you for the incredible help. Do you have any additional suggestion to improve readability or performance?
I understand that the only place where we are not doing explicit bit manipulation is the bias adjustment with eieee += 767. Is there an equivalent adjustment in terms of bit shifts?
Just a minor comment which may or may not be obvious to everyone.
If you have two (reasonably fast) implementations for this conversion and want to verify that they are doing the same thing, you can just run through all the 2^32 bit patterns.
You might find my package VaxData.jl helpful as a reference, since it sounds like you have similar goals supporting (reading, converting, etc) historic floating point representations.
The function ibm_to_ieee_64 is located in the file fp.c that is included in the zip file 3rd_Ed__NumbersAndComputers.zip in the repo, and the repo license file specifies the MIT license, so I think you are safe to publish the code.
For future reference: the solution I marked in this thread is ~10x faster than the one shared in the book (checked with BenchmarkTools.jl).
Thank you everyone for the help and for the additional links. This thread will certainly help future readers interested in this old floating point type.
Even though this issue has been marked as solved, Iâm going to add one more bitâŚ
In May 2024, Cleve Moler, the famous numerical analyst and founder of Matlab, posted a blog entry on this same topic. He included a historical review, annotated sample Matlab codes, and a link to more robust Matlab codes (easily converted to Julia) that convert to/from IBM System/360 hexadecimal format, in both single and double precision. I obtained his permission to use these codes in other open-source software. So that is another good source of information on this topic.
