sfsqrt.c

  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}