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/sfrem.c $Revision: 1.1 $ 26 * 27 * Purpose: 28 * Single Precision Floating-point Remainder 29 * 30 * External Interfaces: 31 * sgl_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 "sgl_float.h" 45 46 /* 47 * Single Precision Floating-point Remainder 48 */ 49 50 int 51 sgl_frem (sgl_floating_point * srcptr1, sgl_floating_point * srcptr2, 52 sgl_floating_point * dstptr, unsigned int *status) 53 { 54 register unsigned int opnd1, opnd2, result; 55 register int opnd1_exponent, opnd2_exponent, dest_exponent, stepcount; 56 register boolean roundup = FALSE; 57 58 opnd1 = *srcptr1; 59 opnd2 = *srcptr2; 60 /* 61 * check first operand for NaN's or infinity 62 */ 63 if ((opnd1_exponent = Sgl_exponent(opnd1)) == SGL_INFINITY_EXPONENT) { 64 if (Sgl_iszero_mantissa(opnd1)) { 65 if (Sgl_isnotnan(opnd2)) { 66 /* invalid since first operand is infinity */ 67 if (Is_invalidtrap_enabled()) 68 return(INVALIDEXCEPTION); 69 Set_invalidflag(); 70 Sgl_makequietnan(result); 71 *dstptr = result; 72 return(NOEXCEPTION); 73 } 74 } 75 else { 76 /* 77 * is NaN; signaling or quiet? 78 */ 79 if (Sgl_isone_signaling(opnd1)) { 80 /* trap if INVALIDTRAP enabled */ 81 if (Is_invalidtrap_enabled()) 82 return(INVALIDEXCEPTION); 83 /* make NaN quiet */ 84 Set_invalidflag(); 85 Sgl_set_quiet(opnd1); 86 } 87 /* 88 * is second operand a signaling NaN? 89 */ 90 else if (Sgl_is_signalingnan(opnd2)) { 91 /* trap if INVALIDTRAP enabled */ 92 if (Is_invalidtrap_enabled()) 93 return(INVALIDEXCEPTION); 94 /* make NaN quiet */ 95 Set_invalidflag(); 96 Sgl_set_quiet(opnd2); 97 *dstptr = opnd2; 98 return(NOEXCEPTION); 99 } 100 /* 101 * return quiet NaN 102 */ 103 *dstptr = opnd1; 104 return(NOEXCEPTION); 105 } 106 } 107 /* 108 * check second operand for NaN's or infinity 109 */ 110 if ((opnd2_exponent = Sgl_exponent(opnd2)) == SGL_INFINITY_EXPONENT) { 111 if (Sgl_iszero_mantissa(opnd2)) { 112 /* 113 * return first operand 114 */ 115 *dstptr = opnd1; 116 return(NOEXCEPTION); 117 } 118 /* 119 * is NaN; signaling or quiet? 120 */ 121 if (Sgl_isone_signaling(opnd2)) { 122 /* trap if INVALIDTRAP enabled */ 123 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 124 /* make NaN quiet */ 125 Set_invalidflag(); 126 Sgl_set_quiet(opnd2); 127 } 128 /* 129 * return quiet NaN 130 */ 131 *dstptr = opnd2; 132 return(NOEXCEPTION); 133 } 134 /* 135 * check second operand for zero 136 */ 137 if (Sgl_iszero_exponentmantissa(opnd2)) { 138 /* invalid since second operand is zero */ 139 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 140 Set_invalidflag(); 141 Sgl_makequietnan(result); 142 *dstptr = result; 143 return(NOEXCEPTION); 144 } 145 146 /* 147 * get sign of result 148 */ 149 result = opnd1; 150 151 /* 152 * check for denormalized operands 153 */ 154 if (opnd1_exponent == 0) { 155 /* check for zero */ 156 if (Sgl_iszero_mantissa(opnd1)) { 157 *dstptr = opnd1; 158 return(NOEXCEPTION); 159 } 160 /* normalize, then continue */ 161 opnd1_exponent = 1; 162 Sgl_normalize(opnd1,opnd1_exponent); 163 } 164 else { 165 Sgl_clear_signexponent_set_hidden(opnd1); 166 } 167 if (opnd2_exponent == 0) { 168 /* normalize, then continue */ 169 opnd2_exponent = 1; 170 Sgl_normalize(opnd2,opnd2_exponent); 171 } 172 else { 173 Sgl_clear_signexponent_set_hidden(opnd2); 174 } 175 176 /* find result exponent and divide step loop count */ 177 dest_exponent = opnd2_exponent - 1; 178 stepcount = opnd1_exponent - opnd2_exponent; 179 180 /* 181 * check for opnd1/opnd2 < 1 182 */ 183 if (stepcount < 0) { 184 /* 185 * check for opnd1/opnd2 > 1/2 186 * 187 * In this case n will round to 1, so 188 * r = opnd1 - opnd2 189 */ 190 if (stepcount == -1 && Sgl_isgreaterthan(opnd1,opnd2)) { 191 Sgl_all(result) = ~Sgl_all(result); /* set sign */ 192 /* align opnd2 with opnd1 */ 193 Sgl_leftshiftby1(opnd2); 194 Sgl_subtract(opnd2,opnd1,opnd2); 195 /* now normalize */ 196 while (Sgl_iszero_hidden(opnd2)) { 197 Sgl_leftshiftby1(opnd2); 198 dest_exponent--; 199 } 200 Sgl_set_exponentmantissa(result,opnd2); 201 goto testforunderflow; 202 } 203 /* 204 * opnd1/opnd2 <= 1/2 205 * 206 * In this case n will round to zero, so 207 * r = opnd1 208 */ 209 Sgl_set_exponentmantissa(result,opnd1); 210 dest_exponent = opnd1_exponent; 211 goto testforunderflow; 212 } 213 214 /* 215 * Generate result 216 * 217 * Do iterative subtract until remainder is less than operand 2. 218 */ 219 while (stepcount-- > 0 && Sgl_all(opnd1)) { 220 if (Sgl_isnotlessthan(opnd1,opnd2)) 221 Sgl_subtract(opnd1,opnd2,opnd1); 222 Sgl_leftshiftby1(opnd1); 223 } 224 /* 225 * Do last subtract, then determine which way to round if remainder 226 * is exactly 1/2 of opnd2 227 */ 228 if (Sgl_isnotlessthan(opnd1,opnd2)) { 229 Sgl_subtract(opnd1,opnd2,opnd1); 230 roundup = TRUE; 231 } 232 if (stepcount > 0 || Sgl_iszero(opnd1)) { 233 /* division is exact, remainder is zero */ 234 Sgl_setzero_exponentmantissa(result); 235 *dstptr = result; 236 return(NOEXCEPTION); 237 } 238 239 /* 240 * Check for cases where opnd1/opnd2 < n 241 * 242 * In this case the result's sign will be opposite that of 243 * opnd1. The mantissa also needs some correction. 244 */ 245 Sgl_leftshiftby1(opnd1); 246 if (Sgl_isgreaterthan(opnd1,opnd2)) { 247 Sgl_invert_sign(result); 248 Sgl_subtract((opnd2<<1),opnd1,opnd1); 249 } 250 /* check for remainder being exactly 1/2 of opnd2 */ 251 else if (Sgl_isequal(opnd1,opnd2) && roundup) { 252 Sgl_invert_sign(result); 253 } 254 255 /* normalize result's mantissa */ 256 while (Sgl_iszero_hidden(opnd1)) { 257 dest_exponent--; 258 Sgl_leftshiftby1(opnd1); 259 } 260 Sgl_set_exponentmantissa(result,opnd1); 261 262 /* 263 * Test for underflow 264 */ 265 testforunderflow: 266 if (dest_exponent <= 0) { 267 /* trap if UNDERFLOWTRAP enabled */ 268 if (Is_underflowtrap_enabled()) { 269 /* 270 * Adjust bias of result 271 */ 272 Sgl_setwrapped_exponent(result,dest_exponent,unfl); 273 *dstptr = result; 274 /* frem is always exact */ 275 return(UNDERFLOWEXCEPTION); 276 } 277 /* 278 * denormalize result or set to signed zero 279 */ 280 if (dest_exponent >= (1 - SGL_P)) { 281 Sgl_rightshift_exponentmantissa(result,1-dest_exponent); 282 } 283 else { 284 Sgl_setzero_exponentmantissa(result); 285 } 286 } 287 else Sgl_set_exponent(result,dest_exponent); 288 *dstptr = result; 289 return(NOEXCEPTION); 290 } 291