QuadFloat.st
author Claus Gittinger <cg@exept.de>
Tue, 25 Jun 2019 14:28:51 +0200
changeset 5050 44fa8672d102
parent 5027 49e378bf72b7
child 5056 0744a9b601c5
permissions -rw-r--r--
#DOCUMENTATION by cg class: SharedQueue comment/format in: #next #nextWithTimeout:

"{ Encoding: utf8 }"

"{ Package: 'stx:libbasic2' }"

"{ NameSpace: Smalltalk }"

LimitedPrecisionReal variableByteSubclass:#QuadFloat
	instanceVariableNames:''
	classVariableNames:'QuadFloatZero QuadFloatOne Pi E Epsilon NaN PositiveInfinity
		NegativeInfinity Halfpi HalfpiNegative'
	poolDictionaries:''
	category:'Magnitude-Numbers'
!

!QuadFloat primitiveDefinitions!
%{
#if __POINTER_SIZE__ == 8
# define SUPPORT_QUADFLOAT
#endif

#if defined(__MSBFIRST__) || defined(__MSB_FIRST__) || defined(MSBFIRST)
# undef LITTLEENDIAN
#else
# define LITTLEENDIAN
#endif

#ifdef SUPPORT_QUADFLOAT
# include <stdint.h>

// The following code is adapted from the Softfloat-3e package,
// which includes the following license/copyright statement:
//
// License for Berkeley SoftFloat Release 3e
//
// John R. Hauser
// 2018 January 20
//
// The following applies to the whole of SoftFloat Release 3e as well as to
// each source file individually.
//
// Copyright 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018 The Regents of the
// University of California.  All rights reserved.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//
//  1. Redistributions of source code must retain the above copyright notice,
//     this list of conditions, and the following disclaimer.
//
//  2. Redistributions in binary form must reproduce the above copyright
//     notice, this list of conditions, and the following disclaimer in the
//     documentation and/or other materials provided with the distribution.
//
//  3. Neither the name of the University nor the names of its contributors
//     may be used to endorse or promote products derived from this software
//     without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS "AS IS", AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, ARE
// DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE FOR ANY
// DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
// (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
// LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
// ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
// THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//

/*----------------------------------------------------------------------------
| Software floating-point underflow tininess-detection mode.
*----------------------------------------------------------------------------*/
enum {
    softfloat_tininess_beforeRounding = 0,
    softfloat_tininess_afterRounding  = 1
};

/*----------------------------------------------------------------------------
| Software floating-point exception flags.
*----------------------------------------------------------------------------*/
enum {
    softfloat_flag_inexact   =  1,
    softfloat_flag_underflow =  2,
    softfloat_flag_overflow  =  4,
    softfloat_flag_infinite  =  8,
    softfloat_flag_invalid   = 16
};
/*----------------------------------------------------------------------------
| Software floating-point rounding mode.  (Mode "odd" is supported only if
| SoftFloat is compiled with macro 'SOFTFLOAT_ROUND_ODD' defined.)
*----------------------------------------------------------------------------*/
enum {
    softfloat_round_near_even   = 0,
    softfloat_round_minMag      = 1,
    softfloat_round_min         = 2,
    softfloat_round_max         = 3,
    softfloat_round_near_maxMag = 4,
    softfloat_round_odd         = 6
};
#define init_detectTininess softfloat_tininess_afterRounding

typedef unsigned char bool;
typedef double        float64_t;

#if defined(LITTLEENDIAN)
struct uint128          { uint64_t v0, v64; };
struct uint64_extra     { uint64_t extra, v; };
struct uint128_extra    { uint64_t extra; struct uint128 v; };
struct extFloat80M      { uint64_t signif; uint16_t signExp; };
#else
struct uint128          { uint64_t v64, v0; };
struct uint64_extra     { uint64_t v, extra; };
struct uint128_extra    { struct uint128 v; uint64_t extra; };
struct extFloat80M      { uint16_t signExp; uint64_t signif; };
#endif

union ui64_f64          { uint64_t ui; float64_t f; };
union ui128_f128        { struct uint128 ui; float128_t f; };
struct exp32_sig128     { int_fast32_t exp; struct uint128 sig; };

typedef struct extFloat80M        extFloat80_t;

/*----------------------------------------------------------------------------
| The bit pattern for a default generated 128-bit floating-point NaN.
*----------------------------------------------------------------------------*/
#define defaultNaNF128UI96_32 0xFFFF8000
#define defaultNaNF128UI64_32 0
#define defaultNaNF128UI32_32 0
#define defaultNaNF128UI0_32  0

// #define defaultNaNF128UI64_64 UINT64_C( 0x7FFF800000000000 )
// #define defaultNaNF128UI0_64  UINT64_C( 0 )

/*----------------------------------------------------------------------------
| The bit pattern for a default generated 128-bit floating-point +INF.
*----------------------------------------------------------------------------*/
#define defaultPInfF128UI96_32 0x7FFF0000
#define defaultPInfF128UI64_32 0
#define defaultPInfF128UI32_32 0
#define defaultPInfF128UI0_32  0
/*----------------------------------------------------------------------------
| The bit pattern for a default generated 128-bit floating-point -INF.
*----------------------------------------------------------------------------*/
#define defaultNInfF128UI96_32 0xFFFF0000
#define defaultNInfF128UI64_32 0
#define defaultNInfF128UI32_32 0
#define defaultNInfF128UI0_32  0

struct commonNaN {
    bool sign;
#ifdef LITTLEENDIAN
    uint64_t v0, v64;
#else
    uint64_t v64, v0;
#endif
};
struct exp16_sig64 { int_fast16_t exp; uint_fast64_t sig; };

#endif // SUPPORT_QUADFLOAT

%}
! !

!QuadFloat primitiveVariables!
%{
#ifdef SUPPORT_QUADFLOAT

uint_fast8_t softfloat_exceptionFlags;
uint_fast8_t softfloat_roundingMode = softfloat_round_near_even;
uint_fast8_t softfloat_detectTininess = init_detectTininess;
uint_fast8_t softfloat_exceptionFlags = 0;

const uint_least8_t softfloat_countLeadingZeros8[256] = {
    8, 7, 6, 6, 5, 5, 5, 5, 4, 4, 4, 4, 4, 4, 4, 4,
    3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
    2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
    2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
};
#endif // SUPPORT_QUADFLOAT

%}
! !

