1 // SPDX-License-Identifier: GPL-2.0-or-later 2 /* 3 * Linux/PA-RISC Project (http://www.parisc-linux.org/) 4 * 5 * Floating-point emulation code 6 * Copyright (C) 2001 Hewlett-Packard (Paul Bame) <bame@debian.org> 7 */ 8 /* 9 * BEGIN_DESC 10 * 11 * File: 12 * @(#) pa/spmath/dfmpy.c $Revision: 1.1 $ 13 * 14 * Purpose: 15 * Double Precision Floating-point Multiply 16 * 17 * External Interfaces: 18 * dbl_fmpy(srcptr1,srcptr2,dstptr,status) 19 * 20 * Internal Interfaces: 21 * 22 * Theory: 23 * <<please update with a overview of the operation of this file>> 24 * 25 * END_DESC 26 */ 27 28 29 #include "float.h" 30 #include "dbl_float.h" 31 32 /* 33 * Double Precision Floating-point Multiply 34 */ 35 36 int 37 dbl_fmpy( 38 dbl_floating_point *srcptr1, 39 dbl_floating_point *srcptr2, 40 dbl_floating_point *dstptr, 41 unsigned int *status) 42 { 43 register unsigned int opnd1p1, opnd1p2, opnd2p1, opnd2p2; 44 register unsigned int opnd3p1, opnd3p2, resultp1, resultp2; 45 register int dest_exponent, count; 46 register boolean inexact = FALSE, guardbit = FALSE, stickybit = FALSE; 47 boolean is_tiny; 48 49 Dbl_copyfromptr(srcptr1,opnd1p1,opnd1p2); 50 Dbl_copyfromptr(srcptr2,opnd2p1,opnd2p2); 51 52 /* 53 * set sign bit of result 54 */ 55 if (Dbl_sign(opnd1p1) ^ Dbl_sign(opnd2p1)) 56 Dbl_setnegativezerop1(resultp1); 57 else Dbl_setzerop1(resultp1); 58 /* 59 * check first operand for NaN's or infinity 60 */ 61 if (Dbl_isinfinity_exponent(opnd1p1)) { 62 if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) { 63 if (Dbl_isnotnan(opnd2p1,opnd2p2)) { 64 if (Dbl_iszero_exponentmantissa(opnd2p1,opnd2p2)) { 65 /* 66 * invalid since operands are infinity 67 * and zero 68 */ 69 if (Is_invalidtrap_enabled()) 70 return(INVALIDEXCEPTION); 71 Set_invalidflag(); 72 Dbl_makequietnan(resultp1,resultp2); 73 Dbl_copytoptr(resultp1,resultp2,dstptr); 74 return(NOEXCEPTION); 75 } 76 /* 77 * return infinity 78 */ 79 Dbl_setinfinity_exponentmantissa(resultp1,resultp2); 80 Dbl_copytoptr(resultp1,resultp2,dstptr); 81 return(NOEXCEPTION); 82 } 83 } 84 else { 85 /* 86 * is NaN; signaling or quiet? 87 */ 88 if (Dbl_isone_signaling(opnd1p1)) { 89 /* trap if INVALIDTRAP enabled */ 90 if (Is_invalidtrap_enabled()) 91 return(INVALIDEXCEPTION); 92 /* make NaN quiet */ 93 Set_invalidflag(); 94 Dbl_set_quiet(opnd1p1); 95 } 96 /* 97 * is second operand a signaling NaN? 98 */ 99 else if (Dbl_is_signalingnan(opnd2p1)) { 100 /* trap if INVALIDTRAP enabled */ 101 if (Is_invalidtrap_enabled()) 102 return(INVALIDEXCEPTION); 103 /* make NaN quiet */ 104 Set_invalidflag(); 105 Dbl_set_quiet(opnd2p1); 106 Dbl_copytoptr(opnd2p1,opnd2p2,dstptr); 107 return(NOEXCEPTION); 108 } 109 /* 110 * return quiet NaN 111 */ 112 Dbl_copytoptr(opnd1p1,opnd1p2,dstptr); 113 return(NOEXCEPTION); 114 } 115 } 116 /* 117 * check second operand for NaN's or infinity 118 */ 119 if (Dbl_isinfinity_exponent(opnd2p1)) { 120 if (Dbl_iszero_mantissa(opnd2p1,opnd2p2)) { 121 if (Dbl_iszero_exponentmantissa(opnd1p1,opnd1p2)) { 122 /* invalid since operands are zero & infinity */ 123 if (Is_invalidtrap_enabled()) 124 return(INVALIDEXCEPTION); 125 Set_invalidflag(); 126 Dbl_makequietnan(opnd2p1,opnd2p2); 127 Dbl_copytoptr(opnd2p1,opnd2p2,dstptr); 128 return(NOEXCEPTION); 129 } 130 /* 131 * return infinity 132 */ 133 Dbl_setinfinity_exponentmantissa(resultp1,resultp2); 134 Dbl_copytoptr(resultp1,resultp2,dstptr); 135 return(NOEXCEPTION); 136 } 137 /* 138 * is NaN; signaling or quiet? 139 */ 140 if (Dbl_isone_signaling(opnd2p1)) { 141 /* trap if INVALIDTRAP enabled */ 142 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 143 /* make NaN quiet */ 144 Set_invalidflag(); 145 Dbl_set_quiet(opnd2p1); 146 } 147 /* 148 * return quiet NaN 149 */ 150 Dbl_copytoptr(opnd2p1,opnd2p2,dstptr); 151 return(NOEXCEPTION); 152 } 153 /* 154 * Generate exponent 155 */ 156 dest_exponent = Dbl_exponent(opnd1p1) + Dbl_exponent(opnd2p1) -DBL_BIAS; 157 158 /* 159 * Generate mantissa 160 */ 161 if (Dbl_isnotzero_exponent(opnd1p1)) { 162 /* set hidden bit */ 163 Dbl_clear_signexponent_set_hidden(opnd1p1); 164 } 165 else { 166 /* check for zero */ 167 if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) { 168 Dbl_setzero_exponentmantissa(resultp1,resultp2); 169 Dbl_copytoptr(resultp1,resultp2,dstptr); 170 return(NOEXCEPTION); 171 } 172 /* is denormalized, adjust exponent */ 173 Dbl_clear_signexponent(opnd1p1); 174 Dbl_leftshiftby1(opnd1p1,opnd1p2); 175 Dbl_normalize(opnd1p1,opnd1p2,dest_exponent); 176 } 177 /* opnd2 needs to have hidden bit set with msb in hidden bit */ 178 if (Dbl_isnotzero_exponent(opnd2p1)) { 179 Dbl_clear_signexponent_set_hidden(opnd2p1); 180 } 181 else { 182 /* check for zero */ 183 if (Dbl_iszero_mantissa(opnd2p1,opnd2p2)) { 184 Dbl_setzero_exponentmantissa(resultp1,resultp2); 185 Dbl_copytoptr(resultp1,resultp2,dstptr); 186 return(NOEXCEPTION); 187 } 188 /* is denormalized; want to normalize */ 189 Dbl_clear_signexponent(opnd2p1); 190 Dbl_leftshiftby1(opnd2p1,opnd2p2); 191 Dbl_normalize(opnd2p1,opnd2p2,dest_exponent); 192 } 193 194 /* Multiply two source mantissas together */ 195 196 /* make room for guard bits */ 197 Dbl_leftshiftby7(opnd2p1,opnd2p2); 198 Dbl_setzero(opnd3p1,opnd3p2); 199 /* 200 * Four bits at a time are inspected in each loop, and a 201 * simple shift and add multiply algorithm is used. 202 */ 203 for (count=1;count<=DBL_P;count+=4) { 204 stickybit |= Dlow4p2(opnd3p2); 205 Dbl_rightshiftby4(opnd3p1,opnd3p2); 206 if (Dbit28p2(opnd1p2)) { 207 /* Twoword_add should be an ADDC followed by an ADD. */ 208 Twoword_add(opnd3p1, opnd3p2, opnd2p1<<3 | opnd2p2>>29, 209 opnd2p2<<3); 210 } 211 if (Dbit29p2(opnd1p2)) { 212 Twoword_add(opnd3p1, opnd3p2, opnd2p1<<2 | opnd2p2>>30, 213 opnd2p2<<2); 214 } 215 if (Dbit30p2(opnd1p2)) { 216 Twoword_add(opnd3p1, opnd3p2, opnd2p1<<1 | opnd2p2>>31, 217 opnd2p2<<1); 218 } 219 if (Dbit31p2(opnd1p2)) { 220 Twoword_add(opnd3p1, opnd3p2, opnd2p1, opnd2p2); 221 } 222 Dbl_rightshiftby4(opnd1p1,opnd1p2); 223 } 224 if (Dbit3p1(opnd3p1)==0) { 225 Dbl_leftshiftby1(opnd3p1,opnd3p2); 226 } 227 else { 228 /* result mantissa >= 2. */ 229 dest_exponent++; 230 } 231 /* check for denormalized result */ 232 while (Dbit3p1(opnd3p1)==0) { 233 Dbl_leftshiftby1(opnd3p1,opnd3p2); 234 dest_exponent--; 235 } 236 /* 237 * check for guard, sticky and inexact bits 238 */ 239 stickybit |= Dallp2(opnd3p2) << 25; 240 guardbit = (Dallp2(opnd3p2) << 24) >> 31; 241 inexact = guardbit | stickybit; 242 243 /* align result mantissa */ 244 Dbl_rightshiftby8(opnd3p1,opnd3p2); 245 246 /* 247 * round result 248 */ 249 if (inexact && (dest_exponent>0 || Is_underflowtrap_enabled())) { 250 Dbl_clear_signexponent(opnd3p1); 251 switch (Rounding_mode()) { 252 case ROUNDPLUS: 253 if (Dbl_iszero_sign(resultp1)) 254 Dbl_increment(opnd3p1,opnd3p2); 255 break; 256 case ROUNDMINUS: 257 if (Dbl_isone_sign(resultp1)) 258 Dbl_increment(opnd3p1,opnd3p2); 259 break; 260 case ROUNDNEAREST: 261 if (guardbit) { 262 if (stickybit || Dbl_isone_lowmantissap2(opnd3p2)) 263 Dbl_increment(opnd3p1,opnd3p2); 264 } 265 } 266 if (Dbl_isone_hidden(opnd3p1)) dest_exponent++; 267 } 268 Dbl_set_mantissa(resultp1,resultp2,opnd3p1,opnd3p2); 269 270 /* 271 * Test for overflow 272 */ 273 if (dest_exponent >= DBL_INFINITY_EXPONENT) { 274 /* trap if OVERFLOWTRAP enabled */ 275 if (Is_overflowtrap_enabled()) { 276 /* 277 * Adjust bias of result 278 */ 279 Dbl_setwrapped_exponent(resultp1,dest_exponent,ovfl); 280 Dbl_copytoptr(resultp1,resultp2,dstptr); 281 if (inexact) 282 if (Is_inexacttrap_enabled()) 283 return (OVERFLOWEXCEPTION | INEXACTEXCEPTION); 284 else Set_inexactflag(); 285 return (OVERFLOWEXCEPTION); 286 } 287 inexact = TRUE; 288 Set_overflowflag(); 289 /* set result to infinity or largest number */ 290 Dbl_setoverflow(resultp1,resultp2); 291 } 292 /* 293 * Test for underflow 294 */ 295 else if (dest_exponent <= 0) { 296 /* trap if UNDERFLOWTRAP enabled */ 297 if (Is_underflowtrap_enabled()) { 298 /* 299 * Adjust bias of result 300 */ 301 Dbl_setwrapped_exponent(resultp1,dest_exponent,unfl); 302 Dbl_copytoptr(resultp1,resultp2,dstptr); 303 if (inexact) 304 if (Is_inexacttrap_enabled()) 305 return (UNDERFLOWEXCEPTION | INEXACTEXCEPTION); 306 else Set_inexactflag(); 307 return (UNDERFLOWEXCEPTION); 308 } 309 310 /* Determine if should set underflow flag */ 311 is_tiny = TRUE; 312 if (dest_exponent == 0 && inexact) { 313 switch (Rounding_mode()) { 314 case ROUNDPLUS: 315 if (Dbl_iszero_sign(resultp1)) { 316 Dbl_increment(opnd3p1,opnd3p2); 317 if (Dbl_isone_hiddenoverflow(opnd3p1)) 318 is_tiny = FALSE; 319 Dbl_decrement(opnd3p1,opnd3p2); 320 } 321 break; 322 case ROUNDMINUS: 323 if (Dbl_isone_sign(resultp1)) { 324 Dbl_increment(opnd3p1,opnd3p2); 325 if (Dbl_isone_hiddenoverflow(opnd3p1)) 326 is_tiny = FALSE; 327 Dbl_decrement(opnd3p1,opnd3p2); 328 } 329 break; 330 case ROUNDNEAREST: 331 if (guardbit && (stickybit || 332 Dbl_isone_lowmantissap2(opnd3p2))) { 333 Dbl_increment(opnd3p1,opnd3p2); 334 if (Dbl_isone_hiddenoverflow(opnd3p1)) 335 is_tiny = FALSE; 336 Dbl_decrement(opnd3p1,opnd3p2); 337 } 338 break; 339 } 340 } 341 342 /* 343 * denormalize result or set to signed zero 344 */ 345 stickybit = inexact; 346 Dbl_denormalize(opnd3p1,opnd3p2,dest_exponent,guardbit, 347 stickybit,inexact); 348 349 /* return zero or smallest number */ 350 if (inexact) { 351 switch (Rounding_mode()) { 352 case ROUNDPLUS: 353 if (Dbl_iszero_sign(resultp1)) { 354 Dbl_increment(opnd3p1,opnd3p2); 355 } 356 break; 357 case ROUNDMINUS: 358 if (Dbl_isone_sign(resultp1)) { 359 Dbl_increment(opnd3p1,opnd3p2); 360 } 361 break; 362 case ROUNDNEAREST: 363 if (guardbit && (stickybit || 364 Dbl_isone_lowmantissap2(opnd3p2))) { 365 Dbl_increment(opnd3p1,opnd3p2); 366 } 367 break; 368 } 369 if (is_tiny) Set_underflowflag(); 370 } 371 Dbl_set_exponentmantissa(resultp1,resultp2,opnd3p1,opnd3p2); 372 } 373 else Dbl_set_exponent(resultp1,dest_exponent); 374 /* check for inexact */ 375 Dbl_copytoptr(resultp1,resultp2,dstptr); 376 if (inexact) { 377 if (Is_inexacttrap_enabled()) return(INEXACTEXCEPTION); 378 else Set_inexactflag(); 379 } 380 return(NOEXCEPTION); 381 } 382