xref: /openbsd-src/sys/arch/hppa/spmath/sfrem.c (revision c2feb25228e11b3c285a2d9400475896f2dd1562)
1*c2feb252Smickey /*	$OpenBSD: sfrem.c,v 1.5 2002/05/07 22:19:30 mickey Exp $	*/
28a472b3eSmickey /*
3*c2feb252Smickey   (c) Copyright 1986 HEWLETT-PACKARD COMPANY
4*c2feb252Smickey   To anyone who acknowledges that this file is provided "AS IS"
5*c2feb252Smickey   without any express or implied warranty:
6*c2feb252Smickey       permission to use, copy, modify, and distribute this file
7*c2feb252Smickey   for any purpose is hereby granted without fee, provided that
8*c2feb252Smickey   the above copyright notice and this notice appears in all
9*c2feb252Smickey   copies, and that the name of Hewlett-Packard Company not be
10*c2feb252Smickey   used in advertising or publicity pertaining to distribution
11*c2feb252Smickey   of the software without specific, written prior permission.
12*c2feb252Smickey   Hewlett-Packard Company makes no representations about the
13*c2feb252Smickey   suitability of this software for any purpose.
148a472b3eSmickey */
15*c2feb252Smickey /* @(#)sfrem.c: Revision: 1.7.88.1 Date: 93/12/07 15:07:10 */
168a472b3eSmickey 
17*c2feb252Smickey #include "float.h"
18*c2feb252Smickey #include "sgl_float.h"
198a472b3eSmickey 
208a472b3eSmickey /*
218a472b3eSmickey  *  Single Precision Floating-point Remainder
228a472b3eSmickey  */
23b94afd46Smickey int
sgl_frem(srcptr1,srcptr2,dstptr,status)248a472b3eSmickey sgl_frem(srcptr1,srcptr2,dstptr,status)
258a472b3eSmickey 	sgl_floating_point *srcptr1, *srcptr2, *dstptr;
268a472b3eSmickey 	unsigned int *status;
278a472b3eSmickey {
288a472b3eSmickey 	register unsigned int opnd1, opnd2, result;
298a472b3eSmickey 	register int opnd1_exponent, opnd2_exponent, dest_exponent, stepcount;
30628bbdd4Smickey 	register int roundup = FALSE;
318a472b3eSmickey 
328a472b3eSmickey 	opnd1 = *srcptr1;
338a472b3eSmickey 	opnd2 = *srcptr2;
348a472b3eSmickey 	/*
358a472b3eSmickey 	 * check first operand for NaN's or infinity
368a472b3eSmickey 	 */
378a472b3eSmickey 	if ((opnd1_exponent = Sgl_exponent(opnd1)) == SGL_INFINITY_EXPONENT) {
388a472b3eSmickey 		if (Sgl_iszero_mantissa(opnd1)) {
398a472b3eSmickey 			if (Sgl_isnotnan(opnd2)) {
408a472b3eSmickey 				/* invalid since first operand is infinity */
418a472b3eSmickey 				if (Is_invalidtrap_enabled())
428a472b3eSmickey 					return(INVALIDEXCEPTION);
438a472b3eSmickey 				Set_invalidflag();
448a472b3eSmickey 				Sgl_makequietnan(result);
458a472b3eSmickey 				*dstptr = result;
468a472b3eSmickey 				return(NOEXCEPTION);
478a472b3eSmickey 			}
488a472b3eSmickey 		}
498a472b3eSmickey 		else {
508a472b3eSmickey 			/*
518a472b3eSmickey 			 * is NaN; signaling or quiet?
528a472b3eSmickey 			 */
538a472b3eSmickey 			if (Sgl_isone_signaling(opnd1)) {
548a472b3eSmickey 				/* trap if INVALIDTRAP enabled */
558a472b3eSmickey 				if (Is_invalidtrap_enabled())
568a472b3eSmickey 					return(INVALIDEXCEPTION);
578a472b3eSmickey 				/* make NaN quiet */
588a472b3eSmickey 				Set_invalidflag();
598a472b3eSmickey 				Sgl_set_quiet(opnd1);
608a472b3eSmickey 			}
618a472b3eSmickey 			/*
628a472b3eSmickey 			 * is second operand a signaling NaN?
638a472b3eSmickey 			 */
648a472b3eSmickey 			else if (Sgl_is_signalingnan(opnd2)) {
658a472b3eSmickey 				/* trap if INVALIDTRAP enabled */
668a472b3eSmickey 				if (Is_invalidtrap_enabled())
678a472b3eSmickey 					return(INVALIDEXCEPTION);
688a472b3eSmickey 				/* make NaN quiet */
698a472b3eSmickey 				Set_invalidflag();
708a472b3eSmickey 				Sgl_set_quiet(opnd2);
718a472b3eSmickey 				*dstptr = opnd2;
728a472b3eSmickey 				return(NOEXCEPTION);
738a472b3eSmickey 			}
748a472b3eSmickey 			/*
758a472b3eSmickey 			 * return quiet NaN
768a472b3eSmickey 			 */
778a472b3eSmickey 			*dstptr = opnd1;
788a472b3eSmickey 			return(NOEXCEPTION);
798a472b3eSmickey 		}
808a472b3eSmickey 	}
818a472b3eSmickey 	/*
828a472b3eSmickey 	 * check second operand for NaN's or infinity
838a472b3eSmickey 	 */
848a472b3eSmickey 	if ((opnd2_exponent = Sgl_exponent(opnd2)) == SGL_INFINITY_EXPONENT) {
858a472b3eSmickey 		if (Sgl_iszero_mantissa(opnd2)) {
868a472b3eSmickey 			/*
878a472b3eSmickey 			 * return first operand
888a472b3eSmickey 			 */
898a472b3eSmickey 			*dstptr = opnd1;
908a472b3eSmickey 			return(NOEXCEPTION);
918a472b3eSmickey 		}
928a472b3eSmickey 		/*
938a472b3eSmickey 		 * is NaN; signaling or quiet?
948a472b3eSmickey 		 */
958a472b3eSmickey 		if (Sgl_isone_signaling(opnd2)) {
968a472b3eSmickey 			/* trap if INVALIDTRAP enabled */
978a472b3eSmickey 			if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
988a472b3eSmickey 			/* make NaN quiet */
998a472b3eSmickey 			Set_invalidflag();
1008a472b3eSmickey 			Sgl_set_quiet(opnd2);
1018a472b3eSmickey 		}
1028a472b3eSmickey 		/*
1038a472b3eSmickey 		 * return quiet NaN
1048a472b3eSmickey 		 */
1058a472b3eSmickey 		*dstptr = opnd2;
1068a472b3eSmickey 		return(NOEXCEPTION);
1078a472b3eSmickey 	}
1088a472b3eSmickey 	/*
1098a472b3eSmickey 	 * check second operand for zero
1108a472b3eSmickey 	 */
1118a472b3eSmickey 	if (Sgl_iszero_exponentmantissa(opnd2)) {
1128a472b3eSmickey 		/* invalid since second operand is zero */
1138a472b3eSmickey 		if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
1148a472b3eSmickey 		Set_invalidflag();
1158a472b3eSmickey 		Sgl_makequietnan(result);
1168a472b3eSmickey 		*dstptr = result;
1178a472b3eSmickey 		return(NOEXCEPTION);
1188a472b3eSmickey 	}
1198a472b3eSmickey 
1208a472b3eSmickey 	/*
1218a472b3eSmickey 	 * get sign of result
1228a472b3eSmickey 	 */
1238a472b3eSmickey 	result = opnd1;
1248a472b3eSmickey 
1258a472b3eSmickey 	/*
1268a472b3eSmickey 	 * check for denormalized operands
1278a472b3eSmickey 	 */
1288a472b3eSmickey 	if (opnd1_exponent == 0) {
1298a472b3eSmickey 		/* check for zero */
1308a472b3eSmickey 		if (Sgl_iszero_mantissa(opnd1)) {
1318a472b3eSmickey 			*dstptr = opnd1;
1328a472b3eSmickey 			return(NOEXCEPTION);
1338a472b3eSmickey 		}
1348a472b3eSmickey 		/* normalize, then continue */
1358a472b3eSmickey 		opnd1_exponent = 1;
1368a472b3eSmickey 		Sgl_normalize(opnd1,opnd1_exponent);
1378a472b3eSmickey 	}
1388a472b3eSmickey 	else {
1398a472b3eSmickey 		Sgl_clear_signexponent_set_hidden(opnd1);
1408a472b3eSmickey 	}
1418a472b3eSmickey 	if (opnd2_exponent == 0) {
1428a472b3eSmickey 		/* normalize, then continue */
1438a472b3eSmickey 		opnd2_exponent = 1;
1448a472b3eSmickey 		Sgl_normalize(opnd2,opnd2_exponent);
1458a472b3eSmickey 	}
1468a472b3eSmickey 	else {
1478a472b3eSmickey 		Sgl_clear_signexponent_set_hidden(opnd2);
1488a472b3eSmickey 	}
1498a472b3eSmickey 
1508a472b3eSmickey 	/* find result exponent and divide step loop count */
1518a472b3eSmickey 	dest_exponent = opnd2_exponent - 1;
1528a472b3eSmickey 	stepcount = opnd1_exponent - opnd2_exponent;
1538a472b3eSmickey 
1548a472b3eSmickey 	/*
1558a472b3eSmickey 	 * check for opnd1/opnd2 < 1
1568a472b3eSmickey 	 */
1578a472b3eSmickey 	if (stepcount < 0) {
1588a472b3eSmickey 		/*
1598a472b3eSmickey 		 * check for opnd1/opnd2 > 1/2
1608a472b3eSmickey 		 *
1618a472b3eSmickey 		 * In this case n will round to 1, so
1628a472b3eSmickey 		 *    r = opnd1 - opnd2
1638a472b3eSmickey 		 */
1648a472b3eSmickey 		if (stepcount == -1 && Sgl_isgreaterthan(opnd1,opnd2)) {
1658a472b3eSmickey 			Sgl_all(result) = ~Sgl_all(result);   /* set sign */
1668a472b3eSmickey 			/* align opnd2 with opnd1 */
1678a472b3eSmickey 			Sgl_leftshiftby1(opnd2);
1688a472b3eSmickey 			Sgl_subtract(opnd2,opnd1,opnd2);
1698a472b3eSmickey 			/* now normalize */
1708a472b3eSmickey 			while (Sgl_iszero_hidden(opnd2)) {
1718a472b3eSmickey 				Sgl_leftshiftby1(opnd2);
1728a472b3eSmickey 				dest_exponent--;
1738a472b3eSmickey 			}
1748a472b3eSmickey 			Sgl_set_exponentmantissa(result,opnd2);
1758a472b3eSmickey 			goto testforunderflow;
1768a472b3eSmickey 		}
1778a472b3eSmickey 		/*
1788a472b3eSmickey 		 * opnd1/opnd2 <= 1/2
1798a472b3eSmickey 		 *
1808a472b3eSmickey 		 * In this case n will round to zero, so
1818a472b3eSmickey 		 *    r = opnd1
1828a472b3eSmickey 		 */
1838a472b3eSmickey 		Sgl_set_exponentmantissa(result,opnd1);
1848a472b3eSmickey 		dest_exponent = opnd1_exponent;
1858a472b3eSmickey 		goto testforunderflow;
1868a472b3eSmickey 	}
1878a472b3eSmickey 
1888a472b3eSmickey 	/*
1898a472b3eSmickey 	 * Generate result
1908a472b3eSmickey 	 *
1918a472b3eSmickey 	 * Do iterative subtract until remainder is less than operand 2.
1928a472b3eSmickey 	 */
1938a472b3eSmickey 	while (stepcount-- > 0 && Sgl_all(opnd1)) {
1948a472b3eSmickey 		if (Sgl_isnotlessthan(opnd1,opnd2))
1958a472b3eSmickey 			Sgl_subtract(opnd1,opnd2,opnd1);
1968a472b3eSmickey 		Sgl_leftshiftby1(opnd1);
1978a472b3eSmickey 	}
1988a472b3eSmickey 	/*
1998a472b3eSmickey 	 * Do last subtract, then determine which way to round if remainder
2008a472b3eSmickey 	 * is exactly 1/2 of opnd2
2018a472b3eSmickey 	 */
2028a472b3eSmickey 	if (Sgl_isnotlessthan(opnd1,opnd2)) {
2038a472b3eSmickey 		Sgl_subtract(opnd1,opnd2,opnd1);
2048a472b3eSmickey 		roundup = TRUE;
2058a472b3eSmickey 	}
2068a472b3eSmickey 	if (stepcount > 0 || Sgl_iszero(opnd1)) {
2078a472b3eSmickey 		/* division is exact, remainder is zero */
2088a472b3eSmickey 		Sgl_setzero_exponentmantissa(result);
2098a472b3eSmickey 		*dstptr = result;
2108a472b3eSmickey 		return(NOEXCEPTION);
2118a472b3eSmickey 	}
2128a472b3eSmickey 
2138a472b3eSmickey 	/*
2148a472b3eSmickey 	 * Check for cases where opnd1/opnd2 < n
2158a472b3eSmickey 	 *
2168a472b3eSmickey 	 * In this case the result's sign will be opposite that of
2178a472b3eSmickey 	 * opnd1.  The mantissa also needs some correction.
2188a472b3eSmickey 	 */
2198a472b3eSmickey 	Sgl_leftshiftby1(opnd1);
2208a472b3eSmickey 	if (Sgl_isgreaterthan(opnd1,opnd2)) {
2218a472b3eSmickey 		Sgl_invert_sign(result);
2228a472b3eSmickey 		Sgl_subtract((opnd2<<1),opnd1,opnd1);
2238a472b3eSmickey 	}
2248a472b3eSmickey 	/* check for remainder being exactly 1/2 of opnd2 */
2258a472b3eSmickey 	else if (Sgl_isequal(opnd1,opnd2) && roundup) {
2268a472b3eSmickey 		Sgl_invert_sign(result);
2278a472b3eSmickey 	}
2288a472b3eSmickey 
2298a472b3eSmickey 	/* normalize result's mantissa */
2308a472b3eSmickey 	while (Sgl_iszero_hidden(opnd1)) {
2318a472b3eSmickey 		dest_exponent--;
2328a472b3eSmickey 		Sgl_leftshiftby1(opnd1);
2338a472b3eSmickey 	}
2348a472b3eSmickey 	Sgl_set_exponentmantissa(result,opnd1);
2358a472b3eSmickey 
2368a472b3eSmickey 	/*
2378a472b3eSmickey 	 * Test for underflow
2388a472b3eSmickey 	 */
2398a472b3eSmickey     testforunderflow:
2408a472b3eSmickey 	if (dest_exponent <= 0) {
2418a472b3eSmickey 		/* trap if UNDERFLOWTRAP enabled */
2428a472b3eSmickey 		if (Is_underflowtrap_enabled()) {
2438a472b3eSmickey 			/*
2448a472b3eSmickey 			 * Adjust bias of result
2458a472b3eSmickey 			 */
2468a472b3eSmickey 			Sgl_setwrapped_exponent(result,dest_exponent,unfl);
2478a472b3eSmickey 			*dstptr = result;
2488a472b3eSmickey 			/* frem is always exact */
2498a472b3eSmickey 			return(UNDERFLOWEXCEPTION);
2508a472b3eSmickey 		}
2518a472b3eSmickey 		/*
2528a472b3eSmickey 		 * denormalize result or set to signed zero
2538a472b3eSmickey 		 */
2548a472b3eSmickey 		if (dest_exponent >= (1 - SGL_P)) {
2558a472b3eSmickey 			Sgl_rightshift_exponentmantissa(result,1-dest_exponent);
256628bbdd4Smickey 		} else {
2578a472b3eSmickey 			Sgl_setzero_exponentmantissa(result);
2588a472b3eSmickey 		}
2598a472b3eSmickey 	}
2608a472b3eSmickey 	else Sgl_set_exponent(result,dest_exponent);
2618a472b3eSmickey 	*dstptr = result;
2628a472b3eSmickey 	return(NOEXCEPTION);
2638a472b3eSmickey }
264