summaryrefslogtreecommitdiff
path: root/arch/parisc/math-emu/sfsqrt.c
diff options
context:
space:
mode:
authorLinus Torvalds <torvalds@ppc970.osdl.org>2005-04-16 15:20:36 -0700
committerLinus Torvalds <torvalds@ppc970.osdl.org>2005-04-16 15:20:36 -0700
commit1da177e4c3f41524e886b7f1b8a0c1fc7321cac2 (patch)
tree0bba044c4ce775e45a88a51686b5d9f90697ea9d /arch/parisc/math-emu/sfsqrt.c
Linux-2.6.12-rc2
Initial git repository build. I'm not bothering with the full history, even though we have it. We can create a separate "historical" git archive of that later if we want to, and in the meantime it's about 3.2GB when imported into git - space that would just make the early git days unnecessarily complicated, when we don't have a lot of good infrastructure for it. Let it rip!
Diffstat (limited to 'arch/parisc/math-emu/sfsqrt.c')
-rw-r--r--arch/parisc/math-emu/sfsqrt.c187
1 files changed, 187 insertions, 0 deletions
diff --git a/arch/parisc/math-emu/sfsqrt.c b/arch/parisc/math-emu/sfsqrt.c
new file mode 100644
index 000000000000..cd3f6db1f6fb
--- /dev/null
+++ b/arch/parisc/math-emu/sfsqrt.c
@@ -0,0 +1,187 @@
+/*
+ * 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/sfsqrt.c $Revision: 1.1 $
+ *
+ * Purpose:
+ * Single Floating-point Square Root
+ *
+ * External Interfaces:
+ * sgl_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 "sgl_float.h"
+
+/*
+ * Single Floating-point Square Root
+ */
+
+/*ARGSUSED*/
+unsigned int
+sgl_fsqrt(
+ sgl_floating_point *srcptr,
+ unsigned int *nullptr,
+ sgl_floating_point *dstptr,
+ unsigned int *status)
+{
+ register unsigned int src, result;
+ register int src_exponent;
+ register unsigned int newbit, sum;
+ register boolean guardbit = FALSE, even_exponent;
+
+ src = *srcptr;
+ /*
+ * check source operand for NaN or infinity
+ */
+ if ((src_exponent = Sgl_exponent(src)) == SGL_INFINITY_EXPONENT) {
+ /*
+ * is signaling NaN?
+ */
+ if (Sgl_isone_signaling(src)) {
+ /* trap if INVALIDTRAP enabled */
+ if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
+ /* make NaN quiet */
+ Set_invalidflag();
+ Sgl_set_quiet(src);
+ }
+ /*
+ * Return quiet NaN or positive infinity.
+ * Fall thru to negative test if negative infinity.
+ */
+ if (Sgl_iszero_sign(src) || Sgl_isnotzero_mantissa(src)) {
+ *dstptr = src;
+ return(NOEXCEPTION);
+ }
+ }
+
+ /*
+ * check for zero source operand
+ */
+ if (Sgl_iszero_exponentmantissa(src)) {
+ *dstptr = src;
+ return(NOEXCEPTION);
+ }
+
+ /*
+ * check for negative source operand
+ */
+ if (Sgl_isone_sign(src)) {
+ /* trap if INVALIDTRAP enabled */
+ if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
+ /* make NaN quiet */
+ Set_invalidflag();
+ Sgl_makequietnan(src);
+ *dstptr = src;
+ return(NOEXCEPTION);
+ }
+
+ /*
+ * Generate result
+ */
+ if (src_exponent > 0) {
+ even_exponent = Sgl_hidden(src);
+ Sgl_clear_signexponent_set_hidden(src);
+ }
+ else {
+ /* normalize operand */
+ Sgl_clear_signexponent(src);
+ src_exponent++;
+ Sgl_normalize(src,src_exponent);
+ even_exponent = src_exponent & 1;
+ }
+ if (even_exponent) {
+ /* exponent is even */
+ /* Add comment here. Explain why odd exponent needs correction */
+ Sgl_leftshiftby1(src);
+ }
+ /*
+ * Add comment here. Explain following algorithm.
+ *
+ * Trust me, it works.
+ *
+ */
+ Sgl_setzero(result);
+ newbit = 1 << SGL_P;
+ while (newbit && Sgl_isnotzero(src)) {
+ Sgl_addition(result,newbit,sum);
+ if(sum <= Sgl_all(src)) {
+ /* update result */
+ Sgl_addition(result,(newbit<<1),result);
+ Sgl_subtract(src,sum,src);
+ }
+ Sgl_rightshiftby1(newbit);
+ Sgl_leftshiftby1(src);
+ }
+ /* correct exponent for pre-shift */
+ if (even_exponent) {
+ Sgl_rightshiftby1(result);
+ }
+
+ /* check for inexact */
+ if (Sgl_isnotzero(src)) {
+ if (!even_exponent && Sgl_islessthan(result,src))
+ Sgl_increment(result);
+ guardbit = Sgl_lowmantissa(result);
+ Sgl_rightshiftby1(result);
+
+ /* now round result */
+ switch (Rounding_mode()) {
+ case ROUNDPLUS:
+ Sgl_increment(result);
+ break;
+ case ROUNDNEAREST:
+ /* stickybit is always true, so guardbit
+ * is enough to determine rounding */
+ if (guardbit) {
+ Sgl_increment(result);
+ }
+ break;
+ }
+ /* increment result exponent by 1 if mantissa overflowed */
+ if (Sgl_isone_hiddenoverflow(result)) src_exponent+=2;
+
+ if (Is_inexacttrap_enabled()) {
+ Sgl_set_exponent(result,
+ ((src_exponent-SGL_BIAS)>>1)+SGL_BIAS);
+ *dstptr = result;
+ return(INEXACTEXCEPTION);
+ }
+ else Set_inexactflag();
+ }
+ else {
+ Sgl_rightshiftby1(result);
+ }
+ Sgl_set_exponent(result,((src_exponent-SGL_BIAS)>>1)+SGL_BIAS);
+ *dstptr = result;
+ return(NOEXCEPTION);
+}