diff options
| author | Samuel Thibault <samuel.thibault@ens-lyon.org> | 2016-08-20 19:14:56 +0200 |
|---|---|---|
| committer | Samuel Thibault <samuel.thibault@ens-lyon.org> | 2016-08-20 19:14:56 +0200 |
| commit | f76453c31593957fec1a99b986bfa5506618b79c (patch) | |
| tree | da353c882fb9b2261c9871bcb9e3876a3e6ed7f6 /sysdeps/ieee754/dbl-64 | |
| parent | 58695b88a9deaecbcf7794760cc333177edaa2b4 (diff) | |
| parent | 78bd7499af46d739ce94410eaeea006e874ca9e5 (diff) | |
Merge tag 'glibc-2.22' into baseline
The GNU C Library
=================
The GNU C Library version 2.22 is now available.
The GNU C Library is used as *the* C library in the GNU system and
in GNU/Linux systems, as well as many other systems that use Linux
as the kernel.
The GNU C Library is primarily designed to be a portable
and high performance C library. It follows all relevant
standards including ISO C11 and POSIX.1-2008. It is also
internationalized and has one of the most complete
internationalization interfaces known.
The GNU C Library webpage is at http://www.gnu.org/software/libc/
Packages for the 2.22 release may be downloaded from:
http://ftpmirror.gnu.org/libc/
http://ftp.gnu.org/gnu/libc/
The mirror list is at http://www.gnu.org/order/ftp.html
NEWS for version 2.22
=====================
* The following bugs are resolved with this release:
438, 4719, 6544, 6792, 11216, 12836, 13028, 13064, 13151, 13152, 14094,
14292, 14841, 14906, 14958, 15319, 15467, 15790, 15969, 16159, 16339,
16350, 16351, 16352, 16353, 16361, 16512, 16526, 16538, 16559, 16560,
16704, 16783, 16850, 17053, 17090, 17195, 17269, 17293, 17322, 17403,
17475, 17523, 17542, 17569, 17581, 17588, 17596, 17620, 17621, 17628,
17631, 17692, 17711, 17715, 17776, 17779, 17792, 17833, 17836, 17841,
17912, 17916, 17930, 17932, 17944, 17949, 17964, 17965, 17967, 17969,
17977, 17978, 17987, 17991, 17996, 17998, 17999, 18007, 18019, 18020,
18029, 18030, 18032, 18034, 18036, 18038, 18039, 18042, 18043, 18046,
18047, 18049, 18068, 18080, 18093, 18100, 18104, 18110, 18111, 18116,
18125, 18128, 18134, 18138, 18185, 18196, 18197, 18206, 18210, 18211,
18217, 18219, 18220, 18221, 18234, 18244, 18245, 18247, 18287, 18319,
18324, 18333, 18346, 18371, 18383, 18397, 18400, 18409, 18410, 18412,
18418, 18422, 18434, 18444, 18457, 18468, 18469, 18470, 18479, 18483,
18495, 18496, 18497, 18498, 18502, 18507, 18508, 18512, 18513, 18519,
18520, 18522, 18527, 18528, 18529, 18530, 18532, 18533, 18534, 18536,
18539, 18540, 18542, 18544, 18545, 18546, 18547, 18549, 18553, 18557,
18558, 18569, 18583, 18585, 18586, 18592, 18593, 18594, 18602, 18612,
18613, 18619, 18633, 18641, 18643, 18648, 18657, 18676, 18694, 18696.
* Cache information can be queried via sysconf() function on s390 e.g. with
_SC_LEVEL1_ICACHE_SIZE as argument.
* A buffer overflow in gethostbyname_r and related functions performing DNS
requests has been fixed. If the NSS functions were called with a
misaligned buffer, the buffer length change due to pointer alignment was
not taken into account. This could result in application crashes or,
potentially arbitrary code execution, using crafted, but syntactically
valid DNS responses. (CVE-2015-1781)
* The time zone file parser has been made more robust against crafted time
zone files, avoiding heap buffer overflows related to the processing of
the tzh_ttisstdcnt and tzh_ttisgmtcnt fields, and a stack overflow due to
large time zone data files. Overly long time zone specifiers in the TZ
variable no longer result in stack overflows and crashes.
* A powerpc and powerpc64 optimization for TLS, similar to TLS descriptors
for LD and GD on x86 and x86-64, has been implemented. You will need
binutils-2.24 or later to enable this optimization.
* Character encoding and ctype tables were updated to Unicode 7.0.0, using
new generator scripts contributed by Pravin Satpute and Mike FABIAN (Red
Hat). These updates cause user visible changes, such as the fix for bug
17998.
* CVE-2014-8121 The NSS backends shared internal state between the getXXent
and getXXbyYY NSS calls for the same database, causing a denial-of-service
condition in some applications.
* Added vector math library named libmvec with the following vectorized x86_64
implementations: cos, cosf, sin, sinf, sincos, sincosf, log, logf, exp, expf,
pow, powf.
The library can be disabled with --disable-mathvec. Use of the functions is
enabled with -fopenmp -ffast-math starting from -O1 for GCC version >= 4.9.0.
Shared library libmvec.so is linked in as needed when using -lm (no need to
specify -lmvec explicitly for not static builds).
Visit <https://sourceware.org/glibc/wiki/libmvec> for detailed information.
* A new fmemopen implementation has been added with the goal of POSIX
compliance. The new implementation fixes the following long-standing
issues: BZ#6544, BZ#11216, BZ#12836, BZ#13151, BZ#13152, and BZ#14292. The
old implementation is still present for use be by existing binaries.
* The 32-bit sparc sigaction ABI was inadvertently broken in the 2.20 and 2.21
releases. It has been fixed to match 2.19 and older, but binaries built
against 2.20 and 2.21 might need to be recompiled. See BZ#18694.
* Port to Native Client running on ARMv7-A (--host=arm-nacl).
Contributed by Roland McGrath (Google).
Contributors
============
This release was made possible by the contributions of many people.
The maintainers are grateful to everyone who has contributed
changes or bug reports. These include:
Adhemerval Zanella
Alan Modra
Alexandre Oliva
Andreas Schwab
Andrew Senkevich
Andriy Rysin
Arjun Shankar
Aurelien Jarno
Benno Schulenberg
Brad Hubbard
Carlos O'Donell
Chris Metcalf
Christian Schmidt
Chung-Lin Tang
Cong Wang
Cyril Hrubis
Daniel Marjamäki
David S. Miller
Dmitry V. Levin
Eric Rannaud
Evangelos Foutras
Feng Gao
Florian Weimer
Gleb Fotengauer-Malinovskiy
H.J. Lu
Igor Zamyatin
J William Piggott
James Cowgill
James Lemke
John David Anglin
Joseph Myers
Kevin Easton
Khem Raj
Leonhard Holz
Mark Wielaard
Marko Myllynen
Martin Galvan
Martin Sebor
Matthew Fortune
Mel Gorman
Mike Frysinger
Miroslav Lichvar
Nathan Lynch
Ondřej Bílka
Paul Eggert
Paul Pluzhnikov
Pavel Kopyl
Pravin Satpute
Rajalakshmi Srinivasaraghavan
Rical Jasan
Richard Henderson
Roland McGrath
Rüdiger Sonderfeld
Samuel Thibault
Siddhesh Poyarekar
Stefan Liebler
Steve Ellcey
Szabolcs Nagy
Torvald Riegel
Tulio Magno Quites Machado Filho
Vincent Bernat
Wilco Dijkstra
Yaakov Selkowitz
Zack Weinberg
Diffstat (limited to 'sysdeps/ieee754/dbl-64')
114 files changed, 944 insertions, 713 deletions
diff --git a/sysdeps/ieee754/dbl-64/Makefile b/sysdeps/ieee754/dbl-64/Makefile index 35f545ff8e..5557c75b45 100644 --- a/sysdeps/ieee754/dbl-64/Makefile +++ b/sysdeps/ieee754/dbl-64/Makefile @@ -2,4 +2,5 @@ ifeq ($(subdir),math) # branred depends on precise IEEE double rounding CFLAGS-branred.c = $(config-cflags-nofma) CFLAGS-e_sqrt.c = $(config-cflags-nofma) +CFLAGS-e_pow.c = $(config-cflags-nofma) endif diff --git a/sysdeps/ieee754/dbl-64/MathLib.h b/sysdeps/ieee754/dbl-64/MathLib.h index 395a57669c..8d8de6fa1b 100644 --- a/sysdeps/ieee754/dbl-64/MathLib.h +++ b/sysdeps/ieee754/dbl-64/MathLib.h @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2014 Free Software Foundation, Inc. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/asincos.tbl b/sysdeps/ieee754/dbl-64/asincos.tbl index 840b61852f..84efa43a2d 100644 --- a/sysdeps/ieee754/dbl-64/asincos.tbl +++ b/sysdeps/ieee754/dbl-64/asincos.tbl @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2014 Free Software Foundation, Inc. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/atnat.h b/sysdeps/ieee754/dbl-64/atnat.h index de612b09c7..a9c265fe41 100644 --- a/sysdeps/ieee754/dbl-64/atnat.h +++ b/sysdeps/ieee754/dbl-64/atnat.h @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2014 Free Software Foundation, Inc. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/atnat2.h b/sysdeps/ieee754/dbl-64/atnat2.h index d5a141deb0..e0d65afd06 100644 --- a/sysdeps/ieee754/dbl-64/atnat2.h +++ b/sysdeps/ieee754/dbl-64/atnat2.h @@ -2,7 +2,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2014 Free Software Foundation, Inc. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/branred.c b/sysdeps/ieee754/dbl-64/branred.c index 9d15fdf4d2..7757a9c833 100644 --- a/sysdeps/ieee754/dbl-64/branred.c +++ b/sysdeps/ieee754/dbl-64/branred.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2014 Free Software Foundation, Inc. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -35,6 +35,7 @@ #include "endian.h" #include "mydefs.h" #include "branred.h" +#include <math.h> #include <math_private.h> #ifndef SECTION @@ -123,7 +124,7 @@ __branred(double x, double *a, double *aa) sum=sum1+sum2; b=b1+b2; - bb = (ABS(b1)>ABS(b2))? (b1-b)+b2 : (b2-b)+b1; + bb = (fabs(b1)>fabs(b2))? (b1-b)+b2 : (b2-b)+b1; if (b > 0.5) {b-=1.0; sum+=1.0;} else if (b < -0.5) diff --git a/sysdeps/ieee754/dbl-64/branred.h b/sysdeps/ieee754/dbl-64/branred.h index 2306e2e023..431bce0c84 100644 --- a/sysdeps/ieee754/dbl-64/branred.h +++ b/sysdeps/ieee754/dbl-64/branred.h @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2014 Free Software Foundation, Inc. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/dbl2mpn.c b/sysdeps/ieee754/dbl-64/dbl2mpn.c index f7d8186098..bb28192883 100644 --- a/sysdeps/ieee754/dbl-64/dbl2mpn.c +++ b/sysdeps/ieee754/dbl-64/dbl2mpn.c @@ -1,4 +1,4 @@ -/* Copyright (C) 1993-2014 Free Software Foundation, Inc. +/* Copyright (C) 1993-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 diff --git a/sysdeps/ieee754/dbl-64/dla.h b/sysdeps/ieee754/dbl-64/dla.h index f90a6f8d2e..641644eeb1 100644 --- a/sysdeps/ieee754/dbl-64/dla.h +++ b/sysdeps/ieee754/dbl-64/dla.h @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2014 Free Software Foundation, Inc. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -17,6 +17,8 @@ * along with this program; if not, see <http://www.gnu.org/licenses/>. */ +#include <math.h> + /***********************************************************************/ /*MODULE_NAME: dla.h */ /* */ @@ -44,7 +46,7 @@ /* z+zz = x+y exactly. */ #define EADD(x,y,z,zz) \ - z=(x)+(y); zz=(ABS(x)>ABS(y)) ? (((x)-(z))+(y)) : (((y)-(z))+(x)); + z=(x)+(y); zz=(fabs(x)>fabs(y)) ? (((x)-(z))+(y)) : (((y)-(z))+(x)); /* Exact subtraction of two single-length floating point numbers, Dekker. */ @@ -52,7 +54,7 @@ /* z+zz = x-y exactly. */ #define ESUB(x,y,z,zz) \ - z=(x)-(y); zz=(ABS(x)>ABS(y)) ? (((x)-(z))-(y)) : ((x)-((y)+(z))); + z=(x)-(y); zz=(fabs(x)>fabs(y)) ? (((x)-(z))-(y)) : ((x)-((y)+(z))); /* Exact multiplication of two single-length floating point numbers, */ @@ -94,7 +96,7 @@ /* storage variables of type double. */ #define ADD2(x, xx, y, yy, z, zz, r, s) \ - r = (x) + (y); s = (ABS (x) > ABS (y)) ? \ + r = (x) + (y); s = (fabs (x) > fabs (y)) ? \ (((((x) - r) + (y)) + (yy)) + (xx)) : \ (((((y) - r) + (x)) + (xx)) + (yy)); \ z = r + s; zz = (r - z) + s; @@ -107,7 +109,7 @@ /* storage variables of type double. */ #define SUB2(x, xx, y, yy, z, zz, r, s) \ - r = (x) - (y); s = (ABS (x) > ABS (y)) ? \ + r = (x) - (y); s = (fabs (x) > fabs (y)) ? \ (((((x) - r) - (y)) - (yy)) + (xx)) : \ ((((x) - ((y) + r)) + (xx)) - (yy)); \ z = r + s; zz = (r - z) + s; @@ -144,16 +146,16 @@ #define ADD2A(x, xx, y, yy, z, zz, r, rr, s, ss, u, uu, w) \ r = (x) + (y); \ - if (ABS (x) > ABS (y)) { rr = ((x) - r) + (y); s = (rr + (yy)) + (xx); } \ + if (fabs (x) > fabs (y)) { rr = ((x) - r) + (y); s = (rr + (yy)) + (xx); } \ else { rr = ((y) - r) + (x); s = (rr + (xx)) + (yy); } \ if (rr != 0.0) { \ z = r + s; zz = (r - z) + s; } \ else { \ - ss = (ABS (xx) > ABS (yy)) ? (((xx) - s) + (yy)) : (((yy) - s) + (xx));\ + ss = (fabs (xx) > fabs (yy)) ? (((xx) - s) + (yy)) : (((yy) - s) + (xx));\ u = r + s; \ - uu = (ABS (r) > ABS (s)) ? ((r - u) + s) : ((s - u) + r); \ + uu = (fabs (r) > fabs (s)) ? ((r - u) + s) : ((s - u) + r); \ w = uu + ss; z = u + w; \ - zz = (ABS (u) > ABS (w)) ? ((u - z) + w) : ((w - z) + u); } + zz = (fabs (u) > fabs (w)) ? ((u - z) + w) : ((w - z) + u); } /* Double-length subtraction, slower but more accurate than SUB2. */ @@ -165,13 +167,13 @@ #define SUB2A(x, xx, y, yy, z, zz, r, rr, s, ss, u, uu, w) \ r = (x) - (y); \ - if (ABS (x) > ABS (y)) { rr = ((x) - r) - (y); s = (rr - (yy)) + (xx); } \ + if (fabs (x) > fabs (y)) { rr = ((x) - r) - (y); s = (rr - (yy)) + (xx); } \ else { rr = (x) - ((y) + r); s = (rr + (xx)) - (yy); } \ if (rr != 0.0) { \ z = r + s; zz = (r - z) + s; } \ else { \ - ss = (ABS (xx) > ABS (yy)) ? (((xx) - s) - (yy)) : ((xx) - ((yy) + s)); \ + ss = (fabs (xx) > fabs (yy)) ? (((xx) - s) - (yy)) : ((xx) - ((yy) + s)); \ u = r + s; \ - uu = (ABS (r) > ABS (s)) ? ((r - u) + s) : ((s - u) + r); \ + uu = (fabs (r) > fabs (s)) ? ((r - u) + s) : ((s - u) + r); \ w = uu + ss; z = u + w; \ - zz = (ABS (u) > ABS (w)) ? ((u - z) + w) : ((w - z) + u); } + zz = (fabs (u) > fabs (w)) ? ((u - z) + w) : ((w - z) + u); } diff --git a/sysdeps/ieee754/dbl-64/doasin.c b/sysdeps/ieee754/dbl-64/doasin.c index 132a3f0113..2372e1bc53 100644 --- a/sysdeps/ieee754/dbl-64/doasin.c +++ b/sysdeps/ieee754/dbl-64/doasin.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001-2014 Free Software Foundation, Inc. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/doasin.h b/sysdeps/ieee754/dbl-64/doasin.h index 39cfcf4aed..ac4c84367e 100644 --- a/sysdeps/ieee754/dbl-64/doasin.h +++ b/sysdeps/ieee754/dbl-64/doasin.h @@ -2,7 +2,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2014 Free Software Foundation, Inc. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/dosincos.c b/sysdeps/ieee754/dbl-64/dosincos.c index de6f57bfff..ed5e2c9027 100644 --- a/sysdeps/ieee754/dbl-64/dosincos.c +++ b/sysdeps/ieee754/dbl-64/dosincos.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001-2014 Free Software Foundation, Inc. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/dosincos.h b/sysdeps/ieee754/dbl-64/dosincos.h index 516afd3977..7d9324f736 100644 --- a/sysdeps/ieee754/dbl-64/dosincos.h +++ b/sysdeps/ieee754/dbl-64/dosincos.h @@ -2,7 +2,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2014 Free Software Foundation, Inc. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/e_asin.c b/sysdeps/ieee754/dbl-64/e_asin.c index 5bb5aeb075..a7684d1078 100644 --- a/sysdeps/ieee754/dbl-64/e_asin.c +++ b/sysdeps/ieee754/dbl-64/e_asin.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001-2014 Free Software Foundation, Inc. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -39,6 +39,8 @@ #include "powtwo.tbl" #include "MathLib.h" #include "uasncs.h" +#include <float.h> +#include <math.h> #include <math_private.h> #ifndef SECTION @@ -67,7 +69,15 @@ __ieee754_asin(double x){ m = u.i[HIGH_HALF]; k = 0x7fffffff&m; /* no sign */ - if (k < 0x3e500000) return x; /* for x->0 => sin(x)=x */ + if (k < 0x3e500000) + { + if (fabs (x) < DBL_MIN) + { + double force_underflow = x * x; + math_force_eval (force_underflow); + } + return x; /* for x->0 => sin(x)=x */ + } /*----------------------2^-26 <= |x| < 2^ -3 -----------------*/ else if (k < 0x3fc00000) { @@ -94,9 +104,9 @@ __ieee754_asin(double x){ __doasin(x,0,w); if (w[0]==(w[0]+1.00000001*w[1])) return w[0]; else { - y=ABS(x); - res=ABS(w[0]); - res1=ABS(w[0]+1.1*w[1]); + y=fabs(x); + res=fabs(w[0]); + res1=fabs(w[0]+1.1*w[1]); return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1); } } @@ -125,11 +135,11 @@ __ieee754_asin(double x){ res1=res+1.1*cor; z=0.5*(res1-res); __dubsin(res,z,w); - z=(w[0]-ABS(x))+w[1]; + z=(w[0]-fabs(x))+w[1]; if (z>1.0e-27) return (m>0)?min(res,res1):-min(res,res1); else if (z<-1.0e-27) return (m>0)?max(res,res1):-max(res,res1); else { - y=ABS(x); + y=fabs(x); return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1); } } @@ -158,11 +168,11 @@ __ieee754_asin(double x){ res1=res+1.1*cor; z=0.5*(res1-res); __dubsin(res,z,w); - z=(w[0]-ABS(x))+w[1]; + z=(w[0]-fabs(x))+w[1]; if (z>1.0e-27) return (m>0)?min(res,res1):-min(res,res1); else if (z<-1.0e-27) return (m>0)?max(res,res1):-max(res,res1); else { - y=ABS(x); + y=fabs(x); return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1); } } @@ -193,11 +203,11 @@ __ieee754_asin(double x){ y=hp0.x-res; z=((hp0.x-y)-res)+(hp1.x-z); __dubcos(y,z,w); - z=(w[0]-ABS(x))+w[1]; + z=(w[0]-fabs(x))+w[1]; if (z>1.0e-27) return (m>0)?min(res,res1):-min(res,res1); else if (z<-1.0e-27) return (m>0)?max(res,res1):-max(res,res1); else { - y=ABS(x); + y=fabs(x); return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1); } } @@ -231,11 +241,11 @@ __ieee754_asin(double x){ z=y+hp1.x; y=(y-z)+hp1.x; __dubcos(z,y,w); - z=(w[0]-ABS(x))+w[1]; + z=(w[0]-fabs(x))+w[1]; if (z>1.0e-27) return (m>0)?min(res,res1):-min(res,res1); else if (z<-1.0e-27) return (m>0)?max(res,res1):-max(res,res1); else { - y=ABS(x); + y=fabs(x); return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1); } } @@ -270,11 +280,11 @@ __ieee754_asin(double x){ z=y+hp1.x; y=(y-z)+hp1.x; __dubcos(z,y,w); - z=(w[0]-ABS(x))+w[1]; + z=(w[0]-fabs(x))+w[1]; if (z>1.0e-27) return (m>0)?min(res,res1):-min(res,res1); else if (z<-1.0e-27) return (m>0)?max(res,res1):-max(res,res1); else { - y=ABS(x); + y=fabs(x); return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1); } } @@ -308,7 +318,7 @@ __ieee754_asin(double x){ cor = (res1-res)+cor; if (res==(res+1.0000001*cor)) return (m>0)?res:-res; else { - y=ABS(x); + y=fabs(x); res1=res+1.1*cor; return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1); } diff --git a/sysdeps/ieee754/dbl-64/e_atan2.c b/sysdeps/ieee754/dbl-64/e_atan2.c index a287ca6656..0c2902855e 100644 --- a/sysdeps/ieee754/dbl-64/e_atan2.c +++ b/sysdeps/ieee754/dbl-64/e_atan2.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001-2014 Free Software Foundation, Inc. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -41,6 +41,9 @@ #include "MathLib.h" #include "uatan.tbl" #include "atnat2.h" +#include <fenv.h> +#include <float.h> +#include <math.h> #include <math_private.h> #include <stap-probe.h> @@ -190,6 +193,7 @@ __ieee754_atan2 (double y, double x) return mhpi.d; } + SET_RESTORE_ROUND (FE_TONEAREST); /* either x/y or y/x is very close to zero */ ax = (x < 0) ? -x : x; ay = (y < 0) ? -y : y; @@ -202,10 +206,18 @@ __ieee754_atan2 (double y, double x) { if (x > 0) { + double ret; if ((z = ay / ax) < TWOM1022) - return normalized (ax, ay, y, z); + ret = normalized (ax, ay, y, z); else - return signArctan2 (y, z); + ret = signArctan2 (y, z); + if (fabs (ret) < DBL_MIN) + { + double vret = ret ? ret : DBL_MIN; + double force_underflow = vret * vret; + math_force_eval (force_underflow); + } + return ret; } else { diff --git a/sysdeps/ieee754/dbl-64/e_atanh.c b/sysdeps/ieee754/dbl-64/e_atanh.c index 21bb9908d1..6b00b800f2 100644 --- a/sysdeps/ieee754/dbl-64/e_atanh.c +++ b/sysdeps/ieee754/dbl-64/e_atanh.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2011-2014 Free Software Foundation, Inc. +/* Copyright (C) 2011-2015 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@gmail.com>, 2011. @@ -35,6 +35,7 @@ */ +#include <float.h> #include <inttypes.h> #include <math.h> #include <math_private.h> @@ -48,16 +49,21 @@ __ieee754_atanh (double x) double t; if (isless (xa, 0.5)) { - if (__builtin_expect (xa < 0x1.0p-28, 0)) + if (__glibc_unlikely (xa < 0x1.0p-28)) { math_force_eval (huge + x); + if (fabs (x) < DBL_MIN) + { + double force_underflow = x * x; + math_force_eval (force_underflow); + } return x; } t = xa + xa; t = 0.5 * __log1p (t + t * xa / (1.0 - xa)); } - else if (__builtin_expect (isless (xa, 1.0), 1)) + else if (__glibc_likely (isless (xa, 1.0))) t = 0.5 * __log1p ((xa + xa) / (1.0 - xa)); else { diff --git a/sysdeps/ieee754/dbl-64/e_cosh.c b/sysdeps/ieee754/dbl-64/e_cosh.c index 6caf943ed0..af3910dd6e 100644 --- a/sysdeps/ieee754/dbl-64/e_cosh.c +++ b/sysdeps/ieee754/dbl-64/e_cosh.c @@ -53,10 +53,10 @@ __ieee754_cosh (double x) /* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */ if (ix < 0x3fd62e43) { + if (ix < 0x3c800000) + return one; /* cosh(tiny) = 1 */ t = __expm1 (fabs (x)); w = one + t; - if (ix < 0x3c800000) - return w; /* cosh(tiny) = 1 */ return one + (t * t) / (w + w); } diff --git a/sysdeps/ieee754/dbl-64/e_exp.c b/sysdeps/ieee754/dbl-64/e_exp.c index 100fc397a7..bb76907f74 100644 --- a/sysdeps/ieee754/dbl-64/e_exp.c +++ b/sysdeps/ieee754/dbl-64/e_exp.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001-2014 Free Software Foundation, Inc. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -32,6 +32,7 @@ /* */ /***************************************************************************/ +#include <math.h> #include "endian.h" #include "uexp.h" #include "mydefs.h" @@ -58,168 +59,178 @@ __ieee754_exp (double x) int4 i, j, m, n, ex; double retval; - SET_RESTORE_ROUND (FE_TONEAREST); + { + SET_RESTORE_ROUND (FE_TONEAREST); - junk1.x = x; - m = junk1.i[HIGH_HALF]; - n = m & hugeint; + junk1.x = x; + m = junk1.i[HIGH_HALF]; + n = m & hugeint; - if (n > smallint && n < bigint) - { - y = x * log2e.x + three51.x; - bexp = y - three51.x; /* multiply the result by 2**bexp */ + if (n > smallint && n < bigint) + { + y = x * log2e.x + three51.x; + bexp = y - three51.x; /* multiply the result by 2**bexp */ - junk1.x = y; + junk1.x = y; - eps = bexp * ln_two2.x; /* x = bexp*ln(2) + t - eps */ - t = x - bexp * ln_two1.x; + eps = bexp * ln_two2.x; /* x = bexp*ln(2) + t - eps */ + t = x - bexp * ln_two1.x; - y = t + three33.x; - base = y - three33.x; /* t rounded to a multiple of 2**-18 */ - junk2.x = y; - del = (t - base) - eps; /* x = bexp*ln(2) + base + del */ - eps = del + del * del * (p3.x * del + p2.x); + y = t + three33.x; + base = y - three33.x; /* t rounded to a multiple of 2**-18 */ + junk2.x = y; + del = (t - base) - eps; /* x = bexp*ln(2) + base + del */ + eps = del + del * del * (p3.x * del + p2.x); - binexp.i[HIGH_HALF] = (junk1.i[LOW_HALF] + 1023) << 20; + binexp.i[HIGH_HALF] = (junk1.i[LOW_HALF] + 1023) << 20; - i = ((junk2.i[LOW_HALF] >> 8) & 0xfffffffe) + 356; - j = (junk2.i[LOW_HALF] & 511) << 1; + i = ((junk2.i[LOW_HALF] >> 8) & 0xfffffffe) + 356; + j = (junk2.i[LOW_HALF] & 511) << 1; - al = coar.x[i] * fine.x[j]; - bet = ((coar.x[i] * fine.x[j + 1] + coar.x[i + 1] * fine.x[j]) - + coar.x[i + 1] * fine.x[j + 1]); + al = coar.x[i] * fine.x[j]; + bet = ((coar.x[i] * fine.x[j + 1] + coar.x[i + 1] * fine.x[j]) + + coar.x[i + 1] * fine.x[j + 1]); - rem = (bet + bet * eps) + al * eps; - res = al + rem; - cor = (al - res) + rem; - if (res == (res + cor * err_0)) - { - retval = res * binexp.x; - goto ret; - } - else - { - retval = __slowexp (x); - goto ret; - } /*if error is over bound */ - } + rem = (bet + bet * eps) + al * eps; + res = al + rem; + cor = (al - res) + rem; + if (res == (res + cor * err_0)) + { + retval = res * binexp.x; + goto ret; + } + else + { + retval = __slowexp (x); + goto ret; + } /*if error is over bound */ + } - if (n <= smallint) - { - retval = 1.0; - goto ret; - } + if (n <= smallint) + { + retval = 1.0; + goto ret; + } - if (n >= badint) - { - if (n > infint) - { - retval = x + x; - goto ret; - } /* x is NaN */ - if (n < infint) - { - retval = (x > 0) ? (hhuge * hhuge) : (tiny * tiny); - goto ret; - } - /* x is finite, cause either overflow or underflow */ - if (junk1.i[LOW_HALF] != 0) - { - retval = x + x; - goto ret; - } /* x is NaN */ - retval = (x > 0) ? inf.x : zero; /* |x| = inf; return either inf or 0 */ - goto ret; - } + if (n >= badint) + { + if (n > infint) + { + retval = x + x; + goto ret; + } /* x is NaN */ + if (n < infint) + { + if (x > 0) + goto ret_huge; + else + goto ret_tiny; + } + /* x is finite, cause either overflow or underflow */ + if (junk1.i[LOW_HALF] != 0) + { + retval = x + x; + goto ret; + } /* x is NaN */ + retval = (x > 0) ? inf.x : zero; /* |x| = inf; return either inf or 0 */ + goto ret; + } - y = x * log2e.x + three51.x; - bexp = y - three51.x; - junk1.x = y; - eps = bexp * ln_two2.x; - t = x - bexp * ln_two1.x; - y = t + three33.x; - base = y - three33.x; - junk2.x = y; - del = (t - base) - eps; - eps = del + del * del * (p3.x * del + p2.x); - i = ((junk2.i[LOW_HALF] >> 8) & 0xfffffffe) + 356; - j = (junk2.i[LOW_HALF] & 511) << 1; - al = coar.x[i] * fine.x[j]; - bet = ((coar.x[i] * fine.x[j + 1] + coar.x[i + 1] * fine.x[j]) - + coar.x[i + 1] * fine.x[j + 1]); - rem = (bet + bet * eps) + al * eps; - res = al + rem; - cor = (al - res) + rem; - if (m >> 31) - { - ex = junk1.i[LOW_HALF]; - if (res < 1.0) - { - res += res; - cor += cor; - ex -= 1; - } - if (ex >= -1022) - { - binexp.i[HIGH_HALF] = (1023 + ex) << 20; - if (res == (res + cor * err_0)) - { - retval = res * binexp.x; - goto ret; - } - else - { - retval = __slowexp (x); - goto check_uflow_ret; - } /*if error is over bound */ - } - ex = -(1022 + ex); - binexp.i[HIGH_HALF] = (1023 - ex) << 20; - res *= binexp.x; - cor *= binexp.x; - eps = 1.0000000001 + err_0 * binexp.x; - t = 1.0 + res; - y = ((1.0 - t) + res) + cor; - res = t + y; - cor = (t - res) + y; - if (res == (res + eps * cor)) - { - binexp.i[HIGH_HALF] = 0x00100000; - retval = (res - 1.0) * binexp.x; - goto check_uflow_ret; - } - else - { - retval = __slowexp (x); - goto check_uflow_ret; - } /* if error is over bound */ - check_uflow_ret: - if (retval < DBL_MIN) - { + y = x * log2e.x + three51.x; + bexp = y - three51.x; + junk1.x = y; + eps = bexp * ln_two2.x; + t = x - bexp * ln_two1.x; + y = t + three33.x; + base = y - three33.x; + junk2.x = y; + del = (t - base) - eps; + eps = del + del * del * (p3.x * del + p2.x); + i = ((junk2.i[LOW_HALF] >> 8) & 0xfffffffe) + 356; + j = (junk2.i[LOW_HALF] & 511) << 1; + al = coar.x[i] * fine.x[j]; + bet = ((coar.x[i] * fine.x[j + 1] + coar.x[i + 1] * fine.x[j]) + + coar.x[i + 1] * fine.x[j + 1]); + rem = (bet + bet * eps) + al * eps; + res = al + rem; + cor = (al - res) + rem; + if (m >> 31) + { + ex = junk1.i[LOW_HALF]; + if (res < 1.0) + { + res += res; + cor += cor; + ex -= 1; + } + if (ex >= -1022) + { + binexp.i[HIGH_HALF] = (1023 + ex) << 20; + if (res == (res + cor * err_0)) + { + retval = res * binexp.x; + goto ret; + } + else + { + retval = __slowexp (x); + goto check_uflow_ret; + } /*if error is over bound */ + } + ex = -(1022 + ex); + binexp.i[HIGH_HALF] = (1023 - ex) << 20; + res *= binexp.x; + cor *= binexp.x; + eps = 1.0000000001 + err_0 * binexp.x; + t = 1.0 + res; + y = ((1.0 - t) + res) + cor; + res = t + y; + cor = (t - res) + y; + if (res == (res + eps * cor)) + { + binexp.i[HIGH_HALF] = 0x00100000; + retval = (res - 1.0) * binexp.x; + goto check_uflow_ret; + } + else + { + retval = __slowexp (x); + goto check_uflow_ret; + } /* if error is over bound */ + check_uflow_ret: + if (retval < DBL_MIN) + { #if FLT_EVAL_METHOD != 0 - volatile + volatile #endif - double force_underflow = tiny * tiny; - math_force_eval (force_underflow); - } - goto ret; - } - else - { - binexp.i[HIGH_HALF] = (junk1.i[LOW_HALF] + 767) << 20; - if (res == (res + cor * err_0)) - { + double force_underflow = tiny * tiny; + math_force_eval (force_underflow); + } + if (retval == 0) + goto ret_tiny; + goto ret; + } + else + { + binexp.i[HIGH_HALF] = (junk1.i[LOW_HALF] + 767) << 20; + if (res == (res + cor * err_0)) retval = res * binexp.x * t256.x; - goto ret; - } - else - { + else retval = __slowexp (x); + if (isinf (retval)) + goto ret_huge; + else goto ret; - } - } + } + } ret: return retval; + + ret_huge: + return hhuge * hhuge; + + ret_tiny: + return tiny * tiny; } #ifndef __ieee754_exp strong_alias (__ieee754_exp, __exp_finite) diff --git a/sysdeps/ieee754/dbl-64/e_exp10.c b/sysdeps/ieee754/dbl-64/e_exp10.c index eac609b71d..b14eaa9833 100644 --- a/sysdeps/ieee754/dbl-64/e_exp10.c +++ b/sysdeps/ieee754/dbl-64/e_exp10.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2012-2014 Free Software Foundation, Inc. +/* Copyright (C) 2012-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 @@ -29,12 +29,14 @@ __ieee754_exp10 (double arg) double arg_high, arg_low; double exp_high, exp_low; - if (!__finite (arg)) + if (!isfinite (arg)) return __ieee754_exp (arg); if (arg < DBL_MIN_10_EXP - DBL_DIG - 10) return DBL_MIN * DBL_MIN; else if (arg > DBL_MAX_10_EXP + 1) return DBL_MAX * DBL_MAX; + else if (fabs (arg) < 0x1p-56) + return 1.0; GET_LOW_WORD (lx, arg); lx &= 0xf8000000; diff --git a/sysdeps/ieee754/dbl-64/e_exp2.c b/sysdeps/ieee754/dbl-64/e_exp2.c index 10e23e2218..30f0a8f529 100644 --- a/sysdeps/ieee754/dbl-64/e_exp2.c +++ b/sysdeps/ieee754/dbl-64/e_exp2.c @@ -1,5 +1,5 @@ /* Double-precision floating point 2^x. - Copyright (C) 1997-2014 Free Software Foundation, Inc. + Copyright (C) 1997-2015 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Geoffrey Keating <geoffk@ozemail.com.au> @@ -43,12 +43,12 @@ __ieee754_exp2 (double x) static const double lomark = (double) (DBL_MIN_EXP - DBL_MANT_DIG - 1); /* Check for usual case. */ - if (__builtin_expect (isless (x, himark), 1)) + if (__glibc_likely (isless (x, himark))) { /* Exceptional cases: */ - if (__builtin_expect (!isgreaterequal (x, lomark), 0)) + if (__glibc_unlikely (!isgreaterequal (x, lomark))) { - if (__isinf (x)) + if (isinf (x)) /* e^-inf == 0, with no error. */ return 0; else @@ -61,6 +61,9 @@ __ieee754_exp2 (double x) double rx, x22, result; union ieee754_double ex2_u, scale_u; + if (fabs (x) < DBL_EPSILON / 4.0) + return 1.0 + x; + { SET_RESTORE_ROUND_NOEX (FE_TONEAREST); @@ -93,7 +96,9 @@ __ieee754_exp2 (double x) /* 3. Compute ex2 = 2^(t/512+e+ex). */ ex2_u.d = exp2_accuratetable[tval & 511]; tval >>= 9; - unsafe = abs (tval) >= -DBL_MIN_EXP - 1; + /* x2 is an integer multiple of 2^-54; avoid intermediate + underflow from the calculation of x22 * x. */ + unsafe = abs (tval) >= -DBL_MIN_EXP - 56; ex2_u.ieee.exponent += tval >> unsafe; scale_u.d = 1.0; scale_u.ieee.exponent += tval - (tval >> unsafe); diff --git a/sysdeps/ieee754/dbl-64/e_fmod.c b/sysdeps/ieee754/dbl-64/e_fmod.c index c83c2aedb2..e82b302200 100644 --- a/sysdeps/ieee754/dbl-64/e_fmod.c +++ b/sysdeps/ieee754/dbl-64/e_fmod.c @@ -45,7 +45,7 @@ __ieee754_fmod (double x, double y) } /* determine ix = ilogb(x) */ - if (__builtin_expect (hx < 0x00100000, 0)) /* subnormal x */ + if (__glibc_unlikely (hx < 0x00100000)) /* subnormal x */ { if (hx == 0) { @@ -62,7 +62,7 @@ __ieee754_fmod (double x, double y) ix = (hx >> 20) - 1023; /* determine iy = ilogb(y) */ - if (__builtin_expect (hy < 0x00100000, 0)) /* subnormal y */ + if (__glibc_unlikely (hy < 0x00100000)) /* subnormal y */ { if (hy == 0) { @@ -79,7 +79,7 @@ __ieee754_fmod (double x, double y) iy = (hy >> 20) - 1023; /* set up {hx,lx}, {hy,ly} and align y to x */ - if (__builtin_expect (ix >= -1022, 1)) + if (__glibc_likely (ix >= -1022)) hx = 0x00100000 | (0x000fffff & hx); else /* subnormal x, shift x to normal */ { @@ -95,7 +95,7 @@ __ieee754_fmod (double x, double y) lx = 0; } } - if (__builtin_expect (iy >= -1022, 1)) + if (__glibc_likely (iy >= -1022)) hy = 0x00100000 | (0x000fffff & hy); else /* subnormal y, shift y to normal */ { @@ -144,7 +144,7 @@ __ieee754_fmod (double x, double y) hx = hx + hx + (lx >> 31); lx = lx + lx; iy -= 1; } - if (__builtin_expect (iy >= -1022, 1)) /* normalize output */ + if (__glibc_likely (iy >= -1022)) /* normalize output */ { hx = ((hx - 0x00100000) | ((iy + 1023) << 20)); INSERT_WORDS (x, hx | sx, lx); diff --git a/sysdeps/ieee754/dbl-64/e_gamma_r.c b/sysdeps/ieee754/dbl-64/e_gamma_r.c index 1c427556bc..adeb61a248 100644 --- a/sysdeps/ieee754/dbl-64/e_gamma_r.c +++ b/sysdeps/ieee754/dbl-64/e_gamma_r.c @@ -1,5 +1,5 @@ /* Implementation of gamma function according to ISO C. - Copyright (C) 1997-2014 Free Software Foundation, Inc. + Copyright (C) 1997-2015 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -104,7 +104,7 @@ gamma_positive (double x, int *exp2_adj) * __ieee754_exp (-x_adj) * __ieee754_sqrt (2 * M_PI / x_adj) / prod); - exp_adj += x_eps * __ieee754_log (x); + exp_adj += x_eps * __ieee754_log (x_adj); double bsum = gamma_coeff[NCOEFF - 1]; double x_adj2 = x_adj * x_adj; for (size_t i = 1; i <= NCOEFF - 1; i++) @@ -119,10 +119,14 @@ __ieee754_gamma_r (double x, int *signgamp) { int32_t hx; u_int32_t lx; +#if FLT_EVAL_METHOD != 0 + volatile +#endif + double ret; EXTRACT_WORDS (hx, lx, x); - if (__builtin_expect (((hx & 0x7fffffff) | lx) == 0, 0)) + if (__glibc_unlikely (((hx & 0x7fffffff) | lx) == 0)) { /* Return value for x == 0 is Inf with divide by zero exception. */ *signgamp = 0; @@ -135,13 +139,13 @@ __ieee754_gamma_r (double x, int *signgamp) *signgamp = 0; return (x - x) / (x - x); } - if (__builtin_expect ((unsigned int) hx == 0xfff00000 && lx == 0, 0)) + if (__glibc_unlikely ((unsigned int) hx == 0xfff00000 && lx == 0)) { /* x == -Inf. According to ISO this is NaN. */ *signgamp = 0; return x - x; } - if (__builtin_expect ((hx & 0x7ff00000) == 0x7ff00000, 0)) + if (__glibc_unlikely ((hx & 0x7ff00000) == 0x7ff00000)) { /* Positive infinity (return positive infinity) or NaN (return NaN). */ @@ -153,36 +157,69 @@ __ieee754_gamma_r (double x, int *signgamp) { /* Overflow. */ *signgamp = 0; - return DBL_MAX * DBL_MAX; + ret = DBL_MAX * DBL_MAX; + return ret; } - else if (x > 0.0) + else { - *signgamp = 0; - int exp2_adj; - double ret = gamma_positive (x, &exp2_adj); - return __scalbn (ret, exp2_adj); + SET_RESTORE_ROUND (FE_TONEAREST); + if (x > 0.0) + { + *signgamp = 0; + int exp2_adj; + double tret = gamma_positive (x, &exp2_adj); + ret = __scalbn (tret, exp2_adj); + } + else if (x >= -DBL_EPSILON / 4.0) + { + *signgamp = 0; + ret = 1.0 / x; + } + else + { + double tx = __trunc (x); + *signgamp = (tx == 2.0 * __trunc (tx / 2.0)) ? -1 : 1; + if (x <= -184.0) + /* Underflow. */ + ret = DBL_MIN * DBL_MIN; + else + { + double frac = tx - x; + if (frac > 0.5) + frac = 1.0 - frac; + double sinpix = (frac <= 0.25 + ? __sin (M_PI * frac) + : __cos (M_PI * (0.5 - frac))); + int exp2_adj; + double tret = M_PI / (-x * sinpix + * gamma_positive (-x, &exp2_adj)); + ret = __scalbn (tret, -exp2_adj); + } + } } - else if (x >= -DBL_EPSILON / 4.0) + if (isinf (ret) && x != 0) { - *signgamp = 0; - return 1.0 / x; + if (*signgamp < 0) + { + ret = -__copysign (DBL_MAX, ret) * DBL_MAX; + ret = -ret; + } + else + ret = __copysign (DBL_MAX, ret) * DBL_MAX; + return ret; } - else + else if (ret == 0) { - double tx = __trunc (x); - *signgamp = (tx == 2.0 * __trunc (tx / 2.0)) ? -1 : 1; - if (x <= -184.0) - /* Underflow. */ - return DBL_MIN * DBL_MIN; - double frac = tx - x; - if (frac > 0.5) - frac = 1.0 - frac; - double sinpix = (frac <= 0.25 - ? __sin (M_PI * frac) - : __cos (M_PI * (0.5 - frac))); - int exp2_adj; - double ret = M_PI / (-x * sinpix * gamma_positive (-x, &exp2_adj)); - return __scalbn (ret, -exp2_adj); + if (*signgamp < 0) + { + ret = -__copysign (DBL_MIN, ret) * DBL_MIN; + ret = -ret; + } + else + ret = __copysign (DBL_MIN, ret) * DBL_MIN; + return ret; } + else + return ret; } strong_alias (__ieee754_gamma_r, __gamma_r_finite) diff --git a/sysdeps/ieee754/dbl-64/e_hypot.c b/sysdeps/ieee754/dbl-64/e_hypot.c index 88242bc4f6..5cbfcbeb48 100644 --- a/sysdeps/ieee754/dbl-64/e_hypot.c +++ b/sysdeps/ieee754/dbl-64/e_hypot.c @@ -70,7 +70,7 @@ __ieee754_hypot (double x, double y) return a + b; } /* x/y > 2**60 */ k = 0; - if (__builtin_expect (ha > 0x5f300000, 0)) /* a>2**500 */ + if (__glibc_unlikely (ha > 0x5f300000)) /* a>2**500 */ { if (ha >= 0x7ff00000) /* Inf or NaN */ { diff --git a/sysdeps/ieee754/dbl-64/e_j0.c b/sysdeps/ieee754/dbl-64/e_j0.c index d165e80925..9f25aa855e 100644 --- a/sysdeps/ieee754/dbl-64/e_j0.c +++ b/sysdeps/ieee754/dbl-64/e_j0.c @@ -305,6 +305,7 @@ pzero (double x) int32_t ix; GET_HIGH_WORD (ix, x); ix &= 0x7fffffff; + /* ix >= 0x40000000 for all calls to this function. */ if (ix >= 0x41b00000) { return one; @@ -321,7 +322,7 @@ pzero (double x) { p = pR3; q = pS3; } - else if (ix >= 0x40000000) + else { p = pR2; q = pS2; } @@ -423,6 +424,7 @@ qzero (double x) int32_t ix; GET_HIGH_WORD (ix, x); ix &= 0x7fffffff; + /* ix >= 0x40000000 for all calls to this function. */ if (ix >= 0x41b00000) { return -.125 / x; @@ -439,7 +441,7 @@ qzero (double x) { p = qR3; q = qS3; } - else if (ix >= 0x40000000) + else { p = qR2; q = qS2; } diff --git a/sysdeps/ieee754/dbl-64/e_j1.c b/sysdeps/ieee754/dbl-64/e_j1.c index ab754c6ee0..26ffdfe282 100644 --- a/sysdeps/ieee754/dbl-64/e_j1.c +++ b/sysdeps/ieee754/dbl-64/e_j1.c @@ -58,6 +58,8 @@ * by method mentioned above. */ +#include <errno.h> +#include <float.h> #include <math.h> #include <math_private.h> @@ -89,7 +91,7 @@ __ieee754_j1 (double x) GET_HIGH_WORD (hx, x); ix = hx & 0x7fffffff; - if (__builtin_expect (ix >= 0x7ff00000, 0)) + if (__glibc_unlikely (ix >= 0x7ff00000)) return one / x; y = fabs (x); if (ix >= 0x40000000) /* |x| >= 2.0 */ @@ -121,10 +123,18 @@ __ieee754_j1 (double x) else return z; } - if (__builtin_expect (ix < 0x3e400000, 0)) /* |x|<2**-27 */ + if (__glibc_unlikely (ix < 0x3e400000)) /* |x|<2**-27 */ { - if (huge + x > one) - return 0.5 * x; /* inexact if x!=0 necessary */ + if (huge + x > one) /* inexact if x!=0 necessary */ + { + double ret = 0.5 * x; + if (fabs (ret) < DBL_MIN) + { + double force_underflow = ret * ret; + math_force_eval (force_underflow); + } + return ret; + } } z = x * x; r1 = z * R[0]; z2 = z * z; @@ -163,12 +173,12 @@ __ieee754_y1 (double x) EXTRACT_WORDS (hx, lx, x); ix = 0x7fffffff & hx; /* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */ - if (__builtin_expect (ix >= 0x7ff00000, 0)) + if (__glibc_unlikely (ix >= 0x7ff00000)) return one / (x + x * x); - if (__builtin_expect ((ix | lx) == 0, 0)) + if (__glibc_unlikely ((ix | lx) == 0)) return -HUGE_VAL + x; /* -inf and overflow exception. */; - if (__builtin_expect (hx < 0, 0)) + if (__glibc_unlikely (hx < 0)) return zero / (zero * x); if (ix >= 0x40000000) /* |x| >= 2.0 */ { @@ -203,9 +213,12 @@ __ieee754_y1 (double x) } return z; } - if (__builtin_expect (ix <= 0x3c900000, 0)) /* x < 2**-54 */ + if (__glibc_unlikely (ix <= 0x3c900000)) /* x < 2**-54 */ { - return (-tpi / x); + z = -tpi / x; + if (isinf (z)) + __set_errno (ERANGE); + return z; } z = x * x; u1 = U0[0] + z * U0[1]; z2 = z * z; @@ -301,6 +314,7 @@ pone (double x) int32_t ix; GET_HIGH_WORD (ix, x); ix &= 0x7fffffff; + /* ix >= 0x40000000 for all calls to this function. */ if (ix >= 0x41b00000) { return one; @@ -317,7 +331,7 @@ pone (double x) { p = pr3; q = ps3; } - else if (ix >= 0x40000000) + else { p = pr2; q = ps2; } @@ -420,6 +434,7 @@ qone (double x) int32_t ix; GET_HIGH_WORD (ix, x); ix &= 0x7fffffff; + /* ix >= 0x40000000 for all calls to this function. */ if (ix >= 0x41b00000) { return .375 / x; @@ -436,7 +451,7 @@ qone (double x) { p = qr3; q = qs3; } - else if (ix >= 0x40000000) + else { p = qr2; q = qs2; } diff --git a/sysdeps/ieee754/dbl-64/e_jn.c b/sysdeps/ieee754/dbl-64/e_jn.c index f48e43a0d9..ccef2dcd80 100644 --- a/sysdeps/ieee754/dbl-64/e_jn.c +++ b/sysdeps/ieee754/dbl-64/e_jn.c @@ -37,6 +37,7 @@ */ #include <errno.h> +#include <float.h> #include <math.h> #include <math_private.h> @@ -51,7 +52,7 @@ double __ieee754_jn (int n, double x) { int32_t i, hx, ix, lx, sgn; - double a, b, temp, di; + double a, b, temp, di, ret; double z, w; /* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x) @@ -60,7 +61,7 @@ __ieee754_jn (int n, double x) EXTRACT_WORDS (hx, lx, x); ix = 0x7fffffff & hx; /* if J(n,NaN) is NaN */ - if (__builtin_expect ((ix | ((u_int32_t) (lx | -lx)) >> 31) > 0x7ff00000, 0)) + if (__glibc_unlikely ((ix | ((u_int32_t) (lx | -lx)) >> 31) > 0x7ff00000)) return x + x; if (n < 0) { @@ -74,14 +75,16 @@ __ieee754_jn (int n, double x) return (__ieee754_j1 (x)); sgn = (n & 1) & (hx >> 31); /* even n -- 0, odd n -- sign(x) */ x = fabs (x); - if (__builtin_expect ((ix | lx) == 0 || ix >= 0x7ff00000, 0)) - /* if x is 0 or inf */ - b = zero; - else if ((double) n <= x) - { - /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */ - if (ix >= 0x52D00000) /* x > 2**302 */ - { /* (x >> n**2) + { + SET_RESTORE_ROUND (FE_TONEAREST); + if (__glibc_unlikely ((ix | lx) == 0 || ix >= 0x7ff00000)) + /* if x is 0 or inf */ + return sgn == 1 ? -zero : zero; + else if ((double) n <= x) + { + /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */ + if (ix >= 0x52D00000) /* x > 2**302 */ + { /* (x >> n**2) * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi) * Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi) * Let s=sin(x), c=cos(x), @@ -94,152 +97,161 @@ __ieee754_jn (int n, double x) * 2 -s+c -c-s * 3 s+c c-s */ - double s; - double c; - __sincos (x, &s, &c); - switch (n & 3) - { - case 0: temp = c + s; break; - case 1: temp = -c + s; break; - case 2: temp = -c - s; break; - case 3: temp = c - s; break; - } - b = invsqrtpi * temp / __ieee754_sqrt (x); - } - else - { - a = __ieee754_j0 (x); - b = __ieee754_j1 (x); - for (i = 1; i < n; i++) - { - temp = b; - b = b * ((double) (i + i) / x) - a; /* avoid underflow */ - a = temp; - } - } - } - else - { - if (ix < 0x3e100000) /* x < 2**-29 */ - { /* x is tiny, return the first Taylor expansion of J(n,x) + double s; + double c; + __sincos (x, &s, &c); + switch (n & 3) + { + case 0: temp = c + s; break; + case 1: temp = -c + s; break; + case 2: temp = -c - s; break; + case 3: temp = c - s; break; + } + b = invsqrtpi * temp / __ieee754_sqrt (x); + } + else + { + a = __ieee754_j0 (x); + b = __ieee754_j1 (x); + for (i = 1; i < n; i++) + { + temp = b; + b = b * ((double) (i + i) / x) - a; /* avoid underflow */ + a = temp; + } + } + } + else + { + if (ix < 0x3e100000) /* x < 2**-29 */ + { /* x is tiny, return the first Taylor expansion of J(n,x) * J(n,x) = 1/n!*(x/2)^n - ... */ - if (n > 33) /* underflow */ - b = zero; - else - { - temp = x * 0.5; b = temp; - for (a = one, i = 2; i <= n; i++) - { - a *= (double) i; /* a = n! */ - b *= temp; /* b = (x/2)^n */ - } - b = b / a; - } - } - else - { - /* use backward recurrence */ - /* x x^2 x^2 - * J(n,x)/J(n-1,x) = ---- ------ ------ ..... - * 2n - 2(n+1) - 2(n+2) - * - * 1 1 1 - * (for large x) = ---- ------ ------ ..... - * 2n 2(n+1) 2(n+2) - * -- - ------ - ------ - - * x x x - * - * Let w = 2n/x and h=2/x, then the above quotient - * is equal to the continued fraction: - * 1 - * = ----------------------- - * 1 - * w - ----------------- - * 1 - * w+h - --------- - * w+2h - ... - * - * To determine how many terms needed, let - * Q(0) = w, Q(1) = w(w+h) - 1, - * Q(k) = (w+k*h)*Q(k-1) - Q(k-2), - * When Q(k) > 1e4 good for single - * When Q(k) > 1e9 good for double - * When Q(k) > 1e17 good for quadruple - */ - /* determine k */ - double t, v; - double q0, q1, h, tmp; int32_t k, m; - w = (n + n) / (double) x; h = 2.0 / (double) x; - q0 = w; z = w + h; q1 = w * z - 1.0; k = 1; - while (q1 < 1.0e9) - { - k += 1; z += h; - tmp = z * q1 - q0; - q0 = q1; - q1 = tmp; - } - m = n + n; - for (t = zero, i = 2 * (n + k); i >= m; i -= 2) - t = one / (i / x - t); - a = t; - b = one; - /* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n) - * Hence, if n*(log(2n/x)) > ... - * single 8.8722839355e+01 - * double 7.09782712893383973096e+02 - * long double 1.1356523406294143949491931077970765006170e+04 - * then recurrent value may overflow and the result is - * likely underflow to zero - */ - tmp = n; - v = two / x; - tmp = tmp * __ieee754_log (fabs (v * tmp)); - if (tmp < 7.09782712893383973096e+02) - { - for (i = n - 1, di = (double) (i + i); i > 0; i--) - { - temp = b; - b *= di; - b = b / x - a; - a = temp; - di -= two; - } - } - else - { - for (i = n - 1, di = (double) (i + i); i > 0; i--) - { - temp = b; - b *= di; - b = b / x - a; - a = temp; - di -= two; - /* scale b to avoid spurious overflow */ - if (b > 1e100) - { - a /= b; - t /= b; - b = one; - } - } - } - /* j0() and j1() suffer enormous loss of precision at and - * near zero; however, we know that their zero points never - * coincide, so just choose the one further away from zero. - */ - z = __ieee754_j0 (x); - w = __ieee754_j1 (x); - if (fabs (z) >= fabs (w)) - b = (t * z / b); - else - b = (t * w / a); - } + if (n > 33) /* underflow */ + b = zero; + else + { + temp = x * 0.5; b = temp; + for (a = one, i = 2; i <= n; i++) + { + a *= (double) i; /* a = n! */ + b *= temp; /* b = (x/2)^n */ + } + b = b / a; + } + } + else + { + /* use backward recurrence */ + /* x x^2 x^2 + * J(n,x)/J(n-1,x) = ---- ------ ------ ..... + * 2n - 2(n+1) - 2(n+2) + * + * 1 1 1 + * (for large x) = ---- ------ ------ ..... + * 2n 2(n+1) 2(n+2) + * -- - ------ - ------ - + * x x x + * + * Let w = 2n/x and h=2/x, then the above quotient + * is equal to the continued fraction: + * 1 + * = ----------------------- + * 1 + * w - ----------------- + * 1 + * w+h - --------- + * w+2h - ... + * + * To determine how many terms needed, let + * Q(0) = w, Q(1) = w(w+h) - 1, + * Q(k) = (w+k*h)*Q(k-1) - Q(k-2), + * When Q(k) > 1e4 good for single + * When Q(k) > 1e9 good for double + * When Q(k) > 1e17 good for quadruple + */ + /* determine k */ + double t, v; + double q0, q1, h, tmp; int32_t k, m; + w = (n + n) / (double) x; h = 2.0 / (double) x; + q0 = w; z = w + h; q1 = w * z - 1.0; k = 1; + while (q1 < 1.0e9) + { + k += 1; z += h; + tmp = z * q1 - q0; + q0 = q1; + q1 = tmp; + } + m = n + n; + for (t = zero, i = 2 * (n + k); i >= m; i -= 2) + t = one / (i / x - t); + a = t; + b = one; + /* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n) + * Hence, if n*(log(2n/x)) > ... + * single 8.8722839355e+01 + * double 7.09782712893383973096e+02 + * long double 1.1356523406294143949491931077970765006170e+04 + * then recurrent value may overflow and the result is + * likely underflow to zero + */ + tmp = n; + v = two / x; + tmp = tmp * __ieee754_log (fabs (v * tmp)); + if (tmp < 7.09782712893383973096e+02) + { + for (i = n - 1, di = (double) (i + i); i > 0; i--) + { + temp = b; + b *= di; + b = b / x - a; + a = temp; + di -= two; + } + } + else + { + for (i = n - 1, di = (double) (i + i); i > 0; i--) + { + temp = b; + b *= di; + b = b / x - a; + a = temp; + di -= two; + /* scale b to avoid spurious overflow */ + if (b > 1e100) + { + a /= b; + t /= b; + b = one; + } + } + } + /* j0() and j1() suffer enormous loss of precision at and + * near zero; however, we know that their zero points never + * coincide, so just choose the one further away from zero. + */ + z = __ieee754_j0 (x); + w = __ieee754_j1 (x); + if (fabs (z) >= fabs (w)) + b = (t * z / b); + else + b = (t * w / a); + } + } + if (sgn == 1) + ret = -b; + else + ret = b; + } + if (ret == 0) + ret = __copysign (DBL_MIN, ret) * DBL_MIN; + else if (fabs (ret) < DBL_MIN) + { + double force_underflow = ret * ret; + math_force_eval (force_underflow); } - if (sgn == 1) - return -b; - else - return b; + return ret; } strong_alias (__ieee754_jn, __jn_finite) @@ -248,17 +260,17 @@ __ieee754_yn (int n, double x) { int32_t i, hx, ix, lx; int32_t sign; - double a, b, temp; + double a, b, temp, ret; EXTRACT_WORDS (hx, lx, x); ix = 0x7fffffff & hx; /* if Y(n,NaN) is NaN */ - if (__builtin_expect ((ix | ((u_int32_t) (lx | -lx)) >> 31) > 0x7ff00000, 0)) + if (__glibc_unlikely ((ix | ((u_int32_t) (lx | -lx)) >> 31) > 0x7ff00000)) return x + x; - if (__builtin_expect ((ix | lx) == 0, 0)) + if (__glibc_unlikely ((ix | lx) == 0)) return -HUGE_VAL + x; /* -inf and overflow exception. */; - if (__builtin_expect (hx < 0, 0)) + if (__glibc_unlikely (hx < 0)) return zero / (zero * x); sign = 1; if (n < 0) @@ -268,57 +280,67 @@ __ieee754_yn (int n, double x) } if (n == 0) return (__ieee754_y0 (x)); - if (n == 1) - return (sign * __ieee754_y1 (x)); - if (__builtin_expect (ix == 0x7ff00000, 0)) - return zero; - if (ix >= 0x52D00000) /* x > 2**302 */ - { /* (x >> n**2) - * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi) - * Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi) - * Let s=sin(x), c=cos(x), - * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then - * - * n sin(xn)*sqt2 cos(xn)*sqt2 - * ---------------------------------- - * 0 s-c c+s - * 1 -s-c -c+s - * 2 -s+c -c-s - * 3 s+c c-s - */ - double c; - double s; - __sincos (x, &s, &c); - switch (n & 3) - { - case 0: temp = s - c; break; - case 1: temp = -s - c; break; - case 2: temp = -s + c; break; - case 3: temp = s + c; break; - } - b = invsqrtpi * temp / __ieee754_sqrt (x); - } - else - { - u_int32_t high; - a = __ieee754_y0 (x); - b = __ieee754_y1 (x); - /* quit if b is -inf */ - GET_HIGH_WORD (high, b); - for (i = 1; i < n && high != 0xfff00000; i++) - { - temp = b; - b = ((double) (i + i) / x) * b - a; - GET_HIGH_WORD (high, b); - a = temp; - } - /* If B is +-Inf, set up errno accordingly. */ - if (!__finite (b)) - __set_errno (ERANGE); - } - if (sign > 0) - return b; - else - return -b; + { + SET_RESTORE_ROUND (FE_TONEAREST); + if (n == 1) + { + ret = sign * __ieee754_y1 (x); + goto out; + } + if (__glibc_unlikely (ix == 0x7ff00000)) + return zero; + if (ix >= 0x52D00000) /* x > 2**302 */ + { /* (x >> n**2) + * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi) + * Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi) + * Let s=sin(x), c=cos(x), + * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then + * + * n sin(xn)*sqt2 cos(xn)*sqt2 + * ---------------------------------- + * 0 s-c c+s + * 1 -s-c -c+s + * 2 -s+c -c-s + * 3 s+c c-s + */ + double c; + double s; + __sincos (x, &s, &c); + switch (n & 3) + { + case 0: temp = s - c; break; + case 1: temp = -s - c; break; + case 2: temp = -s + c; break; + case 3: temp = s + c; break; + } + b = invsqrtpi * temp / __ieee754_sqrt (x); + } + else + { + u_int32_t high; + a = __ieee754_y0 (x); + b = __ieee754_y1 (x); + /* quit if b is -inf */ + GET_HIGH_WORD (high, b); + for (i = 1; i < n && high != 0xfff00000; i++) + { + temp = b; + b = ((double) (i + i) / x) * b - a; + GET_HIGH_WORD (high, b); + a = temp; + } + /* If B is +-Inf, set up errno accordingly. */ + if (!isfinite (b)) + __set_errno (ERANGE); + } + if (sign > 0) + ret = b; + else + ret = -b; + } + out: + if (isinf (ret)) + ret = __copysign (DBL_MAX, ret) * DBL_MAX; + return ret; } strong_alias (__ieee754_yn, __yn_finite) diff --git a/sysdeps/ieee754/dbl-64/e_lgamma_r.c b/sysdeps/ieee754/dbl-64/e_lgamma_r.c index af7d06c82f..fc6f594d62 100644 --- a/sysdeps/ieee754/dbl-64/e_lgamma_r.c +++ b/sysdeps/ieee754/dbl-64/e_lgamma_r.c @@ -77,6 +77,7 @@ * */ +#include <libc-internal.h> #include <math.h> #include <math_private.h> @@ -294,7 +295,18 @@ __ieee754_lgamma_r(double x, int *signgamp) } else /* 2**58 <= x <= inf */ r = x*(__ieee754_log(x)-one); + /* NADJ is set for negative arguments but not otherwise, + resulting in warnings that it may be used uninitialized + although in the cases where it is used it has always been + set. */ + DIAG_PUSH_NEEDS_COMMENT; +#if __GNUC_PREREQ (4, 7) + DIAG_IGNORE_NEEDS_COMMENT (4.9, "-Wmaybe-uninitialized"); +#else + DIAG_IGNORE_NEEDS_COMMENT (4.9, "-Wuninitialized"); +#endif if(hx<0) r = nadj - r; + DIAG_POP_NEEDS_COMMENT; return r; } strong_alias (__ieee754_lgamma_r, __lgamma_r_finite) diff --git a/sysdeps/ieee754/dbl-64/e_log.c b/sysdeps/ieee754/dbl-64/e_log.c index 0b2889cb3b..a9117fb4c5 100644 --- a/sysdeps/ieee754/dbl-64/e_log.c +++ b/sysdeps/ieee754/dbl-64/e_log.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001-2014 Free Software Foundation, Inc. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -38,6 +38,7 @@ #include <dla.h> #include "mpa.h" #include "MathLib.h" +#include <math.h> #include <math_private.h> #include <stap-probe.h> @@ -77,25 +78,29 @@ __ieee754_log (double x) ux = num.i[HIGH_HALF]; dx = num.i[LOW_HALF]; n = 0; - if (__builtin_expect (ux < 0x00100000, 0)) + if (__glibc_unlikely (ux < 0x00100000)) { - if (__builtin_expect (((ux & 0x7fffffff) | dx) == 0, 0)) + if (__glibc_unlikely (((ux & 0x7fffffff) | dx) == 0)) return MHALF / 0.0; /* return -INF */ - if (__builtin_expect (ux < 0, 0)) + if (__glibc_unlikely (ux < 0)) return (x - x) / 0.0; /* return NaN */ n -= 54; x *= two54.d; /* scale x */ num.d = x; } - if (__builtin_expect (ux >= 0x7ff00000, 0)) + if (__glibc_unlikely (ux >= 0x7ff00000)) return x + x; /* INF or NaN */ /* Regular values of x */ w = x - 1; - if (__builtin_expect (ABS (w) > U03, 1)) + if (__glibc_likely (fabs (w) > U03)) goto case_03; + /* log (1) is +0 in all rounding modes. */ + if (w == 0.0) + return 0.0; + /*--- Stage I, the case abs(x-1) < 0.03 */ t8 = MHALF * w; diff --git a/sysdeps/ieee754/dbl-64/e_log10.c b/sysdeps/ieee754/dbl-64/e_log10.c index c3d465a4a9..8548ee3942 100644 --- a/sysdeps/ieee754/dbl-64/e_log10.c +++ b/sysdeps/ieee754/dbl-64/e_log10.c @@ -63,15 +63,15 @@ __ieee754_log10 (double x) k = 0; if (hx < 0x00100000) { /* x < 2**-1022 */ - if (__builtin_expect (((hx & 0x7fffffff) | lx) == 0, 0)) + if (__glibc_unlikely (((hx & 0x7fffffff) | lx) == 0)) return -two54 / (x - x); /* log(+-0)=-inf */ - if (__builtin_expect (hx < 0, 0)) + if (__glibc_unlikely (hx < 0)) return (x - x) / (x - x); /* log(-#) = NaN */ k -= 54; x *= two54; /* subnormal number, scale up x */ GET_HIGH_WORD (hx, x); } - if (__builtin_expect (hx >= 0x7ff00000, 0)) + if (__glibc_unlikely (hx >= 0x7ff00000)) return x + x; k += (hx >> 20) - 1023; i = ((u_int32_t) k & 0x80000000) >> 31; diff --git a/sysdeps/ieee754/dbl-64/e_log2.c b/sysdeps/ieee754/dbl-64/e_log2.c index 890a4a2bd0..997d7cefc8 100644 --- a/sysdeps/ieee754/dbl-64/e_log2.c +++ b/sysdeps/ieee754/dbl-64/e_log2.c @@ -81,15 +81,15 @@ __ieee754_log2 (double x) k = 0; if (hx < 0x00100000) { /* x < 2**-1022 */ - if (__builtin_expect (((hx & 0x7fffffff) | lx) == 0, 0)) + if (__glibc_unlikely (((hx & 0x7fffffff) | lx) == 0)) return -two54 / (x - x); /* log(+-0)=-inf */ - if (__builtin_expect (hx < 0, 0)) + if (__glibc_unlikely (hx < 0)) return (x - x) / (x - x); /* log(-#) = NaN */ k -= 54; x *= two54; /* subnormal number, scale up x */ GET_HIGH_WORD (hx, x); } - if (__builtin_expect (hx >= 0x7ff00000, 0)) + if (__glibc_unlikely (hx >= 0x7ff00000)) return x + x; k += (hx >> 20) - 1023; hx &= 0x000fffff; diff --git a/sysdeps/ieee754/dbl-64/e_pow.c b/sysdeps/ieee754/dbl-64/e_pow.c index 1c5f4b311e..3c027fe7f5 100644 --- a/sysdeps/ieee754/dbl-64/e_pow.c +++ b/sysdeps/ieee754/dbl-64/e_pow.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001-2014 Free Software Foundation, Inc. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -34,6 +34,7 @@ /* round to nearest mode of IEEE 754 standard. */ /* */ /***************************************************************************/ +#include <math.h> #include "endian.h" #include "upow.h" #include <dla.h> @@ -91,27 +92,33 @@ __ieee754_pow (double x, double y) { /* if y<-1 or y>1 */ double retval; - SET_RESTORE_ROUND (FE_TONEAREST); + { + SET_RESTORE_ROUND (FE_TONEAREST); - /* Avoid internal underflow for tiny y. The exact value of y does - not matter if |y| <= 2**-64. */ - if (ABS (y) < 0x1p-64) - y = y < 0 ? -0x1p-64 : 0x1p-64; - z = log1 (x, &aa, &error); /* x^y =e^(y log (X)) */ - t = y * CN; - y1 = t - (t - y); - y2 = y - y1; - t = z * CN; - a1 = t - (t - z); - a2 = (z - a1) + aa; - a = y1 * a1; - aa = y2 * a1 + y * a2; - a1 = a + aa; - a2 = (a - a1) + aa; - error = error * ABS (y); - t = __exp1 (a1, a2, 1.9e16 * error); /* return -10 or 0 if wasn't computed exactly */ - retval = (t > 0) ? t : power1 (x, y); + /* Avoid internal underflow for tiny y. The exact value of y does + not matter if |y| <= 2**-64. */ + if (fabs (y) < 0x1p-64) + y = y < 0 ? -0x1p-64 : 0x1p-64; + z = log1 (x, &aa, &error); /* x^y =e^(y log (X)) */ + t = y * CN; + y1 = t - (t - y); + y2 = y - y1; + t = z * CN; + a1 = t - (t - z); + a2 = (z - a1) + aa; + a = y1 * a1; + aa = y2 * a1 + y * a2; + a1 = a + aa; + a2 = (a - a1) + aa; + error = error * fabs (y); + t = __exp1 (a1, a2, 1.9e16 * error); /* return -10 or 0 if wasn't computed exactly */ + retval = (t > 0) ? t : power1 (x, y); + } + if (isinf (retval)) + retval = huge * huge; + else if (retval == 0) + retval = tiny * tiny; return retval; } @@ -120,7 +127,7 @@ __ieee754_pow (double x, double y) if (((v.i[HIGH_HALF] & 0x7fffffff) == 0x7ff00000 && v.i[LOW_HALF] != 0) || (v.i[HIGH_HALF] & 0x7fffffff) > 0x7ff00000) /* NaN */ return y; - if (ABS (y) > 1.0e20) + if (fabs (y) > 1.0e20) return (y > 0) ? 0 : 1.0 / 0.0; k = checkint (y); if (k == -1) @@ -164,7 +171,21 @@ __ieee754_pow (double x, double y) return y < 0 ? 0.0 : INF.x; } /* if y even or odd */ - return (k == 1) ? __ieee754_pow (-x, y) : -__ieee754_pow (-x, y); + if (k == 1) + return __ieee754_pow (-x, y); + else + { + double retval; + { + SET_RESTORE_ROUND (FE_TONEAREST); + retval = -__ieee754_pow (-x, y); + } + if (isinf (retval)) + retval = -huge * huge; + else if (retval == 0) + retval = -tiny * tiny; + return retval; + } } /* x>0 */ @@ -211,7 +232,7 @@ power1 (double x, double y) aa = ((y1 * a1 - a) + y1 * a2 + y2 * a1) + y2 * a2 + aa * y; a1 = a + aa; a2 = (a - a1) + aa; - error = error * ABS (y); + error = error * fabs (y); t = __exp1 (a1, a2, 1.9e16 * error); return (t >= 0) ? t : __slowpow (x, y, z); } @@ -271,7 +292,7 @@ log1 (double x, double *delta, double *error) * (r7 + t * r8))))) - 0.5 * t2 * (t + t1)); res = e1 + e2; - *error = 1.0e-21 * ABS (t); + *error = 1.0e-21 * fabs (t); *delta = (e1 - res) + e2; return res; } /* |x-1| < 1.5*2**-10 */ @@ -377,7 +398,7 @@ my_log2 (double x, double *delta, double *error) e2 = ((((t - e1) + z) + zz) + t * t * t * (ss3 + t * (s4 + t * (s5 + t * (s6 + t * (s7 + t * s8)))))); res = e1 + e2; - *error = 1.0e-25 * ABS (t); + *error = 1.0e-25 * fabs (t); *delta = (e1 - res) + e2; return res; } diff --git a/sysdeps/ieee754/dbl-64/e_remainder.c b/sysdeps/ieee754/dbl-64/e_remainder.c index 28689614f0..7f3a02b3c2 100644 --- a/sysdeps/ieee754/dbl-64/e_remainder.c +++ b/sysdeps/ieee754/dbl-64/e_remainder.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001-2014 Free Software Foundation, Inc. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -33,6 +33,7 @@ #include "mydefs.h" #include "urem.h" #include "MathLib.h" +#include <math.h> #include <math_private.h> /**************************************************************************/ @@ -66,7 +67,7 @@ __ieee754_remainder (double x, double y) return (xx != 0) ? xx : ((x > 0) ? ZERO.x : nZERO.x); else { - if (ABS (xx) > 0.5 * t.x) + if (fabs (xx) > 0.5 * t.x) return (z > d) ? xx - t.x : xx + t.x; else return xx; @@ -98,10 +99,10 @@ __ieee754_remainder (double x, double y) z = u.x * r.x; d = (z + big.x) - big.x; u.x = (u.x - d * w.x) - d * ww.x; - if (ABS (u.x) < 0.5 * t.x) + if (fabs (u.x) < 0.5 * t.x) return (u.x != 0) ? u.x : ((x > 0) ? ZERO.x : nZERO.x); else - if (ABS (u.x) > 0.5 * t.x) + if (fabs (u.x) > 0.5 * t.x) return (d > z) ? u.x + t.x : u.x - t.x; else { @@ -114,7 +115,7 @@ __ieee754_remainder (double x, double y) { if (kx < 0x7fe00000 && ky < 0x7ff00000 && (ky > 0 || t.i[LOW_HALF] != 0)) { - y = ABS (y) * t128.x; + y = fabs (y) * t128.x; z = __ieee754_remainder (x, y) * t128.x; z = __ieee754_remainder (z, y) * tm128.x; return z; @@ -124,10 +125,10 @@ __ieee754_remainder (double x, double y) if ((kx & 0x7ff00000) == 0x7fe00000 && ky < 0x7ff00000 && (ky > 0 || t.i[LOW_HALF] != 0)) { - y = ABS (y); + y = fabs (y); z = 2.0 * __ieee754_remainder (0.5 * x, y); - d = ABS (z); - if (d <= ABS (d - y)) + d = fabs (z); + if (d <= fabs (d - y)) return z; else return (z > 0) ? z - y : z + y; diff --git a/sysdeps/ieee754/dbl-64/e_sinh.c b/sysdeps/ieee754/dbl-64/e_sinh.c index 851b510aaa..4ff28bf85d 100644 --- a/sysdeps/ieee754/dbl-64/e_sinh.c +++ b/sysdeps/ieee754/dbl-64/e_sinh.c @@ -49,7 +49,7 @@ __ieee754_sinh (double x) ix = jx & 0x7fffffff; /* x is INF or NaN */ - if (__builtin_expect (ix >= 0x7ff00000, 0)) + if (__glibc_unlikely (ix >= 0x7ff00000)) return x + x; h = 0.5; @@ -58,7 +58,7 @@ __ieee754_sinh (double x) /* |x| in [0,22], return sign(x)*0.5*(E+E/(E+1))) */ if (ix < 0x40360000) /* |x|<22 */ { - if (__builtin_expect (ix < 0x3e300000, 0)) /* |x|<2**-28 */ + if (__glibc_unlikely (ix < 0x3e300000)) /* |x|<2**-28 */ if (shuge + x > one) return x; /* sinh(tiny) = tiny with inexact */ diff --git a/sysdeps/ieee754/dbl-64/e_sqrt.c b/sysdeps/ieee754/dbl-64/e_sqrt.c index f095d6cb3c..fff6d148fe 100644 --- a/sysdeps/ieee754/dbl-64/e_sqrt.c +++ b/sysdeps/ieee754/dbl-64/e_sqrt.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001-2014 Free Software Foundation, Inc. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -66,7 +66,7 @@ __ieee754_sqrt (double x) /*----------------- 2^-1022 <= | x |< 2^1024 -----------------*/ if (k > 0x000fffff && k < 0x7ff00000) { - int rm = fegetround (); + int rm = __fegetround (); fenv_t env; libc_feholdexcept_setround (&env, FE_TONEAREST); double ret; diff --git a/sysdeps/ieee754/dbl-64/gamma_product.c b/sysdeps/ieee754/dbl-64/gamma_product.c index 511cb33561..6c0d6210d9 100644 --- a/sysdeps/ieee754/dbl-64/gamma_product.c +++ b/sysdeps/ieee754/dbl-64/gamma_product.c @@ -1,5 +1,5 @@ /* Compute a product of X, X+1, ..., with an error estimate. - Copyright (C) 2013-2014 Free Software Foundation, Inc. + Copyright (C) 2013-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 diff --git a/sysdeps/ieee754/dbl-64/gamma_productf.c b/sysdeps/ieee754/dbl-64/gamma_productf.c index 60c6f7bceb..df59e1beec 100644 --- a/sysdeps/ieee754/dbl-64/gamma_productf.c +++ b/sysdeps/ieee754/dbl-64/gamma_productf.c @@ -1,5 +1,5 @@ /* Compute a product of X, X+1, ..., with an error estimate. - Copyright (C) 2013-2014 Free Software Foundation, Inc. + Copyright (C) 2013-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 diff --git a/sysdeps/ieee754/dbl-64/halfulp.c b/sysdeps/ieee754/dbl-64/halfulp.c index 68d613412d..028cb748f6 100644 --- a/sysdeps/ieee754/dbl-64/halfulp.c +++ b/sysdeps/ieee754/dbl-64/halfulp.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001-2014 Free Software Foundation, Inc. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/mpa-arch.h b/sysdeps/ieee754/dbl-64/mpa-arch.h index 4a9f45ecc8..0940f0bb57 100644 --- a/sysdeps/ieee754/dbl-64/mpa-arch.h +++ b/sysdeps/ieee754/dbl-64/mpa-arch.h @@ -1,5 +1,5 @@ /* Overridable constants and operations. - Copyright (C) 2013-2014 Free Software Foundation, Inc. + Copyright (C) 2013-2015 Free Software Foundation, Inc. This program is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/mpa.c b/sysdeps/ieee754/dbl-64/mpa.c index 9f7f44fa48..7b52da91d5 100644 --- a/sysdeps/ieee754/dbl-64/mpa.c +++ b/sysdeps/ieee754/dbl-64/mpa.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001-2014 Free Software Foundation, Inc. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -50,8 +50,8 @@ #endif #ifndef NO__CONST -const mp_no mpone = { 1, { 1.0, 1.0 } }; -const mp_no mptwo = { 1, { 1.0, 2.0 } }; +const mp_no __mpone = { 1, { 1.0, 1.0 } }; +const mp_no __mptwo = { 1, { 1.0, 2.0 } }; #endif #ifndef NO___ACR @@ -119,7 +119,8 @@ __cpy (const mp_no *x, mp_no *y, int p) #ifndef NO___MP_DBL /* Convert a multiple precision number *X into a double precision - number *Y, normalized case (|x| >= 2**(-1022))). */ + number *Y, normalized case (|x| >= 2**(-1022))). X has precision + P, which is positive. */ static void norm (const mp_no *x, double *y, int p) { @@ -135,7 +136,7 @@ norm (const mp_no *x, double *y, int p) c = X[1] + R * X[2]; else if (p == 3) c = X[1] + R * (X[2] + R * X[3]); - else if (p == 4) + else /* p == 4. */ c = (X[1] + R * X[2]) + R * R * (X[3] + R * X[4]); } else @@ -877,7 +878,7 @@ __inv (const mp_no *x, mp_no *y, int p) { __cpy (y, &w, p); __mul (x, &w, y, p); - __sub (&mptwo, y, &z, p); + __sub (&__mptwo, y, &z, p); __mul (&w, &z, y, p); } } diff --git a/sysdeps/ieee754/dbl-64/mpa.h b/sysdeps/ieee754/dbl-64/mpa.h index bf1ad873d1..47dd6c457c 100644 --- a/sysdeps/ieee754/dbl-64/mpa.h +++ b/sysdeps/ieee754/dbl-64/mpa.h @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2014 Free Software Foundation, Inc. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -71,8 +71,8 @@ typedef union double d; } number; -extern const mp_no mpone; -extern const mp_no mptwo; +extern const mp_no __mpone; +extern const mp_no __mptwo; #define X x->d #define Y y->d @@ -81,8 +81,6 @@ extern const mp_no mptwo; #define EY y->e #define EZ z->e -#define ABS(x) ((x) < 0 ? -(x) : (x)) - #ifndef RADIXI # define RADIXI 0x1.0p-24 /* 2^-24 */ #endif diff --git a/sysdeps/ieee754/dbl-64/mpatan.c b/sysdeps/ieee754/dbl-64/mpatan.c index 49ecdd2d85..d9ac3d62b1 100644 --- a/sysdeps/ieee754/dbl-64/mpatan.c +++ b/sysdeps/ieee754/dbl-64/mpatan.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001-2014 Free Software Foundation, Inc. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -32,6 +32,7 @@ #include "endian.h" #include "mpa.h" +#include <math.h> #ifndef SECTION # define SECTION @@ -65,7 +66,7 @@ __mpatan (mp_no *x, mp_no *y, int p) else { __mp_dbl (x, &dx, p); - dx = ABS (dx); + dx = fabs (dx); for (m = 6; m > 0; m--) { if (dx > __atan_xm[m].d) @@ -83,10 +84,10 @@ __mpatan (mp_no *x, mp_no *y, int p) { for (i = 0; i < m; i++) { - __add (&mpone, &mpsm, &mpt1, p); + __add (&__mpone, &mpsm, &mpt1, p); __mpsqrt (&mpt1, &mpt2, p); __add (&mpt2, &mpt2, &mpt1, p); - __add (&mptwo, &mpsm, &mpt2, p); + __add (&__mptwo, &mpsm, &mpt2, p); __add (&mpt1, &mpt2, &mpt3, p); __dvd (&mpsm, &mpt3, &mpt1, p); __cpy (&mpt1, &mpsm, p); diff --git a/sysdeps/ieee754/dbl-64/mpatan.h b/sysdeps/ieee754/dbl-64/mpatan.h index 4ac4541ff4..2ceb966998 100644 --- a/sysdeps/ieee754/dbl-64/mpatan.h +++ b/sysdeps/ieee754/dbl-64/mpatan.h @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2014 Free Software Foundation, Inc. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/mpatan2.c b/sysdeps/ieee754/dbl-64/mpatan2.c index 5249492666..78f01d9839 100644 --- a/sysdeps/ieee754/dbl-64/mpatan2.c +++ b/sysdeps/ieee754/dbl-64/mpatan2.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001-2014 Free Software Foundation, Inc. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -52,7 +52,7 @@ __mpatan2 (mp_no *y, mp_no *x, mp_no *z, int p) __mul (&mpt1, &mpt1, &mpt2, p); if (mpt1.d[0] != 0) mpt1.d[0] = 1; - __add (&mpt2, &mpone, &mpt3, p); + __add (&mpt2, &__mpone, &mpt3, p); __mpsqrt (&mpt3, &mpt2, p); __add (&mpt1, &mpt2, &mpt3, p); mpt3.d[0] = Y[0]; diff --git a/sysdeps/ieee754/dbl-64/mpexp.c b/sysdeps/ieee754/dbl-64/mpexp.c index 0096afb836..b5a1def93d 100644 --- a/sysdeps/ieee754/dbl-64/mpexp.c +++ b/sysdeps/ieee754/dbl-64/mpexp.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001-2014 Free Software Foundation, Inc. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -140,7 +140,7 @@ __mpexp (mp_no *x, mp_no *y, int p) } __dbl_mp (kf, &mpk, p); __dvd (&mpt2, &mpk, &mpt1, p); - __add (&mpone, &mpt1, &mpt2, p); + __add (&__mpone, &mpt1, &mpt2, p); /* Raise polynomial value to the power of 2**m. Put result in y. */ for (k = 0, j = 0; k < m;) diff --git a/sysdeps/ieee754/dbl-64/mplog.c b/sysdeps/ieee754/dbl-64/mplog.c index 75adac9df6..53e1b57621 100644 --- a/sysdeps/ieee754/dbl-64/mplog.c +++ b/sysdeps/ieee754/dbl-64/mplog.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001-2014 Free Software Foundation, Inc. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -58,7 +58,7 @@ __mplog (mp_no *x, mp_no *y, int p) mpt1.d[0] = -mpt1.d[0]; __mpexp (&mpt1, &mpt2, p); __mul (x, &mpt2, &mpt1, p); - __sub (&mpt1, &mpone, &mpt2, p); + __sub (&mpt1, &__mpone, &mpt2, p); __add (y, &mpt2, &mpt1, p); __cpy (&mpt1, y, p); } diff --git a/sysdeps/ieee754/dbl-64/mpn2dbl.c b/sysdeps/ieee754/dbl-64/mpn2dbl.c index 4e8a764326..c0cf8e238d 100644 --- a/sysdeps/ieee754/dbl-64/mpn2dbl.c +++ b/sysdeps/ieee754/dbl-64/mpn2dbl.c @@ -1,4 +1,4 @@ -/* Copyright (C) 1995-2014 Free Software Foundation, Inc. +/* Copyright (C) 1995-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 diff --git a/sysdeps/ieee754/dbl-64/mpsqrt.c b/sysdeps/ieee754/dbl-64/mpsqrt.c index 86418237bb..f3ef55eb7f 100644 --- a/sysdeps/ieee754/dbl-64/mpsqrt.c +++ b/sysdeps/ieee754/dbl-64/mpsqrt.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001-2014 Free Software Foundation, Inc. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/mpsqrt.h b/sysdeps/ieee754/dbl-64/mpsqrt.h index 9cf725b021..e322398a7a 100644 --- a/sysdeps/ieee754/dbl-64/mpsqrt.h +++ b/sysdeps/ieee754/dbl-64/mpsqrt.h @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2014 Free Software Foundation, Inc. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/mptan.c b/sysdeps/ieee754/dbl-64/mptan.c index aa16782c2c..0d58d9c786 100644 --- a/sysdeps/ieee754/dbl-64/mptan.c +++ b/sysdeps/ieee754/dbl-64/mptan.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001-2014 Free Software Foundation, Inc. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/mydefs.h b/sysdeps/ieee754/dbl-64/mydefs.h index b53dd94764..31585142b7 100644 --- a/sysdeps/ieee754/dbl-64/mydefs.h +++ b/sysdeps/ieee754/dbl-64/mydefs.h @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2014 Free Software Foundation, Inc. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -30,7 +30,6 @@ typedef int int4; typedef union { int4 i[2]; double x; } mynumber; -#define ABS(x) (((x) > 0) ? (x) : -(x)) #define max(x, y) (((y) > (x)) ? (y) : (x)) #define min(x, y) (((y) < (x)) ? (y) : (x)) #endif diff --git a/sysdeps/ieee754/dbl-64/powtwo.tbl b/sysdeps/ieee754/dbl-64/powtwo.tbl index d1d207155f..2406ee5d65 100644 --- a/sysdeps/ieee754/dbl-64/powtwo.tbl +++ b/sysdeps/ieee754/dbl-64/powtwo.tbl @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2014 Free Software Foundation, Inc. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/root.tbl b/sysdeps/ieee754/dbl-64/root.tbl index f52ab96c27..cb0c8b11ec 100644 --- a/sysdeps/ieee754/dbl-64/root.tbl +++ b/sysdeps/ieee754/dbl-64/root.tbl @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2014 Free Software Foundation, Inc. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/s_asinh.c b/sysdeps/ieee754/dbl-64/s_asinh.c index 5500746848..ebe471015b 100644 --- a/sysdeps/ieee754/dbl-64/s_asinh.c +++ b/sysdeps/ieee754/dbl-64/s_asinh.c @@ -21,6 +21,7 @@ * := sign(x)*log1p(|x| + x^2/(1 + sqrt(1+x^2))) */ +#include <float.h> #include <math.h> #include <math_private.h> @@ -36,12 +37,17 @@ __asinh (double x) int32_t hx, ix; GET_HIGH_WORD (hx, x); ix = hx & 0x7fffffff; - if (__builtin_expect (ix < 0x3e300000, 0)) /* |x|<2**-28 */ + if (__glibc_unlikely (ix < 0x3e300000)) /* |x|<2**-28 */ { + if (fabs (x) < DBL_MIN) + { + double force_underflow = x * x; + math_force_eval (force_underflow); + } if (huge + x > one) return x; /* return x inexact except 0 */ } - if (__builtin_expect (ix > 0x41b00000, 0)) /* |x| > 2**28 */ + if (__glibc_unlikely (ix > 0x41b00000)) /* |x| > 2**28 */ { if (ix >= 0x7ff00000) return x + x; /* x is inf or NaN */ diff --git a/sysdeps/ieee754/dbl-64/s_atan.c b/sysdeps/ieee754/dbl-64/s_atan.c index 0aa145cade..5035ae87bc 100644 --- a/sysdeps/ieee754/dbl-64/s_atan.c +++ b/sysdeps/ieee754/dbl-64/s_atan.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001-2014 Free Software Foundation, Inc. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -41,7 +41,10 @@ #include "MathLib.h" #include "uatan.tbl" #include "atnat.h" +#include <fenv.h> +#include <float.h> #include <math.h> +#include <math_private.h> #include <stap-probe.h> void __mpatan (mp_no *, mp_no *, int); /* see definition in mpatan.c */ @@ -79,13 +82,21 @@ atan (double x) return x + x; /* Regular values of x, including denormals +-0 and +-INF */ + SET_RESTORE_ROUND (FE_TONEAREST); u = (x < 0) ? -x : x; if (u < C) { if (u < B) { if (u < A) - return x; + { + if (u < DBL_MIN) + { + double force_underflow = x * x; + math_force_eval (force_underflow); + } + return x; + } else { /* A <= u < B */ v = x * x; diff --git a/sysdeps/ieee754/dbl-64/s_cbrt.c b/sysdeps/ieee754/dbl-64/s_cbrt.c index ada34b6632..a728968a1d 100644 --- a/sysdeps/ieee754/dbl-64/s_cbrt.c +++ b/sysdeps/ieee754/dbl-64/s_cbrt.c @@ -1,5 +1,5 @@ /* Compute cubic root of double value. - Copyright (C) 1997-2014 Free Software Foundation, Inc. + Copyright (C) 1997-2015 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Dirk Alboth <dirka@uni-paderborn.de> and Ulrich Drepper <drepper@cygnus.com>, 1997. diff --git a/sysdeps/ieee754/dbl-64/s_erf.c b/sysdeps/ieee754/dbl-64/s_erf.c index 3f37397224..ea0a73424e 100644 --- a/sysdeps/ieee754/dbl-64/s_erf.c +++ b/sysdeps/ieee754/dbl-64/s_erf.c @@ -128,7 +128,6 @@ static const double * Coefficients for approximation to erf on [0,0.84375] */ efx = 1.28379167095512586316e-01, /* 0x3FC06EBA, 0x8214DB69 */ - efx8 = 1.02703333676410069053e+00, /* 0x3FF06EBA, 0x8214DB69 */ pp[] = { 1.28379167095512558561e-01, /* 0x3FC06EBA, 0x8214DB68 */ -3.25042107247001499370e-01, /* 0xBFD4CD7D, 0x691CB913 */ -2.84817495755985104766e-02, /* 0xBF9D2A51, 0xDBD7194F */ @@ -211,7 +210,16 @@ __erf (double x) if (ix < 0x3e300000) /* |x|<2**-28 */ { if (ix < 0x00800000) - return 0.125 * (8.0 * x + efx8 * x); /*avoid underflow */ + { + /* Avoid spurious underflow. */ + double ret = 0.0625 * (16.0 * x + (16.0 * efx) * x); + if (fabs (ret) < DBL_MIN) + { + double force_underflow = ret * ret; + math_force_eval (force_underflow); + } + return ret; + } return x + efx * x; } z = x * x; diff --git a/sysdeps/ieee754/dbl-64/s_expm1.c b/sysdeps/ieee754/dbl-64/s_expm1.c index 9c9bb34559..41ef63a786 100644 --- a/sysdeps/ieee754/dbl-64/s_expm1.c +++ b/sysdeps/ieee754/dbl-64/s_expm1.c @@ -109,6 +109,7 @@ */ #include <errno.h> +#include <float.h> #include <math.h> #include <math_private.h> #define one Q[0] @@ -194,6 +195,11 @@ __expm1 (double x) } else if (hx < 0x3c900000) /* when |x|<2**-54, return x */ { + if (fabs (x) < DBL_MIN) + { + double force_underflow = x * x; + math_force_eval (force_underflow); + } t = huge + x; /* return x with inexact flags when x!=0 */ return x - (t - (huge + x)); } diff --git a/sysdeps/ieee754/dbl-64/s_fabs.c b/sysdeps/ieee754/dbl-64/s_fabs.c index c82c4210ed..73c09a269e 100644 --- a/sysdeps/ieee754/dbl-64/s_fabs.c +++ b/sysdeps/ieee754/dbl-64/s_fabs.c @@ -19,15 +19,11 @@ static char rcsid[] = "$NetBSD: s_fabs.c,v 1.7 1995/05/10 20:47:13 jtc Exp $"; */ #include <math.h> -#include <math_private.h> double __fabs (double x) { - u_int32_t high; - GET_HIGH_WORD (high, x); - SET_HIGH_WORD (x, high & 0x7fffffff); - return x; + return __builtin_fabs (x); } weak_alias (__fabs, fabs) #ifdef NO_LONG_DOUBLE diff --git a/sysdeps/ieee754/dbl-64/s_fma.c b/sysdeps/ieee754/dbl-64/s_fma.c index cfaa22d184..716b41273d 100644 --- a/sysdeps/ieee754/dbl-64/s_fma.c +++ b/sysdeps/ieee754/dbl-64/s_fma.c @@ -1,5 +1,5 @@ /* Compute x * y + z as ternary operation. - Copyright (C) 2010-2014 Free Software Foundation, Inc. + Copyright (C) 2010-2015 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Jakub Jelinek <jakub@redhat.com>, 2010. @@ -174,7 +174,7 @@ __fma (double x, double y, double z) } /* Ensure correct sign of exact 0 + 0. */ - if (__builtin_expect ((x == 0 || y == 0) && z == 0, 0)) + if (__glibc_unlikely ((x == 0 || y == 0) && z == 0)) return x * y + z; fenv_t env; @@ -198,16 +198,17 @@ __fma (double x, double y, double z) t1 = m1 - t1; t2 = z - t2; double a2 = t1 + t2; + /* Ensure the arithmetic is not scheduled after feclearexcept call. */ + math_force_eval (m2); + math_force_eval (a2); feclearexcept (FE_INEXACT); - /* If the result is an exact zero, ensure it has the correct - sign. */ + /* If the result is an exact zero, ensure it has the correct sign. */ if (a1 == 0 && m2 == 0) { libc_feupdateenv (&env); - /* Ensure that round-to-nearest value of z + m1 is not - reused. */ - asm volatile ("" : "=m" (z) : "m" (z)); + /* Ensure that round-to-nearest value of z + m1 is not reused. */ + z = math_opt_barrier (z); return z + m1; } @@ -216,7 +217,7 @@ __fma (double x, double y, double z) /* Perform m2 + a2 addition with round to odd. */ u.d = a2 + m2; - if (__builtin_expect (adjust < 0, 0)) + if (__glibc_unlikely (adjust < 0)) { if ((u.ieee.mantissa1 & 1) == 0) u.ieee.mantissa1 |= libc_fetestexcept (FE_INEXACT) != 0; @@ -228,14 +229,14 @@ __fma (double x, double y, double z) /* Reset rounding mode and test for inexact simultaneously. */ int j = libc_feupdateenv_test (&env, FE_INEXACT) != 0; - if (__builtin_expect (adjust == 0, 1)) + if (__glibc_likely (adjust == 0)) { if ((u.ieee.mantissa1 & 1) == 0 && u.ieee.exponent != 0x7ff) u.ieee.mantissa1 |= j; /* Result is a1 + u.d. */ return a1 + u.d; } - else if (__builtin_expect (adjust > 0, 1)) + else if (__glibc_likely (adjust > 0)) { if ((u.ieee.mantissa1 & 1) == 0 && u.ieee.exponent != 0x7ff) u.ieee.mantissa1 |= j; diff --git a/sysdeps/ieee754/dbl-64/s_fmaf.c b/sysdeps/ieee754/dbl-64/s_fmaf.c index 9a062f0bb7..86b8662dbf 100644 --- a/sysdeps/ieee754/dbl-64/s_fmaf.c +++ b/sysdeps/ieee754/dbl-64/s_fmaf.c @@ -1,5 +1,5 @@ /* Compute x * y + z as ternary operation. - Copyright (C) 2010-2014 Free Software Foundation, Inc. + Copyright (C) 2010-2015 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Jakub Jelinek <jakub@redhat.com>, 2010. diff --git a/sysdeps/ieee754/dbl-64/s_fpclassify.c b/sysdeps/ieee754/dbl-64/s_fpclassify.c index d875ca4db3..2e8751277d 100644 --- a/sysdeps/ieee754/dbl-64/s_fpclassify.c +++ b/sysdeps/ieee754/dbl-64/s_fpclassify.c @@ -1,5 +1,5 @@ /* Return classification value corresponding to argument. - Copyright (C) 1997-2014 Free Software Foundation, Inc. + Copyright (C) 1997-2015 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. diff --git a/sysdeps/ieee754/dbl-64/s_issignaling.c b/sysdeps/ieee754/dbl-64/s_issignaling.c index 4631ff9c59..925d85fb21 100644 --- a/sysdeps/ieee754/dbl-64/s_issignaling.c +++ b/sysdeps/ieee754/dbl-64/s_issignaling.c @@ -1,5 +1,5 @@ /* Test for signaling NaN. - Copyright (C) 2013-2014 Free Software Foundation, Inc. + Copyright (C) 2013-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 diff --git a/sysdeps/ieee754/dbl-64/s_llrint.c b/sysdeps/ieee754/dbl-64/s_llrint.c index 3bcca28988..62acb67c38 100644 --- a/sysdeps/ieee754/dbl-64/s_llrint.c +++ b/sysdeps/ieee754/dbl-64/s_llrint.c @@ -1,6 +1,6 @@ /* Round argument to nearest integral value according to current rounding direction. - Copyright (C) 1997-2014 Free Software Foundation, Inc. + Copyright (C) 1997-2015 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. diff --git a/sysdeps/ieee754/dbl-64/s_llround.c b/sysdeps/ieee754/dbl-64/s_llround.c index e1cc1770c8..c8b3c190b5 100644 --- a/sysdeps/ieee754/dbl-64/s_llround.c +++ b/sysdeps/ieee754/dbl-64/s_llround.c @@ -1,5 +1,5 @@ /* Round double value to long long int. - Copyright (C) 1997-2014 Free Software Foundation, Inc. + Copyright (C) 1997-2015 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. diff --git a/sysdeps/ieee754/dbl-64/s_log1p.c b/sysdeps/ieee754/dbl-64/s_log1p.c index ea1dc6ca76..cff555b0aa 100644 --- a/sysdeps/ieee754/dbl-64/s_log1p.c +++ b/sysdeps/ieee754/dbl-64/s_log1p.c @@ -78,6 +78,7 @@ * See HP-15C Advanced Functions Handbook, p.193. */ +#include <float.h> #include <math.h> #include <math_private.h> @@ -107,18 +108,25 @@ __log1p (double x) k = 1; if (hx < 0x3FDA827A) /* x < 0.41422 */ { - if (__builtin_expect (ax >= 0x3ff00000, 0)) /* x <= -1.0 */ + if (__glibc_unlikely (ax >= 0x3ff00000)) /* x <= -1.0 */ { if (x == -1.0) - return -two54 / (x - x); /* log1p(-1)=+inf */ + return -two54 / zero; /* log1p(-1)=-inf */ else return (x - x) / (x - x); /* log1p(x<-1)=NaN */ } - if (__builtin_expect (ax < 0x3e200000, 0)) /* |x| < 2**-29 */ + if (__glibc_unlikely (ax < 0x3e200000)) /* |x| < 2**-29 */ { math_force_eval (two54 + x); /* raise inexact */ if (ax < 0x3c900000) /* |x| < 2**-54 */ - return x; + { + if (fabs (x) < DBL_MIN) + { + double force_underflow = x * x; + math_force_eval (force_underflow); + } + return x; + } else return x - x * x * 0.5; } @@ -127,7 +135,7 @@ __log1p (double x) k = 0; f = x; hu = 1; } /* -0.2929<x<0.41422 */ } - else if (__builtin_expect (hx >= 0x7ff00000, 0)) + else if (__glibc_unlikely (hx >= 0x7ff00000)) return x + x; if (k != 0) { @@ -189,8 +197,3 @@ __log1p (double x) else return k * ln2_hi - ((hfsq - (s * (hfsq + R) + (k * ln2_lo + c))) - f); } -weak_alias (__log1p, log1p) -#ifdef NO_LONG_DOUBLE -strong_alias (__log1p, __log1pl) -weak_alias (__log1p, log1pl) -#endif diff --git a/sysdeps/ieee754/dbl-64/s_logb.c b/sysdeps/ieee754/dbl-64/s_logb.c index c065826dd2..7a6c49abf5 100644 --- a/sysdeps/ieee754/dbl-64/s_logb.c +++ b/sysdeps/ieee754/dbl-64/s_logb.c @@ -30,7 +30,7 @@ __logb (double x) return -1.0 / fabs (x); if (ix >= 0x7ff00000) return x * x; - if (__builtin_expect ((rix = ix >> 20) == 0, 0)) + if (__glibc_unlikely ((rix = ix >> 20) == 0)) { /* POSIX specifies that denormal number is treated as though it were normalized. */ diff --git a/sysdeps/ieee754/dbl-64/s_lrint.c b/sysdeps/ieee754/dbl-64/s_lrint.c index 3ab988473a..9cce37f3cf 100644 --- a/sysdeps/ieee754/dbl-64/s_lrint.c +++ b/sysdeps/ieee754/dbl-64/s_lrint.c @@ -1,6 +1,6 @@ /* Round argument to nearest integral value according to current rounding direction. - Copyright (C) 1997-2014 Free Software Foundation, Inc. + Copyright (C) 1997-2015 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. diff --git a/sysdeps/ieee754/dbl-64/s_lround.c b/sysdeps/ieee754/dbl-64/s_lround.c index ae7f178b98..bdc838a676 100644 --- a/sysdeps/ieee754/dbl-64/s_lround.c +++ b/sysdeps/ieee754/dbl-64/s_lround.c @@ -1,5 +1,5 @@ /* Round double value to long int. - Copyright (C) 1997-2014 Free Software Foundation, Inc. + Copyright (C) 1997-2015 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. diff --git a/sysdeps/ieee754/dbl-64/s_modf.c b/sysdeps/ieee754/dbl-64/s_modf.c index 1dce6381ae..0a1e13008f 100644 --- a/sysdeps/ieee754/dbl-64/s_modf.c +++ b/sysdeps/ieee754/dbl-64/s_modf.c @@ -54,7 +54,7 @@ __modf (double x, double *iptr) } } } - else if (__builtin_expect (j0 > 51, 0)) /* no fraction part */ + else if (__glibc_unlikely (j0 > 51)) /* no fraction part */ { *iptr = x * one; /* We must handle NaNs separately. */ diff --git a/sysdeps/ieee754/dbl-64/s_remquo.c b/sysdeps/ieee754/dbl-64/s_remquo.c index 04b1d81bb3..b081d26b0d 100644 --- a/sysdeps/ieee754/dbl-64/s_remquo.c +++ b/sysdeps/ieee754/dbl-64/s_remquo.c @@ -1,5 +1,5 @@ /* Compute remainder and a congruent to the quotient. - Copyright (C) 1997-2014 Free Software Foundation, Inc. + Copyright (C) 1997-2015 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -60,12 +60,12 @@ __remquo (double x, double y, int *quo) y = fabs (y); cquo = 0; - if (x >= 4 * y) + if (hy <= 0x7fcfffff && x >= 4 * y) { x -= 4 * y; cquo += 4; } - if (x >= 2 * y) + if (hy <= 0x7fdfffff && x >= 2 * y) { x -= 2 * y; cquo += 2; @@ -101,6 +101,9 @@ __remquo (double x, double y, int *quo) *quo = qs ? -cquo : cquo; + /* Ensure correct sign of zero result in round-downward mode. */ + if (x == 0.0) + x = 0.0; if (sx) x = -x; return x; diff --git a/sysdeps/ieee754/dbl-64/s_round.c b/sysdeps/ieee754/dbl-64/s_round.c index 7eacdf3dc1..770f8e6386 100644 --- a/sysdeps/ieee754/dbl-64/s_round.c +++ b/sysdeps/ieee754/dbl-64/s_round.c @@ -1,5 +1,5 @@ /* Round double to integer away from zero. - Copyright (C) 1997-2014 Free Software Foundation, Inc. + Copyright (C) 1997-2015 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. diff --git a/sysdeps/ieee754/dbl-64/s_scalbln.c b/sysdeps/ieee754/dbl-64/s_scalbln.c index 6402927d23..32cd12e3b0 100644 --- a/sysdeps/ieee754/dbl-64/s_scalbln.c +++ b/sysdeps/ieee754/dbl-64/s_scalbln.c @@ -31,7 +31,7 @@ __scalbln (double x, long int n) int32_t k, hx, lx; EXTRACT_WORDS (hx, lx, x); k = (hx & 0x7ff00000) >> 20; /* extract exponent */ - if (__builtin_expect (k == 0, 0)) /* 0 or subnormal x */ + if (__glibc_unlikely (k == 0)) /* 0 or subnormal x */ { if ((lx | (hx & 0x7fffffff)) == 0) return x; /* +-0 */ @@ -39,16 +39,16 @@ __scalbln (double x, long int n) GET_HIGH_WORD (hx, x); k = ((hx & 0x7ff00000) >> 20) - 54; } - if (__builtin_expect (k == 0x7ff, 0)) + if (__glibc_unlikely (k == 0x7ff)) return x + x; /* NaN or Inf */ - if (__builtin_expect (n < -50000, 0)) + if (__glibc_unlikely (n < -50000)) return tiny * __copysign (tiny, x); /*underflow*/ - if (__builtin_expect (n > 50000 || k + n > 0x7fe, 0)) + if (__glibc_unlikely (n > 50000 || k + n > 0x7fe)) return huge * __copysign (huge, x); /* overflow */ /* Now k and n are bounded we know that k = k+n does not overflow. */ k = k + n; - if (__builtin_expect (k > 0, 1)) /* normal result */ + if (__glibc_likely (k > 0)) /* normal result */ { SET_HIGH_WORD (x, (hx & 0x800fffff) | (k << 20)); return x; } @@ -58,8 +58,6 @@ __scalbln (double x, long int n) SET_HIGH_WORD (x, (hx & 0x800fffff) | (k << 20)); return x * twom54; } -weak_alias (__scalbln, scalbln) #ifdef NO_LONG_DOUBLE strong_alias (__scalbln, __scalblnl) -weak_alias (__scalbln, scalblnl) #endif diff --git a/sysdeps/ieee754/dbl-64/s_scalbn.c b/sysdeps/ieee754/dbl-64/s_scalbn.c index 6e7d5ad217..0f58034df5 100644 --- a/sysdeps/ieee754/dbl-64/s_scalbn.c +++ b/sysdeps/ieee754/dbl-64/s_scalbn.c @@ -31,7 +31,7 @@ __scalbn (double x, int n) int32_t k, hx, lx; EXTRACT_WORDS (hx, lx, x); k = (hx & 0x7ff00000) >> 20; /* extract exponent */ - if (__builtin_expect (k == 0, 0)) /* 0 or subnormal x */ + if (__glibc_unlikely (k == 0)) /* 0 or subnormal x */ { if ((lx | (hx & 0x7fffffff)) == 0) return x; /* +-0 */ @@ -39,16 +39,16 @@ __scalbn (double x, int n) GET_HIGH_WORD (hx, x); k = ((hx & 0x7ff00000) >> 20) - 54; } - if (__builtin_expect (k == 0x7ff, 0)) + if (__glibc_unlikely (k == 0x7ff)) return x + x; /* NaN or Inf */ - if (__builtin_expect (n < -50000, 0)) + if (__glibc_unlikely (n < -50000)) return tiny * __copysign (tiny, x); /*underflow*/ - if (__builtin_expect (n > 50000 || k + n > 0x7fe, 0)) + if (__glibc_unlikely (n > 50000 || k + n > 0x7fe)) return huge * __copysign (huge, x); /* overflow */ /* Now k and n are bounded we know that k = k+n does not overflow. */ k = k + n; - if (__builtin_expect (k > 0, 1)) /* normal result */ + if (__glibc_likely (k > 0)) /* normal result */ { SET_HIGH_WORD (x, (hx & 0x800fffff) | (k << 20)); return x; } diff --git a/sysdeps/ieee754/dbl-64/s_signbit.c b/sysdeps/ieee754/dbl-64/s_signbit.c index 89bc8b6d67..764f11a2d2 100644 --- a/sysdeps/ieee754/dbl-64/s_signbit.c +++ b/sysdeps/ieee754/dbl-64/s_signbit.c @@ -1,5 +1,5 @@ /* Return nonzero value if number is negative. - Copyright (C) 1997-2014 Free Software Foundation, Inc. + Copyright (C) 1997-2015 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. diff --git a/sysdeps/ieee754/dbl-64/s_sin.c b/sysdeps/ieee754/dbl-64/s_sin.c index 6105e9fbdf..eff120e88d 100644 --- a/sysdeps/ieee754/dbl-64/s_sin.c +++ b/sysdeps/ieee754/dbl-64/s_sin.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001-2014 Free Software Foundation, Inc. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -48,10 +48,12 @@ #include <errno.h> +#include <float.h> #include "endian.h" #include "mydefs.h" #include "usncs.h" #include "MathLib.h" +#include <math.h> #include <math_private.h> #include <fenv.h> @@ -294,7 +296,14 @@ __sin (double x) m = u.i[HIGH_HALF]; k = 0x7fffffff & m; /* no sign */ if (k < 0x3e500000) /* if x->0 =>sin(x)=x */ - retval = x; + { + if (fabs (x) < DBL_MIN) + { + double force_underflow = x * x; + math_force_eval (force_underflow); + } + retval = x; + } /*---------------------------- 2^-26 < |x|< 0.25 ----------------------*/ else if (k < 0x3fd00000) { @@ -355,7 +364,7 @@ __sin (double x) da = xn * mp3; a = y - da; da = (y - a) - da; - eps = ABS (x) * 1.2e-30; + eps = fabs (x) * 1.2e-30; switch (n) { /* quarter of unit circle */ @@ -447,19 +456,21 @@ __sin (double x) } else { + double t; if (a > 0) { m = 1; + t = a; db = da; } else { m = 0; - a = -a; + t = -a; db = -da; } - u.x = big + a; - y = a - (u.x - big); + u.x = big + t; + y = t - (u.x - big); res = do_sin (u, y, db, &cor); cor = (cor > 0) ? 1.035 * cor + eps : 1.035 * cor - eps; retval = ((res == res + cor) ? ((m) ? res : -res) @@ -528,7 +539,7 @@ __cos (double x) else if (k < 0x3feb6000) { /* 2^-27 < |x| < 0.855469 */ - y = ABS (x); + y = fabs (x); u.x = big + y; y = y - (u.x - big); res = do_cos (u, y, &cor); @@ -537,7 +548,7 @@ __cos (double x) else if (k < 0x400368fd) { /* 0.855469 <|x|<2.426265 */ ; - y = hp0 - ABS (x); + y = hp0 - fabs (x); a = y + hp1; da = (y - a) + hp1; xx = a * a; @@ -580,7 +591,7 @@ __cos (double x) da = xn * mp3; a = y - da; da = (y - a) - da; - eps = ABS (x) * 1.2e-30; + eps = fabs (x) * 1.2e-30; switch (n) { @@ -671,19 +682,21 @@ __cos (double x) } else { + double t; if (a > 0) { m = 1; + t = a; db = da; } else { m = 0; - a = -a; + t = -a; db = -da; } - u.x = big + a; - y = a - (u.x - big); + u.x = big + t; + y = t - (u.x - big); res = do_sin (u, y, db, &cor); cor = (cor > 0) ? 1.035 * cor + eps : 1.035 * cor - eps; retval = ((res == res + cor) ? ((m) ? res : -res) @@ -737,7 +750,7 @@ slow (double x) return res; else { - __dubsin (ABS (x), 0, w); + __dubsin (fabs (x), 0, w); if (w[0] == w[0] + 1.000000001 * w[1]) return (x > 0) ? w[0] : -w[0]; else @@ -756,7 +769,7 @@ slow1 (double x) { mynumber u; double w[2], y, cor, res; - y = ABS (x); + y = fabs (x); u.x = big + y; y = y - (u.x - big); res = do_sin_slow (u, y, 0, 0, &cor); @@ -764,7 +777,7 @@ slow1 (double x) return (x > 0) ? res : -res; else { - __dubsin (ABS (x), 0, w); + __dubsin (fabs (x), 0, w); if (w[0] == w[0] + 1.000000005 * w[1]) return (x > 0) ? w[0] : -w[0]; else @@ -783,7 +796,7 @@ slow2 (double x) mynumber u; double w[2], y, y1, y2, cor, res, del; - y = ABS (x); + y = fabs (x); y = hp0 - y; if (y >= 0) { @@ -802,7 +815,7 @@ slow2 (double x) return (x > 0) ? res : -res; else { - y = ABS (x) - hp0; + y = fabs (x) - hp0; y1 = y - hp1; y2 = (y - y1) - hp1; __docos (y1, y2, w); @@ -830,9 +843,9 @@ sloww (double x, double dx, double orig) int4 n; res = TAYLOR_SLOW (x, dx, cor); if (cor > 0) - cor = 1.0005 * cor + ABS (orig) * 3.1e-30; + cor = 1.0005 * cor + fabs (orig) * 3.1e-30; else - cor = 1.0005 * cor - ABS (orig) * 3.1e-30; + cor = 1.0005 * cor - fabs (orig) * 3.1e-30; if (res == res + cor) return res; @@ -840,9 +853,9 @@ sloww (double x, double dx, double orig) { (x > 0) ? __dubsin (x, dx, w) : __dubsin (-x, -dx, w); if (w[1] > 0) - cor = 1.000000001 * w[1] + ABS (orig) * 1.1e-30; + cor = 1.000000001 * w[1] + fabs (orig) * 1.1e-30; else - cor = 1.000000001 * w[1] - ABS (orig) * 1.1e-30; + cor = 1.000000001 * w[1] - fabs (orig) * 1.1e-30; if (w[0] == w[0] + cor) return (x > 0) ? w[0] : -w[0]; @@ -866,9 +879,9 @@ sloww (double x, double dx, double orig) } (a > 0) ? __dubsin (a, da, w) : __dubsin (-a, -da, w); if (w[1] > 0) - cor = 1.000000001 * w[1] + ABS (orig) * 1.1e-40; + cor = 1.000000001 * w[1] + fabs (orig) * 1.1e-40; else - cor = 1.000000001 * w[1] - ABS (orig) * 1.1e-40; + cor = 1.000000001 * w[1] - fabs (orig) * 1.1e-40; if (w[0] == w[0] + cor) return (a > 0) ? w[0] : -w[0]; @@ -894,7 +907,7 @@ sloww1 (double x, double dx, double orig, int m) u.x = big + x; y = x - (u.x - big); - res = do_sin_slow (u, y, dx, 3.1e-30 * ABS (orig), &cor); + res = do_sin_slow (u, y, dx, 3.1e-30 * fabs (orig), &cor); if (res == res + cor) return (m > 0) ? res : -res; @@ -903,9 +916,9 @@ sloww1 (double x, double dx, double orig, int m) __dubsin (x, dx, w); if (w[1] > 0) - cor = 1.000000005 * w[1] + 1.1e-30 * ABS (orig); + cor = 1.000000005 * w[1] + 1.1e-30 * fabs (orig); else - cor = 1.000000005 * w[1] - 1.1e-30 * ABS (orig); + cor = 1.000000005 * w[1] - 1.1e-30 * fabs (orig); if (w[0] == w[0] + cor) return (m > 0) ? w[0] : -w[0]; @@ -930,7 +943,7 @@ sloww2 (double x, double dx, double orig, int n) u.x = big + x; y = x - (u.x - big); - res = do_cos_slow (u, y, dx, 3.1e-30 * ABS (orig), &cor); + res = do_cos_slow (u, y, dx, 3.1e-30 * fabs (orig), &cor); if (res == res + cor) return (n & 2) ? -res : res; @@ -939,9 +952,9 @@ sloww2 (double x, double dx, double orig, int n) __docos (x, dx, w); if (w[1] > 0) - cor = 1.000000005 * w[1] + 1.1e-30 * ABS (orig); + cor = 1.000000005 * w[1] + 1.1e-30 * fabs (orig); else - cor = 1.000000005 * w[1] - 1.1e-30 * ABS (orig); + cor = 1.000000005 * w[1] - 1.1e-30 * fabs (orig); if (w[0] == w[0] + cor) return (n & 2) ? -w[0] : w[0]; @@ -996,7 +1009,7 @@ bsloww1 (double x, double dx, double orig, int n) mynumber u; double w[2], y, cor, res; - y = ABS (x); + y = fabs (x); u.x = big + y; y = y - (u.x - big); dx = (x > 0) ? dx : -dx; @@ -1005,7 +1018,7 @@ bsloww1 (double x, double dx, double orig, int n) return (x > 0) ? res : -res; else { - __dubsin (ABS (x), dx, w); + __dubsin (fabs (x), dx, w); if (w[1] > 0) cor = 1.000000005 * w[1] + 1.1e-24; @@ -1033,7 +1046,7 @@ bsloww2 (double x, double dx, double orig, int n) mynumber u; double w[2], y, cor, res; - y = ABS (x); + y = fabs (x); u.x = big + y; y = y - (u.x - big); dx = (x > 0) ? dx : -dx; @@ -1042,7 +1055,7 @@ bsloww2 (double x, double dx, double orig, int n) return (n & 2) ? -res : res; else { - __docos (ABS (x), dx, w); + __docos (fabs (x), dx, w); if (w[1] > 0) cor = 1.000000005 * w[1] + 1.1e-24; @@ -1068,7 +1081,7 @@ cslow2 (double x) mynumber u; double w[2], y, cor, res; - y = ABS (x); + y = fabs (x); u.x = big + y; y = y - (u.x - big); res = do_cos_slow (u, y, 0, 0, &cor); @@ -1076,7 +1089,7 @@ cslow2 (double x) return res; else { - y = ABS (x); + y = fabs (x); __docos (y, 0, w); if (w[0] == w[0] + 1.000000005 * w[1]) return w[0]; @@ -1105,9 +1118,9 @@ csloww (double x, double dx, double orig) res = TAYLOR_SLOW (x, dx, cor); if (cor > 0) - cor = 1.0005 * cor + ABS (orig) * 3.1e-30; + cor = 1.0005 * cor + fabs (orig) * 3.1e-30; else - cor = 1.0005 * cor - ABS (orig) * 3.1e-30; + cor = 1.0005 * cor - fabs (orig) * 3.1e-30; if (res == res + cor) return res; @@ -1116,9 +1129,9 @@ csloww (double x, double dx, double orig) (x > 0) ? __dubsin (x, dx, w) : __dubsin (-x, -dx, w); if (w[1] > 0) - cor = 1.000000001 * w[1] + ABS (orig) * 1.1e-30; + cor = 1.000000001 * w[1] + fabs (orig) * 1.1e-30; else - cor = 1.000000001 * w[1] - ABS (orig) * 1.1e-30; + cor = 1.000000001 * w[1] - fabs (orig) * 1.1e-30; if (w[0] == w[0] + cor) return (x > 0) ? w[0] : -w[0]; @@ -1143,9 +1156,9 @@ csloww (double x, double dx, double orig) (a > 0) ? __dubsin (a, da, w) : __dubsin (-a, -da, w); if (w[1] > 0) - cor = 1.000000001 * w[1] + ABS (orig) * 1.1e-40; + cor = 1.000000001 * w[1] + fabs (orig) * 1.1e-40; else - cor = 1.000000001 * w[1] - ABS (orig) * 1.1e-40; + cor = 1.000000001 * w[1] - fabs (orig) * 1.1e-40; if (w[0] == w[0] + cor) return (a > 0) ? w[0] : -w[0]; @@ -1171,7 +1184,7 @@ csloww1 (double x, double dx, double orig, int m) u.x = big + x; y = x - (u.x - big); - res = do_sin_slow (u, y, dx, 3.1e-30 * ABS (orig), &cor); + res = do_sin_slow (u, y, dx, 3.1e-30 * fabs (orig), &cor); if (res == res + cor) return (m > 0) ? res : -res; @@ -1179,9 +1192,9 @@ csloww1 (double x, double dx, double orig, int m) { __dubsin (x, dx, w); if (w[1] > 0) - cor = 1.000000005 * w[1] + 1.1e-30 * ABS (orig); + cor = 1.000000005 * w[1] + 1.1e-30 * fabs (orig); else - cor = 1.000000005 * w[1] - 1.1e-30 * ABS (orig); + cor = 1.000000005 * w[1] - 1.1e-30 * fabs (orig); if (w[0] == w[0] + cor) return (m > 0) ? w[0] : -w[0]; else @@ -1206,7 +1219,7 @@ csloww2 (double x, double dx, double orig, int n) u.x = big + x; y = x - (u.x - big); - res = do_cos_slow (u, y, dx, 3.1e-30 * ABS (orig), &cor); + res = do_cos_slow (u, y, dx, 3.1e-30 * fabs (orig), &cor); if (res == res + cor) return (n) ? -res : res; @@ -1214,9 +1227,9 @@ csloww2 (double x, double dx, double orig, int n) { __docos (x, dx, w); if (w[1] > 0) - cor = 1.000000005 * w[1] + 1.1e-30 * ABS (orig); + cor = 1.000000005 * w[1] + 1.1e-30 * fabs (orig); else - cor = 1.000000005 * w[1] - 1.1e-30 * ABS (orig); + cor = 1.000000005 * w[1] - 1.1e-30 * fabs (orig); if (w[0] == w[0] + cor) return (n) ? -w[0] : w[0]; else diff --git a/sysdeps/ieee754/dbl-64/s_sincos.c b/sysdeps/ieee754/dbl-64/s_sincos.c index e4371069d4..c8a99991cc 100644 --- a/sysdeps/ieee754/dbl-64/s_sincos.c +++ b/sysdeps/ieee754/dbl-64/s_sincos.c @@ -1,5 +1,5 @@ /* Compute sine and cosine of argument. - Copyright (C) 1997-2014 Free Software Foundation, Inc. + Copyright (C) 1997-2015 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -17,6 +17,7 @@ License along with the GNU C Library; if not, see <http://www.gnu.org/licenses/>. */ +#include <errno.h> #include <math.h> #include <math_private.h> @@ -36,6 +37,8 @@ __sincos (double x, double *sinx, double *cosx) { /* sin(Inf or NaN) is NaN */ *sinx = *cosx = x - x; + if (__isinf_ns (x)) + __set_errno (EDOM); } else { diff --git a/sysdeps/ieee754/dbl-64/s_tan.c b/sysdeps/ieee754/dbl-64/s_tan.c index a72530dadc..dcb4aca2de 100644 --- a/sysdeps/ieee754/dbl-64/s_tan.c +++ b/sysdeps/ieee754/dbl-64/s_tan.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001-2014 Free Software Foundation, Inc. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/s_trunc.c b/sysdeps/ieee754/dbl-64/s_trunc.c index bd777f6860..ac6c7597b5 100644 --- a/sysdeps/ieee754/dbl-64/s_trunc.c +++ b/sysdeps/ieee754/dbl-64/s_trunc.c @@ -1,5 +1,5 @@ /* Truncate argument to nearest integral value not larger than the argument. - Copyright (C) 1997-2014 Free Software Foundation, Inc. + Copyright (C) 1997-2015 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. diff --git a/sysdeps/ieee754/dbl-64/sincos32.c b/sysdeps/ieee754/dbl-64/sincos32.c index 6b2fa878a4..d942a1f29f 100644 --- a/sysdeps/ieee754/dbl-64/sincos32.c +++ b/sysdeps/ieee754/dbl-64/sincos32.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001-2014 Free Software Foundation, Inc. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -42,6 +42,7 @@ #include "endian.h" #include "mpa.h" #include "sincos32.h" +#include <math.h> #include <math_private.h> #include <stap-probe.h> @@ -118,11 +119,11 @@ __c32 (mp_no *x, mp_no *y, mp_no *z, int p) __mul (&c, &s, &t, p); __sub (&s, &t, &t1, p); __add (&t1, &t1, &s, p); - __sub (&mptwo, &c, &t1, p); + __sub (&__mptwo, &c, &t1, p); __mul (&t1, &c, &t2, p); __add (&t2, &t2, &c, p); } - __sub (&mpone, &c, y, p); + __sub (&__mpone, &c, y, p); __cpy (&s, z, p); } @@ -318,7 +319,7 @@ __mpranred (double x, mp_no *y, int p) int i, k, n; mp_no a, b, c; - if (ABS (x) < 2.8e14) + if (fabs (x) < 2.8e14) { t = (x * hpinv.d + toint.d); xn = t - toint.d; @@ -352,7 +353,7 @@ __mpranred (double x, mp_no *y, int p) if (c.d[1] >= HALFRAD) { t += 1.0; - __sub (&c, &mpone, &b, p); + __sub (&c, &__mpone, &b, p); __mul (&b, &hp, y, p); } else diff --git a/sysdeps/ieee754/dbl-64/sincos32.h b/sysdeps/ieee754/dbl-64/sincos32.h index e289b165c2..ee70b08743 100644 --- a/sysdeps/ieee754/dbl-64/sincos32.h +++ b/sysdeps/ieee754/dbl-64/sincos32.h @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2014 Free Software Foundation, Inc. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/sincostab.c b/sysdeps/ieee754/dbl-64/sincostab.c index e3e24e0f84..44cb525f42 100644 --- a/sysdeps/ieee754/dbl-64/sincostab.c +++ b/sysdeps/ieee754/dbl-64/sincostab.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2014 Free Software Foundation, Inc. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/slowexp.c b/sysdeps/ieee754/dbl-64/slowexp.c index e32742d5d4..7950111412 100644 --- a/sysdeps/ieee754/dbl-64/slowexp.c +++ b/sysdeps/ieee754/dbl-64/slowexp.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001-2014 Free Software Foundation, Inc. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/slowpow.c b/sysdeps/ieee754/dbl-64/slowpow.c index b2e403308e..117225ea46 100644 --- a/sysdeps/ieee754/dbl-64/slowpow.c +++ b/sysdeps/ieee754/dbl-64/slowpow.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001-2014 Free Software Foundation, Inc. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/t_exp.c b/sysdeps/ieee754/dbl-64/t_exp.c index f32fa57a47..8d5c7b81a3 100644 --- a/sysdeps/ieee754/dbl-64/t_exp.c +++ b/sysdeps/ieee754/dbl-64/t_exp.c @@ -1,5 +1,5 @@ /* Accurate tables for exp(). - Copyright (C) 1998-2014 Free Software Foundation, Inc. + Copyright (C) 1998-2015 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Geoffrey Keating <geoffk@ozemail.com.au> diff --git a/sysdeps/ieee754/dbl-64/uasncs.h b/sysdeps/ieee754/dbl-64/uasncs.h index db56a53542..c60e93c0a2 100644 --- a/sysdeps/ieee754/dbl-64/uasncs.h +++ b/sysdeps/ieee754/dbl-64/uasncs.h @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2014 Free Software Foundation, Inc. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/uatan.tbl b/sysdeps/ieee754/dbl-64/uatan.tbl index 24e72e7105..f1337fce21 100644 --- a/sysdeps/ieee754/dbl-64/uatan.tbl +++ b/sysdeps/ieee754/dbl-64/uatan.tbl @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2014 Free Software Foundation, Inc. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/uexp.h b/sysdeps/ieee754/dbl-64/uexp.h index dd59451e3a..6817eaf074 100644 --- a/sysdeps/ieee754/dbl-64/uexp.h +++ b/sysdeps/ieee754/dbl-64/uexp.h @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2014 Free Software Foundation, Inc. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/uexp.tbl b/sysdeps/ieee754/dbl-64/uexp.tbl index 90ff55de41..859e908c5d 100644 --- a/sysdeps/ieee754/dbl-64/uexp.tbl +++ b/sysdeps/ieee754/dbl-64/uexp.tbl @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2014 Free Software Foundation, Inc. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/ulog.h b/sysdeps/ieee754/dbl-64/ulog.h index 70389b4cdb..ce4db49024 100644 --- a/sysdeps/ieee754/dbl-64/ulog.h +++ b/sysdeps/ieee754/dbl-64/ulog.h @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2014 Free Software Foundation, Inc. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/ulog.tbl b/sysdeps/ieee754/dbl-64/ulog.tbl index 454fb61b22..e5d2fdd34b 100644 --- a/sysdeps/ieee754/dbl-64/ulog.tbl +++ b/sysdeps/ieee754/dbl-64/ulog.tbl @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2014 Free Software Foundation, Inc. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/upow.h b/sysdeps/ieee754/dbl-64/upow.h index f26240652e..c8569a9456 100644 --- a/sysdeps/ieee754/dbl-64/upow.h +++ b/sysdeps/ieee754/dbl-64/upow.h @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2014 Free Software Foundation, Inc. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/upow.tbl b/sysdeps/ieee754/dbl-64/upow.tbl index f5e7660d91..1d506b5ad7 100644 --- a/sysdeps/ieee754/dbl-64/upow.tbl +++ b/sysdeps/ieee754/dbl-64/upow.tbl @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2014 Free Software Foundation, Inc. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/urem.h b/sysdeps/ieee754/dbl-64/urem.h index 3bd2883be7..f612c135f8 100644 --- a/sysdeps/ieee754/dbl-64/urem.h +++ b/sysdeps/ieee754/dbl-64/urem.h @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2014 Free Software Foundation, Inc. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/uroot.h b/sysdeps/ieee754/dbl-64/uroot.h index 2db841ffa9..4bbcc3bac7 100644 --- a/sysdeps/ieee754/dbl-64/uroot.h +++ b/sysdeps/ieee754/dbl-64/uroot.h @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2014 Free Software Foundation, Inc. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/usncs.h b/sysdeps/ieee754/dbl-64/usncs.h index b3296ca9c8..e3a8193c40 100644 --- a/sysdeps/ieee754/dbl-64/usncs.h +++ b/sysdeps/ieee754/dbl-64/usncs.h @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2014 Free Software Foundation, Inc. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/utan.h b/sysdeps/ieee754/dbl-64/utan.h index 374132be5f..efa13ca5cf 100644 --- a/sysdeps/ieee754/dbl-64/utan.h +++ b/sysdeps/ieee754/dbl-64/utan.h @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2014 Free Software Foundation, Inc. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/utan.tbl b/sysdeps/ieee754/dbl-64/utan.tbl index 25c59c2759..c088588229 100644 --- a/sysdeps/ieee754/dbl-64/utan.tbl +++ b/sysdeps/ieee754/dbl-64/utan.tbl @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2014 Free Software Foundation, Inc. + * Copyright (C) 2001-2015 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/w_exp.c b/sysdeps/ieee754/dbl-64/w_exp.c index 3c2c3e5f40..61b2dbbd0f 100644 --- a/sysdeps/ieee754/dbl-64/w_exp.c +++ b/sysdeps/ieee754/dbl-64/w_exp.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2011-2014 Free Software Foundation, Inc. +/* Copyright (C) 2011-2015 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@gmail.com>, 2011. @@ -24,9 +24,9 @@ double __exp (double x) { double z = __ieee754_exp (x); - if (__builtin_expect (!__finite (z) || z == 0, 0) - && __finite (x) && _LIB_VERSION != _IEEE_) - return __kernel_standard (x, x, 6 + !!__signbit (x)); + if (__builtin_expect (!isfinite (z) || z == 0, 0) + && isfinite (x) && _LIB_VERSION != _IEEE_) + return __kernel_standard (x, x, 6 + !!signbit (x)); return z; } diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/e_acosh.c b/sysdeps/ieee754/dbl-64/wordsize-64/e_acosh.c index 26268f2498..ccccdaf106 100644 --- a/sysdeps/ieee754/dbl-64/wordsize-64/e_acosh.c +++ b/sysdeps/ieee754/dbl-64/wordsize-64/e_acosh.c @@ -39,7 +39,7 @@ __ieee754_acosh (double x) if (hx > INT64_C (0x4000000000000000)) { - if (__builtin_expect (hx >= INT64_C (0x41b0000000000000), 0)) + if (__glibc_unlikely (hx >= INT64_C (0x41b0000000000000))) { /* x > 2**28 */ if (hx >= INT64_C (0x7ff0000000000000)) @@ -53,13 +53,13 @@ __ieee754_acosh (double x) double t = x * x; return __ieee754_log (2.0 * x - one / (x + __ieee754_sqrt (t - one))); } - else if (__builtin_expect (hx > INT64_C (0x3ff0000000000000), 1)) + else if (__glibc_likely (hx > INT64_C (0x3ff0000000000000))) { /* 1<x<2 */ double t = x - one; return __log1p (t + __ieee754_sqrt (2.0 * t + t * t)); } - else if (__builtin_expect (hx == INT64_C (0x3ff0000000000000), 1)) + else if (__glibc_likely (hx == INT64_C (0x3ff0000000000000))) return 0.0; /* acosh(1) = 0 */ else /* x < 1 */ return (x - x) / (x - x); diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/e_cosh.c b/sysdeps/ieee754/dbl-64/wordsize-64/e_cosh.c index 84593521cc..fca80b13f9 100644 --- a/sysdeps/ieee754/dbl-64/wordsize-64/e_cosh.c +++ b/sysdeps/ieee754/dbl-64/wordsize-64/e_cosh.c @@ -50,9 +50,10 @@ __ieee754_cosh (double x) if (ix < 0x40360000) { /* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */ if(ix<0x3fd62e43) { + if (ix<0x3c800000) /* cosh(tiny) = 1 */ + return one; t = __expm1(fabs(x)); w = one+t; - if (ix<0x3c800000) return w; /* cosh(tiny) = 1 */ return one+(t*t)/(w+w); } diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/e_log10.c b/sysdeps/ieee754/dbl-64/wordsize-64/e_log10.c index dcb7b58a1b..4f5a81669e 100644 --- a/sysdeps/ieee754/dbl-64/wordsize-64/e_log10.c +++ b/sysdeps/ieee754/dbl-64/wordsize-64/e_log10.c @@ -64,16 +64,16 @@ __ieee754_log10 (double x) k = 0; if (hx < INT64_C(0x0010000000000000)) { /* x < 2**-1022 */ - if (__builtin_expect ((hx & UINT64_C(0x7fffffffffffffff)) == 0, 0)) + if (__glibc_unlikely ((hx & UINT64_C(0x7fffffffffffffff)) == 0)) return -two54 / (x - x); /* log(+-0)=-inf */ - if (__builtin_expect (hx < 0, 0)) + if (__glibc_unlikely (hx < 0)) return (x - x) / (x - x); /* log(-#) = NaN */ k -= 54; x *= two54; /* subnormal number, scale up x */ EXTRACT_WORDS64 (hx, x); } /* scale up resulted in a NaN number */ - if (__builtin_expect (hx >= UINT64_C(0x7ff0000000000000), 0)) + if (__glibc_unlikely (hx >= UINT64_C(0x7ff0000000000000))) return x + x; k += (hx >> 52) - 1023; i = ((uint64_t) k & UINT64_C(0x8000000000000000)) >> 63; diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/e_log2.c b/sysdeps/ieee754/dbl-64/wordsize-64/e_log2.c index 6dc7b7d217..5ccb78cf03 100644 --- a/sysdeps/ieee754/dbl-64/wordsize-64/e_log2.c +++ b/sysdeps/ieee754/dbl-64/wordsize-64/e_log2.c @@ -80,15 +80,15 @@ __ieee754_log2 (double x) k = 0; if (hx < INT64_C(0x0010000000000000)) { /* x < 2**-1022 */ - if (__builtin_expect ((hx & UINT64_C(0x7fffffffffffffff)) == 0, 0)) + if (__glibc_unlikely ((hx & UINT64_C(0x7fffffffffffffff)) == 0)) return -two54 / (x - x); /* log(+-0)=-inf */ - if (__builtin_expect (hx < 0, 0)) + if (__glibc_unlikely (hx < 0)) return (x - x) / (x - x); /* log(-#) = NaN */ k -= 54; x *= two54; /* subnormal number, scale up x */ EXTRACT_WORDS64 (hx, x); } - if (__builtin_expect (hx >= UINT64_C(0x7ff0000000000000), 0)) + if (__glibc_unlikely (hx >= UINT64_C(0x7ff0000000000000))) return x + x; k += (hx >> 52) - 1023; hx &= UINT64_C(0x000fffffffffffff); diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_floor.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_floor.c index 3823a53ad4..0dbadb15eb 100644 --- a/sysdeps/ieee754/dbl-64/wordsize-64/s_floor.c +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_floor.c @@ -1,5 +1,5 @@ /* Round double to integer away from zero. - Copyright (C) 2011-2014 Free Software Foundation, Inc. + Copyright (C) 2011-2015 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 2011. diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_frexp.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_frexp.c index d03e33ee1f..f1fa3154ee 100644 --- a/sysdeps/ieee754/dbl-64/wordsize-64/s_frexp.c +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_frexp.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2011-2014 Free Software Foundation, Inc. +/* Copyright (C) 2011-2015 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@gmail.com>, 2011. @@ -39,11 +39,11 @@ __frexp (double x, int *eptr) int32_t ex = 0x7ff & (ix >> 52); int e = 0; - if (__builtin_expect (ex != 0x7ff && x != 0.0, 1)) + if (__glibc_likely (ex != 0x7ff && x != 0.0)) { /* Not zero and finite. */ e = ex - 1022; - if (__builtin_expect (ex == 0, 0)) + if (__glibc_unlikely (ex == 0)) { /* Subnormal. */ x *= 0x1p54; diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_issignaling.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_issignaling.c index 3b8be8430c..67d77d5732 100644 --- a/sysdeps/ieee754/dbl-64/wordsize-64/s_issignaling.c +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_issignaling.c @@ -1,5 +1,5 @@ /* Test for signaling NaN. - Copyright (C) 2013-2014 Free Software Foundation, Inc. + Copyright (C) 2013-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 diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_llround.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_llround.c index da180fef55..8bb52524b7 100644 --- a/sysdeps/ieee754/dbl-64/wordsize-64/s_llround.c +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_llround.c @@ -1,5 +1,5 @@ /* Round double value to long long int. - Copyright (C) 1997-2014 Free Software Foundation, Inc. + Copyright (C) 1997-2015 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -21,6 +21,7 @@ #define __lround __hidden___lround #include <math.h> +#include <sysdep.h> #include <math_private.h> @@ -64,16 +65,31 @@ __llround (double x) weak_alias (__llround, llround) #ifdef NO_LONG_DOUBLE -strong_alias (__llround, __lroundl) -weak_alias (__llround, lroundl) +strong_alias (__llround, __llroundl) +weak_alias (__llround, llroundl) #endif -/* long has the same width as long long on 64-bit machines. */ +/* long has the same width as long long on LP64 machines, so use an alias. + If building for ILP32 on a machine with 64-bit registers, however, + use a cast if necessary. */ #undef lround #undef __lround +#if !defined (_LP64) && REGISTER_CAST_INT32_TO_INT64 +long int +__lround (double x) +{ + return __llround (x); +} +weak_alias (__lround, lround) +# ifdef NO_LONG_DOUBLE +strong_alias (__lround, __lroundl) +weak_alias (__lround, lroundl) +# endif +#else strong_alias (__llround, __lround) weak_alias (__llround, lround) -#ifdef NO_LONG_DOUBLE -strong_alias (__llround, __llroundl) -weak_alias (__llround, llroundl) +# ifdef NO_LONG_DOUBLE +strong_alias (__llround, __lroundl) +weak_alias (__llround, lroundl) +# endif #endif diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_logb.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_logb.c index e51c849592..1d20e93e6d 100644 --- a/sysdeps/ieee754/dbl-64/wordsize-64/s_logb.c +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_logb.c @@ -1,5 +1,5 @@ /* Compute radix independent exponent. - Copyright (C) 2011-2014 Free Software Foundation, Inc. + Copyright (C) 2011-2015 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@gmail.com>, 2011. @@ -34,7 +34,7 @@ __logb (double x) ex = ix >> 52; if (ex == 0x7ff) return x * x; - if (__builtin_expect (ex == 0, 0)) + if (__glibc_unlikely (ex == 0)) { int m = __builtin_clzll (ix); ex -= m - 12; diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_nearbyint.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_nearbyint.c index a58a6202ef..8293819981 100644 --- a/sysdeps/ieee754/dbl-64/wordsize-64/s_nearbyint.c +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_nearbyint.c @@ -46,7 +46,7 @@ __nearbyint(double x) double t = w-TWO52[sx]; math_opt_barrier(t); libc_fesetenv (&env); - return copysign(t, x); + return __copysign (t, x); } } else { if(j0==0x400) return x+x; /* inf or NaN */ diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_remquo.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_remquo.c index b22503f5c7..304e9d0fb4 100644 --- a/sysdeps/ieee754/dbl-64/wordsize-64/s_remquo.c +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_remquo.c @@ -1,5 +1,5 @@ /* Compute remainder and a congruent to the quotient. - Copyright (C) 1997-2014 Free Software Foundation, Inc. + Copyright (C) 1997-2015 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -40,7 +40,7 @@ __remquo (double x, double y, int *quo) hx &= UINT64_C(0x7fffffffffffffff); /* Purge off exception values. */ - if (__builtin_expect (hy == 0, 0)) + if (__glibc_unlikely (hy == 0)) return (x * y) / (x * y); /* y = 0 */ if (__builtin_expect (hx >= UINT64_C(0x7ff0000000000000) /* x not finite */ || hy > UINT64_C(0x7ff0000000000000), 0))/* y is NaN */ @@ -49,22 +49,22 @@ __remquo (double x, double y, int *quo) if (hy <= UINT64_C(0x7fbfffffffffffff)) x = __ieee754_fmod (x, 8 * y); /* now x < 8y */ - if (__builtin_expect (hx == hy, 0)) + if (__glibc_unlikely (hx == hy)) { *quo = qs ? -1 : 1; return zero * x; } - INSERT_WORDS64 (x, hx); + x = fabs (x); INSERT_WORDS64 (y, hy); cquo = 0; - if (x >= 4 * y) + if (hy <= UINT64_C(0x7fcfffffffffffff) && x >= 4 * y) { x -= 4 * y; cquo += 4; } - if (x >= 2 * y) + if (hy <= UINT64_C(0x7fdfffffffffffff) && x >= 2 * y) { x -= 2 * y; cquo += 2; @@ -100,6 +100,9 @@ __remquo (double x, double y, int *quo) *quo = qs ? -cquo : cquo; + /* Ensure correct sign of zero result in round-downward mode. */ + if (x == 0.0) + x = 0.0; if (sx) x = -x; return x; diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_round.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_round.c index 684858c065..565a4622d4 100644 --- a/sysdeps/ieee754/dbl-64/wordsize-64/s_round.c +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_round.c @@ -1,5 +1,5 @@ /* Round double to integer away from zero. - Copyright (C) 1997-2014 Free Software Foundation, Inc. + Copyright (C) 1997-2015 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -32,7 +32,7 @@ __round (double x) EXTRACT_WORDS64 (i0, x); j0 = ((i0 >> 52) & 0x7ff) - 0x3ff; - if (__builtin_expect (j0 < 52, 1)) + if (__glibc_likely (j0 < 52)) { if (j0 < 0) { diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_scalbln.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_scalbln.c index c00db68ddf..8dce51e928 100644 --- a/sysdeps/ieee754/dbl-64/wordsize-64/s_scalbln.c +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_scalbln.c @@ -55,8 +55,6 @@ __scalbln (double x, long int n) INSERT_WORDS64(x,(ix&INT64_C(0x800fffffffffffff))|(k<<52)); return x*twom54; } -weak_alias (__scalbln, scalbln) #ifdef NO_LONG_DOUBLE strong_alias (__scalbln, __scalblnl) -weak_alias (__scalbln, scalblnl) #endif diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_trunc.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_trunc.c index 768e36fdd4..3810f8134d 100644 --- a/sysdeps/ieee754/dbl-64/wordsize-64/s_trunc.c +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_trunc.c @@ -1,5 +1,5 @@ /* Truncate argument to nearest integral value not larger than the argument. - Copyright (C) 1997-2014 Free Software Foundation, Inc. + Copyright (C) 1997-2015 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. diff --git a/sysdeps/ieee754/dbl-64/x2y2m1.c b/sysdeps/ieee754/dbl-64/x2y2m1.c index d39738fc20..c96dae59e1 100644 --- a/sysdeps/ieee754/dbl-64/x2y2m1.c +++ b/sysdeps/ieee754/dbl-64/x2y2m1.c @@ -1,5 +1,5 @@ /* Compute x^2 + y^2 - 1, without large cancellation error. - Copyright (C) 2012-2014 Free Software Foundation, Inc. + Copyright (C) 2012-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 diff --git a/sysdeps/ieee754/dbl-64/x2y2m1f.c b/sysdeps/ieee754/dbl-64/x2y2m1f.c index fd80b9f512..43a8acf8a8 100644 --- a/sysdeps/ieee754/dbl-64/x2y2m1f.c +++ b/sysdeps/ieee754/dbl-64/x2y2m1f.c @@ -1,5 +1,5 @@ /* Compute x^2 + y^2 - 1, without large cancellation error. - Copyright (C) 2012-2014 Free Software Foundation, Inc. + Copyright (C) 2012-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 |