cpl_vax.cpp

/******************************************************************************
 *
 * Project:  CPL
 * Purpose:  Convert between VAX and IEEE floating point formats
 * Author:   Frank Warmerdam, warmerdam@pobox.com
 *
 ******************************************************************************
 * Copyright (c) 2000, Avenza Systems Inc, http://www.avenza.com/
 *
 * SPDX-License-Identifier: MIT
 ****************************************************************************/

#include "cpl_port.h"
#include "cpl_vax.h"

namespace
{
typedef struct dbl
{
    // cppcheck-suppress unusedStructMember
    GUInt32 hi;
    // cppcheck-suppress unusedStructMember
    GUInt32 lo;
} double64_t;
}  // namespace

/************************************************************************/
/*                          CPLVaxToIEEEDouble()                        */
/************************************************************************/

void CPLVaxToIEEEDouble(void *dbl)

{
    double64_t dt;
    GUInt32 sign;
    int exponent;
    GUInt32 rndbits;

    /* -------------------------------------------------------------------- */
    /*      Arrange the VAX double so that it may be accessed by a          */
    /*      double64_t structure, (two GUInt32s).                           */
    /* -------------------------------------------------------------------- */
    {
        const unsigned char *src = static_cast<const unsigned char *>(dbl);
        unsigned char dest[8];
#ifdef CPL_LSB
        dest[2] = src[0];
        dest[3] = src[1];
        dest[0] = src[2];
        dest[1] = src[3];
        dest[6] = src[4];
        dest[7] = src[5];
        dest[4] = src[6];
        dest[5] = src[7];
#else
        dest[1] = src[0];
        dest[0] = src[1];
        dest[3] = src[2];
        dest[2] = src[3];
        dest[5] = src[4];
        dest[4] = src[5];
        dest[7] = src[6];
        dest[6] = src[7];
#endif
        memcpy(&dt, dest, 8);
    }

    /* -------------------------------------------------------------------- */
    /*      Save the sign of the double                                     */
    /* -------------------------------------------------------------------- */
    sign = dt.hi & 0x80000000;

    /* -------------------------------------------------------------------- */
    /*      Adjust the exponent so that we may work with it                 */
    /* -------------------------------------------------------------------- */
    exponent = (dt.hi >> 23) & 0x000000ff;

    if (exponent)
        exponent = exponent - 129 + 1023;

    /* -------------------------------------------------------------------- */
    /*      Save the bits that we are discarding so we can round properly   */
    /* -------------------------------------------------------------------- */
    rndbits = dt.lo & 0x00000007;

    dt.lo = dt.lo >> 3;
    dt.lo = (dt.lo & 0x1fffffff) | (dt.hi << 29);

    if (rndbits)
        dt.lo = dt.lo | 0x00000001;

    /* -------------------------------------------------------------------- */
    /*      Shift the hi-order int over 3 and insert the exponent and sign  */
    /* -------------------------------------------------------------------- */
    dt.hi = dt.hi >> 3;
    dt.hi = dt.hi & 0x000fffff;
    dt.hi = dt.hi | (static_cast<GUInt32>(exponent) << 20) | sign;

#ifdef CPL_LSB
    /* -------------------------------------------------------------------- */
    /*      Change the number to a byte swapped format                      */
    /* -------------------------------------------------------------------- */
    const unsigned char *src = reinterpret_cast<const unsigned char *>(&dt);
    unsigned char *dest = static_cast<unsigned char *>(dbl);

    memcpy(dest + 0, src + 4, 4);
    memcpy(dest + 4, src + 0, 4);
#else
    memcpy(dbl, &dt, 8);
#endif
}

/************************************************************************/
/*                         CPLIEEEToVaxDouble()                         */
/************************************************************************/

void CPLIEEEToVaxDouble(void *dbl)

