| 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/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 | } |