"{ Package: 'stx:libbasic2' }"

"{ NameSpace: Smalltalk }"

LimitedPrecisionReal variableByteSubclass:#QuadFloat

	instanceVariableNames:''

	classVariableNames:'QuadFloatZero QuadFloatOne Pi E Epsilon NaN PositiveInfinity

		NegativeInfinity Halfpi HalfpiNegative Phi'

	poolDictionaries:''

	category:'Magnitude-Numbers'

!

!QuadFloat primitiveDefinitions!

%{

#include <math.h>

extern float128_t STX_addq(float128_t, float128_t, int);

extern float128_t STX_mulq(float128_t, float128_t);

extern float128_t STX_divq(float128_t, float128_t);

extern float128_t STX_negq(float128_t);

extern float128_t STX_absq(float128_t);

extern float128_t STX_floorq(float128_t);

extern float128_t STX_frexpq(float128_t, int*);

extern float128_t STX_ceilq(float128_t);

extern float128_t STX_logq(float128_t);

extern float128_t STX_log10q(float128_t);

extern float128_t STX_log2q(float128_t);

extern float128_t STX_expq(float128_t);

extern float128_t STX_sinq(float128_t);

extern float128_t STX_cosq(float128_t);

extern float128_t STX_tanq(float128_t);

extern float128_t STX_sinhq(float128_t);

extern float128_t STX_coshq(float128_t);

extern float128_t STX_tanhq(float128_t);

extern float128_t STX_asinq(float128_t);

extern float128_t STX_acosq(float128_t);

extern float128_t STX_atanq(float128_t);

extern float128_t STX_asinhq(float128_t);

extern float128_t STX_acoshq(float128_t);

extern float128_t STX_atanhq(float128_t);

extern float128_t STX_qZero;

extern float128_t STX_dbltoq(double);

extern float128_t STX_inttoq(long);

extern double STX_qtodbl(float128_t);

extern int STX_isNanq(float128_t*);

extern int STX_prcmpq(float128_t*, float128_t*);

#define STX_isfiniteq(q)    (!STX_isNanq(&(q)) && !STX_isInfq(&(q)))

#define STX_eqq(x1, x2)     (STX_prcmpq (&(x1), &(x2)) == 0)

#define STX_neqq(x1, x2)    (STX_prcmpq (&(x1), &(x2)) != 0)

#define STX_gtq(x1, x2)     (STX_prcmpq (&(x1), &(x2)) > 0)

#define STX_geq(x1, x2)     (STX_prcmpq (&(x1), &(x2)) >= 0)

#define STX_ltq(x1, x2)     (STX_prcmpq (&(x1), &(x2)) < 0)

#define STX_leq(x1, x2)     (STX_prcmpq (&(x1), &(x2)) <= 0)

#ifdef __win32__

/*

 * no finite(x) ?

 * no isnan(x) ?

 */

# ifndef isnan

#  define isnan(x)      \

	((((unsigned int *)(&x))[0] == 0x00000000) && \

	 (((unsigned int *)(&x))[1] == 0xFFF80000))

# endif

# ifndef isPositiveInfinity

#  define isPositiveInfinity(x) \

	((((unsigned int *)(&x))[0] == 0x00000000) && \

	 (((unsigned int *)(&x))[1] == 0x7FF00000))

# endif

# ifndef isNegativeInfinity

#  define isNegativeInfinity(x) \

	((((unsigned int *)(&x))[0] == 0x00000000) && \

	 (((unsigned int *)(&x))[1] == 0xFFF00000))

# endif

# ifndef isinf

#  define isinf(x) \

	((((unsigned int *)(&x))[0] == 0x00000000) && \

	 ((((unsigned int *)(&x))[1] & 0x7FF00000) == 0x7FF00000))

# endif

# ifndef isfinite

#  define isfinite(x) (!isinf(x) && !isnan(x))

# endif

#else // not win32

#endif

%}

! !

!QuadFloat primitiveVariables!

%{

%}

! !

!QuadFloat primitiveFunctions!

%{

%}

! !

!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 = STX_qZero;

	__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)) {

	double f = __floatVal(aFloat);

	float128_t qf;

	qf = STX_dbltoq (f);

	__qMKQFLOAT(newFloat, qf);   /* OBJECT ALLOCATION */

	RETURN (newFloat);

    }

%}.

    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;

	qf = STX_inttoq( (long)iVal );

	__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 xSUPPORT_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 xSUPPORT_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 xSUPPORT_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 xSUPPORT_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"

!

phi

    "return the constant phi as quadFloat"

    Phi isNil ifTrue:[

	"/ phiDigits has enough digits for 128bit IEEE quads

	"/ do not use as a literal constant here - we cannot depend on the underlying C-compiler here...

	Phi  := self readFrom:(Number phiDigits)

    ].

    ^ Phi

!

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.

     This is an IEEE 128bit quadfloat, where 15 bits are available in the exponent:

    "

    ^ 15