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