!QuadFloat primitiveFunctions!
%{
#ifdef SUPPORT_QUADFLOAT

#if defined(LITTLEENDIAN) || defined(__LSBFIRST__)
#define wordIncr 1
#define indexWord( total, n ) (n)
#define indexWordHi( total ) ((total) - 1)
#define indexWordLo( total ) 0
#define indexMultiword( total, m, n ) (n)
#define indexMultiwordHi( total, n ) ((total) - (n))
#define indexMultiwordLo( total, n ) 0
#define indexMultiwordHiBut( total, n ) (n)
#define indexMultiwordLoBut( total, n ) 0
#define INIT_UINTM4( v3, v2, v1, v0 ) { v0, v1, v2, v3 }
#else
#define wordIncr -1
#define indexWord( total, n ) ((total) - 1 - (n))
#define indexWordHi( total ) 0
#define indexWordLo( total ) ((total) - 1)
#define indexMultiword( total, m, n ) ((total) - 1 - (m))
#define indexMultiwordHi( total, n ) 0
#define indexMultiwordLo( total, n ) ((total) - (n))
#define indexMultiwordHiBut( total, n ) 0
#define indexMultiwordLoBut( total, n ) (n)
#define INIT_UINTM4( v3, v2, v1, v0 ) { v3, v2, v1, v0 }
#endif

#define signF64UI( a ) ((bool) ((uint64_t) (a)>>63))
#define expF64UI( a ) ((int_fast16_t) ((a)>>52) & 0x7FF)
#define fracF64UI( a ) ((a) & UINT64_C( 0x000FFFFFFFFFFFFF ))
#define softfloat_f64UIToCommonNaN( uiA, zPtr ) if ( ! ((uiA) & UINT64_C( 0x0008000000000000 )) ) softfloat_raiseFlags( softfloat_flag_invalid )

#define signF128UI64( a64 ) ((bool) ((uint64_t) (a64)>>63))
#define expF128UI64( a64 ) ((int_fast32_t) ((a64)>>48) & 0x7FFF)
#define fracF128UI64( a64 ) ((a64) & UINT64_C( 0x0000FFFFFFFFFFFF ))
#define packToF128UI64( sign, exp, sig64 ) (((uint_fast64_t) (sign)<<63) + ((uint_fast64_t) (exp)<<48) + (sig64))

#define signF128UI96( a96 ) ((bool) ((uint32_t) (a96)>>31))
#define expF128UI96( a96 ) ((int32_t) ((a96)>>16) & 0x7FFF)
#define fracF128UI96( a96 ) ((a96) & 0x0000FFFF)
#define packToF128UI96( sign, exp, sig96 ) (((uint32_t) (sign)<<31) + ((uint32_t) (exp)<<16) + (sig96))

#define isNaNF128UI( a64, a0 ) (((~(a64) & UINT64_C( 0x7FFF000000000000 )) == 0) && (a0 || ((a64) & UINT64_C( 0x0000FFFFFFFFFFFF ))))
#define softfloat_isSigNaNF128UI( uiA64, uiA0 ) ((((uiA64) & UINT64_C( 0x7FFF800000000000 )) == UINT64_C( 0x7FFF000000000000 )) && ((uiA0) || ((uiA64) & UINT64_C( 0x00007FFFFFFFFFFF ))))
#define softfloat_approxRecip32_1( a ) ((uint32_t) (UINT64_C( 0x7FFFFFFFFFFFFFFF ) / (uint32_t) (a)))
#define isInfF128UI( a64, a0 ) (((a64) & UINT64_C(0x7FFF800000000000)) == UINT64_C(0x7FFF000000000000))

#define signExtF80UI64( a64 ) ((bool) ((uint16_t) (a64)>>15))
#define expExtF80UI64( a64 ) ((a64) & 0x7FFF)
#define packToExtF80UI64( sign, exp ) ((uint_fast16_t) (sign)<<15 | (exp))

#define isNaNExtF80UI( a64, a0 ) ((((a64) & 0x7FFF) == 0x7FFF) && ((a0) & UINT64_C( 0x7FFFFFFFFFFFFFFF )))

/*----------------------------------------------------------------------------
| This function or macro is the same as 'softfloat_shortShiftLeftM' with
| 'size_words' = 5 (N = 160).
*----------------------------------------------------------------------------*/
#define softfloat_shortShiftLeft160M( aPtr, dist, zPtr ) softfloat_shortShiftLeftM( 5, aPtr, dist, zPtr )
/*----------------------------------------------------------------------------
| This function or macro is the same as 'softfloat_shiftRightJamM' with
| 'size_words' = 5 (N = 160).
*----------------------------------------------------------------------------*/
#define softfloat_shiftRightJam160M( aPtr, dist, zPtr ) softfloat_shiftRightJamM( 5, aPtr, dist, zPtr )
/*----------------------------------------------------------------------------
| This function or macro is the same as 'softfloat_shortShiftLeftM' with
| 'size_words' = 4 (N = 128).
*----------------------------------------------------------------------------*/
#define softfloat_shortShiftLeft128M( aPtr, dist, zPtr ) softfloat_shortShiftLeftM( 4, aPtr, dist, zPtr )
/*----------------------------------------------------------------------------
| This function or macro is the same as 'softfloat_shortShiftLeftM' with
| 'size_words' = 3 (N = 96).
*----------------------------------------------------------------------------*/
#define softfloat_shortShiftLeft96M( aPtr, dist, zPtr ) softfloat_shortShiftLeftM( 3, aPtr, dist, zPtr )

static inline void
softfloat_raiseFlags( uint_fast8_t flags ) {
    softfloat_exceptionFlags |= flags;
}

#if 1
static inline void
softfloat_commonNaNToF128M( uint32_t *zWPtr )
{
    zWPtr[indexWord( 4, 3 )] = defaultNaNF128UI96_32;
    zWPtr[indexWord( 4, 2 )] = defaultNaNF128UI64_32;
    zWPtr[indexWord( 4, 1 )] = defaultNaNF128UI32_32;
    zWPtr[indexWord( 4, 0 )] = defaultNaNF128UI0_32;
}
#else
# define softfloat_commonNaNToF128M( zWPtr ) \
{ \
    (zWPtr)[indexWord( 4, 3 )] = defaultNaNF128UI96_32; \
    (zWPtr)[indexWord( 4, 2 )] = defaultNaNF128UI64_32; \
    (zWPtr)[indexWord( 4, 1 )] = defaultNaNF128UI32_32; \
    (zWPtr)[indexWord( 4, 0 )] = defaultNaNF128UI0_32; \
}
#endif

/*----------------------------------------------------------------------------
| Assuming the unsigned integer formed from concatenating 'uiA64' and 'uiA0'
| has the bit pattern of an 80-bit extended floating-point NaN, converts
| this NaN to the common NaN form, and stores the resulting common NaN at the
| location pointed to by 'zPtr'.  If the NaN is a signaling NaN, the invalid
| exception is raised.
*----------------------------------------------------------------------------*/
#define softfloat_extF80UIToCommonNaN( uiA64, uiA0, zPtr ) if ( ! ((uiA0) & UINT64_C( 0x4000000000000000 )) ) softfloat_raiseFlags( softfloat_flag_invalid )

static inline void
softfloat_invalidF128M( uint32_t *zWPtr ) {
    softfloat_raiseFlags( softfloat_flag_invalid );
    softfloat_commonNaNToF128M( zWPtr );
}

static bool
f128M_isSignalingNaN( const float128_t *aPtr )
{
    const uint32_t *aWPtr;
    uint32_t uiA96;

    aWPtr = (const uint32_t *) aPtr;
    uiA96 = aWPtr[indexWordHi( 4 )];
    if ( (uiA96 & 0x7FFF8000) != 0x7FFF0000 ) return 0;
    return
	((uiA96 & 0x00007FFF) != 0)
	    || ((aWPtr[indexWord( 4, 2 )] | aWPtr[indexWord( 4, 1 )]
		     | aWPtr[indexWord( 4, 0 )])
		    != 0);
}

static void
softfloat_shortShiftRightJamM(
     uint_fast8_t size_words,
     const uint32_t *aPtr,
     uint_fast8_t dist,
     uint32_t *zPtr
 )
{
    uint_fast8_t uNegDist;
    unsigned int index, lastIndex;
    uint32_t partWordZ, wordA;

    uNegDist = -dist;
    index = indexWordLo( size_words );
    lastIndex = indexWordHi( size_words );
    wordA = aPtr[index];
    partWordZ = wordA>>dist;
    if ( partWordZ<<dist != wordA ) partWordZ |= 1;
    while ( index != lastIndex ) {
	wordA = aPtr[index + wordIncr];
	zPtr[index] = wordA<<(uNegDist & 31) | partWordZ;
	index += wordIncr;
	partWordZ = wordA>>dist;
    }
    zPtr[index] = partWordZ;

}

void
softfloat_shiftRightJamM(
     uint_fast8_t size_words,
     const uint32_t *aPtr,
     uint32_t dist,
     uint32_t *zPtr
 )
{
    uint32_t wordJam, wordDist, *ptr;
    uint_fast8_t i, innerDist;

    wordJam = 0;
    wordDist = dist>>5;
    if ( wordDist ) {
	if ( size_words < wordDist ) wordDist = size_words;
	ptr = (uint32_t *) (aPtr + indexMultiwordLo( size_words, wordDist ));
	i = wordDist;
	do {
	    wordJam = *ptr++;
	    if ( wordJam ) break;
	    --i;
	} while ( i );
	ptr = zPtr;
    }
    if ( wordDist < size_words ) {
	aPtr += indexMultiwordHiBut( size_words, wordDist );
	innerDist = dist & 31;
	if ( innerDist ) {
	    softfloat_shortShiftRightJamM(
		size_words - wordDist,
		aPtr,
		innerDist,
		(zPtr + indexMultiwordLoBut( size_words, wordDist ))
	    );
	    if ( ! wordDist ) goto wordJam;
	} else {
	    aPtr += indexWordLo( size_words - wordDist );
	    ptr = zPtr + indexWordLo( size_words );
	    for ( i = size_words - wordDist; i; --i ) {
		*ptr = *aPtr;
		aPtr += wordIncr;
		ptr += wordIncr;
	    }
	}
	ptr = zPtr + indexMultiwordHi( size_words, wordDist );
    }
    do {
	*ptr++ = 0;
	--wordDist;
    } while ( wordDist );
 wordJam:
    if ( wordJam ) zPtr[indexWordLo( size_words )] |= 1;

}

void
 softfloat_shortShiftLeftM(
     uint_fast8_t size_words,
     const uint32_t *aPtr,
     uint_fast8_t dist,
     uint32_t *zPtr
 )
{
    uint_fast8_t uNegDist;
    unsigned int index, lastIndex;
    uint32_t partWordZ, wordA;

    uNegDist = -dist;
    index = indexWordHi( size_words );
    lastIndex = indexWordLo( size_words );
    partWordZ = aPtr[index]<<dist;
    while ( index != lastIndex ) {
	wordA = aPtr[index - wordIncr];
	zPtr[index] = partWordZ | wordA>>(uNegDist & 31);
	index -= wordIncr;
	partWordZ = wordA<<dist;
    }
    zPtr[index] = partWordZ;
}

static void
softfloat_shiftLeftM(
     uint_fast8_t size_words,
     const uint32_t *aPtr,
     uint32_t dist,
     uint32_t *zPtr
 )
{
    uint32_t wordDist;
    uint_fast8_t innerDist;
    uint32_t *destPtr;
    uint_fast8_t i;

    wordDist = dist>>5;
    if ( wordDist < size_words ) {
	aPtr += indexMultiwordLoBut( size_words, wordDist );
	innerDist = dist & 31;
	if ( innerDist ) {
	    softfloat_shortShiftLeftM(
		size_words - wordDist,
		aPtr,
		innerDist,
		zPtr + indexMultiwordHiBut( size_words, wordDist )
	    );
	    if ( ! wordDist ) return;
	} else {
	    aPtr += indexWordHi( size_words - wordDist );
	    destPtr = zPtr + indexWordHi( size_words );
	    for ( i = size_words - wordDist; i; --i ) {
		*destPtr = *aPtr;
		aPtr -= wordIncr;
		destPtr -= wordIncr;
	    }
	}
	zPtr += indexMultiwordLo( size_words, wordDist );
    } else {
	wordDist = size_words;
    }
    do {
	*zPtr++ = 0;
	--wordDist;
    } while ( wordDist );
}

/*----------------------------------------------------------------------------
| This function or macro is the same as 'softfloat_shiftLeftM' with
| 'size_words' = 4 (N = 128).
*----------------------------------------------------------------------------*/
#define softfloat_shiftLeft128M( aPtr, dist, zPtr ) softfloat_shiftLeftM( 4, aPtr, dist, zPtr )

bool
softfloat_isNaNF128M( const uint32_t *aWPtr )
{
    uint32_t uiA96;

    uiA96 = aWPtr[indexWordHi( 4 )];
    if ( (~uiA96 & 0x7FFF0000) != 0 ) return 0;
    return
	((uiA96 & 0x0000FFFF) != 0)
	    || ((aWPtr[indexWord( 4, 2 )] | aWPtr[indexWord( 4, 1 )]
		     | aWPtr[indexWord( 4, 0 )])
		    != 0);
}

static inline uint_fast8_t
softfloat_countLeadingZeros32( uint32_t a ) {
    return a ? __builtin_clz( a ) : 32;
}

static inline bool
softfloat_eq128( uint64_t a64, uint64_t a0, uint64_t b64, uint64_t b0 ){
    return (a64 == b64) && (a0 == b0);
}

static inline bool
softfloat_lt128( uint64_t a64, uint64_t a0, uint64_t b64, uint64_t b0 ){
    return (a64 < b64) || ((a64 == b64) && (a0 < b0));
}

static inline bool
softfloat_le128( uint64_t a64, uint64_t a0, uint64_t b64, uint64_t b0 ){
    return (a64 < b64) || ((a64 == b64) && (a0 <= b0));
}

void
softfloat_propagateNaNF128M(
     const uint32_t *aWPtr, const uint32_t *bWPtr, uint32_t *zWPtr )
{
    bool isSigNaNA;
    const uint32_t *ptr;
    bool isSigNaNB;
    uint32_t uiA96, uiB96, wordMagA, wordMagB;

    isSigNaNA = f128M_isSignalingNaN( (const float128_t *) aWPtr );
    ptr = aWPtr;
    if ( ! bWPtr ) {
	if ( isSigNaNA ) softfloat_raiseFlags( softfloat_flag_invalid );
	goto copy;
    }
    isSigNaNB = f128M_isSignalingNaN( (const float128_t *) bWPtr );
    if ( isSigNaNA | isSigNaNB ) {
	softfloat_raiseFlags( softfloat_flag_invalid );
	if ( isSigNaNA ) {
	    if ( isSigNaNB ) goto returnLargerUIMag;
	    if ( softfloat_isNaNF128M( bWPtr ) ) goto copyB;
	    goto copy;
	} else {
	    if ( softfloat_isNaNF128M( aWPtr ) ) goto copy;
	    goto copyB;
	}
    }
 returnLargerUIMag:
    uiA96 = aWPtr[indexWordHi( 4 )];
    uiB96 = bWPtr[indexWordHi( 4 )];
    wordMagA = uiA96 & 0x7FFFFFFF;
    wordMagB = uiB96 & 0x7FFFFFFF;
    if ( wordMagA < wordMagB ) goto copyB;
    if ( wordMagB < wordMagA ) goto copy;
    wordMagA = aWPtr[indexWord( 4, 2 )];
    wordMagB = bWPtr[indexWord( 4, 2 )];
    if ( wordMagA < wordMagB ) goto copyB;
    if ( wordMagB < wordMagA ) goto copy;
    wordMagA = aWPtr[indexWord( 4, 1 )];
    wordMagB = bWPtr[indexWord( 4, 1 )];
    if ( wordMagA < wordMagB ) goto copyB;
    if ( wordMagB < wordMagA ) goto copy;
    wordMagA = aWPtr[indexWord( 4, 0 )];
    wordMagB = bWPtr[indexWord( 4, 0 )];
    if ( wordMagA < wordMagB ) goto copyB;
    if ( wordMagB < wordMagA ) goto copy;
    if ( uiA96 < uiB96 ) goto copy;
 copyB:
    ptr = bWPtr;
 copy:
    zWPtr[indexWordHi( 4 )] = ptr[indexWordHi( 4 )] | 0x00008000;
    zWPtr[indexWord( 4, 2 )] = ptr[indexWord( 4, 2 )];
    zWPtr[indexWord( 4, 1 )] = ptr[indexWord( 4, 1 )];
    zWPtr[indexWord( 4, 0 )] = ptr[indexWord( 4, 0 )];
}

static inline uint_fast8_t
softfloat_countLeadingZeros64( uint64_t a )
{
    uint_fast8_t count;
    uint32_t a32;

    count = 0;
    a32 = a>>32;
    if ( ! a32 ) {
	count = 32;
	a32 = a;
    }
    /*------------------------------------------------------------------------
    | From here, result is current count + count leading zeros of `a32'.
    *------------------------------------------------------------------------*/
    if ( a32 < 0x10000 ) {
	count += 16;
	a32 <<= 16;
    }
    if ( a32 < 0x1000000 ) {
	count += 8;
	a32 <<= 8;
    }
    count += softfloat_countLeadingZeros8[a32>>24];
    return count;
}

static inline struct exp16_sig64
softfloat_normSubnormalF64Sig( uint_fast64_t sig )
{
    int_fast8_t shiftDist;
    struct exp16_sig64 z;

    shiftDist = softfloat_countLeadingZeros64( sig ) - 11;
    z.exp = 1 - shiftDist;
    z.sig = sig<<shiftDist;
    return z;
}

static inline struct uint128_extra
softfloat_shiftRightJam128Extra(
     uint64_t a64, uint64_t a0, uint64_t extra, uint_fast32_t dist )
{
    uint_fast8_t u8NegDist;
    struct uint128_extra z;

    u8NegDist = -dist;
    if ( dist < 64 ) {
	z.v.v64 = a64>>dist;
	z.v.v0 = a64<<(u8NegDist & 63) | a0>>dist;
	z.extra = a0<<(u8NegDist & 63);
    } else {
	z.v.v64 = 0;
	if ( dist == 64 ) {
	    z.v.v0 = a64;
	    z.extra = a0;
	} else {
	    extra |= a0;
	    if ( dist < 128 ) {
		z.v.v0 = a64>>(dist & 63);
		z.extra = a64<<(u8NegDist & 63);
	    } else {
		z.v.v0 = 0;
		z.extra = (dist == 128) ? a64 : (a64 != 0);
	    }
	}
    }
    z.extra |= (extra != 0);
    return z;
}

static inline struct uint128_extra
softfloat_shortShiftRightJam128Extra(
     uint64_t a64, uint64_t a0, uint64_t extra, uint_fast8_t dist )
{
    uint_fast8_t negDist = -dist;
    struct uint128_extra z;
    z.v.v64 = a64>>dist;
    z.v.v0 = a64<<(negDist & 63) | a0>>dist;
    z.extra = a0<<(negDist & 63) | (extra != 0);
    return z;
}

static inline struct uint128
softfloat_shiftRightJam128( uint64_t a64, uint64_t a0, uint_fast32_t dist )
{
    uint_fast8_t u8NegDist;
    struct uint128 z;

    if ( dist < 64 ) {
	u8NegDist = -dist;
	z.v64 = a64>>dist;
	z.v0 =
	    a64<<(u8NegDist & 63) | a0>>dist
		| ((uint64_t) (a0<<(u8NegDist & 63)) != 0);
    } else {
	z.v64 = 0;
	z.v0 =
	    (dist < 127)
		? a64>>(dist & 63)
		      | (((a64 & (((uint_fast64_t) 1<<(dist & 63)) - 1)) | a0)
			     != 0)
		: ((a64 | a0) != 0);
    }
    return z;
}

static inline struct uint128
softfloat_shortShiftLeft128( uint64_t a64, uint64_t a0, uint_fast8_t dist )
{
    struct uint128 z;
    z.v64 = a64<<dist | a0>>(-dist & 63);
    z.v0 = a0<<dist;
    return z;
}

static int_fast8_t
softfloat_compare128M( const uint32_t *aPtr, const uint32_t *bPtr )
{
    unsigned int index, lastIndex;
    uint32_t wordA, wordB;

    index = indexWordHi( 4 );
    lastIndex = indexWordLo( 4 );
    for (;;) {
	wordA = aPtr[index];
	wordB = bPtr[index];
	if ( wordA != wordB ) return (wordA < wordB) ? -1 : 1;
	if ( index == lastIndex ) break;
	index -= wordIncr;
    }
    return 0;
}

static void
softfloat_roundPackMToF128M(
     bool sign, int32_t exp, uint32_t *extSigPtr, uint32_t *zWPtr )
{
    uint_fast8_t roundingMode;
    bool roundNearEven;
    uint32_t sigExtra;
    bool doIncrement, isTiny;
    static const uint32_t maxSig[4] =
	INIT_UINTM4( 0x0001FFFF, 0xFFFFFFFF, 0xFFFFFFFF, 0xFFFFFFFF );
    uint32_t ui, uj;

    /*------------------------------------------------------------------------
    *------------------------------------------------------------------------*/
    roundingMode = softfloat_roundingMode;
    roundNearEven = (roundingMode == softfloat_round_near_even);
    sigExtra = extSigPtr[indexWordLo( 5 )];
    doIncrement = (0x80000000 <= sigExtra);
    if ( ! roundNearEven && (roundingMode != softfloat_round_near_maxMag) ) {
	doIncrement =
	    (roundingMode
		 == (sign ? softfloat_round_min : softfloat_round_max))
		&& sigExtra;
    }
    /*------------------------------------------------------------------------
    *------------------------------------------------------------------------*/
    if ( 0x7FFD <= (uint32_t) exp ) {
	if ( exp < 0 ) {
	    /*----------------------------------------------------------------
	    *----------------------------------------------------------------*/
	    isTiny =
		   (softfloat_detectTininess
			== softfloat_tininess_beforeRounding)
		|| (exp < -1)
		|| ! doIncrement
		|| (softfloat_compare128M(
			extSigPtr + indexMultiwordHi( 5, 4 ), maxSig )
			< 0);
	    softfloat_shiftRightJam160M( extSigPtr, -exp, extSigPtr );
	    exp = 0;
	    sigExtra = extSigPtr[indexWordLo( 5 )];
	    if ( isTiny && sigExtra ) {
		softfloat_raiseFlags( softfloat_flag_underflow );
	    }
	    doIncrement = (0x80000000 <= sigExtra);
	    if (
		   ! roundNearEven
		&& (roundingMode != softfloat_round_near_maxMag)
	    ) {
		doIncrement =
		    (roundingMode
			 == (sign ? softfloat_round_min : softfloat_round_max))
			&& sigExtra;
	    }
	} else if (
	       (0x7FFD < exp)
	    || ((exp == 0x7FFD) && doIncrement
		    && (softfloat_compare128M(
			    extSigPtr + indexMultiwordHi( 5, 4 ), maxSig )
			    == 0))
	) {
	    /*----------------------------------------------------------------
	    *----------------------------------------------------------------*/
	    softfloat_raiseFlags(
		softfloat_flag_overflow | softfloat_flag_inexact );
	    if (
		   roundNearEven
		|| (roundingMode == softfloat_round_near_maxMag)
		|| (roundingMode
			== (sign ? softfloat_round_min : softfloat_round_max))
	    ) {
		ui = packToF128UI96( sign, 0x7FFF, 0 );
		uj = 0;
	    } else {
		ui = packToF128UI96( sign, 0x7FFE, 0x0000FFFF );
		uj = 0xFFFFFFFF;
	    }
	    zWPtr[indexWordHi( 4 )] = ui;
	    zWPtr[indexWord( 4, 2 )] = uj;
	    zWPtr[indexWord( 4, 1 )] = uj;
	    zWPtr[indexWord( 4, 0 )] = uj;
	    return;
	}
    }
    /*------------------------------------------------------------------------
    *------------------------------------------------------------------------*/
    uj = extSigPtr[indexWord( 5, 1 )];
    if ( sigExtra ) {
	softfloat_exceptionFlags |= softfloat_flag_inexact;
#ifdef SOFTFLOAT_ROUND_ODD
	if ( roundingMode == softfloat_round_odd ) {
	    uj |= 1;
	    goto noIncrementPackReturn;
	}
#endif
    }
    if ( doIncrement ) {
	++uj;
	if ( uj ) {
	    if ( ! (sigExtra & 0x7FFFFFFF) && roundNearEven ) uj &= ~1;
	    zWPtr[indexWord( 4, 2 )] = extSigPtr[indexWord( 5, 3 )];
	    zWPtr[indexWord( 4, 1 )] = extSigPtr[indexWord( 5, 2 )];
	    zWPtr[indexWord( 4, 0 )] = uj;
	    ui = extSigPtr[indexWordHi( 5 )];
	} else {
	    zWPtr[indexWord( 4, 0 )] = uj;
	    ui = extSigPtr[indexWord( 5, 2 )] + 1;
	    zWPtr[indexWord( 4, 1 )] = ui;
	    uj = extSigPtr[indexWord( 5, 3 )];
	    if ( ui ) {
		zWPtr[indexWord( 4, 2 )] = uj;
		ui = extSigPtr[indexWordHi( 5 )];
	    } else {
		++uj;
		zWPtr[indexWord( 4, 2 )] = uj;
		ui = extSigPtr[indexWordHi( 5 )];
		if ( ! uj ) ++ui;
	    }
	}
    } else {
 noIncrementPackReturn:
	zWPtr[indexWord( 4, 0 )] = uj;
	ui = extSigPtr[indexWord( 5, 2 )];
	zWPtr[indexWord( 4, 1 )] = ui;
	uj |= ui;
	ui = extSigPtr[indexWord( 5, 3 )];
	zWPtr[indexWord( 4, 2 )] = ui;
	uj |= ui;
	ui = extSigPtr[indexWordHi( 5 )];
	uj |= ui;
	if ( ! uj ) exp = 0;
    }
    zWPtr[indexWordHi( 4 )] = packToF128UI96( sign, exp, ui );
}

static int
softfloat_shiftNormSigF128M(
     const uint32_t *wPtr, uint_fast8_t shiftDist, uint32_t *sigPtr )
{
    uint32_t wordSig;
    int32_t exp;
    uint32_t leadingBit;

    wordSig = wPtr[indexWordHi( 4 )];
    exp = expF128UI96( wordSig );
    if ( exp ) {
	softfloat_shortShiftLeft128M( wPtr, shiftDist, sigPtr );
	leadingBit = 0x00010000<<shiftDist;
	sigPtr[indexWordHi( 4 )] =
	    (sigPtr[indexWordHi( 4 )] & (leadingBit - 1)) | leadingBit;
    } else {
	exp = 16;
	wordSig &= 0x7FFFFFFF;
	if ( ! wordSig ) {
	    exp = -16;
	    wordSig = wPtr[indexWord( 4, 2 )];
	    if ( ! wordSig ) {
		exp = -48;
		wordSig = wPtr[indexWord( 4, 1 )];
		if ( ! wordSig ) {
		    wordSig = wPtr[indexWord( 4, 0 )];
		    if ( ! wordSig ) return -128;
		    exp = -80;
		}
	    }
	}
	exp -= softfloat_countLeadingZeros32( wordSig );
	softfloat_shiftLeft128M( wPtr, 1 - exp + shiftDist, sigPtr );
    }
    return exp;
}

bool
softfloat_tryPropagateNaNF128M(
     const uint32_t *aWPtr, const uint32_t *bWPtr, uint32_t *zWPtr )
{

    if ( softfloat_isNaNF128M( aWPtr ) || softfloat_isNaNF128M( bWPtr ) ) {
	softfloat_propagateNaNF128M( aWPtr, bWPtr, zWPtr );
	return 1;
    }
    return 0;
}

static inline struct uint128
softfloat_add128( uint64_t a64, uint64_t a0, uint64_t b64, uint64_t b0 )
{
    struct uint128 z;
    z.v0 = a0 + b0;
    z.v64 = a64 + b64 + (z.v0 < a0);
    return z;
}

static inline struct uint128
softfloat_sub128( uint64_t a64, uint64_t a0, uint64_t b64, uint64_t b0 )
{
    struct uint128 z;
    z.v0 = a0 - b0;
    z.v64 = a64 - b64;
    z.v64 -= (a0 < b0);
    return z;
}

static inline struct uint128
softfloat_mul128By32( uint64_t a64, uint64_t a0, uint32_t b )
{
    union { unsigned __int128 ui; struct uint128 s; } uZ;
    uZ.ui = ((unsigned __int128) a64<<64 | a0) * b;
    return uZ.s;
}

static void
softfloat_mul128MTo256M(
     const uint32_t *aPtr, const uint32_t *bPtr, uint32_t *zPtr )
{
    uint32_t *lastZPtr, wordB;
    uint64_t dwordProd;
    uint32_t wordZ;
    uint_fast8_t carry;

    bPtr += indexWordLo( 4 );
    lastZPtr = zPtr + indexMultiwordHi( 8, 5 );
    zPtr += indexMultiwordLo( 8, 5 );
    wordB = *bPtr;
    dwordProd = (uint64_t) aPtr[indexWord( 4, 0 )] * wordB;
    zPtr[indexWord( 5, 0 )] = dwordProd;
    dwordProd = (uint64_t) aPtr[indexWord( 4, 1 )] * wordB + (dwordProd>>32);
    zPtr[indexWord( 5, 1 )] = dwordProd;
    dwordProd = (uint64_t) aPtr[indexWord( 4, 2 )] * wordB + (dwordProd>>32);
    zPtr[indexWord( 5, 2 )] = dwordProd;
    dwordProd = (uint64_t) aPtr[indexWord( 4, 3 )] * wordB + (dwordProd>>32);
    zPtr[indexWord( 5, 3 )] = dwordProd;
    zPtr[indexWord( 5, 4 )] = dwordProd>>32;
    do {
	bPtr += wordIncr;
	zPtr += wordIncr;
	wordB = *bPtr;
	dwordProd = (uint64_t) aPtr[indexWord( 4, 0 )] * wordB;
	wordZ = zPtr[indexWord( 5, 0 )] + (uint32_t) dwordProd;
	zPtr[indexWord( 5, 0 )] = wordZ;
	carry = (wordZ < (uint32_t) dwordProd);
	dwordProd =
	    (uint64_t) aPtr[indexWord( 4, 1 )] * wordB + (dwordProd>>32);
	wordZ = zPtr[indexWord( 5, 1 )] + (uint32_t) dwordProd + carry;
	zPtr[indexWord( 5, 1 )] = wordZ;
	if ( wordZ != (uint32_t) dwordProd ) {
	    carry = (wordZ < (uint32_t) dwordProd);
	}
	dwordProd =
	    (uint64_t) aPtr[indexWord( 4, 2 )] * wordB + (dwordProd>>32);
	wordZ = zPtr[indexWord( 5, 2 )] + (uint32_t) dwordProd + carry;
	zPtr[indexWord( 5, 2 )] = wordZ;
	if ( wordZ != (uint32_t) dwordProd ) {
	    carry = (wordZ < (uint32_t) dwordProd);
	}
	dwordProd =
	    (uint64_t) aPtr[indexWord( 4, 3 )] * wordB + (dwordProd>>32);
	wordZ = zPtr[indexWord( 5, 3 )] + (uint32_t) dwordProd + carry;
	zPtr[indexWord( 5, 3 )] = wordZ;
	if ( wordZ != (uint32_t) dwordProd ) {
	    carry = (wordZ < (uint32_t) dwordProd);
	}
	zPtr[indexWord( 5, 4 )] = (dwordProd>>32) + carry;
    } while ( zPtr != lastZPtr );
}

float128_t
softfloat_roundPackToF128(
     bool sign,
     int_fast32_t exp,
     uint_fast64_t sig64,
     uint_fast64_t sig0,
     uint_fast64_t sigExtra
 )
{
    uint_fast8_t roundingMode;
    bool roundNearEven, doIncrement, isTiny;
    struct uint128_extra sig128Extra;
    uint_fast64_t uiZ64, uiZ0;
    struct uint128 sig128;
    union ui128_f128 uZ;

    /*------------------------------------------------------------------------
    *------------------------------------------------------------------------*/
    roundingMode = softfloat_roundingMode;
    roundNearEven = (roundingMode == softfloat_round_near_even);
    doIncrement = (UINT64_C( 0x8000000000000000 ) <= sigExtra);
    if ( ! roundNearEven && (roundingMode != softfloat_round_near_maxMag) ) {
	doIncrement =
	    (roundingMode
		 == (sign ? softfloat_round_min : softfloat_round_max))
		&& sigExtra;
    }
    /*------------------------------------------------------------------------
    *------------------------------------------------------------------------*/
    if ( 0x7FFD <= (uint32_t) exp ) {
	if ( exp < 0 ) {
	    /*----------------------------------------------------------------
	    *----------------------------------------------------------------*/
	    isTiny =
		   (softfloat_detectTininess
			== softfloat_tininess_beforeRounding)
		|| (exp < -1)
		|| ! doIncrement
		|| softfloat_lt128(
		       sig64,
		       sig0,
		       UINT64_C( 0x0001FFFFFFFFFFFF ),
		       UINT64_C( 0xFFFFFFFFFFFFFFFF )
		   );
	    sig128Extra =
		softfloat_shiftRightJam128Extra( sig64, sig0, sigExtra, -exp );
	    sig64 = sig128Extra.v.v64;
	    sig0  = sig128Extra.v.v0;
	    sigExtra = sig128Extra.extra;
	    exp = 0;
	    if ( isTiny && sigExtra ) {
		softfloat_raiseFlags( softfloat_flag_underflow );
	    }
	    doIncrement = (UINT64_C( 0x8000000000000000 ) <= sigExtra);
	    if (
		   ! roundNearEven
		&& (roundingMode != softfloat_round_near_maxMag)
	    ) {
		doIncrement =
		    (roundingMode
			 == (sign ? softfloat_round_min : softfloat_round_max))
			&& sigExtra;
	    }
	} else if (
	       (0x7FFD < exp)
	    || ((exp == 0x7FFD)
		    && softfloat_eq128(
			   sig64,
			   sig0,
			   UINT64_C( 0x0001FFFFFFFFFFFF ),
			   UINT64_C( 0xFFFFFFFFFFFFFFFF )
		       )
		    && doIncrement)
	) {
	    /*----------------------------------------------------------------
	    *----------------------------------------------------------------*/
	    softfloat_raiseFlags(
		softfloat_flag_overflow | softfloat_flag_inexact );
	    if (
		   roundNearEven
		|| (roundingMode == softfloat_round_near_maxMag)
		|| (roundingMode
			== (sign ? softfloat_round_min : softfloat_round_max))
	    ) {
		uiZ64 = packToF128UI64( sign, 0x7FFF, 0 );
		uiZ0  = 0;
	    } else {
		uiZ64 = packToF128UI64( sign, 0x7FFE, UINT64_C( 0x0000FFFFFFFFFFFF ) );
		uiZ0 = UINT64_C( 0xFFFFFFFFFFFFFFFF );
	    }
	    goto uiZ;
	}
    }
    /*------------------------------------------------------------------------
    *------------------------------------------------------------------------*/
    if ( sigExtra ) {
	softfloat_exceptionFlags |= softfloat_flag_inexact;
#ifdef SOFTFLOAT_ROUND_ODD
	if ( roundingMode == softfloat_round_odd ) {
	    sig0 |= 1;
	    goto packReturn;
	}
#endif
    }
    if ( doIncrement ) {
	sig128 = softfloat_add128( sig64, sig0, 0, 1 );
	sig64 = sig128.v64;
	sig0 =
	    sig128.v0
		& ~(uint64_t)
		       (! (sigExtra & UINT64_C( 0x7FFFFFFFFFFFFFFF ))
			    & roundNearEven);
    } else {
	if ( ! (sig64 | sig0) ) exp = 0;
    }
    /*------------------------------------------------------------------------
    *------------------------------------------------------------------------*/
 packReturn:
    uiZ64 = packToF128UI64( sign, exp, sig64 );
    uiZ0  = sig0;
 uiZ:
    uZ.ui.v64 = uiZ64;
    uZ.ui.v0  = uiZ0;
    return uZ.f;
}

float128_t
softfloat_normRoundPackToF128(
     bool sign, int_fast32_t exp, uint_fast64_t sig64, uint_fast64_t sig0 )
{
    int_fast8_t shiftDist;
    struct uint128 sig128;
    union ui128_f128 uZ;
    uint_fast64_t sigExtra;
    struct uint128_extra sig128Extra;

    if ( ! sig64 ) {
	exp -= 64;
	sig64 = sig0;
	sig0 = 0;
    }
    shiftDist = softfloat_countLeadingZeros64( sig64 ) - 15;
    exp -= shiftDist;
    if ( 0 <= shiftDist ) {
	if ( shiftDist ) {
	    sig128 = softfloat_shortShiftLeft128( sig64, sig0, shiftDist );
	    sig64 = sig128.v64;
	    sig0  = sig128.v0;
	}
	if ( (uint32_t) exp < 0x7FFD ) {
	    uZ.ui.v64 = packToF128UI64( sign, sig64 | sig0 ? exp : 0, sig64 );
	    uZ.ui.v0  = sig0;
	    return uZ.f;
	}
	sigExtra = 0;
    } else {
	sig128Extra =
	    softfloat_shortShiftRightJam128Extra( sig64, sig0, 0, -shiftDist );
	sig64 = sig128Extra.v.v64;
	sig0  = sig128Extra.v.v0;
	sigExtra = sig128Extra.extra;
    }
    return softfloat_roundPackToF128( sign, exp, sig64, sig0, sigExtra );
}

struct exp32_sig128
softfloat_normSubnormalF128Sig( uint_fast64_t sig64, uint_fast64_t sig0 )
{
    int_fast8_t shiftDist;
    struct exp32_sig128 z;

    if ( ! sig64 ) {
	shiftDist = softfloat_countLeadingZeros64( sig0 ) - 15;
	z.exp = -63 - shiftDist;
	if ( shiftDist < 0 ) {
	    z.sig.v64 = sig0>>-shiftDist;
	    z.sig.v0  = sig0<<(shiftDist & 63);
	} else {
	    z.sig.v64 = sig0<<shiftDist;
	    z.sig.v0  = 0;
	}
    } else {
	shiftDist = softfloat_countLeadingZeros64( sig64 ) - 15;
	z.exp = 1 - shiftDist;
	z.sig = softfloat_shortShiftLeft128( sig64, sig0, shiftDist );
    }
    return z;
}

struct uint128
softfloat_propagateNaNF128UI(
     uint_fast64_t uiA64,
     uint_fast64_t uiA0,
     uint_fast64_t uiB64,
     uint_fast64_t uiB0
 )
{
    bool isSigNaNA, isSigNaNB;
    uint_fast64_t uiNonsigA64, uiNonsigB64, uiMagA64, uiMagB64;
    struct uint128 uiZ;

    /*------------------------------------------------------------------------
    *------------------------------------------------------------------------*/
    isSigNaNA = softfloat_isSigNaNF128UI( uiA64, uiA0 );
    isSigNaNB = softfloat_isSigNaNF128UI( uiB64, uiB0 );
    /*------------------------------------------------------------------------
    | Make NaNs non-signaling.
    *------------------------------------------------------------------------*/
    uiNonsigA64 = uiA64 | UINT64_C( 0x0000800000000000 );
    uiNonsigB64 = uiB64 | UINT64_C( 0x0000800000000000 );
    /*------------------------------------------------------------------------
    *------------------------------------------------------------------------*/
    if ( isSigNaNA | isSigNaNB ) {
	softfloat_raiseFlags( softfloat_flag_invalid );
	if ( isSigNaNA ) {
	    if ( isSigNaNB ) goto returnLargerMag;
	    if ( isNaNF128UI( uiB64, uiB0 ) ) goto returnB;
	    goto returnA;
	} else {
	    if ( isNaNF128UI( uiA64, uiA0 ) ) goto returnA;
	    goto returnB;
	}
    }
 returnLargerMag:
    uiMagA64 = uiA64 & UINT64_C( 0x7FFFFFFFFFFFFFFF );
    uiMagB64 = uiB64 & UINT64_C( 0x7FFFFFFFFFFFFFFF );
    if ( uiMagA64 < uiMagB64 ) goto returnB;
    if ( uiMagB64 < uiMagA64 ) goto returnA;
    if ( uiA0 < uiB0 ) goto returnB;
    if ( uiB0 < uiA0 ) goto returnA;
    if ( uiNonsigA64 < uiNonsigB64 ) goto returnA;
 returnB:
    uiZ.v64 = uiNonsigB64;
    uiZ.v0  = uiB0;
    return uiZ;
 returnA:
    uiZ.v64 = uiNonsigA64;
    uiZ.v0  = uiA0;
    return uiZ;

}

#if __POINTER_SIZE__ == 4

void
i32_to_f128M( int32_t a, float128_t *zPtr )
{
    uint32_t *zWPtr;
    uint32_t uiZ96, uiZ64;
    bool sign;
    uint32_t absA;
    int_fast8_t shiftDist;
    uint64_t normAbsA;

    zWPtr = (uint32_t *) zPtr;
    uiZ96 = 0;
    uiZ64 = 0;
    if ( a ) {
	sign = (a < 0);
	absA = sign ? -(uint32_t) a : (uint32_t) a;
	shiftDist = softfloat_countLeadingZeros32( absA ) + 17;
	normAbsA = (uint64_t) absA<<shiftDist;
	uiZ96 = packToF128UI96( sign, 0x402E - shiftDist, normAbsA>>32 );
	uiZ64 = normAbsA;
    }
    zWPtr[indexWord( 4, 3 )] = uiZ96;
    zWPtr[indexWord( 4, 2 )] = uiZ64;
    zWPtr[indexWord( 4, 1 )] = 0;
    zWPtr[indexWord( 4, 0 )] = 0;
}

#else

void
i64_to_f128M( int64_t a, float128_t *zPtr )
{
    uint32_t *zWPtr;
    uint32_t uiZ96, uiZ64;
    bool sign;
    uint64_t absA;
    uint_fast8_t shiftDist;
    uint32_t *ptr;

    zWPtr = (uint32_t *) zPtr;
    uiZ96 = 0;
    uiZ64 = 0;
    zWPtr[indexWord( 4, 1 )] = 0;
    zWPtr[indexWord( 4, 0 )] = 0;
    if ( a ) {
	sign = (a < 0);
	absA = sign ? -(uint64_t) a : (uint64_t) a;
	shiftDist = softfloat_countLeadingZeros64( absA ) + 17;
	if ( shiftDist < 32 ) {
	    ptr = zWPtr + indexMultiwordHi( 4, 3 );
	    ptr[indexWord( 3, 2 )] = 0;
	    ptr[indexWord( 3, 1 )] = absA>>32;
	    ptr[indexWord( 3, 0 )] = absA;
	    softfloat_shortShiftLeft96M( ptr, shiftDist, ptr );
	    ptr[indexWordHi( 3 )] =
		packToF128UI96(
		    sign, 0x404E - shiftDist, ptr[indexWordHi( 3 )] );
	    return;
	}
	absA <<= shiftDist - 32;
	uiZ96 = packToF128UI96( sign, 0x404E - shiftDist, absA>>32 );
	uiZ64 = absA;
    }
    zWPtr[indexWord( 4, 3 )] = uiZ96;
    zWPtr[indexWord( 4, 2 )] = uiZ64;
}

#endif

void
f64_to_f128M( float64_t a, float128_t *zPtr )
{
    uint32_t *zWPtr;
    union ui64_f64 uA;
    uint64_t uiA;
    bool sign;
    int_fast16_t exp;
    uint64_t frac;
    struct commonNaN commonNaN;
    uint32_t uiZ96;
    struct exp16_sig64 normExpSig;

    /*------------------------------------------------------------------------
    *------------------------------------------------------------------------*/
    zWPtr = (uint32_t *) zPtr;
    /*------------------------------------------------------------------------
    *------------------------------------------------------------------------*/
    uA.f = a;
    uiA = uA.ui;
    sign = signF64UI( uiA );
    exp  = expF64UI( uiA );
    frac = fracF64UI( uiA );
    /*------------------------------------------------------------------------
    *------------------------------------------------------------------------*/
    zWPtr[indexWord( 4, 0 )] = 0;
    if ( exp == 0x7FF ) {
	if ( frac ) {
	    // NaN
	    softfloat_commonNaNToF128M( zWPtr );
	    return;
	}
	uiZ96 = packToF128UI96( sign, 0x7FFF, 0 );
	goto uiZ;
    }
    /*------------------------------------------------------------------------
    *------------------------------------------------------------------------*/
    if ( ! exp ) {
	if ( ! frac ) {
	    uiZ96 = packToF128UI96( sign, 0, 0 );
	    goto uiZ;
	}
	normExpSig = softfloat_normSubnormalF64Sig( frac );
	exp = normExpSig.exp - 1;
	frac = normExpSig.sig;
    }
    /*------------------------------------------------------------------------
    *------------------------------------------------------------------------*/
    zWPtr[indexWord( 4, 1 )] = (uint32_t) frac<<28;
    frac >>= 4;
    zWPtr[indexWordHi( 4 )] = packToF128UI96( sign, exp + 0x3C00, frac>>32 );
    zWPtr[indexWord( 4, 2 )] = frac;
    return;
    /*------------------------------------------------------------------------
    *------------------------------------------------------------------------*/
 uiZ:
    zWPtr[indexWord( 4, 3 )] = uiZ96;
    zWPtr[indexWord( 4, 2 )] = 0;
    zWPtr[indexWord( 4, 1 )] = 0;

}

float128_t
extF80_to_f128( extFloat80_t a )
{
    union { struct extFloat80M s; extFloat80_t f; } uA;
    uint_fast16_t uiA64;
    uint_fast64_t uiA0;
    uint_fast16_t exp;
    uint_fast64_t frac;
    struct commonNaN commonNaN;
    struct uint128 uiZ;
    bool sign;
    struct uint128 frac128;
    union ui128_f128 uZ;

    uA.f = a;
    uiA64 = uA.s.signExp;
    uiA0  = uA.s.signif;
    exp = expExtF80UI64( uiA64 );
    frac = uiA0 & UINT64_C( 0x7FFFFFFFFFFFFFFF );
    if ( (exp == 0x7FFF) && frac ) {
	softfloat_extF80UIToCommonNaN( uiA64, uiA0, &commonNaN );
	softfloat_commonNaNToF128M( (uint32_t*) &uiZ );
    } else {
	sign = signExtF80UI64( uiA64 );
	frac128 = softfloat_shortShiftLeft128( 0, frac, 49 );
	uiZ.v64 = packToF128UI64( sign, exp, frac128.v64 );
	uiZ.v0  = frac128.v0;
    }
    uZ.ui = uiZ;
    return uZ.f;

}

static float128_t
softfloat_addMagsF128(
     uint64_t uiA64,
     uint64_t uiA0,
     uint64_t uiB64,
     uint64_t uiB0,
     bool signZ
 )
{
    int32_t expA;
    struct uint128 sigA;
    int32_t expB;
    struct uint128 sigB;
    int32_t expDiff;
    struct uint128 uiZ, sigZ;
    int32_t expZ;
    uint64_t sigZExtra;
    struct uint128_extra sig128Extra;
    union ui128_f128 uZ;

    expA = expF128UI64( uiA64 );
    sigA.v64 = fracF128UI64( uiA64 );
    sigA.v0  = uiA0;
    expB = expF128UI64( uiB64 );
    sigB.v64 = fracF128UI64( uiB64 );
    sigB.v0  = uiB0;
    expDiff = expA - expB;
    if ( ! expDiff ) {
	if ( expA == 0x7FFF ) {
	    if ( sigA.v64 | sigA.v0 | sigB.v64 | sigB.v0 ) goto propagateNaN;
	    uiZ.v64 = uiA64;
	    uiZ.v0  = uiA0;
	    goto uiZ;
	}
	sigZ = softfloat_add128( sigA.v64, sigA.v0, sigB.v64, sigB.v0 );
	if ( ! expA ) {
	    uiZ.v64 = packToF128UI64( signZ, 0, sigZ.v64 );
	    uiZ.v0  = sigZ.v0;
	    goto uiZ;
	}
	expZ = expA;
	sigZ.v64 |= UINT64_C( 0x0002000000000000 );
	sigZExtra = 0;
	goto shiftRight1;
    }
    if ( expDiff < 0 ) {
	if ( expB == 0x7FFF ) {
	    if ( sigB.v64 | sigB.v0 ) goto propagateNaN;
	    uiZ.v64 = packToF128UI64( signZ, 0x7FFF, 0 );
	    uiZ.v0  = 0;
	    goto uiZ;
	}
	expZ = expB;
	if ( expA ) {
	    sigA.v64 |= UINT64_C( 0x0001000000000000 );
	} else {
	    ++expDiff;
	    sigZExtra = 0;
	    if ( ! expDiff ) goto newlyAligned;
	}
	sig128Extra =
	    softfloat_shiftRightJam128Extra( sigA.v64, sigA.v0, 0, -expDiff );
	sigA = sig128Extra.v;
	sigZExtra = sig128Extra.extra;
    } else {
	if ( expA == 0x7FFF ) {
	    if ( sigA.v64 | sigA.v0 ) goto propagateNaN;
	    uiZ.v64 = uiA64;
	    uiZ.v0  = uiA0;
	    goto uiZ;
	}
	expZ = expA;
	if ( expB ) {
	    sigB.v64 |= UINT64_C( 0x0001000000000000 );
	} else {
	    --expDiff;
	    sigZExtra = 0;
	    if ( ! expDiff ) goto newlyAligned;
	}
	sig128Extra =
	    softfloat_shiftRightJam128Extra( sigB.v64, sigB.v0, 0, expDiff );
	sigB = sig128Extra.v;
	sigZExtra = sig128Extra.extra;
    }
 newlyAligned:
    sigZ =
	softfloat_add128(
	    sigA.v64 | UINT64_C( 0x0001000000000000 ),
	    sigA.v0,
	    sigB.v64,
	    sigB.v0
	);
    --expZ;
    if ( sigZ.v64 < UINT64_C( 0x0002000000000000 ) ) goto roundAndPack;
    ++expZ;
 shiftRight1:
    sig128Extra =
	softfloat_shortShiftRightJam128Extra(
	    sigZ.v64, sigZ.v0, sigZExtra, 1 );
    sigZ = sig128Extra.v;
    sigZExtra = sig128Extra.extra;
 roundAndPack:
    return
	softfloat_roundPackToF128( signZ, expZ, sigZ.v64, sigZ.v0, sigZExtra );
 propagateNaN:
    uiZ = softfloat_propagateNaNF128UI( uiA64, uiA0, uiB64, uiB0 );
 uiZ:
    uZ.ui = uiZ;
    return uZ.f;

}

float128_t
softfloat_subMagsF128(
     uint_fast64_t uiA64,
     uint_fast64_t uiA0,
     uint_fast64_t uiB64,
     uint_fast64_t uiB0,
     bool signZ
 )
{
    int_fast32_t expA;
    struct uint128 sigA;
    int_fast32_t expB;
    struct uint128 sigB, sigZ;
    int_fast32_t expDiff, expZ;
    struct uint128 uiZ;
    union ui128_f128 uZ;

    expA = expF128UI64( uiA64 );
    sigA.v64 = fracF128UI64( uiA64 );
    sigA.v0  = uiA0;
    expB = expF128UI64( uiB64 );
    sigB.v64 = fracF128UI64( uiB64 );
    sigB.v0  = uiB0;
    sigA = softfloat_shortShiftLeft128( sigA.v64, sigA.v0, 4 );
    sigB = softfloat_shortShiftLeft128( sigB.v64, sigB.v0, 4 );
    expDiff = expA - expB;
    if ( 0 < expDiff ) goto expABigger;
    if ( expDiff < 0 ) goto expBBigger;
    if ( expA == 0x7FFF ) {
	if ( sigA.v64 | sigA.v0 | sigB.v64 | sigB.v0 ) goto propagateNaN;
	softfloat_raiseFlags( softfloat_flag_invalid );
	softfloat_commonNaNToF128M( (uint32_t*)&uiZ );
	// uiZ.v64 = defaultNaNF128UI64;
	// uiZ.v0  = defaultNaNF128UI0;
	goto uiZ;
    }
    expZ = expA;
    if ( ! expZ ) expZ = 1;
    if ( sigB.v64 < sigA.v64 ) goto aBigger;
    if ( sigA.v64 < sigB.v64 ) goto bBigger;
    if ( sigB.v0 < sigA.v0 ) goto aBigger;
    if ( sigA.v0 < sigB.v0 ) goto bBigger;
    uiZ.v64 =
	packToF128UI64(
	    (softfloat_roundingMode == softfloat_round_min), 0, 0 );
    uiZ.v0 = 0;
    goto uiZ;
 expBBigger:
    if ( expB == 0x7FFF ) {
	if ( sigB.v64 | sigB.v0 ) goto propagateNaN;
	uiZ.v64 = packToF128UI64( signZ ^ 1, 0x7FFF, 0 );
	uiZ.v0  = 0;
	goto uiZ;
    }
    if ( expA ) {
	sigA.v64 |= UINT64_C( 0x0010000000000000 );
    } else {
	++expDiff;
	if ( ! expDiff ) goto newlyAlignedBBigger;
    }
    sigA = softfloat_shiftRightJam128( sigA.v64, sigA.v0, -expDiff );
 newlyAlignedBBigger:
    expZ = expB;
    sigB.v64 |= UINT64_C( 0x0010000000000000 );
 bBigger:
    signZ = ! signZ;
    sigZ = softfloat_sub128( sigB.v64, sigB.v0, sigA.v64, sigA.v0 );
    goto normRoundPack;
 expABigger:
    if ( expA == 0x7FFF ) {
	if ( sigA.v64 | sigA.v0 ) goto propagateNaN;
	uiZ.v64 = uiA64;
	uiZ.v0  = uiA0;
	goto uiZ;
    }
    if ( expB ) {
	sigB.v64 |= UINT64_C( 0x0010000000000000 );
    } else {
	--expDiff;
	if ( ! expDiff ) goto newlyAlignedABigger;
    }
    sigB = softfloat_shiftRightJam128( sigB.v64, sigB.v0, expDiff );
 newlyAlignedABigger:
    expZ = expA;
    sigA.v64 |= UINT64_C( 0x0010000000000000 );
 aBigger:
    sigZ = softfloat_sub128( sigA.v64, sigA.v0, sigB.v64, sigB.v0 );
 normRoundPack:
    return softfloat_normRoundPackToF128( signZ, expZ - 5, sigZ.v64, sigZ.v0 );
 propagateNaN:
    uiZ = softfloat_propagateNaNF128UI( uiA64, uiA0, uiB64, uiB0 );
 uiZ:
    uZ.ui = uiZ;
    return uZ.f;
}

void
f128M_add( const float128_t *aPtr, const float128_t *bPtr, float128_t *zPtr )
{
    const uint64_t *aWPtr, *bWPtr;
    uint64_t uiA64, uiA0;
    int signA;
    uint64_t uiB64, uiB0;
    int signB;

    aWPtr = (const uint64_t *) aPtr;
    bWPtr = (const uint64_t *) bPtr;
    uiA64 = aWPtr[indexWord( 2, 1 )];
    uiA0  = aWPtr[indexWord( 2, 0 )];
    signA = signF128UI64( uiA64 );
    uiB64 = bWPtr[indexWord( 2, 1 )];
    uiB0  = bWPtr[indexWord( 2, 0 )];
    signB = signF128UI64( uiB64 );

    if ( signA == signB ) {
	*zPtr = softfloat_addMagsF128( uiA64, uiA0, uiB64, uiB0, signA );
    } else {
	*zPtr = softfloat_subMagsF128( uiA64, uiA0, uiB64, uiB0, signA );
    }
}

void
f128M_sub( const float128_t *aPtr, const float128_t *bPtr, float128_t *zPtr )
{
    const uint64_t *aWPtr, *bWPtr;
    uint_fast64_t uiA64, uiA0;
    bool signA;
    uint_fast64_t uiB64, uiB0;
    bool signB;

    aWPtr = (const uint64_t *) aPtr;
    bWPtr = (const uint64_t *) bPtr;
    uiA64 = aWPtr[indexWord( 2, 1 )];
    uiA0  = aWPtr[indexWord( 2, 0 )];
    signA = signF128UI64( uiA64 );
    uiB64 = bWPtr[indexWord( 2, 1 )];
    uiB0  = bWPtr[indexWord( 2, 0 )];
    signB = signF128UI64( uiB64 );
    if ( signA == signB ) {
	*zPtr = softfloat_subMagsF128( uiA64, uiA0, uiB64, uiB0, signA );
    } else {
	*zPtr = softfloat_addMagsF128( uiA64, uiA0, uiB64, uiB0, signA );
    }
}

void
f128M_mul( const float128_t *aPtr, const float128_t *bPtr, float128_t *zPtr )
{
    const uint32_t *aWPtr, *bWPtr;
    uint32_t *zWPtr;
    uint32_t uiA96;
    int32_t expA;
    uint32_t uiB96;
    int32_t expB;
    bool signZ;
    const uint32_t *ptr;
    uint32_t uiZ96, sigA[4];
    uint_fast8_t shiftDist;
    uint32_t sigB[4];
    int32_t expZ;
    uint32_t sigProd[8], *extSigZPtr;

    /*------------------------------------------------------------------------
    *------------------------------------------------------------------------*/
    aWPtr = (const uint32_t *) aPtr;
    bWPtr = (const uint32_t *) bPtr;
    zWPtr = (uint32_t *) zPtr;
    /*------------------------------------------------------------------------
    *------------------------------------------------------------------------*/
    uiA96 = aWPtr[indexWordHi( 4 )];
    expA = expF128UI96( uiA96 );
    uiB96 = bWPtr[indexWordHi( 4 )];
    expB = expF128UI96( uiB96 );
    signZ = signF128UI96( uiA96 ) ^ signF128UI96( uiB96 );
    /*------------------------------------------------------------------------
    *------------------------------------------------------------------------*/
    if ( (expA == 0x7FFF) || (expB == 0x7FFF) ) {
	if ( softfloat_tryPropagateNaNF128M( aWPtr, bWPtr, zWPtr ) ) return;
	ptr = aWPtr;
	if ( ! expA ) goto possiblyInvalid;
	if ( ! expB ) {
	    ptr = bWPtr;
     possiblyInvalid:
	    if (
		! fracF128UI96( ptr[indexWordHi( 4 )] )
		    && ! (ptr[indexWord( 4, 2 )] | ptr[indexWord( 4, 1 )]
			      | ptr[indexWord( 4, 0 )])
	    ) {
		softfloat_invalidF128M( zWPtr );
		return;
	    }
	}
	uiZ96 = packToF128UI96( signZ, 0x7FFF, 0 );
	goto uiZ96;
    }
    /*------------------------------------------------------------------------
    *------------------------------------------------------------------------*/
    if ( expA ) {
	sigA[indexWordHi( 4 )] = fracF128UI96( uiA96 ) | 0x00010000;
	sigA[indexWord( 4, 2 )] = aWPtr[indexWord( 4, 2 )];
	sigA[indexWord( 4, 1 )] = aWPtr[indexWord( 4, 1 )];
	sigA[indexWord( 4, 0 )] = aWPtr[indexWord( 4, 0 )];
    } else {
	expA = softfloat_shiftNormSigF128M( aWPtr, 0, sigA );
	if ( expA == -128 ) goto zero;
    }
    if ( expB ) {
	sigB[indexWordHi( 4 )] = fracF128UI96( uiB96 ) | 0x00010000;
	sigB[indexWord( 4, 2 )] = bWPtr[indexWord( 4, 2 )];
	sigB[indexWord( 4, 1 )] = bWPtr[indexWord( 4, 1 )];
	sigB[indexWord( 4, 0 )] = bWPtr[indexWord( 4, 0 )];
    } else {
	expB = softfloat_shiftNormSigF128M( bWPtr, 0, sigB );
	if ( expB == -128 ) goto zero;
    }
    /*------------------------------------------------------------------------
    *------------------------------------------------------------------------*/
    expZ = expA + expB - 0x4000;
    softfloat_mul128MTo256M( sigA, sigB, sigProd );
    if (
	sigProd[indexWord( 8, 2 )]
	    || (sigProd[indexWord( 8, 1 )] | sigProd[indexWord( 8, 0 )])
    ) {
	sigProd[indexWord( 8, 3 )] |= 1;
    }
    extSigZPtr = &sigProd[indexMultiwordHi( 8, 5 )];
    shiftDist = 16;
    if ( extSigZPtr[indexWordHi( 5 )] & 2 ) {
	++expZ;
	shiftDist = 15;
    }
    softfloat_shortShiftLeft160M( extSigZPtr, shiftDist, extSigZPtr );
    softfloat_roundPackMToF128M( signZ, expZ, extSigZPtr, zWPtr );
    return;
    /*------------------------------------------------------------------------
    *------------------------------------------------------------------------*/
 zero:
    uiZ96 = packToF128UI96( signZ, 0, 0 );
 uiZ96:
    zWPtr[indexWordHi( 4 )] = uiZ96;
    zWPtr[indexWord( 4, 2 )] = 0;
    zWPtr[indexWord( 4, 1 )] = 0;
    zWPtr[indexWord( 4, 0 )] = 0;
}

float128_t
f128_div( float128_t a, float128_t b )
{
    union ui128_f128 uA;
    uint_fast64_t uiA64, uiA0;
    bool signA;
    int_fast32_t expA;
    struct uint128 sigA;
    union ui128_f128 uB;
    uint_fast64_t uiB64, uiB0;
    bool signB;
    int_fast32_t expB;
    struct uint128 sigB;
    bool signZ;
    struct exp32_sig128 normExpSig;
    int_fast32_t expZ;
    struct uint128 rem;
    uint_fast32_t recip32;
    int ix;
    uint_fast64_t q64;
    uint_fast32_t q;
    struct uint128 term;
    uint_fast32_t qs[3];
    uint_fast64_t sigZExtra;
    struct uint128 sigZ, uiZ;
    union ui128_f128 uZ;

    /*------------------------------------------------------------------------
    *------------------------------------------------------------------------*/
    uA.f = a;
    uiA64 = uA.ui.v64;
    uiA0  = uA.ui.v0;
    signA = signF128UI64( uiA64 );
    expA  = expF128UI64( uiA64 );
    sigA.v64 = fracF128UI64( uiA64 );
    sigA.v0  = uiA0;
    uB.f = b;
    uiB64 = uB.ui.v64;
    uiB0  = uB.ui.v0;
    signB = signF128UI64( uiB64 );
    expB  = expF128UI64( uiB64 );
    sigB.v64 = fracF128UI64( uiB64 );
    sigB.v0  = uiB0;
    signZ = signA ^ signB;
    /*------------------------------------------------------------------------
    *------------------------------------------------------------------------*/
    if ( expA == 0x7FFF ) {
	if ( sigA.v64 | sigA.v0 ) goto propagateNaN;
	if ( expB == 0x7FFF ) {
	    if ( sigB.v64 | sigB.v0 ) goto propagateNaN;
	    goto invalid;
	}
	goto infinity;
    }
    if ( expB == 0x7FFF ) {
	if ( sigB.v64 | sigB.v0 ) goto propagateNaN;
	goto zero;
    }
    /*------------------------------------------------------------------------
    *------------------------------------------------------------------------*/
    if ( ! expB ) {
	if ( ! (sigB.v64 | sigB.v0) ) {
	    if ( ! (expA | sigA.v64 | sigA.v0) ) goto invalid;
	    softfloat_raiseFlags( softfloat_flag_infinite );
	    goto infinity;
	}
	normExpSig = softfloat_normSubnormalF128Sig( sigB.v64, sigB.v0 );
	expB = normExpSig.exp;
	sigB = normExpSig.sig;
    }
    if ( ! expA ) {
	if ( ! (sigA.v64 | sigA.v0) ) goto zero;
	normExpSig = softfloat_normSubnormalF128Sig( sigA.v64, sigA.v0 );
	expA = normExpSig.exp;
	sigA = normExpSig.sig;
    }
    /*------------------------------------------------------------------------
    *------------------------------------------------------------------------*/
    expZ = expA - expB + 0x3FFE;
    sigA.v64 |= UINT64_C( 0x0001000000000000 );
    sigB.v64 |= UINT64_C( 0x0001000000000000 );
    rem = sigA;
    if ( softfloat_lt128( sigA.v64, sigA.v0, sigB.v64, sigB.v0 ) ) {
	--expZ;
	rem = softfloat_add128( sigA.v64, sigA.v0, sigA.v64, sigA.v0 );
    }
    recip32 = softfloat_approxRecip32_1( sigB.v64>>17 );
    ix = 3;
    for (;;) {
	q64 = (uint_fast64_t) (uint32_t) (rem.v64>>19) * recip32;
	q = (q64 + 0x80000000)>>32;
	--ix;
	if ( ix < 0 ) break;
	rem = softfloat_shortShiftLeft128( rem.v64, rem.v0, 29 );
	term = softfloat_mul128By32( sigB.v64, sigB.v0, q );
	rem = softfloat_sub128( rem.v64, rem.v0, term.v64, term.v0 );
	if ( rem.v64 & UINT64_C( 0x8000000000000000 ) ) {
	    --q;
	    rem = softfloat_add128( rem.v64, rem.v0, sigB.v64, sigB.v0 );
	}
	qs[ix] = q;
    }
    /*------------------------------------------------------------------------
    *------------------------------------------------------------------------*/
    if ( ((q + 1) & 7) < 2 ) {
	rem = softfloat_shortShiftLeft128( rem.v64, rem.v0, 29 );
	term = softfloat_mul128By32( sigB.v64, sigB.v0, q );
	rem = softfloat_sub128( rem.v64, rem.v0, term.v64, term.v0 );
	if ( rem.v64 & UINT64_C( 0x8000000000000000 ) ) {
	    --q;
	    rem = softfloat_add128( rem.v64, rem.v0, sigB.v64, sigB.v0 );
	} else if ( softfloat_le128( sigB.v64, sigB.v0, rem.v64, rem.v0 ) ) {
	    ++q;
	    rem = softfloat_sub128( rem.v64, rem.v0, sigB.v64, sigB.v0 );
	}
	if ( rem.v64 | rem.v0 ) q |= 1;
    }
    /*------------------------------------------------------------------------
    *------------------------------------------------------------------------*/
    sigZExtra = (uint64_t) ((uint_fast64_t) q<<60);
    term = softfloat_shortShiftLeft128( 0, qs[1], 54 );
    sigZ =
	softfloat_add128(
	    (uint_fast64_t) qs[2]<<19, ((uint_fast64_t) qs[0]<<25) + (q>>4),
	    term.v64, term.v0
	);
    return
	softfloat_roundPackToF128( signZ, expZ, sigZ.v64, sigZ.v0, sigZExtra );
    /*------------------------------------------------------------------------
    *------------------------------------------------------------------------*/
 propagateNaN:
    uiZ = softfloat_propagateNaNF128UI( uiA64, uiA0, uiB64, uiB0 );
    goto uiZ;
    /*------------------------------------------------------------------------
    *------------------------------------------------------------------------*/
 invalid:
    softfloat_raiseFlags( softfloat_flag_invalid );
    softfloat_commonNaNToF128M( (uint32_t*)&uiZ );
    // uiZ.v64 = defaultNaNF128UI64;
    // uiZ.v0  = defaultNaNF128UI0;
    goto uiZ;
    /*------------------------------------------------------------------------
    *------------------------------------------------------------------------*/
 infinity:
    uiZ.v64 = packToF128UI64( signZ, 0x7FFF, 0 );
    goto uiZ0;
    /*------------------------------------------------------------------------
    *------------------------------------------------------------------------*/
 zero:
    uiZ.v64 = packToF128UI64( signZ, 0, 0 );
 uiZ0:
    uiZ.v0 = 0;
 uiZ:
    uZ.ui = uiZ;
    return uZ.f;
}

static float128_t
f128_rem( float128_t a, float128_t b )
{
    union ui128_f128 uA;
    uint_fast64_t uiA64, uiA0;
    bool signA;
    int_fast32_t expA;
    struct uint128 sigA;
    union ui128_f128 uB;
    uint_fast64_t uiB64, uiB0;
    int_fast32_t expB;
    struct uint128 sigB;
    struct exp32_sig128 normExpSig;
    struct uint128 rem;
    int_fast32_t expDiff;
    uint_fast32_t q, recip32;
    uint_fast64_t q64;
    struct uint128 term, altRem, meanRem;
    bool signRem;
    struct uint128 uiZ;
    union ui128_f128 uZ;

    /*------------------------------------------------------------------------
    *------------------------------------------------------------------------*/
    uA.f = a;
    uiA64 = uA.ui.v64;
    uiA0  = uA.ui.v0;
    signA = signF128UI64( uiA64 );
    expA  = expF128UI64( uiA64 );
    sigA.v64 = fracF128UI64( uiA64 );
    sigA.v0  = uiA0;
    uB.f = b;
    uiB64 = uB.ui.v64;
    uiB0  = uB.ui.v0;
    expB  = expF128UI64( uiB64 );
    sigB.v64 = fracF128UI64( uiB64 );
    sigB.v0  = uiB0;
    /*------------------------------------------------------------------------
    *------------------------------------------------------------------------*/
    if ( expA == 0x7FFF ) {
	if (
	    (sigA.v64 | sigA.v0) || ((expB == 0x7FFF) && (sigB.v64 | sigB.v0))
	) {
	    goto propagateNaN;
	}
	goto invalid;
    }
    if ( expB == 0x7FFF ) {
	if ( sigB.v64 | sigB.v0 ) goto propagateNaN;
	return a;
    }
    /*------------------------------------------------------------------------
    *------------------------------------------------------------------------*/
    if ( ! expB ) {
	if ( ! (sigB.v64 | sigB.v0) ) goto invalid;
	normExpSig = softfloat_normSubnormalF128Sig( sigB.v64, sigB.v0 );
	expB = normExpSig.exp;
	sigB = normExpSig.sig;
    }
    if ( ! expA ) {
	if ( ! (sigA.v64 | sigA.v0) ) return a;
	normExpSig = softfloat_normSubnormalF128Sig( sigA.v64, sigA.v0 );
	expA = normExpSig.exp;
	sigA = normExpSig.sig;
    }
    /*------------------------------------------------------------------------
    *------------------------------------------------------------------------*/
    sigA.v64 |= UINT64_C( 0x0001000000000000 );
    sigB.v64 |= UINT64_C( 0x0001000000000000 );
    rem = sigA;
    expDiff = expA - expB;
    if ( expDiff < 1 ) {
	if ( expDiff < -1 ) return a;
	if ( expDiff ) {
	    --expB;
	    sigB = softfloat_add128( sigB.v64, sigB.v0, sigB.v64, sigB.v0 );
	    q = 0;
	} else {
	    q = softfloat_le128( sigB.v64, sigB.v0, rem.v64, rem.v0 );
	    if ( q ) {
		rem = softfloat_sub128( rem.v64, rem.v0, sigB.v64, sigB.v0 );
	    }
	}
    } else {
	recip32 = softfloat_approxRecip32_1( sigB.v64>>17 );
	expDiff -= 30;
	for (;;) {
	    q64 = (uint_fast64_t) (uint32_t) (rem.v64>>19) * recip32;
	    if ( expDiff < 0 ) break;
	    q = (q64 + 0x80000000)>>32;
	    rem = softfloat_shortShiftLeft128( rem.v64, rem.v0, 29 );
	    term = softfloat_mul128By32( sigB.v64, sigB.v0, q );
	    rem = softfloat_sub128( rem.v64, rem.v0, term.v64, term.v0 );
	    if ( rem.v64 & UINT64_C( 0x8000000000000000 ) ) {
		rem = softfloat_add128( rem.v64, rem.v0, sigB.v64, sigB.v0 );
	    }
	    expDiff -= 29;
	}
	/*--------------------------------------------------------------------
	| (`expDiff' cannot be less than -29 here.)
	*--------------------------------------------------------------------*/
	q = (uint32_t) (q64>>32)>>(~expDiff & 31);
	rem = softfloat_shortShiftLeft128( rem.v64, rem.v0, expDiff + 30 );
	term = softfloat_mul128By32( sigB.v64, sigB.v0, q );
	rem = softfloat_sub128( rem.v64, rem.v0, term.v64, term.v0 );
	if ( rem.v64 & UINT64_C( 0x8000000000000000 ) ) {
	    altRem = softfloat_add128( rem.v64, rem.v0, sigB.v64, sigB.v0 );
	    goto selectRem;
	}
    }
    /*------------------------------------------------------------------------
    *------------------------------------------------------------------------*/
    do {
	altRem = rem;
	++q;
	rem = softfloat_sub128( rem.v64, rem.v0, sigB.v64, sigB.v0 );
    } while ( ! (rem.v64 & UINT64_C( 0x8000000000000000 )) );
 selectRem:
    meanRem = softfloat_add128( rem.v64, rem.v0, altRem.v64, altRem.v0 );
    if (
	(meanRem.v64 & UINT64_C( 0x8000000000000000 ))
	    || (! (meanRem.v64 | meanRem.v0) && (q & 1))
    ) {
	rem = altRem;
    }
    signRem = signA;
    if ( rem.v64 & UINT64_C( 0x8000000000000000 ) ) {
	signRem = ! signRem;
	rem = softfloat_sub128( 0, 0, rem.v64, rem.v0 );
    }
    return softfloat_normRoundPackToF128( signRem, expB - 1, rem.v64, rem.v0 );
    /*------------------------------------------------------------------------
    *------------------------------------------------------------------------*/
 propagateNaN:
    uiZ = softfloat_propagateNaNF128UI( uiA64, uiA0, uiB64, uiB0 );
    goto uiZ;
    /*------------------------------------------------------------------------
    *------------------------------------------------------------------------*/
 invalid:
    softfloat_raiseFlags( softfloat_flag_invalid );
    softfloat_commonNaNToF128M( (uint32_t*)&uiZ );
    // uiZ.v64 = defaultNaNF128UI64;
    // uiZ.v0  = defaultNaNF128UI0;
 uiZ:
    uZ.ui = uiZ;
    return uZ.f;
}

bool
f128_isNan( float128_t a) {
    int_fast32_t expA;
    struct uint128 sigA;
    uint_fast64_t uiA64, uiA0;
    union ui128_f128 uA;

    uA.f = a;
    uiA64 = uA.ui.v64;
    uiA0  = uA.ui.v0;
    return isNaNF128UI( uiA64, uiA0 );
}

bool
f128_isInf( float128_t a) {
    int_fast32_t expA;
    struct uint128 sigA;
    uint_fast64_t uiA64, uiA0;
    union ui128_f128 uA;

    uA.f = a;
    uiA64 = uA.ui.v64;
    uiA0  = uA.ui.v0;
    return isInfF128UI( uiA64, uiA0 );
//    expA  = expF128UI64( uiA64 );
//    sigA.v64 = fracF128UI64( uiA64 );
//    sigA.v0  = uiA0;
//
//    if ( expA == 0x7FFF ) {
//        if ( sigA.v64 == 0 ) return 1;
//    }
//    return 0;
}

static bool
f128M_eq( float128_t* aPtr, float128_t* bPtr )
{
    union ui128_f128 uA;
    uint_fast64_t uiA64, uiA0;
    union ui128_f128 uB;
    uint_fast64_t uiB64, uiB0;

    uA.f = *aPtr;
    uiA64 = uA.ui.v64;
    uiA0  = uA.ui.v0;
    uB.f = *bPtr;
    uiB64 = uB.ui.v64;
    uiB0  = uB.ui.v0;

    // fprintf(stderr, "isNaNF128UI(a): %d\n", isNaNF128UI( uiA64, uiA0 ));
    // fprintf(stderr, "isNaNF128UI(a): %d\n", isNaNF128UI( uiB64, uiB0 ));

    if ( isNaNF128UI( uiA64, uiA0 ) || isNaNF128UI( uiB64, uiB0 ) ) {
	if (
	       softfloat_isSigNaNF128UI( uiA64, uiA0 )
	    || softfloat_isSigNaNF128UI( uiB64, uiB0 )
	) {
	    softfloat_raiseFlags( softfloat_flag_invalid );
	}
	return 0;
    }

    // fprintf(stderr, "uiA0 == uiB0: %d\n", (uiA0 == uiB0));
    // fprintf(stderr, "uiA64 == uiB64: %d\n", (uiA64 == uiB64));

    return
	   (uiA0 == uiB0)
	&& (   (uiA64 == uiB64)
	    || (! uiA0 && ! ((uiA64 | uiB64) & UINT64_C( 0x7FFFFFFFFFFFFFFF )))
	   );
}

static bool
f128M_lt( float128_t* aPtr, float128_t* bPtr )
{
    union ui128_f128 uA;
    uint_fast64_t uiA64, uiA0;
    union ui128_f128 uB;
    uint_fast64_t uiB64, uiB0;
    bool signA, signB;

    uA.f = *aPtr;
    uiA64 = uA.ui.v64;
    uiA0  = uA.ui.v0;
    uB.f = *bPtr;
    uiB64 = uB.ui.v64;
    uiB0  = uB.ui.v0;
    if ( isNaNF128UI( uiA64, uiA0 ) || isNaNF128UI( uiB64, uiB0 ) ) {
	softfloat_raiseFlags( softfloat_flag_invalid );
	return 0;
    }
    signA = signF128UI64( uiA64 );
    signB = signF128UI64( uiB64 );
    return
	(signA != signB)
	    ? signA
		  && (((uiA64 | uiB64) & UINT64_C( 0x7FFFFFFFFFFFFFFF ))
			  | uiA0 | uiB0)
	    : ((uiA64 != uiB64) || (uiA0 != uiB0))
		  && (signA ^ softfloat_lt128( uiA64, uiA0, uiB64, uiB0 ));
}

static inline void
f128M_div( const float128_t *aPtr, const float128_t *bPtr, float128_t *zPtr )
{
    *zPtr = f128_div( *aPtr, *bPtr );
}

static inline void
f128M_rem( const float128_t *aPtr, const float128_t *bPtr, float128_t *zPtr )
{
    *zPtr = f128_rem( *aPtr, *bPtr );
}

static inline float128_t
f128_negate( float128_t a) {
    struct uint128 sigA;
    uint_fast64_t uiA64, uiA0;
    union ui128_f128 uA;

    uA.f = a;
    uiA64 = uA.ui.v64;
    uiA0  = uA.ui.v0;
    if (isNaNF128UI( uiA64, uiA0 )) return a;
    uA.ui.v64 ^= UINT64_C(0x8000000000000000);
    return uA.f;
}

static inline void
f128M_negate( const float128_t *aPtr, float128_t *zPtr )
{
    *zPtr = f128_negate( *aPtr );
}

#endif // SUPPORT_QUADFLOAT

%}
! !

