1 /*
2 * CDDL HEADER START
3 *
4 * The contents of this file are subject to the terms of the
5 * Common Development and Distribution License, Version 1.0 only
6 * (the "License"). You may not use this file except in compliance
7 * with the License.
8 *
9 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
10 * or http://www.opensolaris.org/os/licensing.
11 * See the License for the specific language governing permissions
12 * and limitations under the License.
13 *
14 * When distributing Covered Code, include this CDDL HEADER in each
15 * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
16 * If applicable, add the following below this CDDL HEADER, with the
17 * fields enclosed by brackets "[]" replaced with your own identifying
18 * information: Portions Copyright [yyyy] [name of copyright owner]
19 *
20 * CDDL HEADER END
21 */
22 /*
23 * Copyright (c) 1994-1997, by Sun Microsystems, Inc.
24 * All rights reserved.
25 */
26
27 /*
28 * This file contains __quad_mag_add and __quad_mag_sub, the core
29 * of the quad precision add and subtract operations.
30 */
31
32 #include "quad.h"
33
34 /*
35 * __quad_mag_add(x, y, z, fsr)
36 *
37 * Sets *z = *x + *y, rounded according to the rounding mode in *fsr,
38 * and updates the current exceptions in *fsr. This routine assumes
39 * *x and *y are finite, with the same sign (i.e., an addition of
40 * magnitudes), |*x| >= |*y|, and *z already has its sign bit set.
41 */
42 void
__quad_mag_add(const union longdouble * x,const union longdouble * y,union longdouble * z,unsigned int * fsr)43 __quad_mag_add(const union longdouble *x, const union longdouble *y,
44 union longdouble *z, unsigned int *fsr)
45 {
46 unsigned int lx, ly, ex, ey, frac2, frac3, frac4;
47 unsigned int round, sticky, carry, rm;
48 int e, uflo;
49
50 /* get the leading significand words and exponents */
51 ex = (x->l.msw & 0x7fffffff) >> 16;
52 lx = x->l.msw & 0xffff;
53 if (ex == 0)
54 ex = 1;
55 else
56 lx |= 0x10000;
57
58 ey = (y->l.msw & 0x7fffffff) >> 16;
59 ly = y->l.msw & 0xffff;
60 if (ey == 0)
61 ey = 1;
62 else
63 ly |= 0x10000;
64
65 /* prenormalize y */
66 e = (int) ex - (int) ey;
67 round = sticky = 0;
68 if (e >= 114) {
69 frac2 = x->l.frac2;
70 frac3 = x->l.frac3;
71 frac4 = x->l.frac4;
72 sticky = ly | y->l.frac2 | y->l.frac3 | y->l.frac4;
73 } else {
74 frac2 = y->l.frac2;
75 frac3 = y->l.frac3;
76 frac4 = y->l.frac4;
77 if (e >= 96) {
78 sticky = frac4 | frac3 | (frac2 & 0x7fffffff);
79 round = frac2 & 0x80000000;
80 frac4 = ly;
81 frac3 = frac2 = ly = 0;
82 e -= 96;
83 } else if (e >= 64) {
84 sticky = frac4 | (frac3 & 0x7fffffff);
85 round = frac3 & 0x80000000;
86 frac4 = frac2;
87 frac3 = ly;
88 frac2 = ly = 0;
89 e -= 64;
90 } else if (e >= 32) {
91 sticky = frac4 & 0x7fffffff;
92 round = frac4 & 0x80000000;
93 frac4 = frac3;
94 frac3 = frac2;
95 frac2 = ly;
96 ly = 0;
97 e -= 32;
98 }
99 if (e) {
100 sticky |= round | (frac4 & ((1 << (e - 1)) - 1));
101 round = frac4 & (1 << (e - 1));
102 frac4 = (frac4 >> e) | (frac3 << (32 - e));
103 frac3 = (frac3 >> e) | (frac2 << (32 - e));
104 frac2 = (frac2 >> e) | (ly << (32 - e));
105 ly >>= e;
106 }
107
108 /* add, propagating carries */
109 frac4 += x->l.frac4;
110 carry = (frac4 < x->l.frac4);
111 frac3 += x->l.frac3;
112 if (carry) {
113 frac3++;
114 carry = (frac3 <= x->l.frac3);
115 } else {
116 carry = (frac3 < x->l.frac3);
117 }
118 frac2 += x->l.frac2;
119 if (carry) {
120 frac2++;
121 carry = (frac2 <= x->l.frac2);
122 } else {
123 carry = (frac2 < x->l.frac2);
124 }
125 lx += ly;
126 if (carry)
127 lx++;
128
129 /* postnormalize */
130 if (lx >= 0x20000) {
131 sticky |= round;
132 round = frac4 & 1;
133 frac4 = (frac4 >> 1) | (frac3 << 31);
134 frac3 = (frac3 >> 1) | (frac2 << 31);
135 frac2 = (frac2 >> 1) | (lx << 31);
136 lx >>= 1;
137 ex++;
138 }
139 }
140
141 /* keep track of whether the result before rounding is tiny */
142 uflo = (lx < 0x10000);
143
144 /* get the rounding mode, fudging directed rounding modes */
145 /* as though the result were positive */
146 rm = *fsr >> 30;
147 if (z->l.msw)
148 rm ^= (rm >> 1);
149
150 /* see if we need to round */
151 if (round | sticky) {
152 *fsr |= FSR_NXC;
153
154 /* round up if necessary */
155 if (rm == FSR_RP || (rm == FSR_RN && round &&
156 (sticky || (frac4 & 1)))) {
157 if (++frac4 == 0)
158 if (++frac3 == 0)
159 if (++frac2 == 0)
160 if (++lx >= 0x20000) {
161 lx >>= 1;
162 ex++;
163 }
164 }
165 }
166
167 /* check for overflow */
168 if (ex >= 0x7fff) {
169 /* store the default overflowed result */
170 *fsr |= FSR_OFC | FSR_NXC;
171 if (rm == FSR_RN || rm == FSR_RP) {
172 z->l.msw |= 0x7fff0000;
173 z->l.frac2 = z->l.frac3 = z->l.frac4 = 0;
174 } else {
175 z->l.msw |= 0x7ffeffff;
176 z->l.frac2 = z->l.frac3 = z->l.frac4 = 0xffffffff;
177 }
178 } else {
179 /* store the result */
180 if (lx >= 0x10000)
181 z->l.msw |= (ex << 16);
182 z->l.msw |= (lx & 0xffff);
183 z->l.frac2 = frac2;
184 z->l.frac3 = frac3;
185 z->l.frac4 = frac4;
186
187 /* if the pre-rounded result was tiny and underflow trapping */
188 /* is enabled, simulate underflow */
189 if (uflo && (*fsr & FSR_UFM))
190 *fsr |= FSR_UFC;
191 }
192 }
193
194 /*
195 * __quad_mag_sub(x, y, z, fsr)
196 *
197 * Sets *z = *x - *y, rounded according to the rounding mode in *fsr,
198 * and updates the current exceptions in *fsr. This routine assumes
199 * *x and *y are finite, with opposite signs (i.e., a subtraction of
200 * magnitudes), |*x| >= |*y|, and *z already has its sign bit set.
201 */
202 void
__quad_mag_sub(const union longdouble * x,const union longdouble * y,union longdouble * z,unsigned int * fsr)203 __quad_mag_sub(const union longdouble *x, const union longdouble *y,
204 union longdouble *z, unsigned int *fsr)
205 {
206 unsigned int lx, ly, ex, ey, frac2, frac3, frac4;
207 unsigned int guard, round, sticky, borrow, rm;
208 int e;
209
210 /* get the leading significand words and exponents */
211 ex = (x->l.msw & 0x7fffffff) >> 16;
212 lx = x->l.msw & 0xffff;
213 if (ex == 0)
214 ex = 1;
215 else
216 lx |= 0x10000;
217
218 ey = (y->l.msw & 0x7fffffff) >> 16;
219 ly = y->l.msw & 0xffff;
220 if (ey == 0)
221 ey = 1;
222 else
223 ly |= 0x10000;
224
225 /* prenormalize y */
226 e = (int) ex - (int) ey;
227 guard = round = sticky = 0;
228 if (e > 114) {
229 sticky = ly | y->l.frac2 | y->l.frac3 | y->l.frac4;
230 ly = frac2 = frac3 = frac4 = 0;
231 } else {
232 frac2 = y->l.frac2;
233 frac3 = y->l.frac3;
234 frac4 = y->l.frac4;
235 if (e >= 96) {
236 sticky = frac4 | frac3 | (frac2 & 0x3fffffff);
237 round = frac2 & 0x40000000;
238 guard = frac2 & 0x80000000;
239 frac4 = ly;
240 frac3 = frac2 = ly = 0;
241 e -= 96;
242 } else if (e >= 64) {
243 sticky = frac4 | (frac3 & 0x3fffffff);
244 round = frac3 & 0x40000000;
245 guard = frac3 & 0x80000000;
246 frac4 = frac2;
247 frac3 = ly;
248 frac2 = ly = 0;
249 e -= 64;
250 } else if (e >= 32) {
251 sticky = frac4 & 0x3fffffff;
252 round = frac4 & 0x40000000;
253 guard = frac4 & 0x80000000;
254 frac4 = frac3;
255 frac3 = frac2;
256 frac2 = ly;
257 ly = 0;
258 e -= 32;
259 }
260 if (e > 1) {
261 sticky |= guard | round |
262 (frac4 & ((1 << (e - 2)) - 1));
263 round = frac4 & (1 << (e - 2));
264 guard = frac4 & (1 << (e - 1));
265 frac4 = (frac4 >> e) | (frac3 << (32 - e));
266 frac3 = (frac3 >> e) | (frac2 << (32 - e));
267 frac2 = (frac2 >> e) | (ly << (32 - e));
268 ly >>= e;
269 } else if (e == 1) {
270 sticky |= round;
271 round = guard;
272 guard = frac4 & 1;
273 frac4 = (frac4 >> 1) | (frac3 << 31);
274 frac3 = (frac3 >> 1) | (frac2 << 31);
275 frac2 = (frac2 >> 1) | (ly << 31);
276 ly >>= 1;
277 }
278 }
279
280 /* complement guard, round, and sticky as need be */
281 if (sticky) {
282 round = !round;
283 guard = !guard;
284 } else if (round) {
285 guard = !guard;
286 }
287 borrow = (guard | round | sticky);
288
289 /* subtract, propagating borrows */
290 frac4 = x->l.frac4 - frac4;
291 if (borrow) {
292 frac4--;
293 borrow = (frac4 >= x->l.frac4);
294 } else {
295 borrow = (frac4 > x->l.frac4);
296 }
297 frac3 = x->l.frac3 - frac3;
298 if (borrow) {
299 frac3--;
300 borrow = (frac3 >= x->l.frac3);
301 } else {
302 borrow = (frac3 > x->l.frac3);
303 }
304 frac2 = x->l.frac2 - frac2;
305 if (borrow) {
306 frac2--;
307 borrow = (frac2 >= x->l.frac2);
308 } else {
309 borrow = (frac2 > x->l.frac2);
310 }
311 lx -= ly;
312 if (borrow)
313 lx--;
314
315 /* get the rounding mode */
316 rm = *fsr >> 30;
317
318 /* handle zero result */
319 if (!(lx | frac2 | frac3 | frac4 | guard)) {
320 z->l.msw = ((rm == FSR_RM)? 0x80000000 : 0);
321 z->l.frac2 = z->l.frac3 = z->l.frac4 = 0;
322 return;
323 }
324
325 /* postnormalize */
326 if (lx < 0x10000) {
327 /* if cancellation occurred or the exponent is 1, */
328 /* the result is exact */
329 if (lx < 0x8000 || ex == 1) {
330 while ((lx | (frac2 & 0xfffe0000)) == 0 && ex > 32) {
331 lx = frac2;
332 frac2 = frac3;
333 frac3 = frac4;
334 frac4 = ((guard)? 0x80000000 : 0);
335 guard = 0;
336 ex -= 32;
337 }
338 while (lx < 0x10000 && ex > 1) {
339 lx = (lx << 1) | (frac2 >> 31);
340 frac2 = (frac2 << 1) | (frac3 >> 31);
341 frac3 = (frac3 << 1) | (frac4 >> 31);
342 frac4 <<= 1;
343 if (guard) {
344 frac4 |= 1;
345 guard = 0;
346 }
347 ex--;
348 }
349 if (lx >= 0x10000)
350 z->l.msw |= (ex << 16);
351 z->l.msw |= (lx & 0xffff);
352 z->l.frac2 = frac2;
353 z->l.frac3 = frac3;
354 z->l.frac4 = frac4;
355
356 /* if the result is tiny and underflow trapping is */
357 /* enabled, simulate underflow */
358 if (lx < 0x10000 && (*fsr & FSR_UFM))
359 *fsr |= FSR_UFC;
360 return;
361 }
362
363 /* otherwise we only borrowed one place */
364 lx = (lx << 1) | (frac2 >> 31);
365 frac2 = (frac2 << 1) | (frac3 >> 31);
366 frac3 = (frac3 << 1) | (frac4 >> 31);
367 frac4 <<= 1;
368 if (guard)
369 frac4 |= 1;
370 ex--;
371 } else {
372 sticky |= round;
373 round = guard;
374 }
375
376 /* fudge directed rounding modes as though the result were positive */
377 if (z->l.msw)
378 rm ^= (rm >> 1);
379
380 /* see if we need to round */
381 if (round | sticky) {
382 *fsr |= FSR_NXC;
383
384 /* round up if necessary */
385 if (rm == FSR_RP || (rm == FSR_RN && round &&
386 (sticky || (frac4 & 1)))) {
387 if (++frac4 == 0)
388 if (++frac3 == 0)
389 if (++frac2 == 0)
390 if (++lx >= 0x20000) {
391 lx >>= 1;
392 ex++;
393 }
394 }
395 }
396
397 /* store the result */
398 z->l.msw |= (ex << 16) | (lx & 0xffff);
399 z->l.frac2 = frac2;
400 z->l.frac3 = frac3;
401 z->l.frac4 = frac4;
402 }
403