Home | History | Annotate | Line # | Download | only in fpe
fpu_rem.c revision 1.12
      1  1.12   isaki /*	$NetBSD: fpu_rem.c,v 1.12 2013/03/26 10:57:13 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.12   isaki __KERNEL_RCSID(0, "$NetBSD: fpu_rem.c,v 1.12 2013/03/26 10:57:13 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.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.1  briggs  *            3.1 If R = Y, go to Step 9.
     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.1  briggs  *       Step 4.  At this point, R = X - QY = MOD(X,Y). Set
     66   1.1  briggs  *                Last_Subtract := false (used in Step 7 below). If
     67   1.1  briggs  *                MOD is requested, go to Step 6.
     68   1.1  briggs  *
     69   1.1  briggs  *       Step 5.  R = MOD(X,Y), but REM(X,Y) is requested.
     70   1.1  briggs  *            5.1 If R < Y/2, then R = MOD(X,Y) = REM(X,Y). Go to
     71   1.1  briggs  *                Step 6.
     72   1.1  briggs  *            5.2 If R > Y/2, then { set Last_Subtract := true,
     73   1.1  briggs  *                Q := Q + 1, Y := signY*Y }. Go to Step 6.
     74   1.1  briggs  *            5.3 This is the tricky case of R = Y/2. If Q is odd,
     75   1.1  briggs  *                then { Q := Q + 1, signX := -signX }.
     76   1.1  briggs  *
     77   1.1  briggs  *       Step 6.  R := signX*R.
     78   1.1  briggs  *
     79   1.1  briggs  *       Step 7.  If Last_Subtract = true, R := R - Y.
     80   1.1  briggs  *
     81   1.1  briggs  *       Step 8.  Return signQ, last 7 bits of Q, and R as required.
     82   1.1  briggs  *
     83   1.1  briggs  *       Step 9.  At this point, R = 2^(-j)*X - Q Y = Y. Thus,
     84   1.1  briggs  *                X = 2^(j)*(Q+1)Y. set Q := 2^(j)*(Q+1),
     85   1.1  briggs  *                R := 0. Return signQ, last 7 bits of Q, and R.
     86  1.10   isaki  */
     87   1.1  briggs 
     88   1.7     dsl static struct fpn * __fpu_modrem(struct fpemu *fe, int modrem);
     89   1.1  briggs 
     90   1.1  briggs static struct fpn *
     91   1.9     dsl __fpu_modrem(struct fpemu *fe, int modrem)
     92   1.1  briggs {
     93  1.10   isaki 	static struct fpn X, Y;
     94  1.10   isaki 	struct fpn *x, *y, *r;
     95  1.12   isaki 	struct fpn r_bkup;
     96  1.10   isaki 	u_int signX, signY, signQ;
     97  1.10   isaki 	int j, k, l, q;
     98  1.10   isaki 	int Last_Subtract;
     99  1.10   isaki 
    100  1.10   isaki 	CPYFPN(&X, &fe->fe_f1);
    101  1.10   isaki 	CPYFPN(&Y, &fe->fe_f2);
    102  1.10   isaki 	x = &X;
    103  1.10   isaki 	y = &Y;
    104  1.10   isaki 	r = &fe->fe_f2;
    105   1.1  briggs 
    106   1.1  briggs 	/*
    107  1.10   isaki 	 * Step 1
    108   1.1  briggs 	 */
    109  1.10   isaki 	signX = x->fp_sign;
    110  1.10   isaki 	signY = y->fp_sign;
    111  1.10   isaki 	signQ = (signX ^ signY);
    112  1.10   isaki 	x->fp_sign = y->fp_sign = 0;
    113   1.1  briggs 
    114  1.10   isaki 	/*
    115  1.10   isaki 	 * Step 2
    116  1.10   isaki 	 */
    117  1.10   isaki 	l = x->fp_exp - y->fp_exp;
    118  1.10   isaki 	k = 0;
    119  1.10   isaki 	q = 0;
    120  1.10   isaki 	if (l >= 0) {
    121  1.10   isaki 		CPYFPN(r, x);
    122  1.10   isaki 		r->fp_exp -= l;
    123  1.10   isaki 		j = l;
    124  1.10   isaki 
    125  1.10   isaki 		/*
    126  1.10   isaki 		 * Step 3
    127  1.10   isaki 		 */
    128  1.11   isaki 		while (y->fp_exp != r->fp_exp ||
    129  1.11   isaki 		       y->fp_mant[0] != r->fp_mant[0] ||
    130  1.10   isaki 		       y->fp_mant[1] != r->fp_mant[1] ||
    131  1.10   isaki 		       y->fp_mant[2] != r->fp_mant[2]) {
    132  1.10   isaki 
    133  1.10   isaki 			/* Step 3.2 */
    134  1.12   isaki 			CPYFPN(&r_bkup, r);
    135  1.12   isaki 			CPYFPN(&fe->fe_f1, r);
    136  1.12   isaki 			CPYFPN(&fe->fe_f2, y);
    137  1.12   isaki 			fe->fe_f2.fp_sign = 1;
    138  1.12   isaki 			r = fpu_add(fe);
    139  1.12   isaki 			if (r->fp_sign == 0) {
    140  1.10   isaki 				q++;
    141  1.12   isaki 			} else {
    142  1.12   isaki 				CPYFPN(r, &r_bkup);
    143  1.10   isaki 			}
    144  1.10   isaki 
    145  1.10   isaki 			/* Step 3.3 */
    146  1.10   isaki 			if (j == 0)
    147  1.10   isaki 				goto Step4;
    148  1.10   isaki 
    149  1.10   isaki 			/* Step 3.4 */
    150  1.10   isaki 			k++;
    151  1.10   isaki 			j--;
    152  1.10   isaki 			q += q;
    153  1.10   isaki 			r->fp_exp++;
    154  1.10   isaki 		}
    155  1.10   isaki 		/* Step 9 */
    156  1.10   isaki 		goto Step9;
    157   1.1  briggs 	}
    158   1.1  briggs  Step4:
    159  1.10   isaki 	Last_Subtract = 0;
    160  1.10   isaki 	if (modrem == 0)
    161  1.10   isaki 		goto Step6;
    162  1.10   isaki 
    163  1.10   isaki 	/*
    164  1.10   isaki 	 * Step 5
    165  1.10   isaki 	 */
    166  1.10   isaki 	/* Step 5.1 */
    167  1.10   isaki 	if (r->fp_exp + 1 < y->fp_exp ||
    168  1.10   isaki 	    (r->fp_exp + 1 == y->fp_exp &&
    169  1.11   isaki 	     (r->fp_mant[0] < y->fp_mant[0] ||
    170  1.11   isaki 	      r->fp_mant[1] < y->fp_mant[1] ||
    171  1.10   isaki 	      r->fp_mant[2] < y->fp_mant[2]))) {
    172  1.10   isaki 		/* if r < y/2 */
    173  1.10   isaki 		goto Step6;
    174  1.10   isaki 	}
    175  1.10   isaki 	/* Step 5.2 */
    176  1.10   isaki 	if (r->fp_exp + 1 != y->fp_exp ||
    177  1.11   isaki 	    r->fp_mant[0] != y->fp_mant[0] ||
    178  1.11   isaki 	    r->fp_mant[1] != y->fp_mant[1] ||
    179  1.10   isaki 	    r->fp_mant[2] != y->fp_mant[2]) {
    180  1.10   isaki 		/* if (!(r < y/2) && !(r == y/2)) */
    181  1.10   isaki 		Last_Subtract = 1;
    182  1.10   isaki 		q++;
    183  1.10   isaki 		y->fp_sign = signY;
    184  1.10   isaki 	} else {
    185  1.10   isaki 		/* Step 5.3 */
    186  1.10   isaki 		/* r == y/2 */
    187  1.10   isaki 		if (q % 2) {
    188  1.10   isaki 			q++;
    189  1.10   isaki 			signX = !signX;
    190  1.10   isaki 		}
    191   1.1  briggs 	}
    192   1.1  briggs 
    193   1.1  briggs  Step6:
    194  1.10   isaki 	r->fp_sign = signX;
    195   1.1  briggs 
    196  1.10   isaki 	/*
    197  1.10   isaki 	 * Step 7
    198  1.10   isaki 	 */
    199  1.10   isaki 	if (Last_Subtract) {
    200  1.10   isaki 		CPYFPN(&fe->fe_f1, r);
    201  1.10   isaki 		CPYFPN(&fe->fe_f2, y);
    202  1.10   isaki 		fe->fe_f2.fp_sign = !y->fp_sign;
    203  1.10   isaki 		r = fpu_add(fe);
    204  1.10   isaki 	}
    205  1.10   isaki 	/*
    206  1.10   isaki 	 * Step 8
    207  1.10   isaki 	 */
    208  1.10   isaki 	q &= 0x7f;
    209  1.10   isaki 	q |= (signQ << 7);
    210  1.10   isaki 	fe->fe_fpframe->fpf_fpsr =
    211   1.1  briggs 	fe->fe_fpsr =
    212  1.11   isaki 	    (fe->fe_fpsr & ~FPSR_QTT) | (q << 16);
    213  1.10   isaki 	return r;
    214   1.3  briggs 
    215   1.3  briggs  Step9:
    216  1.10   isaki 	fe->fe_f1.fp_class = FPC_ZERO;
    217  1.10   isaki 	q++;
    218  1.10   isaki 	q &= 0x7f;
    219  1.10   isaki 	q |= (signQ << 7);
    220  1.10   isaki 	fe->fe_fpframe->fpf_fpsr =
    221   1.3  briggs 	fe->fe_fpsr =
    222  1.11   isaki 	    (fe->fe_fpsr & ~FPSR_QTT) | (q << 16);
    223  1.10   isaki 	return &fe->fe_f1;
    224   1.1  briggs }
    225   1.1  briggs 
    226   1.1  briggs struct fpn *
    227   1.8     dsl fpu_rem(struct fpemu *fe)
    228   1.1  briggs {
    229  1.10   isaki 	return __fpu_modrem(fe, 1);
    230   1.1  briggs }
    231   1.1  briggs 
    232   1.1  briggs struct fpn *
    233   1.8     dsl fpu_mod(struct fpemu *fe)
    234   1.1  briggs {
    235  1.10   isaki 	return __fpu_modrem(fe, 0);
    236   1.1  briggs }
    237