VAX F_FLOAT and D_FLOAT to IEEE T_FLOAT and S_FLOAT (double)

Not quite a nice first post to introduce myself or so, but a nice uncommented Python code that converts VAX floating point formats (8 bytes D_FLOAT and 4 bytes F_FLOAT) to the widespread IEEE float formats (8 bytes T_FLOAT, aka "double" and the regular 4 bytes S_FLOAT).

It has been heavily inspired from the libvaxdata library (especially for the constants).

Some notes about the platform used:

• We assume here that it is little endian.
• Both 32 and 64 bits platforms are supported, thanks to the bitandnot function (equivalent to "v1 & ~v2" in C) and the use of arithmetics instead of bitwise operators in some cases.
  

OK, the code is not documented but still, the variable names are quite explicit. 


It is nicely documented in the following document: http://pubs.usgs.gov/of/2005/1424/of2005-1424_v1.2.pdf

Refresher about the floating point representations (source: Code Project article):


VAX single precision floating point: 

SEF :       S        EEEEEEEE        FFFFFFF        FFFFFFFF        FFFFFFFF
bits :      1        2      9        10                                   32
bytes :     byte2           byte1                   byte4           byte3

VAX double precision floating point: 

SEF :       S        EEEEEEEE        FFFFFFF        FFFFFFFF        FFFFFFFF
bits :      1        2      9        10                                   32
bytes :     byte2           byte1                   byte4           byte3
frctn.:                    L2                     L1

IEEE single precision floating point: 

SEF :    S        EEEEEEEE        FFFFFFF        FFFFFFFF        FFFFFFFF
bits :   1        2      9        10                                   32
bytes :  byte1           byte2                   byte3           byte4

IEEE double precision floating point:

SEF:   S     EEEEEEE EEEE  FFFF FFFFFFFF FFFFFFFF FFFFFFFF FFFFFFFF FFFFFFFF FFFFFFFF
bits:  1     2          12 13                                                      64
bytes: byte1         byte2      byte3    byte4    byte5    byte6    byte7    byte8
frctn.:                    L1                     L2

General encoding:

V = (-1)^S * M * A^(E - B)
M = C + F

IEEE single float :     A = 2   B = 127   C = 1
IEEE double float :     A = 2   B = 1023  C = 1
IBM  single float :     A = 16  B = 64    C = 0
IBM  double float :     A = 16  B = 64    C = 0
VAX  single float :     A = 2   B = 128   C = 0.5
VAX  double float :     A = 2   B = 128   C = 0.5

And hey, small trick: note that the 4 first bytes of a VAX D_FLOAT represent a VAX F_FLOAT. Easy to switch from one format to another, simply ignore the second half of your 8 bytes array.


import struct

SIGN_BIT = 0x80000000

VAX_D_SIGN_BIT = SIGN_BIT
VAX_D_EXPONENT_MASK = 0x7F800000
VAX_D_EXPONENT_SIZE = 8
VAX_D_EXPONENT_BIAS = (1 << ( VAX_D_EXPONENT_SIZE - 1))
VAX_D_MANTISSA_MASK = 0x007FFFFF
VAX_D_MANTISSA_SIZE = 23
VAX_D_HIDDEN_BIT    = (1 << VAX_D_MANTISSA_SIZE)

VAX_F_MANTISSA_MASK = 0x007FFFFF
VAX_F_MANTISSA_SIZE = 23
VAX_F_HIDDEN_BIT = (1 << VAX_F_MANTISSA_SIZE)
VAX_F_EXPONENT_SIZE = 8
VAX_F_EXPONENT_BIAS = (1 << (VAX_F_EXPONENT_SIZE - 1))
VAX_F_EXPONENT_MASK = 0x7F800000

IEEE_T_MANTISSA_SIZE = 20
IEEE_T_EXPONENT_SIZE = 11
IEEE_T_EXPONENT_BIAS = ((1 << (IEEE_T_EXPONENT_SIZE - 1)) - 1)

IEEE_S_EXPONENT_SIZE = 8
IEEE_S_EXPONENT_BIAS = ((1 << (IEEE_S_EXPONENT_SIZE - 1)) - 1)
IEEE_S_MANTISSA_SIZE = 23

def bitandnot(v1, v2):
    return (((v1 / 65536) & (0xffff - (v2 / 65536))) * 65536
          + ((v1 % 65536) & (0xffff - (v2 % 65536))))

def vd8tod(stream):
    EXPONENT_RIGHT_SHIFT = (VAX_D_MANTISSA_SIZE - IEEE_T_MANTISSA_SIZE)
    EXPONENT_ADJUSTMENT = (1 + VAX_D_EXPONENT_BIAS - IEEE_T_EXPONENT_BIAS)
    IN_PLACE_EXPONENT_ADJUSTMENT = (EXPONENT_ADJUSTMENT << IEEE_T_MANTISSA_SIZE)

    v1a, v1b, v2a, v2b = struct.unpack('HHHH', stream)

    v1 = v1b + 65536 * v1a
    v2 = v2b + 65536 * v2a

    if (v1 & VAX_D_EXPONENT_MASK) == 0:
        if (v1 & SIGN_BIT) == SIGN_BIT:
            raise
        res = '\x00' * 8
    else:
        i1 = (((v1 & SIGN_BIT) | (bitandnot(v1, SIGN_BIT) >> EXPONENT_RIGHT_SHIFT)) - IN_PLACE_EXPONENT_ADJUSTMENT )
        i2 = (((v1 << (32 - EXPONENT_RIGHT_SHIFT)) | (v2 >> EXPONENT_RIGHT_SHIFT))) & 0xffff

        res = struct.pack('LL', i2, i1);

    return struct.unpack('d', res)[0];

def vr4tof(stream):
    EXPONENT_ADJUSTMENT = (1 + VAX_F_EXPONENT_BIAS - IEEE_S_EXPONENT_BIAS)
    IN_PLACE_EXPONENT_ADJUSTMENT = (EXPONENT_ADJUSTMENT << IEEE_S_MANTISSA_SIZE)

    v1a, v1b = struct.unpack('HH', stream)
    v1 = v1b + 65536 * v1a

    e = v1 & VAX_F_EXPONENT_MASK

    if e == 0:
        if v1 & SIGN_BIT == SIGN_BIT:
            raise
        res = '\x00' * 4
    else:
        e >>= VAX_F_MANTISSA_SIZE
        e -= EXPONENT_ADJUSTMENT

        if e > 0:
            res = v1 - IN_PLACE_EXPONENT_ADJUSTMENT
        else:
            res = (v1 & SIGN_BIT) | ((VAX_F_HIDDEN_BIT | (v1 & VAX_F_MANTISSA_MASK)) >> (1 - e))

        res = struct.pack('L', res)

    return struct.unpack('f', res)[0]

if __name__ == "__main__":
    stream = '\x40\xc1\x00\x00\x00\x00\x00\x00' # -3.0
    print vd8tod(stream)
    stream = '\x9b\x41\x33\x33\x00\x00\x00\x00' # approx 4.85
    print vd8tod(stream)
    print vr4tof(stream[:4])