This is the mail archive of the
glibc-cvs@sourceware.org
mailing list for the glibc project.
GNU C Library master sources branch master updated. glibc-2.21-557-g6471190
- From: wilco at sourceware dot org
- To: glibc-cvs at sourceware dot org
- Date: 6 Jul 2015 11:56:00 -0000
- Subject: GNU C Library master sources branch master updated. glibc-2.21-557-g6471190
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "GNU C Library master sources".
The branch, master has been updated
via 6471190491435c6ee819232f3eb71ad9638563cc (commit)
via a1b85e1d6c2fbd815818d50ba7060dd3771888e3 (commit)
from 2b1d7148e3664eeb177ae2fc91bf282d75da7623 (commit)
Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.
- Log -----------------------------------------------------------------
http://sourceware.org/git/gitweb.cgi?p=glibc.git;a=commitdiff;h=6471190491435c6ee819232f3eb71ad9638563cc
commit 6471190491435c6ee819232f3eb71ad9638563cc
Author: Wilco Dijkstra <wdijkstr@arm.com>
Date: Mon Jul 6 12:43:31 2015 +0100
Inline __ieee754_sqrt and __ieee754_sqrtf. Also add external definitions.
diff --git a/ChangeLog b/ChangeLog
index f057e9d..0f5d45e 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,5 +1,14 @@
2015-07-06 Wilco Dijkstra <wdijkstr@arm.com>
+ * sysdeps/aarch64/fpu/math_private.h (__ieee754_sqrt):
+ New function. (__ieee754_sqrtf): New function.
+ * sysdeps/aarch64/fpu/e_sqrt.c (__ieee754_sqrt):
+ New function.
+ * sysdeps/aarch64/fpu/e_sqrtf.c (__ieee754_sqrtf):
+ New function.
+
+2015-07-06 Wilco Dijkstra <wdijkstr@arm.com>
+
* sysdeps/ieee754/support.c: Remove unused file.
2015-07-06 Szabolcs Nagy <szabolcs.nagy@arm.com>
diff --git a/sysdeps/aarch64/fpu/e_sqrt.c b/sysdeps/aarch64/fpu/e_sqrt.c
new file mode 100644
index 0000000..4f11ca2
--- /dev/null
+++ b/sysdeps/aarch64/fpu/e_sqrt.c
@@ -0,0 +1,28 @@
+/* Square root of floating point number.
+ Copyright (C) 2015 Free Software Foundation, Inc.
+ 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 Lesser General Public
+ License as published by the Free Software Foundation; either
+ version 2.1 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
+ Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public
+ License along with the GNU C Library; if not, see
+ <http://www.gnu.org/licenses/>. */
+
+#include <math_private.h>
+
+double
+__ieee754_sqrt (double d)
+{
+ double res;
+ asm ("fsqrt %d0, %d1" : "=w" (res) : "w" (d));
+ return res;
+}
+strong_alias (__ieee754_sqrt, __sqrt_finite)
diff --git a/sysdeps/aarch64/fpu/e_sqrtf.c b/sysdeps/aarch64/fpu/e_sqrtf.c
new file mode 100644
index 0000000..a2e99e1
--- /dev/null
+++ b/sysdeps/aarch64/fpu/e_sqrtf.c
@@ -0,0 +1,28 @@
+/* Single-precision floating point square root.
+ Copyright (C) 2015 Free Software Foundation, Inc.
+ 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 Lesser General Public
+ License as published by the Free Software Foundation; either
+ version 2.1 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
+ Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public
+ License along with the GNU C Library; if not, see
+ <http://www.gnu.org/licenses/>. */
+
+#include <math_private.h>
+
+float
+__ieee754_sqrtf (float s)
+{
+ float res;
+ asm ("fsqrt %s0, %s1" : "=w" (res) : "w" (s));
+ return res;
+}
+strong_alias (__ieee754_sqrtf, __sqrtf_finite)
diff --git a/sysdeps/aarch64/fpu/math_private.h b/sysdeps/aarch64/fpu/math_private.h
index 52a6ad9..b3c2509 100644
--- a/sysdeps/aarch64/fpu/math_private.h
+++ b/sysdeps/aarch64/fpu/math_private.h
@@ -22,6 +22,22 @@
#include <fenv.h>
#include <fpu_control.h>
+extern __always_inline double
+__ieee754_sqrt (double d)
+{
+ double res;
+ asm __volatile__ ("fsqrt %d0, %d1" : "=w" (res) : "w" (d));
+ return res;
+}
+
+extern __always_inline float
+__ieee754_sqrtf (float s)
+{
+ float res;
+ asm __volatile__ ("fsqrt %s0, %s1" : "=w" (res) : "w" (s));
+ return res;
+}
+
static __always_inline void
libc_feholdexcept_aarch64 (fenv_t *envp)
{
http://sourceware.org/git/gitweb.cgi?p=glibc.git;a=commitdiff;h=a1b85e1d6c2fbd815818d50ba7060dd3771888e3
commit a1b85e1d6c2fbd815818d50ba7060dd3771888e3
Author: Wilco Dijkstra <wdijkstr@arm.com>
Date: Thu Jun 4 13:11:43 2015 +0100
Remove unused file sysdeps/ieee754/support.c
diff --git a/ChangeLog b/ChangeLog
index 1c9a42f..f057e9d 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,7 @@
+2015-07-06 Wilco Dijkstra <wdijkstr@arm.com>
+
+ * sysdeps/ieee754/support.c: Remove unused file.
+
2015-07-06 Szabolcs Nagy <szabolcs.nagy@arm.com>
* sysdeps/unix/sysv/linux/aarch64/getcontext.S (__getcontext): Use q
diff --git a/sysdeps/ieee754/support.c b/sysdeps/ieee754/support.c
deleted file mode 100644
index 00476c0..0000000
--- a/sysdeps/ieee754/support.c
+++ /dev/null
@@ -1,520 +0,0 @@
-/*
- * Copyright (c) 1985, 1993
- * The Regents of the University of California. All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- * notice, this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- * notice, this list of conditions and the following disclaimer in the
- * documentation and/or other materials provided with the distribution.
- * 4. Neither the name of the University nor the names of its contributors
- * may be used to endorse or promote products derived from this software
- * without specific prior written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- */
-
-#ifndef lint
-static char sccsid[] = "@(#)support.c 8.1 (Berkeley) 6/4/93";
-#endif /* not lint */
-
-/*
- * Some IEEE standard 754 recommended functions and remainder and sqrt for
- * supporting the C elementary functions.
- ******************************************************************************
- * WARNING:
- * These codes are developed (in double) to support the C elementary
- * functions temporarily. They are not universal, and some of them are very
- * slow (in particular, drem and sqrt is extremely inefficient). Each
- * computer system should have its implementation of these functions using
- * its own assembler.
- ******************************************************************************
- *
- * IEEE 754 required operations:
- * drem(x,p)
- * returns x REM y = x - [x/y]*y , where [x/y] is the integer
- * nearest x/y; in half way case, choose the even one.
- * sqrt(x)
- * returns the square root of x correctly rounded according to
- * the rounding mod.
- *
- * IEEE 754 recommended functions:
- * (a) copysign(x,y)
- * returns x with the sign of y.
- * (b) scalb(x,N)
- * returns x * (2**N), for integer values N.
- * (c) logb(x)
- * returns the unbiased exponent of x, a signed integer in
- * double precision, except that logb(0) is -INF, logb(INF)
- * is +INF, and logb(NAN) is that NAN.
- * (d) finite(x)
- * returns the value TRUE if -INF < x < +INF and returns
- * FALSE otherwise.
- *
- *
- * CODED IN C BY K.C. NG, 11/25/84;
- * REVISED BY K.C. NG on 1/22/85, 2/13/85, 3/24/85.
- */
-
-#include "mathimpl.h"
-
-#if defined(vax)||defined(tahoe) /* VAX D format */
-#include <errno.h>
- static const unsigned short msign=0x7fff , mexp =0x7f80 ;
- static const short prep1=57, gap=7, bias=129 ;
- static const double novf=1.7E38, nunf=3.0E-39, zero=0.0 ;
-#else /* defined(vax)||defined(tahoe) */
- static const unsigned short msign=0x7fff, mexp =0x7ff0 ;
- static const short prep1=54, gap=4, bias=1023 ;
- static const double novf=1.7E308, nunf=3.0E-308,zero=0.0;
-#endif /* defined(vax)||defined(tahoe) */
-
-double scalb(x,N)
-double x; int N;
-{
- int k;
-
-#ifdef national
- unsigned short *px=(unsigned short *) &x + 3;
-#else /* national */
- unsigned short *px=(unsigned short *) &x;
-#endif /* national */
-
- if( x == zero ) return(x);
-
-#if defined(vax)||defined(tahoe)
- if( (k= *px & mexp ) != ~msign ) {
- if (N < -260)
- return(nunf*nunf);
- else if (N > 260) {
- return(copysign(infnan(ERANGE),x));
- }
-#else /* defined(vax)||defined(tahoe) */
- if( (k= *px & mexp ) != mexp ) {
- if( N<-2100) return(nunf*nunf); else if(N>2100) return(novf+novf);
- if( k == 0 ) {
- x *= scalb(1.0,(int)prep1); N -= prep1; return(scalb(x,N));}
-#endif /* defined(vax)||defined(tahoe) */
-
- if((k = (k>>gap)+ N) > 0 )
- if( k < (mexp>>gap) ) *px = (*px&~mexp) | (k<<gap);
- else x=novf+novf; /* overflow */
- else
- if( k > -prep1 )
- /* gradual underflow */
- {*px=(*px&~mexp)|(short)(1<<gap); x *= scalb(1.0,k-1);}
- else
- return(nunf*nunf);
- }
- return(x);
-}
-
-
-double copysign(x,y)
-double x,y;
-{
-#ifdef national
- unsigned short *px=(unsigned short *) &x+3,
- *py=(unsigned short *) &y+3;
-#else /* national */
- unsigned short *px=(unsigned short *) &x,
- *py=(unsigned short *) &y;
-#endif /* national */
-
-#if defined(vax)||defined(tahoe)
- if ( (*px & mexp) == 0 ) return(x);
-#endif /* defined(vax)||defined(tahoe) */
-
- *px = ( *px & msign ) | ( *py & ~msign );
- return(x);
-}
-
-double logb(x)
-double x;
-{
-
-#ifdef national
- short *px=(short *) &x+3, k;
-#else /* national */
- short *px=(short *) &x, k;
-#endif /* national */
-
-#if defined(vax)||defined(tahoe)
- return (int)(((*px&mexp)>>gap)-bias);
-#else /* defined(vax)||defined(tahoe) */
- if( (k= *px & mexp ) != mexp )
- if ( k != 0 )
- return ( (k>>gap) - bias );
- else if( x != zero)
- return ( -1022.0 );
- else
- return(-(1.0/zero));
- else if(x != x)
- return(x);
- else
- {*px &= msign; return(x);}
-#endif /* defined(vax)||defined(tahoe) */
-}
-
-finite(x)
-double x;
-{
-#if defined(vax)||defined(tahoe)
- return(1);
-#else /* defined(vax)||defined(tahoe) */
-#ifdef national
- return( (*((short *) &x+3 ) & mexp ) != mexp );
-#else /* national */
- return( (*((short *) &x ) & mexp ) != mexp );
-#endif /* national */
-#endif /* defined(vax)||defined(tahoe) */
-}
-
-double drem(x,p)
-double x,p;
-{
- short sign;
- double hp,dp,tmp;
- unsigned short k;
-#ifdef national
- unsigned short
- *px=(unsigned short *) &x +3,
- *pp=(unsigned short *) &p +3,
- *pd=(unsigned short *) &dp +3,
- *pt=(unsigned short *) &tmp+3;
-#else /* national */
- unsigned short
- *px=(unsigned short *) &x ,
- *pp=(unsigned short *) &p ,
- *pd=(unsigned short *) &dp ,
- *pt=(unsigned short *) &tmp;
-#endif /* national */
-
- *pp &= msign ;
-
-#if defined(vax)||defined(tahoe)
- if( ( *px & mexp ) == ~msign ) /* is x a reserved operand? */
-#else /* defined(vax)||defined(tahoe) */
- if( ( *px & mexp ) == mexp )
-#endif /* defined(vax)||defined(tahoe) */
- return (x-p)-(x-p); /* create nan if x is inf */
- if (p == zero) {
-#if defined(vax)||defined(tahoe)
- return(infnan(EDOM));
-#else /* defined(vax)||defined(tahoe) */
- return zero/zero;
-#endif /* defined(vax)||defined(tahoe) */
- }
-
-#if defined(vax)||defined(tahoe)
- if( ( *pp & mexp ) == ~msign ) /* is p a reserved operand? */
-#else /* defined(vax)||defined(tahoe) */
- if( ( *pp & mexp ) == mexp )
-#endif /* defined(vax)||defined(tahoe) */
- { if (p != p) return p; else return x;}
-
- else if ( ((*pp & mexp)>>gap) <= 1 )
- /* subnormal p, or almost subnormal p */
- { double b; b=scalb(1.0,(int)prep1);
- p *= b; x = drem(x,p); x *= b; return(drem(x,p)/b);}
- else if ( p >= novf/2)
- { p /= 2 ; x /= 2; return(drem(x,p)*2);}
- else
- {
- dp=p+p; hp=p/2;
- sign= *px & ~msign ;
- *px &= msign ;
- while ( x > dp )
- {
- k=(*px & mexp) - (*pd & mexp) ;
- tmp = dp ;
- *pt += k ;
-
-#if defined(vax)||defined(tahoe)
- if( x < tmp ) *pt -= 128 ;
-#else /* defined(vax)||defined(tahoe) */
- if( x < tmp ) *pt -= 16 ;
-#endif /* defined(vax)||defined(tahoe) */
-
- x -= tmp ;
- }
- if ( x > hp )
- { x -= p ; if ( x >= hp ) x -= p ; }
-
-#if defined(vax)||defined(tahoe)
- if (x)
-#endif /* defined(vax)||defined(tahoe) */
- *px ^= sign;
- return( x);
-
- }
-}
-
-
-double sqrt(x)
-double x;
-{
- double q,s,b,r;
- double t;
- double const zero=0.0;
- int m,n,i;
-#if defined(vax)||defined(tahoe)
- int k=54;
-#else /* defined(vax)||defined(tahoe) */
- int k=51;
-#endif /* defined(vax)||defined(tahoe) */
-
- /* sqrt(NaN) is NaN, sqrt(+-0) = +-0 */
- if(x!=x||x==zero) return(x);
-
- /* sqrt(negative) is invalid */
- if(x<zero) {
-#if defined(vax)||defined(tahoe)
- return (infnan(EDOM)); /* NaN */
-#else /* defined(vax)||defined(tahoe) */
- return(zero/zero);
-#endif /* defined(vax)||defined(tahoe) */
- }
-
- /* sqrt(INF) is INF */
- if(!finite(x)) return(x);
-
- /* scale x to [1,4) */
- n=logb(x);
- x=scalb(x,-n);
- if((m=logb(x))!=0) x=scalb(x,-m); /* subnormal number */
- m += n;
- n = m/2;
- if((n+n)!=m) {x *= 2; m -=1; n=m/2;}
-
- /* generate sqrt(x) bit by bit (accumulating in q) */
- q=1.0; s=4.0; x -= 1.0; r=1;
- for(i=1;i<=k;i++) {
- t=s+1; x *= 4; r /= 2;
- if(t<=x) {
- s=t+t+2, x -= t; q += r;}
- else
- s *= 2;
- }
-
- /* generate the last bit and determine the final rounding */
- r/=2; x *= 4;
- if(x==zero) goto end; 100+r; /* trigger inexact flag */
- if(s<x) {
- q+=r; x -=s; s += 2; s *= 2; x *= 4;
- t = (x-s)-5;
- b=1.0+3*r/4; if(b==1.0) goto end; /* b==1 : Round-to-zero */
- b=1.0+r/4; if(b>1.0) t=1; /* b>1 : Round-to-(+INF) */
- if(t>=0) q+=r; } /* else: Round-to-nearest */
- else {
- s *= 2; x *= 4;
- t = (x-s)-1;
- b=1.0+3*r/4; if(b==1.0) goto end;
- b=1.0+r/4; if(b>1.0) t=1;
- if(t>=0) q+=r; }
-
-end: return(scalb(q,n));
-}
-
-#if 0
-/* DREM(X,Y)
- * RETURN X REM Y =X-N*Y, N=[X/Y] ROUNDED (ROUNDED TO EVEN IN THE HALF WAY CASE)
- * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
- * INTENDED FOR ASSEMBLY LANGUAGE
- * CODED IN C BY K.C. NG, 3/23/85, 4/8/85.
- *
- * Warning: this code should not get compiled in unless ALL of
- * the following machine-dependent routines are supplied.
- *
- * Required machine dependent functions (not on a VAX):
- * swapINX(i): save inexact flag and reset it to "i"
- * swapENI(e): save inexact enable and reset it to "e"
- */
-
-double drem(x,y)
-double x,y;
-{
-
-#ifdef national /* order of words in floating point number */
- static const n0=3,n1=2,n2=1,n3=0;
-#else /* VAX, SUN, ZILOG, TAHOE */
- static const n0=0,n1=1,n2=2,n3=3;
-#endif
-
- static const unsigned short mexp =0x7ff0, m25 =0x0190, m57 =0x0390;
- static const double zero=0.0;
- double hy,y1,t,t1;
- short k;
- long n;
- int i,e;
- unsigned short xexp,yexp, *px =(unsigned short *) &x ,
- nx,nf, *py =(unsigned short *) &y ,
- sign, *pt =(unsigned short *) &t ,
- *pt1 =(unsigned short *) &t1 ;
-
- xexp = px[n0] & mexp ; /* exponent of x */
- yexp = py[n0] & mexp ; /* exponent of y */
- sign = px[n0] &0x8000; /* sign of x */
-
-/* return NaN if x is NaN, or y is NaN, or x is INF, or y is zero */
- if(x!=x) return(x); if(y!=y) return(y); /* x or y is NaN */
- if( xexp == mexp ) return(zero/zero); /* x is INF */
- if(y==zero) return(y/y);
-
-/* save the inexact flag and inexact enable in i and e respectively
- * and reset them to zero
- */
- i=swapINX(0); e=swapENI(0);
-
-/* subnormal number */
- nx=0;
- if(yexp==0) {t=1.0,pt[n0]+=m57; y*=t; nx=m57;}
-
-/* if y is tiny (biased exponent <= 57), scale up y to y*2**57 */
- if( yexp <= m57 ) {py[n0]+=m57; nx+=m57; yexp+=m57;}
-
- nf=nx;
- py[n0] &= 0x7fff;
- px[n0] &= 0x7fff;
-
-/* mask off the least significant 27 bits of y */
- t=y; pt[n3]=0; pt[n2]&=0xf800; y1=t;
-
-/* LOOP: argument reduction on x whenever x > y */
-loop:
- while ( x > y )
- {
- t=y;
- t1=y1;
- xexp=px[n0]&mexp; /* exponent of x */
- k=xexp-yexp-m25;
- if(k>0) /* if x/y >= 2**26, scale up y so that x/y < 2**26 */
- {pt[n0]+=k;pt1[n0]+=k;}
- n=x/t; x=(x-n*t1)-n*(t-t1);
- }
- /* end while (x > y) */
-
- if(nx!=0) {t=1.0; pt[n0]+=nx; x*=t; nx=0; goto loop;}
-
-/* final adjustment */
-
- hy=y/2.0;
- if(x>hy||((x==hy)&&n%2==1)) x-=y;
- px[n0] ^= sign;
- if(nf!=0) { t=1.0; pt[n0]-=nf; x*=t;}
-
-/* restore inexact flag and inexact enable */
- swapINX(i); swapENI(e);
-
- return(x);
-}
-#endif
-
-#if 0
-/* SQRT
- * RETURN CORRECTLY ROUNDED (ACCORDING TO THE ROUNDING MODE) SQRT
- * FOR IEEE DOUBLE PRECISION ONLY, INTENDED FOR ASSEMBLY LANGUAGE
- * CODED IN C BY K.C. NG, 3/22/85.
- *
- * Warning: this code should not get compiled in unless ALL of
- * the following machine-dependent routines are supplied.
- *
- * Required machine dependent functions:
- * swapINX(i) ...return the status of INEXACT flag and reset it to "i"
- * swapRM(r) ...return the current Rounding Mode and reset it to "r"
- * swapENI(e) ...return the status of inexact enable and reset it to "e"
- * addc(t) ...perform t=t+1 regarding t as a 64 bit unsigned integer
- * subc(t) ...perform t=t-1 regarding t as a 64 bit unsigned integer
- */
-
-static const unsigned long table[] = {
-0, 1204, 3062, 5746, 9193, 13348, 18162, 23592, 29598, 36145, 43202, 50740,
-58733, 67158, 75992, 85215, 83599, 71378, 60428, 50647, 41945, 34246, 27478,
-21581, 16499, 12183, 8588, 5674, 3403, 1742, 661, 130, };
-
-double newsqrt(x)
-double x;
-{
- double y,z,t,addc(),subc()
- double const b54=134217728.*134217728.; /* b54=2**54 */
- long mx,scalx;
- long const mexp=0x7ff00000;
- int i,j,r,e,swapINX(),swapRM(),swapENI();
- unsigned long *py=(unsigned long *) &y ,
- *pt=(unsigned long *) &t ,
- *px=(unsigned long *) &x ;
-#ifdef national /* ordering of word in a floating point number */
- const int n0=1, n1=0;
-#else
- const int n0=0, n1=1;
-#endif
-/* Rounding Mode: RN ...round-to-nearest
- * RZ ...round-towards 0
- * RP ...round-towards +INF
- * RM ...round-towards -INF
- */
- const int RN=0,RZ=1,RP=2,RM=3;
- /* machine dependent: work on a Zilog Z8070
- * and a National 32081 & 16081
- */
-
-/* exceptions */
- if(x!=x||x==0.0) return(x); /* sqrt(NaN) is NaN, sqrt(+-0) = +-0 */
- if(x<0) return((x-x)/(x-x)); /* sqrt(negative) is invalid */
- if((mx=px[n0]&mexp)==mexp) return(x); /* sqrt(+INF) is +INF */
-
-/* save, reset, initialize */
- e=swapENI(0); /* ...save and reset the inexact enable */
- i=swapINX(0); /* ...save INEXACT flag */
- r=swapRM(RN); /* ...save and reset the Rounding Mode to RN */
- scalx=0;
-
-/* subnormal number, scale up x to x*2**54 */
- if(mx==0) {x *= b54 ; scalx-=0x01b00000;}
-
-/* scale x to avoid intermediate over/underflow:
- * if (x > 2**512) x=x/2**512; if (x < 2**-512) x=x*2**512 */
- if(mx>0x5ff00000) {px[n0] -= 0x20000000; scalx+= 0x10000000;}
- if(mx<0x1ff00000) {px[n0] += 0x20000000; scalx-= 0x10000000;}
-
-/* magic initial approximation to almost 8 sig. bits */
- py[n0]=(px[n0]>>1)+0x1ff80000;
- py[n0]=py[n0]-table[(py[n0]>>15)&31];
-
-/* Heron's rule once with correction to improve y to almost 18 sig. bits */
- t=x/y; y=y+t; py[n0]=py[n0]-0x00100006; py[n1]=0;
-
-/* triple to almost 56 sig. bits; now y approx. sqrt(x) to within 1 ulp */
- t=y*y; z=t; pt[n0]+=0x00100000; t+=z; z=(x-z)*y;
- t=z/(t+x) ; pt[n0]+=0x00100000; y+=t;
-
-/* twiddle last bit to force y correctly rounded */
- swapRM(RZ); /* ...set Rounding Mode to round-toward-zero */
- swapINX(0); /* ...clear INEXACT flag */
- swapENI(e); /* ...restore inexact enable status */
- t=x/y; /* ...chopped quotient, possibly inexact */
- j=swapINX(i); /* ...read and restore inexact flag */
- if(j==0) { if(t==y) goto end; else t=subc(t); } /* ...t=t-ulp */
- b54+0.1; /* ..trigger inexact flag, sqrt(x) is inexact */
- if(r==RN) t=addc(t); /* ...t=t+ulp */
- else if(r==RP) { t=addc(t);y=addc(y);}/* ...t=t+ulp;y=y+ulp; */
- y=y+t; /* ...chopped sum */
- py[n0]=py[n0]-0x00100000; /* ...correctly rounded sqrt(x) */
-end: py[n0]=py[n0]+scalx; /* ...scale back y */
- swapRM(r); /* ...restore Rounding Mode */
- return(y);
-}
-#endif
-----------------------------------------------------------------------
Summary of changes:
ChangeLog | 13 +
sysdeps/{x86_64 => aarch64}/fpu/e_sqrt.c | 9 +-
sysdeps/{s390 => aarch64}/fpu/e_sqrtf.c | 9 +-
sysdeps/aarch64/fpu/math_private.h | 16 +
sysdeps/ieee754/support.c | 520 ------------------------------
5 files changed, 36 insertions(+), 531 deletions(-)
copy sysdeps/{x86_64 => aarch64}/fpu/e_sqrt.c (84%)
copy sysdeps/{s390 => aarch64}/fpu/e_sqrtf.c (81%)
delete mode 100644 sysdeps/ieee754/support.c
hooks/post-receive
--
GNU C Library master sources