1 /*============================================================================
2 This source file is an extension to the SoftFloat IEC/IEEE Floating-point
3 Arithmetic Package, Release 2b, written for Bochs (x86 achitecture simulator)
4 floating point emulation.
5 
6 THIS SOFTWARE IS DISTRIBUTED AS IS, FOR FREE.  Although reasonable effort has
7 been made to avoid it, THIS SOFTWARE MAY CONTAIN FAULTS THAT WILL AT TIMES
8 RESULT IN INCORRECT BEHAVIOR.  USE OF THIS SOFTWARE IS RESTRICTED TO PERSONS
9 AND ORGANIZATIONS WHO CAN AND WILL TAKE FULL RESPONSIBILITY FOR ALL LOSSES,
10 COSTS, OR OTHER PROBLEMS THEY INCUR DUE TO THE SOFTWARE, AND WHO FURTHERMORE
11 EFFECTIVELY INDEMNIFY JOHN HAUSER AND THE INTERNATIONAL COMPUTER SCIENCE
12 INSTITUTE (possibly via similar legal warning) AGAINST ALL LOSSES, COSTS, OR
13 OTHER PROBLEMS INCURRED BY THEIR CUSTOMERS AND CLIENTS DUE TO THE SOFTWARE.
14 
15 Derivative works are acceptable, even for commercial purposes, so long as
16 (1) the source code for the derivative work includes prominent notice that
17 the work is derivative, and (2) the source code includes prominent notice with
18 these four paragraphs for those parts of this code that are retained.
19 =============================================================================*/
20 
21 /*============================================================================
22  * Written for Bochs (x86 achitecture simulator) by
23  *            Stanislav Shwartsman [sshwarts at sourceforge net]
24  * ==========================================================================*/
25 
26 #include "softfloatx80.h"
27 #include "softfloat-round-pack.h"
28 #define USE_estimateDiv128To64
29 #include "softfloat-macros.h"
30 
31 /* executes single exponent reduction cycle */
remainder_kernel(Bit64u aSig0,Bit64u bSig,int expDiff,Bit64u * zSig0,Bit64u * zSig1)32 static Bit64u remainder_kernel(Bit64u aSig0, Bit64u bSig, int expDiff, Bit64u *zSig0, Bit64u *zSig1)
33 {
34     Bit64u term0, term1;
35     Bit64u aSig1 = 0;
36 
37     shortShift128Left(aSig1, aSig0, expDiff, &aSig1, &aSig0);
38     Bit64u q = estimateDiv128To64(aSig1, aSig0, bSig);
39     mul64To128(bSig, q, &term0, &term1);
40     sub128(aSig1, aSig0, term0, term1, zSig1, zSig0);
41     while ((Bit64s)(*zSig1) < 0) {
42         --q;
43         add128(*zSig1, *zSig0, 0, bSig, zSig1, zSig0);
44     }
45     return q;
46 }
47 
do_fprem(floatx80 a,floatx80 b,floatx80 & r,Bit64u & q,int rounding_mode,float_status_t & status)48 static int do_fprem(floatx80 a, floatx80 b, floatx80 &r, Bit64u &q, int rounding_mode, float_status_t &status)
49 {
50     Bit32s aExp, bExp, zExp, expDiff;
51     Bit64u aSig0, aSig1, bSig;
52     int aSign;
53     q = 0;
54 
55     // handle unsupported extended double-precision floating encodings
56     if (floatx80_is_unsupported(a) || floatx80_is_unsupported(b))
57     {
58         float_raise(status, float_flag_invalid);
59         r = floatx80_default_nan;
60         return -1;
61     }
62 
63     aSig0 = extractFloatx80Frac(a);
64     aExp = extractFloatx80Exp(a);
65     aSign = extractFloatx80Sign(a);
66     bSig = extractFloatx80Frac(b);
67     bExp = extractFloatx80Exp(b);
68 
69     if (aExp == 0x7FFF) {
70         if ((Bit64u) (aSig0<<1) || ((bExp == 0x7FFF) && (Bit64u) (bSig<<1))) {
71             r = propagateFloatx80NaN(a, b, status);
72             return -1;
73         }
74         float_raise(status, float_flag_invalid);
75         r = floatx80_default_nan;
76         return -1;
77     }
78     if (bExp == 0x7FFF) {
79         if ((Bit64u) (bSig<<1)) {
80             r = propagateFloatx80NaN(a, b, status);
81             return -1;
82         }
83         if (aExp == 0 && aSig0) {
84             float_raise(status, float_flag_denormal);
85             normalizeFloatx80Subnormal(aSig0, &aExp, &aSig0);
86             r = (a.fraction & BX_CONST64(0x8000000000000000)) ?
87                     packFloatx80(aSign, aExp, aSig0) : a;
88             return 0;
89         }
90         r = a;
91         return 0;
92 
93     }
94     if (bExp == 0) {
95         if (bSig == 0) {
96             float_raise(status, float_flag_invalid);
97             r = floatx80_default_nan;
98             return -1;
99         }
100         float_raise(status, float_flag_denormal);
101         normalizeFloatx80Subnormal(bSig, &bExp, &bSig);
102     }
103     if (aExp == 0) {
104         if (aSig0 == 0) {
105             r = a;
106             return 0;
107         }
108         float_raise(status, float_flag_denormal);
109         normalizeFloatx80Subnormal(aSig0, &aExp, &aSig0);
110     }
111     expDiff = aExp - bExp;
112     aSig1 = 0;
113 
114     Bit32u overflow = 0;
115 
116     if (expDiff >= 64) {
117         int n = (expDiff & 0x1f) | 0x20;
118         remainder_kernel(aSig0, bSig, n, &aSig0, &aSig1);
119         zExp = aExp - n;
120         overflow = 1;
121     }
122     else {
123         zExp = bExp;
124 
125         if (expDiff < 0) {
126             if (expDiff < -1) {
127                r = (a.fraction & BX_CONST64(0x8000000000000000)) ?
128                     packFloatx80(aSign, aExp, aSig0) : a;
129                return 0;
130             }
131             shift128Right(aSig0, 0, 1, &aSig0, &aSig1);
132             expDiff = 0;
133         }
134 
135         if (expDiff > 0) {
136             q = remainder_kernel(aSig0, bSig, expDiff, &aSig0, &aSig1);
137         }
138         else {
139             if (bSig <= aSig0) {
140                aSig0 -= bSig;
141                q = 1;
142             }
143         }
144 
145         if (rounding_mode == float_round_nearest_even)
146         {
147             Bit64u term0, term1;
148             shift128Right(bSig, 0, 1, &term0, &term1);
149 
150             if (! lt128(aSig0, aSig1, term0, term1))
151             {
152                int lt = lt128(term0, term1, aSig0, aSig1);
153                int eq = eq128(aSig0, aSig1, term0, term1);
154 
155                if ((eq && (q & 1)) || lt) {
156                   aSign = !aSign;
157                   ++q;
158                }
159                if (lt) sub128(bSig, 0, aSig0, aSig1, &aSig0, &aSig1);
160             }
161         }
162     }
163 
164     r = normalizeRoundAndPackFloatx80(80, aSign, zExp, aSig0, aSig1, status);
165     return overflow;
166 }
167 
168 /*----------------------------------------------------------------------------
169 | Returns the remainder of the extended double-precision floating-point value
170 | `a' with respect to the corresponding value `b'.  The operation is performed
171 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
172 *----------------------------------------------------------------------------*/
173 
floatx80_ieee754_remainder(floatx80 a,floatx80 b,floatx80 & r,Bit64u & q,float_status_t & status)174 int floatx80_ieee754_remainder(floatx80 a, floatx80 b, floatx80 &r, Bit64u &q, float_status_t &status)
175 {
176     return do_fprem(a, b, r, q, float_round_nearest_even, status);
177 }
178 
179 /*----------------------------------------------------------------------------
180 | Returns the remainder of the extended double-precision floating-point value
181 | `a' with  respect to  the corresponding value `b'. Unlike previous function
182 | the  function  does not compute  the remainder  specified  in  the IEC/IEEE
183 | Standard  for Binary  Floating-Point  Arithmetic.  This  function  operates
184 | differently  from the  previous  function in  the way  that it  rounds  the
185 | quotient of 'a' divided by 'b' to an integer.
186 *----------------------------------------------------------------------------*/
187 
floatx80_remainder(floatx80 a,floatx80 b,floatx80 & r,Bit64u & q,float_status_t & status)188 int floatx80_remainder(floatx80 a, floatx80 b, floatx80 &r, Bit64u &q, float_status_t &status)
189 {
190     return do_fprem(a, b, r, q, float_round_to_zero, status);
191 }
192