!QuadFloat class methodsFor:'documentation'!

documentation
"
    QuadFloats represent rational numbers with limited precision
    and are mapped to IEEE quadruple precision format (128bit),
    also called binary128.

    If the underlying cpu supports them natively, the machine format (long double) is
    used. Otherwise, a software emulation is done, which is much slower.
    Thus only use them, if you really need the additional precision;
    if not, use Float (which are doubles) or LongFloats which usually have IEEE extended precision (80bit).

    QuadFloats give you definite 128 bit quadruple floats,
    thus, code using quadFloats is guaranteed to be portable from one architecture to another.

    Representation:
	    128bit quadruple IEEE floats (16bytes);
	    112 bit mantissa,
	    16 bit exponent,
	    34 decimal digits (approx.)

    On Sparc CPUs, this is a native supported type (long double) and fast;
    on x86 CPUs, this is emulated and slow.

    Mixed mode arithmetic:
	quadFloat op anyFloat    -> quadFloat
	anyFloat op quadFloat    -> quadFloat

    Range and precision of storage formats: see LimitedPrecisionReal >> documentation

    [author:]
	Claus Gittinger

    [see also:]
	Number
	Float ShortFloat LongFloat Fraction FixedPoint Integer Complex
	FloatArray DoubleArray
	https://en.wikipedia.org/wiki/Extended_precision
