"{ 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 )]