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