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