mirror of git://sourceware.org/git/glibc.git
* sysdeps/alpha/fpu/e_sqrt.c: Use the asm version when the input is
a finite non-denormal, deferring to the full IEEE version otherwise.
This commit is contained in:
Richard Henderson
parent
d0c425dbc5
commit
9b1370b857
|
|
@ -1,4 +1,4 @@
|
||||||
/* Copyright (C) 1996, 1997 Free Software Foundation, Inc.
|
/* Copyright (C) 1996, 1997, 1998 Free Software Foundation, Inc.
|
||||||
Contributed by David Mosberger (davidm@cs.arizona.edu).
|
Contributed by David Mosberger (davidm@cs.arizona.edu).
|
||||||
|
|
||||||
This file is part of the GNU C Library.
|
This file is part of the GNU C Library.
|
||||||
|
|
@ -18,17 +18,16 @@
|
||||||
write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
|
write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
|
||||||
Boston, MA 02111-1307, USA. */
|
Boston, MA 02111-1307, USA. */
|
||||||
|
|
||||||
|
|
||||||
|
#if !defined(_IEEE_FP_INEXACT)
|
||||||
|
|
||||||
/*
|
/*
|
||||||
* We have three versions, depending on how exact we need the results.
|
* This version is much faster than generic sqrt implementation, but
|
||||||
|
* it doesn't handle the inexact flag. It doesn't handle exceptional
|
||||||
|
* values either, but will defer to the full ieee754_sqrt routine which
|
||||||
|
* can.
|
||||||
*/
|
*/
|
||||||
|
|
||||||
#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. */
|
/* Careful with rearranging this without consulting the assembly below. */
|
||||||
const static struct sqrt_data_struct {
|
const static struct sqrt_data_struct {
|
||||||
unsigned long dn, up, half, almost_three_half;
|
unsigned long dn, up, half, almost_three_half;
|
||||||
|
|
@ -54,112 +53,6 @@ const static struct sqrt_data_struct {
|
||||||
0x1527f,0x1334a,0x11051,0xe951, 0xbe01, 0x8e0d, 0x5924, 0x1edd }
|
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 ("\
|
asm ("\
|
||||||
/* Define offsets into the structure defined in C above. */
|
/* Define offsets into the structure defined in C above. */
|
||||||
$DN = 0*8
|
$DN = 0*8
|
||||||
|
|
@ -174,7 +67,7 @@ asm ("\
|
||||||
$Y = 8
|
$Y = 8
|
||||||
|
|
||||||
.text
|
.text
|
||||||
.align 3
|
.align 5
|
||||||
.globl __ieee754_sqrt
|
.globl __ieee754_sqrt
|
||||||
.ent __ieee754_sqrt
|
.ent __ieee754_sqrt
|
||||||
__ieee754_sqrt:
|
__ieee754_sqrt:
|
||||||
|
|
@ -187,72 +80,86 @@ __ieee754_sqrt:
|
||||||
#endif
|
#endif
|
||||||
" .prologue 1
|
" .prologue 1
|
||||||
|
|
||||||
stt $f16, $K($sp)
|
.align 4
|
||||||
lda $4, sqrt_data # load base address into t3
|
stt $f16, $K($sp) # e0 :
|
||||||
fblt $f16, $negative
|
mult $f31, $f31, $f31 # .. fm :
|
||||||
|
lda $4, sqrt_data # e0 :
|
||||||
/* Compute initial guess. */
|
fblt $f16, $fixup # .. fa :
|
||||||
|
|
||||||
.align 3
|
|
||||||
|
|
||||||
ldah $2, 0x5fe8 # e0 :
|
ldah $2, 0x5fe8 # e0 :
|
||||||
ldq $3, $K($sp) # .. e1 :
|
ldq $3, $K($sp) # .. e1 :
|
||||||
ldt $f12, $HALF($4) # e0 :
|
ldt $f12, $HALF($4) # e0 :
|
||||||
ldt $f18, $ALMOST_THREE_HALF($4) # .. e1 :
|
ldt $f18, $ALMOST_THREE_HALF($4) # .. e1 :
|
||||||
srl $3, 33, $1 # e0 :
|
|
||||||
mult $f16, $f12, $f11 # .. fm : $f11 = x * 0.5
|
sll $3, 52, $5 # e0 :
|
||||||
subl $2, $1, $2 # e0 :
|
lda $6, 0x7fd # .. e1 :
|
||||||
|
fnop # .. fa :
|
||||||
|
fnop # .. fm :
|
||||||
|
|
||||||
|
subq $5, 1, $5 # e1 :
|
||||||
|
srl $3, 33, $1 # .. e0 :
|
||||||
|
cmpule $5, $6, $5 # e0 :
|
||||||
|
beq $5, $fixup # .. e1 :
|
||||||
|
|
||||||
|
mult $f16, $f12, $f11 # fm : $f11 = x * 0.5
|
||||||
|
subl $2, $1, $2 # .. e0 :
|
||||||
addt $f12, $f12, $f17 # .. fa : $f17 = 1.0
|
addt $f12, $f12, $f17 # .. fa : $f17 = 1.0
|
||||||
srl $2, 12, $1 # e0 :
|
srl $2, 12, $1 # e0 :
|
||||||
and $1, 0xfc, $1 # .. e1 :
|
|
||||||
addq $1, $4, $1 # e0 :
|
and $1, 0xfc, $1 # e0 :
|
||||||
ldl $1, $T2($1) # .. e1 :
|
addq $1, $4, $1 # e1 :
|
||||||
addt $f12, $f17, $f15 # fa : $f15 = 1.5
|
ldl $1, $T2($1) # e0 :
|
||||||
subl $2, $1, $2 # .. e1 :
|
addt $f12, $f17, $f15 # .. fa : $f15 = 1.5
|
||||||
sll $2, 32, $2 # e0 :
|
|
||||||
|
subl $2, $1, $2 # e0 :
|
||||||
ldt $f14, $DN($4) # .. e1 :
|
ldt $f14, $DN($4) # .. e1 :
|
||||||
|
sll $2, 32, $2 # e0 :
|
||||||
stq $2, $Y($sp) # e0 :
|
stq $2, $Y($sp) # e0 :
|
||||||
nop # .. e1 : avoid pipe flash
|
|
||||||
nop # e0 :
|
|
||||||
ldt $f13, $Y($sp) # .. e1 :
|
|
||||||
|
|
||||||
mult/su $f11, $f13, $f10 # fm : $f10 = (x * 0.5) * y
|
ldt $f13, $Y($sp) # e0 :
|
||||||
mult $f10, $f13, $f10 # fm : $f10 = ((x * 0.5) * y) * y
|
mult/su $f11, $f13, $f10 # fm 2: $f10 = (x * 0.5) * y
|
||||||
subt $f15, $f10, $f1 # fa : $f1 = (1.5 - 0.5*x*y*y)
|
mult $f10, $f13, $f10 # fm 4: $f10 = ((x * 0.5) * y) * y
|
||||||
mult $f13, $f1, $f13 # fm : yp = y*(1.5 - 0.5*x*y*y)
|
subt $f15, $f10, $f1 # fa 4: $f1 = (1.5 - 0.5*x*y*y)
|
||||||
mult/su $f11, $f13, $f1 # fm : $f11 = x * 0.5 * yp
|
|
||||||
mult $f1, $f13, $f11 # fm : $f11 = (x * 0.5 * yp) * yp
|
mult $f13, $f1, $f13 # fm 4: yp = y*(1.5 - 0.5*x*y*y)
|
||||||
subt $f18, $f11, $f1 # fa : $f1= (1.5-2^-30) - 0.5*x*yp*yp
|
mult/su $f11, $f13, $f1 # fm 4: $f11 = x * 0.5 * yp
|
||||||
mult $f13, $f1, $f13 # fm : ypp = $f13 = yp*$f1
|
mult $f1, $f13, $f11 # fm 4: $f11 = (x * 0.5 * yp) * yp
|
||||||
subt $f15, $f12, $f1 # fa : $f1 = (1.5 - 0.5)
|
subt $f18, $f11, $f1 # fa 4: $f1= (1.5-2^-30) - 0.5*x*yp*yp
|
||||||
ldt $f15, $UP($4) # .. e1 :
|
|
||||||
mult/su $f16, $f13, $f10 # fm : z = $f10 = x * ypp
|
mult $f13, $f1, $f13 # fm 4: ypp = $f13 = yp*$f1
|
||||||
mult $f10, $f13, $f11 # fm : $f11 = z*ypp
|
subt $f15, $f12, $f1 # .. fa : $f1 = (1.5 - 0.5)
|
||||||
|
ldt $f15, $UP($4) # .. e0 :
|
||||||
|
mult/su $f16, $f13, $f10 # fm 4: z = $f10 = x * ypp
|
||||||
|
|
||||||
|
mult $f10, $f13, $f11 # fm 4: $f11 = z*ypp
|
||||||
mult $f10, $f12, $f12 # fm : $f12 = z*0.5
|
mult $f10, $f12, $f12 # fm : $f12 = z*0.5
|
||||||
subt $f1, $f11, $f1 # .. fa : $f1 = 1 - z*ypp
|
subt $f1, $f11, $f1 # fa 4: $f1 = 1 - z*ypp
|
||||||
mult $f12, $f1, $f12 # fm : $f12 = z*0.5*(1 - z*ypp)
|
mult $f12, $f1, $f12 # fm 4: $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
|
addt $f10, $f12, $f0 # fa 4: zp=res= z + z*0.5*(1 - z*ypp)
|
||||||
|
mult/c $f0, $f14, $f12 # fm 4: zmi = zp * DN
|
||||||
mult/c $f0, $f15, $f11 # fm : zpl = zp * UP
|
mult/c $f0, $f15, $f11 # fm : zpl = zp * UP
|
||||||
mult/c $f0, $f12, $f1 # fm : $f1 = zp * zmi
|
mult/c $f0, $f12, $f1 # fm : $f1 = zp * zmi
|
||||||
|
|
||||||
mult/c $f0, $f11, $f15 # fm : $f15 = zp * zpl
|
mult/c $f0, $f11, $f15 # fm : $f15 = zp * zpl
|
||||||
|
subt/su $f1, $f16, $f13 # .. fa : y1 = zp*zmi - x
|
||||||
|
subt/su $f15, $f16, $f14 # fa 4: y2 = zp*zpl - x
|
||||||
|
fcmovge $f13, $f12, $f0 # fa 3: res = (y1 >= 0) ? zmi : res
|
||||||
|
|
||||||
subt/su $f1, $f16, $f13 # fa : y1 = zp*zmi - x
|
fcmovlt $f14, $f11, $f0 # fa 4: res = (y2 < 0) ? zpl : res
|
||||||
subt/su $f15, $f16, $f14 # fa : y2 = zp*zpl - x
|
addq $sp, 16, $sp # .. e0 :
|
||||||
|
|
||||||
fcmovge $f13, $f12, $f0 # res = (y1 >= 0) ? zmi : res
|
|
||||||
fcmovlt $f14, $f11, $f0 # res = (y2 < 0) ? zpl : res
|
|
||||||
|
|
||||||
addq $sp, 16, $sp # e0 :
|
|
||||||
ret # .. e1 :
|
ret # .. e1 :
|
||||||
|
|
||||||
$negative:
|
.align 4
|
||||||
ldt $f0, $NAN($4)
|
$fixup:
|
||||||
addq $sp, 16, $sp
|
addq $sp, 16, $sp
|
||||||
ret
|
br "ASM_ALPHA_NG_SYMBOL_PREFIX"__full_ieee754_sqrt..ng
|
||||||
|
|
||||||
.end __ieee754_sqrt");
|
.end __ieee754_sqrt");
|
||||||
|
|
||||||
#endif /* _IEEE_FP */
|
static double __full_ieee754_sqrt(double) __attribute__((unused));
|
||||||
#endif /* _IEEE_FP && _IEEE_FP_INEXACT */
|
#define __ieee754_sqrt __full_ieee754_sqrt
|
||||||
|
|
||||||
|
#endif /* _IEEE_FP_INEXACT */
|
||||||
|
|
||||||
|
#include <sysdeps/libm-ieee754/e_sqrt.c>
|
||||||
|
|
|
||||||
Loading…
Reference in New Issue