{
    double64_t dt;

#ifdef CPL_LSB
    {
        const GByte *src = static_cast<const GByte *>(dbl);
        GByte dest[8];

        dest[0] = src[4];
        dest[1] = src[5];
        dest[2] = src[6];
        dest[3] = src[7];
        dest[4] = src[0];
        dest[5] = src[1];
        dest[6] = src[2];
        dest[7] = src[3];
        memcpy(&dt, dest, 8);
    }
#else
    memcpy(&dt, dbl, 8);
#endif

    GInt32 sign = dt.hi & 0x80000000;
    GInt32 exponent = dt.hi >> 20;
    exponent = exponent & 0x000007ff;

    /* -------------------------------------------------------------------- */
    /*      An exponent of zero means a zero value.                         */
    /* -------------------------------------------------------------------- */
    if (exponent)
        exponent = exponent - 1023 + 129;

    /* -------------------------------------------------------------------- */
    /*      In the case of overflow, return the largest number we can       */
    /* -------------------------------------------------------------------- */
    if (exponent > 255)
    {
        GByte dest[8];

        if (sign)
            dest[1] = 0xff;
        else
            dest[1] = 0x7f;

        dest[0] = 0xff;
        dest[2] = 0xff;
        dest[3] = 0xff;
        dest[4] = 0xff;
        dest[5] = 0xff;
        dest[6] = 0xff;
        dest[7] = 0xff;
        memcpy(dbl, dest, 8);

        return;
    }

    /* -------------------------------------------------------------------- */
    /*      In the case of of underflow return zero                         */
    /* -------------------------------------------------------------------- */
    else if ((exponent < 0) || (exponent == 0 && sign == 0))
    {
        memset(dbl, 0, 8);

        return;
    }
    else
    {
        /* --------------------------------------------------------------------
         */
        /*          Shift the fraction 3 bits left and set the exponent and
         * sign*/
        /* --------------------------------------------------------------------
         */
        dt.hi = dt.hi << 3;
        dt.hi = dt.hi | (dt.lo >> 29);
        dt.hi = dt.hi & 0x007fffff;
        dt.hi = dt.hi | (exponent << 23) | sign;

        dt.lo = dt.lo << 3;
    }

    /* -------------------------------------------------------------------- */
    /*      Convert the double back to VAX format                           */
    /* -------------------------------------------------------------------- */
    const GByte *src = reinterpret_cast<GByte *>(&dt);

#ifdef CPL_LSB
    GByte *dest = static_cast<GByte *>(dbl);
    memcpy(dest + 2, src + 0, 2);
    memcpy(dest + 0, src + 2, 2);
    memcpy(dest + 6, src + 4, 2);
    memcpy(dest + 4, src + 6, 2);
#else
    GByte dest[8];
    dest[1] = src[0];
    dest[0] = src[1];
    dest[3] = src[2];
    dest[2] = src[3];
    dest[5] = src[4];
    dest[4] = src[5];
    dest[7] = src[6];
    dest[6] = src[7];
    memcpy(dbl, dest, 8);
#endif
}

//////////////////////////////////////////////////////////////////////////
/// Below code is adapted from Public Domain VICAR project
/// https://github.com/nasa/VICAR/blob/master/vos/rtl/source/conv_vax_ieee_r.c
//////////////////////////////////////////////////////////////////////////

static void real_byte_swap(const unsigned char from[4], unsigned char to[4])
{
    to[0] = from[1];
    to[1] = from[0];
    to[2] = from[3];
    to[3] = from[2];
}

/* Shift x[1]..x[3] right one bit by bytes, don't bother with x[0] */
#define SHIFT_RIGHT(x)                                                         \
    {                                                                          \
        x[3] = ((x[3] >> 1) & 0x7F) | ((x[2] << 7) & 0x80);                    \
        x[2] = ((x[2] >> 1) & 0x7F) | ((x[1] << 7) & 0x80);                    \
        x[1] = (x[1] >> 1) & 0x7F;                                             \
    }

/* Shift x[1]..x[3] left one bit by bytes, don't bother with x[0] */
#define SHIFT_LEFT(x)                                                          \
    {                                                                          \
        x[1] = ((x[1] << 1) & 0xFE) | ((x[2] >> 7) & 0x01);                    \
        x[2] = ((x[2] << 1) & 0xFE) | ((x[3] >> 7) & 0x01);                    \
        x[3] = (x[3] << 1) & 0xFE;                                             \
    }

/************************************************************************/
/* Convert between IEEE and Vax single-precision floating point.        */
/* Both formats are represented as:                                     */
/* (-1)^s * f * 2^(e-bias)                                              */
/* where s is the sign bit, f is the mantissa (see below), e is the     */
/* exponent, and bias is the exponent bias (see below).                 */
/* There is an assumed leading 1 on the mantissa (except for IEEE       */
/* denormalized numbers), but the placement of the binary point varies. */
/*                                                                      */
/* IEEE format:    seeeeeee efffffff 8*f 8*f                            */
/*        where e is exponent with bias of 127 and f is of the          */
/*        form 1.fffff...                                               */
/* Special cases:                                                       */
/*    e=255, f!=0:        NaN (Not a Number)                            */
/*    e=255, f=0:        Infinity (+/- depending on s)                  */
/*    e=0, f!=0:        Denormalized numbers, of the form               */
/*                (-1)^s * (0.ffff) * 2^(-126)                          */
/*    e=0, f=0:        Zero (can be +/-)                                */
/*                                                                      */
/* VAX format:    seeeeeee efffffff 8*f 8*f                             */
/*        where e is exponent with bias of 128 and f is of the          */
/*        form .1fffff...                                               */
/* Byte swapping: Note that the above format is the logical format,     */
/*        which can be represented as bytes SE1 E2F1 F2 F3.             */
/*        The actual order in memory is E2F1 SE1 F3 F2 (which is        */
/*        two half-word swaps, NOT a full-word swap).                   */
/* Special cases:                                                       */
/*    e=0, s=0:        Zero (no +/-)                                    */
/*    e=0, s=1:        Invalid, causes Reserved Operand error           */
/*                                                                      */
/* The same code works on all byte-order machines because only byte     */
/* operations are performed.  It could perhaps be done more efficiently */
/* on a longword basis, but then the code would be byte-order dependent.*/
/* MAKE SURE any mods will work on either byte order!!!                 */
/************************************************************************/

