Linus Torvalds | 1da177e | 2005-04-16 15:20:36 -0700 | [diff] [blame] | 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 | } |