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