"
! !

!QuadFloat class methodsFor:'instance creation'!

basicNew
    "return a new quadFloat - here we return 0.0
     - QuadFloats are usually NOT created this way ...
     Its implemented here to allow things like binary store & load
     of quadFloats.
     (but it is not a good idea to store the bits of a float - the reader might have a
      totally different representation - so floats should be
      binary stored in a device independent format)."

%{  /* NOCONTEXT */
#ifdef SUPPORT_QUADFLOAT
    OBJ newFloat;

    if (sizeof(long double) == sizeof(float128_t)) {
	__qMKLFLOAT(newFloat, 0.0);   /* OBJECT ALLOCATION */
    } else {
	float128_t qf;
	f64_to_f128M(0.0, &qf);
	__qMKQFLOAT(newFloat, qf);   /* OBJECT ALLOCATION */
    }
    RETURN (newFloat);
#endif /* SUPPORT_QUADFLOAT */
%}.
    self error:'QuadFloats not supported on this patform'

    "Created: / 06-06-2019 / 17:18:58 / Claus Gittinger"
!

fromFloat:aFloat
    "return a new quadFloat, given a float value"

%{  /* NOCONTEXT */
#ifdef SUPPORT_QUADFLOAT
    OBJ newFloat;

    if (__isFloatLike(aFloat)) {
	float64_t f = __floatVal(aFloat);
	float128_t qf;

	f64_to_f128M(f, &qf);
	__qMKQFLOAT(newFloat, qf);   /* OBJECT ALLOCATION */
	RETURN (newFloat);
    }
#endif /* SUPPORT_QUADFLOAT */
%}.
    aFloat isFloat ifTrue:[
	self errorUnsupported
    ].
    ArgumentError raise

    "
     QuadFloat fromFloat:123.0
     123.0 asQuadFloat
     123 asQuadFloat
    "

    "Created: / 06-06-2019 / 18:01:03 / Claus Gittinger"
