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