1 // SPDX-License-Identifier: GPL-2.0-or-later 2 /* 3 * Linux/PA-RISC Project (http://www.parisc-linux.org/) 4 * 5 * Floating-point emulation code 6 * Copyright (C) 2001 Hewlett-Packard (Paul Bame) <bame@debian.org> 7 */ 8 /* 9 * BEGIN_DESC 10 * 11 * File: 12 * @(#) pa/spmath/dfrem.c $Revision: 1.1 $ 13 * 14 * Purpose: 15 * Double Precision Floating-point Remainder 16 * 17 * External Interfaces: 18 * dbl_frem(srcptr1,srcptr2,dstptr,status) 19 * 20 * Internal Interfaces: 21 * 22 * Theory: 23 * <<please update with a overview of the operation of this file>> 24 * 25 * END_DESC 26 */ 27 28 29 30 #include "float.h" 31 #include "dbl_float.h" 32 33 /* 34 * Double Precision Floating-point Remainder 35 */ 36 37 int 38 dbl_frem (dbl_floating_point * srcptr1, dbl_floating_point * srcptr2, 39 dbl_floating_point * dstptr, unsigned int *status) 40 { 41 register unsigned int opnd1p1, opnd1p2, opnd2p1, opnd2p2; 42 register unsigned int resultp1, resultp2; 43 register int opnd1_exponent, opnd2_exponent, dest_exponent, stepcount; 44 register boolean roundup = FALSE; 45 46 Dbl_copyfromptr(srcptr1,opnd1p1,opnd1p2); 47 Dbl_copyfromptr(srcptr2,opnd2p1,opnd2p2); 48 /* 49 * check first operand for NaN's or infinity 50 */ 51 if ((opnd1_exponent = Dbl_exponent(opnd1p1)) == DBL_INFINITY_EXPONENT) { 52 if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) { 53 if (Dbl_isnotnan(opnd2p1,opnd2p2)) { 54 /* invalid since first operand is infinity */ 55 if (Is_invalidtrap_enabled()) 56 return(INVALIDEXCEPTION); 57 Set_invalidflag(); 58 Dbl_makequietnan(resultp1,resultp2); 59 Dbl_copytoptr(resultp1,resultp2,dstptr); 60 return(NOEXCEPTION); 61 } 62 } 63 else { 64 /* 65 * is NaN; signaling or quiet? 66 */ 67 if (Dbl_isone_signaling(opnd1p1)) { 68 /* trap if INVALIDTRAP enabled */ 69 if (Is_invalidtrap_enabled()) 70 return(INVALIDEXCEPTION); 71 /* make NaN quiet */ 72 Set_invalidflag(); 73 Dbl_set_quiet(opnd1p1); 74 } 75 /* 76 * is second operand a signaling NaN? 77 */ 78 else if (Dbl_is_signalingnan(opnd2p1)) { 79 /* trap if INVALIDTRAP enabled */ 80 if (Is_invalidtrap_enabled()) 81 return(INVALIDEXCEPTION); 82 /* make NaN quiet */ 83 Set_invalidflag(); 84 Dbl_set_quiet(opnd2p1); 85 Dbl_copytoptr(opnd2p1,opnd2p2,dstptr); 86 return(NOEXCEPTION); 87 } 88 /* 89 * return quiet NaN 90 */ 91 Dbl_copytoptr(opnd1p1,opnd1p2,dstptr); 92 return(NOEXCEPTION); 93 } 94 } 95 /* 96 * check second operand for NaN's or infinity 97 */ 98 if ((opnd2_exponent = Dbl_exponent(opnd2p1)) == DBL_INFINITY_EXPONENT) { 99 if (Dbl_iszero_mantissa(opnd2p1,opnd2p2)) { 100 /* 101 * return first operand 102 */ 103 Dbl_copytoptr(opnd1p1,opnd1p2,dstptr); 104 return(NOEXCEPTION); 105 } 106 /* 107 * is NaN; signaling or quiet? 108 */ 109 if (Dbl_isone_signaling(opnd2p1)) { 110 /* trap if INVALIDTRAP enabled */ 111 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 112 /* make NaN quiet */ 113 Set_invalidflag(); 114 Dbl_set_quiet(opnd2p1); 115 } 116 /* 117 * return quiet NaN 118 */ 119 Dbl_copytoptr(opnd2p1,opnd2p2,dstptr); 120 return(NOEXCEPTION); 121 } 122 /* 123 * check second operand for zero 124 */ 125 if (Dbl_iszero_exponentmantissa(opnd2p1,opnd2p2)) { 126 /* invalid since second operand is zero */ 127 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 128 Set_invalidflag(); 129 Dbl_makequietnan(resultp1,resultp2); 130 Dbl_copytoptr(resultp1,resultp2,dstptr); 131 return(NOEXCEPTION); 132 } 133 134 /* 135 * get sign of result 136 */ 137 resultp1 = opnd1p1; 138 139 /* 140 * check for denormalized operands 141 */ 142 if (opnd1_exponent == 0) { 143 /* check for zero */ 144 if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) { 145 Dbl_copytoptr(opnd1p1,opnd1p2,dstptr); 146 return(NOEXCEPTION); 147 } 148 /* normalize, then continue */ 149 opnd1_exponent = 1; 150 Dbl_normalize(opnd1p1,opnd1p2,opnd1_exponent); 151 } 152 else { 153 Dbl_clear_signexponent_set_hidden(opnd1p1); 154 } 155 if (opnd2_exponent == 0) { 156 /* normalize, then continue */ 157 opnd2_exponent = 1; 158 Dbl_normalize(opnd2p1,opnd2p2,opnd2_exponent); 159 } 160 else { 161 Dbl_clear_signexponent_set_hidden(opnd2p1); 162 } 163 164 /* find result exponent and divide step loop count */ 165 dest_exponent = opnd2_exponent - 1; 166 stepcount = opnd1_exponent - opnd2_exponent; 167 168 /* 169 * check for opnd1/opnd2 < 1 170 */ 171 if (stepcount < 0) { 172 /* 173 * check for opnd1/opnd2 > 1/2 174 * 175 * In this case n will round to 1, so 176 * r = opnd1 - opnd2 177 */ 178 if (stepcount == -1 && 179 Dbl_isgreaterthan(opnd1p1,opnd1p2,opnd2p1,opnd2p2)) { 180 /* set sign */ 181 Dbl_allp1(resultp1) = ~Dbl_allp1(resultp1); 182 /* align opnd2 with opnd1 */ 183 Dbl_leftshiftby1(opnd2p1,opnd2p2); 184 Dbl_subtract(opnd2p1,opnd2p2,opnd1p1,opnd1p2, 185 opnd2p1,opnd2p2); 186 /* now normalize */ 187 while (Dbl_iszero_hidden(opnd2p1)) { 188 Dbl_leftshiftby1(opnd2p1,opnd2p2); 189 dest_exponent--; 190 } 191 Dbl_set_exponentmantissa(resultp1,resultp2,opnd2p1,opnd2p2); 192 goto testforunderflow; 193 } 194 /* 195 * opnd1/opnd2 <= 1/2 196 * 197 * In this case n will round to zero, so 198 * r = opnd1 199 */ 200 Dbl_set_exponentmantissa(resultp1,resultp2,opnd1p1,opnd1p2); 201 dest_exponent = opnd1_exponent; 202 goto testforunderflow; 203 } 204 205 /* 206 * Generate result 207 * 208 * Do iterative subtract until remainder is less than operand 2. 209 */ 210 while (stepcount-- > 0 && (Dbl_allp1(opnd1p1) || Dbl_allp2(opnd1p2))) { 211 if (Dbl_isnotlessthan(opnd1p1,opnd1p2,opnd2p1,opnd2p2)) { 212 Dbl_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2,opnd1p1,opnd1p2); 213 } 214 Dbl_leftshiftby1(opnd1p1,opnd1p2); 215 } 216 /* 217 * Do last subtract, then determine which way to round if remainder 218 * is exactly 1/2 of opnd2 219 */ 220 if (Dbl_isnotlessthan(opnd1p1,opnd1p2,opnd2p1,opnd2p2)) { 221 Dbl_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2,opnd1p1,opnd1p2); 222 roundup = TRUE; 223 } 224 if (stepcount > 0 || Dbl_iszero(opnd1p1,opnd1p2)) { 225 /* division is exact, remainder is zero */ 226 Dbl_setzero_exponentmantissa(resultp1,resultp2); 227 Dbl_copytoptr(resultp1,resultp2,dstptr); 228 return(NOEXCEPTION); 229 } 230 231 /* 232 * Check for cases where opnd1/opnd2 < n 233 * 234 * In this case the result's sign will be opposite that of 235 * opnd1. The mantissa also needs some correction. 236 */ 237 Dbl_leftshiftby1(opnd1p1,opnd1p2); 238 if (Dbl_isgreaterthan(opnd1p1,opnd1p2,opnd2p1,opnd2p2)) { 239 Dbl_invert_sign(resultp1); 240 Dbl_leftshiftby1(opnd2p1,opnd2p2); 241 Dbl_subtract(opnd2p1,opnd2p2,opnd1p1,opnd1p2,opnd1p1,opnd1p2); 242 } 243 /* check for remainder being exactly 1/2 of opnd2 */ 244 else if (Dbl_isequal(opnd1p1,opnd1p2,opnd2p1,opnd2p2) && roundup) { 245 Dbl_invert_sign(resultp1); 246 } 247 248 /* normalize result's mantissa */ 249 while (Dbl_iszero_hidden(opnd1p1)) { 250 dest_exponent--; 251 Dbl_leftshiftby1(opnd1p1,opnd1p2); 252 } 253 Dbl_set_exponentmantissa(resultp1,resultp2,opnd1p1,opnd1p2); 254 255 /* 256 * Test for underflow 257 */ 258 testforunderflow: 259 if (dest_exponent <= 0) { 260 /* trap if UNDERFLOWTRAP enabled */ 261 if (Is_underflowtrap_enabled()) { 262 /* 263 * Adjust bias of result 264 */ 265 Dbl_setwrapped_exponent(resultp1,dest_exponent,unfl); 266 /* frem is always exact */ 267 Dbl_copytoptr(resultp1,resultp2,dstptr); 268 return(UNDERFLOWEXCEPTION); 269 } 270 /* 271 * denormalize result or set to signed zero 272 */ 273 if (dest_exponent >= (1 - DBL_P)) { 274 Dbl_rightshift_exponentmantissa(resultp1,resultp2, 275 1-dest_exponent); 276 } 277 else { 278 Dbl_setzero_exponentmantissa(resultp1,resultp2); 279 } 280 } 281 else Dbl_set_exponent(resultp1,dest_exponent); 282 Dbl_copytoptr(resultp1,resultp2,dstptr); 283 return(NOEXCEPTION); 284 } 285