!

fromInteger:anInteger
    "return a new quadFloat, given an integer value"

%{  /* NOCONTEXT */
#ifdef SUPPORT_QUADFLOAT
    OBJ newFloat;

    if (__isSmallInteger(anInteger)) {
	INT iVal = __intVal(anInteger);
	float128_t qf;

#if __POINTER_SIZE__ == 4
	i32_to_f128M( iVal, &qf );
else
	i64_to_f128M( iVal, &qf );
#endif
	__qMKQFLOAT(newFloat, qf);   /* OBJECT ALLOCATION */
	RETURN (newFloat);
    }
#endif /* SUPPORT_QUADFLOAT */
%}.
    ^ super fromInteger:anInteger

    "
     QuadFloat fromInteger:123
     123 asQuadFloat
    "
!

fromLongFloat:aFloat
    "return a new quadFloat, given a long float value"

%{  /* NOCONTEXT */
#ifdef SUPPORT_QUADFLOAT
    OBJ newFloat;
    union {
	LONGFLOAT_t lf;         // is long double
	extFloat80_t ef;        // is 80bit ext
	float128_t qf;          // is 128bit quad
    } u;

    if (__isLongFloat(aFloat)) {
	u.lf = __longFloatVal(aFloat);

	if (sizeof(LONGFLOAT_t) == 16) {
	    // longFloat is already 128 bits in size (sparc)
	    __qMKQFLOAT(newFloat, u.qf);   /* OBJECT ALLOCATION */
	    RETURN (newFloat);
	}
	if (sizeof(LONGFLOAT_t) < 16) {
	    // assume 80bit extended float format (amd64, x86_64)
	    u.qf = extF80_to_f128( u.ef);
	    __qMKQFLOAT(newFloat, u.qf);   /* OBJECT ALLOCATION */
	    RETURN (newFloat);
	}
	// fall into error case
    }
#endif /* SUPPORT_QUADFLOAT */
%}.
    aFloat isLongFloat ifTrue:[
	self errorUnsupported
    ].
    ArgumentError raise

    "
     QuadFloat fromLongFloat:123.0 asLongFloat
    "
