xref: /openbmc/linux/arch/parisc/math-emu/dfrem.c (revision 660662f8)
1660662f8SThomas Gleixner // SPDX-License-Identifier: GPL-2.0-or-later
21da177e4SLinus Torvalds /*
31da177e4SLinus Torvalds  * Linux/PA-RISC Project (http://www.parisc-linux.org/)
41da177e4SLinus Torvalds  *
51da177e4SLinus Torvalds  * Floating-point emulation code
61da177e4SLinus Torvalds  *  Copyright (C) 2001 Hewlett-Packard (Paul Bame) <bame@debian.org>
71da177e4SLinus Torvalds  */
81da177e4SLinus Torvalds /*
91da177e4SLinus Torvalds  * BEGIN_DESC
101da177e4SLinus Torvalds  *
111da177e4SLinus Torvalds  *  File:
121da177e4SLinus Torvalds  *	@(#)	pa/spmath/dfrem.c		$Revision: 1.1 $
131da177e4SLinus Torvalds  *
141da177e4SLinus Torvalds  *  Purpose:
151da177e4SLinus Torvalds  *	Double Precision Floating-point Remainder
161da177e4SLinus Torvalds  *
171da177e4SLinus Torvalds  *  External Interfaces:
181da177e4SLinus Torvalds  *	dbl_frem(srcptr1,srcptr2,dstptr,status)
191da177e4SLinus Torvalds  *
201da177e4SLinus Torvalds  *  Internal Interfaces:
211da177e4SLinus Torvalds  *
221da177e4SLinus Torvalds  *  Theory:
231da177e4SLinus Torvalds  *	<<please update with a overview of the operation of this file>>
241da177e4SLinus Torvalds  *
251da177e4SLinus Torvalds  * END_DESC
261da177e4SLinus Torvalds */
271da177e4SLinus Torvalds 
281da177e4SLinus Torvalds 
291da177e4SLinus Torvalds 
301da177e4SLinus Torvalds #include "float.h"
311da177e4SLinus Torvalds #include "dbl_float.h"
321da177e4SLinus Torvalds 
331da177e4SLinus Torvalds /*
341da177e4SLinus Torvalds  *  Double Precision Floating-point Remainder
351da177e4SLinus Torvalds  */
361da177e4SLinus Torvalds 
371da177e4SLinus Torvalds int
dbl_frem(dbl_floating_point * srcptr1,dbl_floating_point * srcptr2,dbl_floating_point * dstptr,unsigned int * status)381da177e4SLinus Torvalds dbl_frem (dbl_floating_point * srcptr1, dbl_floating_point * srcptr2,
391da177e4SLinus Torvalds 	  dbl_floating_point * dstptr, unsigned int *status)
401da177e4SLinus Torvalds {
411da177e4SLinus Torvalds 	register unsigned int opnd1p1, opnd1p2, opnd2p1, opnd2p2;
421da177e4SLinus Torvalds 	register unsigned int resultp1, resultp2;
431da177e4SLinus Torvalds 	register int opnd1_exponent, opnd2_exponent, dest_exponent, stepcount;
441da177e4SLinus Torvalds 	register boolean roundup = FALSE;
451da177e4SLinus Torvalds 
461da177e4SLinus Torvalds 	Dbl_copyfromptr(srcptr1,opnd1p1,opnd1p2);
471da177e4SLinus Torvalds 	Dbl_copyfromptr(srcptr2,opnd2p1,opnd2p2);
481da177e4SLinus Torvalds 	/*
491da177e4SLinus Torvalds 	 * check first operand for NaN's or infinity
501da177e4SLinus Torvalds 	 */
511da177e4SLinus Torvalds 	if ((opnd1_exponent = Dbl_exponent(opnd1p1)) == DBL_INFINITY_EXPONENT) {
521da177e4SLinus Torvalds 		if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) {
531da177e4SLinus Torvalds 			if (Dbl_isnotnan(opnd2p1,opnd2p2)) {
541da177e4SLinus Torvalds 				/* invalid since first operand is infinity */
551da177e4SLinus Torvalds 				if (Is_invalidtrap_enabled())
561da177e4SLinus Torvalds                                 	return(INVALIDEXCEPTION);
571da177e4SLinus Torvalds                                 Set_invalidflag();
581da177e4SLinus Torvalds                                 Dbl_makequietnan(resultp1,resultp2);
591da177e4SLinus Torvalds 				Dbl_copytoptr(resultp1,resultp2,dstptr);
601da177e4SLinus Torvalds 				return(NOEXCEPTION);
611da177e4SLinus Torvalds 			}
621da177e4SLinus Torvalds 		}
631da177e4SLinus Torvalds 		else {
641da177e4SLinus Torvalds                 	/*
651da177e4SLinus Torvalds                  	 * is NaN; signaling or quiet?
661da177e4SLinus Torvalds                  	 */
671da177e4SLinus Torvalds                 	if (Dbl_isone_signaling(opnd1p1)) {
681da177e4SLinus Torvalds                         	/* trap if INVALIDTRAP enabled */
691da177e4SLinus Torvalds                         	if (Is_invalidtrap_enabled())
701da177e4SLinus Torvalds                             		return(INVALIDEXCEPTION);
711da177e4SLinus Torvalds                         	/* make NaN quiet */
721da177e4SLinus Torvalds                         	Set_invalidflag();
731da177e4SLinus Torvalds                         	Dbl_set_quiet(opnd1p1);
741da177e4SLinus Torvalds                 	}
751da177e4SLinus Torvalds 			/*
761da177e4SLinus Torvalds 			 * is second operand a signaling NaN?
771da177e4SLinus Torvalds 			 */
781da177e4SLinus Torvalds 			else if (Dbl_is_signalingnan(opnd2p1)) {
791da177e4SLinus Torvalds                         	/* trap if INVALIDTRAP enabled */
801da177e4SLinus Torvalds                         	if (Is_invalidtrap_enabled())
811da177e4SLinus Torvalds                             		return(INVALIDEXCEPTION);
821da177e4SLinus Torvalds                         	/* make NaN quiet */
831da177e4SLinus Torvalds                         	Set_invalidflag();
841da177e4SLinus Torvalds                         	Dbl_set_quiet(opnd2p1);
851da177e4SLinus Torvalds 				Dbl_copytoptr(opnd2p1,opnd2p2,dstptr);
861da177e4SLinus Torvalds                 		return(NOEXCEPTION);
871da177e4SLinus Torvalds 			}
881da177e4SLinus Torvalds                 	/*
891da177e4SLinus Torvalds                  	 * return quiet NaN
901da177e4SLinus Torvalds                  	 */
911da177e4SLinus Torvalds 			Dbl_copytoptr(opnd1p1,opnd1p2,dstptr);
921da177e4SLinus Torvalds                 	return(NOEXCEPTION);
931da177e4SLinus Torvalds 		}
941da177e4SLinus Torvalds 	}
951da177e4SLinus Torvalds 	/*
961da177e4SLinus Torvalds 	 * check second operand for NaN's or infinity
971da177e4SLinus Torvalds 	 */
981da177e4SLinus Torvalds 	if ((opnd2_exponent = Dbl_exponent(opnd2p1)) == DBL_INFINITY_EXPONENT) {
991da177e4SLinus Torvalds 		if (Dbl_iszero_mantissa(opnd2p1,opnd2p2)) {
1001da177e4SLinus Torvalds 			/*
1011da177e4SLinus Torvalds 			 * return first operand
1021da177e4SLinus Torvalds 			 */
1031da177e4SLinus Torvalds 			Dbl_copytoptr(opnd1p1,opnd1p2,dstptr);
1041da177e4SLinus Torvalds 			return(NOEXCEPTION);
1051da177e4SLinus Torvalds 		}
1061da177e4SLinus Torvalds                 /*
1071da177e4SLinus Torvalds                  * is NaN; signaling or quiet?
1081da177e4SLinus Torvalds                  */
1091da177e4SLinus Torvalds                 if (Dbl_isone_signaling(opnd2p1)) {
1101da177e4SLinus Torvalds                         /* trap if INVALIDTRAP enabled */
1111da177e4SLinus Torvalds                         if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
1121da177e4SLinus Torvalds                         /* make NaN quiet */
1131da177e4SLinus Torvalds                         Set_invalidflag();
1141da177e4SLinus Torvalds                         Dbl_set_quiet(opnd2p1);
1151da177e4SLinus Torvalds                 }
1161da177e4SLinus Torvalds                 /*
1171da177e4SLinus Torvalds                  * return quiet NaN
1181da177e4SLinus Torvalds                  */
1191da177e4SLinus Torvalds 		Dbl_copytoptr(opnd2p1,opnd2p2,dstptr);
1201da177e4SLinus Torvalds                 return(NOEXCEPTION);
1211da177e4SLinus Torvalds 	}
1221da177e4SLinus Torvalds 	/*
1231da177e4SLinus Torvalds 	 * check second operand for zero
1241da177e4SLinus Torvalds 	 */
1251da177e4SLinus Torvalds 	if (Dbl_iszero_exponentmantissa(opnd2p1,opnd2p2)) {
1261da177e4SLinus Torvalds 		/* invalid since second operand is zero */
1271da177e4SLinus Torvalds 		if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
1281da177e4SLinus Torvalds                 Set_invalidflag();
1291da177e4SLinus Torvalds                 Dbl_makequietnan(resultp1,resultp2);
1301da177e4SLinus Torvalds 		Dbl_copytoptr(resultp1,resultp2,dstptr);
1311da177e4SLinus Torvalds 		return(NOEXCEPTION);
1321da177e4SLinus Torvalds 	}
1331da177e4SLinus Torvalds 
1341da177e4SLinus Torvalds 	/*
1351da177e4SLinus Torvalds 	 * get sign of result
1361da177e4SLinus Torvalds 	 */
1371da177e4SLinus Torvalds 	resultp1 = opnd1p1;
1381da177e4SLinus Torvalds 
1391da177e4SLinus Torvalds 	/*
1401da177e4SLinus Torvalds 	 * check for denormalized operands
1411da177e4SLinus Torvalds 	 */
1421da177e4SLinus Torvalds 	if (opnd1_exponent == 0) {
1431da177e4SLinus Torvalds 		/* check for zero */
1441da177e4SLinus Torvalds 		if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) {
1451da177e4SLinus Torvalds 			Dbl_copytoptr(opnd1p1,opnd1p2,dstptr);
1461da177e4SLinus Torvalds 			return(NOEXCEPTION);
1471da177e4SLinus Torvalds 		}
1481da177e4SLinus Torvalds 		/* normalize, then continue */
1491da177e4SLinus Torvalds 		opnd1_exponent = 1;
1501da177e4SLinus Torvalds 		Dbl_normalize(opnd1p1,opnd1p2,opnd1_exponent);
1511da177e4SLinus Torvalds 	}
1521da177e4SLinus Torvalds 	else {
1531da177e4SLinus Torvalds 		Dbl_clear_signexponent_set_hidden(opnd1p1);
1541da177e4SLinus Torvalds 	}
1551da177e4SLinus Torvalds 	if (opnd2_exponent == 0) {
1561da177e4SLinus Torvalds 		/* normalize, then continue */
1571da177e4SLinus Torvalds 		opnd2_exponent = 1;
1581da177e4SLinus Torvalds 		Dbl_normalize(opnd2p1,opnd2p2,opnd2_exponent);
1591da177e4SLinus Torvalds 	}
1601da177e4SLinus Torvalds 	else {
1611da177e4SLinus Torvalds 		Dbl_clear_signexponent_set_hidden(opnd2p1);
1621da177e4SLinus Torvalds 	}
1631da177e4SLinus Torvalds 
1641da177e4SLinus Torvalds 	/* find result exponent and divide step loop count */
1651da177e4SLinus Torvalds 	dest_exponent = opnd2_exponent - 1;
1661da177e4SLinus Torvalds 	stepcount = opnd1_exponent - opnd2_exponent;
1671da177e4SLinus Torvalds 
1681da177e4SLinus Torvalds 	/*
1691da177e4SLinus Torvalds 	 * check for opnd1/opnd2 < 1
1701da177e4SLinus Torvalds 	 */
1711da177e4SLinus Torvalds 	if (stepcount < 0) {
1721da177e4SLinus Torvalds 		/*
1731da177e4SLinus Torvalds 		 * check for opnd1/opnd2 > 1/2
1741da177e4SLinus Torvalds 		 *
1751da177e4SLinus Torvalds 		 * In this case n will round to 1, so
1761da177e4SLinus Torvalds 		 *    r = opnd1 - opnd2
1771da177e4SLinus Torvalds 		 */
1781da177e4SLinus Torvalds 		if (stepcount == -1 &&
1791da177e4SLinus Torvalds 		    Dbl_isgreaterthan(opnd1p1,opnd1p2,opnd2p1,opnd2p2)) {
1801da177e4SLinus Torvalds 			/* set sign */
1811da177e4SLinus Torvalds 			Dbl_allp1(resultp1) = ~Dbl_allp1(resultp1);
1821da177e4SLinus Torvalds 			/* align opnd2 with opnd1 */
1831da177e4SLinus Torvalds 			Dbl_leftshiftby1(opnd2p1,opnd2p2);
1841da177e4SLinus Torvalds 			Dbl_subtract(opnd2p1,opnd2p2,opnd1p1,opnd1p2,
1851da177e4SLinus Torvalds 			 opnd2p1,opnd2p2);
1861da177e4SLinus Torvalds 			/* now normalize */
1871da177e4SLinus Torvalds                 	while (Dbl_iszero_hidden(opnd2p1)) {
1881da177e4SLinus Torvalds                         	Dbl_leftshiftby1(opnd2p1,opnd2p2);
1891da177e4SLinus Torvalds                         	dest_exponent--;
1901da177e4SLinus Torvalds 			}
1911da177e4SLinus Torvalds 			Dbl_set_exponentmantissa(resultp1,resultp2,opnd2p1,opnd2p2);
1921da177e4SLinus Torvalds 			goto testforunderflow;
1931da177e4SLinus Torvalds 		}
1941da177e4SLinus Torvalds 		/*
1951da177e4SLinus Torvalds 		 * opnd1/opnd2 <= 1/2
1961da177e4SLinus Torvalds 		 *
1971da177e4SLinus Torvalds 		 * In this case n will round to zero, so
1981da177e4SLinus Torvalds 		 *    r = opnd1
1991da177e4SLinus Torvalds 		 */
2001da177e4SLinus Torvalds 		Dbl_set_exponentmantissa(resultp1,resultp2,opnd1p1,opnd1p2);
2011da177e4SLinus Torvalds 		dest_exponent = opnd1_exponent;
2021da177e4SLinus Torvalds 		goto testforunderflow;
2031da177e4SLinus Torvalds 	}
2041da177e4SLinus Torvalds 
2051da177e4SLinus Torvalds 	/*
2061da177e4SLinus Torvalds 	 * Generate result
2071da177e4SLinus Torvalds 	 *
2081da177e4SLinus Torvalds 	 * Do iterative subtract until remainder is less than operand 2.
2091da177e4SLinus Torvalds 	 */
2101da177e4SLinus Torvalds 	while (stepcount-- > 0 && (Dbl_allp1(opnd1p1) || Dbl_allp2(opnd1p2))) {
2111da177e4SLinus Torvalds 		if (Dbl_isnotlessthan(opnd1p1,opnd1p2,opnd2p1,opnd2p2)) {
2121da177e4SLinus Torvalds 			Dbl_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2,opnd1p1,opnd1p2);
2131da177e4SLinus Torvalds 		}
2141da177e4SLinus Torvalds 		Dbl_leftshiftby1(opnd1p1,opnd1p2);
2151da177e4SLinus Torvalds 	}
2161da177e4SLinus Torvalds 	/*
2171da177e4SLinus Torvalds 	 * Do last subtract, then determine which way to round if remainder
2181da177e4SLinus Torvalds 	 * is exactly 1/2 of opnd2
2191da177e4SLinus Torvalds 	 */
2201da177e4SLinus Torvalds 	if (Dbl_isnotlessthan(opnd1p1,opnd1p2,opnd2p1,opnd2p2)) {
2211da177e4SLinus Torvalds 		Dbl_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2,opnd1p1,opnd1p2);
2221da177e4SLinus Torvalds 		roundup = TRUE;
2231da177e4SLinus Torvalds 	}
2241da177e4SLinus Torvalds 	if (stepcount > 0 || Dbl_iszero(opnd1p1,opnd1p2)) {
2251da177e4SLinus Torvalds 		/* division is exact, remainder is zero */
2261da177e4SLinus Torvalds 		Dbl_setzero_exponentmantissa(resultp1,resultp2);
2271da177e4SLinus Torvalds 		Dbl_copytoptr(resultp1,resultp2,dstptr);
2281da177e4SLinus Torvalds 		return(NOEXCEPTION);
2291da177e4SLinus Torvalds 	}
2301da177e4SLinus Torvalds 
2311da177e4SLinus Torvalds 	/*
2321da177e4SLinus Torvalds 	 * Check for cases where opnd1/opnd2 < n
2331da177e4SLinus Torvalds 	 *
2341da177e4SLinus Torvalds 	 * In this case the result's sign will be opposite that of
2351da177e4SLinus Torvalds 	 * opnd1.  The mantissa also needs some correction.
2361da177e4SLinus Torvalds 	 */
2371da177e4SLinus Torvalds 	Dbl_leftshiftby1(opnd1p1,opnd1p2);
2381da177e4SLinus Torvalds 	if (Dbl_isgreaterthan(opnd1p1,opnd1p2,opnd2p1,opnd2p2)) {
2391da177e4SLinus Torvalds 		Dbl_invert_sign(resultp1);
2401da177e4SLinus Torvalds 		Dbl_leftshiftby1(opnd2p1,opnd2p2);
2411da177e4SLinus Torvalds 		Dbl_subtract(opnd2p1,opnd2p2,opnd1p1,opnd1p2,opnd1p1,opnd1p2);
2421da177e4SLinus Torvalds 	}
2431da177e4SLinus Torvalds 	/* check for remainder being exactly 1/2 of opnd2 */
2441da177e4SLinus Torvalds 	else if (Dbl_isequal(opnd1p1,opnd1p2,opnd2p1,opnd2p2) && roundup) {
2451da177e4SLinus Torvalds 		Dbl_invert_sign(resultp1);
2461da177e4SLinus Torvalds 	}
2471da177e4SLinus Torvalds 
2481da177e4SLinus Torvalds 	/* normalize result's mantissa */
2491da177e4SLinus Torvalds         while (Dbl_iszero_hidden(opnd1p1)) {
2501da177e4SLinus Torvalds                 dest_exponent--;
2511da177e4SLinus Torvalds                 Dbl_leftshiftby1(opnd1p1,opnd1p2);
2521da177e4SLinus Torvalds         }
2531da177e4SLinus Torvalds 	Dbl_set_exponentmantissa(resultp1,resultp2,opnd1p1,opnd1p2);
2541da177e4SLinus Torvalds 
2551da177e4SLinus Torvalds         /*
2561da177e4SLinus Torvalds          * Test for underflow
2571da177e4SLinus Torvalds          */
2581da177e4SLinus Torvalds     testforunderflow:
2591da177e4SLinus Torvalds 	if (dest_exponent <= 0) {
2601da177e4SLinus Torvalds                 /* trap if UNDERFLOWTRAP enabled */
2611da177e4SLinus Torvalds                 if (Is_underflowtrap_enabled()) {
2621da177e4SLinus Torvalds                         /*
2631da177e4SLinus Torvalds                          * Adjust bias of result
2641da177e4SLinus Torvalds                          */
2651da177e4SLinus Torvalds                         Dbl_setwrapped_exponent(resultp1,dest_exponent,unfl);
2661da177e4SLinus Torvalds 			/* frem is always exact */
2671da177e4SLinus Torvalds 			Dbl_copytoptr(resultp1,resultp2,dstptr);
2681da177e4SLinus Torvalds 			return(UNDERFLOWEXCEPTION);
2691da177e4SLinus Torvalds                 }
2701da177e4SLinus Torvalds                 /*
2711da177e4SLinus Torvalds                  * denormalize result or set to signed zero
2721da177e4SLinus Torvalds                  */
2731da177e4SLinus Torvalds                 if (dest_exponent >= (1 - DBL_P)) {
2741da177e4SLinus Torvalds 			Dbl_rightshift_exponentmantissa(resultp1,resultp2,
2751da177e4SLinus Torvalds 			 1-dest_exponent);
2761da177e4SLinus Torvalds                 }
2771da177e4SLinus Torvalds                 else {
2781da177e4SLinus Torvalds 			Dbl_setzero_exponentmantissa(resultp1,resultp2);
2791da177e4SLinus Torvalds 		}
2801da177e4SLinus Torvalds 	}
2811da177e4SLinus Torvalds 	else Dbl_set_exponent(resultp1,dest_exponent);
2821da177e4SLinus Torvalds 	Dbl_copytoptr(resultp1,resultp2,dstptr);
2831da177e4SLinus Torvalds 	return(NOEXCEPTION);
2841da177e4SLinus Torvalds }
285