Home | History | Annotate | Line # | Download | only in fpe
fpu_rem.c revision 1.16.10.1
      1  1.16.10.1   skrll /*	$NetBSD: fpu_rem.c,v 1.16.10.1 2015/04/06 15:17:58 skrll 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.16.10.1   skrll __KERNEL_RCSID(0, "$NetBSD: fpu_rem.c,v 1.16.10.1 2015/04/06 15:17:58 skrll 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.16.10.1   skrll  *            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.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.16   isaki  *       Step 4.  R := signX*R.
     66       1.16   isaki  *
     67  1.16.10.1   skrll  *       Step 5.  If MOD is requested, go to Step 7.
     68       1.16   isaki  *
     69       1.16   isaki  *       Step 6.  Now, R = MOD(X,Y), convert to REM(X,Y) is requested.
     70       1.16   isaki  *                Do banker's rounding.
     71       1.16   isaki  *                If   abs(R) > Y/2
     72       1.16   isaki  *                 || (abs(R) == Y/2 && Q % 2 == 1) then
     73       1.16   isaki  *                 { Q := Q + 1, R := R - signX * Y }.
     74       1.16   isaki  *
     75       1.16   isaki  *       Step 7.  Return signQ, last 7 bits of Q, and R as required.
     76       1.10   isaki  */
     77        1.1  briggs 
     78       1.14   isaki static struct fpn * __fpu_modrem(struct fpemu *fe, int is_mod);
     79       1.16   isaki static int abscmp3(const struct fpn *a, const struct fpn *b);
     80       1.16   isaki 
     81       1.16   isaki /* Absolute FORTRAN Compare */
     82       1.16   isaki static int
     83       1.16   isaki abscmp3(const struct fpn *a, const struct fpn *b)
     84       1.16   isaki {
     85       1.16   isaki 	int i;
     86       1.16   isaki 
     87       1.16   isaki 	if (a->fp_exp < b->fp_exp) {
     88       1.16   isaki 		return -1;
     89       1.16   isaki 	} else if (a->fp_exp > b->fp_exp) {
     90       1.16   isaki 		return 1;
     91       1.16   isaki 	} else {
     92       1.16   isaki 		for (i = 0; i < 3; i++) {
     93       1.16   isaki 			if (a->fp_mant[i] < b->fp_mant[i])
     94       1.16   isaki 				return -1;
     95       1.16   isaki 			else if (a->fp_mant[i] > b->fp_mant[i])
     96       1.16   isaki 				return 1;
     97       1.16   isaki 		}
     98       1.16   isaki 	}
     99       1.16   isaki 	return 0;
    100       1.16   isaki }
    101        1.1  briggs 
    102        1.1  briggs static struct fpn *
    103       1.14   isaki __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.13   isaki 	uint32_t signX, signY, signQ;
    108       1.10   isaki 	int j, k, l, q;
    109       1.16   isaki 	int cmp;
    110       1.16   isaki 
    111       1.16   isaki 	if (ISNAN(&fe->fe_f1) || ISNAN(&fe->fe_f2))
    112       1.16   isaki 		return fpu_newnan(fe);
    113       1.16   isaki 	if (ISINF(&fe->fe_f1) || ISZERO(&fe->fe_f2))
    114       1.16   isaki 		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.16   isaki 	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.16   isaki 	/* Special treatment that just return input value but Q is necessary */
    132       1.16   isaki 	if (ISZERO(x) || ISINF(y)) {
    133       1.16   isaki 		r = &fe->fe_f1;
    134       1.16   isaki 		goto Step7;
    135       1.16   isaki 	}
    136       1.16   isaki 
    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.15   isaki 	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.16   isaki 		for (;;) {
    151       1.16   isaki 			cmp = abscmp3(r, y);
    152       1.16   isaki 
    153       1.16   isaki 			/* Step 3.1 */
    154       1.16   isaki 			if (cmp == 0)
    155       1.16   isaki 				break;
    156       1.10   isaki 
    157       1.10   isaki 			/* Step 3.2 */
    158       1.16   isaki 			if (cmp > 0) {
    159       1.16   isaki 				CPYFPN(&fe->fe_f1, r);
    160       1.16   isaki 				CPYFPN(&fe->fe_f2, y);
    161       1.16   isaki 				fe->fe_f2.fp_sign = 1;
    162       1.16   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.16   isaki 		/* R == Y */
    177       1.16   isaki 		q++;
    178       1.16   isaki 		r->fp_class = FPC_ZERO;
    179       1.16   isaki 		goto Step7;
    180        1.1  briggs 	}
    181        1.1  briggs  Step4:
    182       1.16   isaki 	r->fp_sign = signX;
    183       1.10   isaki 
    184       1.10   isaki 	/*
    185       1.10   isaki 	 * Step 5
    186       1.10   isaki 	 */
    187       1.16   isaki 	if (is_mod)
    188       1.16   isaki 		goto Step7;
    189        1.1  briggs 
    190       1.10   isaki 	/*
    191       1.16   isaki 	 * Step 6
    192       1.10   isaki 	 */
    193       1.16   isaki 	/* y = y / 2 */
    194       1.16   isaki 	y->fp_exp--;
    195       1.16   isaki 	/* abscmp3 ignore sign */
    196       1.16   isaki 	cmp = abscmp3(r, y);
    197       1.16   isaki 	/* revert y */
    198       1.16   isaki 	y->fp_exp++;
    199       1.16   isaki 
    200       1.16   isaki 	if (cmp > 0 || (cmp == 0 && q % 2)) {
    201       1.16   isaki 		q++;
    202       1.10   isaki 		CPYFPN(&fe->fe_f1, r);
    203       1.10   isaki 		CPYFPN(&fe->fe_f2, y);
    204       1.16   isaki 		fe->fe_f2.fp_sign = !signX;
    205       1.10   isaki 		r = fpu_add(fe);
    206       1.10   isaki 	}
    207       1.16   isaki 
    208       1.10   isaki 	/*
    209       1.16   isaki 	 * Step 7
    210       1.10   isaki 	 */
    211       1.16   isaki  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.14   isaki 	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.14   isaki 	return __fpu_modrem(fe, 1);
    230        1.1  briggs }
    231