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*/
37unsigned int
38sgl_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