!

fromShortFloat:aShortFloat
    "return a new quadFloat, given a float value"

    ^ self fromFloat:(aShortFloat asFloat)

    "
     QuadFloat fromShortFloat:123.0 asShortFloat
    "

    "Created: / 08-06-2019 / 03:28:37 / Claus Gittinger"
! !

!QuadFloat class methodsFor:'coercing & converting'!

coerce:aNumber
    "convert the argument aNumber into an instance of the receiver's class and return it."

    ^ aNumber asQuadFloat.

    "Created: / 06-06-2019 / 16:51:01 / Claus Gittinger"
! !

!QuadFloat class methodsFor:'constants'!

NaN
    "return a quadFloat which represents not-a-Number (i.e. an invalid number)"

    |nan|

    NaN isNil ifTrue:[
%{  /* NOCONTEXT */
#ifdef SUPPORT_QUADFLOAT
	{
	    OBJ newFloat;
	    float128_t qf;

	    softfloat_commonNaNToF128M( (uint32_t*)(&qf) );
	    __qMKQFLOAT(newFloat, qf);   /* OBJECT ALLOCATION */
	    nan = newFloat;
	}
#endif /* SUPPORT_QUADFLOAT */
%}.
	NaN := nan
    ].
    ^ NaN
!

e
    "return the constant e as quadFloat"

    E isNil ifTrue:[
	"/ eDigits has enough digits for 128bit IEEE quads
	"/ do not use as a literal constant here - we cannot depend on the underlying C-compiler here...
	E  := self readFrom:(Number eDigits)
    ].
    ^ E

    "Created: / 06-06-2019 / 17:01:54 / Claus Gittinger"
