aboutsummaryrefslogtreecommitdiffstats
path: root/arch/parisc/math-emu/dfadd.c
diff options
context:
space:
mode:
Diffstat (limited to 'arch/parisc/math-emu/dfadd.c')
-rw-r--r--arch/parisc/math-emu/dfadd.c524
1 files changed, 524 insertions, 0 deletions
diff --git a/arch/parisc/math-emu/dfadd.c b/arch/parisc/math-emu/dfadd.c
new file mode 100644
index 000000000000..e147d7d3b0f4
--- /dev/null
+++ b/arch/parisc/math-emu/dfadd.c
@@ -0,0 +1,524 @@
1/*
2 * Linux/PA-RISC Project (http://www.parisc-linux.org/)
3 *
4 * Floating-point emulation code
5 * Copyright (C) 2001 Hewlett-Packard (Paul Bame) <bame@debian.org>
6 *
7 * This program is free software; you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation; either version 2, or (at your option)
10 * any later version.
11 *
12 * This program is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with this program; if not, write to the Free Software
19 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 */
21/*
22 * BEGIN_DESC
23 *
24 * File:
25 * @(#) pa/spmath/dfadd.c $Revision: 1.1 $
26 *
27 * Purpose:
28 * Double_add: add two double precision values.
29 *
30 * External Interfaces:
31 * dbl_fadd(leftptr, rightptr, dstptr, status)
32 *
33 * Internal Interfaces:
34 *
35 * Theory:
36 * <<please update with a overview of the operation of this file>>
37 *
38 * END_DESC
39*/
40
41
42#include "float.h"
43#include "dbl_float.h"
44
45/*
46 * Double_add: add two double precision values.
47 */
48dbl_fadd(
49 dbl_floating_point *leftptr,
50 dbl_floating_point *rightptr,
51 dbl_floating_point *dstptr,
52 unsigned int *status)
53{
54 register unsigned int signless_upper_left, signless_upper_right, save;
55 register unsigned int leftp1, leftp2, rightp1, rightp2, extent;
56 register unsigned int resultp1 = 0, resultp2 = 0;
57
58 register int result_exponent, right_exponent, diff_exponent;
59 register int sign_save, jumpsize;
60 register boolean inexact = FALSE;
61 register boolean underflowtrap;
62
63 /* Create local copies of the numbers */
64 Dbl_copyfromptr(leftptr,leftp1,leftp2);
65 Dbl_copyfromptr(rightptr,rightp1,rightp2);
66
67 /* A zero "save" helps discover equal operands (for later), *
68 * and is used in swapping operands (if needed). */
69 Dbl_xortointp1(leftp1,rightp1,/*to*/save);
70
71 /*
72 * check first operand for NaN's or infinity
73 */
74 if ((result_exponent = Dbl_exponent(leftp1)) == DBL_INFINITY_EXPONENT)
75 {
76 if (Dbl_iszero_mantissa(leftp1,leftp2))
77 {
78 if (Dbl_isnotnan(rightp1,rightp2))
79 {
80 if (Dbl_isinfinity(rightp1,rightp2) && save!=0)
81 {
82 /*
83 * invalid since operands are opposite signed infinity's
84 */
85 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
86 Set_invalidflag();
87 Dbl_makequietnan(resultp1,resultp2);
88 Dbl_copytoptr(resultp1,resultp2,dstptr);
89 return(NOEXCEPTION);
90 }
91 /*
92 * return infinity
93 */
94 Dbl_copytoptr(leftp1,leftp2,dstptr);
95 return(NOEXCEPTION);
96 }
97 }
98 else
99 {
100 /*
101 * is NaN; signaling or quiet?
102 */
103 if (Dbl_isone_signaling(leftp1))
104 {
105 /* trap if INVALIDTRAP enabled */
106 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
107 /* make NaN quiet */
108 Set_invalidflag();
109 Dbl_set_quiet(leftp1);
110 }
111 /*
112 * is second operand a signaling NaN?
113 */
114 else if (Dbl_is_signalingnan(rightp1))
115 {
116 /* trap if INVALIDTRAP enabled */
117 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
118 /* make NaN quiet */
119 Set_invalidflag();
120 Dbl_set_quiet(rightp1);
121 Dbl_copytoptr(rightp1,rightp2,dstptr);
122 return(NOEXCEPTION);
123 }
124 /*
125 * return quiet NaN
126 */
127 Dbl_copytoptr(leftp1,leftp2,dstptr);
128 return(NOEXCEPTION);
129 }
130 } /* End left NaN or Infinity processing */
131 /*
132 * check second operand for NaN's or infinity
133 */
134 if (Dbl_isinfinity_exponent(rightp1))
135 {
136 if (Dbl_iszero_mantissa(rightp1,rightp2))
137 {
138 /* return infinity */
139 Dbl_copytoptr(rightp1,rightp2,dstptr);
140 return(NOEXCEPTION);
141 }
142 /*
143 * is NaN; signaling or quiet?
144 */
145 if (Dbl_isone_signaling(rightp1))
146 {
147 /* trap if INVALIDTRAP enabled */
148 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
149 /* make NaN quiet */
150 Set_invalidflag();
151 Dbl_set_quiet(rightp1);
152 }
153 /*
154 * return quiet NaN
155 */
156 Dbl_copytoptr(rightp1,rightp2,dstptr);
157 return(NOEXCEPTION);
158 } /* End right NaN or Infinity processing */
159
160 /* Invariant: Must be dealing with finite numbers */
161
162 /* Compare operands by removing the sign */
163 Dbl_copytoint_exponentmantissap1(leftp1,signless_upper_left);
164 Dbl_copytoint_exponentmantissap1(rightp1,signless_upper_right);
165
166 /* sign difference selects add or sub operation. */
167 if(Dbl_ismagnitudeless(leftp2,rightp2,signless_upper_left,signless_upper_right))
168 {
169 /* Set the left operand to the larger one by XOR swap *
170 * First finish the first word using "save" */
171 Dbl_xorfromintp1(save,rightp1,/*to*/rightp1);
172 Dbl_xorfromintp1(save,leftp1,/*to*/leftp1);
173 Dbl_swap_lower(leftp2,rightp2);
174 result_exponent = Dbl_exponent(leftp1);
175 }
176 /* Invariant: left is not smaller than right. */
177
178 if((right_exponent = Dbl_exponent(rightp1)) == 0)
179 {
180 /* Denormalized operands. First look for zeroes */
181 if(Dbl_iszero_mantissa(rightp1,rightp2))
182 {
183 /* right is zero */
184 if(Dbl_iszero_exponentmantissa(leftp1,leftp2))
185 {
186 /* Both operands are zeros */
187 if(Is_rounding_mode(ROUNDMINUS))
188 {
189 Dbl_or_signs(leftp1,/*with*/rightp1);
190 }
191 else
192 {
193 Dbl_and_signs(leftp1,/*with*/rightp1);
194 }
195 }
196 else
197 {
198 /* Left is not a zero and must be the result. Trapped
199 * underflows are signaled if left is denormalized. Result
200 * is always exact. */
201 if( (result_exponent == 0) && Is_underflowtrap_enabled() )
202 {
203 /* need to normalize results mantissa */
204 sign_save = Dbl_signextendedsign(leftp1);
205 Dbl_leftshiftby1(leftp1,leftp2);
206 Dbl_normalize(leftp1,leftp2,result_exponent);
207 Dbl_set_sign(leftp1,/*using*/sign_save);
208 Dbl_setwrapped_exponent(leftp1,result_exponent,unfl);
209 Dbl_copytoptr(leftp1,leftp2,dstptr);
210 /* inexact = FALSE */
211 return(UNDERFLOWEXCEPTION);
212 }
213 }
214 Dbl_copytoptr(leftp1,leftp2,dstptr);
215 return(NOEXCEPTION);
216 }
217
218 /* Neither are zeroes */
219 Dbl_clear_sign(rightp1); /* Exponent is already cleared */
220 if(result_exponent == 0 )
221 {
222 /* Both operands are denormalized. The result must be exact
223 * and is simply calculated. A sum could become normalized and a
224 * difference could cancel to a true zero. */
225 if( (/*signed*/int) save < 0 )
226 {
227 Dbl_subtract(leftp1,leftp2,/*minus*/rightp1,rightp2,
228 /*into*/resultp1,resultp2);
229 if(Dbl_iszero_mantissa(resultp1,resultp2))
230 {
231 if(Is_rounding_mode(ROUNDMINUS))
232 {
233 Dbl_setone_sign(resultp1);
234 }
235 else
236 {
237 Dbl_setzero_sign(resultp1);
238 }
239 Dbl_copytoptr(resultp1,resultp2,dstptr);
240 return(NOEXCEPTION);
241 }
242 }
243 else
244 {
245 Dbl_addition(leftp1,leftp2,rightp1,rightp2,
246 /*into*/resultp1,resultp2);
247 if(Dbl_isone_hidden(resultp1))
248 {
249 Dbl_copytoptr(resultp1,resultp2,dstptr);
250 return(NOEXCEPTION);
251 }
252 }
253 if(Is_underflowtrap_enabled())
254 {
255 /* need to normalize result */
256 sign_save = Dbl_signextendedsign(resultp1);
257 Dbl_leftshiftby1(resultp1,resultp2);
258 Dbl_normalize(resultp1,resultp2,result_exponent);
259 Dbl_set_sign(resultp1,/*using*/sign_save);
260 Dbl_setwrapped_exponent(resultp1,result_exponent,unfl);
261 Dbl_copytoptr(resultp1,resultp2,dstptr);
262 /* inexact = FALSE */
263 return(UNDERFLOWEXCEPTION);
264 }
265 Dbl_copytoptr(resultp1,resultp2,dstptr);
266 return(NOEXCEPTION);
267 }
268 right_exponent = 1; /* Set exponent to reflect different bias
269 * with denomalized numbers. */
270 }
271 else
272 {
273 Dbl_clear_signexponent_set_hidden(rightp1);
274 }
275 Dbl_clear_exponent_set_hidden(leftp1);
276 diff_exponent = result_exponent - right_exponent;
277
278 /*
279 * Special case alignment of operands that would force alignment
280 * beyond the extent of the extension. A further optimization
281 * could special case this but only reduces the path length for this
282 * infrequent case.
283 */
284 if(diff_exponent > DBL_THRESHOLD)
285 {
286 diff_exponent = DBL_THRESHOLD;
287 }
288
289 /* Align right operand by shifting to right */
290 Dbl_right_align(/*operand*/rightp1,rightp2,/*shifted by*/diff_exponent,
291 /*and lower to*/extent);
292
293 /* Treat sum and difference of the operands separately. */
294 if( (/*signed*/int) save < 0 )
295 {
296 /*
297 * Difference of the two operands. Their can be no overflow. A
298 * borrow can occur out of the hidden bit and force a post
299 * normalization phase.
300 */
301 Dbl_subtract_withextension(leftp1,leftp2,/*minus*/rightp1,rightp2,
302 /*with*/extent,/*into*/resultp1,resultp2);
303 if(Dbl_iszero_hidden(resultp1))
304 {
305 /* Handle normalization */
306 /* A straight foward algorithm would now shift the result
307 * and extension left until the hidden bit becomes one. Not
308 * all of the extension bits need participate in the shift.
309 * Only the two most significant bits (round and guard) are
310 * needed. If only a single shift is needed then the guard
311 * bit becomes a significant low order bit and the extension
312 * must participate in the rounding. If more than a single
313 * shift is needed, then all bits to the right of the guard
314 * bit are zeros, and the guard bit may or may not be zero. */
315 sign_save = Dbl_signextendedsign(resultp1);
316 Dbl_leftshiftby1_withextent(resultp1,resultp2,extent,resultp1,resultp2);
317
318 /* Need to check for a zero result. The sign and exponent
319 * fields have already been zeroed. The more efficient test
320 * of the full object can be used.
321 */
322 if(Dbl_iszero(resultp1,resultp2))
323 /* Must have been "x-x" or "x+(-x)". */
324 {
325 if(Is_rounding_mode(ROUNDMINUS)) Dbl_setone_sign(resultp1);
326 Dbl_copytoptr(resultp1,resultp2,dstptr);
327 return(NOEXCEPTION);
328 }
329 result_exponent--;
330 /* Look to see if normalization is finished. */
331 if(Dbl_isone_hidden(resultp1))
332 {
333 if(result_exponent==0)
334 {
335 /* Denormalized, exponent should be zero. Left operand *
336 * was normalized, so extent (guard, round) was zero */
337 goto underflow;
338 }
339 else
340 {
341 /* No further normalization is needed. */
342 Dbl_set_sign(resultp1,/*using*/sign_save);
343 Ext_leftshiftby1(extent);
344 goto round;
345 }
346 }
347
348 /* Check for denormalized, exponent should be zero. Left *
349 * operand was normalized, so extent (guard, round) was zero */
350 if(!(underflowtrap = Is_underflowtrap_enabled()) &&
351 result_exponent==0) goto underflow;
352
353 /* Shift extension to complete one bit of normalization and
354 * update exponent. */
355 Ext_leftshiftby1(extent);
356
357 /* Discover first one bit to determine shift amount. Use a
358 * modified binary search. We have already shifted the result
359 * one position right and still not found a one so the remainder
360 * of the extension must be zero and simplifies rounding. */
361 /* Scan bytes */
362 while(Dbl_iszero_hiddenhigh7mantissa(resultp1))
363 {
364 Dbl_leftshiftby8(resultp1,resultp2);
365 if((result_exponent -= 8) <= 0 && !underflowtrap)
366 goto underflow;
367 }
368 /* Now narrow it down to the nibble */
369 if(Dbl_iszero_hiddenhigh3mantissa(resultp1))
370 {
371 /* The lower nibble contains the normalizing one */
372 Dbl_leftshiftby4(resultp1,resultp2);
373 if((result_exponent -= 4) <= 0 && !underflowtrap)
374 goto underflow;
375 }
376 /* Select case were first bit is set (already normalized)
377 * otherwise select the proper shift. */
378 if((jumpsize = Dbl_hiddenhigh3mantissa(resultp1)) > 7)
379 {
380 /* Already normalized */
381 if(result_exponent <= 0) goto underflow;
382 Dbl_set_sign(resultp1,/*using*/sign_save);
383 Dbl_set_exponent(resultp1,/*using*/result_exponent);
384 Dbl_copytoptr(resultp1,resultp2,dstptr);
385 return(NOEXCEPTION);
386 }
387 Dbl_sethigh4bits(resultp1,/*using*/sign_save);
388 switch(jumpsize)
389 {
390 case 1:
391 {
392 Dbl_leftshiftby3(resultp1,resultp2);
393 result_exponent -= 3;
394 break;
395 }
396 case 2:
397 case 3:
398 {
399 Dbl_leftshiftby2(resultp1,resultp2);
400 result_exponent -= 2;
401 break;
402 }
403 case 4:
404 case 5:
405 case 6:
406 case 7:
407 {
408 Dbl_leftshiftby1(resultp1,resultp2);
409 result_exponent -= 1;
410 break;
411 }
412 }
413 if(result_exponent > 0)
414 {
415 Dbl_set_exponent(resultp1,/*using*/result_exponent);
416 Dbl_copytoptr(resultp1,resultp2,dstptr);
417 return(NOEXCEPTION); /* Sign bit is already set */
418 }
419 /* Fixup potential underflows */
420 underflow:
421 if(Is_underflowtrap_enabled())
422 {
423 Dbl_set_sign(resultp1,sign_save);
424 Dbl_setwrapped_exponent(resultp1,result_exponent,unfl);
425 Dbl_copytoptr(resultp1,resultp2,dstptr);
426 /* inexact = FALSE */
427 return(UNDERFLOWEXCEPTION);
428 }
429 /*
430 * Since we cannot get an inexact denormalized result,
431 * we can now return.
432 */
433 Dbl_fix_overshift(resultp1,resultp2,(1-result_exponent),extent);
434 Dbl_clear_signexponent(resultp1);
435 Dbl_set_sign(resultp1,sign_save);
436 Dbl_copytoptr(resultp1,resultp2,dstptr);
437 return(NOEXCEPTION);
438 } /* end if(hidden...)... */
439 /* Fall through and round */
440 } /* end if(save < 0)... */
441 else
442 {
443 /* Add magnitudes */
444 Dbl_addition(leftp1,leftp2,rightp1,rightp2,/*to*/resultp1,resultp2);
445 if(Dbl_isone_hiddenoverflow(resultp1))
446 {
447 /* Prenormalization required. */
448 Dbl_rightshiftby1_withextent(resultp2,extent,extent);
449 Dbl_arithrightshiftby1(resultp1,resultp2);
450 result_exponent++;
451 } /* end if hiddenoverflow... */
452 } /* end else ...add magnitudes... */
453
454 /* Round the result. If the extension is all zeros,then the result is
455 * exact. Otherwise round in the correct direction. No underflow is
456 * possible. If a postnormalization is necessary, then the mantissa is
457 * all zeros so no shift is needed. */
458 round:
459 if(Ext_isnotzero(extent))
460 {
461 inexact = TRUE;
462 switch(Rounding_mode())
463 {
464 case ROUNDNEAREST: /* The default. */
465 if(Ext_isone_sign(extent))
466 {
467 /* at least 1/2 ulp */
468 if(Ext_isnotzero_lower(extent) ||
469 Dbl_isone_lowmantissap2(resultp2))
470 {
471 /* either exactly half way and odd or more than 1/2ulp */
472 Dbl_increment(resultp1,resultp2);
473 }
474 }
475 break;
476
477 case ROUNDPLUS:
478 if(Dbl_iszero_sign(resultp1))
479 {
480 /* Round up positive results */
481 Dbl_increment(resultp1,resultp2);
482 }
483 break;
484
485 case ROUNDMINUS:
486 if(Dbl_isone_sign(resultp1))
487 {
488 /* Round down negative results */
489 Dbl_increment(resultp1,resultp2);
490 }
491
492 case ROUNDZERO:;
493 /* truncate is simple */
494 } /* end switch... */
495 if(Dbl_isone_hiddenoverflow(resultp1)) result_exponent++;
496 }
497 if(result_exponent == DBL_INFINITY_EXPONENT)
498 {
499 /* Overflow */
500 if(Is_overflowtrap_enabled())
501 {
502 Dbl_setwrapped_exponent(resultp1,result_exponent,ovfl);
503 Dbl_copytoptr(resultp1,resultp2,dstptr);
504 if (inexact)
505 if (Is_inexacttrap_enabled())
506 return(OVERFLOWEXCEPTION | INEXACTEXCEPTION);
507 else Set_inexactflag();
508 return(OVERFLOWEXCEPTION);
509 }
510 else
511 {
512 inexact = TRUE;
513 Set_overflowflag();
514 Dbl_setoverflow(resultp1,resultp2);
515 }
516 }
517 else Dbl_set_exponent(resultp1,result_exponent);
518 Dbl_copytoptr(resultp1,resultp2,dstptr);
519 if(inexact)
520 if(Is_inexacttrap_enabled())
521 return(INEXACTEXCEPTION);
522 else Set_inexactflag();
523 return(NOEXCEPTION);
524}
generated by cgit v1.2.2 (git 2.25.0) at 2024-12-25 12:38:02 -0500