Home | History | Annotate | Line # | Download | only in fpe
      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