![[BACK]](/icons/back.gif){kind=icon}
Return to dfrem.c CVS log ![[TXT]](/icons/text.gif){kind=icon}
![[DIR]](/icons/dir.gif){kind=icon}
Up to [local] / sys / arch / hppa / spmath|
Return to dfrem.c CVS log | Up to [local] / sys / arch / hppa / spmath |
1.1 nbrk 1: /* $OpenBSD: dfrem.c,v 1.5 2002/05/07 22:19:30 mickey Exp $ */
2: /*
3: (c) Copyright 1986 HEWLETT-PACKARD COMPANY
4: To anyone who acknowledges that this file is provided "AS IS"
5: without any express or implied warranty:
6: permission to use, copy, modify, and distribute this file
7: for any purpose is hereby granted without fee, provided that
8: the above copyright notice and this notice appears in all
9: copies, and that the name of Hewlett-Packard Company not be
10: used in advertising or publicity pertaining to distribution
11: of the software without specific, written prior permission.
12: Hewlett-Packard Company makes no representations about the
13: suitability of this software for any purpose.
14: */
15: /* @(#)dfrem.c: Revision: 1.7.88.1 Date: 93/12/07 15:05:43 */
16:
17: #include "float.h"
18: #include "dbl_float.h"
19:
20: /*
21: * Double Precision Floating-point Remainder
22: */
23: int
24: dbl_frem(srcptr1,srcptr2,dstptr,status)
25: dbl_floating_point *srcptr1, *srcptr2, *dstptr;
26: unsigned int *status;
27: {
28: register unsigned int opnd1p1, opnd1p2, opnd2p1, opnd2p2;
29: register unsigned int resultp1, resultp2;
30: register int opnd1_exponent, opnd2_exponent, dest_exponent, stepcount;
31: register int roundup = FALSE;
32:
33: Dbl_copyfromptr(srcptr1,opnd1p1,opnd1p2);
34: Dbl_copyfromptr(srcptr2,opnd2p1,opnd2p2);
35: /*
36: * check first operand for NaN's or infinity
37: */
38: if ((opnd1_exponent = Dbl_exponent(opnd1p1)) == DBL_INFINITY_EXPONENT) {
39: if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) {
40: if (Dbl_isnotnan(opnd2p1,opnd2p2)) {
41: /* invalid since first operand is infinity */
42: if (Is_invalidtrap_enabled())
43: return(INVALIDEXCEPTION);
44: Set_invalidflag();
45: Dbl_makequietnan(resultp1,resultp2);
46: Dbl_copytoptr(resultp1,resultp2,dstptr);
47: return(NOEXCEPTION);
48: }
49: }
50: else {
51: /*
52: * is NaN; signaling or quiet?
53: */
54: if (Dbl_isone_signaling(opnd1p1)) {
55: /* trap if INVALIDTRAP enabled */
56: if (Is_invalidtrap_enabled())
57: return(INVALIDEXCEPTION);
58: /* make NaN quiet */
59: Set_invalidflag();
60: Dbl_set_quiet(opnd1p1);
61: }
62: /*
63: * is second operand a signaling NaN?
64: */
65: else if (Dbl_is_signalingnan(opnd2p1)) {
66: /* trap if INVALIDTRAP enabled */
67: if (Is_invalidtrap_enabled())
68: return(INVALIDEXCEPTION);
69: /* make NaN quiet */
70: Set_invalidflag();
71: Dbl_set_quiet(opnd2p1);
72: Dbl_copytoptr(opnd2p1,opnd2p2,dstptr);
73: return(NOEXCEPTION);
74: }
75: /*
76: * return quiet NaN
77: */
78: Dbl_copytoptr(opnd1p1,opnd1p2,dstptr);
79: return(NOEXCEPTION);
80: }
81: }
82: /*
83: * check second operand for NaN's or infinity
84: */
85: if ((opnd2_exponent = Dbl_exponent(opnd2p1)) == DBL_INFINITY_EXPONENT) {
86: if (Dbl_iszero_mantissa(opnd2p1,opnd2p2)) {
87: /*
88: * return first operand
89: */
90: Dbl_copytoptr(opnd1p1,opnd1p2,dstptr);
91: return(NOEXCEPTION);
92: }
93: /*
94: * is NaN; signaling or quiet?
95: */
96: if (Dbl_isone_signaling(opnd2p1)) {
97: /* trap if INVALIDTRAP enabled */
98: if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
99: /* make NaN quiet */
100: Set_invalidflag();
101: Dbl_set_quiet(opnd2p1);
102: }
103: /*
104: * return quiet NaN
105: */
106: Dbl_copytoptr(opnd2p1,opnd2p2,dstptr);
107: return(NOEXCEPTION);
108: }
109: /*
110: * check second operand for zero
111: */
112: if (Dbl_iszero_exponentmantissa(opnd2p1,opnd2p2)) {
113: /* invalid since second operand is zero */
114: if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
115: Set_invalidflag();
116: Dbl_makequietnan(resultp1,resultp2);
117: Dbl_copytoptr(resultp1,resultp2,dstptr);
118: return(NOEXCEPTION);
119: }
120:
121: /*
122: * get sign of result
123: */
124: resultp1 = opnd1p1;
125:
126: /*
127: * check for denormalized operands
128: */
129: if (opnd1_exponent == 0) {
130: /* check for zero */
131: if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) {
132: Dbl_copytoptr(opnd1p1,opnd1p2,dstptr);
133: return(NOEXCEPTION);
134: }
135: /* normalize, then continue */
136: opnd1_exponent = 1;
137: Dbl_normalize(opnd1p1,opnd1p2,opnd1_exponent);
138: }
139: else {
140: Dbl_clear_signexponent_set_hidden(opnd1p1);
141: }
142: if (opnd2_exponent == 0) {
143: /* normalize, then continue */
144: opnd2_exponent = 1;
145: Dbl_normalize(opnd2p1,opnd2p2,opnd2_exponent);
146: }
147: else {
148: Dbl_clear_signexponent_set_hidden(opnd2p1);
149: }
150:
151: /* find result exponent and divide step loop count */
152: dest_exponent = opnd2_exponent - 1;
153: stepcount = opnd1_exponent - opnd2_exponent;
154:
155: /*
156: * check for opnd1/opnd2 < 1
157: */
158: if (stepcount < 0) {
159: /*
160: * check for opnd1/opnd2 > 1/2
161: *
162: * In this case n will round to 1, so
163: * r = opnd1 - opnd2
164: */
165: if (stepcount == -1 &&
166: Dbl_isgreaterthan(opnd1p1,opnd1p2,opnd2p1,opnd2p2)) {
167: /* set sign */
168: Dbl_allp1(resultp1) = ~Dbl_allp1(resultp1);
169: /* align opnd2 with opnd1 */
170: Dbl_leftshiftby1(opnd2p1,opnd2p2);
171: Dbl_subtract(opnd2p1,opnd2p2,opnd1p1,opnd1p2,
172: opnd2p1,opnd2p2);
173: /* now normalize */
174: while (Dbl_iszero_hidden(opnd2p1)) {
175: Dbl_leftshiftby1(opnd2p1,opnd2p2);
176: dest_exponent--;
177: }
178: Dbl_set_exponentmantissa(resultp1,resultp2,opnd2p1,opnd2p2);
179: goto testforunderflow;
180: }
181: /*
182: * opnd1/opnd2 <= 1/2
183: *
184: * In this case n will round to zero, so
185: * r = opnd1
186: */
187: Dbl_set_exponentmantissa(resultp1,resultp2,opnd1p1,opnd1p2);
188: dest_exponent = opnd1_exponent;
189: goto testforunderflow;
190: }
191:
192: /*
193: * Generate result
194: *
195: * Do iterative subtract until remainder is less than operand 2.
196: */
197: while (stepcount-- > 0 && (Dbl_allp1(opnd1p1) || Dbl_allp2(opnd1p2))) {
198: if (Dbl_isnotlessthan(opnd1p1,opnd1p2,opnd2p1,opnd2p2)) {
199: Dbl_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2,opnd1p1,opnd1p2);
200: }
201: Dbl_leftshiftby1(opnd1p1,opnd1p2);
202: }
203: /*
204: * Do last subtract, then determine which way to round if remainder
205: * is exactly 1/2 of opnd2
206: */
207: if (Dbl_isnotlessthan(opnd1p1,opnd1p2,opnd2p1,opnd2p2)) {
208: Dbl_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2,opnd1p1,opnd1p2);
209: roundup = TRUE;
210: }
211: if (stepcount > 0 || Dbl_iszero(opnd1p1,opnd1p2)) {
212: /* division is exact, remainder is zero */
213: Dbl_setzero_exponentmantissa(resultp1,resultp2);
214: Dbl_copytoptr(resultp1,resultp2,dstptr);
215: return(NOEXCEPTION);
216: }
217:
218: /*
219: * Check for cases where opnd1/opnd2 < n
220: *
221: * In this case the result's sign will be opposite that of
222: * opnd1. The mantissa also needs some correction.
223: */
224: Dbl_leftshiftby1(opnd1p1,opnd1p2);
225: if (Dbl_isgreaterthan(opnd1p1,opnd1p2,opnd2p1,opnd2p2)) {
226: Dbl_invert_sign(resultp1);
227: Dbl_leftshiftby1(opnd2p1,opnd2p2);
228: Dbl_subtract(opnd2p1,opnd2p2,opnd1p1,opnd1p2,opnd1p1,opnd1p2);
229: }
230: /* check for remainder being exactly 1/2 of opnd2 */
231: else if (Dbl_isequal(opnd1p1,opnd1p2,opnd2p1,opnd2p2) && roundup) {
232: Dbl_invert_sign(resultp1);
233: }
234:
235: /* normalize result's mantissa */
236: while (Dbl_iszero_hidden(opnd1p1)) {
237: dest_exponent--;
238: Dbl_leftshiftby1(opnd1p1,opnd1p2);
239: }
240: Dbl_set_exponentmantissa(resultp1,resultp2,opnd1p1,opnd1p2);
241:
242: /*
243: * Test for underflow
244: */
245: testforunderflow:
246: if (dest_exponent <= 0) {
247: /* trap if UNDERFLOWTRAP enabled */
248: if (Is_underflowtrap_enabled()) {
249: /*
250: * Adjust bias of result
251: */
252: Dbl_setwrapped_exponent(resultp1,dest_exponent,unfl);
253: /* frem is always exact */
254: Dbl_copytoptr(resultp1,resultp2,dstptr);
255: return(UNDERFLOWEXCEPTION);
256: }
257: /*
258: * denormalize result or set to signed zero
259: */
260: if (dest_exponent >= (1 - DBL_P)) {
261: Dbl_rightshift_exponentmantissa(resultp1,resultp2,
262: 1-dest_exponent);
263: }
264: else {
265: Dbl_setzero_exponentmantissa(resultp1,resultp2);
266: }
267: }
268: else Dbl_set_exponent(resultp1,dest_exponent);
269: Dbl_copytoptr(resultp1,resultp2,dstptr);
270: return(NOEXCEPTION);
271: }