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/sfmpy.c $Revision: 1.1 $ 26 * 27 * Purpose: 28 * Single Precision Floating-point Multiply 29 * 30 * External Interfaces: 31 * sgl_fmpy(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 #include "float.h" 43 #include "sgl_float.h" 44 45 /* 46 * Single Precision Floating-point Multiply 47 */ 48 49 int 50 sgl_fmpy( 51 sgl_floating_point *srcptr1, 52 sgl_floating_point *srcptr2, 53 sgl_floating_point *dstptr, 54 unsigned int *status) 55 { 56 register unsigned int opnd1, opnd2, opnd3, result; 57 register int dest_exponent, count; 58 register boolean inexact = FALSE, guardbit = FALSE, stickybit = FALSE; 59 boolean is_tiny; 60 61 opnd1 = *srcptr1; 62 opnd2 = *srcptr2; 63 /* 64 * set sign bit of result 65 */ 66 if (Sgl_sign(opnd1) ^ Sgl_sign(opnd2)) Sgl_setnegativezero(result); 67 else Sgl_setzero(result); 68 /* 69 * check first operand for NaN's or infinity 70 */ 71 if (Sgl_isinfinity_exponent(opnd1)) { 72 if (Sgl_iszero_mantissa(opnd1)) { 73 if (Sgl_isnotnan(opnd2)) { 74 if (Sgl_iszero_exponentmantissa(opnd2)) { 75 /* 76 * invalid since operands are infinity 77 * and zero 78 */ 79 if (Is_invalidtrap_enabled()) 80 return(INVALIDEXCEPTION); 81 Set_invalidflag(); 82 Sgl_makequietnan(result); 83 *dstptr = result; 84 return(NOEXCEPTION); 85 } 86 /* 87 * return infinity 88 */ 89 Sgl_setinfinity_exponentmantissa(result); 90 *dstptr = result; 91 return(NOEXCEPTION); 92 } 93 } 94 else { 95 /* 96 * is NaN; signaling or quiet? 97 */ 98 if (Sgl_isone_signaling(opnd1)) { 99 /* trap if INVALIDTRAP enabled */ 100 if (Is_invalidtrap_enabled()) 101 return(INVALIDEXCEPTION); 102 /* make NaN quiet */ 103 Set_invalidflag(); 104 Sgl_set_quiet(opnd1); 105 } 106 /* 107 * is second operand a signaling NaN? 108 */ 109 else if (Sgl_is_signalingnan(opnd2)) { 110 /* trap if INVALIDTRAP enabled */ 111 if (Is_invalidtrap_enabled()) 112 return(INVALIDEXCEPTION); 113 /* make NaN quiet */ 114 Set_invalidflag(); 115 Sgl_set_quiet(opnd2); 116 *dstptr = opnd2; 117 return(NOEXCEPTION); 118 } 119 /* 120 * return quiet NaN 121 */ 122 *dstptr = opnd1; 123 return(NOEXCEPTION); 124 } 125 } 126 /* 127 * check second operand for NaN's or infinity 128 */ 129 if (Sgl_isinfinity_exponent(opnd2)) { 130 if (Sgl_iszero_mantissa(opnd2)) { 131 if (Sgl_iszero_exponentmantissa(opnd1)) { 132 /* invalid since operands are zero & infinity */ 133 if (Is_invalidtrap_enabled()) 134 return(INVALIDEXCEPTION); 135 Set_invalidflag(); 136 Sgl_makequietnan(opnd2); 137 *dstptr = opnd2; 138 return(NOEXCEPTION); 139 } 140 /* 141 * return infinity 142 */ 143 Sgl_setinfinity_exponentmantissa(result); 144 *dstptr = result; 145 return(NOEXCEPTION); 146 } 147 /* 148 * is NaN; signaling or quiet? 149 */ 150 if (Sgl_isone_signaling(opnd2)) { 151 /* trap if INVALIDTRAP enabled */ 152 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 153 154 /* make NaN quiet */ 155 Set_invalidflag(); 156 Sgl_set_quiet(opnd2); 157 } 158 /* 159 * return quiet NaN 160 */ 161 *dstptr = opnd2; 162 return(NOEXCEPTION); 163 } 164 /* 165 * Generate exponent 166 */ 167 dest_exponent = Sgl_exponent(opnd1) + Sgl_exponent(opnd2) - SGL_BIAS; 168 169 /* 170 * Generate mantissa 171 */ 172 if (Sgl_isnotzero_exponent(opnd1)) { 173 /* set hidden bit */ 174 Sgl_clear_signexponent_set_hidden(opnd1); 175 } 176 else { 177 /* check for zero */ 178 if (Sgl_iszero_mantissa(opnd1)) { 179 Sgl_setzero_exponentmantissa(result); 180 *dstptr = result; 181 return(NOEXCEPTION); 182 } 183 /* is denormalized, adjust exponent */ 184 Sgl_clear_signexponent(opnd1); 185 Sgl_leftshiftby1(opnd1); 186 Sgl_normalize(opnd1,dest_exponent); 187 } 188 /* opnd2 needs to have hidden bit set with msb in hidden bit */ 189 if (Sgl_isnotzero_exponent(opnd2)) { 190 Sgl_clear_signexponent_set_hidden(opnd2); 191 } 192 else { 193 /* check for zero */ 194 if (Sgl_iszero_mantissa(opnd2)) { 195 Sgl_setzero_exponentmantissa(result); 196 *dstptr = result; 197 return(NOEXCEPTION); 198 } 199 /* is denormalized; want to normalize */ 200 Sgl_clear_signexponent(opnd2); 201 Sgl_leftshiftby1(opnd2); 202 Sgl_normalize(opnd2,dest_exponent); 203 } 204 205 /* Multiply two source mantissas together */ 206 207 Sgl_leftshiftby4(opnd2); /* make room for guard bits */ 208 Sgl_setzero(opnd3); 209 /* 210 * Four bits at a time are inspected in each loop, and a 211 * simple shift and add multiply algorithm is used. 212 */ 213 for (count=1;count<SGL_P;count+=4) { 214 stickybit |= Slow4(opnd3); 215 Sgl_rightshiftby4(opnd3); 216 if (Sbit28(opnd1)) Sall(opnd3) += (Sall(opnd2) << 3); 217 if (Sbit29(opnd1)) Sall(opnd3) += (Sall(opnd2) << 2); 218 if (Sbit30(opnd1)) Sall(opnd3) += (Sall(opnd2) << 1); 219 if (Sbit31(opnd1)) Sall(opnd3) += Sall(opnd2); 220 Sgl_rightshiftby4(opnd1); 221 } 222 /* make sure result is left-justified */ 223 if (Sgl_iszero_sign(opnd3)) { 224 Sgl_leftshiftby1(opnd3); 225 } 226 else { 227 /* result mantissa >= 2. */ 228 dest_exponent++; 229 } 230 /* check for denormalized result */ 231 while (Sgl_iszero_sign(opnd3)) { 232 Sgl_leftshiftby1(opnd3); 233 dest_exponent--; 234 } 235 /* 236 * check for guard, sticky and inexact bits 237 */ 238 stickybit |= Sgl_all(opnd3) << (SGL_BITLENGTH - SGL_EXP_LENGTH + 1); 239 guardbit = Sbit24(opnd3); 240 inexact = guardbit | stickybit; 241 242 /* re-align mantissa */ 243 Sgl_rightshiftby8(opnd3); 244 245 /* 246 * round result 247 */ 248 if (inexact && (dest_exponent>0 || Is_underflowtrap_enabled())) { 249 Sgl_clear_signexponent(opnd3); 250 switch (Rounding_mode()) { 251 case ROUNDPLUS: 252 if (Sgl_iszero_sign(result)) 253 Sgl_increment(opnd3); 254 break; 255 case ROUNDMINUS: 256 if (Sgl_isone_sign(result)) 257 Sgl_increment(opnd3); 258 break; 259 case ROUNDNEAREST: 260 if (guardbit) { 261 if (stickybit || Sgl_isone_lowmantissa(opnd3)) 262 Sgl_increment(opnd3); 263 } 264 } 265 if (Sgl_isone_hidden(opnd3)) dest_exponent++; 266 } 267 Sgl_set_mantissa(result,opnd3); 268 269 /* 270 * Test for overflow 271 */ 272 if (dest_exponent >= SGL_INFINITY_EXPONENT) { 273 /* trap if OVERFLOWTRAP enabled */ 274 if (Is_overflowtrap_enabled()) { 275 /* 276 * Adjust bias of result 277 */ 278 Sgl_setwrapped_exponent(result,dest_exponent,ovfl); 279 *dstptr = result; 280 if (inexact) 281 if (Is_inexacttrap_enabled()) 282 return(OVERFLOWEXCEPTION | INEXACTEXCEPTION); 283 else Set_inexactflag(); 284 return(OVERFLOWEXCEPTION); 285 } 286 inexact = TRUE; 287 Set_overflowflag(); 288 /* set result to infinity or largest number */ 289 Sgl_setoverflow(result); 290 } 291 /* 292 * Test for underflow 293 */ 294 else if (dest_exponent <= 0) { 295 /* trap if UNDERFLOWTRAP enabled */ 296 if (Is_underflowtrap_enabled()) { 297 /* 298 * Adjust bias of result 299 */ 300 Sgl_setwrapped_exponent(result,dest_exponent,unfl); 301 *dstptr = result; 302 if (inexact) 303 if (Is_inexacttrap_enabled()) 304 return(UNDERFLOWEXCEPTION | INEXACTEXCEPTION); 305 else Set_inexactflag(); 306 return(UNDERFLOWEXCEPTION); 307 } 308 309 /* Determine if should set underflow flag */ 310 is_tiny = TRUE; 311 if (dest_exponent == 0 && inexact) { 312 switch (Rounding_mode()) { 313 case ROUNDPLUS: 314 if (Sgl_iszero_sign(result)) { 315 Sgl_increment(opnd3); 316 if (Sgl_isone_hiddenoverflow(opnd3)) 317 is_tiny = FALSE; 318 Sgl_decrement(opnd3); 319 } 320 break; 321 case ROUNDMINUS: 322 if (Sgl_isone_sign(result)) { 323 Sgl_increment(opnd3); 324 if (Sgl_isone_hiddenoverflow(opnd3)) 325 is_tiny = FALSE; 326 Sgl_decrement(opnd3); 327 } 328 break; 329 case ROUNDNEAREST: 330 if (guardbit && (stickybit || 331 Sgl_isone_lowmantissa(opnd3))) { 332 Sgl_increment(opnd3); 333 if (Sgl_isone_hiddenoverflow(opnd3)) 334 is_tiny = FALSE; 335 Sgl_decrement(opnd3); 336 } 337 break; 338 } 339 } 340 341 /* 342 * denormalize result or set to signed zero 343 */ 344 stickybit = inexact; 345 Sgl_denormalize(opnd3,dest_exponent,guardbit,stickybit,inexact); 346 347 /* return zero or smallest number */ 348 if (inexact) { 349 switch (Rounding_mode()) { 350 case ROUNDPLUS: 351 if (Sgl_iszero_sign(result)) { 352 Sgl_increment(opnd3); 353 } 354 break; 355 case ROUNDMINUS: 356 if (Sgl_isone_sign(result)) { 357 Sgl_increment(opnd3); 358 } 359 break; 360 case ROUNDNEAREST: 361 if (guardbit && (stickybit || 362 Sgl_isone_lowmantissa(opnd3))) { 363 Sgl_increment(opnd3); 364 } 365 break; 366 } 367 if (is_tiny) Set_underflowflag(); 368 } 369 Sgl_set_exponentmantissa(result,opnd3); 370 } 371 else Sgl_set_exponent(result,dest_exponent); 372 *dstptr = result; 373 374 /* check for inexact */ 375 if (inexact) { 376 if (Is_inexacttrap_enabled()) return(INEXACTEXCEPTION); 377 else Set_inexactflag(); 378 } 379 return(NOEXCEPTION); 380 } 381