!

infinity
    "return a quadFloat which represents +INF"

    |inf|

    PositiveInfinity isNil ifTrue:[
%{  /* NOCONTEXT */
#ifdef SUPPORT_QUADFLOAT
	{
	    OBJ newFloat;
	    struct uint128 uiZ;
	    union ui128_f128 uZ;
	    float128_t qf;

	    uiZ.v64 = packToF128UI64( 0, 0x7FFF, 0 );
	    uiZ.v0 = 0;
	    uZ.ui = uiZ;
	    qf = uZ.f;
	    __qMKQFLOAT(newFloat, qf);   /* OBJECT ALLOCATION */
	    inf = newFloat;
	}
#endif /* SUPPORT_QUADFLOAT */
%}.
	PositiveInfinity := inf
    ].
    ^ PositiveInfinity

    "Created: / 08-06-2019 / 14:05:26 / Claus Gittinger"
!

negativeInfinity
    "return a quadFloat which represents -INF"

    |inf|

    NegativeInfinity isNil ifTrue:[
%{  /* NOCONTEXT */
#ifdef SUPPORT_QUADFLOAT
	{
	    OBJ newFloat;
	    struct uint128 uiZ;
	    union ui128_f128 uZ;
	    float128_t qf;

	    uiZ.v64 = packToF128UI64( 1, 0x7FFF, 0 );
	    uiZ.v0 = 0;
	    uZ.ui = uiZ;
	    qf = uZ.f;
	    __qMKQFLOAT(newFloat, qf);   /* OBJECT ALLOCATION */
	    inf = newFloat;
	}
#endif /* SUPPORT_QUADFLOAT */
%}.
	NegativeInfinity := inf
    ].
    ^ NegativeInfinity

    "Created: / 08-06-2019 / 14:05:50 / Claus Gittinger"
!

