Home | History | Annotate | Line # | Download | only in fpe
fpu_rem.c revision 1.11.2.1
      1  1.11.2.1    yamt /*	$NetBSD: fpu_rem.c,v 1.11.2.1 2014/05/22 11:39:55 yamt 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.11.2.1    yamt __KERNEL_RCSID(0, "$NetBSD: fpu_rem.c,v 1.11.2.1 2014/05/22 11:39:55 yamt 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.1  briggs  *       Step 2.  Set L := expo(X)-expo(Y), k := 0, 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.11.2.1    yamt  *            3.1 If R = Y, then { Q := Q + 1, R := 0, go to Step 8. }
     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.1  briggs  *            3.4 k := k + 1, j := j - 1, Q := 2Q, R := 2R. Go to
     63       1.1  briggs  *                Step 3.1.
     64       1.1  briggs  *
     65  1.11.2.1    yamt  *       Step 4.  R := signX*R.
     66  1.11.2.1    yamt  *
     67  1.11.2.1    yamt  *       Step 5.  If MOD is requested, go to Step .
     68  1.11.2.1    yamt  *
     69  1.11.2.1    yamt  *       Step 6.  Now, R = MOD(X,Y), convert to REM(X,Y) is requested.
     70  1.11.2.1    yamt  *                Do banker's rounding.
     71  1.11.2.1    yamt  *                If   abs(R) > Y/2
     72  1.11.2.1    yamt  *                 || (abs(R) == Y/2 && Q % 2 == 1) then
     73  1.11.2.1    yamt  *                 { Q := Q + 1, R := R - signX * Y }.
     74  1.11.2.1    yamt  *
     75  1.11.2.1    yamt  *       Step 7.  Return signQ, last 7 bits of Q, and R as required.
     76      1.10   isaki  */
     77       1.1  briggs 
     78  1.11.2.1    yamt static struct fpn * __fpu_modrem(struct fpemu *fe, int is_mod);
     79  1.11.2.1    yamt static int abscmp3(const struct fpn *a, const struct fpn *b);
     80  1.11.2.1    yamt 
     81  1.11.2.1    yamt /* Absolute FORTRAN Compare */
     82  1.11.2.1    yamt static int
     83  1.11.2.1    yamt abscmp3(const struct fpn *a, const struct fpn *b)
     84  1.11.2.1    yamt {
     85  1.11.2.1    yamt 	int i;
     86  1.11.2.1    yamt 
     87  1.11.2.1    yamt 	if (a->fp_exp < b->fp_exp) {
     88  1.11.2.1    yamt 		return -1;
     89  1.11.2.1    yamt 	} else if (a->fp_exp > b->fp_exp) {
     90  1.11.2.1    yamt 		return 1;
     91  1.11.2.1    yamt 	} else {
     92  1.11.2.1    yamt 		for (i = 0; i < 3; i++) {
     93  1.11.2.1    yamt 			if (a->fp_mant[i] < b->fp_mant[i])
     94  1.11.2.1    yamt 				return -1;
     95  1.11.2.1    yamt 			else if (a->fp_mant[i] > b->fp_mant[i])
     96  1.11.2.1    yamt 				return 1;
     97  1.11.2.1    yamt 		}
     98  1.11.2.1    yamt 	}
     99  1.11.2.1    yamt 	return 0;
    100  1.11.2.1    yamt }
    101       1.1  briggs 
    102       1.1  briggs static struct fpn *
    103  1.11.2.1    yamt __fpu_modrem(struct fpemu *fe, int is_mod)
    104       1.1  briggs {
    105      1.10   isaki 	static struct fpn X, Y;
    106      1.10   isaki 	struct fpn *x, *y, *r;
    107  1.11.2.1    yamt 	uint32_t signX, signY, signQ;
    108      1.10   isaki 	int j, k, l, q;
    109  1.11.2.1    yamt 	int cmp;
    110  1.11.2.1    yamt 
    111  1.11.2.1    yamt 	if (ISNAN(&fe->fe_f1) || ISNAN(&fe->fe_f2))
    112  1.11.2.1    yamt 		return fpu_newnan(fe);
    113  1.11.2.1    yamt 	if (ISINF(&fe->fe_f1) || ISZERO(&fe->fe_f2))
    114  1.11.2.1    yamt 		return fpu_newnan(fe);
    115      1.10   isaki 
    116      1.10   isaki 	CPYFPN(&X, &fe->fe_f1);
    117      1.10   isaki 	CPYFPN(&Y, &fe->fe_f2);
    118      1.10   isaki 	x = &X;
    119      1.10   isaki 	y = &Y;
    120  1.11.2.1    yamt 	q = 0;
    121      1.10   isaki 	r = &fe->fe_f2;
    122       1.1  briggs 
    123       1.1  briggs 	/*
    124      1.10   isaki 	 * Step 1
    125       1.1  briggs 	 */
    126      1.10   isaki 	signX = x->fp_sign;
    127      1.10   isaki 	signY = y->fp_sign;
    128      1.10   isaki 	signQ = (signX ^ signY);
    129      1.10   isaki 	x->fp_sign = y->fp_sign = 0;
    130       1.1  briggs 
    131  1.11.2.1    yamt 	/* Special treatment that just return input value but Q is necessary */
    132  1.11.2.1    yamt 	if (ISZERO(x) || ISINF(y)) {
    133  1.11.2.1    yamt 		r = &fe->fe_f1;
    134  1.11.2.1    yamt 		goto Step7;
    135  1.11.2.1    yamt 	}
    136  1.11.2.1    yamt 
    137      1.10   isaki 	/*
    138      1.10   isaki 	 * Step 2
    139      1.10   isaki 	 */
    140      1.10   isaki 	l = x->fp_exp - y->fp_exp;
    141      1.10   isaki 	k = 0;
    142  1.11.2.1    yamt 	CPYFPN(r, x);
    143      1.10   isaki 	if (l >= 0) {
    144      1.10   isaki 		r->fp_exp -= l;
    145      1.10   isaki 		j = l;
    146      1.10   isaki 
    147      1.10   isaki 		/*
    148      1.10   isaki 		 * Step 3
    149      1.10   isaki 		 */
    150  1.11.2.1    yamt 		for (;;) {
    151  1.11.2.1    yamt 			cmp = abscmp3(r, y);
    152  1.11.2.1    yamt 
    153  1.11.2.1    yamt 			/* Step 3.1 */
    154  1.11.2.1    yamt 			if (cmp == 0)
    155  1.11.2.1    yamt 				break;
    156      1.10   isaki 
    157      1.10   isaki 			/* Step 3.2 */
    158  1.11.2.1    yamt 			if (cmp > 0) {
    159      1.10   isaki 				CPYFPN(&fe->fe_f1, r);
    160      1.10   isaki 				CPYFPN(&fe->fe_f2, y);
    161      1.10   isaki 				fe->fe_f2.fp_sign = 1;
    162      1.10   isaki 				r = fpu_add(fe);
    163      1.10   isaki 				q++;
    164      1.10   isaki 			}
    165      1.10   isaki 
    166      1.10   isaki 			/* Step 3.3 */
    167      1.10   isaki 			if (j == 0)
    168      1.10   isaki 				goto Step4;
    169      1.10   isaki 
    170      1.10   isaki 			/* Step 3.4 */
    171      1.10   isaki 			k++;
    172      1.10   isaki 			j--;
    173      1.10   isaki 			q += q;
    174      1.10   isaki 			r->fp_exp++;
    175      1.10   isaki 		}
    176  1.11.2.1    yamt 		/* R == Y */
    177  1.11.2.1    yamt 		q++;
    178  1.11.2.1    yamt 		r->fp_class = FPC_ZERO;
    179  1.11.2.1    yamt 		goto Step7;
    180       1.1  briggs 	}
    181       1.1  briggs  Step4:
    182  1.11.2.1    yamt 	r->fp_sign = signX;
    183      1.10   isaki 
    184      1.10   isaki 	/*
    185      1.10   isaki 	 * Step 5
    186      1.10   isaki 	 */
    187  1.11.2.1    yamt 	if (is_mod)
    188  1.11.2.1    yamt 		goto Step7;
    189       1.1  briggs 
    190      1.10   isaki 	/*
    191  1.11.2.1    yamt 	 * Step 6
    192      1.10   isaki 	 */
    193  1.11.2.1    yamt 	/* y = y / 2 */
    194  1.11.2.1    yamt 	y->fp_exp--;
    195  1.11.2.1    yamt 	/* abscmp3 ignore sign */
    196  1.11.2.1    yamt 	cmp = abscmp3(r, y);
    197  1.11.2.1    yamt 	/* revert y */
    198  1.11.2.1    yamt 	y->fp_exp++;
    199  1.11.2.1    yamt 
    200  1.11.2.1    yamt 	if (cmp > 0 || (cmp == 0 && q % 2)) {
    201  1.11.2.1    yamt 		q++;
    202      1.10   isaki 		CPYFPN(&fe->fe_f1, r);
    203      1.10   isaki 		CPYFPN(&fe->fe_f2, y);
    204  1.11.2.1    yamt 		fe->fe_f2.fp_sign = !signX;
    205      1.10   isaki 		r = fpu_add(fe);
    206      1.10   isaki 	}
    207  1.11.2.1    yamt 
    208      1.10   isaki 	/*
    209  1.11.2.1    yamt 	 * Step 7
    210      1.10   isaki 	 */
    211  1.11.2.1    yamt  Step7:
    212      1.10   isaki 	q &= 0x7f;
    213      1.10   isaki 	q |= (signQ << 7);
    214      1.10   isaki 	fe->fe_fpframe->fpf_fpsr =
    215       1.1  briggs 	fe->fe_fpsr =
    216      1.11   isaki 	    (fe->fe_fpsr & ~FPSR_QTT) | (q << 16);
    217      1.10   isaki 	return r;
    218       1.1  briggs }
    219       1.1  briggs 
    220       1.1  briggs struct fpn *
    221       1.8     dsl fpu_rem(struct fpemu *fe)
    222       1.1  briggs {
    223  1.11.2.1    yamt 	return __fpu_modrem(fe, 0);
    224       1.1  briggs }
    225       1.1  briggs 
    226       1.1  briggs struct fpn *
    227       1.8     dsl fpu_mod(struct fpemu *fe)
    228       1.1  briggs {
    229  1.11.2.1    yamt 	return __fpu_modrem(fe, 1);
    230       1.1  briggs }
    231