cachepc-linux

Fork of AMDESE/linux with modifications for CachePC side-channel attack
git clone https://git.sinitax.com/sinitax/cachepc-linux
Log | Files | Refs | README | LICENSE | sfeed.txt

sfsqrt.c (4163B)

      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}