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/sfsqrt.c $Revision: 1.1 $ 26 * 27 * Purpose: 28 * Single Floating-point Square Root 29 * 30 * External Interfaces: 31 * sgl_fsqrt(srcptr,nullptr,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 Floating-point Square Root 47 */ 48 49 /*ARGSUSED*/ 50 unsigned int 51 sgl_fsqrt( 52 sgl_floating_point *srcptr, 53 unsigned int *nullptr, 54 sgl_floating_point *dstptr, 55 unsigned int *status) 56 { 57 register unsigned int src, result; 58 register int src_exponent; 59 register unsigned int newbit, sum; 60 register boolean guardbit = FALSE, even_exponent; 61 62 src = *srcptr; 63 /* 64 * check source operand for NaN or infinity 65 */ 66 if ((src_exponent = Sgl_exponent(src)) == SGL_INFINITY_EXPONENT) { 67 /* 68 * is signaling NaN? 69 */ 70 if (Sgl_isone_signaling(src)) { 71 /* trap if INVALIDTRAP enabled */ 72 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 73 /* make NaN quiet */ 74 Set_invalidflag(); 75 Sgl_set_quiet(src); 76 } 77 /* 78 * Return quiet NaN or positive infinity. 79 * Fall through to negative test if negative infinity. 80 */ 81 if (Sgl_iszero_sign(src) || Sgl_isnotzero_mantissa(src)) { 82 *dstptr = src; 83 return(NOEXCEPTION); 84 } 85 } 86 87 /* 88 * check for zero source operand 89 */ 90 if (Sgl_iszero_exponentmantissa(src)) { 91 *dstptr = src; 92 return(NOEXCEPTION); 93 } 94 95 /* 96 * check for negative source operand 97 */ 98 if (Sgl_isone_sign(src)) { 99 /* trap if INVALIDTRAP enabled */ 100 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 101 /* make NaN quiet */ 102 Set_invalidflag(); 103 Sgl_makequietnan(src); 104 *dstptr = src; 105 return(NOEXCEPTION); 106 } 107 108 /* 109 * Generate result 110 */ 111 if (src_exponent > 0) { 112 even_exponent = Sgl_hidden(src); 113 Sgl_clear_signexponent_set_hidden(src); 114 } 115 else { 116 /* normalize operand */ 117 Sgl_clear_signexponent(src); 118 src_exponent++; 119 Sgl_normalize(src,src_exponent); 120 even_exponent = src_exponent & 1; 121 } 122 if (even_exponent) { 123 /* exponent is even */ 124 /* Add comment here. Explain why odd exponent needs correction */ 125 Sgl_leftshiftby1(src); 126 } 127 /* 128 * Add comment here. Explain following algorithm. 129 * 130 * Trust me, it works. 131 * 132 */ 133 Sgl_setzero(result); 134 newbit = 1 << SGL_P; 135 while (newbit && Sgl_isnotzero(src)) { 136 Sgl_addition(result,newbit,sum); 137 if(sum <= Sgl_all(src)) { 138 /* update result */ 139 Sgl_addition(result,(newbit<<1),result); 140 Sgl_subtract(src,sum,src); 141 } 142 Sgl_rightshiftby1(newbit); 143 Sgl_leftshiftby1(src); 144 } 145 /* correct exponent for pre-shift */ 146 if (even_exponent) { 147 Sgl_rightshiftby1(result); 148 } 149 150 /* check for inexact */ 151 if (Sgl_isnotzero(src)) { 152 if (!even_exponent && Sgl_islessthan(result,src)) 153 Sgl_increment(result); 154 guardbit = Sgl_lowmantissa(result); 155 Sgl_rightshiftby1(result); 156 157 /* now round result */ 158 switch (Rounding_mode()) { 159 case ROUNDPLUS: 160 Sgl_increment(result); 161 break; 162 case ROUNDNEAREST: 163 /* stickybit is always true, so guardbit 164 * is enough to determine rounding */ 165 if (guardbit) { 166 Sgl_increment(result); 167 } 168 break; 169 } 170 /* increment result exponent by 1 if mantissa overflowed */ 171 if (Sgl_isone_hiddenoverflow(result)) src_exponent+=2; 172 173 if (Is_inexacttrap_enabled()) { 174 Sgl_set_exponent(result, 175 ((src_exponent-SGL_BIAS)>>1)+SGL_BIAS); 176 *dstptr = result; 177 return(INEXACTEXCEPTION); 178 } 179 else Set_inexactflag(); 180 } 181 else { 182 Sgl_rightshiftby1(result); 183 } 184 Sgl_set_exponent(result,((src_exponent-SGL_BIAS)>>1)+SGL_BIAS); 185 *dstptr = result; 186 return(NOEXCEPTION); 187 } 188