QuadFloat.st
changeset 5269 83163a7005f1
parent 5268 d3a464884f60
child 5276 b81062a163d4
equal deleted inserted replaced
5268:d3a464884f60 5269:83163a7005f1
     5 "{ NameSpace: Smalltalk }"
     5 "{ NameSpace: Smalltalk }"
     6 
     6 
     7 LimitedPrecisionReal variableByteSubclass:#QuadFloat
     7 LimitedPrecisionReal variableByteSubclass:#QuadFloat
     8 	instanceVariableNames:''
     8 	instanceVariableNames:''
     9 	classVariableNames:'QuadFloatZero QuadFloatOne Pi E Epsilon NaN PositiveInfinity
     9 	classVariableNames:'QuadFloatZero QuadFloatOne Pi E Epsilon NaN PositiveInfinity
    10 		NegativeInfinity Halfpi HalfpiNegative Phi DefaultPrintFormat
    10 		NegativeInfinity Halfpi HalfpiNegative Phi'
    11 		DefaultPrintfFormat'
       
    12 	poolDictionaries:''
    11 	poolDictionaries:''
    13 	category:'Magnitude-Numbers'
    12 	category:'Magnitude-Numbers'
    14 !
    13 !
    15 
    14 
    16 !QuadFloat primitiveDefinitions!
    15 !QuadFloat primitiveDefinitions!
    17 %{
    16 %{
    18 #define SUPPORT_QUADFLOAT
    17 #if __POINTER_SIZE__ == 8
    19 
    18 # define SUPPORT_QUADFLOAT
    20 #include <math.h>
       
    21 
       
    22 extern float128_t STX_qadd(float128_t, float128_t, int);
       
    23 extern float128_t STX_qmul(float128_t, float128_t);
       
    24 extern float128_t STX_qdiv(float128_t, float128_t);
       
    25 extern float128_t STX_qneg(float128_t);
       
    26 extern float128_t STX_qabs(float128_t);
       
    27 extern float128_t STX_qfloor(float128_t);
       
    28 extern float128_t STX_qfrexp(float128_t, int*);
       
    29 extern float128_t STX_qceil(float128_t);
       
    30 extern float128_t STX_qlog(float128_t);
       
    31 extern float128_t STX_qlog10(float128_t);
       
    32 extern float128_t STX_qlog2(float128_t);
       
    33 extern float128_t STX_qexp(float128_t);
       
    34 extern float128_t STX_qsin(float128_t);
       
    35 extern float128_t STX_qcos(float128_t);
       
    36 extern float128_t STX_qtan(float128_t);
       
    37 extern float128_t STX_qsinh(float128_t);
       
    38 extern float128_t STX_qcosh(float128_t);
       
    39 extern float128_t STX_qtanh(float128_t);
       
    40 extern float128_t STX_qasin(float128_t);
       
    41 extern float128_t STX_qacos(float128_t);
       
    42 extern float128_t STX_qatan(float128_t);
       
    43 extern float128_t STX_qasinh(float128_t);
       
    44 extern float128_t STX_qacosh(float128_t);
       
    45 extern float128_t STX_qatanh(float128_t);
       
    46 extern float128_t STX_qZero;
       
    47 extern float128_t STX_dbltoq(double);
       
    48 extern float128_t STX_inttoq(long);
       
    49 extern double STX_qtodbl(float128_t);
       
    50 extern int STX_qisNan(float128_t*);
       
    51 extern int STX_qprcmp(float128_t*, float128_t*);
       
    52 
       
    53 #define STX_qisfinite(q)    (!STX_qisNan(&(q)) && !STX_qisInf(&(q)))
       
    54 #define STX_qeq(x1, x2)     (STX_qprcmp (&(x1), &(x2)) == 0)
       
    55 #define STX_qneq(x1, x2)    (STX_qprcmp (&(x1), &(x2)) != 0)
       
    56 #define STX_qgt(x1, x2)     (STX_qprcmp (&(x1), &(x2)) > 0)
       
    57 #define STX_qge(x1, x2)     (STX_qprcmp (&(x1), &(x2)) >= 0)
       
    58 #define STX_qlt(x1, x2)     (STX_qprcmp (&(x1), &(x2)) < 0)
       
    59 #define STX_qle(x1, x2)     (STX_qprcmp (&(x1), &(x2)) <= 0)
       
    60 
       
    61 #ifdef __win32__
       
    62 /*
       
    63  * no finite(x) ?
       
    64  * no isnan(x) ?
       
    65  */
       
    66 # ifndef isnan
       
    67 #  define isnan(x)      \
       
    68 	((((unsigned int *)(&x))[0] == 0x00000000) && \
       
    69 	 (((unsigned int *)(&x))[1] == 0xFFF80000))
       
    70 # endif
       
    71 
       
    72 # ifndef isPositiveInfinity
       
    73 #  define isPositiveInfinity(x) \
       
    74 	((((unsigned int *)(&x))[0] == 0x00000000) && \
       
    75 	 (((unsigned int *)(&x))[1] == 0x7FF00000))
       
    76 # endif
       
    77 
       
    78 # ifndef isNegativeInfinity
       
    79 #  define isNegativeInfinity(x) \
       
    80 	((((unsigned int *)(&x))[0] == 0x00000000) && \
       
    81 	 (((unsigned int *)(&x))[1] == 0xFFF00000))
       
    82 # endif
       
    83 
       
    84 # ifndef isinf
       
    85 #  define isinf(x) \
       
    86 	((((unsigned int *)(&x))[0] == 0x00000000) && \
       
    87 	 ((((unsigned int *)(&x))[1] & 0x7FF00000) == 0x7FF00000))
       
    88 # endif
       
    89 
       
    90 # ifndef isfinite
       
    91 #  define isfinite(x) (!isinf(x) && !isnan(x))
       
    92 # endif
       
    93 
       
    94 #else // not win32
       
    95 #endif
    19 #endif
       
    20 
       
    21 #if defined(__MSBFIRST__) || defined(__MSB_FIRST__) || defined(MSBFIRST)
       
    22 # undef LITTLEENDIAN
       
    23 #else
       
    24 # define LITTLEENDIAN
       
    25 #endif
       
    26 
       
    27 #ifdef SUPPORT_QUADFLOAT
       
    28 # include <stdint.h>
       
    29 
       
    30 // The following code is adapted from the Softfloat-3e package,
       
    31 // which includes the following license/copyright statement:
       
    32 //
       
    33 // License for Berkeley SoftFloat Release 3e
       
    34 //
       
    35 // John R. Hauser
       
    36 // 2018 January 20
       
    37 //
       
    38 // The following applies to the whole of SoftFloat Release 3e as well as to
       
    39 // each source file individually.
       
    40 //
       
    41 // Copyright 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018 The Regents of the
       
    42 // University of California.  All rights reserved.
       
    43 //
       
    44 // Redistribution and use in source and binary forms, with or without
       
    45 // modification, are permitted provided that the following conditions are met:
       
    46 //
       
    47 //  1. Redistributions of source code must retain the above copyright notice,
       
    48 //     this list of conditions, and the following disclaimer.
       
    49 //
       
    50 //  2. Redistributions in binary form must reproduce the above copyright
       
    51 //     notice, this list of conditions, and the following disclaimer in the
       
    52 //     documentation and/or other materials provided with the distribution.
       
    53 //
       
    54 //  3. Neither the name of the University nor the names of its contributors
       
    55 //     may be used to endorse or promote products derived from this software
       
    56 //     without specific prior written permission.
       
    57 //
       
    58 // THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS "AS IS", AND ANY
       
    59 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
       
    60 // WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, ARE
       
    61 // DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE FOR ANY
       
    62 // DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
       
    63 // (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
       
    64 // LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
       
    65 // ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
       
    66 // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
       
    67 // THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
       
    68 //
       
    69 
       
    70 /*----------------------------------------------------------------------------
       
    71 | Software floating-point underflow tininess-detection mode.
       
    72 *----------------------------------------------------------------------------*/
       
    73 enum {
       
    74     softfloat_tininess_beforeRounding = 0,
       
    75     softfloat_tininess_afterRounding  = 1
       
    76 };
       
    77 
       
    78 /*----------------------------------------------------------------------------
       
    79 | Software floating-point exception flags.
       
    80 *----------------------------------------------------------------------------*/
       
    81 enum {
       
    82     softfloat_flag_inexact   =  1,
       
    83     softfloat_flag_underflow =  2,
       
    84     softfloat_flag_overflow  =  4,
       
    85     softfloat_flag_infinite  =  8,
       
    86     softfloat_flag_invalid   = 16
       
    87 };
       
    88 /*----------------------------------------------------------------------------
       
    89 | Software floating-point rounding mode.  (Mode "odd" is supported only if
       
    90 | SoftFloat is compiled with macro 'SOFTFLOAT_ROUND_ODD' defined.)
       
    91 *----------------------------------------------------------------------------*/
       
    92 enum {
       
    93     softfloat_round_near_even   = 0,
       
    94     softfloat_round_minMag      = 1,
       
    95     softfloat_round_min         = 2,
       
    96     softfloat_round_max         = 3,
       
    97     softfloat_round_near_maxMag = 4,
       
    98     softfloat_round_odd         = 6
       
    99 };
       
   100 #define init_detectTininess softfloat_tininess_afterRounding
       
   101 
       
   102 typedef unsigned char bool;
       
   103 typedef double        float64_t;
       
   104 
       
   105 #if defined(LITTLEENDIAN)
       
   106 struct uint128          { uint64_t v0, v64; };
       
   107 struct uint64_extra     { uint64_t extra, v; };
       
   108 struct uint128_extra    { uint64_t extra; struct uint128 v; };
       
   109 struct extFloat80M      { uint64_t signif; uint16_t signExp; };
       
   110 #else
       
   111 struct uint128          { uint64_t v64, v0; };
       
   112 struct uint64_extra     { uint64_t v, extra; };
       
   113 struct uint128_extra    { struct uint128 v; uint64_t extra; };
       
   114 struct extFloat80M      { uint16_t signExp; uint64_t signif; };
       
   115 #endif
       
   116 
       
   117 union ui64_f64          { uint64_t ui; float64_t f; };
       
   118 union ui128_f128        { struct uint128 ui; float128_t f; };
       
   119 struct exp32_sig128     { int_fast32_t exp; struct uint128 sig; };
       
   120 
       
   121 typedef struct extFloat80M        extFloat80_t;
       
   122 
       
   123 /*----------------------------------------------------------------------------
       
   124 | The bit pattern for a default generated 128-bit floating-point NaN.
       
   125 *----------------------------------------------------------------------------*/
       
   126 #define defaultNaNF128UI96_32 0xFFFF8000
       
   127 #define defaultNaNF128UI64_32 0
       
   128 #define defaultNaNF128UI32_32 0
       
   129 #define defaultNaNF128UI0_32  0
       
   130 
       
   131 // #define defaultNaNF128UI64_64 UINT64_C( 0x7FFF800000000000 )
       
   132 // #define defaultNaNF128UI0_64  UINT64_C( 0 )
       
   133 
       
   134 /*----------------------------------------------------------------------------
       
   135 | The bit pattern for a default generated 128-bit floating-point +INF.
       
   136 *----------------------------------------------------------------------------*/
       
   137 #define defaultPInfF128UI96_32 0x7FFF0000
       
   138 #define defaultPInfF128UI64_32 0
       
   139 #define defaultPInfF128UI32_32 0
       
   140 #define defaultPInfF128UI0_32  0
       
   141 /*----------------------------------------------------------------------------
       
   142 | The bit pattern for a default generated 128-bit floating-point -INF.
       
   143 *----------------------------------------------------------------------------*/
       
   144 #define defaultNInfF128UI96_32 0xFFFF0000
       
   145 #define defaultNInfF128UI64_32 0
       
   146 #define defaultNInfF128UI32_32 0
       
   147 #define defaultNInfF128UI0_32  0
       
   148 
       
   149 struct commonNaN {
       
   150     bool sign;
       
   151 #ifdef LITTLEENDIAN
       
   152     uint64_t v0, v64;
       
   153 #else
       
   154     uint64_t v64, v0;
       
   155 #endif
       
   156 };
       
   157 struct exp16_sig64 { int_fast16_t exp; uint_fast64_t sig; };
       
   158 
       
   159 #endif // SUPPORT_QUADFLOAT
    96 
   160 
    97 %}
   161 %}
    98 ! !
   162 ! !
    99 
   163 
   100 !QuadFloat primitiveVariables!
   164 !QuadFloat primitiveVariables!
   101 %{
   165 %{
       
   166 #ifdef SUPPORT_QUADFLOAT
       
   167 
       
   168 uint_fast8_t softfloat_exceptionFlags;
       
   169 uint_fast8_t softfloat_roundingMode = softfloat_round_near_even;
       
   170 uint_fast8_t softfloat_detectTininess = init_detectTininess;
       
   171 uint_fast8_t softfloat_exceptionFlags = 0;
       
   172 
       
   173 const uint_least8_t softfloat_countLeadingZeros8[256] = {
       
   174     8, 7, 6, 6, 5, 5, 5, 5, 4, 4, 4, 4, 4, 4, 4, 4,
       
   175     3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
       
   176     2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
       
   177     2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
       
   178     1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       
   179     1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       
   180     1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       
   181     1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       
   182     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       
   183     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       
   184     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       
   185     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       
   186     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       
   187     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       
   188     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       
   189     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
       
   190 };
       
   191 #endif // SUPPORT_QUADFLOAT
       
   192 
   102 %}
   193 %}
   103 ! !
   194 ! !
   104 
   195 
   105 !QuadFloat primitiveFunctions!
   196 !QuadFloat primitiveFunctions!
   106 %{
   197 %{
       
   198 #ifdef SUPPORT_QUADFLOAT
       
   199 
       
   200 #if defined(LITTLEENDIAN) || defined(__LSBFIRST__)
       
   201 #define wordIncr 1
       
   202 #define indexWord( total, n ) (n)
       
   203 #define indexWordHi( total ) ((total) - 1)
       
   204 #define indexWordLo( total ) 0
       
   205 #define indexMultiword( total, m, n ) (n)
       
   206 #define indexMultiwordHi( total, n ) ((total) - (n))
       
   207 #define indexMultiwordLo( total, n ) 0
       
   208 #define indexMultiwordHiBut( total, n ) (n)
       
   209 #define indexMultiwordLoBut( total, n ) 0
       
   210 #define INIT_UINTM4( v3, v2, v1, v0 ) { v0, v1, v2, v3 }
       
   211 #else
       
   212 #define wordIncr -1
       
   213 #define indexWord( total, n ) ((total) - 1 - (n))
       
   214 #define indexWordHi( total ) 0
       
   215 #define indexWordLo( total ) ((total) - 1)
       
   216 #define indexMultiword( total, m, n ) ((total) - 1 - (m))
       
   217 #define indexMultiwordHi( total, n ) 0
       
   218 #define indexMultiwordLo( total, n ) ((total) - (n))
       
   219 #define indexMultiwordHiBut( total, n ) 0
       
   220 #define indexMultiwordLoBut( total, n ) (n)
       
   221 #define INIT_UINTM4( v3, v2, v1, v0 ) { v3, v2, v1, v0 }
       
   222 #endif
       
   223 
       
   224 #define signF64UI( a ) ((bool) ((uint64_t) (a)>>63))
       
   225 #define expF64UI( a ) ((int_fast16_t) ((a)>>52) & 0x7FF)
       
   226 #define fracF64UI( a ) ((a) & UINT64_C( 0x000FFFFFFFFFFFFF ))
       
   227 #define softfloat_f64UIToCommonNaN( uiA, zPtr ) if ( ! ((uiA) & UINT64_C( 0x0008000000000000 )) ) softfloat_raiseFlags( softfloat_flag_invalid )
       
   228 
       
   229 #define signF128UI64( a64 ) ((bool) ((uint64_t) (a64)>>63))
       
   230 #define expF128UI64( a64 ) ((int_fast32_t) ((a64)>>48) & 0x7FFF)
       
   231 #define fracF128UI64( a64 ) ((a64) & UINT64_C( 0x0000FFFFFFFFFFFF ))
       
   232 #define packToF128UI64( sign, exp, sig64 ) (((uint_fast64_t) (sign)<<63) + ((uint_fast64_t) (exp)<<48) + (sig64))
       
   233 
       
   234 #define signF128UI96( a96 ) ((bool) ((uint32_t) (a96)>>31))
       
   235 #define expF128UI96( a96 ) ((int32_t) ((a96)>>16) & 0x7FFF)
       
   236 #define fracF128UI96( a96 ) ((a96) & 0x0000FFFF)
       
   237 #define packToF128UI96( sign, exp, sig96 ) (((uint32_t) (sign)<<31) + ((uint32_t) (exp)<<16) + (sig96))
       
   238 
       
   239 #define isNaNF128UI( a64, a0 ) (((~(a64) & UINT64_C( 0x7FFF000000000000 )) == 0) && (a0 || ((a64) & UINT64_C( 0x0000FFFFFFFFFFFF ))))
       
   240 #define softfloat_isSigNaNF128UI( uiA64, uiA0 ) ((((uiA64) & UINT64_C( 0x7FFF800000000000 )) == UINT64_C( 0x7FFF000000000000 )) && ((uiA0) || ((uiA64) & UINT64_C( 0x00007FFFFFFFFFFF ))))
       
   241 #define softfloat_approxRecip32_1( a ) ((uint32_t) (UINT64_C( 0x7FFFFFFFFFFFFFFF ) / (uint32_t) (a)))
       
   242 #define isInfF128UI( a64, a0 ) (((a64) & UINT64_C(0x7FFF800000000000)) == UINT64_C(0x7FFF000000000000))
       
   243 
       
   244 #define signExtF80UI64( a64 ) ((bool) ((uint16_t) (a64)>>15))
       
   245 #define expExtF80UI64( a64 ) ((a64) & 0x7FFF)
       
   246 #define packToExtF80UI64( sign, exp ) ((uint_fast16_t) (sign)<<15 | (exp))
       
   247 
       
   248 #define isNaNExtF80UI( a64, a0 ) ((((a64) & 0x7FFF) == 0x7FFF) && ((a0) & UINT64_C( 0x7FFFFFFFFFFFFFFF )))
       
   249 
       
   250 /*----------------------------------------------------------------------------
       
   251 | This function or macro is the same as 'softfloat_shortShiftLeftM' with
       
   252 | 'size_words' = 5 (N = 160).
       
   253 *----------------------------------------------------------------------------*/
       
   254 #define softfloat_shortShiftLeft160M( aPtr, dist, zPtr ) softfloat_shortShiftLeftM( 5, aPtr, dist, zPtr )
       
   255 /*----------------------------------------------------------------------------
       
   256 | This function or macro is the same as 'softfloat_shiftRightJamM' with
       
   257 | 'size_words' = 5 (N = 160).
       
   258 *----------------------------------------------------------------------------*/
       
   259 #define softfloat_shiftRightJam160M( aPtr, dist, zPtr ) softfloat_shiftRightJamM( 5, aPtr, dist, zPtr )
       
   260 /*----------------------------------------------------------------------------
       
   261 | This function or macro is the same as 'softfloat_shortShiftLeftM' with
       
   262 | 'size_words' = 4 (N = 128).
       
   263 *----------------------------------------------------------------------------*/
       
   264 #define softfloat_shortShiftLeft128M( aPtr, dist, zPtr ) softfloat_shortShiftLeftM( 4, aPtr, dist, zPtr )
       
   265 /*----------------------------------------------------------------------------
       
   266 | This function or macro is the same as 'softfloat_shortShiftLeftM' with
       
   267 | 'size_words' = 3 (N = 96).
       
   268 *----------------------------------------------------------------------------*/
       
   269 #define softfloat_shortShiftLeft96M( aPtr, dist, zPtr ) softfloat_shortShiftLeftM( 3, aPtr, dist, zPtr )
       
   270 
       
   271 static inline void
       
   272 softfloat_raiseFlags( uint_fast8_t flags ) {
       
   273     softfloat_exceptionFlags |= flags;
       
   274 }
       
   275 
       
   276 #if 1
       
   277 static inline void
       
   278 softfloat_commonNaNToF128M( uint32_t *zWPtr )
       
   279 {
       
   280     zWPtr[indexWord( 4, 3 )] = defaultNaNF128UI96_32;
       
   281     zWPtr[indexWord( 4, 2 )] = defaultNaNF128UI64_32;
       
   282     zWPtr[indexWord( 4, 1 )] = defaultNaNF128UI32_32;
       
   283     zWPtr[indexWord( 4, 0 )] = defaultNaNF128UI0_32;
       
   284 }
       
   285 #else
       
   286 # define softfloat_commonNaNToF128M( zWPtr ) \
       
   287 { \
       
   288     (zWPtr)[indexWord( 4, 3 )] = defaultNaNF128UI96_32; \
       
   289     (zWPtr)[indexWord( 4, 2 )] = defaultNaNF128UI64_32; \
       
   290     (zWPtr)[indexWord( 4, 1 )] = defaultNaNF128UI32_32; \
       
   291     (zWPtr)[indexWord( 4, 0 )] = defaultNaNF128UI0_32; \
       
   292 }
       
   293 #endif
       
   294 
       
   295 /*----------------------------------------------------------------------------
       
   296 | Assuming the unsigned integer formed from concatenating 'uiA64' and 'uiA0'
       
   297 | has the bit pattern of an 80-bit extended floating-point NaN, converts
       
   298 | this NaN to the common NaN form, and stores the resulting common NaN at the
       
   299 | location pointed to by 'zPtr'.  If the NaN is a signaling NaN, the invalid
       
   300 | exception is raised.
       
   301 *----------------------------------------------------------------------------*/
       
   302 #define softfloat_extF80UIToCommonNaN( uiA64, uiA0, zPtr ) if ( ! ((uiA0) & UINT64_C( 0x4000000000000000 )) ) softfloat_raiseFlags( softfloat_flag_invalid )
       
   303 
       
   304 static inline void
       
   305 softfloat_invalidF128M( uint32_t *zWPtr ) {
       
   306     softfloat_raiseFlags( softfloat_flag_invalid );
       
   307     softfloat_commonNaNToF128M( zWPtr );
       
   308 }
       
   309 
       
   310 static bool
       
   311 f128M_isSignalingNaN( const float128_t *aPtr )
       
   312 {
       
   313     const uint32_t *aWPtr;
       
   314     uint32_t uiA96;
       
   315 
       
   316     aWPtr = (const uint32_t *) aPtr;
       
   317     uiA96 = aWPtr[indexWordHi( 4 )];
       
   318     if ( (uiA96 & 0x7FFF8000) != 0x7FFF0000 ) return 0;
       
   319     return
       
   320 	((uiA96 & 0x00007FFF) != 0)
       
   321 	    || ((aWPtr[indexWord( 4, 2 )] | aWPtr[indexWord( 4, 1 )]
       
   322 		     | aWPtr[indexWord( 4, 0 )])
       
   323 		    != 0);
       
   324 }
       
   325 
       
   326 static void
       
   327 softfloat_shortShiftRightJamM(
       
   328      uint_fast8_t size_words,
       
   329      const uint32_t *aPtr,
       
   330      uint_fast8_t dist,
       
   331      uint32_t *zPtr
       
   332  )
       
   333 {
       
   334     uint_fast8_t uNegDist;
       
   335     unsigned int index, lastIndex;
       
   336     uint32_t partWordZ, wordA;
       
   337 
       
   338     uNegDist = -dist;
       
   339     index = indexWordLo( size_words );
       
   340     lastIndex = indexWordHi( size_words );
       
   341     wordA = aPtr[index];
       
   342     partWordZ = wordA>>dist;
       
   343     if ( partWordZ<<dist != wordA ) partWordZ |= 1;
       
   344     while ( index != lastIndex ) {
       
   345 	wordA = aPtr[index + wordIncr];
       
   346 	zPtr[index] = wordA<<(uNegDist & 31) | partWordZ;
       
   347 	index += wordIncr;
       
   348 	partWordZ = wordA>>dist;
       
   349     }
       
   350     zPtr[index] = partWordZ;
       
   351 
       
   352 }
       
   353 
       
   354 void
       
   355 softfloat_shiftRightJamM(
       
   356      uint_fast8_t size_words,
       
   357      const uint32_t *aPtr,
       
   358      uint32_t dist,
       
   359      uint32_t *zPtr
       
   360  )
       
   361 {
       
   362     uint32_t wordJam, wordDist, *ptr;
       
   363     uint_fast8_t i, innerDist;
       
   364 
       
   365     wordJam = 0;
       
   366     wordDist = dist>>5;
       
   367     if ( wordDist ) {
       
   368 	if ( size_words < wordDist ) wordDist = size_words;
       
   369 	ptr = (uint32_t *) (aPtr + indexMultiwordLo( size_words, wordDist ));
       
   370 	i = wordDist;
       
   371 	do {
       
   372 	    wordJam = *ptr++;
       
   373 	    if ( wordJam ) break;
       
   374 	    --i;
       
   375 	} while ( i );
       
   376 	ptr = zPtr;
       
   377     }
       
   378     if ( wordDist < size_words ) {
       
   379 	aPtr += indexMultiwordHiBut( size_words, wordDist );
       
   380 	innerDist = dist & 31;
       
   381 	if ( innerDist ) {
       
   382 	    softfloat_shortShiftRightJamM(
       
   383 		size_words - wordDist,
       
   384 		aPtr,
       
   385 		innerDist,
       
   386 		(zPtr + indexMultiwordLoBut( size_words, wordDist ))
       
   387 	    );
       
   388 	    if ( ! wordDist ) goto wordJam;
       
   389 	} else {
       
   390 	    aPtr += indexWordLo( size_words - wordDist );
       
   391 	    ptr = zPtr + indexWordLo( size_words );
       
   392 	    for ( i = size_words - wordDist; i; --i ) {
       
   393 		*ptr = *aPtr;
       
   394 		aPtr += wordIncr;
       
   395 		ptr += wordIncr;
       
   396 	    }
       
   397 	}
       
   398 	ptr = zPtr + indexMultiwordHi( size_words, wordDist );
       
   399     }
       
   400     do {
       
   401 	*ptr++ = 0;
       
   402 	--wordDist;
       
   403     } while ( wordDist );
       
   404  wordJam:
       
   405     if ( wordJam ) zPtr[indexWordLo( size_words )] |= 1;
       
   406 
       
   407 }
       
   408 
       
   409 void
       
   410  softfloat_shortShiftLeftM(
       
   411      uint_fast8_t size_words,
       
   412      const uint32_t *aPtr,
       
   413      uint_fast8_t dist,
       
   414      uint32_t *zPtr
       
   415  )
       
   416 {
       
   417     uint_fast8_t uNegDist;
       
   418     unsigned int index, lastIndex;
       
   419     uint32_t partWordZ, wordA;
       
   420 
       
   421     uNegDist = -dist;
       
   422     index = indexWordHi( size_words );
       
   423     lastIndex = indexWordLo( size_words );
       
   424     partWordZ = aPtr[index]<<dist;
       
   425     while ( index != lastIndex ) {
       
   426 	wordA = aPtr[index - wordIncr];
       
   427 	zPtr[index] = partWordZ | wordA>>(uNegDist & 31);
       
   428 	index -= wordIncr;
       
   429 	partWordZ = wordA<<dist;
       
   430     }
       
   431     zPtr[index] = partWordZ;
       
   432 }
       
   433 
       
   434 static void
       
   435 softfloat_shiftLeftM(
       
   436      uint_fast8_t size_words,
       
   437      const uint32_t *aPtr,
       
   438      uint32_t dist,
       
   439      uint32_t *zPtr
       
   440  )
       
   441 {
       
   442     uint32_t wordDist;
       
   443     uint_fast8_t innerDist;
       
   444     uint32_t *destPtr;
       
   445     uint_fast8_t i;
       
   446 
       
   447     wordDist = dist>>5;
       
   448     if ( wordDist < size_words ) {
       
   449 	aPtr += indexMultiwordLoBut( size_words, wordDist );
       
   450 	innerDist = dist & 31;
       
   451 	if ( innerDist ) {
       
   452 	    softfloat_shortShiftLeftM(
       
   453 		size_words - wordDist,
       
   454 		aPtr,
       
   455 		innerDist,
       
   456 		zPtr + indexMultiwordHiBut( size_words, wordDist )
       
   457 	    );
       
   458 	    if ( ! wordDist ) return;
       
   459 	} else {
       
   460 	    aPtr += indexWordHi( size_words - wordDist );
       
   461 	    destPtr = zPtr + indexWordHi( size_words );
       
   462 	    for ( i = size_words - wordDist; i; --i ) {
       
   463 		*destPtr = *aPtr;
       
   464 		aPtr -= wordIncr;
       
   465 		destPtr -= wordIncr;
       
   466 	    }
       
   467 	}
       
   468 	zPtr += indexMultiwordLo( size_words, wordDist );
       
   469     } else {
       
   470 	wordDist = size_words;
       
   471     }
       
   472     do {
       
   473 	*zPtr++ = 0;
       
   474 	--wordDist;
       
   475     } while ( wordDist );
       
   476 }
       
   477 
       
   478 /*----------------------------------------------------------------------------
       
   479 | This function or macro is the same as 'softfloat_shiftLeftM' with
       
   480 | 'size_words' = 4 (N = 128).
       
   481 *----------------------------------------------------------------------------*/
       
   482 #define softfloat_shiftLeft128M( aPtr, dist, zPtr ) softfloat_shiftLeftM( 4, aPtr, dist, zPtr )
       
   483 
       
   484 bool
       
   485 softfloat_isNaNF128M( const uint32_t *aWPtr )
       
   486 {
       
   487     uint32_t uiA96;
       
   488 
       
   489     uiA96 = aWPtr[indexWordHi( 4 )];
       
   490     if ( (~uiA96 & 0x7FFF0000) != 0 ) return 0;
       
   491     return
       
   492 	((uiA96 & 0x0000FFFF) != 0)
       
   493 	    || ((aWPtr[indexWord( 4, 2 )] | aWPtr[indexWord( 4, 1 )]
       
   494 		     | aWPtr[indexWord( 4, 0 )])
       
   495 		    != 0);
       
   496 }
       
   497 
       
   498 static inline uint_fast8_t
       
   499 softfloat_countLeadingZeros32( uint32_t a ) {
       
   500     return a ? __builtin_clz( a ) : 32;
       
   501 }
       
   502 
       
   503 static inline bool
       
   504 softfloat_eq128( uint64_t a64, uint64_t a0, uint64_t b64, uint64_t b0 ){
       
   505     return (a64 == b64) && (a0 == b0);
       
   506 }
       
   507 
       
   508 static inline bool
       
   509 softfloat_lt128( uint64_t a64, uint64_t a0, uint64_t b64, uint64_t b0 ){
       
   510     return (a64 < b64) || ((a64 == b64) && (a0 < b0));
       
   511 }
       
   512 
       
   513 static inline bool
       
   514 softfloat_le128( uint64_t a64, uint64_t a0, uint64_t b64, uint64_t b0 ){
       
   515     return (a64 < b64) || ((a64 == b64) && (a0 <= b0));
       
   516 }
       
   517 
       
   518 void
       
   519 softfloat_propagateNaNF128M(
       
   520      const uint32_t *aWPtr, const uint32_t *bWPtr, uint32_t *zWPtr )
       
   521 {
       
   522     bool isSigNaNA;
       
   523     const uint32_t *ptr;
       
   524     bool isSigNaNB;
       
   525     uint32_t uiA96, uiB96, wordMagA, wordMagB;
       
   526 
       
   527     isSigNaNA = f128M_isSignalingNaN( (const float128_t *) aWPtr );
       
   528     ptr = aWPtr;
       
   529     if ( ! bWPtr ) {
       
   530 	if ( isSigNaNA ) softfloat_raiseFlags( softfloat_flag_invalid );
       
   531 	goto copy;
       
   532     }
       
   533     isSigNaNB = f128M_isSignalingNaN( (const float128_t *) bWPtr );
       
   534     if ( isSigNaNA | isSigNaNB ) {
       
   535 	softfloat_raiseFlags( softfloat_flag_invalid );
       
   536 	if ( isSigNaNA ) {
       
   537 	    if ( isSigNaNB ) goto returnLargerUIMag;
       
   538 	    if ( softfloat_isNaNF128M( bWPtr ) ) goto copyB;
       
   539 	    goto copy;
       
   540 	} else {
       
   541 	    if ( softfloat_isNaNF128M( aWPtr ) ) goto copy;
       
   542 	    goto copyB;
       
   543 	}
       
   544     }
       
   545  returnLargerUIMag:
       
   546     uiA96 = aWPtr[indexWordHi( 4 )];
       
   547     uiB96 = bWPtr[indexWordHi( 4 )];
       
   548     wordMagA = uiA96 & 0x7FFFFFFF;
       
   549     wordMagB = uiB96 & 0x7FFFFFFF;
       
   550     if ( wordMagA < wordMagB ) goto copyB;
       
   551     if ( wordMagB < wordMagA ) goto copy;
       
   552     wordMagA = aWPtr[indexWord( 4, 2 )];
       
   553     wordMagB = bWPtr[indexWord( 4, 2 )];
       
   554     if ( wordMagA < wordMagB ) goto copyB;
       
   555     if ( wordMagB < wordMagA ) goto copy;
       
   556     wordMagA = aWPtr[indexWord( 4, 1 )];
       
   557     wordMagB = bWPtr[indexWord( 4, 1 )];
       
   558     if ( wordMagA < wordMagB ) goto copyB;
       
   559     if ( wordMagB < wordMagA ) goto copy;
       
   560     wordMagA = aWPtr[indexWord( 4, 0 )];
       
   561     wordMagB = bWPtr[indexWord( 4, 0 )];
       
   562     if ( wordMagA < wordMagB ) goto copyB;
       
   563     if ( wordMagB < wordMagA ) goto copy;
       
   564     if ( uiA96 < uiB96 ) goto copy;
       
   565  copyB:
       
   566     ptr = bWPtr;
       
   567  copy:
       
   568     zWPtr[indexWordHi( 4 )] = ptr[indexWordHi( 4 )] | 0x00008000;
       
   569     zWPtr[indexWord( 4, 2 )] = ptr[indexWord( 4, 2 )];
       
   570     zWPtr[indexWord( 4, 1 )] = ptr[indexWord( 4, 1 )];
       
   571     zWPtr[indexWord( 4, 0 )] = ptr[indexWord( 4, 0 )];
       
   572 }
       
   573 
       
   574 static inline uint_fast8_t
       
   575 softfloat_countLeadingZeros64( uint64_t a )
       
   576 {
       
   577     uint_fast8_t count;
       
   578     uint32_t a32;
       
   579 
       
   580     count = 0;
       
   581     a32 = a>>32;
       
   582     if ( ! a32 ) {
       
   583 	count = 32;
       
   584 	a32 = a;
       
   585     }
       
   586     /*------------------------------------------------------------------------
       
   587     | From here, result is current count + count leading zeros of `a32'.
       
   588     *------------------------------------------------------------------------*/
       
   589     if ( a32 < 0x10000 ) {
       
   590 	count += 16;
       
   591 	a32 <<= 16;
       
   592     }
       
   593     if ( a32 < 0x1000000 ) {
       
   594 	count += 8;
       
   595 	a32 <<= 8;
       
   596     }
       
   597     count += softfloat_countLeadingZeros8[a32>>24];
       
   598     return count;
       
   599 }
       
   600 
       
   601 static inline struct exp16_sig64
       
   602 softfloat_normSubnormalF64Sig( uint_fast64_t sig )
       
   603 {
       
   604     int_fast8_t shiftDist;
       
   605     struct exp16_sig64 z;
       
   606 
       
   607     shiftDist = softfloat_countLeadingZeros64( sig ) - 11;
       
   608     z.exp = 1 - shiftDist;
       
   609     z.sig = sig<<shiftDist;
       
   610     return z;
       
   611 }
       
   612 
       
   613 static inline struct uint128_extra
       
   614 softfloat_shiftRightJam128Extra(
       
   615      uint64_t a64, uint64_t a0, uint64_t extra, uint_fast32_t dist )
       
   616 {
       
   617     uint_fast8_t u8NegDist;
       
   618     struct uint128_extra z;
       
   619 
       
   620     u8NegDist = -dist;
       
   621     if ( dist < 64 ) {
       
   622 	z.v.v64 = a64>>dist;
       
   623 	z.v.v0 = a64<<(u8NegDist & 63) | a0>>dist;
       
   624 	z.extra = a0<<(u8NegDist & 63);
       
   625     } else {
       
   626 	z.v.v64 = 0;
       
   627 	if ( dist == 64 ) {
       
   628 	    z.v.v0 = a64;
       
   629 	    z.extra = a0;
       
   630 	} else {
       
   631 	    extra |= a0;
       
   632 	    if ( dist < 128 ) {
       
   633 		z.v.v0 = a64>>(dist & 63);
       
   634 		z.extra = a64<<(u8NegDist & 63);
       
   635 	    } else {
       
   636 		z.v.v0 = 0;
       
   637 		z.extra = (dist == 128) ? a64 : (a64 != 0);
       
   638 	    }
       
   639 	}
       
   640     }
       
   641     z.extra |= (extra != 0);
       
   642     return z;
       
   643 }
       
   644 
       
   645 static inline struct uint128_extra
       
   646 softfloat_shortShiftRightJam128Extra(
       
   647      uint64_t a64, uint64_t a0, uint64_t extra, uint_fast8_t dist )
       
   648 {
       
   649     uint_fast8_t negDist = -dist;
       
   650     struct uint128_extra z;
       
   651     z.v.v64 = a64>>dist;
       
   652     z.v.v0 = a64<<(negDist & 63) | a0>>dist;
       
   653     z.extra = a0<<(negDist & 63) | (extra != 0);
       
   654     return z;
       
   655 }
       
   656 
       
   657 static inline struct uint128
       
   658 softfloat_shiftRightJam128( uint64_t a64, uint64_t a0, uint_fast32_t dist )
       
   659 {
       
   660     uint_fast8_t u8NegDist;
       
   661     struct uint128 z;
       
   662 
       
   663     if ( dist < 64 ) {
       
   664 	u8NegDist = -dist;
       
   665 	z.v64 = a64>>dist;
       
   666 	z.v0 =
       
   667 	    a64<<(u8NegDist & 63) | a0>>dist
       
   668 		| ((uint64_t) (a0<<(u8NegDist & 63)) != 0);
       
   669     } else {
       
   670 	z.v64 = 0;
       
   671 	z.v0 =
       
   672 	    (dist < 127)
       
   673 		? a64>>(dist & 63)
       
   674 		      | (((a64 & (((uint_fast64_t) 1<<(dist & 63)) - 1)) | a0)
       
   675 			     != 0)
       
   676 		: ((a64 | a0) != 0);
       
   677     }
       
   678     return z;
       
   679 }
       
   680 
       
   681 static inline struct uint128
       
   682 softfloat_shortShiftLeft128( uint64_t a64, uint64_t a0, uint_fast8_t dist )
       
   683 {
       
   684     struct uint128 z;
       
   685     z.v64 = a64<<dist | a0>>(-dist & 63);
       
   686     z.v0 = a0<<dist;
       
   687     return z;
       
   688 }
       
   689 
       
   690 static int_fast8_t
       
   691 softfloat_compare128M( const uint32_t *aPtr, const uint32_t *bPtr )
       
   692 {
       
   693     unsigned int index, lastIndex;
       
   694     uint32_t wordA, wordB;
       
   695 
       
   696     index = indexWordHi( 4 );
       
   697     lastIndex = indexWordLo( 4 );
       
   698     for (;;) {
       
   699 	wordA = aPtr[index];
       
   700 	wordB = bPtr[index];
       
   701 	if ( wordA != wordB ) return (wordA < wordB) ? -1 : 1;
       
   702 	if ( index == lastIndex ) break;
       
   703 	index -= wordIncr;
       
   704     }
       
   705     return 0;
       
   706 }
       
   707 
       
   708 static void
       
   709 softfloat_roundPackMToF128M(
       
   710      bool sign, int32_t exp, uint32_t *extSigPtr, uint32_t *zWPtr )
       
   711 {
       
   712     uint_fast8_t roundingMode;
       
   713     bool roundNearEven;
       
   714     uint32_t sigExtra;
       
   715     bool doIncrement, isTiny;
       
   716     static const uint32_t maxSig[4] =
       
   717 	INIT_UINTM4( 0x0001FFFF, 0xFFFFFFFF, 0xFFFFFFFF, 0xFFFFFFFF );
       
   718     uint32_t ui, uj;
       
   719 
       
   720     /*------------------------------------------------------------------------
       
   721     *------------------------------------------------------------------------*/
       
   722     roundingMode = softfloat_roundingMode;
       
   723     roundNearEven = (roundingMode == softfloat_round_near_even);
       
   724     sigExtra = extSigPtr[indexWordLo( 5 )];
       
   725     doIncrement = (0x80000000 <= sigExtra);
       
   726     if ( ! roundNearEven && (roundingMode != softfloat_round_near_maxMag) ) {
       
   727 	doIncrement =
       
   728 	    (roundingMode
       
   729 		 == (sign ? softfloat_round_min : softfloat_round_max))
       
   730 		&& sigExtra;
       
   731     }
       
   732     /*------------------------------------------------------------------------
       
   733     *------------------------------------------------------------------------*/
       
   734     if ( 0x7FFD <= (uint32_t) exp ) {
       
   735 	if ( exp < 0 ) {
       
   736 	    /*----------------------------------------------------------------
       
   737 	    *----------------------------------------------------------------*/
       
   738 	    isTiny =
       
   739 		   (softfloat_detectTininess
       
   740 			== softfloat_tininess_beforeRounding)
       
   741 		|| (exp < -1)
       
   742 		|| ! doIncrement
       
   743 		|| (softfloat_compare128M(
       
   744 			extSigPtr + indexMultiwordHi( 5, 4 ), maxSig )
       
   745 			< 0);
       
   746 	    softfloat_shiftRightJam160M( extSigPtr, -exp, extSigPtr );
       
   747 	    exp = 0;
       
   748 	    sigExtra = extSigPtr[indexWordLo( 5 )];
       
   749 	    if ( isTiny && sigExtra ) {
       
   750 		softfloat_raiseFlags( softfloat_flag_underflow );
       
   751 	    }
       
   752 	    doIncrement = (0x80000000 <= sigExtra);
       
   753 	    if (
       
   754 		   ! roundNearEven
       
   755 		&& (roundingMode != softfloat_round_near_maxMag)
       
   756 	    ) {
       
   757 		doIncrement =
       
   758 		    (roundingMode
       
   759 			 == (sign ? softfloat_round_min : softfloat_round_max))
       
   760 			&& sigExtra;
       
   761 	    }
       
   762 	} else if (
       
   763 	       (0x7FFD < exp)
       
   764 	    || ((exp == 0x7FFD) && doIncrement
       
   765 		    && (softfloat_compare128M(
       
   766 			    extSigPtr + indexMultiwordHi( 5, 4 ), maxSig )
       
   767 			    == 0))
       
   768 	) {
       
   769 	    /*----------------------------------------------------------------
       
   770 	    *----------------------------------------------------------------*/
       
   771 	    softfloat_raiseFlags(
       
   772 		softfloat_flag_overflow | softfloat_flag_inexact );
       
   773 	    if (
       
   774 		   roundNearEven
       
   775 		|| (roundingMode == softfloat_round_near_maxMag)
       
   776 		|| (roundingMode
       
   777 			== (sign ? softfloat_round_min : softfloat_round_max))
       
   778 	    ) {
       
   779 		ui = packToF128UI96( sign, 0x7FFF, 0 );
       
   780 		uj = 0;
       
   781 	    } else {
       
   782 		ui = packToF128UI96( sign, 0x7FFE, 0x0000FFFF );
       
   783 		uj = 0xFFFFFFFF;
       
   784 	    }
       
   785 	    zWPtr[indexWordHi( 4 )] = ui;
       
   786 	    zWPtr[indexWord( 4, 2 )] = uj;
       
   787 	    zWPtr[indexWord( 4, 1 )] = uj;
       
   788 	    zWPtr[indexWord( 4, 0 )] = uj;
       
   789 	    return;
       
   790 	}
       
   791     }
       
   792     /*------------------------------------------------------------------------
       
   793     *------------------------------------------------------------------------*/
       
   794     uj = extSigPtr[indexWord( 5, 1 )];
       
   795     if ( sigExtra ) {
       
   796 	softfloat_exceptionFlags |= softfloat_flag_inexact;
       
   797 #ifdef SOFTFLOAT_ROUND_ODD
       
   798 	if ( roundingMode == softfloat_round_odd ) {
       
   799 	    uj |= 1;
       
   800 	    goto noIncrementPackReturn;
       
   801 	}
       
   802 #endif
       
   803     }
       
   804     if ( doIncrement ) {
       
   805 	++uj;
       
   806 	if ( uj ) {
       
   807 	    if ( ! (sigExtra & 0x7FFFFFFF) && roundNearEven ) uj &= ~1;
       
   808 	    zWPtr[indexWord( 4, 2 )] = extSigPtr[indexWord( 5, 3 )];
       
   809 	    zWPtr[indexWord( 4, 1 )] = extSigPtr[indexWord( 5, 2 )];
       
   810 	    zWPtr[indexWord( 4, 0 )] = uj;
       
   811 	    ui = extSigPtr[indexWordHi( 5 )];
       
   812 	} else {
       
   813 	    zWPtr[indexWord( 4, 0 )] = uj;
       
   814 	    ui = extSigPtr[indexWord( 5, 2 )] + 1;
       
   815 	    zWPtr[indexWord( 4, 1 )] = ui;
       
   816 	    uj = extSigPtr[indexWord( 5, 3 )];
       
   817 	    if ( ui ) {
       
   818 		zWPtr[indexWord( 4, 2 )] = uj;
       
   819 		ui = extSigPtr[indexWordHi( 5 )];
       
   820 	    } else {
       
   821 		++uj;
       
   822 		zWPtr[indexWord( 4, 2 )] = uj;
       
   823 		ui = extSigPtr[indexWordHi( 5 )];
       
   824 		if ( ! uj ) ++ui;
       
   825 	    }
       
   826 	}
       
   827     } else {
       
   828  noIncrementPackReturn:
       
   829 	zWPtr[indexWord( 4, 0 )] = uj;
       
   830 	ui = extSigPtr[indexWord( 5, 2 )];
       
   831 	zWPtr[indexWord( 4, 1 )] = ui;
       
   832 	uj |= ui;
       
   833 	ui = extSigPtr[indexWord( 5, 3 )];
       
   834 	zWPtr[indexWord( 4, 2 )] = ui;
       
   835 	uj |= ui;
       
   836 	ui = extSigPtr[indexWordHi( 5 )];
       
   837 	uj |= ui;
       
   838 	if ( ! uj ) exp = 0;
       
   839     }
       
   840     zWPtr[indexWordHi( 4 )] = packToF128UI96( sign, exp, ui );
       
   841 }
       
   842 
       
   843 static int
       
   844 softfloat_shiftNormSigF128M(
       
   845      const uint32_t *wPtr, uint_fast8_t shiftDist, uint32_t *sigPtr )
       
   846 {
       
   847     uint32_t wordSig;
       
   848     int32_t exp;
       
   849     uint32_t leadingBit;
       
   850 
       
   851     wordSig = wPtr[indexWordHi( 4 )];
       
   852     exp = expF128UI96( wordSig );
       
   853     if ( exp ) {
       
   854 	softfloat_shortShiftLeft128M( wPtr, shiftDist, sigPtr );
       
   855 	leadingBit = 0x00010000<<shiftDist;
       
   856 	sigPtr[indexWordHi( 4 )] =
       
   857 	    (sigPtr[indexWordHi( 4 )] & (leadingBit - 1)) | leadingBit;
       
   858     } else {
       
   859 	exp = 16;
       
   860 	wordSig &= 0x7FFFFFFF;
       
   861 	if ( ! wordSig ) {
       
   862 	    exp = -16;
       
   863 	    wordSig = wPtr[indexWord( 4, 2 )];
       
   864 	    if ( ! wordSig ) {
       
   865 		exp = -48;
       
   866 		wordSig = wPtr[indexWord( 4, 1 )];
       
   867 		if ( ! wordSig ) {
       
   868 		    wordSig = wPtr[indexWord( 4, 0 )];
       
   869 		    if ( ! wordSig ) return -128;
       
   870 		    exp = -80;
       
   871 		}
       
   872 	    }
       
   873 	}
       
   874 	exp -= softfloat_countLeadingZeros32( wordSig );
       
   875 	softfloat_shiftLeft128M( wPtr, 1 - exp + shiftDist, sigPtr );
       
   876     }
       
   877     return exp;
       
   878 }
       
   879 
       
   880 bool
       
   881 softfloat_tryPropagateNaNF128M(
       
   882      const uint32_t *aWPtr, const uint32_t *bWPtr, uint32_t *zWPtr )
       
   883 {
       
   884 
       
   885     if ( softfloat_isNaNF128M( aWPtr ) || softfloat_isNaNF128M( bWPtr ) ) {
       
   886 	softfloat_propagateNaNF128M( aWPtr, bWPtr, zWPtr );
       
   887 	return 1;
       
   888     }
       
   889     return 0;
       
   890 }
       
   891 
       
   892 static inline struct uint128
       
   893 softfloat_add128( uint64_t a64, uint64_t a0, uint64_t b64, uint64_t b0 )
       
   894 {
       
   895     struct uint128 z;
       
   896     z.v0 = a0 + b0;
       
   897     z.v64 = a64 + b64 + (z.v0 < a0);
       
   898     return z;
       
   899 }
       
   900 
       
   901 static inline struct uint128
       
   902 softfloat_sub128( uint64_t a64, uint64_t a0, uint64_t b64, uint64_t b0 )
       
   903 {
       
   904     struct uint128 z;
       
   905     z.v0 = a0 - b0;
       
   906     z.v64 = a64 - b64;
       
   907     z.v64 -= (a0 < b0);
       
   908     return z;
       
   909 }
       
   910 
       
   911 static inline struct uint128
       
   912 softfloat_mul128By32( uint64_t a64, uint64_t a0, uint32_t b )
       
   913 {
       
   914     union { unsigned __int128 ui; struct uint128 s; } uZ;
       
   915     uZ.ui = ((unsigned __int128) a64<<64 | a0) * b;
       
   916     return uZ.s;
       
   917 }
       
   918 
       
   919 static void
       
   920 softfloat_mul128MTo256M(
       
   921      const uint32_t *aPtr, const uint32_t *bPtr, uint32_t *zPtr )
       
   922 {
       
   923     uint32_t *lastZPtr, wordB;
       
   924     uint64_t dwordProd;
       
   925     uint32_t wordZ;
       
   926     uint_fast8_t carry;
       
   927 
       
   928     bPtr += indexWordLo( 4 );
       
   929     lastZPtr = zPtr + indexMultiwordHi( 8, 5 );
       
   930     zPtr += indexMultiwordLo( 8, 5 );
       
   931     wordB = *bPtr;
       
   932     dwordProd = (uint64_t) aPtr[indexWord( 4, 0 )] * wordB;
       
   933     zPtr[indexWord( 5, 0 )] = dwordProd;
       
   934     dwordProd = (uint64_t) aPtr[indexWord( 4, 1 )] * wordB + (dwordProd>>32);
       
   935     zPtr[indexWord( 5, 1 )] = dwordProd;
       
   936     dwordProd = (uint64_t) aPtr[indexWord( 4, 2 )] * wordB + (dwordProd>>32);
       
   937     zPtr[indexWord( 5, 2 )] = dwordProd;
       
   938     dwordProd = (uint64_t) aPtr[indexWord( 4, 3 )] * wordB + (dwordProd>>32);
       
   939     zPtr[indexWord( 5, 3 )] = dwordProd;
       
   940     zPtr[indexWord( 5, 4 )] = dwordProd>>32;
       
   941     do {
       
   942 	bPtr += wordIncr;
       
   943 	zPtr += wordIncr;
       
   944 	wordB = *bPtr;
       
   945 	dwordProd = (uint64_t) aPtr[indexWord( 4, 0 )] * wordB;
       
   946 	wordZ = zPtr[indexWord( 5, 0 )] + (uint32_t) dwordProd;
       
   947 	zPtr[indexWord( 5, 0 )] = wordZ;
       
   948 	carry = (wordZ < (uint32_t) dwordProd);
       
   949 	dwordProd =
       
   950 	    (uint64_t) aPtr[indexWord( 4, 1 )] * wordB + (dwordProd>>32);
       
   951 	wordZ = zPtr[indexWord( 5, 1 )] + (uint32_t) dwordProd + carry;
       
   952 	zPtr[indexWord( 5, 1 )] = wordZ;
       
   953 	if ( wordZ != (uint32_t) dwordProd ) {
       
   954 	    carry = (wordZ < (uint32_t) dwordProd);
       
   955 	}
       
   956 	dwordProd =
       
   957 	    (uint64_t) aPtr[indexWord( 4, 2 )] * wordB + (dwordProd>>32);
       
   958 	wordZ = zPtr[indexWord( 5, 2 )] + (uint32_t) dwordProd + carry;
       
   959 	zPtr[indexWord( 5, 2 )] = wordZ;
       
   960 	if ( wordZ != (uint32_t) dwordProd ) {
       
   961 	    carry = (wordZ < (uint32_t) dwordProd);
       
   962 	}
       
   963 	dwordProd =
       
   964 	    (uint64_t) aPtr[indexWord( 4, 3 )] * wordB + (dwordProd>>32);
       
   965 	wordZ = zPtr[indexWord( 5, 3 )] + (uint32_t) dwordProd + carry;
       
   966 	zPtr[indexWord( 5, 3 )] = wordZ;
       
   967 	if ( wordZ != (uint32_t) dwordProd ) {
       
   968 	    carry = (wordZ < (uint32_t) dwordProd);
       
   969 	}
       
   970 	zPtr[indexWord( 5, 4 )] = (dwordProd>>32) + carry;
       
   971     } while ( zPtr != lastZPtr );
       
   972 }
       
   973 
       
   974 float128_t
       
   975 softfloat_roundPackToF128(
       
   976      bool sign,
       
   977      int_fast32_t exp,
       
   978      uint_fast64_t sig64,
       
   979      uint_fast64_t sig0,
       
   980      uint_fast64_t sigExtra
       
   981  )
       
   982 {
       
   983     uint_fast8_t roundingMode;
       
   984     bool roundNearEven, doIncrement, isTiny;
       
   985     struct uint128_extra sig128Extra;
       
   986     uint_fast64_t uiZ64, uiZ0;
       
   987     struct uint128 sig128;
       
   988     union ui128_f128 uZ;
       
   989 
       
   990     /*------------------------------------------------------------------------
       
   991     *------------------------------------------------------------------------*/
       
   992     roundingMode = softfloat_roundingMode;
       
   993     roundNearEven = (roundingMode == softfloat_round_near_even);
       
   994     doIncrement = (UINT64_C( 0x8000000000000000 ) <= sigExtra);
       
   995     if ( ! roundNearEven && (roundingMode != softfloat_round_near_maxMag) ) {
       
   996 	doIncrement =
       
   997 	    (roundingMode
       
   998 		 == (sign ? softfloat_round_min : softfloat_round_max))
       
   999 		&& sigExtra;
       
  1000     }
       
  1001     /*------------------------------------------------------------------------
       
  1002     *------------------------------------------------------------------------*/
       
  1003     if ( 0x7FFD <= (uint32_t) exp ) {
       
  1004 	if ( exp < 0 ) {
       
  1005 	    /*----------------------------------------------------------------
       
  1006 	    *----------------------------------------------------------------*/
       
  1007 	    isTiny =
       
  1008 		   (softfloat_detectTininess
       
  1009 			== softfloat_tininess_beforeRounding)
       
  1010 		|| (exp < -1)
       
  1011 		|| ! doIncrement
       
  1012 		|| softfloat_lt128(
       
  1013 		       sig64,
       
  1014 		       sig0,
       
  1015 		       UINT64_C( 0x0001FFFFFFFFFFFF ),
       
  1016 		       UINT64_C( 0xFFFFFFFFFFFFFFFF )
       
  1017 		   );
       
  1018 	    sig128Extra =
       
  1019 		softfloat_shiftRightJam128Extra( sig64, sig0, sigExtra, -exp );
       
  1020 	    sig64 = sig128Extra.v.v64;
       
  1021 	    sig0  = sig128Extra.v.v0;
       
  1022 	    sigExtra = sig128Extra.extra;
       
  1023 	    exp = 0;
       
  1024 	    if ( isTiny && sigExtra ) {
       
  1025 		softfloat_raiseFlags( softfloat_flag_underflow );
       
  1026 	    }
       
  1027 	    doIncrement = (UINT64_C( 0x8000000000000000 ) <= sigExtra);
       
  1028 	    if (
       
  1029 		   ! roundNearEven
       
  1030 		&& (roundingMode != softfloat_round_near_maxMag)
       
  1031 	    ) {
       
  1032 		doIncrement =
       
  1033 		    (roundingMode
       
  1034 			 == (sign ? softfloat_round_min : softfloat_round_max))
       
  1035 			&& sigExtra;
       
  1036 	    }
       
  1037 	} else if (
       
  1038 	       (0x7FFD < exp)
       
  1039 	    || ((exp == 0x7FFD)
       
  1040 		    && softfloat_eq128(
       
  1041 			   sig64,
       
  1042 			   sig0,
       
  1043 			   UINT64_C( 0x0001FFFFFFFFFFFF ),
       
  1044 			   UINT64_C( 0xFFFFFFFFFFFFFFFF )
       
  1045 		       )
       
  1046 		    && doIncrement)
       
  1047 	) {
       
  1048 	    /*----------------------------------------------------------------
       
  1049 	    *----------------------------------------------------------------*/
       
  1050 	    softfloat_raiseFlags(
       
  1051 		softfloat_flag_overflow | softfloat_flag_inexact );
       
  1052 	    if (
       
  1053 		   roundNearEven
       
  1054 		|| (roundingMode == softfloat_round_near_maxMag)
       
  1055 		|| (roundingMode
       
  1056 			== (sign ? softfloat_round_min : softfloat_round_max))
       
  1057 	    ) {
       
  1058 		uiZ64 = packToF128UI64( sign, 0x7FFF, 0 );
       
  1059 		uiZ0  = 0;
       
  1060 	    } else {
       
  1061 		uiZ64 = packToF128UI64( sign, 0x7FFE, UINT64_C( 0x0000FFFFFFFFFFFF ) );
       
  1062 		uiZ0 = UINT64_C( 0xFFFFFFFFFFFFFFFF );
       
  1063 	    }
       
  1064 	    goto uiZ;
       
  1065 	}
       
  1066     }
       
  1067     /*------------------------------------------------------------------------
       
  1068     *------------------------------------------------------------------------*/
       
  1069     if ( sigExtra ) {
       
  1070 	softfloat_exceptionFlags |= softfloat_flag_inexact;
       
  1071 #ifdef SOFTFLOAT_ROUND_ODD
       
  1072 	if ( roundingMode == softfloat_round_odd ) {
       
  1073 	    sig0 |= 1;
       
  1074 	    goto packReturn;
       
  1075 	}
       
  1076 #endif
       
  1077     }
       
  1078     if ( doIncrement ) {
       
  1079 	sig128 = softfloat_add128( sig64, sig0, 0, 1 );
       
  1080 	sig64 = sig128.v64;
       
  1081 	sig0 =
       
  1082 	    sig128.v0
       
  1083 		& ~(uint64_t)
       
  1084 		       (! (sigExtra & UINT64_C( 0x7FFFFFFFFFFFFFFF ))
       
  1085 			    & roundNearEven);
       
  1086     } else {
       
  1087 	if ( ! (sig64 | sig0) ) exp = 0;
       
  1088     }
       
  1089     /*------------------------------------------------------------------------
       
  1090     *------------------------------------------------------------------------*/
       
  1091  packReturn:
       
  1092     uiZ64 = packToF128UI64( sign, exp, sig64 );
       
  1093     uiZ0  = sig0;
       
  1094  uiZ:
       
  1095     uZ.ui.v64 = uiZ64;
       
  1096     uZ.ui.v0  = uiZ0;
       
  1097     return uZ.f;
       
  1098 }
       
  1099 
       
  1100 float128_t
       
  1101 softfloat_normRoundPackToF128(
       
  1102      bool sign, int_fast32_t exp, uint_fast64_t sig64, uint_fast64_t sig0 )
       
  1103 {
       
  1104     int_fast8_t shiftDist;
       
  1105     struct uint128 sig128;
       
  1106     union ui128_f128 uZ;
       
  1107     uint_fast64_t sigExtra;
       
  1108     struct uint128_extra sig128Extra;
       
  1109 
       
  1110     if ( ! sig64 ) {
       
  1111 	exp -= 64;
       
  1112 	sig64 = sig0;
       
  1113 	sig0 = 0;
       
  1114     }
       
  1115     shiftDist = softfloat_countLeadingZeros64( sig64 ) - 15;
       
  1116     exp -= shiftDist;
       
  1117     if ( 0 <= shiftDist ) {
       
  1118 	if ( shiftDist ) {
       
  1119 	    sig128 = softfloat_shortShiftLeft128( sig64, sig0, shiftDist );
       
  1120 	    sig64 = sig128.v64;
       
  1121 	    sig0  = sig128.v0;
       
  1122 	}
       
  1123 	if ( (uint32_t) exp < 0x7FFD ) {
       
  1124 	    uZ.ui.v64 = packToF128UI64( sign, sig64 | sig0 ? exp : 0, sig64 );
       
  1125 	    uZ.ui.v0  = sig0;
       
  1126 	    return uZ.f;
       
  1127 	}
       
  1128 	sigExtra = 0;
       
  1129     } else {
       
  1130 	sig128Extra =
       
  1131 	    softfloat_shortShiftRightJam128Extra( sig64, sig0, 0, -shiftDist );
       
  1132 	sig64 = sig128Extra.v.v64;
       
  1133 	sig0  = sig128Extra.v.v0;
       
  1134 	sigExtra = sig128Extra.extra;
       
  1135     }
       
  1136     return softfloat_roundPackToF128( sign, exp, sig64, sig0, sigExtra );
       
  1137 }
       
  1138 
       
  1139 struct exp32_sig128
       
  1140 softfloat_normSubnormalF128Sig( uint_fast64_t sig64, uint_fast64_t sig0 )
       
  1141 {
       
  1142     int_fast8_t shiftDist;
       
  1143     struct exp32_sig128 z;
       
  1144 
       
  1145     if ( ! sig64 ) {
       
  1146 	shiftDist = softfloat_countLeadingZeros64( sig0 ) - 15;
       
  1147 	z.exp = -63 - shiftDist;
       
  1148 	if ( shiftDist < 0 ) {
       
  1149 	    z.sig.v64 = sig0>>-shiftDist;
       
  1150 	    z.sig.v0  = sig0<<(shiftDist & 63);
       
  1151 	} else {
       
  1152 	    z.sig.v64 = sig0<<shiftDist;
       
  1153 	    z.sig.v0  = 0;
       
  1154 	}
       
  1155     } else {
       
  1156 	shiftDist = softfloat_countLeadingZeros64( sig64 ) - 15;
       
  1157 	z.exp = 1 - shiftDist;
       
  1158 	z.sig = softfloat_shortShiftLeft128( sig64, sig0, shiftDist );
       
  1159     }
       
  1160     return z;
       
  1161 }
       
  1162 
       
  1163 struct uint128
       
  1164 softfloat_propagateNaNF128UI(
       
  1165      uint_fast64_t uiA64,
       
  1166      uint_fast64_t uiA0,
       
  1167      uint_fast64_t uiB64,
       
  1168      uint_fast64_t uiB0
       
  1169  )
       
  1170 {
       
  1171     bool isSigNaNA, isSigNaNB;
       
  1172     uint_fast64_t uiNonsigA64, uiNonsigB64, uiMagA64, uiMagB64;
       
  1173     struct uint128 uiZ;
       
  1174 
       
  1175     /*------------------------------------------------------------------------
       
  1176     *------------------------------------------------------------------------*/
       
  1177     isSigNaNA = softfloat_isSigNaNF128UI( uiA64, uiA0 );
       
  1178     isSigNaNB = softfloat_isSigNaNF128UI( uiB64, uiB0 );
       
  1179     /*------------------------------------------------------------------------
       
  1180     | Make NaNs non-signaling.
       
  1181     *------------------------------------------------------------------------*/
       
  1182     uiNonsigA64 = uiA64 | UINT64_C( 0x0000800000000000 );
       
  1183     uiNonsigB64 = uiB64 | UINT64_C( 0x0000800000000000 );
       
  1184     /*------------------------------------------------------------------------
       
  1185     *------------------------------------------------------------------------*/
       
  1186     if ( isSigNaNA | isSigNaNB ) {
       
  1187 	softfloat_raiseFlags( softfloat_flag_invalid );
       
  1188 	if ( isSigNaNA ) {
       
  1189 	    if ( isSigNaNB ) goto returnLargerMag;
       
  1190 	    if ( isNaNF128UI( uiB64, uiB0 ) ) goto returnB;
       
  1191 	    goto returnA;
       
  1192 	} else {
       
  1193 	    if ( isNaNF128UI( uiA64, uiA0 ) ) goto returnA;
       
  1194 	    goto returnB;
       
  1195 	}
       
  1196     }
       
  1197  returnLargerMag:
       
  1198     uiMagA64 = uiA64 & UINT64_C( 0x7FFFFFFFFFFFFFFF );
       
  1199     uiMagB64 = uiB64 & UINT64_C( 0x7FFFFFFFFFFFFFFF );
       
  1200     if ( uiMagA64 < uiMagB64 ) goto returnB;
       
  1201     if ( uiMagB64 < uiMagA64 ) goto returnA;
       
  1202     if ( uiA0 < uiB0 ) goto returnB;
       
  1203     if ( uiB0 < uiA0 ) goto returnA;
       
  1204     if ( uiNonsigA64 < uiNonsigB64 ) goto returnA;
       
  1205  returnB:
       
  1206     uiZ.v64 = uiNonsigB64;
       
  1207     uiZ.v0  = uiB0;
       
  1208     return uiZ;
       
  1209  returnA:
       
  1210     uiZ.v64 = uiNonsigA64;
       
  1211     uiZ.v0  = uiA0;
       
  1212     return uiZ;
       
  1213 
       
  1214 }
       
  1215 
       
  1216 #if __POINTER_SIZE__ == 4
       
  1217 
       
  1218 void
       
  1219 i32_to_f128M( int32_t a, float128_t *zPtr )
       
  1220 {
       
  1221     uint32_t *zWPtr;
       
  1222     uint32_t uiZ96, uiZ64;
       
  1223     bool sign;
       
  1224     uint32_t absA;
       
  1225     int_fast8_t shiftDist;
       
  1226     uint64_t normAbsA;
       
  1227 
       
  1228     zWPtr = (uint32_t *) zPtr;
       
  1229     uiZ96 = 0;
       
  1230     uiZ64 = 0;
       
  1231     if ( a ) {
       
  1232 	sign = (a < 0);
       
  1233 	absA = sign ? -(uint32_t) a : (uint32_t) a;
       
  1234 	shiftDist = softfloat_countLeadingZeros32( absA ) + 17;
       
  1235 	normAbsA = (uint64_t) absA<<shiftDist;
       
  1236 	uiZ96 = packToF128UI96( sign, 0x402E - shiftDist, normAbsA>>32 );
       
  1237 	uiZ64 = normAbsA;
       
  1238     }
       
  1239     zWPtr[indexWord( 4, 3 )] = uiZ96;
       
  1240     zWPtr[indexWord( 4, 2 )] = uiZ64;
       
  1241     zWPtr[indexWord( 4, 1 )] = 0;
       
  1242     zWPtr[indexWord( 4, 0 )] = 0;
       
  1243 }
       
  1244 
       
  1245 #else
       
  1246 
       
  1247 void
       
  1248 i64_to_f128M( int64_t a, float128_t *zPtr )
       
  1249 {
       
  1250     uint32_t *zWPtr;
       
  1251     uint32_t uiZ96, uiZ64;
       
  1252     bool sign;
       
  1253     uint64_t absA;
       
  1254     uint_fast8_t shiftDist;
       
  1255     uint32_t *ptr;
       
  1256 
       
  1257     zWPtr = (uint32_t *) zPtr;
       
  1258     uiZ96 = 0;
       
  1259     uiZ64 = 0;
       
  1260     zWPtr[indexWord( 4, 1 )] = 0;
       
  1261     zWPtr[indexWord( 4, 0 )] = 0;
       
  1262     if ( a ) {
       
  1263 	sign = (a < 0);
       
  1264 	absA = sign ? -(uint64_t) a : (uint64_t) a;
       
  1265 	shiftDist = softfloat_countLeadingZeros64( absA ) + 17;
       
  1266 	if ( shiftDist < 32 ) {
       
  1267 	    ptr = zWPtr + indexMultiwordHi( 4, 3 );
       
  1268 	    ptr[indexWord( 3, 2 )] = 0;
       
  1269 	    ptr[indexWord( 3, 1 )] = absA>>32;
       
  1270 	    ptr[indexWord( 3, 0 )] = absA;
       
  1271 	    softfloat_shortShiftLeft96M( ptr, shiftDist, ptr );
       
  1272 	    ptr[indexWordHi( 3 )] =
       
  1273 		packToF128UI96(
       
  1274 		    sign, 0x404E - shiftDist, ptr[indexWordHi( 3 )] );
       
  1275 	    return;
       
  1276 	}
       
  1277 	absA <<= shiftDist - 32;
       
  1278 	uiZ96 = packToF128UI96( sign, 0x404E - shiftDist, absA>>32 );
       
  1279 	uiZ64 = absA;
       
  1280     }
       
  1281     zWPtr[indexWord( 4, 3 )] = uiZ96;
       
  1282     zWPtr[indexWord( 4, 2 )] = uiZ64;
       
  1283 }
       
  1284 
       
  1285 #endif
       
  1286 
       
  1287 void
       
  1288 f64_to_f128M( float64_t a, float128_t *zPtr )
       
  1289 {
       
  1290     uint32_t *zWPtr;
       
  1291     union ui64_f64 uA;
       
  1292     uint64_t uiA;
       
  1293     bool sign;
       
  1294     int_fast16_t exp;
       
  1295     uint64_t frac;
       
  1296     struct commonNaN commonNaN;
       
  1297     uint32_t uiZ96;
       
  1298     struct exp16_sig64 normExpSig;
       
  1299 
       
  1300     /*------------------------------------------------------------------------
       
  1301     *------------------------------------------------------------------------*/
       
  1302     zWPtr = (uint32_t *) zPtr;
       
  1303     /*------------------------------------------------------------------------
       
  1304     *------------------------------------------------------------------------*/
       
  1305     uA.f = a;
       
  1306     uiA = uA.ui;
       
  1307     sign = signF64UI( uiA );
       
  1308     exp  = expF64UI( uiA );
       
  1309     frac = fracF64UI( uiA );
       
  1310     /*------------------------------------------------------------------------
       
  1311     *------------------------------------------------------------------------*/
       
  1312     zWPtr[indexWord( 4, 0 )] = 0;
       
  1313     if ( exp == 0x7FF ) {
       
  1314 	if ( frac ) {
       
  1315 	    // NaN
       
  1316 	    softfloat_commonNaNToF128M( zWPtr );
       
  1317 	    return;
       
  1318 	}
       
  1319 	uiZ96 = packToF128UI96( sign, 0x7FFF, 0 );
       
  1320 	goto uiZ;
       
  1321     }
       
  1322     /*------------------------------------------------------------------------
       
  1323     *------------------------------------------------------------------------*/
       
  1324     if ( ! exp ) {
       
  1325 	if ( ! frac ) {
       
  1326 	    uiZ96 = packToF128UI96( sign, 0, 0 );
       
  1327 	    goto uiZ;
       
  1328 	}
       
  1329 	normExpSig = softfloat_normSubnormalF64Sig( frac );
       
  1330 	exp = normExpSig.exp - 1;
       
  1331 	frac = normExpSig.sig;
       
  1332     }
       
  1333     /*------------------------------------------------------------------------
       
  1334     *------------------------------------------------------------------------*/
       
  1335     zWPtr[indexWord( 4, 1 )] = (uint32_t) frac<<28;
       
  1336     frac >>= 4;
       
  1337     zWPtr[indexWordHi( 4 )] = packToF128UI96( sign, exp + 0x3C00, frac>>32 );
       
  1338     zWPtr[indexWord( 4, 2 )] = frac;
       
  1339     return;
       
  1340     /*------------------------------------------------------------------------
       
  1341     *------------------------------------------------------------------------*/
       
  1342  uiZ:
       
  1343     zWPtr[indexWord( 4, 3 )] = uiZ96;
       
  1344     zWPtr[indexWord( 4, 2 )] = 0;
       
  1345     zWPtr[indexWord( 4, 1 )] = 0;
       
  1346 
       
  1347 }
       
  1348 
       
  1349 float128_t
       
  1350 extF80_to_f128( extFloat80_t a )
       
  1351 {
       
  1352     union { struct extFloat80M s; extFloat80_t f; } uA;
       
  1353     uint_fast16_t uiA64;
       
  1354     uint_fast64_t uiA0;
       
  1355     uint_fast16_t exp;
       
  1356     uint_fast64_t frac;
       
  1357     struct commonNaN commonNaN;
       
  1358     struct uint128 uiZ;
       
  1359     bool sign;
       
  1360     struct uint128 frac128;
       
  1361     union ui128_f128 uZ;
       
  1362 
       
  1363     uA.f = a;
       
  1364     uiA64 = uA.s.signExp;
       
  1365     uiA0  = uA.s.signif;
       
  1366     exp = expExtF80UI64( uiA64 );
       
  1367     frac = uiA0 & UINT64_C( 0x7FFFFFFFFFFFFFFF );
       
  1368     if ( (exp == 0x7FFF) && frac ) {
       
  1369 	softfloat_extF80UIToCommonNaN( uiA64, uiA0, &commonNaN );
       
  1370 	softfloat_commonNaNToF128M( (uint32_t*) &uiZ );
       
  1371     } else {
       
  1372 	sign = signExtF80UI64( uiA64 );
       
  1373 	frac128 = softfloat_shortShiftLeft128( 0, frac, 49 );
       
  1374 	uiZ.v64 = packToF128UI64( sign, exp, frac128.v64 );
       
  1375 	uiZ.v0  = frac128.v0;
       
  1376     }
       
  1377     uZ.ui = uiZ;
       
  1378     return uZ.f;
       
  1379 
       
  1380 }
       
  1381 
       
  1382 static float128_t
       
  1383 softfloat_addMagsF128(
       
  1384      uint64_t uiA64,
       
  1385      uint64_t uiA0,
       
  1386      uint64_t uiB64,
       
  1387      uint64_t uiB0,
       
  1388      bool signZ
       
  1389  )
       
  1390 {
       
  1391     int32_t expA;
       
  1392     struct uint128 sigA;
       
  1393     int32_t expB;
       
  1394     struct uint128 sigB;
       
  1395     int32_t expDiff;
       
  1396     struct uint128 uiZ, sigZ;
       
  1397     int32_t expZ;
       
  1398     uint64_t sigZExtra;
       
  1399     struct uint128_extra sig128Extra;
       
  1400     union ui128_f128 uZ;
       
  1401 
       
  1402     expA = expF128UI64( uiA64 );
       
  1403     sigA.v64 = fracF128UI64( uiA64 );
       
  1404     sigA.v0  = uiA0;
       
  1405     expB = expF128UI64( uiB64 );
       
  1406     sigB.v64 = fracF128UI64( uiB64 );
       
  1407     sigB.v0  = uiB0;
       
  1408     expDiff = expA - expB;
       
  1409     if ( ! expDiff ) {
       
  1410 	if ( expA == 0x7FFF ) {
       
  1411 	    if ( sigA.v64 | sigA.v0 | sigB.v64 | sigB.v0 ) goto propagateNaN;
       
  1412 	    uiZ.v64 = uiA64;
       
  1413 	    uiZ.v0  = uiA0;
       
  1414 	    goto uiZ;
       
  1415 	}
       
  1416 	sigZ = softfloat_add128( sigA.v64, sigA.v0, sigB.v64, sigB.v0 );
       
  1417 	if ( ! expA ) {
       
  1418 	    uiZ.v64 = packToF128UI64( signZ, 0, sigZ.v64 );
       
  1419 	    uiZ.v0  = sigZ.v0;
       
  1420 	    goto uiZ;
       
  1421 	}
       
  1422 	expZ = expA;
       
  1423 	sigZ.v64 |= UINT64_C( 0x0002000000000000 );
       
  1424 	sigZExtra = 0;
       
  1425 	goto shiftRight1;
       
  1426     }
       
  1427     if ( expDiff < 0 ) {
       
  1428 	if ( expB == 0x7FFF ) {
       
  1429 	    if ( sigB.v64 | sigB.v0 ) goto propagateNaN;
       
  1430 	    uiZ.v64 = packToF128UI64( signZ, 0x7FFF, 0 );
       
  1431 	    uiZ.v0  = 0;
       
  1432 	    goto uiZ;
       
  1433 	}
       
  1434 	expZ = expB;
       
  1435 	if ( expA ) {
       
  1436 	    sigA.v64 |= UINT64_C( 0x0001000000000000 );
       
  1437 	} else {
       
  1438 	    ++expDiff;
       
  1439 	    sigZExtra = 0;
       
  1440 	    if ( ! expDiff ) goto newlyAligned;
       
  1441 	}
       
  1442 	sig128Extra =
       
  1443 	    softfloat_shiftRightJam128Extra( sigA.v64, sigA.v0, 0, -expDiff );
       
  1444 	sigA = sig128Extra.v;
       
  1445 	sigZExtra = sig128Extra.extra;
       
  1446     } else {
       
  1447 	if ( expA == 0x7FFF ) {
       
  1448 	    if ( sigA.v64 | sigA.v0 ) goto propagateNaN;
       
  1449 	    uiZ.v64 = uiA64;
       
  1450 	    uiZ.v0  = uiA0;
       
  1451 	    goto uiZ;
       
  1452 	}
       
  1453 	expZ = expA;
       
  1454 	if ( expB ) {
       
  1455 	    sigB.v64 |= UINT64_C( 0x0001000000000000 );
       
  1456 	} else {
       
  1457 	    --expDiff;
       
  1458 	    sigZExtra = 0;
       
  1459 	    if ( ! expDiff ) goto newlyAligned;
       
  1460 	}
       
  1461 	sig128Extra =
       
  1462 	    softfloat_shiftRightJam128Extra( sigB.v64, sigB.v0, 0, expDiff );
       
  1463 	sigB = sig128Extra.v;
       
  1464 	sigZExtra = sig128Extra.extra;
       
  1465     }
       
  1466  newlyAligned:
       
  1467     sigZ =
       
  1468 	softfloat_add128(
       
  1469 	    sigA.v64 | UINT64_C( 0x0001000000000000 ),
       
  1470 	    sigA.v0,
       
  1471 	    sigB.v64,
       
  1472 	    sigB.v0
       
  1473 	);
       
  1474     --expZ;
       
  1475     if ( sigZ.v64 < UINT64_C( 0x0002000000000000 ) ) goto roundAndPack;
       
  1476     ++expZ;
       
  1477  shiftRight1:
       
  1478     sig128Extra =
       
  1479 	softfloat_shortShiftRightJam128Extra(
       
  1480 	    sigZ.v64, sigZ.v0, sigZExtra, 1 );
       
  1481     sigZ = sig128Extra.v;
       
  1482     sigZExtra = sig128Extra.extra;
       
  1483  roundAndPack:
       
  1484     return
       
  1485 	softfloat_roundPackToF128( signZ, expZ, sigZ.v64, sigZ.v0, sigZExtra );
       
  1486  propagateNaN:
       
  1487     uiZ = softfloat_propagateNaNF128UI( uiA64, uiA0, uiB64, uiB0 );
       
  1488  uiZ:
       
  1489     uZ.ui = uiZ;
       
  1490     return uZ.f;
       
  1491 
       
  1492 }
       
  1493 
       
  1494 float128_t
       
  1495 softfloat_subMagsF128(
       
  1496      uint_fast64_t uiA64,
       
  1497      uint_fast64_t uiA0,
       
  1498      uint_fast64_t uiB64,
       
  1499      uint_fast64_t uiB0,
       
  1500      bool signZ
       
  1501  )
       
  1502 {
       
  1503     int_fast32_t expA;
       
  1504     struct uint128 sigA;
       
  1505     int_fast32_t expB;
       
  1506     struct uint128 sigB, sigZ;
       
  1507     int_fast32_t expDiff, expZ;
       
  1508     struct uint128 uiZ;
       
  1509     union ui128_f128 uZ;
       
  1510 
       
  1511     expA = expF128UI64( uiA64 );
       
  1512     sigA.v64 = fracF128UI64( uiA64 );
       
  1513     sigA.v0  = uiA0;
       
  1514     expB = expF128UI64( uiB64 );
       
  1515     sigB.v64 = fracF128UI64( uiB64 );
       
  1516     sigB.v0  = uiB0;
       
  1517     sigA = softfloat_shortShiftLeft128( sigA.v64, sigA.v0, 4 );
       
  1518     sigB = softfloat_shortShiftLeft128( sigB.v64, sigB.v0, 4 );
       
  1519     expDiff = expA - expB;
       
  1520     if ( 0 < expDiff ) goto expABigger;
       
  1521     if ( expDiff < 0 ) goto expBBigger;
       
  1522     if ( expA == 0x7FFF ) {
       
  1523 	if ( sigA.v64 | sigA.v0 | sigB.v64 | sigB.v0 ) goto propagateNaN;
       
  1524 	softfloat_raiseFlags( softfloat_flag_invalid );
       
  1525 	softfloat_commonNaNToF128M( (uint32_t*)&uiZ );
       
  1526 	// uiZ.v64 = defaultNaNF128UI64;
       
  1527 	// uiZ.v0  = defaultNaNF128UI0;
       
  1528 	goto uiZ;
       
  1529     }
       
  1530     expZ = expA;
       
  1531     if ( ! expZ ) expZ = 1;
       
  1532     if ( sigB.v64 < sigA.v64 ) goto aBigger;
       
  1533     if ( sigA.v64 < sigB.v64 ) goto bBigger;
       
  1534     if ( sigB.v0 < sigA.v0 ) goto aBigger;
       
  1535     if ( sigA.v0 < sigB.v0 ) goto bBigger;
       
  1536     uiZ.v64 =
       
  1537 	packToF128UI64(
       
  1538 	    (softfloat_roundingMode == softfloat_round_min), 0, 0 );
       
  1539     uiZ.v0 = 0;
       
  1540     goto uiZ;
       
  1541  expBBigger:
       
  1542     if ( expB == 0x7FFF ) {
       
  1543 	if ( sigB.v64 | sigB.v0 ) goto propagateNaN;
       
  1544 	uiZ.v64 = packToF128UI64( signZ ^ 1, 0x7FFF, 0 );
       
  1545 	uiZ.v0  = 0;
       
  1546 	goto uiZ;
       
  1547     }
       
  1548     if ( expA ) {
       
  1549 	sigA.v64 |= UINT64_C( 0x0010000000000000 );
       
  1550     } else {
       
  1551 	++expDiff;
       
  1552 	if ( ! expDiff ) goto newlyAlignedBBigger;
       
  1553     }
       
  1554     sigA = softfloat_shiftRightJam128( sigA.v64, sigA.v0, -expDiff );
       
  1555  newlyAlignedBBigger:
       
  1556     expZ = expB;
       
  1557     sigB.v64 |= UINT64_C( 0x0010000000000000 );
       
  1558  bBigger:
       
  1559     signZ = ! signZ;
       
  1560     sigZ = softfloat_sub128( sigB.v64, sigB.v0, sigA.v64, sigA.v0 );
       
  1561     goto normRoundPack;
       
  1562  expABigger:
       
  1563     if ( expA == 0x7FFF ) {
       
  1564 	if ( sigA.v64 | sigA.v0 ) goto propagateNaN;
       
  1565 	uiZ.v64 = uiA64;
       
  1566 	uiZ.v0  = uiA0;
       
  1567 	goto uiZ;
       
  1568     }
       
  1569     if ( expB ) {
       
  1570 	sigB.v64 |= UINT64_C( 0x0010000000000000 );
       
  1571     } else {
       
  1572 	--expDiff;
       
  1573 	if ( ! expDiff ) goto newlyAlignedABigger;
       
  1574     }
       
  1575     sigB = softfloat_shiftRightJam128( sigB.v64, sigB.v0, expDiff );
       
  1576  newlyAlignedABigger:
       
  1577     expZ = expA;
       
  1578     sigA.v64 |= UINT64_C( 0x0010000000000000 );
       
  1579  aBigger:
       
  1580     sigZ = softfloat_sub128( sigA.v64, sigA.v0, sigB.v64, sigB.v0 );
       
  1581  normRoundPack:
       
  1582     return softfloat_normRoundPackToF128( signZ, expZ - 5, sigZ.v64, sigZ.v0 );
       
  1583  propagateNaN:
       
  1584     uiZ = softfloat_propagateNaNF128UI( uiA64, uiA0, uiB64, uiB0 );
       
  1585  uiZ:
       
  1586     uZ.ui = uiZ;
       
  1587     return uZ.f;
       
  1588 }
       
  1589 
       
  1590 void
       
  1591 f128M_add( const float128_t *aPtr, const float128_t *bPtr, float128_t *zPtr )
       
  1592 {
       
  1593     const uint64_t *aWPtr, *bWPtr;
       
  1594     uint64_t uiA64, uiA0;
       
  1595     int signA;
       
  1596     uint64_t uiB64, uiB0;
       
  1597     int signB;
       
  1598 
       
  1599     aWPtr = (const uint64_t *) aPtr;
       
  1600     bWPtr = (const uint64_t *) bPtr;
       
  1601     uiA64 = aWPtr[indexWord( 2, 1 )];
       
  1602     uiA0  = aWPtr[indexWord( 2, 0 )];
       
  1603     signA = signF128UI64( uiA64 );
       
  1604     uiB64 = bWPtr[indexWord( 2, 1 )];
       
  1605     uiB0  = bWPtr[indexWord( 2, 0 )];
       
  1606     signB = signF128UI64( uiB64 );
       
  1607 
       
  1608     if ( signA == signB ) {
       
  1609 	*zPtr = softfloat_addMagsF128( uiA64, uiA0, uiB64, uiB0, signA );
       
  1610     } else {
       
  1611 	*zPtr = softfloat_subMagsF128( uiA64, uiA0, uiB64, uiB0, signA );
       
  1612     }
       
  1613 }
       
  1614 
       
  1615 void
       
  1616 f128M_sub( const float128_t *aPtr, const float128_t *bPtr, float128_t *zPtr )
       
  1617 {
       
  1618     const uint64_t *aWPtr, *bWPtr;
       
  1619     uint_fast64_t uiA64, uiA0;
       
  1620     bool signA;
       
  1621     uint_fast64_t uiB64, uiB0;
       
  1622     bool signB;
       
  1623 
       
  1624     aWPtr = (const uint64_t *) aPtr;
       
  1625     bWPtr = (const uint64_t *) bPtr;
       
  1626     uiA64 = aWPtr[indexWord( 2, 1 )];
       
  1627     uiA0  = aWPtr[indexWord( 2, 0 )];
       
  1628     signA = signF128UI64( uiA64 );
       
  1629     uiB64 = bWPtr[indexWord( 2, 1 )];
       
  1630     uiB0  = bWPtr[indexWord( 2, 0 )];
       
  1631     signB = signF128UI64( uiB64 );
       
  1632     if ( signA == signB ) {
       
  1633 	*zPtr = softfloat_subMagsF128( uiA64, uiA0, uiB64, uiB0, signA );
       
  1634     } else {
       
  1635 	*zPtr = softfloat_addMagsF128( uiA64, uiA0, uiB64, uiB0, signA );
       
  1636     }
       
  1637 }
       
  1638 
       
  1639 void
       
  1640 f128M_mul( const float128_t *aPtr, const float128_t *bPtr, float128_t *zPtr )
       
  1641 {
       
  1642     const uint32_t *aWPtr, *bWPtr;
       
  1643     uint32_t *zWPtr;
       
  1644     uint32_t uiA96;
       
  1645     int32_t expA;
       
  1646     uint32_t uiB96;
       
  1647     int32_t expB;
       
  1648     bool signZ;
       
  1649     const uint32_t *ptr;
       
  1650     uint32_t uiZ96, sigA[4];
       
  1651     uint_fast8_t shiftDist;
       
  1652     uint32_t sigB[4];
       
  1653     int32_t expZ;
       
  1654     uint32_t sigProd[8], *extSigZPtr;
       
  1655 
       
  1656     /*------------------------------------------------------------------------
       
  1657     *------------------------------------------------------------------------*/
       
  1658     aWPtr = (const uint32_t *) aPtr;
       
  1659     bWPtr = (const uint32_t *) bPtr;
       
  1660     zWPtr = (uint32_t *) zPtr;
       
  1661     /*------------------------------------------------------------------------
       
  1662     *------------------------------------------------------------------------*/
       
  1663     uiA96 = aWPtr[indexWordHi( 4 )];
       
  1664     expA = expF128UI96( uiA96 );
       
  1665     uiB96 = bWPtr[indexWordHi( 4 )];
       
  1666     expB = expF128UI96( uiB96 );
       
  1667     signZ = signF128UI96( uiA96 ) ^ signF128UI96( uiB96 );
       
  1668     /*------------------------------------------------------------------------
       
  1669     *------------------------------------------------------------------------*/
       
  1670     if ( (expA == 0x7FFF) || (expB == 0x7FFF) ) {
       
  1671 	if ( softfloat_tryPropagateNaNF128M( aWPtr, bWPtr, zWPtr ) ) return;
       
  1672 	ptr = aWPtr;
       
  1673 	if ( ! expA ) goto possiblyInvalid;
       
  1674 	if ( ! expB ) {
       
  1675 	    ptr = bWPtr;
       
  1676      possiblyInvalid:
       
  1677 	    if (
       
  1678 		! fracF128UI96( ptr[indexWordHi( 4 )] )
       
  1679 		    && ! (ptr[indexWord( 4, 2 )] | ptr[indexWord( 4, 1 )]
       
  1680 			      | ptr[indexWord( 4, 0 )])
       
  1681 	    ) {
       
  1682 		softfloat_invalidF128M( zWPtr );
       
  1683 		return;
       
  1684 	    }
       
  1685 	}
       
  1686 	uiZ96 = packToF128UI96( signZ, 0x7FFF, 0 );
       
  1687 	goto uiZ96;
       
  1688     }
       
  1689     /*------------------------------------------------------------------------
       
  1690     *------------------------------------------------------------------------*/
       
  1691     if ( expA ) {
       
  1692 	sigA[indexWordHi( 4 )] = fracF128UI96( uiA96 ) | 0x00010000;
       
  1693 	sigA[indexWord( 4, 2 )] = aWPtr[indexWord( 4, 2 )];
       
  1694 	sigA[indexWord( 4, 1 )] = aWPtr[indexWord( 4, 1 )];
       
  1695 	sigA[indexWord( 4, 0 )] = aWPtr[indexWord( 4, 0 )];
       
  1696     } else {
       
  1697 	expA = softfloat_shiftNormSigF128M( aWPtr, 0, sigA );
       
  1698 	if ( expA == -128 ) goto zero;
       
  1699     }
       
  1700     if ( expB ) {
       
  1701 	sigB[indexWordHi( 4 )] = fracF128UI96( uiB96 ) | 0x00010000;
       
  1702 	sigB[indexWord( 4, 2 )] = bWPtr[indexWord( 4, 2 )];
       
  1703 	sigB[indexWord( 4, 1 )] = bWPtr[indexWord( 4, 1 )];
       
  1704 	sigB[indexWord( 4, 0 )] = bWPtr[indexWord( 4, 0 )];
       
  1705     } else {
       
  1706 	expB = softfloat_shiftNormSigF128M( bWPtr, 0, sigB );
       
  1707 	if ( expB == -128 ) goto zero;
       
  1708     }
       
  1709     /*------------------------------------------------------------------------
       
  1710     *------------------------------------------------------------------------*/
       
  1711     expZ = expA + expB - 0x4000;
       
  1712     softfloat_mul128MTo256M( sigA, sigB, sigProd );
       
  1713     if (
       
  1714 	sigProd[indexWord( 8, 2 )]
       
  1715 	    || (sigProd[indexWord( 8, 1 )] | sigProd[indexWord( 8, 0 )])
       
  1716     ) {
       
  1717 	sigProd[indexWord( 8, 3 )] |= 1;
       
  1718     }
       
  1719     extSigZPtr = &sigProd[indexMultiwordHi( 8, 5 )];
       
  1720     shiftDist = 16;
       
  1721     if ( extSigZPtr[indexWordHi( 5 )] & 2 ) {
       
  1722 	++expZ;
       
  1723 	shiftDist = 15;
       
  1724     }
       
  1725     softfloat_shortShiftLeft160M( extSigZPtr, shiftDist, extSigZPtr );
       
  1726     softfloat_roundPackMToF128M( signZ, expZ, extSigZPtr, zWPtr );
       
  1727     return;
       
  1728     /*------------------------------------------------------------------------
       
  1729     *------------------------------------------------------------------------*/
       
  1730  zero:
       
  1731     uiZ96 = packToF128UI96( signZ, 0, 0 );
       
  1732  uiZ96:
       
  1733     zWPtr[indexWordHi( 4 )] = uiZ96;
       
  1734     zWPtr[indexWord( 4, 2 )] = 0;
       
  1735     zWPtr[indexWord( 4, 1 )] = 0;
       
  1736     zWPtr[indexWord( 4, 0 )] = 0;
       
  1737 }
       
  1738 
       
  1739 float128_t
       
  1740 f128_div( float128_t a, float128_t b )
       
  1741 {
       
  1742     union ui128_f128 uA;
       
  1743     uint_fast64_t uiA64, uiA0;
       
  1744     bool signA;
       
  1745     int_fast32_t expA;
       
  1746     struct uint128 sigA;
       
  1747     union ui128_f128 uB;
       
  1748     uint_fast64_t uiB64, uiB0;
       
  1749     bool signB;
       
  1750     int_fast32_t expB;
       
  1751     struct uint128 sigB;
       
  1752     bool signZ;
       
  1753     struct exp32_sig128 normExpSig;
       
  1754     int_fast32_t expZ;
       
  1755     struct uint128 rem;
       
  1756     uint_fast32_t recip32;
       
  1757     int ix;
       
  1758     uint_fast64_t q64;
       
  1759     uint_fast32_t q;
       
  1760     struct uint128 term;
       
  1761     uint_fast32_t qs[3];
       
  1762     uint_fast64_t sigZExtra;
       
  1763     struct uint128 sigZ, uiZ;
       
  1764     union ui128_f128 uZ;
       
  1765 
       
  1766     /*------------------------------------------------------------------------
       
  1767     *------------------------------------------------------------------------*/
       
  1768     uA.f = a;
       
  1769     uiA64 = uA.ui.v64;
       
  1770     uiA0  = uA.ui.v0;
       
  1771     signA = signF128UI64( uiA64 );
       
  1772     expA  = expF128UI64( uiA64 );
       
  1773     sigA.v64 = fracF128UI64( uiA64 );
       
  1774     sigA.v0  = uiA0;
       
  1775     uB.f = b;
       
  1776     uiB64 = uB.ui.v64;
       
  1777     uiB0  = uB.ui.v0;
       
  1778     signB = signF128UI64( uiB64 );
       
  1779     expB  = expF128UI64( uiB64 );
       
  1780     sigB.v64 = fracF128UI64( uiB64 );
       
  1781     sigB.v0  = uiB0;
       
  1782     signZ = signA ^ signB;
       
  1783     /*------------------------------------------------------------------------
       
  1784     *------------------------------------------------------------------------*/
       
  1785     if ( expA == 0x7FFF ) {
       
  1786 	if ( sigA.v64 | sigA.v0 ) goto propagateNaN;
       
  1787 	if ( expB == 0x7FFF ) {
       
  1788 	    if ( sigB.v64 | sigB.v0 ) goto propagateNaN;
       
  1789 	    goto invalid;
       
  1790 	}
       
  1791 	goto infinity;
       
  1792     }
       
  1793     if ( expB == 0x7FFF ) {
       
  1794 	if ( sigB.v64 | sigB.v0 ) goto propagateNaN;
       
  1795 	goto zero;
       
  1796     }
       
  1797     /*------------------------------------------------------------------------
       
  1798     *------------------------------------------------------------------------*/
       
  1799     if ( ! expB ) {
       
  1800 	if ( ! (sigB.v64 | sigB.v0) ) {
       
  1801 	    if ( ! (expA | sigA.v64 | sigA.v0) ) goto invalid;
       
  1802 	    softfloat_raiseFlags( softfloat_flag_infinite );
       
  1803 	    goto infinity;
       
  1804 	}
       
  1805 	normExpSig = softfloat_normSubnormalF128Sig( sigB.v64, sigB.v0 );
       
  1806 	expB = normExpSig.exp;
       
  1807 	sigB = normExpSig.sig;
       
  1808     }
       
  1809     if ( ! expA ) {
       
  1810 	if ( ! (sigA.v64 | sigA.v0) ) goto zero;
       
  1811 	normExpSig = softfloat_normSubnormalF128Sig( sigA.v64, sigA.v0 );
       
  1812 	expA = normExpSig.exp;
       
  1813 	sigA = normExpSig.sig;
       
  1814     }
       
  1815     /*------------------------------------------------------------------------
       
  1816     *------------------------------------------------------------------------*/
       
  1817     expZ = expA - expB + 0x3FFE;
       
  1818     sigA.v64 |= UINT64_C( 0x0001000000000000 );
       
  1819     sigB.v64 |= UINT64_C( 0x0001000000000000 );
       
  1820     rem = sigA;
       
  1821     if ( softfloat_lt128( sigA.v64, sigA.v0, sigB.v64, sigB.v0 ) ) {
       
  1822 	--expZ;
       
  1823 	rem = softfloat_add128( sigA.v64, sigA.v0, sigA.v64, sigA.v0 );
       
  1824     }
       
  1825     recip32 = softfloat_approxRecip32_1( sigB.v64>>17 );
       
  1826     ix = 3;
       
  1827     for (;;) {
       
  1828 	q64 = (uint_fast64_t) (uint32_t) (rem.v64>>19) * recip32;
       
  1829 	q = (q64 + 0x80000000)>>32;
       
  1830 	--ix;
       
  1831 	if ( ix < 0 ) break;
       
  1832 	rem = softfloat_shortShiftLeft128( rem.v64, rem.v0, 29 );
       
  1833 	term = softfloat_mul128By32( sigB.v64, sigB.v0, q );
       
  1834 	rem = softfloat_sub128( rem.v64, rem.v0, term.v64, term.v0 );
       
  1835 	if ( rem.v64 & UINT64_C( 0x8000000000000000 ) ) {
       
  1836 	    --q;
       
  1837 	    rem = softfloat_add128( rem.v64, rem.v0, sigB.v64, sigB.v0 );
       
  1838 	}
       
  1839 	qs[ix] = q;
       
  1840     }
       
  1841     /*------------------------------------------------------------------------
       
  1842     *------------------------------------------------------------------------*/
       
  1843     if ( ((q + 1) & 7) < 2 ) {
       
  1844 	rem = softfloat_shortShiftLeft128( rem.v64, rem.v0, 29 );
       
  1845 	term = softfloat_mul128By32( sigB.v64, sigB.v0, q );
       
  1846 	rem = softfloat_sub128( rem.v64, rem.v0, term.v64, term.v0 );
       
  1847 	if ( rem.v64 & UINT64_C( 0x8000000000000000 ) ) {
       
  1848 	    --q;
       
  1849 	    rem = softfloat_add128( rem.v64, rem.v0, sigB.v64, sigB.v0 );
       
  1850 	} else if ( softfloat_le128( sigB.v64, sigB.v0, rem.v64, rem.v0 ) ) {
       
  1851 	    ++q;
       
  1852 	    rem = softfloat_sub128( rem.v64, rem.v0, sigB.v64, sigB.v0 );
       
  1853 	}
       
  1854 	if ( rem.v64 | rem.v0 ) q |= 1;
       
  1855     }
       
  1856     /*------------------------------------------------------------------------
       
  1857     *------------------------------------------------------------------------*/
       
  1858     sigZExtra = (uint64_t) ((uint_fast64_t) q<<60);
       
  1859     term = softfloat_shortShiftLeft128( 0, qs[1], 54 );
       
  1860     sigZ =
       
  1861 	softfloat_add128(
       
  1862 	    (uint_fast64_t) qs[2]<<19, ((uint_fast64_t) qs[0]<<25) + (q>>4),
       
  1863 	    term.v64, term.v0
       
  1864 	);
       
  1865     return
       
  1866 	softfloat_roundPackToF128( signZ, expZ, sigZ.v64, sigZ.v0, sigZExtra );
       
  1867     /*------------------------------------------------------------------------
       
  1868     *------------------------------------------------------------------------*/
       
  1869  propagateNaN:
       
  1870     uiZ = softfloat_propagateNaNF128UI( uiA64, uiA0, uiB64, uiB0 );
       
  1871     goto uiZ;
       
  1872     /*------------------------------------------------------------------------
       
  1873     *------------------------------------------------------------------------*/
       
  1874  invalid:
       
  1875     softfloat_raiseFlags( softfloat_flag_invalid );
       
  1876     softfloat_commonNaNToF128M( (uint32_t*)&uiZ );
       
  1877     // uiZ.v64 = defaultNaNF128UI64;
       
  1878     // uiZ.v0  = defaultNaNF128UI0;
       
  1879     goto uiZ;
       
  1880     /*------------------------------------------------------------------------
       
  1881     *------------------------------------------------------------------------*/
       
  1882  infinity:
       
  1883     uiZ.v64 = packToF128UI64( signZ, 0x7FFF, 0 );
       
  1884     goto uiZ0;
       
  1885     /*------------------------------------------------------------------------
       
  1886     *------------------------------------------------------------------------*/
       
  1887  zero:
       
  1888     uiZ.v64 = packToF128UI64( signZ, 0, 0 );
       
  1889  uiZ0:
       
  1890     uiZ.v0 = 0;
       
  1891  uiZ:
       
  1892     uZ.ui = uiZ;
       
  1893     return uZ.f;
       
  1894 }
       
  1895 
       
  1896 static float128_t
       
  1897 f128_rem( float128_t a, float128_t b )
       
  1898 {
       
  1899     union ui128_f128 uA;
       
  1900     uint_fast64_t uiA64, uiA0;
       
  1901     bool signA;
       
  1902     int_fast32_t expA;
       
  1903     struct uint128 sigA;
       
  1904     union ui128_f128 uB;
       
  1905     uint_fast64_t uiB64, uiB0;
       
  1906     int_fast32_t expB;
       
  1907     struct uint128 sigB;
       
  1908     struct exp32_sig128 normExpSig;
       
  1909     struct uint128 rem;
       
  1910     int_fast32_t expDiff;
       
  1911     uint_fast32_t q, recip32;
       
  1912     uint_fast64_t q64;
       
  1913     struct uint128 term, altRem, meanRem;
       
  1914     bool signRem;
       
  1915     struct uint128 uiZ;
       
  1916     union ui128_f128 uZ;
       
  1917 
       
  1918     /*------------------------------------------------------------------------
       
  1919     *------------------------------------------------------------------------*/
       
  1920     uA.f = a;
       
  1921     uiA64 = uA.ui.v64;
       
  1922     uiA0  = uA.ui.v0;
       
  1923     signA = signF128UI64( uiA64 );
       
  1924     expA  = expF128UI64( uiA64 );
       
  1925     sigA.v64 = fracF128UI64( uiA64 );
       
  1926     sigA.v0  = uiA0;
       
  1927     uB.f = b;
       
  1928     uiB64 = uB.ui.v64;
       
  1929     uiB0  = uB.ui.v0;
       
  1930     expB  = expF128UI64( uiB64 );
       
  1931     sigB.v64 = fracF128UI64( uiB64 );
       
  1932     sigB.v0  = uiB0;
       
  1933     /*------------------------------------------------------------------------
       
  1934     *------------------------------------------------------------------------*/
       
  1935     if ( expA == 0x7FFF ) {
       
  1936 	if (
       
  1937 	    (sigA.v64 | sigA.v0) || ((expB == 0x7FFF) && (sigB.v64 | sigB.v0))
       
  1938 	) {
       
  1939 	    goto propagateNaN;
       
  1940 	}
       
  1941 	goto invalid;
       
  1942     }
       
  1943     if ( expB == 0x7FFF ) {
       
  1944 	if ( sigB.v64 | sigB.v0 ) goto propagateNaN;
       
  1945 	return a;
       
  1946     }
       
  1947     /*------------------------------------------------------------------------
       
  1948     *------------------------------------------------------------------------*/
       
  1949     if ( ! expB ) {
       
  1950 	if ( ! (sigB.v64 | sigB.v0) ) goto invalid;
       
  1951 	normExpSig = softfloat_normSubnormalF128Sig( sigB.v64, sigB.v0 );
       
  1952 	expB = normExpSig.exp;
       
  1953 	sigB = normExpSig.sig;
       
  1954     }
       
  1955     if ( ! expA ) {
       
  1956 	if ( ! (sigA.v64 | sigA.v0) ) return a;
       
  1957 	normExpSig = softfloat_normSubnormalF128Sig( sigA.v64, sigA.v0 );
       
  1958 	expA = normExpSig.exp;
       
  1959 	sigA = normExpSig.sig;
       
  1960     }
       
  1961     /*------------------------------------------------------------------------
       
  1962     *------------------------------------------------------------------------*/
       
  1963     sigA.v64 |= UINT64_C( 0x0001000000000000 );
       
  1964     sigB.v64 |= UINT64_C( 0x0001000000000000 );
       
  1965     rem = sigA;
       
  1966     expDiff = expA - expB;
       
  1967     if ( expDiff < 1 ) {
       
  1968 	if ( expDiff < -1 ) return a;
       
  1969 	if ( expDiff ) {
       
  1970 	    --expB;
       
  1971 	    sigB = softfloat_add128( sigB.v64, sigB.v0, sigB.v64, sigB.v0 );
       
  1972 	    q = 0;
       
  1973 	} else {
       
  1974 	    q = softfloat_le128( sigB.v64, sigB.v0, rem.v64, rem.v0 );
       
  1975 	    if ( q ) {
       
  1976 		rem = softfloat_sub128( rem.v64, rem.v0, sigB.v64, sigB.v0 );
       
  1977 	    }
       
  1978 	}
       
  1979     } else {
       
  1980 	recip32 = softfloat_approxRecip32_1( sigB.v64>>17 );
       
  1981 	expDiff -= 30;
       
  1982 	for (;;) {
       
  1983 	    q64 = (uint_fast64_t) (uint32_t) (rem.v64>>19) * recip32;
       
  1984 	    if ( expDiff < 0 ) break;
       
  1985 	    q = (q64 + 0x80000000)>>32;
       
  1986 	    rem = softfloat_shortShiftLeft128( rem.v64, rem.v0, 29 );
       
  1987 	    term = softfloat_mul128By32( sigB.v64, sigB.v0, q );
       
  1988 	    rem = softfloat_sub128( rem.v64, rem.v0, term.v64, term.v0 );
       
  1989 	    if ( rem.v64 & UINT64_C( 0x8000000000000000 ) ) {
       
  1990 		rem = softfloat_add128( rem.v64, rem.v0, sigB.v64, sigB.v0 );
       
  1991 	    }
       
  1992 	    expDiff -= 29;
       
  1993 	}
       
  1994 	/*--------------------------------------------------------------------
       
  1995 	| (`expDiff' cannot be less than -29 here.)
       
  1996 	*--------------------------------------------------------------------*/
       
  1997 	q = (uint32_t) (q64>>32)>>(~expDiff & 31);
       
  1998 	rem = softfloat_shortShiftLeft128( rem.v64, rem.v0, expDiff + 30 );
       
  1999 	term = softfloat_mul128By32( sigB.v64, sigB.v0, q );
       
  2000 	rem = softfloat_sub128( rem.v64, rem.v0, term.v64, term.v0 );
       
  2001 	if ( rem.v64 & UINT64_C( 0x8000000000000000 ) ) {
       
  2002 	    altRem = softfloat_add128( rem.v64, rem.v0, sigB.v64, sigB.v0 );
       
  2003 	    goto selectRem;
       
  2004 	}
       
  2005     }
       
  2006     /*------------------------------------------------------------------------
       
  2007     *------------------------------------------------------------------------*/
       
  2008     do {
       
  2009 	altRem = rem;
       
  2010 	++q;
       
  2011 	rem = softfloat_sub128( rem.v64, rem.v0, sigB.v64, sigB.v0 );
       
  2012     } while ( ! (rem.v64 & UINT64_C( 0x8000000000000000 )) );
       
  2013  selectRem:
       
  2014     meanRem = softfloat_add128( rem.v64, rem.v0, altRem.v64, altRem.v0 );
       
  2015     if (
       
  2016 	(meanRem.v64 & UINT64_C( 0x8000000000000000 ))
       
  2017 	    || (! (meanRem.v64 | meanRem.v0) && (q & 1))
       
  2018     ) {
       
  2019 	rem = altRem;
       
  2020     }
       
  2021     signRem = signA;
       
  2022     if ( rem.v64 & UINT64_C( 0x8000000000000000 ) ) {
       
  2023 	signRem = ! signRem;
       
  2024 	rem = softfloat_sub128( 0, 0, rem.v64, rem.v0 );
       
  2025     }
       
  2026     return softfloat_normRoundPackToF128( signRem, expB - 1, rem.v64, rem.v0 );
       
  2027     /*------------------------------------------------------------------------
       
  2028     *------------------------------------------------------------------------*/
       
  2029  propagateNaN:
       
  2030     uiZ = softfloat_propagateNaNF128UI( uiA64, uiA0, uiB64, uiB0 );
       
  2031     goto uiZ;
       
  2032     /*------------------------------------------------------------------------
       
  2033     *------------------------------------------------------------------------*/
       
  2034  invalid:
       
  2035     softfloat_raiseFlags( softfloat_flag_invalid );
       
  2036     softfloat_commonNaNToF128M( (uint32_t*)&uiZ );
       
  2037     // uiZ.v64 = defaultNaNF128UI64;
       
  2038     // uiZ.v0  = defaultNaNF128UI0;
       
  2039  uiZ:
       
  2040     uZ.ui = uiZ;
       
  2041     return uZ.f;
       
  2042 }
       
  2043 
       
  2044 bool
       
  2045 f128_isNan( float128_t a) {
       
  2046     int_fast32_t expA;
       
  2047     struct uint128 sigA;
       
  2048     uint_fast64_t uiA64, uiA0;
       
  2049     union ui128_f128 uA;
       
  2050 
       
  2051     uA.f = a;
       
  2052     uiA64 = uA.ui.v64;
       
  2053     uiA0  = uA.ui.v0;
       
  2054     return isNaNF128UI( uiA64, uiA0 );
       
  2055 }
       
  2056 
       
  2057 bool
       
  2058 f128_isInf( float128_t a) {
       
  2059     int_fast32_t expA;
       
  2060     struct uint128 sigA;
       
  2061     uint_fast64_t uiA64, uiA0;
       
  2062     union ui128_f128 uA;
       
  2063 
       
  2064     uA.f = a;
       
  2065     uiA64 = uA.ui.v64;
       
  2066     uiA0  = uA.ui.v0;
       
  2067     return isInfF128UI( uiA64, uiA0 );
       
  2068 //    expA  = expF128UI64( uiA64 );
       
  2069 //    sigA.v64 = fracF128UI64( uiA64 );
       
  2070 //    sigA.v0  = uiA0;
       
  2071 //
       
  2072 //    if ( expA == 0x7FFF ) {
       
  2073 //        if ( sigA.v64 == 0 ) return 1;
       
  2074 //    }
       
  2075 //    return 0;
       
  2076 }
       
  2077 
       
  2078 static bool
       
  2079 f128M_eq( float128_t* aPtr, float128_t* bPtr )
       
  2080 {
       
  2081     union ui128_f128 uA;
       
  2082     uint_fast64_t uiA64, uiA0;
       
  2083     union ui128_f128 uB;
       
  2084     uint_fast64_t uiB64, uiB0;
       
  2085 
       
  2086     uA.f = *aPtr;
       
  2087     uiA64 = uA.ui.v64;
       
  2088     uiA0  = uA.ui.v0;
       
  2089     uB.f = *bPtr;
       
  2090     uiB64 = uB.ui.v64;
       
  2091     uiB0  = uB.ui.v0;
       
  2092 
       
  2093     // fprintf(stderr, "isNaNF128UI(a): %d\n", isNaNF128UI( uiA64, uiA0 ));
       
  2094     // fprintf(stderr, "isNaNF128UI(a): %d\n", isNaNF128UI( uiB64, uiB0 ));
       
  2095 
       
  2096     if ( isNaNF128UI( uiA64, uiA0 ) || isNaNF128UI( uiB64, uiB0 ) ) {
       
  2097 	if (
       
  2098 	       softfloat_isSigNaNF128UI( uiA64, uiA0 )
       
  2099 	    || softfloat_isSigNaNF128UI( uiB64, uiB0 )
       
  2100 	) {
       
  2101 	    softfloat_raiseFlags( softfloat_flag_invalid );
       
  2102 	}
       
  2103 	return 0;
       
  2104     }
       
  2105 
       
  2106     // fprintf(stderr, "uiA0 == uiB0: %d\n", (uiA0 == uiB0));
       
  2107     // fprintf(stderr, "uiA64 == uiB64: %d\n", (uiA64 == uiB64));
       
  2108 
       
  2109     return
       
  2110 	   (uiA0 == uiB0)
       
  2111 	&& (   (uiA64 == uiB64)
       
  2112 	    || (! uiA0 && ! ((uiA64 | uiB64) & UINT64_C( 0x7FFFFFFFFFFFFFFF )))
       
  2113 	   );
       
  2114 }
       
  2115 
       
  2116 static bool
       
  2117 f128M_lt( float128_t* aPtr, float128_t* bPtr )
       
  2118 {
       
  2119     union ui128_f128 uA;
       
  2120     uint_fast64_t uiA64, uiA0;
       
  2121     union ui128_f128 uB;
       
  2122     uint_fast64_t uiB64, uiB0;
       
  2123     bool signA, signB;
       
  2124 
       
  2125     uA.f = *aPtr;
       
  2126     uiA64 = uA.ui.v64;
       
  2127     uiA0  = uA.ui.v0;
       
  2128     uB.f = *bPtr;
       
  2129     uiB64 = uB.ui.v64;
       
  2130     uiB0  = uB.ui.v0;
       
  2131     if ( isNaNF128UI( uiA64, uiA0 ) || isNaNF128UI( uiB64, uiB0 ) ) {
       
  2132 	softfloat_raiseFlags( softfloat_flag_invalid );
       
  2133 	return 0;
       
  2134     }
       
  2135     signA = signF128UI64( uiA64 );
       
  2136     signB = signF128UI64( uiB64 );
       
  2137     return
       
  2138 	(signA != signB)
       
  2139 	    ? signA
       
  2140 		  && (((uiA64 | uiB64) & UINT64_C( 0x7FFFFFFFFFFFFFFF ))
       
  2141 			  | uiA0 | uiB0)
       
  2142 	    : ((uiA64 != uiB64) || (uiA0 != uiB0))
       
  2143 		  && (signA ^ softfloat_lt128( uiA64, uiA0, uiB64, uiB0 ));
       
  2144 }
       
  2145 
       
  2146 static inline void
       
  2147 f128M_div( const float128_t *aPtr, const float128_t *bPtr, float128_t *zPtr )
       
  2148 {
       
  2149     *zPtr = f128_div( *aPtr, *bPtr );
       
  2150 }
       
  2151 
       
  2152 static inline void
       
  2153 f128M_rem( const float128_t *aPtr, const float128_t *bPtr, float128_t *zPtr )
       
  2154 {
       
  2155     *zPtr = f128_rem( *aPtr, *bPtr );
       
  2156 }
       
  2157 
       
  2158 static inline float128_t
       
  2159 f128_negate( float128_t a) {
       
  2160     struct uint128 sigA;
       
  2161     uint_fast64_t uiA64, uiA0;
       
  2162     union ui128_f128 uA;
       
  2163 
       
  2164     uA.f = a;
       
  2165     uiA64 = uA.ui.v64;
       
  2166     uiA0  = uA.ui.v0;
       
  2167     if (isNaNF128UI( uiA64, uiA0 )) return a;
       
  2168     uA.ui.v64 ^= UINT64_C(0x8000000000000000);
       
  2169     return uA.f;
       
  2170 }
       
  2171 
       
  2172 static inline void
       
  2173 f128M_negate( const float128_t *aPtr, float128_t *zPtr )
       
  2174 {
       
  2175     *zPtr = f128_negate( *aPtr );
       
  2176 }
       
  2177 
       
  2178 #endif // SUPPORT_QUADFLOAT
       
  2179 
   107 %}
  2180 %}
   108 ! !
  2181 ! !
   109 
  2182 
   110 !QuadFloat class methodsFor:'documentation'!
  2183 !QuadFloat class methodsFor:'documentation'!
   111 
  2184 
   165     OBJ newFloat;
  2238     OBJ newFloat;
   166 
  2239 
   167     if (sizeof(long double) == sizeof(float128_t)) {
  2240     if (sizeof(long double) == sizeof(float128_t)) {
   168 	__qMKLFLOAT(newFloat, 0.0);   /* OBJECT ALLOCATION */
  2241 	__qMKLFLOAT(newFloat, 0.0);   /* OBJECT ALLOCATION */
   169     } else {
  2242     } else {
   170 	float128_t qf = STX_qZero;
  2243 	float128_t qf;
       
  2244 	f64_to_f128M(0.0, &qf);
   171 	__qMKQFLOAT(newFloat, qf);   /* OBJECT ALLOCATION */
  2245 	__qMKQFLOAT(newFloat, qf);   /* OBJECT ALLOCATION */
   172     }
  2246     }
   173     RETURN (newFloat);
  2247     RETURN (newFloat);
   174 #endif /* SUPPORT_QUADFLOAT */
  2248 #endif /* SUPPORT_QUADFLOAT */
   175 %}.
  2249 %}.
   184 %{  /* NOCONTEXT */
  2258 %{  /* NOCONTEXT */
   185 #ifdef SUPPORT_QUADFLOAT
  2259 #ifdef SUPPORT_QUADFLOAT
   186     OBJ newFloat;
  2260     OBJ newFloat;
   187 
  2261 
   188     if (__isFloatLike(aFloat)) {
  2262     if (__isFloatLike(aFloat)) {
   189 	double f = __floatVal(aFloat);
  2263 	float64_t f = __floatVal(aFloat);
   190 	float128_t qf;
  2264 	float128_t qf;
   191 
  2265 
   192 	qf = STX_dbltoq (f);
  2266 	f64_to_f128M(f, &qf);
   193 	__qMKQFLOAT(newFloat, qf);   /* OBJECT ALLOCATION */
  2267 	__qMKQFLOAT(newFloat, qf);   /* OBJECT ALLOCATION */
   194 	RETURN (newFloat);
  2268 	RETURN (newFloat);
   195     }
  2269     }
   196 #endif /* SUPPORT_QUADFLOAT */
  2270 #endif /* SUPPORT_QUADFLOAT */
   197 %}.
  2271 %}.
   218 
  2292 
   219     if (__isSmallInteger(anInteger)) {
  2293     if (__isSmallInteger(anInteger)) {
   220 	INT iVal = __intVal(anInteger);
  2294 	INT iVal = __intVal(anInteger);
   221 	float128_t qf;
  2295 	float128_t qf;
   222 
  2296 
   223 	qf = STX_inttoq( (long)iVal );
  2297 #if __POINTER_SIZE__ == 4
       
  2298 	i32_to_f128M( iVal, &qf );
       
  2299 else
       
  2300 	i64_to_f128M( iVal, &qf );
       
  2301 #endif
   224 	__qMKQFLOAT(newFloat, qf);   /* OBJECT ALLOCATION */
  2302 	__qMKQFLOAT(newFloat, qf);   /* OBJECT ALLOCATION */
   225 	RETURN (newFloat);
  2303 	RETURN (newFloat);
   226     }
  2304     }
   227 #endif /* SUPPORT_QUADFLOAT */
  2305 #endif /* SUPPORT_QUADFLOAT */
   228 %}.
  2306 %}.
   236 
  2314 
   237 fromLongFloat:aFloat
  2315 fromLongFloat:aFloat
   238     "return a new quadFloat, given a long float value"
  2316     "return a new quadFloat, given a long float value"
   239 
  2317 
   240 %{  /* NOCONTEXT */
  2318 %{  /* NOCONTEXT */
   241 #ifdef xSUPPORT_QUADFLOAT
  2319 #ifdef SUPPORT_QUADFLOAT
   242     OBJ newFloat;
  2320     OBJ newFloat;
   243     union {
  2321     union {
   244 	LONGFLOAT_t lf;         // is long double
  2322 	LONGFLOAT_t lf;         // is long double
   245 	extFloat80_t ef;        // is 80bit ext
  2323 	extFloat80_t ef;        // is 80bit ext
   246 	float128_t qf;          // is 128bit quad
  2324 	float128_t qf;          // is 128bit quad
   284     "
  2362     "
   285 
  2363 
   286     "Created: / 08-06-2019 / 03:28:37 / Claus Gittinger"
  2364     "Created: / 08-06-2019 / 03:28:37 / Claus Gittinger"
   287 ! !
  2365 ! !
   288 
  2366 
   289 !QuadFloat class methodsFor:'accessing'!
       
   290 
       
   291 defaultPrintFormat
       
   292     "/ by default, I will print 19 digits
       
   293     "/  ShortFloat pi   -> 3.141593  
       
   294     "/  Float pi        -> 3.14159265358979  
       
   295     "/  LongFloat pi    -> 3.141592653589793239  
       
   296     "/  QuadFloat pi    -> 3.141592653589793239  
       
   297 
       
   298     ^ DefaultPrintFormat
       
   299 ! !
       
   300 
       
   301 !QuadFloat class methodsFor:'class initialization'!
       
   302 
       
   303 initialize
       
   304     DefaultPrintFormat := '.34'.  "/ 34 valid digits
       
   305     DefaultPrintfFormat := '%34f'.
       
   306 
       
   307     "
       
   308      self initialize
       
   309     "
       
   310 
       
   311     "
       
   312      DefaultPrintFormat := '.19'.
       
   313      QuadFloat pi printString.
       
   314 
       
   315      DefaultPrintFormat := '.9'.
       
   316      QuadFloat pi printString.
       
   317 
       
   318      DefaultPrintFormat := '.6'.
       
   319      QuadFloat pi printString.
       
   320     "
       
   321 ! !
       
   322 
       
   323 !QuadFloat class methodsFor:'coercing & converting'!
  2367 !QuadFloat class methodsFor:'coercing & converting'!
   324 
  2368 
   325 coerce:aNumber
  2369 coerce:aNumber
   326     "convert the argument aNumber into an instance of the receiver's class and return it."
  2370     "convert the argument aNumber into an instance of the receiver's class and return it."
   327 
  2371 
   337 
  2381 
   338     |nan|
  2382     |nan|
   339 
  2383 
   340     NaN isNil ifTrue:[
  2384     NaN isNil ifTrue:[
   341 %{  /* NOCONTEXT */
  2385 %{  /* NOCONTEXT */
   342 #ifdef xSUPPORT_QUADFLOAT
  2386 #ifdef SUPPORT_QUADFLOAT
   343 	{
  2387 	{
   344 	    OBJ newFloat;
  2388 	    OBJ newFloat;
   345 	    float128_t qf;
  2389 	    float128_t qf;
   346 
  2390 
   347 	    softfloat_commonNaNToF128M( (uint32_t*)(&qf) );
  2391 	    softfloat_commonNaNToF128M( (uint32_t*)(&qf) );
   373 
  2417 
   374     |inf|
  2418     |inf|
   375 
  2419 
   376     PositiveInfinity isNil ifTrue:[
  2420     PositiveInfinity isNil ifTrue:[
   377 %{  /* NOCONTEXT */
  2421 %{  /* NOCONTEXT */
   378 #ifdef xSUPPORT_QUADFLOAT
  2422 #ifdef SUPPORT_QUADFLOAT
   379 	{
  2423 	{
   380 	    OBJ newFloat;
  2424 	    OBJ newFloat;
   381 	    struct uint128 uiZ;
  2425 	    struct uint128 uiZ;
   382 	    union ui128_f128 uZ;
  2426 	    union ui128_f128 uZ;
   383 	    float128_t qf;
  2427 	    float128_t qf;
   403 
  2447 
   404     |inf|
  2448     |inf|
   405 
  2449 
   406     NegativeInfinity isNil ifTrue:[
  2450     NegativeInfinity isNil ifTrue:[
   407 %{  /* NOCONTEXT */
  2451 %{  /* NOCONTEXT */
   408 #ifdef xSUPPORT_QUADFLOAT
  2452 #ifdef SUPPORT_QUADFLOAT
   409 	{
  2453 	{
   410 	    OBJ newFloat;
  2454 	    OBJ newFloat;
   411 	    struct uint128 uiZ;
  2455 	    struct uint128 uiZ;
   412 	    union ui128_f128 uZ;
  2456 	    union ui128_f128 uZ;
   413 	    float128_t qf;
  2457 	    float128_t qf;
   430 
  2474 
   431 phi
  2475 phi
   432     "return the constant phi as quadFloat"
  2476     "return the constant phi as quadFloat"
   433 
  2477 
   434     Phi isNil ifTrue:[
  2478     Phi isNil ifTrue:[
   435 	"/ phiDigits has enough digits for 128bit IEEE quads
  2479         "/ phiDigits has enough digits for 128bit IEEE quads
   436 	"/ do not use as a literal constant here - we cannot depend on the underlying C-compiler here...
  2480         "/ do not use as a literal constant here - we cannot depend on the underlying C-compiler here...
   437 	Phi  := self readFrom:(Number phiDigits)
  2481         Phi  := self readFrom:(Number phiDigits)
   438     ].
  2482     ].
   439     ^ Phi
  2483     ^ Phi
   440 
  2484 
   441 
  2485 
   442 !
  2486 !
   566 
  2610 
   567 * aNumber
  2611 * aNumber
   568     "return the product of the receiver and the argument."
  2612     "return the product of the receiver and the argument."
   569 
  2613 
   570     aNumber class == QuadFloat ifTrue:[
  2614     aNumber class == QuadFloat ifTrue:[
   571 	^ aNumber productFromQuadFloat:self
  2615         ^ aNumber productFromQuadFloat:self
   572     ].
  2616     ].
   573 
  2617 
   574     thisContext isReallyRecursive ifTrue:[self error].
  2618     thisContext isReallyRecursive ifTrue:[self error].
   575     ^ aNumber productFromQuadFloat:self
  2619     ^ aNumber productFromQuadFloat:self
   576 !
  2620 !
   597 	^ ZeroDivide raiseRequestWith:thisContext.
  2641 	^ ZeroDivide raiseRequestWith:thisContext.
   598     ].
  2642     ].
   599     ^ aNumber quotientFromQuadFloat:self
  2643     ^ aNumber quotientFromQuadFloat:self
   600 !
  2644 !
   601 
  2645 
   602 abs
  2646 negated
   603     "return the absolute value of the receiver
  2647     "return the receiver negated"
   604      reimplemented here for speed"
       
   605 
  2648 
   606 %{  /* NOCONTEXT */
  2649 %{  /* NOCONTEXT */
   607 #ifdef SUPPORT_QUADFLOAT
  2650 #ifdef SUPPORT_QUADFLOAT
   608     OBJ newFloat;
  2651     OBJ newFloat;
   609     float128_t result, myVal;
  2652     float128_t result, myVal;
   610 
  2653 
   611     myVal = __quadFloatVal(self);
  2654     myVal = __quadFloatVal(self);
   612     result = STX_qabs(myVal);
  2655     f128M_negate( &myVal, &result );
   613     __qMKQFLOAT(newFloat, result);
       
   614     RETURN ( newFloat );
       
   615 #endif
       
   616 %}.
       
   617 !
       
   618 
       
   619 ceiling
       
   620     "return the smallest integer which is greater or equal to the receiver."
       
   621 
       
   622     |val|
       
   623 
       
   624 %{
       
   625     float128_t qVal;
       
   626     float128_t qMinInt;
       
   627     float128_t qMaxInt;
       
   628 
       
   629     qVal = __quadFloatVal(self);
       
   630     qVal = STX_qceil(qVal);
       
   631 
       
   632     qMinInt = STX_dbltoq((double)_MIN_INT);
       
   633     qMaxInt = STX_dbltoq((double)_MAX_INT);
       
   634     if (STX_qge(qVal, qMinInt) && STX_qle(qVal, qMaxInt)) {
       
   635 	double dVal = STX_qtodbl(qVal);
       
   636 	RETURN ( __mkSmallInteger( (INT) dVal ) );
       
   637     }
       
   638     __qMKQFLOAT(val, qVal);
       
   639 %}.
       
   640     ^ val asInteger
       
   641 
       
   642     "
       
   643      0.5 asQuadFloat ceiling
       
   644      0.5 asQuadFloat ceilingAsFloat
       
   645      -0.5 asQuadFloat ceiling
       
   646      -0.5 asQuadFloat ceilingAsFloat
       
   647     "
       
   648 !
       
   649 
       
   650 ceilingAsFloat
       
   651     "return the smallest integer-valued float greater or equal to the receiver.
       
   652      This is much like #ceiling, but avoids a (possibly expensive) conversion
       
   653      of the result to an integer.
       
   654      It may be useful, if the result is to be further used in another float-operation."
       
   655 
       
   656 %{  /* NOCONTEXT */
       
   657 #ifdef SUPPORT_QUADFLOAT
       
   658     OBJ newFloat;
       
   659     float128_t result, myVal;
       
   660 
       
   661     myVal = __quadFloatVal(self);
       
   662     result = STX_qceil(myVal);
       
   663     __qMKQFLOAT(newFloat, result);
       
   664     RETURN ( newFloat );
       
   665 #endif
       
   666 %}.
       
   667 !
       
   668 
       
   669 cos
       
   670     "return the cosine of the receiver (interpreted as radians)"
       
   671 
       
   672 %{  /* NOCONTEXT */
       
   673 #ifdef SUPPORT_QUADFLOAT
       
   674     OBJ newFloat;
       
   675     float128_t result, myVal;
       
   676 
       
   677     myVal = __quadFloatVal(self);
       
   678     result = STX_qcos(myVal);
       
   679     __qMKQFLOAT(newFloat, result);
       
   680     RETURN ( newFloat );
       
   681 #endif
       
   682 %}.
       
   683 !
       
   684 
       
   685 cosh
       
   686     "return the hyperbolic cosine of the receiver (interpreted as radians)"
       
   687 
       
   688 %{  /* NOCONTEXT */
       
   689 #ifdef SUPPORT_QUADFLOAT
       
   690     OBJ newFloat;
       
   691     float128_t result, myVal;
       
   692 
       
   693     myVal = __quadFloatVal(self);
       
   694     result = STX_qcosh(myVal);
       
   695     __qMKQFLOAT(newFloat, result);
       
   696     RETURN ( newFloat );
       
   697 #endif
       
   698 %}.
       
   699 !
       
   700 
       
   701 exp
       
   702     "return e raised to the power of the receiver"
       
   703 
       
   704 %{  /* NOCONTEXT */
       
   705 #ifdef SUPPORT_QUADFLOAT
       
   706     OBJ newFloat;
       
   707     float128_t result, myVal;
       
   708 
       
   709     myVal = __quadFloatVal(self);
       
   710     result = STX_qexp(myVal);
       
   711     __qMKQFLOAT(newFloat, result);
       
   712     RETURN ( newFloat );
       
   713 #endif
       
   714 %}.
       
   715 !
       
   716 
       
   717 exponent
       
   718     "extract a normalized float's exponent.
       
   719      The returned value depends on the float-representation of
       
   720      the underlying machine and is therefore highly unportable.
       
   721      This is not for general use.
       
   722      This assumes that the mantissa is normalized to
       
   723      0.5 .. 1.0 and the float's value is: mantissa * 2^exp"
       
   724 
       
   725 %{  /* NOCONTEXT */
       
   726     float128_t myVal;
       
   727     float128_t frac;
       
   728     int exp;
       
   729 
       
   730     myVal = __quadFloatVal(self);
       
   731 #if 1
       
   732     // should we?
       
   733     if (! (STX_qisNan(&myVal) || STX_qisInf(&myVal)))
       
   734 #endif
       
   735     {
       
   736 	frac = STX_qfrexp(myVal, &exp);
       
   737 	RETURN (__mkSmallInteger(exp));
       
   738     }
       
   739 %}.
       
   740     ^ super exponent
       
   741 
       
   742     "
       
   743      1.0 exponent
       
   744      1.0 asQuadFloat exponent
       
   745      2.0 exponent
       
   746      2.0 asQuadFloat exponent
       
   747      3.0 exponent
       
   748      3.0 asQuadFloat exponent
       
   749      4.0 exponent
       
   750      4.0 asQuadFloat exponent
       
   751      0.5 exponent
       
   752      0.5 asQuadFloat exponent
       
   753      0.4 exponent
       
   754      0.4 asQuadFloat exponent
       
   755      0.25 exponent
       
   756      0.25 asQuadFloat exponent
       
   757      0.2 exponent
       
   758      0.2 asQuadFloat exponent
       
   759      0.00000011111 exponent
       
   760      0.00000011111 asQuadFloat exponent
       
   761      0.0 exponent
       
   762      0.0 asQuadFloat exponent
       
   763 
       
   764      1e1000 exponent -> error (INF)
       
   765      1e1000 asQuadFloat exponent -> error (INF)
       
   766     "
       
   767 !
       
   768 
       
   769 floor
       
   770     "return the integer nearest the receiver towards negative infinity."
       
   771 
       
   772     |val|
       
   773 
       
   774 %{
       
   775     float128_t qVal;
       
   776     float128_t qMinInt;
       
   777     float128_t qMaxInt;
       
   778 
       
   779     qVal = __quadFloatVal(self);
       
   780     qVal = STX_qfloor(qVal);
       
   781 
       
   782     qMinInt = STX_dbltoq((double)_MIN_INT);
       
   783     qMaxInt = STX_dbltoq((double)_MAX_INT);
       
   784     if (STX_qge(qVal, qMinInt) && STX_qle(qVal, qMaxInt)) {
       
   785 	double dVal = STX_qtodbl(qVal);
       
   786 	RETURN ( __mkSmallInteger( (INT) dVal ) );
       
   787     }
       
   788     __qMKQFLOAT(val, qVal);
       
   789 %}.
       
   790     ^ val asInteger
       
   791 
       
   792     "
       
   793      0.5 asQuadFloat floor
       
   794      0.5 asQuadFloat floorAsFloat
       
   795      -0.5 asQuadFloat floor
       
   796      -0.5 asQuadFloat floorAsFloat
       
   797     "
       
   798 !
       
   799 
       
   800 floorAsFloat
       
   801     "return the integer nearest the receiver towards negative infinity as a float.
       
   802      This is much like #floor, but avoids a (possibly expensive) conversion
       
   803      of the result to an integer.
       
   804      It may be useful, if the result is to be further used in another float-operation."
       
   805 
       
   806 %{  /* NOCONTEXT */
       
   807 #ifdef SUPPORT_QUADFLOAT
       
   808     OBJ newFloat;
       
   809     float128_t result, myVal;
       
   810 
       
   811     myVal = __quadFloatVal(self);
       
   812     result = STX_qfloor(myVal);
       
   813     __qMKQFLOAT(newFloat, result);
       
   814     RETURN ( newFloat );
       
   815 #endif
       
   816 %}.
       
   817 !
       
   818 
       
   819 ln
       
   820     "return natural logarithm of the receiver."
       
   821 
       
   822 %{  /* NOCONTEXT */
       
   823 #ifdef SUPPORT_QUADFLOAT
       
   824     OBJ newFloat;
       
   825     float128_t result, myVal;
       
   826 
       
   827     myVal = __quadFloatVal(self);
       
   828     result = STX_qlog(myVal);
       
   829     __qMKQFLOAT(newFloat, result);
       
   830     RETURN ( newFloat );
       
   831 #endif
       
   832 %}.
       
   833 !
       
   834 
       
   835 log
       
   836     "return log base 10 of the receiver.
       
   837      Alias for log:10."
       
   838 
       
   839 %{  /* NOCONTEXT */
       
   840 #ifdef SUPPORT_QUADFLOAT
       
   841     OBJ newFloat;
       
   842     float128_t result, myVal;
       
   843 
       
   844     myVal = __quadFloatVal(self);
       
   845     result = STX_qlog10(myVal);
       
   846     __qMKQFLOAT(newFloat, result);
       
   847     RETURN ( newFloat );
       
   848 #endif
       
   849 %}.
       
   850 !
       
   851 
       
   852 log2
       
   853     "return logarithm dualis of the receiver."
       
   854 
       
   855 %{  /* NOCONTEXT */
       
   856 #ifdef SUPPORT_QUADFLOAT
       
   857     OBJ newFloat;
       
   858     float128_t result, myVal;
       
   859 
       
   860     myVal = __quadFloatVal(self);
       
   861     result = STX_qlog2(myVal);
       
   862     __qMKQFLOAT(newFloat, result);
       
   863     RETURN ( newFloat );
       
   864 #endif
       
   865 %}.
       
   866 !
       
   867 
       
   868 mantissa
       
   869     "extract a normalized float's mantissa.
       
   870      The returned value depends on the float-representation of
       
   871      the underlying machine and is therefore highly unportable.
       
   872      This is not for general use.
       
   873      This assumes that the mantissa is normalized to
       
   874      0.5 .. 1.0 and the float's value is mantissa * 2^exp"
       
   875 
       
   876 %{  /* NOCONTEXT */
       
   877     float128_t myVal;
       
   878     float128_t frac;
       
   879     int exp;
       
   880     OBJ newFloat;
       
   881 
       
   882     myVal = __quadFloatVal(self);
       
   883     // ouch: math libs seem to not care for NaN here;
       
   884 #if 1
       
   885     // should we?
       
   886     if (! (STX_qisNan(&myVal) || STX_qisInf(&myVal)))
       
   887 #endif
       
   888     {
       
   889 	frac = STX_qfrexp(myVal, &exp);
       
   890 	__qMKQFLOAT(newFloat, frac);
       
   891 	RETURN ( newFloat );
       
   892     }
       
   893 %}.
       
   894     ^ super mantissa
       
   895 
       
   896     "
       
   897      1.0 exponent
       
   898      1.0 asQuadFloat exponent
       
   899      1.0 mantissa
       
   900      1.0 asQuadFloat mantissa
       
   901 
       
   902      0.25 exponent
       
   903      0.25 asQuadFloat exponent
       
   904      0.25 mantissa
       
   905      0.25 asQuadFloat mantissa
       
   906 
       
   907      0.00000011111 exponent
       
   908      0.00000011111 mantissa
       
   909 
       
   910      1e1000 mantissa
       
   911     "
       
   912 
       
   913     "Modified: / 20-06-2017 / 11:37:13 / cg"
       
   914     "Modified (comment): / 26-05-2019 / 03:12:55 / Claus Gittinger"
       
   915 !
       
   916 
       
   917 negated
       
   918     "return the receiver negated"
       
   919 
       
   920 %{  /* NOCONTEXT */
       
   921 #ifdef SUPPORT_QUADFLOAT
       
   922     OBJ newFloat;
       
   923     float128_t result, myVal;
       
   924 
       
   925     myVal = __quadFloatVal(self);
       
   926     result = STX_qneg(myVal);
       
   927     __qMKQFLOAT(newFloat, result);
  2656     __qMKQFLOAT(newFloat, result);
   928     RETURN ( newFloat );
  2657     RETURN ( newFloat );
   929 #endif
  2658 #endif
   930 %}.
  2659 %}.
   931 !
  2660 !
   938 	 No, you shalt not divide by zero
  2667 	 No, you shalt not divide by zero
   939 	"
  2668 	"
   940 	^ ZeroDivide raiseRequestWith:thisContext.
  2669 	^ ZeroDivide raiseRequestWith:thisContext.
   941     ].
  2670     ].
   942     ^ aNumber remainderFromLongFloat:self
  2671     ^ aNumber remainderFromLongFloat:self
   943 !
       
   944 
       
   945 sin
       
   946     "return the sine of the receiver (interpreted as radians)"
       
   947 
       
   948 %{  /* NOCONTEXT */
       
   949 #ifdef SUPPORT_QUADFLOAT
       
   950     OBJ newFloat;
       
   951     float128_t result, myVal;
       
   952 
       
   953     myVal = __quadFloatVal(self);
       
   954     result = STX_qsin(myVal);
       
   955     __qMKQFLOAT(newFloat, result);
       
   956     RETURN ( newFloat );
       
   957 #endif
       
   958 %}.
       
   959 !
       
   960 
       
   961 sinh
       
   962     "return the hyperbolic sine of the receiver (interpreted as radians)"
       
   963 
       
   964 %{  /* NOCONTEXT */
       
   965 #ifdef SUPPORT_QUADFLOAT
       
   966     OBJ newFloat;
       
   967     float128_t result, myVal;
       
   968 
       
   969     myVal = __quadFloatVal(self);
       
   970     result = STX_qsinh(myVal);
       
   971     __qMKQFLOAT(newFloat, result);
       
   972     RETURN ( newFloat );
       
   973 #endif
       
   974 %}.
       
   975 !
       
   976 
       
   977 tan
       
   978     "return the tangent of the receiver (interpreted as radians)"
       
   979 
       
   980 %{  /* NOCONTEXT */
       
   981 #ifdef SUPPORT_QUADFLOAT
       
   982     OBJ newFloat;
       
   983     float128_t result, myVal;
       
   984 
       
   985     myVal = __quadFloatVal(self);
       
   986     result = STX_qtan(myVal);
       
   987     __qMKQFLOAT(newFloat, result);
       
   988     RETURN ( newFloat );
       
   989 #endif
       
   990 %}.
       
   991 !
       
   992 
       
   993 tanh
       
   994     "return the hyperbolic tangent of the receiver (interpreted as radians)"
       
   995 
       
   996 %{  /* NOCONTEXT */
       
   997 #ifdef SUPPORT_QUADFLOAT
       
   998     OBJ newFloat;
       
   999     float128_t result, myVal;
       
  1000 
       
  1001     myVal = __quadFloatVal(self);
       
  1002     result = STX_qtanh(myVal);
       
  1003     __qMKQFLOAT(newFloat, result);
       
  1004     RETURN ( newFloat );
       
  1005 #endif
       
  1006 %}.
       
  1007 ! !
  2672 ! !
  1008 
  2673 
  1009 !QuadFloat methodsFor:'coercing & converting'!
  2674 !QuadFloat methodsFor:'coercing & converting'!
  1010 
       
  1011 asFloat
       
  1012     "return a Float with same value as the receiver.
       
  1013      Raises an error if the receiver exceeds the float range."
       
  1014 
       
  1015 %{  /* NOCONTEXT */
       
  1016 
       
  1017     OBJ newFloat;
       
  1018     float128_t qVal = __quadFloatVal(self);
       
  1019     double dVal = STX_qtodbl(qVal);
       
  1020 
       
  1021     if (isfinite(dVal) || !STX_qisfinite(qVal)) {
       
  1022 	__qMKFLOAT(newFloat, dVal);
       
  1023 	RETURN ( newFloat );
       
  1024     }
       
  1025 %}.
       
  1026     "
       
  1027      value out of range
       
  1028      if you need -INF for a zero receiver, try Number trapInfinity:[...]
       
  1029     "
       
  1030     ^ self class
       
  1031 	raise:#infiniteResultSignal
       
  1032 	receiver:self
       
  1033 	selector:#asFloat
       
  1034 	arguments:#()
       
  1035 	errorString:'receiver is out of the double-precision float range'
       
  1036 
       
  1037     "
       
  1038      1.0 asQuadFloat asFloat
       
  1039     "
       
  1040 !
       
  1041 
  2675 
  1042 generality
  2676 generality
  1043     "return the generality value - see ArithmeticValue>>retry:coercing:"
  2677     "return the generality value - see ArithmeticValue>>retry:coercing:"
  1044 
  2678 
  1045     ^ 93
  2679     ^ 93
  1115     OBJ newFloat;
  2749     OBJ newFloat;
  1116     float128_t result, myVal, argVal;
  2750     float128_t result, myVal, argVal;
  1117 
  2751 
  1118     myVal = __quadFloatVal(self);
  2752     myVal = __quadFloatVal(self);
  1119     argVal = __quadFloatVal(aQuadFloat);
  2753     argVal = __quadFloatVal(aQuadFloat);
  1120     result = STX_qadd( argVal, myVal, 1 );
  2754     f128M_sub( &myVal, &argVal, &result );
  1121     __qMKQFLOAT(newFloat, result);
  2755     __qMKQFLOAT(newFloat, result);
  1122     RETURN ( newFloat );
  2756     RETURN ( newFloat );
  1123 #endif // SUPPORT_QUADFLOAT
  2757 #endif // SUPPORT_QUADFLOAT
  1124 %}.
  2758 %}.
  1125     self errorUnsupported
  2759     self errorUnsupported
  1133     OBJ newFloat;
  2767     OBJ newFloat;
  1134     float128_t result, myVal, argVal;
  2768     float128_t result, myVal, argVal;
  1135 
  2769 
  1136     myVal = __quadFloatVal(self);
  2770     myVal = __quadFloatVal(self);
  1137     argVal = __quadFloatVal(aQuadFloat);
  2771     argVal = __quadFloatVal(aQuadFloat);
  1138     RETURN (STX_qeq(argVal, myVal) ? true : false);
  2772     RETURN (f128M_eq( &argVal, &myVal ) ? true : false);
  1139 #endif // SUPPORT_QUADFLOAT
  2773 #endif // SUPPORT_QUADFLOAT
  1140 %}.
  2774 %}.
  1141     self errorUnsupported
  2775     self errorUnsupported
  1142 
  2776 
  1143     "Modified: / 08-06-2019 / 13:31:48 / Claus Gittinger"
  2777     "Modified: / 08-06-2019 / 13:31:48 / Claus Gittinger"
  1151     OBJ newFloat;
  2785     OBJ newFloat;
  1152     float128_t result, myVal, argVal;
  2786     float128_t result, myVal, argVal;
  1153 
  2787 
  1154     myVal = __quadFloatVal(self);
  2788     myVal = __quadFloatVal(self);
  1155     argVal = __quadFloatVal(aQuadFloat);
  2789     argVal = __quadFloatVal(aQuadFloat);
  1156     RETURN (STX_qlt(argVal, myVal) ? true : false);
  2790     RETURN (f128M_lt( &argVal, &myVal ) ? true : false);
  1157 #endif // SUPPORT_QUADFLOAT
  2791 #endif // SUPPORT_QUADFLOAT
  1158 %}.
  2792 %}.
  1159     self errorUnsupported
  2793     self errorUnsupported
  1160 !
  2794 !
  1161 
  2795 
  1167     OBJ newFloat;
  2801     OBJ newFloat;
  1168     float128_t result, myVal, argVal;
  2802     float128_t result, myVal, argVal;
  1169 
  2803 
  1170     myVal = __quadFloatVal(self);
  2804     myVal = __quadFloatVal(self);
  1171     argVal = __quadFloatVal(aQuadFloat);
  2805     argVal = __quadFloatVal(aQuadFloat);
  1172     result = STX_qmul(myVal, argVal);
  2806     f128M_mul( &myVal, &argVal, &result );
  1173     __qMKQFLOAT(newFloat, result);
  2807     __qMKQFLOAT(newFloat, result);
  1174     RETURN ( newFloat );
  2808     RETURN ( newFloat );
  1175 #endif // SUPPORT_QUADFLOAT
  2809 #endif // SUPPORT_QUADFLOAT
  1176 %}.
  2810 %}.
  1177     self errorUnsupported
  2811     self errorUnsupported
  1185     OBJ newFloat;
  2819     OBJ newFloat;
  1186     float128_t result, myVal, argVal;
  2820     float128_t result, myVal, argVal;
  1187 
  2821 
  1188     myVal = __quadFloatVal(self);
  2822     myVal = __quadFloatVal(self);
  1189     argVal = __quadFloatVal(aQuadFloat);
  2823     argVal = __quadFloatVal(aQuadFloat);
  1190     result = STX_qdiv(argVal, myVal);
  2824     f128M_div( &myVal, &argVal, &result );
  1191     __qMKQFLOAT(newFloat, result);
  2825     __qMKQFLOAT(newFloat, result);
  1192     RETURN ( newFloat );
  2826     RETURN ( newFloat );
  1193 #endif // SUPPORT_QUADFLOAT
  2827 #endif // SUPPORT_QUADFLOAT
  1194 %}.
  2828 %}.
  1195     self errorUnsupported
  2829     self errorUnsupported
  1203     OBJ newFloat;
  2837     OBJ newFloat;
  1204     float128_t result, myVal, argVal;
  2838     float128_t result, myVal, argVal;
  1205 
  2839 
  1206     myVal = __quadFloatVal(self);
  2840     myVal = __quadFloatVal(self);
  1207     argVal = __quadFloatVal(aQuadFloat);
  2841     argVal = __quadFloatVal(aQuadFloat);
  1208     result = STX_qadd( myVal, argVal, 0 );
  2842     f128M_add( &myVal, &argVal, &result );
  1209     __qMKQFLOAT(newFloat, result);
  2843     __qMKQFLOAT(newFloat, result);
  1210     RETURN ( newFloat );
  2844     RETURN ( newFloat );
  1211 #endif // SUPPORT_QUADFLOAT
  2845 #endif // SUPPORT_QUADFLOAT
  1212 %}.
  2846 %}.
  1213     self errorUnsupported
  2847     self errorUnsupported
  1222 ! !
  2856 ! !
  1223 
  2857 
  1224 !QuadFloat methodsFor:'printing'!
  2858 !QuadFloat methodsFor:'printing'!
  1225 
  2859 
  1226 printOn:aStream
  2860 printOn:aStream
  1227     PrintfScanf printf:'%e' on:aStream argument:self
  2861     |mantissa exponent|
  1228 
  2862 
  1229 "/    |mantissa exponent|
  2863     mantissa := self mantissa.
  1230 "/
  2864     exponent := self exponent.
  1231 "/    mantissa := self mantissa.
  2865 
  1232 "/    exponent := self exponent.
  2866     self exponent == 0 ifTrue:[
  1233 "/
  2867 	mantissa printOn:aStream.
  1234 "/    self exponent == 0 ifTrue:[
  2868 	aStream nextPutAll:'.0'.
  1235 "/        mantissa printOn:aStream.
  2869 	^ self
  1236 "/        aStream nextPutAll:'.0'.
  2870     ].
  1237 "/        ^ self
  2871     mantissa == 0 ifTrue:[
  1238 "/    ].
  2872 	"/ a zero mantissa is impossible - except for zero and a few others
  1239 "/    mantissa == 0 ifTrue:[
  2873 	exponent == 0 ifTrue:[ aStream nextPutAll:'0.0'. ^ self].
  1240 "/        "/ a zero mantissa is impossible - except for zero and a few others
  2874 	self == NaN ifTrue:[ aStream nextPutAll:'NAN'. ^ self ].
  1241 "/        exponent == 0 ifTrue:[ aStream nextPutAll:'0.0'. ^ self].
  2875 	self == NegativeInfinity ifTrue:[ aStream nextPutAll:'-INF'. ^ self].
  1242 "/        self == NaN ifTrue:[ aStream nextPutAll:'NAN'. ^ self ].
  2876 	self == PositiveInfinity ifTrue:[ aStream nextPutAll:'INF'. ^ self].
  1243 "/        self == NegativeInfinity ifTrue:[ aStream nextPutAll:'-INF'. ^ self].
  2877 
  1244 "/        self == PositiveInfinity ifTrue:[ aStream nextPutAll:'INF'. ^ self].
  2878 	self error:'invalid largeFloat' mayProceed:true.
  1245 "/
  2879 	aStream nextPutAll:'Invalid'. ^ self.
  1246 "/        self error:'invalid largeFloat' mayProceed:true.
  2880     ].
  1247 "/        aStream nextPutAll:'Invalid'. ^ self.
  2881 
  1248 "/    ].
  2882     exponent >= 0 ifTrue:[
  1249 "/
  2883 	(mantissa bitShift:exponent) printOn:aStream.
  1250 "/    exponent >= 0 ifTrue:[
  2884 	aStream nextPutAll:'.0'.
  1251 "/        (mantissa bitShift:exponent) printOn:aStream.
  2885 	^ self
  1252 "/        aStream nextPutAll:'.0'.
  2886     ].
  1253 "/        ^ self
  2887     ((mantissa / (1 bitShift:exponent negated)) asFixedPoint:6) printOn:aStream.
  1254 "/    ].
       
  1255 "/    ((mantissa / (1 bitShift:exponent negated)) asFixedPoint:6) printOn:aStream.
       
  1256 
  2888 
  1257     "Created: / 11-06-2019 / 00:13:00 / Claus Gittinger"
  2889     "Created: / 11-06-2019 / 00:13:00 / Claus Gittinger"
  1258 ! !
  2890 ! !
  1259 
  2891 
  1260 !QuadFloat methodsFor:'queries'!
  2892 !QuadFloat methodsFor:'queries'!
  1265 %{  /* NOCONTEXT */
  2897 %{  /* NOCONTEXT */
  1266 #ifdef SUPPORT_QUADFLOAT
  2898 #ifdef SUPPORT_QUADFLOAT
  1267     float128_t myVal;
  2899     float128_t myVal;
  1268 
  2900 
  1269     myVal = __quadFloatVal(self);
  2901     myVal = __quadFloatVal(self);
  1270     RETURN (STX_qisfinite(myVal) ? true : false);
  2902     RETURN ((f128_isInf(myVal) || f128_isNan(myVal)) ? false : true);
  1271 #endif // SUPPORT_QUADFLOAT
  2903 #endif // SUPPORT_QUADFLOAT
  1272 %}.
  2904 %}.
  1273 
  2905 
  1274     "
  2906     "
  1275 	1.0 isFinite
  2907         1.0 isFinite
  1276 	self NaN isFinite
  2908         self NaN isFinite
  1277 	self infinity isFinite
  2909         self infinity isFinite
  1278 	self negativeInfinity isFinite
  2910         self negativeInfinity isFinite
  1279 	(0.0 uncheckedDivide: 0.0) isFinite
  2911         (0.0 uncheckedDivide: 0.0) isFinite
  1280 	(1.0 uncheckedDivide: 0.0) isFinite
  2912         (1.0 uncheckedDivide: 0.0) isFinite
  1281     "
  2913     "
  1282 !
  2914 !
  1283 
  2915 
  1284 isInfinite
  2916 isInfinite
  1285     "return true, if the receiver is an infinite float (+Inf or -Inf)."
  2917     "return true, if the receiver is an infinite float (+Inf or -Inf)."
  1287 %{  /* NOCONTEXT */
  2919 %{  /* NOCONTEXT */
  1288 #ifdef SUPPORT_QUADFLOAT
  2920 #ifdef SUPPORT_QUADFLOAT
  1289     float128_t myVal;
  2921     float128_t myVal;
  1290 
  2922 
  1291     myVal = __quadFloatVal(self);
  2923     myVal = __quadFloatVal(self);
  1292     RETURN (STX_qisInf(&myVal) ? true : false);
  2924     RETURN (f128_isInf(myVal) ? true : false);
  1293 #endif // SUPPORT_QUADFLOAT
  2925 #endif // SUPPORT_QUADFLOAT
  1294 %}.
  2926 %}.
  1295 
  2927 
  1296     "
  2928     "
  1297 	1.0 asQuadFloat isFinite
  2929         1.0 isFinite
  1298 	1.0 asQuadFloat isInfinite
  2930         1.0 isInfinite
  1299 	self NaN isFinite
  2931         self NaN isFinite
  1300 	self NaN isInfinite
  2932         self NaN isInfinite
  1301 	self infinity isFinite
  2933         self infinity isFinite
  1302 	self infinity isInfinite
  2934         self infinity isInfinite
  1303 	self negativeInfinity isFinite
  2935         self negativeInfinity isFinite
  1304 	self negativeInfinity isInfinite
  2936         self negativeInfinity isInfinite
  1305 	(0.0 uncheckedDivide: 0.0) isFinite
  2937         (0.0 uncheckedDivide: 0.0) isFinite
  1306 	(1.0 uncheckedDivide: 0.0) isFinite
  2938         (1.0 uncheckedDivide: 0.0) isFinite
  1307     "
  2939     "
  1308 !
  2940 !
  1309 
  2941 
  1310 isNaN
  2942 isNaN
  1311     "return true, if the receiver is an invalid float (NaN - not a number).
  2943     "return true, if the receiver is an invalid float (NaN - not a number).
  1315 %{  /* NOCONTEXT */
  2947 %{  /* NOCONTEXT */
  1316 #ifdef SUPPORT_QUADFLOAT
  2948 #ifdef SUPPORT_QUADFLOAT
  1317     float128_t myVal;
  2949     float128_t myVal;
  1318 
  2950 
  1319     myVal = __quadFloatVal(self);
  2951     myVal = __quadFloatVal(self);
  1320     RETURN (STX_qisNan(&myVal) ? true : false);
  2952     RETURN (f128_isNan(myVal) ? true : false);
  1321 #endif // SUPPORT_QUADFLOAT
  2953 #endif // SUPPORT_QUADFLOAT
  1322 %}.
  2954 %}.
  1323 
  2955 
  1324     "
  2956     "
  1325 	1.0 asQuadFloat isFinite
  2957         1.0 isFinite
  1326 	self NaN isFinite
  2958         self NaN isFinite
  1327 	self infinity isFinite
  2959         self infinity isFinite
  1328 	(0.0 uncheckedDivide: 0.0) isFinite
  2960         (0.0 uncheckedDivide: 0.0) isFinite
  1329 	(1.0 uncheckedDivide: 0.0) isFinite
  2961         (1.0 uncheckedDivide: 0.0) isFinite
  1330     "
  2962     "
  1331 ! !
       
  1332 
       
  1333 !QuadFloat methodsFor:'testing'!
       
  1334 
       
  1335 isQuadFloat
       
  1336     "return true, if the receiver is some kind of quad floating point number (iee quad precision)"
       
  1337 
       
  1338     ^ true
       
  1339 ! !
  2963 ! !
  1340 
  2964 
  1341 !QuadFloat class methodsFor:'documentation'!
  2965 !QuadFloat class methodsFor:'documentation'!
  1342 
  2966 
  1343 version_CVS
  2967 version_CVS
  1344     ^ '$Header$'
  2968     ^ '$Header$'
  1345 ! !
  2969 ! !
  1346 
  2970 
  1347 
       
  1348 QuadFloat initialize!