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/dfsqrt.c $Revision: 1.1 $ 26 * 27 * Purpose: 28 * Double Floating-point Square Root 29 * 30 * External Interfaces: 31 * dbl_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 "dbl_float.h" 44 45 /* 46 * Double Floating-point Square Root 47 */ 48 49 /*ARGSUSED*/ 50 unsigned int 51 dbl_fsqrt( 52 dbl_floating_point *srcptr, 53 unsigned int *nullptr, 54 dbl_floating_point *dstptr, 55 unsigned int *status) 56 { 57 register unsigned int srcp1, srcp2, resultp1, resultp2; 58 register unsigned int newbitp1, newbitp2, sump1, sump2; 59 register int src_exponent; 60 register boolean guardbit = FALSE, even_exponent; 61 62 Dbl_copyfromptr(srcptr,srcp1,srcp2); 63 /* 64 * check source operand for NaN or infinity 65 */ 66 if ((src_exponent = Dbl_exponent(srcp1)) == DBL_INFINITY_EXPONENT) { 67 /* 68 * is signaling NaN? 69 */ 70 if (Dbl_isone_signaling(srcp1)) { 71 /* trap if INVALIDTRAP enabled */ 72 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 73 /* make NaN quiet */ 74 Set_invalidflag(); 75 Dbl_set_quiet(srcp1); 76 } 77 /* 78 * Return quiet NaN or positive infinity. 79 * Fall through to negative test if negative infinity. 80 */ 81 if (Dbl_iszero_sign(srcp1) || 82 Dbl_isnotzero_mantissa(srcp1,srcp2)) { 83 Dbl_copytoptr(srcp1,srcp2,dstptr); 84 return(NOEXCEPTION); 85 } 86 } 87 88 /* 89 * check for zero source operand 90 */ 91 if (Dbl_iszero_exponentmantissa(srcp1,srcp2)) { 92 Dbl_copytoptr(srcp1,srcp2,dstptr); 93 return(NOEXCEPTION); 94 } 95 96 /* 97 * check for negative source operand 98 */ 99 if (Dbl_isone_sign(srcp1)) { 100 /* trap if INVALIDTRAP enabled */ 101 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 102 /* make NaN quiet */ 103 Set_invalidflag(); 104 Dbl_makequietnan(srcp1,srcp2); 105 Dbl_copytoptr(srcp1,srcp2,dstptr); 106 return(NOEXCEPTION); 107 } 108 109 /* 110 * Generate result 111 */ 112 if (src_exponent > 0) { 113 even_exponent = Dbl_hidden(srcp1); 114 Dbl_clear_signexponent_set_hidden(srcp1); 115 } 116 else { 117 /* normalize operand */ 118 Dbl_clear_signexponent(srcp1); 119 src_exponent++; 120 Dbl_normalize(srcp1,srcp2,src_exponent); 121 even_exponent = src_exponent & 1; 122 } 123 if (even_exponent) { 124 /* exponent is even */ 125 /* Add comment here. Explain why odd exponent needs correction */ 126 Dbl_leftshiftby1(srcp1,srcp2); 127 } 128 /* 129 * Add comment here. Explain following algorithm. 130 * 131 * Trust me, it works. 132 * 133 */ 134 Dbl_setzero(resultp1,resultp2); 135 Dbl_allp1(newbitp1) = 1 << (DBL_P - 32); 136 Dbl_setzero_mantissap2(newbitp2); 137 while (Dbl_isnotzero(newbitp1,newbitp2) && Dbl_isnotzero(srcp1,srcp2)) { 138 Dbl_addition(resultp1,resultp2,newbitp1,newbitp2,sump1,sump2); 139 if(Dbl_isnotgreaterthan(sump1,sump2,srcp1,srcp2)) { 140 Dbl_leftshiftby1(newbitp1,newbitp2); 141 /* update result */ 142 Dbl_addition(resultp1,resultp2,newbitp1,newbitp2, 143 resultp1,resultp2); 144 Dbl_subtract(srcp1,srcp2,sump1,sump2,srcp1,srcp2); 145 Dbl_rightshiftby2(newbitp1,newbitp2); 146 } 147 else { 148 Dbl_rightshiftby1(newbitp1,newbitp2); 149 } 150 Dbl_leftshiftby1(srcp1,srcp2); 151 } 152 /* correct exponent for pre-shift */ 153 if (even_exponent) { 154 Dbl_rightshiftby1(resultp1,resultp2); 155 } 156 157 /* check for inexact */ 158 if (Dbl_isnotzero(srcp1,srcp2)) { 159 if (!even_exponent && Dbl_islessthan(resultp1,resultp2,srcp1,srcp2)) { 160 Dbl_increment(resultp1,resultp2); 161 } 162 guardbit = Dbl_lowmantissap2(resultp2); 163 Dbl_rightshiftby1(resultp1,resultp2); 164 165 /* now round result */ 166 switch (Rounding_mode()) { 167 case ROUNDPLUS: 168 Dbl_increment(resultp1,resultp2); 169 break; 170 case ROUNDNEAREST: 171 /* stickybit is always true, so guardbit 172 * is enough to determine rounding */ 173 if (guardbit) { 174 Dbl_increment(resultp1,resultp2); 175 } 176 break; 177 } 178 /* increment result exponent by 1 if mantissa overflowed */ 179 if (Dbl_isone_hiddenoverflow(resultp1)) src_exponent+=2; 180 181 if (Is_inexacttrap_enabled()) { 182 Dbl_set_exponent(resultp1, 183 ((src_exponent-DBL_BIAS)>>1)+DBL_BIAS); 184 Dbl_copytoptr(resultp1,resultp2,dstptr); 185 return(INEXACTEXCEPTION); 186 } 187 else Set_inexactflag(); 188 } 189 else { 190 Dbl_rightshiftby1(resultp1,resultp2); 191 } 192 Dbl_set_exponent(resultp1,((src_exponent-DBL_BIAS)>>1)+DBL_BIAS); 193 Dbl_copytoptr(resultp1,resultp2,dstptr); 194 return(NOEXCEPTION); 195 } 196