![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to fpu_div.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [Isaki's NoNo m68k/m88k emulator] / nono / fpe|
Return to fpu_div.c CVS log | Up to [Isaki's NoNo m68k/m88k emulator] / nono / fpe |
1.1 root 1: /* $NetBSD: fpu_div.c,v 1.10 2014/01/01 05:23:40 isaki Exp $ */
2:
3: /*
4: * Copyright (c) 1992, 1993
5: * The Regents of the University of California. All rights reserved.
6: *
7: * This software was developed by the Computer Systems Engineering group
8: * at Lawrence Berkeley Laboratory under DARPA contract BG 91-66 and
9: * contributed to Berkeley.
10: *
11: * All advertising materials mentioning features or use of this software
12: * must display the following acknowledgement:
13: * This product includes software developed by the University of
14: * California, Lawrence Berkeley Laboratory.
15: *
16: * Redistribution and use in source and binary forms, with or without
17: * modification, are permitted provided that the following conditions
18: * are met:
19: * 1. Redistributions of source code must retain the above copyright
20: * notice, this list of conditions and the following disclaimer.
21: * 2. Redistributions in binary form must reproduce the above copyright
22: * notice, this list of conditions and the following disclaimer in the
23: * documentation and/or other materials provided with the distribution.
24: * 3. Neither the name of the University nor the names of its contributors
25: * may be used to endorse or promote products derived from this software
26: * without specific prior written permission.
27: *
28: * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
29: * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
30: * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
31: * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
32: * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
33: * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
34: * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
35: * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
36: * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
37: * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
38: * SUCH DAMAGE.
39: *
40: * @(#)fpu_div.c 8.1 (Berkeley) 6/11/93
41: */
42:
43: /*
44: * Perform an FPU divide (return x / y).
45: */
46:
47: #include "fpu_arith.h"
48: #include "fpu_emulate.h"
49:
50: /*
51: * Division of normal numbers is done as follows:
52: *
53: * x and y are floating point numbers, i.e., in the form 1.bbbb * 2^e.
54: * If X and Y are the mantissas (1.bbbb's), the quotient is then:
55: *
56: * q = (X / Y) * 2^((x exponent) - (y exponent))
57: *
58: * Since X and Y are both in [1.0,2.0), the quotient's mantissa (X / Y)
59: * will be in [0.5,2.0). Moreover, it will be less than 1.0 if and only
60: * if X < Y. In that case, it will have to be shifted left one bit to
61: * become a normal number, and the exponent decremented. Thus, the
62: * desired exponent is:
63: *
64: * left_shift = x->fp_mant < y->fp_mant;
65: * result_exp = x->fp_exp - y->fp_exp - left_shift;
66: *
67: * The quotient mantissa X/Y can then be computed one bit at a time
68: * using the following algorithm:
69: *
70: * Q = 0; -- Initial quotient.
71: * R = X; -- Initial remainder,
72: * if (left_shift) -- but fixed up in advance.
73: * R *= 2;
74: * for (bit = FP_NMANT; --bit >= 0; R *= 2) {
75: * if (R >= Y) {
76: * Q |= 1 << bit;
77: * R -= Y;
78: * }
79: * }
80: *
81: * The subtraction R -= Y always removes the uppermost bit from R (and
82: * can sometimes remove additional lower-order 1 bits); this proof is
83: * left to the reader.
84: *
85: * This loop correctly calculates the guard and round bits since they are
86: * included in the expanded internal representation. The sticky bit
87: * is to be set if and only if any other bits beyond guard and round
88: * would be set. From the above it is obvious that this is true if and
89: * only if the remainder R is nonzero when the loop terminates.
90: *
91: * Examining the loop above, we can see that the quotient Q is built
92: * one bit at a time ``from the top down''. This means that we can
93: * dispense with the multi-word arithmetic and just build it one word
94: * at a time, writing each result word when it is done.
95: *
96: * Furthermore, since X and Y are both in [1.0,2.0), we know that,
97: * initially, R >= Y. (Recall that, if X < Y, R is set to X * 2 and
98: * is therefore at in [2.0,4.0).) Thus Q is sure to have bit FP_NMANT-1
99: * set, and R can be set initially to either X - Y (when X >= Y) or
100: * 2X - Y (when X < Y). In addition, comparing R and Y is difficult,
101: * so we will simply calculate R - Y and see if that underflows.
102: * This leads to the following revised version of the algorithm:
103: *
104: * R = X;
105: * bit = FP_1;
106: * D = R - Y;
107: * if (D >= 0) {
108: * result_exp = x->fp_exp - y->fp_exp;
109: * R = D;
110: * q = bit;
111: * bit >>= 1;
112: * } else {
113: * result_exp = x->fp_exp - y->fp_exp - 1;
114: * q = 0;
115: * }
116: * R <<= 1;
117: * do {
118: * D = R - Y;
119: * if (D >= 0) {
120: * q |= bit;
121: * R = D;
122: * }
123: * R <<= 1;
124: * } while ((bit >>= 1) != 0);
125: * Q[0] = q;
126: * for (i = 1; i < 4; i++) {
127: * q = 0, bit = 1 << 31;
128: * do {
129: * D = R - Y;
130: * if (D >= 0) {
131: * q |= bit;
132: * R = D;
133: * }
134: * R <<= 1;
135: * } while ((bit >>= 1) != 0);
136: * Q[i] = q;
137: * }
138: *
139: * This can be refined just a bit further by moving the `R <<= 1'
140: * calculations to the front of the do-loops and eliding the first one.
141: * The process can be terminated immediately whenever R becomes 0, but
142: * this is relatively rare, and we do not bother.
143: */
144:
145: struct fpn *
146: fpu_div(struct fpemu *fe)
147: {
148: struct fpn *x = &fe->fe_f1, *y = &fe->fe_f2;
149: uint32_t q, bit;
150: uint32_t r0, r1, r2, d0, d1, d2, y0, y1, y2;
151: FPU_DECL_CARRY
152:
153: fe->fe_fpsr &= ~FPSR_EXCP; /* clear all exceptions */
154:
155: /*
156: * Since divide is not commutative, we cannot just use ORDER.
157: * Check either operand for NaN first; if there is at least one,
158: * order the signalling one (if only one) onto the right, then
159: * return it. Otherwise we have the following cases:
160: *
161: * Inf / Inf = NaN, plus NV exception
162: * Inf / num = Inf
163: * Inf / 0 = Inf
164: * 0 / Inf = 0
165: * 0 / num = 0
166: * 0 / 0 = NaN, plus NV exception
167: * num / Inf = 0
168: * num / num = num (do the divide)
169: * num / 0 = Inf, plus DZ exception
170: */
171: if (ISNAN(x) || ISNAN(y)) {
172: ORDER(x, y);
173: return (y);
174: }
175: if (ISINF(x) || ISZERO(x)) {
176: if (x->fp_class == y->fp_class)
177: return (fpu_newnan(fe));
178: /* all results at this point use XOR of operand signs */
179: x->fp_sign ^= y->fp_sign;
180: return (x);
181: }
182:
183: /* all results at this point use XOR of operand signs */
184: x->fp_sign ^= y->fp_sign;
185: if (ISINF(y)) {
186: x->fp_class = FPC_ZERO;
187: return (x);
188: }
189: if (ISZERO(y)) {
190: fe->fe_fpsr |= FPSR_DZ;
191: x->fp_class = FPC_INF;
192: return (x);
193: }
194:
195: /*
196: * Macros for the divide. See comments at top for algorithm.
197: * Note that we expand R, D, and Y here.
198: */
199:
200: #define SUBTRACT /* D = R - Y */ \
201: FPU_SUBS(d2, r2, y2); \
202: FPU_SUBCS(d1, r1, y1); FPU_SUBC(d0, r0, y0)
203:
204: #define NONNEGATIVE /* D >= 0 */ \
205: ((int)d0 >= 0)
206:
207: #ifdef FPU_SHL1_BY_ADD
208: #define SHL1 /* R <<= 1 */ \
209: FPU_ADDS(r2, r2, r2); \
210: FPU_ADDCS(r1, r1, r1); FPU_ADDC(r0, r0, r0)
211: #else
212: #define SHL1 \
213: r0 = (r0 << 1) | (r1 >> 31), r1 = (r1 << 1) | (r2 >> 31), \
214: r2 <<= 1
215: #endif
216:
217: #define LOOP /* do ... while (bit >>= 1) */ \
218: do { \
219: SHL1; \
220: SUBTRACT; \
221: if (NONNEGATIVE) { \
222: q |= bit; \
223: r0 = d0, r1 = d1, r2 = d2; \
224: } \
225: } while ((bit >>= 1) != 0)
226:
227: #define WORD(r, i) /* calculate r->fp_mant[i] */ \
228: q = 0; \
229: bit = 1 << 31; \
230: LOOP; \
231: (x)->fp_mant[i] = q
232:
233: /* Setup. Note that we put our result in x. */
234: r0 = x->fp_mant[0];
235: r1 = x->fp_mant[1];
236: r2 = x->fp_mant[2];
237: y0 = y->fp_mant[0];
238: y1 = y->fp_mant[1];
239: y2 = y->fp_mant[2];
240:
241: bit = FP_1;
242: SUBTRACT;
243: if (NONNEGATIVE) {
244: x->fp_exp -= y->fp_exp;
245: r0 = d0, r1 = d1, r2 = d2;
246: q = bit;
247: bit >>= 1;
248: } else {
249: x->fp_exp -= y->fp_exp + 1;
250: q = 0;
251: }
252: LOOP;
253: x->fp_mant[0] = q;
254: WORD(x, 1);
255: WORD(x, 2);
256: x->fp_sticky = r0 | r1 | r2;
257:
258: return (x);
259: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.