^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/dfsqrt.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) * Double 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) * dbl_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 "dbl_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) * Double 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) dbl_fsqrt(
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 39) dbl_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) dbl_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 srcp1, srcp2, resultp1, resultp2;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 45) register unsigned int newbitp1, newbitp2, sump1, sump2;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 46) register int src_exponent;
^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) Dbl_copyfromptr(srcptr,srcp1,srcp2);
^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 = Dbl_exponent(srcp1)) == DBL_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 (Dbl_isone_signaling(srcp1)) {
^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) Dbl_set_quiet(srcp1);
^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 (Dbl_iszero_sign(srcp1) ||
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 69) Dbl_isnotzero_mantissa(srcp1,srcp2)) {
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 70) Dbl_copytoptr(srcp1,srcp2,dstptr);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 71) return(NOEXCEPTION);
^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) /*
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 76) * check for zero source operand
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 77) */
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 78) if (Dbl_iszero_exponentmantissa(srcp1,srcp2)) {
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 79) Dbl_copytoptr(srcp1,srcp2,dstptr);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 80) return(NOEXCEPTION);
^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) /*
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 84) * check for negative source operand
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 85) */
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 86) if (Dbl_isone_sign(srcp1)) {
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 87) /* trap if INVALIDTRAP enabled */
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 88) if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 89) /* make NaN quiet */
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 90) Set_invalidflag();
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 91) Dbl_makequietnan(srcp1,srcp2);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 92) Dbl_copytoptr(srcp1,srcp2,dstptr);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 93) return(NOEXCEPTION);
^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) /*
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 97) * Generate result
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 98) */
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 99) if (src_exponent > 0) {
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 100) even_exponent = Dbl_hidden(srcp1);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 101) Dbl_clear_signexponent_set_hidden(srcp1);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 102) }
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 103) else {
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 104) /* normalize operand */
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 105) Dbl_clear_signexponent(srcp1);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 106) src_exponent++;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 107) Dbl_normalize(srcp1,srcp2,src_exponent);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 108) even_exponent = src_exponent & 1;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 109) }
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 110) if (even_exponent) {
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 111) /* exponent is even */
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 112) /* Add comment here. Explain why odd exponent needs correction */
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 113) Dbl_leftshiftby1(srcp1,srcp2);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 114) }
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 115) /*
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 116) * Add comment here. Explain following algorithm.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 117) *
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 118) * Trust me, it works.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 119) *
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 120) */
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 121) Dbl_setzero(resultp1,resultp2);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 122) Dbl_allp1(newbitp1) = 1 << (DBL_P - 32);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 123) Dbl_setzero_mantissap2(newbitp2);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 124) while (Dbl_isnotzero(newbitp1,newbitp2) && Dbl_isnotzero(srcp1,srcp2)) {
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 125) Dbl_addition(resultp1,resultp2,newbitp1,newbitp2,sump1,sump2);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 126) if(Dbl_isnotgreaterthan(sump1,sump2,srcp1,srcp2)) {
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 127) Dbl_leftshiftby1(newbitp1,newbitp2);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 128) /* update result */
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 129) Dbl_addition(resultp1,resultp2,newbitp1,newbitp2,
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 130) resultp1,resultp2);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 131) Dbl_subtract(srcp1,srcp2,sump1,sump2,srcp1,srcp2);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 132) Dbl_rightshiftby2(newbitp1,newbitp2);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 133) }
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 134) else {
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 135) Dbl_rightshiftby1(newbitp1,newbitp2);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 136) }
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 137) Dbl_leftshiftby1(srcp1,srcp2);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 138) }
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 139) /* correct exponent for pre-shift */
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 140) if (even_exponent) {
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 141) Dbl_rightshiftby1(resultp1,resultp2);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 142) }
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 143)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 144) /* check for inexact */
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 145) if (Dbl_isnotzero(srcp1,srcp2)) {
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 146) if (!even_exponent && Dbl_islessthan(resultp1,resultp2,srcp1,srcp2)) {
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 147) Dbl_increment(resultp1,resultp2);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 148) }
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 149) guardbit = Dbl_lowmantissap2(resultp2);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 150) Dbl_rightshiftby1(resultp1,resultp2);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 151)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 152) /* now round result */
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 153) switch (Rounding_mode()) {
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 154) case ROUNDPLUS:
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 155) Dbl_increment(resultp1,resultp2);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 156) break;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 157) case ROUNDNEAREST:
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 158) /* stickybit is always true, so guardbit
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 159) * is enough to determine rounding */
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 160) if (guardbit) {
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 161) Dbl_increment(resultp1,resultp2);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 162) }
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 163) break;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 164) }
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 165) /* increment result exponent by 1 if mantissa overflowed */
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 166) if (Dbl_isone_hiddenoverflow(resultp1)) src_exponent+=2;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 167)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 168) if (Is_inexacttrap_enabled()) {
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 169) Dbl_set_exponent(resultp1,
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 170) ((src_exponent-DBL_BIAS)>>1)+DBL_BIAS);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 171) Dbl_copytoptr(resultp1,resultp2,dstptr);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 172) return(INEXACTEXCEPTION);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 173) }
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 174) else Set_inexactflag();
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 175) }
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 176) else {
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 177) Dbl_rightshiftby1(resultp1,resultp2);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 178) }
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 179) Dbl_set_exponent(resultp1,((src_exponent-DBL_BIAS)>>1)+DBL_BIAS);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 180) Dbl_copytoptr(resultp1,resultp2,dstptr);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 181) return(NOEXCEPTION);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 182) }
Orange Pi5 kernel
Deprecated Linux kernel 5.10.110 for OrangePi 5/5B/5+ boards
3 Commits
0 Branches
0 Tags