pi
    "return the constant pi as quadFloat"

    Pi isNil ifTrue:[
	"/ piDigits has enough digits for 128bit IEEE quads
	"/ do not use as a literal constant here - we cannot depend on the underlying C-compiler here...
	Pi  := self readFrom:(Number piDigits)
    ].
    ^ Pi

    "Created: / 06-06-2019 / 17:09:51 / Claus Gittinger"
!

unity
    "return the neutral element for multiplication (1.0) as QuadFloat"

    QuadFloatOne isNil ifTrue:[
	QuadFloatOne := 1.0 asQuadFloat.
    ].
    ^ QuadFloatOne

    "Created: / 07-06-2019 / 03:26:38 / Claus Gittinger"
!

zero
    "return the neutral element for addition (0.0) as QuadFloat"

    QuadFloatZero isNil ifTrue:[
	QuadFloatZero := 0.0 asQuadFloat
    ].
    ^ QuadFloatZero

    "Created: / 07-06-2019 / 09:22:56 / Claus Gittinger"
! !

!QuadFloat class methodsFor:'error reportng'!

errorUnsupported
    self error:'QuadFloats not supported on this patform'

    "Created: / 07-06-2019 / 02:44:39 / Claus Gittinger"
! !

!QuadFloat class methodsFor:'queries'!

epsilon
    "return the maximum relative spacing of instances of mySelf
     (i.e. the value-delta of the least significant bit)"

    Epsilon isNil ifTrue:[
	Epsilon := self computeEpsilon.
    ].
    ^ Epsilon

    "
     self epsilon
    "

    "Created: / 10-06-2019 / 21:21:18 / Claus Gittinger"
!

exponentCharacter
    "return the character used to print between mantissa an exponent.
     Also used by the scanner when reading numbers."

    ^ $Q

    "Created: / 10-06-2019 / 21:28:04 / Claus Gittinger"
!

numBitsInExponent
    "answer the number of bits in the exponent
     the hidden bit is not counted here:

     This is an 128bit quadfloat,
	   where 15 bits are available in the exponent:
	seeeeeee eeeeeeee mmmmmmmm mmmmmmmm mmmmmmmm mmmmmmmm mmmmmmmm mmmmmmmm mmmmmmmm mmmmmmmm mmmmmmmm mmmmmmmm mmmmmmmm mmmmmmmm mmmmmmmm mmmmmmmm...
    "

    ^ 15

    "
     1.0 class numBitsInExponent -> 11
     1.0 asShortFloat class numBitsInExponent -> 8
     1.0 asLongFloat class numBitsInExponent -> 15
     1.0 asQuadFloat class numBitsInExponent -> 15
    "

    "Created: / 11-06-2019 / 00:14:55 / Claus Gittinger"
!

numBitsInMantissa
    "answer the number of bits in the mantissa
     the hidden bit is not counted here:

     This is an 128bit quadfloat,
	   where 112 bits are available in the mantissa:
	seeeeeee eeeeeeee mmmmmmmm mmmmmmmm mmmmmmmm mmmmmmmm mmmmmmmm mmmmmmmm mmmmmmmm mmmmmmmm mmmmmmmm mmmmmmmm mmmmmmmm mmmmmmmm mmmmmmmm mmmmmmmm...
    "

    ^ 112

    "
     1.0 class numBitsInMantissa
     1.0 asShortFloat class numBitsInMantissa
     1.0 asLongFloat class numBitsInMantissa
     1.0 asQuadFloat class numBitsInMantissa
    "

    "Created: / 07-06-2019 / 03:24:20 / Claus Gittinger"
! !

!QuadFloat methodsFor:'arithmetic'!

* aNumber
    "return the product of the receiver and the argument."

thisContext isRecursive ifTrue:[self halt].
    ^ aNumber productFromQuadFloat:self
!

+ aNumber
    "return the sum of the receiver and the argument, aNumber"

    ^ aNumber sumFromQuadFloat:self
!

- aNumber
    "return the difference of the receiver and the argument, aNumber"

    ^ aNumber differenceFromQuadFloat:self
!

/ aNumber
    "return the quotient of the receiver and the argument, aNumber"

    aNumber isZero ifTrue:[
	"
	 No, you shalt not divide by zero
	"
	^ ZeroDivide raiseRequestWith:thisContext.
    ].
    ^ aNumber quotientFromQuadFloat:self
!

negated
    "return the receiver negated"

%{  /* NOCONTEXT */
#ifdef SUPPORT_QUADFLOAT
    OBJ newFloat;
    float128_t result, myVal;

    myVal = __quadFloatVal(self);
    f128M_negate( &myVal, &result );
    __qMKQFLOAT(newFloat, result);
    RETURN ( newFloat );
#endif
%}.
!

rem: aNumber
    "return the floating point remainder of the receiver and the argument, aNumber"

    aNumber isZero ifTrue:[
	"
	 No, you shalt not divide by zero
	"
	^ ZeroDivide raiseRequestWith:thisContext.
    ].
    ^ aNumber remainderFromLongFloat:self
! !

!QuadFloat methodsFor:'coercing & converting'!

generality
    "return the generality value - see ArithmeticValue>>retry:coercing:"

    ^ 93

    "Created: / 07-06-2019 / 09:30:58 / Claus Gittinger"
! !

!QuadFloat methodsFor:'comparing'!

< aNumber
    "return true, if the argument is greater"

    ^ aNumber lessFromQuadFloat:self

    "Created: / 07-06-2019 / 09:25:47 / Claus Gittinger"
!

= aNumber
    "return true, if the argument represents the same numeric value
     as the receiver, false otherwise"

    ^ aNumber equalFromQuadFloat:self

    "Created: / 07-06-2019 / 09:25:27 / Claus Gittinger"
!

hash
    "return a number for hashing; redefined, since floats compare
     by numeric value (i.e. 3.0 = 3), therefore 3.0 hash must be the same
     as 3 hash."

    |i|

    (self >= SmallInteger minVal and:[self <= SmallInteger maxVal]) ifTrue:[
	i := self asInteger.
	self = i ifTrue:[
	    ^ i hash
	].
    ].

    ^ self asFloat hash

    "
     1.2345 hash
     1.2345 asShortFloat hash
     1.2345 asLongFloat hash
     1.2345 asQuadFloat hash

     1.0 hash
     1.0 asShortFloat hash
     1.0 asLongFloat hash
     1.0 asQuadFloat hash

     0.5 asShortFloat hash
     0.5 asShortFloat hash
     0.5 asLongFloat hash
     0.5 asQuadFloat hash

     0.25 asShortFloat hash
     0.25 asShortFloat hash
     0.25 asLongFloat hash
     0.25 asQuadFloat hash
    "

    "Created: / 07-06-2019 / 09:28:07 / Claus Gittinger"
! !

!QuadFloat methodsFor:'double dispatching'!

differenceFromQuadFloat:aQuadFloat
%{
#ifdef SUPPORT_QUADFLOAT
    OBJ newFloat;
    float128_t result, myVal, argVal;

    myVal = __quadFloatVal(self);
    argVal = __quadFloatVal(aQuadFloat);
    f128M_sub( &myVal, &argVal, &result );
    __qMKQFLOAT(newFloat, result);
    RETURN ( newFloat );
#endif // SUPPORT_QUADFLOAT
%}.
    self errorUnsupported
!

equalFromQuadFloat:aQuadFloat
    "sent when aQuadFloat does not know how to compare agaist the receiver, self"

%{
#ifdef SUPPORT_QUADFLOAT
    OBJ newFloat;
    float128_t result, myVal, argVal;

    myVal = __quadFloatVal(self);
    argVal = __quadFloatVal(aQuadFloat);
    RETURN (f128M_eq( &argVal, &myVal ) ? true : false);
#endif // SUPPORT_QUADFLOAT
%}.
    self errorUnsupported

    "Modified: / 08-06-2019 / 13:31:48 / Claus Gittinger"
!

lessFromQuadFloat:aQuadFloat
    "sent when aQuadFloat does not know how to compare agaist the receiver, self"

%{
#ifdef SUPPORT_QUADFLOAT
    OBJ newFloat;
    float128_t result, myVal, argVal;

    myVal = __quadFloatVal(self);
    argVal = __quadFloatVal(aQuadFloat);
    RETURN (f128M_lt( &argVal, &myVal ) ? true : false);
#endif // SUPPORT_QUADFLOAT
%}.
    self errorUnsupported
!

productFromQuadFloat:aQuadFloat
    "sent when aQuadFloat does not know how to multiply the receiver, self"

%{
#ifdef SUPPORT_QUADFLOAT
    OBJ newFloat;
    float128_t result, myVal, argVal;

    myVal = __quadFloatVal(self);
    argVal = __quadFloatVal(aQuadFloat);
    f128M_mul( &myVal, &argVal, &result );
    __qMKQFLOAT(newFloat, result);
    RETURN ( newFloat );
#endif // SUPPORT_QUADFLOAT
%}.
    self errorUnsupported
!

quotientFromQuadFloat:aQuadFloat
    "sent when aQuadFloat does not know how to multiply the receiver, self"

%{
#ifdef SUPPORT_QUADFLOAT
    OBJ newFloat;
    float128_t result, myVal, argVal;

    myVal = __quadFloatVal(self);
    argVal = __quadFloatVal(aQuadFloat);
    f128M_div( &myVal, &argVal, &result );
    __qMKQFLOAT(newFloat, result);
    RETURN ( newFloat );
#endif // SUPPORT_QUADFLOAT
%}.
    self errorUnsupported
!

sumFromQuadFloat:aQuadFloat
    "sent when aQuadFloat does not know how to add the receiver, self"

%{
#ifdef SUPPORT_QUADFLOAT
    OBJ newFloat;
    float128_t result, myVal, argVal;

    myVal = __quadFloatVal(self);
    argVal = __quadFloatVal(aQuadFloat);
    f128M_add( &myVal, &argVal, &result );
    __qMKQFLOAT(newFloat, result);
    RETURN ( newFloat );
#endif // SUPPORT_QUADFLOAT
%}.
    self errorUnsupported
! !

!QuadFloat methodsFor:'error reportng'!

errorUnsupported
    self class errorUnsupported

    "Modified: / 07-06-2019 / 02:44:51 / Claus Gittinger"
! !

!QuadFloat methodsFor:'printing'!

printOn:aStream
    |mantissa exponent|

    mantissa := self mantissa.
    exponent := self exponent.

    self exponent == 0 ifTrue:[
	mantissa printOn:aStream.
	aStream nextPutAll:'.0'.
	^ self
    ].
    mantissa == 0 ifTrue:[
	"/ a zero mantissa is impossible - except for zero and a few others
	exponent == 0 ifTrue:[ aStream nextPutAll:'0.0'. ^ self].
	self == NaN ifTrue:[ aStream nextPutAll:'NAN'. ^ self ].
	self == NegativeInfinity ifTrue:[ aStream nextPutAll:'-INF'. ^ self].
	self == PositiveInfinity ifTrue:[ aStream nextPutAll:'INF'. ^ self].

	self error:'invalid largeFloat' mayProceed:true.
	aStream nextPutAll:'Invalid'. ^ self.
    ].

    exponent >= 0 ifTrue:[
	(mantissa bitShift:exponent) printOn:aStream.
	aStream nextPutAll:'.0'.
	^ self
    ].
    ((mantissa / (1 bitShift:exponent negated)) asFixedPoint:6) printOn:aStream.

    "Created: / 11-06-2019 / 00:13:00 / Claus Gittinger"
! !

!QuadFloat methodsFor:'queries'!

isFinite
    "return true, if the receiver is a finite float
     i.e. not NaN and not infinite."

%{  /* NOCONTEXT */
#ifdef SUPPORT_QUADFLOAT
    float128_t myVal;

    myVal = __quadFloatVal(self);
    RETURN ((f128_isInf(myVal) || f128_isNan(myVal)) ? false : true);
#endif // SUPPORT_QUADFLOAT
%}.

    "
	1.0 isFinite
	self NaN isFinite
	self infinity isFinite
	self negativeInfinity isFinite
	(0.0 uncheckedDivide: 0.0) isFinite
	(1.0 uncheckedDivide: 0.0) isFinite
    "
!

isInfinite
    "return true, if the receiver is an infinite float.
     i.e. +Inf or -Inf."

%{  /* NOCONTEXT */
#ifdef SUPPORT_QUADFLOAT
    float128_t myVal;

    myVal = __quadFloatVal(self);
    RETURN (f128_isInf(myVal) ? true : false);
#endif // SUPPORT_QUADFLOAT
%}.

    "
	1.0 isFinite
	1.0 isInfinite
	self NaN isFinite
	self NaN isInfinite
	self infinity isFinite
	self infinity isInfinite
	self negativeInfinity isFinite
	self negativeInfinity isInfinite
	(0.0 uncheckedDivide: 0.0) isFinite
	(1.0 uncheckedDivide: 0.0) isFinite
    "
!

isNaN
    "return true, if the receiver is a finite float
     i.e. not NaN and not infinite."

%{  /* NOCONTEXT */
#ifdef SUPPORT_QUADFLOAT
    float128_t myVal;

    myVal = __quadFloatVal(self);
    RETURN (f128_isNan(myVal) ? true : false);
#endif // SUPPORT_QUADFLOAT
%}.

    "
	1.0 isFinite
	self NaN isFinite
	self infinity isFinite
	(0.0 uncheckedDivide: 0.0) isFinite
	(1.0 uncheckedDivide: 0.0) isFinite
    "
! !

!QuadFloat class methodsFor:'documentation'!

version_CVS
    ^ '$Header$'
! !