1 /* $OpenBSD: sfsqrt.c,v 1.7 2003/04/10 17:27:59 mickey Exp $ */ 2 /* 3 (c) Copyright 1986 HEWLETT-PACKARD COMPANY 4 To anyone who acknowledges that this file is provided "AS IS" 5 without any express or implied warranty: 6 permission to use, copy, modify, and distribute this file 7 for any purpose is hereby granted without fee, provided that 8 the above copyright notice and this notice appears in all 9 copies, and that the name of Hewlett-Packard Company not be 10 used in advertising or publicity pertaining to distribution 11 of the software without specific, written prior permission. 12 Hewlett-Packard Company makes no representations about the 13 suitability of this software for any purpose. 14 */ 15 /* @(#)sfsqrt.c: Revision: 1.9.88.1 Date: 93/12/07 15:07:13 */ 16 17 #include "float.h" 18 #include "sgl_float.h" 19 20 /* 21 * Single Floating-point Square Root 22 */ 23 24 /*ARGSUSED*/ 25 int 26 sgl_fsqrt(srcptr, null, dstptr, status) 27 sgl_floating_point *srcptr, *null, *dstptr; 28 unsigned int *status; 29 { 30 register unsigned int src, result; 31 register int src_exponent, newbit, sum; 32 register int guardbit = FALSE, even_exponent; 33 34 src = *srcptr; 35 /* 36 * check source operand for NaN or infinity 37 */ 38 if ((src_exponent = Sgl_exponent(src)) == SGL_INFINITY_EXPONENT) { 39 /* 40 * is signaling NaN? 41 */ 42 if (Sgl_isone_signaling(src)) { 43 /* trap if INVALIDTRAP enabled */ 44 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 45 /* make NaN quiet */ 46 Set_invalidflag(); 47 Sgl_set_quiet(src); 48 } 49 /* 50 * Return quiet NaN or positive infinity. 51 * Fall thru to negative test if negative infinity. 52 */ 53 if (Sgl_iszero_sign(src) || Sgl_isnotzero_mantissa(src)) { 54 *dstptr = src; 55 return(NOEXCEPTION); 56 } 57 } 58 59 /* 60 * check for zero source operand 61 */ 62 if (Sgl_iszero_exponentmantissa(src)) { 63 *dstptr = src; 64 return(NOEXCEPTION); 65 } 66 67 /* 68 * check for negative source operand 69 */ 70 if (Sgl_isone_sign(src)) { 71 /* trap if INVALIDTRAP enabled */ 72 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 73 /* make NaN quiet */ 74 Set_invalidflag(); 75 Sgl_makequietnan(src); 76 *dstptr = src; 77 return(NOEXCEPTION); 78 } 79 80 /* 81 * Generate result 82 */ 83 if (src_exponent > 0) { 84 even_exponent = Sgl_hidden(src); 85 Sgl_clear_signexponent_set_hidden(src); 86 } 87 else { 88 /* normalize operand */ 89 Sgl_clear_signexponent(src); 90 src_exponent++; 91 Sgl_normalize(src,src_exponent); 92 even_exponent = src_exponent & 1; 93 } 94 if (even_exponent) { 95 /* exponent is even */ 96 /* Add comment here. Explain why odd exponent needs correction */ 97 Sgl_leftshiftby1(src); 98 } 99 /* 100 * Add comment here. Explain following algorithm. 101 * 102 * Trust me, it works. 103 * 104 */ 105 Sgl_setzero(result); 106 newbit = 1 << SGL_P; 107 while (newbit && Sgl_isnotzero(src)) { 108 Sgl_addition(result,newbit,sum); 109 if(sum <= Sgl_all(src)) { 110 /* update result */ 111 Sgl_addition(result,(newbit<<1),result); 112 Sgl_subtract(src,sum,src); 113 } 114 Sgl_rightshiftby1(newbit); 115 Sgl_leftshiftby1(src); 116 } 117 /* correct exponent for pre-shift */ 118 if (even_exponent) { 119 Sgl_rightshiftby1(result); 120 } 121 122 /* check for inexact */ 123 if (Sgl_isnotzero(src)) { 124 if (!even_exponent & Sgl_islessthan(result,src)) 125 Sgl_increment(result); 126 guardbit = Sgl_lowmantissa(result); 127 Sgl_rightshiftby1(result); 128 129 /* now round result */ 130 switch (Rounding_mode()) { 131 case ROUNDPLUS: 132 Sgl_increment(result); 133 break; 134 case ROUNDNEAREST: 135 /* stickybit is always true, so guardbit 136 * is enough to determine rounding */ 137 if (guardbit) { 138 Sgl_increment(result); 139 } 140 break; 141 } 142 /* increment result exponent by 1 if mantissa overflowed */ 143 if (Sgl_isone_hiddenoverflow(result)) src_exponent+=2; 144 145 if (Is_inexacttrap_enabled()) { 146 Sgl_set_exponent(result, 147 ((src_exponent-SGL_BIAS)>>1)+SGL_BIAS); 148 *dstptr = result; 149 return(INEXACTEXCEPTION); 150 } 151 else Set_inexactflag(); 152 } 153 else { 154 Sgl_rightshiftby1(result); 155 } 156 Sgl_set_exponent(result,((src_exponent-SGL_BIAS)>>1)+SGL_BIAS); 157 *dstptr = result; 158 return(NOEXCEPTION); 159 } 160