1
2 /*============================================================================
3
4 This C source file is part of the SoftFloat IEEE Floating-Point Arithmetic
5 Package, Release 3a, by John R. Hauser.
6
7 Copyright 2011, 2012, 2013, 2014 The Regents of the University of California.
8 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
extF80_rem(extFloat80_t a,extFloat80_t b)44 extFloat80_t extF80_rem( extFloat80_t a, extFloat80_t b )
45 {
46 union { struct extFloat80M s; extFloat80_t f; } uA;
47 uint_fast16_t uiA64;
48 uint_fast64_t uiA0;
49 bool signA;
50 int_fast32_t expA;
51 uint_fast64_t sigA;
52 union { struct extFloat80M s; extFloat80_t f; } uB;
53 uint_fast16_t uiB64;
54 uint_fast64_t uiB0;
55 bool signB;
56 int_fast32_t expB;
57 uint_fast64_t sigB;
58 struct exp32_sig64 normExpSig;
59 int_fast32_t expDiff;
60 struct uint128 rem, shiftedSigB;
61 uint_fast32_t q, recip32;
62 uint_fast64_t q64;
63 struct uint128 term, altRem, meanRem;
64 bool signRem;
65 struct uint128 uiZ;
66 uint_fast16_t uiZ64;
67 uint_fast64_t uiZ0;
68 union { struct extFloat80M s; extFloat80_t f; } uZ;
69
70 /*------------------------------------------------------------------------
71 *------------------------------------------------------------------------*/
72 uA.f = a;
73 uiA64 = uA.s.signExp;
74 uiA0 = uA.s.signif;
75 signA = signExtF80UI64( uiA64 );
76 expA = expExtF80UI64( uiA64 );
77 sigA = uiA0;
78 uB.f = b;
79 uiB64 = uB.s.signExp;
80 uiB0 = uB.s.signif;
81 signB = signExtF80UI64( uiB64 );
82 expB = expExtF80UI64( uiB64 );
83 sigB = uiB0;
84 /*------------------------------------------------------------------------
85 *------------------------------------------------------------------------*/
86 if ( expA == 0x7FFF ) {
87 if (
88 (sigA & UINT64_C( 0x7FFFFFFFFFFFFFFF ))
89 || ((expB == 0x7FFF) && (sigB & UINT64_C( 0x7FFFFFFFFFFFFFFF )))
90 ) {
91 goto propagateNaN;
92 }
93 goto invalid;
94 }
95 if ( expB == 0x7FFF ) {
96 if ( sigB & UINT64_C( 0x7FFFFFFFFFFFFFFF ) ) goto propagateNaN;
97 /*--------------------------------------------------------------------
98 | Argument b is an infinity. Doubling `expB' is an easy way to ensure
99 | that `expDiff' later is less than -1, which will result in returning
100 | a canonicalized version of argument a.
101 *--------------------------------------------------------------------*/
102 expB += expB;
103 }
104 /*------------------------------------------------------------------------
105 *------------------------------------------------------------------------*/
106 if ( ! expB ) expB = 1;
107 if ( ! (sigB & UINT64_C( 0x8000000000000000 )) ) {
108 if ( ! sigB ) goto invalid;
109 normExpSig = softfloat_normSubnormalExtF80Sig( sigB );
110 expB += normExpSig.exp;
111 sigB = normExpSig.sig;
112 }
113 if ( ! expA ) expA = 1;
114 if ( ! (sigA & UINT64_C( 0x8000000000000000 )) ) {
115 if ( ! sigA ) {
116 expA = 0;
117 goto copyA;
118 }
119 normExpSig = softfloat_normSubnormalExtF80Sig( sigA );
120 expA += normExpSig.exp;
121 sigA = normExpSig.sig;
122 }
123 /*------------------------------------------------------------------------
124 *------------------------------------------------------------------------*/
125 expDiff = expA - expB;
126 if ( expDiff < -1 ) goto copyA;
127 rem = softfloat_shortShiftLeft128( 0, sigA, 32 );
128 shiftedSigB = softfloat_shortShiftLeft128( 0, sigB, 32 );
129 if ( expDiff < 1 ) {
130 if ( expDiff ) {
131 --expB;
132 shiftedSigB = softfloat_shortShiftLeft128( 0, sigB, 33 );
133 q = 0;
134 } else {
135 q = (sigB <= sigA);
136 if ( q ) {
137 rem =
138 softfloat_sub128(
139 rem.v64, rem.v0, shiftedSigB.v64, shiftedSigB.v0 );
140 }
141 }
142 } else {
143 recip32 = softfloat_approxRecip32_1( sigB>>32 );
144 expDiff -= 30;
145 for (;;) {
146 q64 = (uint_fast64_t) (uint32_t) (rem.v64>>2) * recip32;
147 if ( expDiff < 0 ) break;
148 q = (q64 + 0x80000000)>>32;
149 rem = softfloat_shortShiftLeft128( rem.v64, rem.v0, 29 );
150 term = softfloat_mul64ByShifted32To128( sigB, q );
151 rem = softfloat_sub128( rem.v64, rem.v0, term.v64, term.v0 );
152 if ( rem.v64 & UINT64_C( 0x8000000000000000 ) ) {
153 rem =
154 softfloat_add128(
155 rem.v64, rem.v0, shiftedSigB.v64, shiftedSigB.v0 );
156 }
157 expDiff -= 29;
158 }
159 /*--------------------------------------------------------------------
160 | (`expDiff' cannot be less than -29 here.)
161 *--------------------------------------------------------------------*/
162 q = (uint32_t) (q64>>32)>>(~expDiff & 31);
163 rem = softfloat_shortShiftLeft128( rem.v64, rem.v0, expDiff + 30 );
164 term = softfloat_mul64ByShifted32To128( sigB, q );
165 rem = softfloat_sub128( rem.v64, rem.v0, term.v64, term.v0 );
166 if ( rem.v64 & UINT64_C( 0x8000000000000000 ) ) {
167 altRem =
168 softfloat_add128(
169 rem.v64, rem.v0, shiftedSigB.v64, shiftedSigB.v0 );
170 goto selectRem;
171 }
172 }
173 /*------------------------------------------------------------------------
174 *------------------------------------------------------------------------*/
175 do {
176 altRem = rem;
177 ++q;
178 rem =
179 softfloat_sub128(
180 rem.v64, rem.v0, shiftedSigB.v64, shiftedSigB.v0 );
181 } while ( ! (rem.v64 & UINT64_C( 0x8000000000000000 )) );
182 selectRem:
183 meanRem = softfloat_add128( rem.v64, rem.v0, altRem.v64, altRem.v0 );
184 if (
185 (meanRem.v64 & UINT64_C( 0x8000000000000000 ))
186 || (! (meanRem.v64 | meanRem.v0) && (q & 1))
187 ) {
188 rem = altRem;
189 }
190 signRem = signA;
191 if ( rem.v64 & UINT64_C( 0x8000000000000000 ) ) {
192 signRem = ! signRem;
193 rem = softfloat_sub128( 0, 0, rem.v64, rem.v0 );
194 }
195 return
196 softfloat_normRoundPackToExtF80(
197 signRem, expB + 32, rem.v64, rem.v0, 80 );
198 /*------------------------------------------------------------------------
199 *------------------------------------------------------------------------*/
200 propagateNaN:
201 uiZ = softfloat_propagateNaNExtF80UI( uiA64, uiA0, uiB64, uiB0 );
202 uiZ64 = uiZ.v64;
203 uiZ0 = uiZ.v0;
204 goto uiZ;
205 /*------------------------------------------------------------------------
206 *------------------------------------------------------------------------*/
207 invalid:
208 softfloat_raiseFlags( softfloat_flag_invalid );
209 uiZ64 = defaultNaNExtF80UI64;
210 uiZ0 = defaultNaNExtF80UI0;
211 goto uiZ;
212 /*------------------------------------------------------------------------
213 *------------------------------------------------------------------------*/
214 copyA:
215 if ( expA < 1 ) {
216 sigA >>= 1 - expA;
217 expA = 0;
218 }
219 uiZ64 = packToExtF80UI64( signA, expA );
220 uiZ0 = sigA;
221 uiZ:
222 uZ.s.signExp = uiZ64;
223 uZ.s.signif = uiZ0;
224 return uZ.f;
225
226 }
227
228