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