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