Home | History | Annotate | Line # | Download | only in spmath
frnd.c revision 1.4.80.1
      1  1.4.80.1      yamt /*	$NetBSD: frnd.c,v 1.4.80.1 2012/04/17 00:06:27 yamt Exp $	*/
      2       1.1  fredette 
      3       1.1  fredette /*	$OpenBSD: frnd.c,v 1.5 2001/03/29 03:58:18 mickey Exp $	*/
      4       1.1  fredette 
      5       1.1  fredette /*
      6       1.1  fredette  * Copyright 1996 1995 by Open Software Foundation, Inc.
      7       1.1  fredette  *              All Rights Reserved
      8       1.1  fredette  *
      9       1.1  fredette  * Permission to use, copy, modify, and distribute this software and
     10       1.1  fredette  * its documentation for any purpose and without fee is hereby granted,
     11       1.1  fredette  * provided that the above copyright notice appears in all copies and
     12       1.1  fredette  * that both the copyright notice and this permission notice appear in
     13       1.1  fredette  * supporting documentation.
     14       1.1  fredette  *
     15       1.1  fredette  * OSF DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE
     16       1.1  fredette  * INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
     17       1.1  fredette  * FOR A PARTICULAR PURPOSE.
     18       1.1  fredette  *
     19       1.1  fredette  * IN NO EVENT SHALL OSF BE LIABLE FOR ANY SPECIAL, INDIRECT, OR
     20       1.1  fredette  * CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM
     21       1.1  fredette  * LOSS OF USE, DATA OR PROFITS, WHETHER IN ACTION OF CONTRACT,
     22       1.1  fredette  * NEGLIGENCE, OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION
     23       1.1  fredette  * WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
     24       1.1  fredette  *
     25       1.1  fredette  */
     26       1.1  fredette /*
     27       1.1  fredette  * pmk1.1
     28       1.1  fredette  */
     29       1.1  fredette /*
     30       1.1  fredette  * (c) Copyright 1986 HEWLETT-PACKARD COMPANY
     31       1.1  fredette  *
     32       1.1  fredette  * To anyone who acknowledges that this file is provided "AS IS"
     33       1.1  fredette  * without any express or implied warranty:
     34       1.1  fredette  *     permission to use, copy, modify, and distribute this file
     35       1.1  fredette  * for any purpose is hereby granted without fee, provided that
     36       1.1  fredette  * the above copyright notice and this notice appears in all
     37       1.1  fredette  * copies, and that the name of Hewlett-Packard Company not be
     38       1.1  fredette  * used in advertising or publicity pertaining to distribution
     39       1.1  fredette  * of the software without specific, written prior permission.
     40       1.1  fredette  * Hewlett-Packard Company makes no representations about the
     41       1.1  fredette  * suitability of this software for any purpose.
     42       1.1  fredette  */
     43       1.2     lukem 
     44       1.2     lukem #include <sys/cdefs.h>
     45  1.4.80.1      yamt __KERNEL_RCSID(0, "$NetBSD: frnd.c,v 1.4.80.1 2012/04/17 00:06:27 yamt Exp $");
     46       1.1  fredette 
     47       1.1  fredette #include "../spmath/float.h"
     48       1.1  fredette #include "../spmath/sgl_float.h"
     49       1.1  fredette #include "../spmath/dbl_float.h"
     50       1.1  fredette #include "../spmath/quad_float.h"
     51       1.1  fredette #include "../spmath/cnv_float.h"
     52       1.1  fredette 
     53       1.1  fredette /*
     54       1.1  fredette  *  Single Floating-point Round to Integer
     55       1.1  fredette  */
     56       1.1  fredette 
     57       1.1  fredette /*ARGSUSED*/
     58       1.1  fredette int
     59  1.4.80.1      yamt sgl_frnd(sgl_floating_point *srcptr, sgl_floating_point *dstptr,
     60  1.4.80.1      yamt     unsigned int *status)
     61       1.1  fredette {
     62       1.1  fredette 	register unsigned int src, result;
     63       1.1  fredette 	register int src_exponent;
     64       1.4   thorpej 	register int inexact = false;
     65       1.1  fredette 
     66       1.1  fredette 	src = *srcptr;
     67       1.1  fredette 	/*
     68       1.1  fredette 	 * check source operand for NaN or infinity
     69       1.1  fredette 	 */
     70       1.1  fredette 	if ((src_exponent = Sgl_exponent(src)) == SGL_INFINITY_EXPONENT) {
     71       1.1  fredette 		/*
     72       1.1  fredette 		 * is signaling NaN?
     73       1.1  fredette 		 */
     74       1.1  fredette 		if (Sgl_isone_signaling(src)) {
     75       1.1  fredette 			/* trap if INVALIDTRAP enabled */
     76       1.1  fredette 			if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
     77       1.1  fredette 			/* make NaN quiet */
     78       1.1  fredette 			Set_invalidflag();
     79       1.1  fredette 			Sgl_set_quiet(src);
     80       1.1  fredette 		}
     81       1.1  fredette 		/*
     82       1.1  fredette 		 * return quiet NaN or infinity
     83       1.1  fredette 		 */
     84       1.1  fredette 		*dstptr = src;
     85       1.1  fredette 		return(NOEXCEPTION);
     86       1.1  fredette 	}
     87       1.1  fredette 	/*
     88       1.1  fredette 	 * Need to round?
     89       1.1  fredette 	 */
     90       1.1  fredette 	if ((src_exponent -= SGL_BIAS) >= SGL_P - 1) {
     91       1.1  fredette 		*dstptr = src;
     92       1.1  fredette 		return(NOEXCEPTION);
     93       1.1  fredette 	}
     94       1.1  fredette 	/*
     95       1.1  fredette 	 * Generate result
     96       1.1  fredette 	 */
     97       1.1  fredette 	if (src_exponent >= 0) {
     98       1.1  fredette 		Sgl_clear_exponent_set_hidden(src);
     99       1.1  fredette 		result = src;
    100       1.1  fredette 		Sgl_rightshift(result,(SGL_P-1) - (src_exponent));
    101       1.1  fredette 		/* check for inexact */
    102       1.1  fredette 		if (Sgl_isinexact_to_fix(src,src_exponent)) {
    103       1.4   thorpej 			inexact = true;
    104       1.1  fredette 			/*  round result  */
    105       1.1  fredette 			switch (Rounding_mode()) {
    106       1.1  fredette 			case ROUNDPLUS:
    107       1.1  fredette 			     if (Sgl_iszero_sign(src)) Sgl_increment(result);
    108       1.1  fredette 			     break;
    109       1.1  fredette 			case ROUNDMINUS:
    110       1.1  fredette 			     if (Sgl_isone_sign(src)) Sgl_increment(result);
    111       1.1  fredette 			     break;
    112       1.1  fredette 			case ROUNDNEAREST:
    113       1.1  fredette 			     if (Sgl_isone_roundbit(src,src_exponent))
    114       1.1  fredette 				if (Sgl_isone_stickybit(src,src_exponent)
    115       1.1  fredette 				|| (Sgl_isone_lowmantissa(result)))
    116       1.1  fredette 					Sgl_increment(result);
    117       1.1  fredette 			}
    118       1.1  fredette 		}
    119       1.1  fredette 		Sgl_leftshift(result,(SGL_P-1) - (src_exponent));
    120       1.1  fredette 		if (Sgl_isone_hiddenoverflow(result))
    121       1.1  fredette 			Sgl_set_exponent(result,src_exponent + (SGL_BIAS+1));
    122       1.1  fredette 		else Sgl_set_exponent(result,src_exponent + SGL_BIAS);
    123       1.1  fredette 	}
    124       1.1  fredette 	else {
    125       1.1  fredette 		result = src;		/* set sign */
    126       1.1  fredette 		Sgl_setzero_exponentmantissa(result);
    127       1.1  fredette 		/* check for inexact */
    128       1.1  fredette 		if (Sgl_isnotzero_exponentmantissa(src)) {
    129       1.4   thorpej 			inexact = true;
    130       1.1  fredette 			/*  round result  */
    131       1.1  fredette 			switch (Rounding_mode()) {
    132       1.1  fredette 			case ROUNDPLUS:
    133       1.1  fredette 			     if (Sgl_iszero_sign(src))
    134       1.1  fredette 				Sgl_set_exponent(result,SGL_BIAS);
    135       1.1  fredette 			     break;
    136       1.1  fredette 			case ROUNDMINUS:
    137       1.1  fredette 			     if (Sgl_isone_sign(src))
    138       1.1  fredette 				Sgl_set_exponent(result,SGL_BIAS);
    139       1.1  fredette 			     break;
    140       1.1  fredette 			case ROUNDNEAREST:
    141       1.1  fredette 			     if (src_exponent == -1)
    142       1.1  fredette 				if (Sgl_isnotzero_mantissa(src))
    143       1.1  fredette 				   Sgl_set_exponent(result,SGL_BIAS);
    144       1.1  fredette 			}
    145       1.1  fredette 		}
    146       1.1  fredette 	}
    147       1.1  fredette 	*dstptr = result;
    148       1.1  fredette 	if (inexact) {
    149       1.1  fredette 		if (Is_inexacttrap_enabled()) return(INEXACTEXCEPTION);
    150       1.1  fredette 		else Set_inexactflag();
    151       1.1  fredette 	}
    152       1.1  fredette 	return(NOEXCEPTION);
    153       1.1  fredette }
    154       1.1  fredette 
    155       1.1  fredette /*
    156       1.1  fredette  *  Double Floating-point Round to Integer
    157       1.1  fredette  */
    158       1.1  fredette 
    159       1.1  fredette /*ARGSUSED*/
    160       1.1  fredette int
    161  1.4.80.1      yamt dbl_frnd(dbl_floating_point *srcptr, dbl_floating_point *dstptr,
    162  1.4.80.1      yamt     unsigned int *status)
    163       1.1  fredette {
    164       1.1  fredette 	register unsigned int srcp1, srcp2, resultp1, resultp2;
    165       1.1  fredette 	register int src_exponent;
    166       1.4   thorpej 	register int inexact = false;
    167       1.1  fredette 
    168       1.1  fredette 	Dbl_copyfromptr(srcptr,srcp1,srcp2);
    169       1.1  fredette 	/*
    170       1.1  fredette 	 * check source operand for NaN or infinity
    171       1.1  fredette 	 */
    172       1.1  fredette 	if ((src_exponent = Dbl_exponent(srcp1)) == DBL_INFINITY_EXPONENT) {
    173       1.1  fredette 		/*
    174       1.1  fredette 		 * is signaling NaN?
    175       1.1  fredette 		 */
    176       1.1  fredette 		if (Dbl_isone_signaling(srcp1)) {
    177       1.1  fredette 			/* trap if INVALIDTRAP enabled */
    178       1.1  fredette 			if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
    179       1.1  fredette 			/* make NaN quiet */
    180       1.1  fredette 			Set_invalidflag();
    181       1.1  fredette 			Dbl_set_quiet(srcp1);
    182       1.1  fredette 		}
    183       1.1  fredette 		/*
    184       1.1  fredette 		 * return quiet NaN or infinity
    185       1.1  fredette 		 */
    186       1.1  fredette 		Dbl_copytoptr(srcp1,srcp2,dstptr);
    187       1.1  fredette 		return(NOEXCEPTION);
    188       1.1  fredette 	}
    189       1.1  fredette 	/*
    190       1.1  fredette 	 * Need to round?
    191       1.1  fredette 	 */
    192       1.1  fredette 	if ((src_exponent -= DBL_BIAS) >= DBL_P - 1) {
    193       1.1  fredette 		Dbl_copytoptr(srcp1,srcp2,dstptr);
    194       1.1  fredette 		return(NOEXCEPTION);
    195       1.1  fredette 	}
    196       1.1  fredette 	/*
    197       1.1  fredette 	 * Generate result
    198       1.1  fredette 	 */
    199       1.1  fredette 	if (src_exponent >= 0) {
    200       1.1  fredette 		Dbl_clear_exponent_set_hidden(srcp1);
    201       1.1  fredette 		resultp1 = srcp1;
    202       1.1  fredette 		resultp2 = srcp2;
    203       1.1  fredette 		Dbl_rightshift(resultp1,resultp2,(DBL_P-1) - (src_exponent));
    204       1.1  fredette 		/* check for inexact */
    205       1.1  fredette 		if (Dbl_isinexact_to_fix(srcp1,srcp2,src_exponent)) {
    206       1.4   thorpej 			inexact = true;
    207       1.1  fredette 			/*  round result  */
    208       1.1  fredette 			switch (Rounding_mode()) {
    209       1.1  fredette 			case ROUNDPLUS:
    210       1.1  fredette 			     if (Dbl_iszero_sign(srcp1))
    211       1.1  fredette 				Dbl_increment(resultp1,resultp2);
    212       1.1  fredette 			     break;
    213       1.1  fredette 			case ROUNDMINUS:
    214       1.1  fredette 			     if (Dbl_isone_sign(srcp1))
    215       1.1  fredette 				Dbl_increment(resultp1,resultp2);
    216       1.1  fredette 			     break;
    217       1.1  fredette 			case ROUNDNEAREST:
    218       1.1  fredette 			     if (Dbl_isone_roundbit(srcp1,srcp2,src_exponent))
    219       1.1  fredette 			      if (Dbl_isone_stickybit(srcp1,srcp2,src_exponent)
    220       1.1  fredette 				  || (Dbl_isone_lowmantissap2(resultp2)))
    221       1.1  fredette 					Dbl_increment(resultp1,resultp2);
    222       1.1  fredette 			}
    223       1.1  fredette 		}
    224       1.1  fredette 		Dbl_leftshift(resultp1,resultp2,(DBL_P-1) - (src_exponent));
    225       1.1  fredette 		if (Dbl_isone_hiddenoverflow(resultp1))
    226       1.1  fredette 			Dbl_set_exponent(resultp1,src_exponent + (DBL_BIAS+1));
    227       1.1  fredette 		else Dbl_set_exponent(resultp1,src_exponent + DBL_BIAS);
    228       1.1  fredette 	}
    229       1.1  fredette 	else {
    230       1.1  fredette 		resultp1 = srcp1;  /* set sign */
    231       1.1  fredette 		Dbl_setzero_exponentmantissa(resultp1,resultp2);
    232       1.1  fredette 		/* check for inexact */
    233       1.1  fredette 		if (Dbl_isnotzero_exponentmantissa(srcp1,srcp2)) {
    234       1.4   thorpej 			inexact = true;
    235       1.1  fredette 			/*  round result  */
    236       1.1  fredette 			switch (Rounding_mode()) {
    237       1.1  fredette 			case ROUNDPLUS:
    238       1.1  fredette 			     if (Dbl_iszero_sign(srcp1))
    239       1.1  fredette 				Dbl_set_exponent(resultp1,DBL_BIAS);
    240       1.1  fredette 			     break;
    241       1.1  fredette 			case ROUNDMINUS:
    242       1.1  fredette 			     if (Dbl_isone_sign(srcp1))
    243       1.1  fredette 				Dbl_set_exponent(resultp1,DBL_BIAS);
    244       1.1  fredette 			     break;
    245       1.1  fredette 			case ROUNDNEAREST:
    246       1.1  fredette 			     if (src_exponent == -1)
    247       1.1  fredette 				if (Dbl_isnotzero_mantissa(srcp1,srcp2))
    248       1.1  fredette 				   Dbl_set_exponent(resultp1,DBL_BIAS);
    249       1.1  fredette 			}
    250       1.1  fredette 		}
    251       1.1  fredette 	}
    252       1.1  fredette 	Dbl_copytoptr(resultp1,resultp2,dstptr);
    253       1.1  fredette 	if (inexact) {
    254       1.1  fredette 		if (Is_inexacttrap_enabled()) return(INEXACTEXCEPTION);
    255       1.1  fredette 		else Set_inexactflag();
    256       1.1  fredette 	}
    257       1.1  fredette 	return(NOEXCEPTION);
    258       1.1  fredette }
    259       1.1  fredette 
    260       1.1  fredette /*ARGSUSED*/
    261       1.1  fredette int
    262  1.4.80.1      yamt quad_frnd(quad_floating_point *srcptr, quad_floating_point *dstptr,
    263  1.4.80.1      yamt     unsigned int *status)
    264       1.1  fredette {
    265       1.1  fredette 	return(UNIMPLEMENTEDEXCEPTION);
    266       1.1  fredette }
    267       1.1  fredette 
    268