glibc/sysdeps/alpha/w_sqrt.S
1997-06-21 01:37:20 +00:00

162 lines
5.0 KiB
ArmAsm

/* Copyright (C) 1996 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by David Mosberger <davidm@cs.arizona.edu>, 1996.
Based on public-domain C source by Linus Torvalds.
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. */
/* This version is much faster than generic sqrt implementation, but
it doesn't handle exceptional values or the inexact flag. Don't use
this if _IEEE_FP or _IEEE_FP_INEXACT is in effect. */
#ifndef _IEEE_FP
#define _ERRNO_H
#include <bits/errno.h>
#include <sysdep.h>
.set noreorder
#ifdef __ELF__
.section .rodata
#else
.rdata
#endif
.align 5 # align to cache line
/* Do all memory accesses relative to sqrtdata. */
sqrtdata:
#define DN 0x00
#define UP 0x08
#define HALF 0x10
#define ALMOST_THREE_HALF 0x18
#define T2 0x20
.quad 0x3fefffffffffffff /* DN = next(1.0) */
.quad 0x3ff0000000000001 /* UP = prev(1.0) */
.quad 0x3fe0000000000000 /* HALF = 0.5 */
.quad 0x3ff7ffffffc00000 /* ALMOST_THREE_HALF = 1.5-2^-30 */
/* table T2: */
.long 0x1500, 0x2ef8, 0x4d67, 0x6b02, 0x87be, 0xa395, 0xbe7a, 0xd866
.long 0xf14a, 0x1091b, 0x11fcd, 0x13552, 0x14999, 0x15c98, 0x16e34, 0x17e5f
.long 0x18d03, 0x19a01, 0x1a545, 0x1ae8a, 0x1b5c4, 0x1bb01, 0x1bfde, 0x1c28d
.long 0x1c2de, 0x1c0db, 0x1ba73, 0x1b11c, 0x1a4b5, 0x1953d, 0x18266, 0x16be0
.long 0x1683e, 0x179d8, 0x18a4d, 0x19992, 0x1a789, 0x1b445, 0x1bf61, 0x1c989
.long 0x1d16d, 0x1d77b, 0x1dddf, 0x1e2ad, 0x1e5bf, 0x1e6e8, 0x1e654, 0x1e3cd
.long 0x1df2a, 0x1d635, 0x1cb16, 0x1be2c, 0x1ae4e, 0x19bde, 0x1868e, 0x16e2e
.long 0x1527f, 0x1334a, 0x11051, 0xe951, 0xbe01, 0x8e0d, 0x5924, 0x1edd
/*
* Stack variables:
*/
#define K 16(sp)
#define Y 24(sp)
#define FSIZE 32
.text
LEAF(__sqrt, FSIZE)
lda sp, -FSIZE(sp)
ldgp gp, .-__sqrt(pv)
stq ra, 0(sp)
#ifdef PROF
lda AT, _mcount
jsr AT, (AT), _mcount
#endif
.prologue 1
stt $f16, K
lda t3, sqrtdata # load base address into t3
fblt $f16, $negative
/* Compute initial guess. */
.align 3
ldah t1, 0x5fe8 # e0 :
ldq t2, K # .. e1 :
ldt $f12, HALF(t3) # e0 :
ldt $f18, ALMOST_THREE_HALF(t3) # .. e1 :
srl t2, 33, t0 # e0 :
mult $f16, $f12, $f11 # .. fm : $f11 = x * 0.5
subl t1, t0, t1 # e0 :
addt $f12, $f12, $f17 # .. fa : $f17 = 1.0
srl t1, 12, t0 # e0 :
and t0, 0xfc, t0 # .. e1 :
addq t0, t3, t0 # e0 :
ldl t0, T2(t0) # .. e1 :
addt $f12, $f17, $f15 # fa : $f15 = 1.5
subl t1, t0, t1 # .. e1 :
sll t1, 32, t1 # e0 :
ldt $f14, DN(t3) # .. e1 :
stq t1, Y # e0 :
ldt $f13, Y # e1 :
addq sp, FSIZE, sp # e0 :
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(t3) # .. 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
ret
$negative:
lda t1, -1
stq t1, K
lda t1, EDOM
stl t1, errno
#ifdef _LIBC_REENTRANT
jsr ra, __errno_location
lda t1, -1
ldq ra, 0(sp)
stl t1, 0(v0)
#endif
ldt $f0, K # res = (double) 0xffffffffffffffff
addq sp, FSIZE, sp
ret
END(__sqrt)
weak_alias(__sqrt, sqrt)
#endif /* !_IEEE_FP */