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