1 
2 /*============================================================================
3 
4 This C source file is part of the SoftFloat IEEE Floating-Point Arithmetic
5 Package, Release 3e, by John R. Hauser.
6 
7 Copyright 2011, 2012, 2013, 2014, 2015, 2016 The Regents of the University of
8 California.  All rights reserved.
9 
10 Redistribution and use in source and binary forms, with or without
11 modification, 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 
24 THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS "AS IS", AND ANY
25 EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
26 WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, ARE
27 DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE FOR ANY
28 DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
29 (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
30 LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
31 ON 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
33 SOFTWARE, 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 
46 float64_t
softfloat_mulAddF64(uint_fast64_t uiA,uint_fast64_t uiB,uint_fast64_t uiC,uint_fast8_t op)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 
246 float64_t
softfloat_mulAddF64(uint_fast64_t uiA,uint_fast64_t uiB,uint_fast64_t uiC,uint_fast8_t op)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