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/dfdiv.c $Revision: 1.1 $ 26 * 27 * Purpose: 28 * Double Precision Floating-point Divide 29 * 30 * External Interfaces: 31 * dbl_fdiv(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 "dbl_float.h" 44 45 /* 46 * Double Precision Floating-point Divide 47 */ 48 49 int 50 dbl_fdiv (dbl_floating_point * srcptr1, dbl_floating_point * srcptr2, 51 dbl_floating_point * dstptr, unsigned int *status) 52 { 53 register unsigned int opnd1p1, opnd1p2, opnd2p1, opnd2p2; 54 register unsigned int opnd3p1, opnd3p2, resultp1, resultp2; 55 register int dest_exponent, count; 56 register boolean inexact = FALSE, guardbit = FALSE, stickybit = FALSE; 57 boolean is_tiny; 58 59 Dbl_copyfromptr(srcptr1,opnd1p1,opnd1p2); 60 Dbl_copyfromptr(srcptr2,opnd2p1,opnd2p2); 61 /* 62 * set sign bit of result 63 */ 64 if (Dbl_sign(opnd1p1) ^ Dbl_sign(opnd2p1)) 65 Dbl_setnegativezerop1(resultp1); 66 else Dbl_setzerop1(resultp1); 67 /* 68 * check first operand for NaN's or infinity 69 */ 70 if (Dbl_isinfinity_exponent(opnd1p1)) { 71 if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) { 72 if (Dbl_isnotnan(opnd2p1,opnd2p2)) { 73 if (Dbl_isinfinity(opnd2p1,opnd2p2)) { 74 /* 75 * invalid since both operands 76 * are infinity 77 */ 78 if (Is_invalidtrap_enabled()) 79 return(INVALIDEXCEPTION); 80 Set_invalidflag(); 81 Dbl_makequietnan(resultp1,resultp2); 82 Dbl_copytoptr(resultp1,resultp2,dstptr); 83 return(NOEXCEPTION); 84 } 85 /* 86 * return infinity 87 */ 88 Dbl_setinfinity_exponentmantissa(resultp1,resultp2); 89 Dbl_copytoptr(resultp1,resultp2,dstptr); 90 return(NOEXCEPTION); 91 } 92 } 93 else { 94 /* 95 * is NaN; signaling or quiet? 96 */ 97 if (Dbl_isone_signaling(opnd1p1)) { 98 /* trap if INVALIDTRAP enabled */ 99 if (Is_invalidtrap_enabled()) 100 return(INVALIDEXCEPTION); 101 /* make NaN quiet */ 102 Set_invalidflag(); 103 Dbl_set_quiet(opnd1p1); 104 } 105 /* 106 * is second operand a signaling NaN? 107 */ 108 else if (Dbl_is_signalingnan(opnd2p1)) { 109 /* trap if INVALIDTRAP enabled */ 110 if (Is_invalidtrap_enabled()) 111 return(INVALIDEXCEPTION); 112 /* make NaN quiet */ 113 Set_invalidflag(); 114 Dbl_set_quiet(opnd2p1); 115 Dbl_copytoptr(opnd2p1,opnd2p2,dstptr); 116 return(NOEXCEPTION); 117 } 118 /* 119 * return quiet NaN 120 */ 121 Dbl_copytoptr(opnd1p1,opnd1p2,dstptr); 122 return(NOEXCEPTION); 123 } 124 } 125 /* 126 * check second operand for NaN's or infinity 127 */ 128 if (Dbl_isinfinity_exponent(opnd2p1)) { 129 if (Dbl_iszero_mantissa(opnd2p1,opnd2p2)) { 130 /* 131 * return zero 132 */ 133 Dbl_setzero_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 * check for division by zero 155 */ 156 if (Dbl_iszero_exponentmantissa(opnd2p1,opnd2p2)) { 157 if (Dbl_iszero_exponentmantissa(opnd1p1,opnd1p2)) { 158 /* invalid since both operands are zero */ 159 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 160 Set_invalidflag(); 161 Dbl_makequietnan(resultp1,resultp2); 162 Dbl_copytoptr(resultp1,resultp2,dstptr); 163 return(NOEXCEPTION); 164 } 165 if (Is_divisionbyzerotrap_enabled()) 166 return(DIVISIONBYZEROEXCEPTION); 167 Set_divisionbyzeroflag(); 168 Dbl_setinfinity_exponentmantissa(resultp1,resultp2); 169 Dbl_copytoptr(resultp1,resultp2,dstptr); 170 return(NOEXCEPTION); 171 } 172 /* 173 * Generate exponent 174 */ 175 dest_exponent = Dbl_exponent(opnd1p1) - Dbl_exponent(opnd2p1) + DBL_BIAS; 176 177 /* 178 * Generate mantissa 179 */ 180 if (Dbl_isnotzero_exponent(opnd1p1)) { 181 /* set hidden bit */ 182 Dbl_clear_signexponent_set_hidden(opnd1p1); 183 } 184 else { 185 /* check for zero */ 186 if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) { 187 Dbl_setzero_exponentmantissa(resultp1,resultp2); 188 Dbl_copytoptr(resultp1,resultp2,dstptr); 189 return(NOEXCEPTION); 190 } 191 /* is denormalized, want to normalize */ 192 Dbl_clear_signexponent(opnd1p1); 193 Dbl_leftshiftby1(opnd1p1,opnd1p2); 194 Dbl_normalize(opnd1p1,opnd1p2,dest_exponent); 195 } 196 /* opnd2 needs to have hidden bit set with msb in hidden bit */ 197 if (Dbl_isnotzero_exponent(opnd2p1)) { 198 Dbl_clear_signexponent_set_hidden(opnd2p1); 199 } 200 else { 201 /* is denormalized; want to normalize */ 202 Dbl_clear_signexponent(opnd2p1); 203 Dbl_leftshiftby1(opnd2p1,opnd2p2); 204 while (Dbl_iszero_hiddenhigh7mantissa(opnd2p1)) { 205 dest_exponent+=8; 206 Dbl_leftshiftby8(opnd2p1,opnd2p2); 207 } 208 if (Dbl_iszero_hiddenhigh3mantissa(opnd2p1)) { 209 dest_exponent+=4; 210 Dbl_leftshiftby4(opnd2p1,opnd2p2); 211 } 212 while (Dbl_iszero_hidden(opnd2p1)) { 213 dest_exponent++; 214 Dbl_leftshiftby1(opnd2p1,opnd2p2); 215 } 216 } 217 218 /* Divide the source mantissas */ 219 220 /* 221 * A non-restoring divide algorithm is used. 222 */ 223 Twoword_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2); 224 Dbl_setzero(opnd3p1,opnd3p2); 225 for (count=1; count <= DBL_P && (opnd1p1 || opnd1p2); count++) { 226 Dbl_leftshiftby1(opnd1p1,opnd1p2); 227 Dbl_leftshiftby1(opnd3p1,opnd3p2); 228 if (Dbl_iszero_sign(opnd1p1)) { 229 Dbl_setone_lowmantissap2(opnd3p2); 230 Twoword_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2); 231 } 232 else { 233 Twoword_add(opnd1p1, opnd1p2, opnd2p1, opnd2p2); 234 } 235 } 236 if (count <= DBL_P) { 237 Dbl_leftshiftby1(opnd3p1,opnd3p2); 238 Dbl_setone_lowmantissap2(opnd3p2); 239 Dbl_leftshift(opnd3p1,opnd3p2,(DBL_P-count)); 240 if (Dbl_iszero_hidden(opnd3p1)) { 241 Dbl_leftshiftby1(opnd3p1,opnd3p2); 242 dest_exponent--; 243 } 244 } 245 else { 246 if (Dbl_iszero_hidden(opnd3p1)) { 247 /* need to get one more bit of result */ 248 Dbl_leftshiftby1(opnd1p1,opnd1p2); 249 Dbl_leftshiftby1(opnd3p1,opnd3p2); 250 if (Dbl_iszero_sign(opnd1p1)) { 251 Dbl_setone_lowmantissap2(opnd3p2); 252 Twoword_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2); 253 } 254 else { 255 Twoword_add(opnd1p1,opnd1p2,opnd2p1,opnd2p2); 256 } 257 dest_exponent--; 258 } 259 if (Dbl_iszero_sign(opnd1p1)) guardbit = TRUE; 260 stickybit = Dbl_allp1(opnd1p1) || Dbl_allp2(opnd1p2); 261 } 262 inexact = guardbit | stickybit; 263 264 /* 265 * round result 266 */ 267 if (inexact && (dest_exponent > 0 || Is_underflowtrap_enabled())) { 268 Dbl_clear_signexponent(opnd3p1); 269 switch (Rounding_mode()) { 270 case ROUNDPLUS: 271 if (Dbl_iszero_sign(resultp1)) 272 Dbl_increment(opnd3p1,opnd3p2); 273 break; 274 case ROUNDMINUS: 275 if (Dbl_isone_sign(resultp1)) 276 Dbl_increment(opnd3p1,opnd3p2); 277 break; 278 case ROUNDNEAREST: 279 if (guardbit && (stickybit || 280 Dbl_isone_lowmantissap2(opnd3p2))) { 281 Dbl_increment(opnd3p1,opnd3p2); 282 } 283 } 284 if (Dbl_isone_hidden(opnd3p1)) dest_exponent++; 285 } 286 Dbl_set_mantissa(resultp1,resultp2,opnd3p1,opnd3p2); 287 288 /* 289 * Test for overflow 290 */ 291 if (dest_exponent >= DBL_INFINITY_EXPONENT) { 292 /* trap if OVERFLOWTRAP enabled */ 293 if (Is_overflowtrap_enabled()) { 294 /* 295 * Adjust bias of result 296 */ 297 Dbl_setwrapped_exponent(resultp1,dest_exponent,ovfl); 298 Dbl_copytoptr(resultp1,resultp2,dstptr); 299 if (inexact) 300 if (Is_inexacttrap_enabled()) 301 return(OVERFLOWEXCEPTION | INEXACTEXCEPTION); 302 else Set_inexactflag(); 303 return(OVERFLOWEXCEPTION); 304 } 305 Set_overflowflag(); 306 /* set result to infinity or largest number */ 307 Dbl_setoverflow(resultp1,resultp2); 308 inexact = TRUE; 309 } 310 /* 311 * Test for underflow 312 */ 313 else if (dest_exponent <= 0) { 314 /* trap if UNDERFLOWTRAP enabled */ 315 if (Is_underflowtrap_enabled()) { 316 /* 317 * Adjust bias of result 318 */ 319 Dbl_setwrapped_exponent(resultp1,dest_exponent,unfl); 320 Dbl_copytoptr(resultp1,resultp2,dstptr); 321 if (inexact) 322 if (Is_inexacttrap_enabled()) 323 return(UNDERFLOWEXCEPTION | INEXACTEXCEPTION); 324 else Set_inexactflag(); 325 return(UNDERFLOWEXCEPTION); 326 } 327 328 /* Determine if should set underflow flag */ 329 is_tiny = TRUE; 330 if (dest_exponent == 0 && inexact) { 331 switch (Rounding_mode()) { 332 case ROUNDPLUS: 333 if (Dbl_iszero_sign(resultp1)) { 334 Dbl_increment(opnd3p1,opnd3p2); 335 if (Dbl_isone_hiddenoverflow(opnd3p1)) 336 is_tiny = FALSE; 337 Dbl_decrement(opnd3p1,opnd3p2); 338 } 339 break; 340 case ROUNDMINUS: 341 if (Dbl_isone_sign(resultp1)) { 342 Dbl_increment(opnd3p1,opnd3p2); 343 if (Dbl_isone_hiddenoverflow(opnd3p1)) 344 is_tiny = FALSE; 345 Dbl_decrement(opnd3p1,opnd3p2); 346 } 347 break; 348 case ROUNDNEAREST: 349 if (guardbit && (stickybit || 350 Dbl_isone_lowmantissap2(opnd3p2))) { 351 Dbl_increment(opnd3p1,opnd3p2); 352 if (Dbl_isone_hiddenoverflow(opnd3p1)) 353 is_tiny = FALSE; 354 Dbl_decrement(opnd3p1,opnd3p2); 355 } 356 break; 357 } 358 } 359 360 /* 361 * denormalize result or set to signed zero 362 */ 363 stickybit = inexact; 364 Dbl_denormalize(opnd3p1,opnd3p2,dest_exponent,guardbit, 365 stickybit,inexact); 366 367 /* return rounded number */ 368 if (inexact) { 369 switch (Rounding_mode()) { 370 case ROUNDPLUS: 371 if (Dbl_iszero_sign(resultp1)) { 372 Dbl_increment(opnd3p1,opnd3p2); 373 } 374 break; 375 case ROUNDMINUS: 376 if (Dbl_isone_sign(resultp1)) { 377 Dbl_increment(opnd3p1,opnd3p2); 378 } 379 break; 380 case ROUNDNEAREST: 381 if (guardbit && (stickybit || 382 Dbl_isone_lowmantissap2(opnd3p2))) { 383 Dbl_increment(opnd3p1,opnd3p2); 384 } 385 break; 386 } 387 if (is_tiny) Set_underflowflag(); 388 } 389 Dbl_set_exponentmantissa(resultp1,resultp2,opnd3p1,opnd3p2); 390 } 391 else Dbl_set_exponent(resultp1,dest_exponent); 392 Dbl_copytoptr(resultp1,resultp2,dstptr); 393 394 /* check for inexact */ 395 if (inexact) { 396 if (Is_inexacttrap_enabled()) return(INEXACTEXCEPTION); 397 else Set_inexactflag(); 398 } 399 return(NOEXCEPTION); 400 } 401