summaryrefslogtreecommitdiff
path: root/arch/parisc/math-emu/dfsqrt.c
diff options
context:
space:
mode:
authorAndré Fabian Silva Delgado <emulatorman@parabola.nu>2015-08-05 17:04:01 -0300
committerAndré Fabian Silva Delgado <emulatorman@parabola.nu>2015-08-05 17:04:01 -0300
commit57f0f512b273f60d52568b8c6b77e17f5636edc0 (patch)
tree5e910f0e82173f4ef4f51111366a3f1299037a7b /arch/parisc/math-emu/dfsqrt.c
Initial import
Diffstat (limited to 'arch/parisc/math-emu/dfsqrt.c')
-rw-r--r--arch/parisc/math-emu/dfsqrt.c195
1 files changed, 195 insertions, 0 deletions
diff --git a/arch/parisc/math-emu/dfsqrt.c b/arch/parisc/math-emu/dfsqrt.c
new file mode 100644
index 000000000..9542c6d28
--- /dev/null
+++ b/arch/parisc/math-emu/dfsqrt.c
@@ -0,0 +1,195 @@
+/*
+ * Linux/PA-RISC Project (http://www.parisc-linux.org/)
+ *
+ * Floating-point emulation code
+ * Copyright (C) 2001 Hewlett-Packard (Paul Bame) <bame@debian.org>
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2, or (at your option)
+ * any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ */
+/*
+ * BEGIN_DESC
+ *
+ * File:
+ * @(#) pa/spmath/dfsqrt.c $Revision: 1.1 $
+ *
+ * Purpose:
+ * Double Floating-point Square Root
+ *
+ * External Interfaces:
+ * dbl_fsqrt(srcptr,nullptr,dstptr,status)
+ *
+ * Internal Interfaces:
+ *
+ * Theory:
+ * <<please update with a overview of the operation of this file>>
+ *
+ * END_DESC
+*/
+
+
+#include "float.h"
+#include "dbl_float.h"
+
+/*
+ * Double Floating-point Square Root
+ */
+
+/*ARGSUSED*/
+unsigned int
+dbl_fsqrt(
+ dbl_floating_point *srcptr,
+ unsigned int *nullptr,
+ dbl_floating_point *dstptr,
+ unsigned int *status)
+{
+ register unsigned int srcp1, srcp2, resultp1, resultp2;
+ register unsigned int newbitp1, newbitp2, sump1, sump2;
+ register int src_exponent;
+ register boolean guardbit = FALSE, even_exponent;
+
+ Dbl_copyfromptr(srcptr,srcp1,srcp2);
+ /*
+ * check source operand for NaN or infinity
+ */
+ if ((src_exponent = Dbl_exponent(srcp1)) == DBL_INFINITY_EXPONENT) {
+ /*
+ * is signaling NaN?
+ */
+ if (Dbl_isone_signaling(srcp1)) {
+ /* trap if INVALIDTRAP enabled */
+ if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
+ /* make NaN quiet */
+ Set_invalidflag();
+ Dbl_set_quiet(srcp1);
+ }
+ /*
+ * Return quiet NaN or positive infinity.
+ * Fall through to negative test if negative infinity.
+ */
+ if (Dbl_iszero_sign(srcp1) ||
+ Dbl_isnotzero_mantissa(srcp1,srcp2)) {
+ Dbl_copytoptr(srcp1,srcp2,dstptr);
+ return(NOEXCEPTION);
+ }
+ }
+
+ /*
+ * check for zero source operand
+ */
+ if (Dbl_iszero_exponentmantissa(srcp1,srcp2)) {
+ Dbl_copytoptr(srcp1,srcp2,dstptr);
+ return(NOEXCEPTION);
+ }
+
+ /*
+ * check for negative source operand
+ */
+ if (Dbl_isone_sign(srcp1)) {
+ /* trap if INVALIDTRAP enabled */
+ if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
+ /* make NaN quiet */
+ Set_invalidflag();
+ Dbl_makequietnan(srcp1,srcp2);
+ Dbl_copytoptr(srcp1,srcp2,dstptr);
+ return(NOEXCEPTION);
+ }
+
+ /*
+ * Generate result
+ */
+ if (src_exponent > 0) {
+ even_exponent = Dbl_hidden(srcp1);
+ Dbl_clear_signexponent_set_hidden(srcp1);
+ }
+ else {
+ /* normalize operand */
+ Dbl_clear_signexponent(srcp1);
+ src_exponent++;
+ Dbl_normalize(srcp1,srcp2,src_exponent);
+ even_exponent = src_exponent & 1;
+ }
+ if (even_exponent) {
+ /* exponent is even */
+ /* Add comment here. Explain why odd exponent needs correction */
+ Dbl_leftshiftby1(srcp1,srcp2);
+ }
+ /*
+ * Add comment here. Explain following algorithm.
+ *
+ * Trust me, it works.
+ *
+ */
+ Dbl_setzero(resultp1,resultp2);
+ Dbl_allp1(newbitp1) = 1 << (DBL_P - 32);
+ Dbl_setzero_mantissap2(newbitp2);
+ while (Dbl_isnotzero(newbitp1,newbitp2) && Dbl_isnotzero(srcp1,srcp2)) {
+ Dbl_addition(resultp1,resultp2,newbitp1,newbitp2,sump1,sump2);
+ if(Dbl_isnotgreaterthan(sump1,sump2,srcp1,srcp2)) {
+ Dbl_leftshiftby1(newbitp1,newbitp2);
+ /* update result */
+ Dbl_addition(resultp1,resultp2,newbitp1,newbitp2,
+ resultp1,resultp2);
+ Dbl_subtract(srcp1,srcp2,sump1,sump2,srcp1,srcp2);
+ Dbl_rightshiftby2(newbitp1,newbitp2);
+ }
+ else {
+ Dbl_rightshiftby1(newbitp1,newbitp2);
+ }
+ Dbl_leftshiftby1(srcp1,srcp2);
+ }
+ /* correct exponent for pre-shift */
+ if (even_exponent) {
+ Dbl_rightshiftby1(resultp1,resultp2);
+ }
+
+ /* check for inexact */
+ if (Dbl_isnotzero(srcp1,srcp2)) {
+ if (!even_exponent && Dbl_islessthan(resultp1,resultp2,srcp1,srcp2)) {
+ Dbl_increment(resultp1,resultp2);
+ }
+ guardbit = Dbl_lowmantissap2(resultp2);
+ Dbl_rightshiftby1(resultp1,resultp2);
+
+ /* now round result */
+ switch (Rounding_mode()) {
+ case ROUNDPLUS:
+ Dbl_increment(resultp1,resultp2);
+ break;
+ case ROUNDNEAREST:
+ /* stickybit is always true, so guardbit
+ * is enough to determine rounding */
+ if (guardbit) {
+ Dbl_increment(resultp1,resultp2);
+ }
+ break;
+ }
+ /* increment result exponent by 1 if mantissa overflowed */
+ if (Dbl_isone_hiddenoverflow(resultp1)) src_exponent+=2;
+
+ if (Is_inexacttrap_enabled()) {
+ Dbl_set_exponent(resultp1,
+ ((src_exponent-DBL_BIAS)>>1)+DBL_BIAS);
+ Dbl_copytoptr(resultp1,resultp2,dstptr);
+ return(INEXACTEXCEPTION);
+ }
+ else Set_inexactflag();
+ }
+ else {
+ Dbl_rightshiftby1(resultp1,resultp2);
+ }
+ Dbl_set_exponent(resultp1,((src_exponent-DBL_BIAS)>>1)+DBL_BIAS);
+ Dbl_copytoptr(resultp1,resultp2,dstptr);
+ return(NOEXCEPTION);
+}