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