[SciPy-User] Conversion from 32bits IEEE floats to IBM floats

[Matthieu Brucher]

Thanks a lot Robert!
I didn't know about the view() method, interesting.

Matthieu

2011/6/9 Robert Kern <robert.kern at gmail.com>

> On Thu, Jun 9, 2011 at 08:17, Matthieu Brucher
> <matthieu.brucher at gmail.com> wrote:
> > Hi,
> > I wondered if anyone had a conversion routine for 32 bits IEEE floats in
> an
> > array to IBM floats (stored in a 4 bytes integers). I have a routine for
> > doing the opposite, but not IEEE->IBM.
> > There are codes in C or other languages, but the trick is starting from
> an
> > array (don't know if it can be reinterpreted as an integer array easily).
>
> This is almost certainly not the most elegant, but it seems to work for me:
>
> import numpy as np
>
> def ibm2ieee(ibm):
>    """ Converts an IBM floating point number into IEEE format. """
>
>    sign = ibm >> 31 & 0x01
>
>    exponent = ibm >> 24 & 0x7f
>
>    mantissa = ibm & 0x00ffffff
>    mantissa = (mantissa * np.float32(1.0)) / pow(2, 24)
>
>    ieee = (1 - 2 * sign) * mantissa * np.power(np.float32(16.0), exponent -
> 64)
>
>    return ieee
>
> def ieee2ibm(ieee):
>    ieee = ieee.astype(np.float32)
>    expmask = 0x7f800000
>    signmask = 0x80000000
>    mantmask = 0x7fffff
>    asint = ieee.view('i4')
>    signbit = asint & signmask
>    exponent = ((asint & expmask) >> 23) - 127
>    # The IBM 7-bit exponent is to the base 16 and the mantissa is presumed
> to
>    # be entirely to the right of the radix point. In contrast, the IEEE
>    # exponent is to the base 2 and there is an assumed 1-bit to the left of
> the
>    # radix point.
>    exp16 = ((exponent+1) // 4)
>    exp_remainder = (exponent+1) % 4
>    exp16 += exp_remainder != 0
>    downshift = np.where(exp_remainder, 4-exp_remainder, 0)
>    ibm_exponent = np.clip(exp16 + 64, 0, 127)
>    expbits = ibm_exponent << 24
>    # Add the implicit initial 1-bit to the 23-bit IEEE mantissa to get the
>    # 24-bit IBM mantissa. Downshift it by the remainder from the exponent's
>    # division by 4. It is allowed to have up to 3 leading 0s.
>    ibm_mantissa = ((asint & mantmask) | 0x800000) >> downshift
>    # Special-case 0.0
>    ibm_mantissa = np.where(ieee, ibm_mantissa, 0)
>    expbits = np.where(ieee, expbits, 0)
>    return signbit | expbits | ibm_mantissa
>
> --
> Robert Kern
>
> "I have come to believe that the whole world is an enigma, a harmless
> enigma that is made terrible by our own mad attempt to interpret it as
> though it had an underlying truth."
>   -- Umberto Eco
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>



-- 
Information System Engineer, Ph.D.
Blog: http://matt.eifelle.com
LinkedIn: http://www.linkedin.com/in/matthieubrucher
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20110610/a98cbdba/attachment.html>