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