summaryrefslogtreecommitdiff
path: root/sysdeps/alpha/fpu/e_sqrt.c
diff options
context:
space:
mode:
authorUlrich Drepper <drepper@redhat.com>1997-09-05 01:43:20 +0000
committerUlrich Drepper <drepper@redhat.com>1997-09-05 01:43:20 +0000
commit26b4d7667169f8db26fd8194b3c498ec58e50f90 (patch)
tree491a6edd9fe8961fe781a50c39517dae8a35321d /sysdeps/alpha/fpu/e_sqrt.c
parent800d775e426b9c0af63f711b79b09bf540c97456 (diff)
1997-09-05 03:15 Ulrich Drepper <drepper@cygnus.com> * Makefile (subdirs): Reorganize order so that nss follows add-ons. * config.make.in (LDAP): Add. * configure.in: Prepare for ldap add-on. * shlib-versions: Add version numbers for ldap NSS module. * nss/Makefile: Add rules to build ldap NSS module. * db/btree/bt_page.c (__bt_free): Update meta data information. Update db from db-2.3.6. * db2/db.h: Updated. * db2/btree/bt_close.c: Updated. * db2/btree/bt_compare.c: Updated. * db2/btree/bt_conv.c: Updated. * db2/btree/bt_cursor.c: Updated. * db2/btree/bt_delete.c: Updated. * db2/btree/bt_put.c: Updated. * db2/btree/bt_rec.c: Updated. * db2/btree/bt_recno.c: Updated. * db2/btree/bt_search.c: Updated. * db2/btree/bt_split.c: Updated. * db2/btree/bt_stat.c: Updated. * db2/common/db_appinit.c: Updated. * db2/common/db_apprec.c: Updated. * db2/common/db_byteorder.c: Updated. * db2/common/db_err.c: Updated. * db2/common/db_region.c: Updated. * db2/db/db.c: Updated. * db2/db/db_conv.c: Updated. * db2/db/db_dup.c: Updated. * db2/db/db_pr.c: Updated. * db2/db/db_ret.c: Updated. * db2/db185/db185.c: Updated. * db2/dbm/dbm.c: Updated. * db2/hash/hash_dup.c: Updated. * db2/include/db_am.h: Updated. * db2/include/db_page.h: Updated. * db2/include/mp.h: Updated. * db2/include/queue.h: Updated. * db2/log/log.c: Updated. * db2/log/log_findckp.c: Updated. * db2/log/log_get.c: Updated. * db2/log/log_rec.c: Updated. * db2/mp/mp_bh.c: Updated. * db2/mp/mp_fopen.c: Updated. * db2/mp/mp_pr.c: Updated. * db2/mp/mp_sync.c: Updated. * db2/mutex/x86.gcc: Updated. * db2/os/db_os_dir.c: Updated. * db2/os/db_os_fid.c: Updated. * db2/progs/db_archive/db_archive.c: Updated. * db2/progs/db_checkpoint/db_checkpoint.c: Updated. * db2/progs/db_deadlock/db_deadlock.c: Updated. * db2/progs/db_dump/db_dump.c: Updated. * db2/progs/db_dump185/db_dump185.c: Updated. * db2/progs/db_load/db_load.c: Updated. * db2/progs/db_printlog/db_printlog.c: Updated. * db2/progs/db_recover/db_recover.c: Updated. * db2/progs/db_stat/db_stat.c: Updated. * db2/txn/txn.c: Updated. * db2/txn/txn_rec.c: Updated. * db2/Makefile: Add rule to build db_printlog. * db2/config.h: Don't include endian.h and don't define WORDS_BIGENDIAN. Only define HAVE_ENDIAN_H. * elf/dl-open.c (_dl_open): Make thread-safe. * elf/dl-close.c (_dl_close): Likewise. * elf/dlclose.c: Pretty print. * elf/dl-version.c (make_string): Use __strdup not strdup. (find_needed): Don't use l_searchlist of _dl_loaded, use the l_next list. * elf/dl-deps.c (_dl_map_object_deps): Remove variable head. * manual/math.texi: Add note about == and != for FP values. * manual/message.texi: Mention problems with relative paths and binstextdomain. * manual/string.texi: Mark strdupa and strndupa correctly as macros. * sunrpc/Makefile (CFLAGS-xkey_prot.c): Define to prevent warnings. * sunrpc/rpc_cout.c: Avoid needless races around return value. * sysdeps/generic/stpcpy.c: Include <config.h> if wanted. Update from tzdata-1997h. * time/africa: Updated. * time/europe: Updated. * time/iso3166.tab: Updated. * time/northamerica: Updated. * time/southamerica: Updated. * time/zone.tab: Updated. 1997-09-04 13:19 Richard Henderson <rth@cygnus.com> * sysdeps/alpha/w_sqrt.S: Removed. * sysdeps/alpha/fpu/e_sqrt.c: New. Obey -mieee and -mieee-with-inexact and build a version that is as fast as possible given the constraint. 1997-08-30 Andreas Jaeger <aj@arthur.rhein-neckar.de> * manual/maint.texi: Document some defaults for configuration parameters. 1997-09-04 15:57 Ulrich Drepper <drepper@cygnus.com> * termios/cfsetspeed.c (cfsetspeed): Change return value to int. * termios/termios.h: Change prototype accordingly. Reported by Andreas Jaeger <aj@arthur.rhein-neckar.de>. 1997-07-26 14:42 H.J. Lu <hjl@gnu.ai.mit.edu> * libio/stdfiles.c (_libio_using_thunks): Define if _G_USING_THUNKS is defined. * libio/_G_config.h: Moved to ... * sysdeps/generic/_G_config.h: ...here. * sysdeps/unix/sysv/linux/_G_config.h: New, Add vtable thunks support. 1997-09-02 Andreas Schwab <schwab@issan.informatik.uni-dortmund.de> * sysdeps/m68k/Makefile (pic-ccflag): Override it for more efficient code (together with a change in gcc). * Makeconfig (pic-ccflag): New variable. (CFLAGS-.os): Use it instead of hardcoding -fPIC. 1997-09-02 Andreas Schwab <schwab@issan.informatik.uni-dortmund.de> * sysdeps/m68k/dl-machine.h (elf_machine_load_address): Use _dl_start instead of a local label, to avoid getting an unreadable label name. 1997-09-02 Andreas Schwab <schwab@issan.informatik.uni-dortmund.de> * nss/Makefile ($(objpfx)libnss_db.so): Use $(dbobjdir) to find libdb.so. * Makeconfig (dbobjdir): Use db2. 1997-08-31 Andreas Schwab <schwab@issan.informatik.uni-dortmund.de> * manual/maint.texi (Porting to Unix): Update description of the implementation of Unix system calls. (Contributors): Use real umlaut in tege's name. 1997-08-30 Andreas Schwab <schwab@issan.informatik.uni-dortmund.de> * sysdeps/libm-ieee754/w_gamma.c: Make compilable with _IEEE_LIBM. * sysdeps/libm-ieee754/w_gammaf.c: Likewise. * sysdeps/libm-ieee754/w_gammal.c: Likewise.
Diffstat (limited to 'sysdeps/alpha/fpu/e_sqrt.c')
-rw-r--r--sysdeps/alpha/fpu/e_sqrt.c256
1 files changed, 256 insertions, 0 deletions
diff --git a/sysdeps/alpha/fpu/e_sqrt.c b/sysdeps/alpha/fpu/e_sqrt.c
new file mode 100644
index 0000000000..76fa015622
--- /dev/null
+++ b/sysdeps/alpha/fpu/e_sqrt.c
@@ -0,0 +1,256 @@
+/* Copyright (C) 1996, 1997 Free Software Foundation, Inc.
+ Contributed by David Mosberger (davidm@cs.arizona.edu).
+
+ This file is part of the GNU C Library.
+
+ The GNU C Library is free software; you can redistribute it and/or
+ modify it under the terms of the GNU Library General Public License as
+ published by the Free Software Foundation; either version 2 of the
+ License, or (at your option) any later version.
+
+ The GNU C Library 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
+ Library General Public License for more details.
+
+ You should have received a copy of the GNU Library General Public
+ License along with the GNU C Library; see the file COPYING.LIB. If not,
+ write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+ Boston, MA 02111-1307, USA. */
+
+/*
+ * We have three versions, depending on how exact we need the results.
+ */
+
+#if defined(_IEEE_FP) && defined(_IEEE_FP_INEXACT)
+
+/* Most demanding: go to the original source. */
+#include <libm-ieee754/e_sqrt.c>
+
+#else
+
+/* Careful with rearranging this without consulting the assembly below. */
+const static struct sqrt_data_struct {
+ unsigned long dn, up, half, almost_three_half;
+ unsigned long one_and_a_half, two_to_minus_30, one, nan;
+ const int T2[64];
+} sqrt_data = {
+ 0x3fefffffffffffff, /* __dn = nextafter(1,-Inf) */
+ 0x3ff0000000000001, /* __up = nextafter(1,+Inf) */
+ 0x3fe0000000000000, /* half */
+ 0x3ff7ffffffc00000, /* almost_three_half = 1.5-2^-30 */
+ 0x3ff8000000000000, /* one_and_a_half */
+ 0x3e10000000000000, /* two_to_minus_30 */
+ 0x3ff0000000000000, /* one */
+ 0xffffffffffffffff, /* nan */
+
+ { 0x1500, 0x2ef8, 0x4d67, 0x6b02, 0x87be, 0xa395, 0xbe7a, 0xd866,
+ 0xf14a, 0x1091b,0x11fcd,0x13552,0x14999,0x15c98,0x16e34,0x17e5f,
+ 0x18d03,0x19a01,0x1a545,0x1ae8a,0x1b5c4,0x1bb01,0x1bfde,0x1c28d,
+ 0x1c2de,0x1c0db,0x1ba73,0x1b11c,0x1a4b5,0x1953d,0x18266,0x16be0,
+ 0x1683e,0x179d8,0x18a4d,0x19992,0x1a789,0x1b445,0x1bf61,0x1c989,
+ 0x1d16d,0x1d77b,0x1dddf,0x1e2ad,0x1e5bf,0x1e6e8,0x1e654,0x1e3cd,
+ 0x1df2a,0x1d635,0x1cb16,0x1be2c,0x1ae4e,0x19bde,0x1868e,0x16e2e,
+ 0x1527f,0x1334a,0x11051,0xe951, 0xbe01, 0x8e0d, 0x5924, 0x1edd }
+};
+
+#ifdef _IEEE_FP
+/*
+ * This version is much faster than the standard one included above,
+ * but it doesn't maintain the inexact flag.
+ */
+
+#define lobits(x) (((unsigned int *)&x)[0])
+#define hibits(x) (((unsigned int *)&x)[1])
+
+static inline double initial_guess(double x, unsigned int k,
+ const struct sqrt_data_struct * const ptr)
+{
+ double ret = 0.0;
+
+ k = 0x5fe80000 - (k >> 1);
+ k = k - ptr->T2[63&(k>>14)];
+ hibits(ret) = k;
+ return ret;
+}
+
+/* up = nextafter(1,+Inf), dn = nextafter(1,-Inf) */
+
+#define __half (ptr->half)
+#define __one_and_a_half (ptr->one_and_a_half)
+#define __two_to_minus_30 (ptr->two_to_minus_30)
+#define __one (ptr->one)
+#define __up (ptr->up)
+#define __dn (ptr->dn)
+#define __Nan (ptr->nan)
+
+#define Double(x) (*(double *)&x)
+
+/* Multiply with chopping rounding.. */
+#define choppedmul(a,b,c) \
+ __asm__("multc %1,%2,%0":"=&f" (c):"f" (a), "f" (b))
+
+double
+__ieee754_sqrt(double x)
+{
+ const struct sqrt_data_struct * const ptr = &sqrt_data;
+ unsigned long k, bits;
+ double y, z, zp, zn;
+ double dn, up, low, high;
+ double half, one_and_a_half, one, two_to_minus_30;
+
+ *(double *)&bits = x;
+ k = bits;
+
+ /* Negative or NaN or Inf */
+ if ((k >> 52) >= 0x7ff)
+ goto special;
+ y = initial_guess(x, k >> 32, ptr);
+ half = Double(__half);
+ one_and_a_half = Double(__one_and_a_half);
+ y = y*(one_and_a_half - half*x*y*y);
+ dn = Double(__dn);
+ two_to_minus_30 = Double(__two_to_minus_30);
+ y = y*((one_and_a_half - two_to_minus_30) - half*x*y*y);
+ up = Double(__up);
+ z = x*y;
+ one = Double(__one);
+ z = z + half*z*(one-z*y);
+
+ choppedmul(z,dn,zp);
+ choppedmul(z,up,zn);
+
+ choppedmul(z,zp,low);
+ low = low - x;
+ choppedmul(z,zn,high);
+ high = high - x;
+
+ /* I can't get gcc to use fcmov's.. */
+ __asm__("fcmovge %2,%3,%0"
+ :"=f" (z)
+ :"0" (z), "f" (low), "f" (zp));
+ __asm__("fcmovlt %2,%3,%0"
+ :"=f" (z)
+ :"0" (z), "f" (high), "f" (zn));
+ return z; /* Argh! gcc jumps to end here */
+
+special:
+ /* throw away sign bit */
+ k <<= 1;
+ /* -0 */
+ if (!k)
+ return x;
+ /* special? */
+ if ((k >> 53) == 0x7ff) {
+ /* NaN? */
+ if (k << 11)
+ return x;
+ /* sqrt(+Inf) = +Inf */
+ if (x > 0)
+ return x;
+ }
+
+ x = Double(__Nan);
+ return x;
+}
+
+#else
+/*
+ * This version is much faster than generic sqrt implementation, but
+ * it doesn't handle exceptional values or the inexact flag.
+ */
+
+asm ("\
+ /* Define offsets into the structure defined in C above. */
+ $DN = 0*8
+ $UP = 1*8
+ $HALF = 2*8
+ $ALMOST_THREE_HALF = 3*8
+ $NAN = 7*8
+ $T2 = 8*8
+
+ /* Stack variables. */
+ $K = 0
+ $Y = 8
+
+ .text
+ .align 3
+ .globl __ieee754_sqrt
+ .ent __ieee754_sqrt
+__ieee754_sqrt:
+ ldgp $29, 0($27)
+ subq $sp, 16, $sp
+ .frame $sp, 16, $26, 0\n"
+#ifdef PROF
+" lda $28, _mcount
+ jsr $28, ($28), _mcount\n"
+#endif
+" .prologue 1
+
+ stt $f16, $K($sp)
+ lda $4, sqrt_data # load base address into t3
+ fblt $f16, $negative
+
+ /* Compute initial guess. */
+
+ .align 3
+
+ ldah $2, 0x5fe8 # e0 :
+ ldq $3, $K($sp) # .. e1 :
+ ldt $f12, $HALF($4) # e0 :
+ ldt $f18, $ALMOST_THREE_HALF($4) # .. e1 :
+ srl $3, 33, $1 # e0 :
+ mult $f16, $f12, $f11 # .. fm : $f11 = x * 0.5
+ subl $2, $1, $2 # e0 :
+ addt $f12, $f12, $f17 # .. fa : $f17 = 1.0
+ srl $2, 12, $1 # e0 :
+ and $1, 0xfc, $1 # .. e1 :
+ addq $1, $4, $1 # e0 :
+ ldl $1, $T2($1) # .. e1 :
+ addt $f12, $f17, $f15 # fa : $f15 = 1.5
+ subl $2, $1, $2 # .. e1 :
+ sll $2, 32, $2 # e0 :
+ ldt $f14, $DN($4) # .. e1 :
+ stq $2, $Y($sp) # e0 :
+ ldt $f13, $Y($sp) # e1 :
+
+ mult $f11, $f13, $f10 # fm : $f10 = (x * 0.5) * y
+ mult $f10, $f13, $f10 # fm : $f10 = ((x * 0.5) * y) * y
+ subt $f15, $f10, $f1 # fa : $f1 = (1.5 - 0.5*x*y*y)
+ mult $f13, $f1, $f13 # fm : yp = y*(1.5 - 0.5*x*y*y)
+ mult $f11, $f13, $f11 # fm : $f11 = x * 0.5 * yp
+ mult $f11, $f13, $f11 # fm : $f11 = (x * 0.5 * yp) * yp
+ subt $f18, $f11, $f1 # fa : $f1= (1.5-2^-30) - 0.5*x*yp*yp
+ mult $f13, $f1, $f13 # fm : ypp = $f13 = yp*$f1
+ subt $f15, $f12, $f1 # fa : $f1 = (1.5 - 0.5)
+ ldt $f15, $UP($4) # .. e1 :
+ mult $f16, $f13, $f10 # fm : z = $f10 = x * ypp
+ mult $f10, $f13, $f11 # fm : $f11 = z*ypp
+ mult $f10, $f12, $f12 # fm : $f12 = z*0.5
+ subt $f1, $f11, $f1 # .. fa : $f1 = 1 - z*ypp
+ mult $f12, $f1, $f12 # fm : $f12 = z*0.5*(1 - z*ypp)
+ addt $f10, $f12, $f0 # fa : zp=res=$f0= z + z*0.5*(1 - z*ypp)
+
+ mult/c $f0, $f14, $f12 # fm : zmi = zp * DN
+ mult/c $f0, $f15, $f11 # fm : zpl = zp * UP
+ mult/c $f0, $f12, $f1 # fm : $f1 = zp * zmi
+ mult/c $f0, $f11, $f15 # fm : $f15 = zp * zpl
+
+ subt $f1, $f16, $f13 # fa : y1 = zp*zmi - x
+ subt $f15, $f16, $f15 # fa : y2 = zp*zpl - x
+
+ fcmovge $f13, $f12, $f0 # res = (y1 >= 0) ? zmi : res
+ fcmovlt $f15, $f11, $f0 # res = (y2 < 0) ? zpl : res
+
+ addq $sp, 16, $sp # e0 :
+ ret # .. e1 :
+
+$negative:
+ ldt $f0, $NAN($4)
+ addq $sp, 16, $sp
+ ret
+
+ .end __ieee754_sqrt");
+
+#endif /* _IEEE_FP */
+#endif /* _IEEE_FP && _IEEE_FP_INEXACT */