1 /* $OpenBSD: sfmpy.c,v 1.5 2002/05/07 22:19:30 mickey Exp $ */ 2 /* 3 (c) Copyright 1986 HEWLETT-PACKARD COMPANY 4 To anyone who acknowledges that this file is provided "AS IS" 5 without any express or implied warranty: 6 permission to use, copy, modify, and distribute this file 7 for any purpose is hereby granted without fee, provided that 8 the above copyright notice and this notice appears in all 9 copies, and that the name of Hewlett-Packard Company not be 10 used in advertising or publicity pertaining to distribution 11 of the software without specific, written prior permission. 12 Hewlett-Packard Company makes no representations about the 13 suitability of this software for any purpose. 14 */ 15 /* @(#)sfmpy.c: Revision: 2.9.88.1 Date: 93/12/07 15:07:07 */ 16 17 #include "float.h" 18 #include "sgl_float.h" 19 20 /* 21 * Single Precision Floating-point Multiply 22 */ 23 int 24 sgl_fmpy(srcptr1,srcptr2,dstptr,status) 25 sgl_floating_point *srcptr1, *srcptr2, *dstptr; 26 unsigned int *status; 27 { 28 register unsigned int opnd1, opnd2, opnd3, result; 29 register int dest_exponent, count; 30 register int inexact = FALSE, guardbit = FALSE, stickybit = FALSE; 31 int is_tiny; 32 33 opnd1 = *srcptr1; 34 opnd2 = *srcptr2; 35 /* 36 * set sign bit of result 37 */ 38 if (Sgl_sign(opnd1) ^ Sgl_sign(opnd2)) Sgl_setnegativezero(result); 39 else Sgl_setzero(result); 40 /* 41 * check first operand for NaN's or infinity 42 */ 43 if (Sgl_isinfinity_exponent(opnd1)) { 44 if (Sgl_iszero_mantissa(opnd1)) { 45 if (Sgl_isnotnan(opnd2)) { 46 if (Sgl_iszero_exponentmantissa(opnd2)) { 47 /* 48 * invalid since operands are infinity 49 * and zero 50 */ 51 if (Is_invalidtrap_enabled()) 52 return(INVALIDEXCEPTION); 53 Set_invalidflag(); 54 Sgl_makequietnan(result); 55 *dstptr = result; 56 return(NOEXCEPTION); 57 } 58 /* 59 * return infinity 60 */ 61 Sgl_setinfinity_exponentmantissa(result); 62 *dstptr = result; 63 return(NOEXCEPTION); 64 } 65 } 66 else { 67 /* 68 * is NaN; signaling or quiet? 69 */ 70 if (Sgl_isone_signaling(opnd1)) { 71 /* trap if INVALIDTRAP enabled */ 72 if (Is_invalidtrap_enabled()) 73 return(INVALIDEXCEPTION); 74 /* make NaN quiet */ 75 Set_invalidflag(); 76 Sgl_set_quiet(opnd1); 77 } 78 /* 79 * is second operand a signaling NaN? 80 */ 81 else if (Sgl_is_signalingnan(opnd2)) { 82 /* trap if INVALIDTRAP enabled */ 83 if (Is_invalidtrap_enabled()) 84 return(INVALIDEXCEPTION); 85 /* make NaN quiet */ 86 Set_invalidflag(); 87 Sgl_set_quiet(opnd2); 88 *dstptr = opnd2; 89 return(NOEXCEPTION); 90 } 91 /* 92 * return quiet NaN 93 */ 94 *dstptr = opnd1; 95 return(NOEXCEPTION); 96 } 97 } 98 /* 99 * check second operand for NaN's or infinity 100 */ 101 if (Sgl_isinfinity_exponent(opnd2)) { 102 if (Sgl_iszero_mantissa(opnd2)) { 103 if (Sgl_iszero_exponentmantissa(opnd1)) { 104 /* invalid since operands are zero & infinity */ 105 if (Is_invalidtrap_enabled()) 106 return(INVALIDEXCEPTION); 107 Set_invalidflag(); 108 Sgl_makequietnan(opnd2); 109 *dstptr = opnd2; 110 return(NOEXCEPTION); 111 } 112 /* 113 * return infinity 114 */ 115 Sgl_setinfinity_exponentmantissa(result); 116 *dstptr = result; 117 return(NOEXCEPTION); 118 } 119 /* 120 * is NaN; signaling or quiet? 121 */ 122 if (Sgl_isone_signaling(opnd2)) { 123 /* trap if INVALIDTRAP enabled */ 124 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 125 126 /* make NaN quiet */ 127 Set_invalidflag(); 128 Sgl_set_quiet(opnd2); 129 } 130 /* 131 * return quiet NaN 132 */ 133 *dstptr = opnd2; 134 return(NOEXCEPTION); 135 } 136 /* 137 * Generate exponent 138 */ 139 dest_exponent = Sgl_exponent(opnd1) + Sgl_exponent(opnd2) - SGL_BIAS; 140 141 /* 142 * Generate mantissa 143 */ 144 if (Sgl_isnotzero_exponent(opnd1)) { 145 /* set hidden bit */ 146 Sgl_clear_signexponent_set_hidden(opnd1); 147 } 148 else { 149 /* check for zero */ 150 if (Sgl_iszero_mantissa(opnd1)) { 151 Sgl_setzero_exponentmantissa(result); 152 *dstptr = result; 153 return(NOEXCEPTION); 154 } 155 /* is denormalized, adjust exponent */ 156 Sgl_clear_signexponent(opnd1); 157 Sgl_leftshiftby1(opnd1); 158 Sgl_normalize(opnd1,dest_exponent); 159 } 160 /* opnd2 needs to have hidden bit set with msb in hidden bit */ 161 if (Sgl_isnotzero_exponent(opnd2)) { 162 Sgl_clear_signexponent_set_hidden(opnd2); 163 } 164 else { 165 /* check for zero */ 166 if (Sgl_iszero_mantissa(opnd2)) { 167 Sgl_setzero_exponentmantissa(result); 168 *dstptr = result; 169 return(NOEXCEPTION); 170 } 171 /* is denormalized; want to normalize */ 172 Sgl_clear_signexponent(opnd2); 173 Sgl_leftshiftby1(opnd2); 174 Sgl_normalize(opnd2,dest_exponent); 175 } 176 177 /* Multiply two source mantissas together */ 178 179 Sgl_leftshiftby4(opnd2); /* make room for guard bits */ 180 Sgl_setzero(opnd3); 181 /* 182 * Four bits at a time are inspected in each loop, and a 183 * simple shift and add multiply algorithm is used. 184 */ 185 for (count=1;count<SGL_P;count+=4) { 186 stickybit |= Slow4(opnd3); 187 Sgl_rightshiftby4(opnd3); 188 if (Sbit28(opnd1)) Sall(opnd3) += (Sall(opnd2) << 3); 189 if (Sbit29(opnd1)) Sall(opnd3) += (Sall(opnd2) << 2); 190 if (Sbit30(opnd1)) Sall(opnd3) += (Sall(opnd2) << 1); 191 if (Sbit31(opnd1)) Sall(opnd3) += Sall(opnd2); 192 Sgl_rightshiftby4(opnd1); 193 } 194 /* make sure result is left-justified */ 195 if (Sgl_iszero_sign(opnd3)) { 196 Sgl_leftshiftby1(opnd3); 197 } 198 else { 199 /* result mantissa >= 2. */ 200 dest_exponent++; 201 } 202 /* check for denormalized result */ 203 while (Sgl_iszero_sign(opnd3)) { 204 Sgl_leftshiftby1(opnd3); 205 dest_exponent--; 206 } 207 /* 208 * check for guard, sticky and inexact bits 209 */ 210 stickybit |= Sgl_all(opnd3) << (SGL_BITLENGTH - SGL_EXP_LENGTH + 1); 211 guardbit = Sbit24(opnd3); 212 inexact = guardbit | stickybit; 213 214 /* re-align mantissa */ 215 Sgl_rightshiftby8(opnd3); 216 217 /* 218 * round result 219 */ 220 if (inexact && (dest_exponent>0 || Is_underflowtrap_enabled())) { 221 Sgl_clear_signexponent(opnd3); 222 switch (Rounding_mode()) { 223 case ROUNDPLUS: 224 if (Sgl_iszero_sign(result)) 225 Sgl_increment(opnd3); 226 break; 227 case ROUNDMINUS: 228 if (Sgl_isone_sign(result)) 229 Sgl_increment(opnd3); 230 break; 231 case ROUNDNEAREST: 232 if (guardbit && 233 (stickybit || Sgl_isone_lowmantissa(opnd3))) 234 Sgl_increment(opnd3); 235 break; 236 } 237 if (Sgl_isone_hidden(opnd3)) dest_exponent++; 238 } 239 Sgl_set_mantissa(result,opnd3); 240 241 /* 242 * Test for overflow 243 */ 244 if (dest_exponent >= SGL_INFINITY_EXPONENT) { 245 /* trap if OVERFLOWTRAP enabled */ 246 if (Is_overflowtrap_enabled()) { 247 /* 248 * Adjust bias of result 249 */ 250 Sgl_setwrapped_exponent(result,dest_exponent,ovfl); 251 *dstptr = result; 252 if (inexact) { 253 if (Is_inexacttrap_enabled()) 254 return(OVERFLOWEXCEPTION | INEXACTEXCEPTION); 255 else Set_inexactflag(); 256 } 257 return(OVERFLOWEXCEPTION); 258 } 259 inexact = TRUE; 260 Set_overflowflag(); 261 /* set result to infinity or largest number */ 262 Sgl_setoverflow(result); 263 } 264 /* 265 * Test for underflow 266 */ 267 else if (dest_exponent <= 0) { 268 /* trap if UNDERFLOWTRAP enabled */ 269 if (Is_underflowtrap_enabled()) { 270 /* 271 * Adjust bias of result 272 */ 273 Sgl_setwrapped_exponent(result,dest_exponent,unfl); 274 *dstptr = result; 275 if (inexact) { 276 if (Is_inexacttrap_enabled()) 277 return(UNDERFLOWEXCEPTION | INEXACTEXCEPTION); 278 else Set_inexactflag(); 279 } 280 return(UNDERFLOWEXCEPTION); 281 } 282 283 /* Determine if should set underflow flag */ 284 is_tiny = TRUE; 285 if (dest_exponent == 0 && inexact) { 286 switch (Rounding_mode()) { 287 case ROUNDPLUS: 288 if (Sgl_iszero_sign(result)) { 289 Sgl_increment(opnd3); 290 if (Sgl_isone_hiddenoverflow(opnd3)) 291 is_tiny = FALSE; 292 Sgl_decrement(opnd3); 293 } 294 break; 295 case ROUNDMINUS: 296 if (Sgl_isone_sign(result)) { 297 Sgl_increment(opnd3); 298 if (Sgl_isone_hiddenoverflow(opnd3)) 299 is_tiny = FALSE; 300 Sgl_decrement(opnd3); 301 } 302 break; 303 case ROUNDNEAREST: 304 if (guardbit && (stickybit || 305 Sgl_isone_lowmantissa(opnd3))) { 306 Sgl_increment(opnd3); 307 if (Sgl_isone_hiddenoverflow(opnd3)) 308 is_tiny = FALSE; 309 Sgl_decrement(opnd3); 310 } 311 break; 312 } 313 } 314 315 /* 316 * denormalize result or set to signed zero 317 */ 318 stickybit = inexact; 319 Sgl_denormalize(opnd3,dest_exponent,guardbit,stickybit,inexact); 320 321 /* return zero or smallest number */ 322 if (inexact) { 323 switch (Rounding_mode()) { 324 case ROUNDPLUS: 325 if (Sgl_iszero_sign(result)) 326 Sgl_increment(opnd3); 327 break; 328 case ROUNDMINUS: 329 if (Sgl_isone_sign(result)) 330 Sgl_increment(opnd3); 331 break; 332 case ROUNDNEAREST: 333 if (guardbit && (stickybit || 334 Sgl_isone_lowmantissa(opnd3))) 335 Sgl_increment(opnd3); 336 break; 337 } 338 if (is_tiny) Set_underflowflag(); 339 } 340 Sgl_set_exponentmantissa(result,opnd3); 341 } 342 else Sgl_set_exponent(result,dest_exponent); 343 *dstptr = result; 344 345 /* check for inexact */ 346 if (inexact) { 347 if (Is_inexacttrap_enabled()) return(INEXACTEXCEPTION); 348 else Set_inexactflag(); 349 } 350 return(NOEXCEPTION); 351 } 352