Function Dragon4_PrintFloat_IEEE_binary128

numpy/core/src/multiarray/dragon4.c:2718–2800

* IEEE binary128 floating-point format * * sign: 1 bit * exponent: 15 bits * mantissa: 112 bits * * Currently binary128 format exists on only a few CPUs, such as on the POWER9 * arch or aarch64. Because of this, this code has not been extensively tested. * I am not sure if the arch also supports uint128, and C does not seem to * support int128 literals. So we use uint64 to do manip

 * support int128 literals. So we use uint64 to do manipulation.
 */
static npy_uint32
Dragon4_PrintFloat_IEEE_binary128(
    Dragon4_Scratch *scratch, FloatVal128 val128, Dragon4_Options *opt)
{
    char *buffer = scratch->repr;
    const npy_uint32 bufferSize = sizeof(scratch->repr);
    BigInt *bigints = scratch->bigints;

    npy_uint32 floatExponent, floatSign;

    npy_uint64 mantissa_hi, mantissa_lo;
    npy_int32 exponent;
    npy_uint32 mantissaBit;
    npy_bool hasUnequalMargins;
    char signbit = '\0';

    mantissa_hi = val128.hi & bitmask_u64(48);
    mantissa_lo = val128.lo;
    floatExponent = (val128.hi >> 48) & bitmask_u32(15);
    floatSign = val128.hi >> 63;

    /* output the sign */
    if (floatSign != 0) {
        signbit = '-';
    }
    else if (opt->sign) {
        signbit = '+';
    }

    /* if this is a special value */
    if (floatExponent == bitmask_u32(15)) {
        npy_uint64 mantissa_zero = mantissa_hi == 0 && mantissa_lo == 0;
        return PrintInfNan(buffer, bufferSize, !mantissa_zero, 16, signbit);
    }
    /* else this is a number */

    /* factor the value into its parts */
    if (floatExponent != 0) {
        /*
         * normal
         * The floating point equation is:
         *  value = (1 + mantissa/2^112) * 2 ^ (exponent-16383)
         * We convert the integer equation by factoring a 2^112 out of the
         * exponent
         *  value = (1 + mantissa/2^112) * 2^112 * 2 ^ (exponent-16383-112)
         *  value = (2^112 + mantissa) * 2 ^ (exponent-16383-112)
         * Because of the implied 1 in front of the mantissa we have 112 bits of
         * precision.
         *   m = (2^112 + mantissa)
         *   e = (exponent-16383+1-112)
         *
         *   Adding 2^112 to the mantissa is the same as adding 2^48 to the hi
         *   64 bit part.
         */
        mantissa_hi         = (1ull << 48) | mantissa_hi;
        /* mantissa_lo is unchanged */
        exponent            = floatExponent - 16383 - 112;
        mantissaBit         = 112;