1
2/*============================================================================
3
4This C source file is part of the SoftFloat IEEE Floating-Point Arithmetic
5Package, Release 3e, by John R. Hauser.
6
7Copyright 2011, 2012, 2013, 2014, 2015, 2016 The Regents of the University of
8California.  All rights reserved.
9
10Redistribution and use in source and binary forms, with or without
11modification, are permitted provided that the following conditions are met:
12
13 1. Redistributions of source code must retain the above copyright notice,
14    this list of conditions, and the following disclaimer.
15
16 2. Redistributions in binary form must reproduce the above copyright notice,
17    this list of conditions, and the following disclaimer in the documentation
18    and/or other materials provided with the distribution.
19
20 3. Neither the name of the University nor the names of its contributors may
21    be used to endorse or promote products derived from this software without
22    specific prior written permission.
23
24THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS "AS IS", AND ANY
25EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
26WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, ARE
27DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE FOR ANY
28DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
29(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
30LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
31ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
32(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
33SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
34
35=============================================================================*/
36
37#include <stdbool.h>
38#include <stdint.h>
39#include "platform.h"
40#include "internals.h"
41#include "specialize.h"
42#include "softfloat.h"
43
44#ifdef SOFTFLOAT_FAST_INT64
45
46float64_t
47 softfloat_mulAddF64(
48     uint_fast64_t uiA, uint_fast64_t uiB, uint_fast64_t uiC, uint_fast8_t op )
49{
50    bool signA;
51    int_fast16_t expA;
52    uint_fast64_t sigA;
53    bool signB;
54    int_fast16_t expB;
55    uint_fast64_t sigB;
56    bool signC;
57    int_fast16_t expC;
58    uint_fast64_t sigC;
59    bool signZ;
60    uint_fast64_t magBits, uiZ;
61    struct exp16_sig64 normExpSig;
62    int_fast16_t expZ;
63    struct uint128 sig128Z;
64    uint_fast64_t sigZ;
65    int_fast16_t expDiff;
66    struct uint128 sig128C;
67    int_fast8_t shiftDist;
68    union ui64_f64 uZ;
69
70    /*------------------------------------------------------------------------
71    *------------------------------------------------------------------------*/
72    signA = signF64UI( uiA );
73    expA  = expF64UI( uiA );
74    sigA  = fracF64UI( uiA );
75    signB = signF64UI( uiB );
76    expB  = expF64UI( uiB );
77    sigB  = fracF64UI( uiB );
78    signC = signF64UI( uiC ) ^ (op == softfloat_mulAdd_subC);
79    expC  = expF64UI( uiC );
80    sigC  = fracF64UI( uiC );
81    signZ = signA ^ signB ^ (op == softfloat_mulAdd_subProd);
82    /*------------------------------------------------------------------------
83    *------------------------------------------------------------------------*/
84    if ( expA == 0x7FF ) {
85        if ( sigA || ((expB == 0x7FF) && sigB) ) goto propagateNaN_ABC;
86        magBits = expB | sigB;
87        goto infProdArg;
88    }
89    if ( expB == 0x7FF ) {
90        if ( sigB ) goto propagateNaN_ABC;
91        magBits = expA | sigA;
92        goto infProdArg;
93    }
94    if ( expC == 0x7FF ) {
95        if ( sigC ) {
96            uiZ = 0;
97            goto propagateNaN_ZC;
98        }
99        uiZ = uiC;
100        goto uiZ;
101    }
102    /*------------------------------------------------------------------------
103    *------------------------------------------------------------------------*/
104    if ( ! expA ) {
105        if ( ! sigA ) goto zeroProd;
106        normExpSig = softfloat_normSubnormalF64Sig( sigA );
107        expA = normExpSig.exp;
108        sigA = normExpSig.sig;
109    }
110    if ( ! expB ) {
111        if ( ! sigB ) goto zeroProd;
112        normExpSig = softfloat_normSubnormalF64Sig( sigB );
113        expB = normExpSig.exp;
114        sigB = normExpSig.sig;
115    }
116    /*------------------------------------------------------------------------
117    *------------------------------------------------------------------------*/
118    expZ = expA + expB - 0x3FE;
119    sigA = (sigA | UINT64_C( 0x0010000000000000 ))<<10;
120    sigB = (sigB | UINT64_C( 0x0010000000000000 ))<<10;
121    sig128Z = softfloat_mul64To128( sigA, sigB );
122    if ( sig128Z.v64 < UINT64_C( 0x2000000000000000 ) ) {
123        --expZ;
124        sig128Z =
125            softfloat_add128(
126                sig128Z.v64, sig128Z.v0, sig128Z.v64, sig128Z.v0 );
127    }
128    if ( ! expC ) {
129        if ( ! sigC ) {
130            --expZ;
131            sigZ = sig128Z.v64<<1 | (sig128Z.v0 != 0);
132            goto roundPack;
133        }
134        normExpSig = softfloat_normSubnormalF64Sig( sigC );
135        expC = normExpSig.exp;
136        sigC = normExpSig.sig;
137    }
138    sigC = (sigC | UINT64_C( 0x0010000000000000 ))<<9;
139    /*------------------------------------------------------------------------
140    *------------------------------------------------------------------------*/
141    expDiff = expZ - expC;
142    if ( expDiff < 0 ) {
143        expZ = expC;
144        if ( (signZ == signC) || (expDiff < -1) ) {
145            sig128Z.v64 = softfloat_shiftRightJam64( sig128Z.v64, -expDiff );
146        } else {
147            sig128Z =
148                softfloat_shortShiftRightJam128( sig128Z.v64, sig128Z.v0, 1 );
149        }
150    } else if ( expDiff ) {
151        sig128C = softfloat_shiftRightJam128( sigC, 0, expDiff );
152    }
153    /*------------------------------------------------------------------------
154    *------------------------------------------------------------------------*/
155    if ( signZ == signC ) {
156        /*--------------------------------------------------------------------
157        *--------------------------------------------------------------------*/
158        if ( expDiff <= 0 ) {
159            sigZ = (sigC + sig128Z.v64) | (sig128Z.v0 != 0);
160        } else {
161            sig128Z =
162                softfloat_add128(
163                    sig128Z.v64, sig128Z.v0, sig128C.v64, sig128C.v0 );
164            sigZ = sig128Z.v64 | (sig128Z.v0 != 0);
165        }
166        if ( sigZ < UINT64_C( 0x4000000000000000 ) ) {
167            --expZ;
168            sigZ <<= 1;
169        }
170    } else {
171        /*--------------------------------------------------------------------
172        *--------------------------------------------------------------------*/
173        if ( expDiff < 0 ) {
174            signZ = signC;
175            sig128Z = softfloat_sub128( sigC, 0, sig128Z.v64, sig128Z.v0 );
176        } else if ( ! expDiff ) {
177            sig128Z.v64 = sig128Z.v64 - sigC;
178            if ( ! (sig128Z.v64 | sig128Z.v0) ) goto completeCancellation;
179            if ( sig128Z.v64 & UINT64_C( 0x8000000000000000 ) ) {
180                signZ = ! signZ;
181                sig128Z = softfloat_sub128( 0, 0, sig128Z.v64, sig128Z.v0 );
182            }
183        } else {
184            sig128Z =
185                softfloat_sub128(
186                    sig128Z.v64, sig128Z.v0, sig128C.v64, sig128C.v0 );
187        }
188        /*--------------------------------------------------------------------
189        *--------------------------------------------------------------------*/
190        if ( ! sig128Z.v64 ) {
191            expZ -= 64;
192            sig128Z.v64 = sig128Z.v0;
193            sig128Z.v0 = 0;
194        }
195        shiftDist = softfloat_countLeadingZeros64( sig128Z.v64 ) - 1;
196        expZ -= shiftDist;
197        if ( shiftDist < 0 ) {
198            sigZ = softfloat_shortShiftRightJam64( sig128Z.v64, -shiftDist );
199        } else {
200            sig128Z =
201                softfloat_shortShiftLeft128(
202                    sig128Z.v64, sig128Z.v0, shiftDist );
203            sigZ = sig128Z.v64;
204        }
205        sigZ |= (sig128Z.v0 != 0);
206    }
207 roundPack:
208    return softfloat_roundPackToF64( signZ, expZ, sigZ );
209    /*------------------------------------------------------------------------
210    *------------------------------------------------------------------------*/
211 propagateNaN_ABC:
212    uiZ = softfloat_propagateNaNF64UI( uiA, uiB );
213    goto propagateNaN_ZC;
214    /*------------------------------------------------------------------------
215    *------------------------------------------------------------------------*/
216 infProdArg:
217    if ( magBits ) {
218        uiZ = packToF64UI( signZ, 0x7FF, 0 );
219        if ( expC != 0x7FF ) goto uiZ;
220        if ( sigC ) goto propagateNaN_ZC;
221        if ( signZ == signC ) goto uiZ;
222    }
223    softfloat_raiseFlags( softfloat_flag_invalid );
224    uiZ = defaultNaNF64UI;
225 propagateNaN_ZC:
226    uiZ = softfloat_propagateNaNF64UI( uiZ, uiC );
227    goto uiZ;
228    /*------------------------------------------------------------------------
229    *------------------------------------------------------------------------*/
230 zeroProd:
231    uiZ = uiC;
232    if ( ! (expC | sigC) && (signZ != signC) ) {
233 completeCancellation:
234        uiZ =
235            packToF64UI(
236                (softfloat_roundingMode == softfloat_round_min), 0, 0 );
237    }
238 uiZ:
239    uZ.ui = uiZ;
240    return uZ.f;
241
242}
243
244#else
245
246float64_t
247 softfloat_mulAddF64(
248     uint_fast64_t uiA, uint_fast64_t uiB, uint_fast64_t uiC, uint_fast8_t op )
249{
250    bool signA;
251    int_fast16_t expA;
252    uint64_t sigA;
253    bool signB;
254    int_fast16_t expB;
255    uint64_t sigB;
256    bool signC;
257    int_fast16_t expC;
258    uint64_t sigC;
259    bool signZ;
260    uint64_t magBits, uiZ;
261    struct exp16_sig64 normExpSig;
262    int_fast16_t expZ;
263    uint32_t sig128Z[4];
264    uint64_t sigZ;
265    int_fast16_t shiftDist, expDiff;
266    uint32_t sig128C[4];
267    union ui64_f64 uZ;
268
269    /*------------------------------------------------------------------------
270    *------------------------------------------------------------------------*/
271    signA = signF64UI( uiA );
272    expA  = expF64UI( uiA );
273    sigA  = fracF64UI( uiA );
274    signB = signF64UI( uiB );
275    expB  = expF64UI( uiB );
276    sigB  = fracF64UI( uiB );
277    signC = signF64UI( uiC ) ^ (op == softfloat_mulAdd_subC);
278    expC  = expF64UI( uiC );
279    sigC  = fracF64UI( uiC );
280    signZ = signA ^ signB ^ (op == softfloat_mulAdd_subProd);
281    /*------------------------------------------------------------------------
282    *------------------------------------------------------------------------*/
283    if ( expA == 0x7FF ) {
284        if ( sigA || ((expB == 0x7FF) && sigB) ) goto propagateNaN_ABC;
285        magBits = expB | sigB;
286        goto infProdArg;
287    }
288    if ( expB == 0x7FF ) {
289        if ( sigB ) goto propagateNaN_ABC;
290        magBits = expA | sigA;
291        goto infProdArg;
292    }
293    if ( expC == 0x7FF ) {
294        if ( sigC ) {
295            uiZ = 0;
296            goto propagateNaN_ZC;
297        }
298        uiZ = uiC;
299        goto uiZ;
300    }
301    /*------------------------------------------------------------------------
302    *------------------------------------------------------------------------*/
303    if ( ! expA ) {
304        if ( ! sigA ) goto zeroProd;
305        normExpSig = softfloat_normSubnormalF64Sig( sigA );
306        expA = normExpSig.exp;
307        sigA = normExpSig.sig;
308    }
309    if ( ! expB ) {
310        if ( ! sigB ) goto zeroProd;
311        normExpSig = softfloat_normSubnormalF64Sig( sigB );
312        expB = normExpSig.exp;
313        sigB = normExpSig.sig;
314    }
315    /*------------------------------------------------------------------------
316    *------------------------------------------------------------------------*/
317    expZ = expA + expB - 0x3FE;
318    sigA = (sigA | UINT64_C( 0x0010000000000000 ))<<10;
319    sigB = (sigB | UINT64_C( 0x0010000000000000 ))<<11;
320    softfloat_mul64To128M( sigA, sigB, sig128Z );
321    sigZ =
322        (uint64_t) sig128Z[indexWord( 4, 3 )]<<32 | sig128Z[indexWord( 4, 2 )];
323    shiftDist = 0;
324    if ( ! (sigZ & UINT64_C( 0x4000000000000000 )) ) {
325        --expZ;
326        shiftDist = -1;
327    }
328    if ( ! expC ) {
329        if ( ! sigC ) {
330            if ( shiftDist ) sigZ <<= 1;
331            goto sigZ;
332        }
333        normExpSig = softfloat_normSubnormalF64Sig( sigC );
334        expC = normExpSig.exp;
335        sigC = normExpSig.sig;
336    }
337    sigC = (sigC | UINT64_C( 0x0010000000000000 ))<<10;
338    /*------------------------------------------------------------------------
339    *------------------------------------------------------------------------*/
340    expDiff = expZ - expC;
341    if ( expDiff < 0 ) {
342        expZ = expC;
343        if ( (signZ == signC) || (expDiff < -1) ) {
344            shiftDist -= expDiff;
345            if ( shiftDist) {
346                sigZ = softfloat_shiftRightJam64( sigZ, shiftDist );
347            }
348        } else {
349            if ( ! shiftDist ) {
350                softfloat_shortShiftRight128M( sig128Z, 1, sig128Z );
351            }
352        }
353    } else {
354        if ( shiftDist ) softfloat_add128M( sig128Z, sig128Z, sig128Z );
355        if ( ! expDiff ) {
356            sigZ =
357                (uint64_t) sig128Z[indexWord( 4, 3 )]<<32
358                    | sig128Z[indexWord( 4, 2 )];
359        } else {
360            sig128C[indexWord( 4, 3 )] = sigC>>32;
361            sig128C[indexWord( 4, 2 )] = sigC;
362            sig128C[indexWord( 4, 1 )] = 0;
363            sig128C[indexWord( 4, 0 )] = 0;
364            softfloat_shiftRightJam128M( sig128C, expDiff, sig128C );
365        }
366    }
367    /*------------------------------------------------------------------------
368    *------------------------------------------------------------------------*/
369    if ( signZ == signC ) {
370        /*--------------------------------------------------------------------
371        *--------------------------------------------------------------------*/
372        if ( expDiff <= 0 ) {
373            sigZ += sigC;
374        } else {
375            softfloat_add128M( sig128Z, sig128C, sig128Z );
376            sigZ =
377                (uint64_t) sig128Z[indexWord( 4, 3 )]<<32
378                    | sig128Z[indexWord( 4, 2 )];
379        }
380        if ( sigZ & UINT64_C( 0x8000000000000000 ) ) {
381            ++expZ;
382            sigZ = softfloat_shortShiftRightJam64( sigZ, 1 );
383        }
384    } else {
385        /*--------------------------------------------------------------------
386        *--------------------------------------------------------------------*/
387        if ( expDiff < 0 ) {
388            signZ = signC;
389            if ( expDiff < -1 ) {
390                sigZ = sigC - sigZ;
391                if (
392                    sig128Z[indexWord( 4, 1 )] || sig128Z[indexWord( 4, 0 )]
393                ) {
394                    sigZ = (sigZ - 1) | 1;
395                }
396                if ( ! (sigZ & UINT64_C( 0x4000000000000000 )) ) {
397                    --expZ;
398                    sigZ <<= 1;
399                }
400                goto roundPack;
401            } else {
402                sig128C[indexWord( 4, 3 )] = sigC>>32;
403                sig128C[indexWord( 4, 2 )] = sigC;
404                sig128C[indexWord( 4, 1 )] = 0;
405                sig128C[indexWord( 4, 0 )] = 0;
406                softfloat_sub128M( sig128C, sig128Z, sig128Z );
407            }
408        } else if ( ! expDiff ) {
409            sigZ -= sigC;
410            if (
411                ! sigZ && ! sig128Z[indexWord( 4, 1 )]
412                    && ! sig128Z[indexWord( 4, 0 )]
413            ) {
414                goto completeCancellation;
415            }
416            sig128Z[indexWord( 4, 3 )] = sigZ>>32;
417            sig128Z[indexWord( 4, 2 )] = sigZ;
418            if ( sigZ & UINT64_C( 0x8000000000000000 ) ) {
419                signZ = ! signZ;
420                softfloat_negX128M( sig128Z );
421            }
422        } else {
423            softfloat_sub128M( sig128Z, sig128C, sig128Z );
424            if ( 1 < expDiff ) {
425                sigZ =
426                    (uint64_t) sig128Z[indexWord( 4, 3 )]<<32
427                        | sig128Z[indexWord( 4, 2 )];
428                if ( ! (sigZ & UINT64_C( 0x4000000000000000 )) ) {
429                    --expZ;
430                    sigZ <<= 1;
431                }
432                goto sigZ;
433            }
434        }
435        /*--------------------------------------------------------------------
436        *--------------------------------------------------------------------*/
437        shiftDist = 0;
438        sigZ =
439            (uint64_t) sig128Z[indexWord( 4, 3 )]<<32
440                | sig128Z[indexWord( 4, 2 )];
441        if ( ! sigZ ) {
442            shiftDist = 64;
443            sigZ =
444                (uint64_t) sig128Z[indexWord( 4, 1 )]<<32
445                    | sig128Z[indexWord( 4, 0 )];
446        }
447        shiftDist += softfloat_countLeadingZeros64( sigZ ) - 1;
448        if ( shiftDist ) {
449            expZ -= shiftDist;
450            softfloat_shiftLeft128M( sig128Z, shiftDist, sig128Z );
451            sigZ =
452                (uint64_t) sig128Z[indexWord( 4, 3 )]<<32
453                    | sig128Z[indexWord( 4, 2 )];
454        }
455    }
456 sigZ:
457    if ( sig128Z[indexWord( 4, 1 )] || sig128Z[indexWord( 4, 0 )] ) sigZ |= 1;
458 roundPack:
459    return softfloat_roundPackToF64( signZ, expZ - 1, sigZ );
460    /*------------------------------------------------------------------------
461    *------------------------------------------------------------------------*/
462 propagateNaN_ABC:
463    uiZ = softfloat_propagateNaNF64UI( uiA, uiB );
464    goto propagateNaN_ZC;
465    /*------------------------------------------------------------------------
466    *------------------------------------------------------------------------*/
467 infProdArg:
468    if ( magBits ) {
469        uiZ = packToF64UI( signZ, 0x7FF, 0 );
470        if ( expC != 0x7FF ) goto uiZ;
471        if ( sigC ) goto propagateNaN_ZC;
472        if ( signZ == signC ) goto uiZ;
473    }
474    softfloat_raiseFlags( softfloat_flag_invalid );
475    uiZ = defaultNaNF64UI;
476 propagateNaN_ZC:
477    uiZ = softfloat_propagateNaNF64UI( uiZ, uiC );
478    goto uiZ;
479    /*------------------------------------------------------------------------
480    *------------------------------------------------------------------------*/
481 zeroProd:
482    uiZ = uiC;
483    if ( ! (expC | sigC) && (signZ != signC) ) {
484 completeCancellation:
485        uiZ =
486            packToF64UI(
487                (softfloat_roundingMode == softfloat_round_min), 0, 0 );
488    }
489 uiZ:
490    uZ.ui = uiZ;
491    return uZ.f;
492
493}
494
495#endif
496
497