1 // SPDX-License-Identifier: GPL-2.0-or-later 2 /* 3 * Linux/PA-RISC Project (http://www.parisc-linux.org/) 4 * 5 * Floating-point emulation code 6 * Copyright (C) 2001 Hewlett-Packard (Paul Bame) <bame@debian.org> 7 */ 8 /* 9 * BEGIN_DESC 10 * 11 * File: 12 * @(#) pa/spmath/sfsqrt.c $Revision: 1.1 $ 13 * 14 * Purpose: 15 * Single Floating-point Square Root 16 * 17 * External Interfaces: 18 * sgl_fsqrt(srcptr,_nullptr,dstptr,status) 19 * 20 * Internal Interfaces: 21 * 22 * Theory: 23 * <<please update with a overview of the operation of this file>> 24 * 25 * END_DESC 26 */ 27 28 29 #include "float.h" 30 #include "sgl_float.h" 31 32 /* 33 * Single Floating-point Square Root 34 */ 35 36 /*ARGSUSED*/ 37 unsigned int 38 sgl_fsqrt( 39 sgl_floating_point *srcptr, 40 unsigned int *_nullptr, 41 sgl_floating_point *dstptr, 42 unsigned int *status) 43 { 44 register unsigned int src, result; 45 register int src_exponent; 46 register unsigned int newbit, sum; 47 register boolean guardbit = FALSE, even_exponent; 48 49 src = *srcptr; 50 /* 51 * check source operand for NaN or infinity 52 */ 53 if ((src_exponent = Sgl_exponent(src)) == SGL_INFINITY_EXPONENT) { 54 /* 55 * is signaling NaN? 56 */ 57 if (Sgl_isone_signaling(src)) { 58 /* trap if INVALIDTRAP enabled */ 59 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 60 /* make NaN quiet */ 61 Set_invalidflag(); 62 Sgl_set_quiet(src); 63 } 64 /* 65 * Return quiet NaN or positive infinity. 66 * Fall through to negative test if negative infinity. 67 */ 68 if (Sgl_iszero_sign(src) || Sgl_isnotzero_mantissa(src)) { 69 *dstptr = src; 70 return(NOEXCEPTION); 71 } 72 } 73 74 /* 75 * check for zero source operand 76 */ 77 if (Sgl_iszero_exponentmantissa(src)) { 78 *dstptr = src; 79 return(NOEXCEPTION); 80 } 81 82 /* 83 * check for negative source operand 84 */ 85 if (Sgl_isone_sign(src)) { 86 /* trap if INVALIDTRAP enabled */ 87 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION); 88 /* make NaN quiet */ 89 Set_invalidflag(); 90 Sgl_makequietnan(src); 91 *dstptr = src; 92 return(NOEXCEPTION); 93 } 94 95 /* 96 * Generate result 97 */ 98 if (src_exponent > 0) { 99 even_exponent = Sgl_hidden(src); 100 Sgl_clear_signexponent_set_hidden(src); 101 } 102 else { 103 /* normalize operand */ 104 Sgl_clear_signexponent(src); 105 src_exponent++; 106 Sgl_normalize(src,src_exponent); 107 even_exponent = src_exponent & 1; 108 } 109 if (even_exponent) { 110 /* exponent is even */ 111 /* Add comment here. Explain why odd exponent needs correction */ 112 Sgl_leftshiftby1(src); 113 } 114 /* 115 * Add comment here. Explain following algorithm. 116 * 117 * Trust me, it works. 118 * 119 */ 120 Sgl_setzero(result); 121 newbit = 1 << SGL_P; 122 while (newbit && Sgl_isnotzero(src)) { 123 Sgl_addition(result,newbit,sum); 124 if(sum <= Sgl_all(src)) { 125 /* update result */ 126 Sgl_addition(result,(newbit<<1),result); 127 Sgl_subtract(src,sum,src); 128 } 129 Sgl_rightshiftby1(newbit); 130 Sgl_leftshiftby1(src); 131 } 132 /* correct exponent for pre-shift */ 133 if (even_exponent) { 134 Sgl_rightshiftby1(result); 135 } 136 137 /* check for inexact */ 138 if (Sgl_isnotzero(src)) { 139 if (!even_exponent && Sgl_islessthan(result,src)) 140 Sgl_increment(result); 141 guardbit = Sgl_lowmantissa(result); 142 Sgl_rightshiftby1(result); 143 144 /* now round result */ 145 switch (Rounding_mode()) { 146 case ROUNDPLUS: 147 Sgl_increment(result); 148 break; 149 case ROUNDNEAREST: 150 /* stickybit is always true, so guardbit 151 * is enough to determine rounding */ 152 if (guardbit) { 153 Sgl_increment(result); 154 } 155 break; 156 } 157 /* increment result exponent by 1 if mantissa overflowed */ 158 if (Sgl_isone_hiddenoverflow(result)) src_exponent+=2; 159 160 if (Is_inexacttrap_enabled()) { 161 Sgl_set_exponent(result, 162 ((src_exponent-SGL_BIAS)>>1)+SGL_BIAS); 163 *dstptr = result; 164 return(INEXACTEXCEPTION); 165 } 166 else Set_inexactflag(); 167 } 168 else { 169 Sgl_rightshiftby1(result); 170 } 171 Sgl_set_exponent(result,((src_exponent-SGL_BIAS)>>1)+SGL_BIAS); 172 *dstptr = result; 173 return(NOEXCEPTION); 174 } 175