/************************************************************************/
/* This routine will convert VAX F floating point values to IEEE        */
/* single precision floating point.                                     */
/************************************************************************/

static void vax_ieee_r(const unsigned char *from, unsigned char *ieee)
{
    unsigned char vaxf[4];
    unsigned char exp;

    real_byte_swap(from, vaxf); /* Put bytes in rational order */
    memcpy(ieee, vaxf, 4);      /* Since most bits are the same */

    exp = ((vaxf[0] << 1) & 0xFE) | ((vaxf[1] >> 7) & 0x01);

    if (exp == 0)
    { /* Zero or invalid pattern */
        if (vaxf[0] & 0x80)
        {                   /* Sign bit set, which is illegal for VAX */
            ieee[0] = 0x7F; /* IEEE NaN */
            ieee[1] = 0xFF;
            ieee[2] = 0xFF;
            ieee[3] = 0xFF;
        }
        else
        { /* Zero */
            ieee[0] = ieee[1] = ieee[2] = ieee[3] = 0;
        }
    }

    else if (exp >= 3)
    { /* Normal case */
        exp -= 2;
        ieee[0] =
            (vaxf[0] & 0x80) | ((exp >> 1) & 0x7F); /* remake sign + exponent */
    } /* Low bit of exp can't change, so don't bother w/it */

    else if (exp == 2)
    {                                      /* Denormalize the number */
        SHIFT_RIGHT(ieee);                 /* Which means shift right 1, */
        ieee[1] = (ieee[1] & 0x3F) | 0x40; /* Add suppressed most signif bit, */
        ieee[0] = vaxf[0] & 0x80; /* and set exponent to 0 (preserving sign) */
    }

    else
    {                      /* Exp==1, denormalize again */
        SHIFT_RIGHT(ieee); /* Like above but shift by 2 */
        SHIFT_RIGHT(ieee);
        ieee[1] = (ieee[1] & 0x1F) | 0x20;
        ieee[0] = vaxf[0] & 0x80;
    }

#ifdef CPL_LSB
    CPL_SWAP32PTR(ieee);
#endif
}

/************************************************************************/
/* This routine will convert IEEE single precision floating point       */
/* values to VAX F floating point.                                      */
/************************************************************************/

static void ieee_vax_r(unsigned char *ieee, unsigned char *to)
{
    unsigned char vaxf[4];
    unsigned char exp;

#ifdef CPL_LSB
    CPL_SWAP32PTR(ieee);
#endif

    memcpy(vaxf, ieee, 4); /* Since most bits are the same */

    exp = ((ieee[0] << 1) & 0xFE) | ((ieee[1] >> 7) & 0x01);

    /* Exponent 255 means NaN or Infinity, exponent 254 is too large for */
    /* VAX notation.  In either case, set to sign * highest possible number */

    if (exp == 255 || exp == 254)
    { /* Infinity or NaN or too big */
        vaxf[0] = 0x7F | (ieee[0] & 0x80);
        vaxf[1] = 0xFF;
        vaxf[2] = 0xFF;
        vaxf[3] = 0xFF;
    }

    else if (exp != 0)
    { /* Normal case */
        exp += 2;
        vaxf[0] =
            (ieee[0] & 0x80) | ((exp >> 1) & 0x7F); /* remake sign + exponent */
    } /* Low bit of exp can't change, so don't bother w/it */

    else
    { /* exp == 0, zero or denormalized number */
        if (ieee[1] == 0 && ieee[2] == 0 && ieee[3] == 0)
        { /* +/- 0 */
            vaxf[0] = vaxf[1] = vaxf[2] = vaxf[3] = 0;
        }
        else
        { /* denormalized number */
            if (ieee[1] & 0x40)
            {                                      /* hi bit set (0.1ffff) */
                SHIFT_LEFT(vaxf);                  /* Renormalize */
                vaxf[1] = vaxf[1] & 0x7F;          /* Set vax exponent to 2 */
                vaxf[0] = (ieee[0] & 0x80) | 0x01; /* sign, exponent==2 */
            }
            else if (ieee[1] & 0x20)
            {                     /* next bit set (0.01ffff) */
                SHIFT_LEFT(vaxf); /* Renormalize */
                SHIFT_LEFT(vaxf);
                vaxf[1] = vaxf[1] | 0x80; /* Set vax exponent to 1 */
                vaxf[0] = ieee[0] & 0x80; /* sign, exponent==1 */
            }
            else
            { /* Number too small for VAX */
                vaxf[0] = vaxf[1] = vaxf[2] = vaxf[3] = 0; /* so set to 0 */
            }
        }
    }

    real_byte_swap(vaxf, to); /* Put bytes in weird VAX order */
}

void CPLVaxToIEEEFloat(void *f)
{
    unsigned char res[4];
    vax_ieee_r(static_cast<const unsigned char *>(f), res);
    memcpy(f, res, 4);
}

void CPLIEEEToVaxFloat(void *f)
{
    unsigned char res[4];
    ieee_vax_r(static_cast<unsigned char *>(f), res);
    memcpy(f, res, 4);
}