1 /* $NetBSD: fpu_rem.c,v 1.17 2015/02/05 12:22:06 isaki Exp $ */ 2 3 /* 4 * Copyright (c) 1995 Ken Nakata 5 * All rights reserved. 6 * 7 * Redistribution and use in source and binary forms, with or without 8 * modification, are permitted provided that the following conditions 9 * are met: 10 * 1. Redistributions of source code must retain the above copyright 11 * notice, this list of conditions and the following disclaimer. 12 * 2. Redistributions in binary form must reproduce the above copyright 13 * notice, this list of conditions and the following disclaimer in the 14 * documentation and/or other materials provided with the distribution. 15 * 3. Neither the name of the author nor the names of its contributors 16 * may be used to endorse or promote products derived from this software 17 * without specific prior written permission. 18 * 19 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND 20 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 21 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 22 * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE 23 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 24 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 25 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 26 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 27 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 28 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 29 * SUCH DAMAGE. 30 * 31 * @(#)fpu_rem.c 10/24/95 32 */ 33 34 #include <sys/cdefs.h> 35 __KERNEL_RCSID(0, "$NetBSD: fpu_rem.c,v 1.17 2015/02/05 12:22:06 isaki Exp $"); 36 37 #include <sys/types.h> 38 #include <sys/signal.h> 39 #include <machine/frame.h> 40 41 #include "fpu_emulate.h" 42 43 /* 44 * ALGORITHM 45 * 46 * Step 1. Save and strip signs of X and Y: signX := sign(X), 47 * signY := sign(Y), X := *X*, Y := *Y*, 48 * signQ := signX EOR signY. Record whether MOD or REM 49 * is requested. 50 * 51 * Step 2. Set L := expo(X)-expo(Y), k := 0, Q := 0. 52 * If (L < 0) then 53 * R := X, go to Step 4. 54 * else 55 * R := 2^(-L)X, j := L. 56 * endif 57 * 58 * Step 3. Perform MOD(X,Y) 59 * 3.1 If R = Y, then { Q := Q + 1, R := 0, go to Step 7. } 60 * 3.2 If R > Y, then { R := R - Y, Q := Q + 1} 61 * 3.3 If j = 0, go to Step 4. 62 * 3.4 k := k + 1, j := j - 1, Q := 2Q, R := 2R. Go to 63 * Step 3.1. 64 * 65 * Step 4. R := signX*R. 66 * 67 * Step 5. If MOD is requested, go to Step 7. 68 * 69 * Step 6. Now, R = MOD(X,Y), convert to REM(X,Y) is requested. 70 * Do banker's rounding. 71 * If abs(R) > Y/2 72 * || (abs(R) == Y/2 && Q % 2 == 1) then 73 * { Q := Q + 1, R := R - signX * Y }. 74 * 75 * Step 7. Return signQ, last 7 bits of Q, and R as required. 76 */ 77 78 static struct fpn * __fpu_modrem(struct fpemu *fe, int is_mod); 79 static int abscmp3(const struct fpn *a, const struct fpn *b); 80 81 /* Absolute FORTRAN Compare */ 82 static int 83 abscmp3(const struct fpn *a, const struct fpn *b) 84 { 85 int i; 86 87 if (a->fp_exp < b->fp_exp) { 88 return -1; 89 } else if (a->fp_exp > b->fp_exp) { 90 return 1; 91 } else { 92 for (i = 0; i < 3; i++) { 93 if (a->fp_mant[i] < b->fp_mant[i]) 94 return -1; 95 else if (a->fp_mant[i] > b->fp_mant[i]) 96 return 1; 97 } 98 } 99 return 0; 100 } 101 102 static struct fpn * 103 __fpu_modrem(struct fpemu *fe, int is_mod) 104 { 105 static struct fpn X, Y; 106 struct fpn *x, *y, *r; 107 uint32_t signX, signY, signQ; 108 int j, k, l, q; 109 int cmp; 110 111 if (ISNAN(&fe->fe_f1) || ISNAN(&fe->fe_f2)) 112 return fpu_newnan(fe); 113 if (ISINF(&fe->fe_f1) || ISZERO(&fe->fe_f2)) 114 return fpu_newnan(fe); 115 116 CPYFPN(&X, &fe->fe_f1); 117 CPYFPN(&Y, &fe->fe_f2); 118 x = &X; 119 y = &Y; 120 q = 0; 121 r = &fe->fe_f2; 122 123 /* 124 * Step 1 125 */ 126 signX = x->fp_sign; 127 signY = y->fp_sign; 128 signQ = (signX ^ signY); 129 x->fp_sign = y->fp_sign = 0; 130 131 /* Special treatment that just return input value but Q is necessary */ 132 if (ISZERO(x) || ISINF(y)) { 133 r = &fe->fe_f1; 134 goto Step7; 135 } 136 137 /* 138 * Step 2 139 */ 140 l = x->fp_exp - y->fp_exp; 141 k = 0; 142 CPYFPN(r, x); 143 if (l >= 0) { 144 r->fp_exp -= l; 145 j = l; 146 147 /* 148 * Step 3 149 */ 150 for (;;) { 151 cmp = abscmp3(r, y); 152 153 /* Step 3.1 */ 154 if (cmp == 0) 155 break; 156 157 /* Step 3.2 */ 158 if (cmp > 0) { 159 CPYFPN(&fe->fe_f1, r); 160 CPYFPN(&fe->fe_f2, y); 161 fe->fe_f2.fp_sign = 1; 162 r = fpu_add(fe); 163 q++; 164 } 165 166 /* Step 3.3 */ 167 if (j == 0) 168 goto Step4; 169 170 /* Step 3.4 */ 171 k++; 172 j--; 173 q += q; 174 r->fp_exp++; 175 } 176 /* R == Y */ 177 q++; 178 r->fp_class = FPC_ZERO; 179 goto Step7; 180 } 181 Step4: 182 r->fp_sign = signX; 183 184 /* 185 * Step 5 186 */ 187 if (is_mod) 188 goto Step7; 189 190 /* 191 * Step 6 192 */ 193 /* y = y / 2 */ 194 y->fp_exp--; 195 /* abscmp3 ignore sign */ 196 cmp = abscmp3(r, y); 197 /* revert y */ 198 y->fp_exp++; 199 200 if (cmp > 0 || (cmp == 0 && q % 2)) { 201 q++; 202 CPYFPN(&fe->fe_f1, r); 203 CPYFPN(&fe->fe_f2, y); 204 fe->fe_f2.fp_sign = !signX; 205 r = fpu_add(fe); 206 } 207 208 /* 209 * Step 7 210 */ 211 Step7: 212 q &= 0x7f; 213 q |= (signQ << 7); 214 fe->fe_fpframe->fpf_fpsr = 215 fe->fe_fpsr = 216 (fe->fe_fpsr & ~FPSR_QTT) | (q << 16); 217 return r; 218 } 219 220 struct fpn * 221 fpu_rem(struct fpemu *fe) 222 { 223 return __fpu_modrem(fe, 0); 224 } 225 226 struct fpn * 227 fpu_mod(struct fpemu *fe) 228 { 229 return __fpu_modrem(fe, 1); 230 } 231