xref: /openbmc/linux/arch/parisc/math-emu/sfsqrt.c (revision 55fd7e02)
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