diff options
author | dim <dim@FreeBSD.org> | 2017-09-26 09:01:56 +0000 |
---|---|---|
committer | Luiz Souza <luiz@netgate.com> | 2018-02-24 03:14:30 -0300 |
commit | e121cea29ef9da649959dcfcfd6e42b212b04136 (patch) | |
tree | 56231c03d772ba0cd91ec0bba5b100e7637e22f6 /lib/msun/src | |
parent | b8bedcfda33a0e7fab7a9fc796a47021422aa2b9 (diff) | |
download | FreeBSD-src-e121cea29ef9da649959dcfcfd6e42b212b04136.zip FreeBSD-src-e121cea29ef9da649959dcfcfd6e42b212b04136.tar.gz |
Synchronize most of libm with head as of r323004. This excludes a few
arch-specific updates for powerpcspe, mips and riscv, for which support
has not been merged yet.
Bump __FreeBSD_version for the addition of cacoshl, cacosl, casinhl,
casinl, catanl, catanhl, sincos, sincosf, and sincosl.
MFC r305382 (by bde):
Add asm versions of fmod(), fmodf() and fmodl() on amd64. Add asm
versions of fmodf() amd fmodl() on i387.
fmod is similar to remainder, and the C versions are 3 to 9 times
slower than the asm versions on x86 for both, but we had the strange
mixture of all 6 variants of remainder in asm and only 1 of 6
variants of fmod in asm.
MFC r305384 (by bde):
Disconnect the "optimized" asm variants of cos(), sin() and tan() from
the build on i386. Leave them in the source tree for regression tests.
The asm functions were always much less accurate (by a factor of more
than 10**18 in the worst case). They were faster on old CPUs. But
with each new generation of CPUs they get relatively slower. The
double precision C version's average advantage is about a factor of 2
on Haswell.
The asm functions were already intentionally avoided in float and long
double precision on i386 and in all precisions on amd64. Float
precision and amd64 give larger advantages to the C version. The long
double precision C code and compilers' understanding of long double
precision are not so good, so the i387 is still slightly faster for
long double precision, except for the unimportant subcase of huge args
where the sub-optimal C code now somehow beats the i387 by about a
factor of 2.
MFC r305385 (by bde):
Oops, the previous i386 version of e_fmodf.S and e_fmodl.S was
actually the amd64 version.
MFC r306409 (by emaste):
libm: fix some unused variable (rcsid) and dangling else warnings
s_{fabs,fmax,logb,scalb}{,f,l}.c may be built elsewhere with a higher
WARNS setting.
Reviewed by: ed
Sponsored by: The FreeBSD Foundation
Differential Revision: https://reviews.freebsd.org/D8061
MFC r306410 (by emaste):
libm: simplify i387 subdir logic with make's :S substitution
MFC r306527 (by emaste):
libm: remove unused variables for LDBL_MANT_DIG != 113
Sponsored by: The FreeBSD Foundation
MFC r306709 (by emaste):
libm: remove unused variables
Sponsored by: The FreeBSD Foundation
MFC r307066 (by br):
Don't use fmaxl/fminl on platforms with no long double support,
use fmax/fmin instead.
This fixes fmaxmin test failure on MIPS64.
Reviewed by: emaste
Sponsored by: DARPA, AFRL
Sponsored by: HEIF5
Differential Revision: https://reviews.freebsd.org/D8216
MFC r308172 (by emaste):
libm: add braces around initialization of subobjects
This cleans up a warning when building libm at higher WARNS levels and
makes the intent more clear. By the C standard the values are assigned
to subobject members in order so this change introduces no functional
change. (6.7.9 20)
Reviewed by: kib
Sponsored by: The FreeBSD Foundation
Differential Revision: https://reviews.freebsd.org/D8333
MFC r313761 (by mmokhi):
Add casinl() cacosl() catanl() casinhl() cacoshl() catanhl() APIs to msun
to improve C11 conformance.
PR: 216850 216851 216852 216856 216857 216858
Submitted by: mmokhi
Reported by: sgk@troutmask.apl.washington.edu
Reviewed by: bde, mat, theraven
Approved by: bde (src committer), mat (mentor)
Differential Revision: https://reviews.freebsd.org/D9491
MFC r313863 (by mmokhi):
Fix building of r313761 on platforms that
`long double` is alias of `double` (MIPS, etc)
PR: 216850 216851 216852 216856 216857 216858
Reported by: emsate
Reviewed by: bde emaste hselasky
Approved by: bde emaste hselasky
Differential Revision: https://reviews.freebsd.org/D9491
MFC r313864 (by mmokhi):
Add documentations related to new APIs of r313761
PR: 216850 216851 216852 216856 216857 216858
Submitted by: sgk@troutmask.apl.washington.edu
Reported by: sgk@troutmask.apl.washington.edu
Reviewed by: bde emaste hselasky
Approved by: bde emaste hselasky
Differential Revision: https://reviews.freebsd.org/D9491
MFC r314950 (by ngie):
Don't expect :test_large_inputs to fail with i386 anymore
Recent changes (maybe a side-effect of the ATF-ification in r314649)
invalidate the failure expectation.
PR: 205446
Sponsored by: Dell EMC Isilon
MFC r317349 (by pfg):
msun: Remove trailing space in Sunsoft copyright statement.
Submittedby: kargl
MFC r319047 (by mmel):
Implement sincos, sincosf, and sincosl.
The primary benefit of these functions is that argument
reduction is done once instead of twice in independent
calls to sin() and cos().
* lib/msun/Makefile:
. Add s_sincos[fl].c to the build.
. Add sincos.3 documentation.
. Add appropriate MLINKS.
* lib/msun/Symbol.map:
. Expose sincos[fl] symbols in dynamic libm.so.
* lib/msun/man/sincos.3:
. Documentation for sincos[fl].
* lib/msun/src/k_sincos.h:
. Kernel for sincos() function. This merges the individual kernels
for sin() and cos(). The merger offered an opportunity to re-arrange
the individual kernels for better performance.
* lib/msun/src/k_sincosf.h:
. Kernel for sincosf() function. This merges the individual kernels
for sinf() and cosf(). The merger offered an opportunity to re-arrange
the individual kernels for better performance.
* lib/msun/src/k_sincosl.h:
. Kernel for sincosl() function. This merges the individual kernels
for sinl() and cosl(). The merger offered an opportunity to re-arrange
the individual kernels for better performance.
* lib/msun/src/math.h:
. Add prototytpes for sincos[fl]().
* lib/msun/src/math_private.h:
. Add RETURNV macros. This is needed to reset fpsetprec on I386
hardware for a function with type void.
* lib/msun/src/s_sincos.c:
. Implementation of sincos() where sin() and cos() were merged into
one routine and possibly re-arranged for better performance.
* lib/msun/src/s_sincosf.c:
. Implementation of sincosf() where sinf() and cosf() were merged into
one routine and possibly re-arranged for better performance.
* lib/msun/src/s_sincosl.c:
. Implementation of sincosl() where sinl() and cosl() were merged into
one routine and possibly re-arranged for better performance.
PR: 215977, 218300
Submitted by: Steven G. Kargl <sgk@troutmask.apl.washington.edu>
Differential Revision: https://reviews.freebsd.org/D10765
MFC r321457 (by ngie):
Mark :reduction as an expected failure
It fails with clang 5.0+.
PR: 220989
Reported by: Jenkins
MFC r322418 (by rlibby):
lib/msun: avoid referring to broken LDBL_MAX
LDBL_MAX is broken on i386:
https://lists.freebsd.org/pipermail/freebsd-numerics/2012-September/000288.html
Gcc has produced +Infinity for LDBL_MAX on i386 and amd64 with -m32
for some time, and newer versions of gcc are now warning that the
"floating constant exceeds range of 'long double'". Avoid this by
referring to half the value of LDBL_MAX instead.
Reviewed by: bde
Approved by: markj (mentor)
Sponsored by: Dell EMC Isilon
MFC r322435 (by rlibby):
Revert r322418, LDBL_MAX_EXP unsuitable for macro pasting on some arches
Either need a different way to spell HALF_LDBL_MAX, or a different way
to spell LDBL_MAX_EXP, or a different approach.
Reported by: ian
MFC r322921 (by ngie):
Revert r321457
It doesn't fail after ^/head@r322855 (the releng_50 clang merge).
PR: 220989
(cherry picked from commit 68896d08c13fabea18690d4c5369de8a5e24ba82)
Diffstat (limited to 'lib/msun/src')
-rw-r--r-- | lib/msun/src/catrig.c | 9 | ||||
-rw-r--r-- | lib/msun/src/catrigl.c | 417 | ||||
-rw-r--r-- | lib/msun/src/e_asin.c | 2 | ||||
-rw-r--r-- | lib/msun/src/e_coshl.c | 2 | ||||
-rw-r--r-- | lib/msun/src/e_lgammaf_r.c | 2 | ||||
-rw-r--r-- | lib/msun/src/e_sinhl.c | 5 | ||||
-rw-r--r-- | lib/msun/src/k_sincos.h | 52 | ||||
-rw-r--r-- | lib/msun/src/k_sincosf.h | 43 | ||||
-rw-r--r-- | lib/msun/src/k_sincosl.h | 134 | ||||
-rw-r--r-- | lib/msun/src/math.h | 3 | ||||
-rw-r--r-- | lib/msun/src/math_private.h | 14 | ||||
-rw-r--r-- | lib/msun/src/s_fabs.c | 5 | ||||
-rw-r--r-- | lib/msun/src/s_fmax.c | 5 | ||||
-rw-r--r-- | lib/msun/src/s_fmin.c | 5 | ||||
-rw-r--r-- | lib/msun/src/s_logbl.c | 5 | ||||
-rw-r--r-- | lib/msun/src/s_scalbn.c | 12 | ||||
-rw-r--r-- | lib/msun/src/s_scalbnf.c | 11 | ||||
-rw-r--r-- | lib/msun/src/s_scalbnl.c | 12 | ||||
-rw-r--r-- | lib/msun/src/s_sincos.c | 80 | ||||
-rw-r--r-- | lib/msun/src/s_sincosf.c | 126 | ||||
-rw-r--r-- | lib/msun/src/s_sincosl.c | 105 | ||||
-rw-r--r-- | lib/msun/src/s_tanhl.c | 2 |
22 files changed, 1023 insertions, 28 deletions
diff --git a/lib/msun/src/catrig.c b/lib/msun/src/catrig.c index f392862..ebfe356 100644 --- a/lib/msun/src/catrig.c +++ b/lib/msun/src/catrig.c @@ -637,3 +637,12 @@ catan(double complex z) return (CMPLX(cimag(w), creal(w))); } + +#if LDBL_MANT_DIG == 53 +__weak_reference(cacosh, cacoshl); +__weak_reference(cacos, cacosl); +__weak_reference(casinh, casinhl); +__weak_reference(casin, casinl); +__weak_reference(catanh, catanhl); +__weak_reference(catan, catanl); +#endif diff --git a/lib/msun/src/catrigl.c b/lib/msun/src/catrigl.c new file mode 100644 index 0000000..b798570 --- /dev/null +++ b/lib/msun/src/catrigl.c @@ -0,0 +1,417 @@ +/*- + * Copyright (c) 2012 Stephen Montgomery-Smith <stephen@FreeBSD.ORG> + * Copyright (c) 2017 Mahdi Mokhtari <mmokhi@FreeBSD.org> + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * 1. Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND + * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL + * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS + * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) + * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT + * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY + * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF + * SUCH DAMAGE. + */ + +/* + * The algorithm is very close to that in "Implementing the complex arcsine + * and arccosine functions using exception handling" by T. E. Hull, Thomas F. + * Fairgrieve, and Ping Tak Peter Tang, published in ACM Transactions on + * Mathematical Software, Volume 23 Issue 3, 1997, Pages 299-335, + * http://dl.acm.org/citation.cfm?id=275324. + * + * See catrig.c for complete comments. + * + * XXX comments were removed automatically, and even short ones on the right + * of statements were removed (all of them), contrary to normal style. Only + * a few comments on the right of declarations remain. + */ + +#include <sys/cdefs.h> +__FBSDID("$FreeBSD$"); + +#include <complex.h> +#include <float.h> + +#include "invtrig.h" +#include "math.h" +#include "math_private.h" + +#undef isinf +#define isinf(x) (fabsl(x) == INFINITY) +#undef isnan +#define isnan(x) ((x) != (x)) +#define raise_inexact() do { volatile float junk __unused = 1 + tiny; } while(0) +#undef signbit +#define signbit(x) (__builtin_signbitl(x)) + +#if LDBL_MAX_EXP != 0x4000 +#error "Unsupported long double format" +#endif + +static const long double +A_crossover = 10, +B_crossover = 0.6417, +FOUR_SQRT_MIN = 0x1p-8189L, +HALF_MAX = 0x1p16383L, +QUARTER_SQRT_MAX = 0x1p8189L, +RECIP_EPSILON = 1 / LDBL_EPSILON, +SQRT_MIN = 0x1p-8191L; + +#if LDBL_MANT_DIG == 64 +static const union IEEEl2bits +um_e = LD80C(0xadf85458a2bb4a9b, 1, 2.71828182845904523536e+0L), +um_ln2 = LD80C(0xb17217f7d1cf79ac, -1, 6.93147180559945309417e-1L); +#define m_e um_e.e +#define m_ln2 um_ln2.e +static const long double +/* The next 2 literals for non-i386. Misrounding them on i386 is harmless. */ +SQRT_3_EPSILON = 5.70316273435758915310e-10, /* 0x9cc470a0490973e8.0p-94 */ +SQRT_6_EPSILON = 8.06549008734932771664e-10; /* 0xddb3d742c265539e.0p-94 */ +#elif LDBL_MANT_DIG == 113 +static const long double +m_e = 2.71828182845904523536028747135266250e0L, /* 0x15bf0a8b1457695355fb8ac404e7a.0p-111 */ +m_ln2 = 6.93147180559945309417232121458176568e-1L, /* 0x162e42fefa39ef35793c7673007e6.0p-113 */ +SQRT_3_EPSILON = 2.40370335797945490975336727199878124e-17, /* 0x1bb67ae8584caa73b25742d7078b8.0p-168 */ +SQRT_6_EPSILON = 3.39934988877629587239082586223300391e-17; /* 0x13988e1409212e7d0321914321a55.0p-167 */ +#else +#error "Unsupported long double format" +#endif + +static const volatile float +tiny = 0x1p-100; + +static long double complex clog_for_large_values(long double complex z); + +static inline long double +f(long double a, long double b, long double hypot_a_b) +{ + if (b < 0) + return ((hypot_a_b - b) / 2); + if (b == 0) + return (a / 2); + return (a * a / (hypot_a_b + b) / 2); +} + +static inline void +do_hard_work(long double x, long double y, long double *rx, int *B_is_usable, + long double *B, long double *sqrt_A2my2, long double *new_y) +{ + long double R, S, A; + long double Am1, Amy; + + R = hypotl(x, y + 1); + S = hypotl(x, y - 1); + + A = (R + S) / 2; + if (A < 1) + A = 1; + + if (A < A_crossover) { + if (y == 1 && x < LDBL_EPSILON * LDBL_EPSILON / 128) { + *rx = sqrtl(x); + } else if (x >= LDBL_EPSILON * fabsl(y - 1)) { + Am1 = f(x, 1 + y, R) + f(x, 1 - y, S); + *rx = log1pl(Am1 + sqrtl(Am1 * (A + 1))); + } else if (y < 1) { + *rx = x / sqrtl((1 - y) * (1 + y)); + } else { + *rx = log1pl((y - 1) + sqrtl((y - 1) * (y + 1))); + } + } else { + *rx = logl(A + sqrtl(A * A - 1)); + } + + *new_y = y; + + if (y < FOUR_SQRT_MIN) { + *B_is_usable = 0; + *sqrt_A2my2 = A * (2 / LDBL_EPSILON); + *new_y = y * (2 / LDBL_EPSILON); + return; + } + + *B = y / A; + *B_is_usable = 1; + + if (*B > B_crossover) { + *B_is_usable = 0; + if (y == 1 && x < LDBL_EPSILON / 128) { + *sqrt_A2my2 = sqrtl(x) * sqrtl((A + y) / 2); + } else if (x >= LDBL_EPSILON * fabsl(y - 1)) { + Amy = f(x, y + 1, R) + f(x, y - 1, S); + *sqrt_A2my2 = sqrtl(Amy * (A + y)); + } else if (y > 1) { + *sqrt_A2my2 = x * (4 / LDBL_EPSILON / LDBL_EPSILON) * y / + sqrtl((y + 1) * (y - 1)); + *new_y = y * (4 / LDBL_EPSILON / LDBL_EPSILON); + } else { + *sqrt_A2my2 = sqrtl((1 - y) * (1 + y)); + } + } +} + +long double complex +casinhl(long double complex z) +{ + long double x, y, ax, ay, rx, ry, B, sqrt_A2my2, new_y; + int B_is_usable; + long double complex w; + + x = creall(z); + y = cimagl(z); + ax = fabsl(x); + ay = fabsl(y); + + if (isnan(x) || isnan(y)) { + if (isinf(x)) + return (CMPLXL(x, y + y)); + if (isinf(y)) + return (CMPLXL(y, x + x)); + if (y == 0) + return (CMPLXL(x + x, y)); + return (CMPLXL(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); + } + + if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) { + if (signbit(x) == 0) + w = clog_for_large_values(z) + m_ln2; + else + w = clog_for_large_values(-z) + m_ln2; + return (CMPLXL(copysignl(creall(w), x), + copysignl(cimagl(w), y))); + } + + if (x == 0 && y == 0) + return (z); + + raise_inexact(); + + if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4) + return (z); + + do_hard_work(ax, ay, &rx, &B_is_usable, &B, &sqrt_A2my2, &new_y); + if (B_is_usable) + ry = asinl(B); + else + ry = atan2l(new_y, sqrt_A2my2); + return (CMPLXL(copysignl(rx, x), copysignl(ry, y))); +} + +long double complex +casinl(long double complex z) +{ + long double complex w; + + w = casinhl(CMPLXL(cimagl(z), creall(z))); + return (CMPLXL(cimagl(w), creall(w))); +} + +long double complex +cacosl(long double complex z) +{ + long double x, y, ax, ay, rx, ry, B, sqrt_A2mx2, new_x; + int sx, sy; + int B_is_usable; + long double complex w; + + x = creall(z); + y = cimagl(z); + sx = signbit(x); + sy = signbit(y); + ax = fabsl(x); + ay = fabsl(y); + + if (isnan(x) || isnan(y)) { + if (isinf(x)) + return (CMPLXL(y + y, -INFINITY)); + if (isinf(y)) + return (CMPLXL(x + x, -y)); + if (x == 0) + return (CMPLXL(pio2_hi + pio2_lo, y + y)); + return (CMPLXL(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); + } + + if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) { + w = clog_for_large_values(z); + rx = fabsl(cimagl(w)); + ry = creall(w) + m_ln2; + if (sy == 0) + ry = -ry; + return (CMPLXL(rx, ry)); + } + + if (x == 1 && y == 0) + return (CMPLXL(0, -y)); + + raise_inexact(); + + if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4) + return (CMPLXL(pio2_hi - (x - pio2_lo), -y)); + + do_hard_work(ay, ax, &ry, &B_is_usable, &B, &sqrt_A2mx2, &new_x); + if (B_is_usable) { + if (sx == 0) + rx = acosl(B); + else + rx = acosl(-B); + } else { + if (sx == 0) + rx = atan2l(sqrt_A2mx2, new_x); + else + rx = atan2l(sqrt_A2mx2, -new_x); + } + if (sy == 0) + ry = -ry; + return (CMPLXL(rx, ry)); +} + +long double complex +cacoshl(long double complex z) +{ + long double complex w; + long double rx, ry; + + w = cacosl(z); + rx = creall(w); + ry = cimagl(w); + if (isnan(rx) && isnan(ry)) + return (CMPLXL(ry, rx)); + if (isnan(rx)) + return (CMPLXL(fabsl(ry), rx)); + if (isnan(ry)) + return (CMPLXL(ry, ry)); + return (CMPLXL(fabsl(ry), copysignl(rx, cimagl(z)))); +} + +static long double complex +clog_for_large_values(long double complex z) +{ + long double x, y; + long double ax, ay, t; + + x = creall(z); + y = cimagl(z); + ax = fabsl(x); + ay = fabsl(y); + if (ax < ay) { + t = ax; + ax = ay; + ay = t; + } + + if (ax > HALF_MAX) + return (CMPLXL(logl(hypotl(x / m_e, y / m_e)) + 1, + atan2l(y, x))); + + if (ax > QUARTER_SQRT_MAX || ay < SQRT_MIN) + return (CMPLXL(logl(hypotl(x, y)), atan2l(y, x))); + + return (CMPLXL(logl(ax * ax + ay * ay) / 2, atan2l(y, x))); +} + +static inline long double +sum_squares(long double x, long double y) +{ + + if (y < SQRT_MIN) + return (x * x); + + return (x * x + y * y); +} + +static inline long double +real_part_reciprocal(long double x, long double y) +{ + long double scale; + uint16_t hx, hy; + int16_t ix, iy; + + GET_LDBL_EXPSIGN(hx, x); + ix = hx & 0x7fff; + GET_LDBL_EXPSIGN(hy, y); + iy = hy & 0x7fff; +#define BIAS (LDBL_MAX_EXP - 1) +#define CUTOFF (LDBL_MANT_DIG / 2 + 1) + if (ix - iy >= CUTOFF || isinf(x)) + return (1 / x); + if (iy - ix >= CUTOFF) + return (x / y / y); + if (ix <= BIAS + LDBL_MAX_EXP / 2 - CUTOFF) + return (x / (x * x + y * y)); + scale = 1; + SET_LDBL_EXPSIGN(scale, 0x7fff - ix); + x *= scale; + y *= scale; + return (x / (x * x + y * y) * scale); +} + +long double complex +catanhl(long double complex z) +{ + long double x, y, ax, ay, rx, ry; + + x = creall(z); + y = cimagl(z); + ax = fabsl(x); + ay = fabsl(y); + + if (y == 0 && ax <= 1) + return (CMPLXL(atanhl(x), y)); + + if (x == 0) + return (CMPLXL(x, atanl(y))); + + if (isnan(x) || isnan(y)) { + if (isinf(x)) + return (CMPLXL(copysignl(0, x), y + y)); + if (isinf(y)) + return (CMPLXL(copysignl(0, x), + copysignl(pio2_hi + pio2_lo, y))); + return (CMPLXL(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); + } + + if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) + return (CMPLXL(real_part_reciprocal(x, y), + copysignl(pio2_hi + pio2_lo, y))); + + if (ax < SQRT_3_EPSILON / 2 && ay < SQRT_3_EPSILON / 2) { + raise_inexact(); + return (z); + } + + if (ax == 1 && ay < LDBL_EPSILON) + rx = (m_ln2 - logl(ay)) / 2; + else + rx = log1pl(4 * ax / sum_squares(ax - 1, ay)) / 4; + + if (ax == 1) + ry = atan2l(2, -ay) / 2; + else if (ay < LDBL_EPSILON) + ry = atan2l(2 * ay, (1 - ax) * (1 + ax)) / 2; + else + ry = atan2l(2 * ay, (1 - ax) * (1 + ax) - ay * ay) / 2; + + return (CMPLXL(copysignl(rx, x), copysignl(ry, y))); +} + +long double complex +catanl(long double complex z) +{ + long double complex w; + + w = catanhl(CMPLXL(cimagl(z), creall(z))); + return (CMPLXL(cimagl(w), creall(w))); +} diff --git a/lib/msun/src/e_asin.c b/lib/msun/src/e_asin.c index 27de207..931b270 100644 --- a/lib/msun/src/e_asin.c +++ b/lib/msun/src/e_asin.c @@ -6,7 +6,7 @@ * * Developed at SunSoft, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice + * software is freely granted, provided that this notice * is preserved. * ==================================================== */ diff --git a/lib/msun/src/e_coshl.c b/lib/msun/src/e_coshl.c index 0a21277..4e3b283 100644 --- a/lib/msun/src/e_coshl.c +++ b/lib/msun/src/e_coshl.c @@ -86,7 +86,9 @@ long double coshl(long double x) { long double hi,lo,x2,x4; +#if LDBL_MANT_DIG == 113 double dx2; +#endif uint16_t ix; GET_LDBL_EXPSIGN(ix,x); diff --git a/lib/msun/src/e_lgammaf_r.c b/lib/msun/src/e_lgammaf_r.c index 9084e18..48346c33 100644 --- a/lib/msun/src/e_lgammaf_r.c +++ b/lib/msun/src/e_lgammaf_r.c @@ -122,7 +122,7 @@ sin_pif(float x) float __ieee754_lgammaf_r(float x, int *signgamp) { - float nadj,p,p1,p2,p3,q,r,t,w,y,z; + float nadj,p,p1,p2,q,r,t,w,y,z; int32_t hx; int i,ix; diff --git a/lib/msun/src/e_sinhl.c b/lib/msun/src/e_sinhl.c index ce7e333..38d3df1 100644 --- a/lib/msun/src/e_sinhl.c +++ b/lib/msun/src/e_sinhl.c @@ -85,7 +85,10 @@ long double sinhl(long double x) { long double hi,lo,x2,x4; - double dx2,s; +#if LDBL_MANT_DIG == 113 + double dx2; +#endif + double s; int16_t ix,jx; GET_LDBL_EXPSIGN(jx,x); diff --git a/lib/msun/src/k_sincos.h b/lib/msun/src/k_sincos.h new file mode 100644 index 0000000..6f03be2 --- /dev/null +++ b/lib/msun/src/k_sincos.h @@ -0,0 +1,52 @@ +/*- + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunSoft, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + * + * k_sin.c and k_cos.c merged by Steven G. Kargl. + */ + +#include <sys/cdefs.h> +__FBSDID("$FreeBSD$"); + +static const double +S1 = -1.66666666666666324348e-01, /* 0xBFC55555, 0x55555549 */ +S2 = 8.33333333332248946124e-03, /* 0x3F811111, 0x1110F8A6 */ +S3 = -1.98412698298579493134e-04, /* 0xBF2A01A0, 0x19C161D5 */ +S4 = 2.75573137070700676789e-06, /* 0x3EC71DE3, 0x57B1FE7D */ +S5 = -2.50507602534068634195e-08, /* 0xBE5AE5E6, 0x8A2B9CEB */ +S6 = 1.58969099521155010221e-10; /* 0x3DE5D93A, 0x5ACFD57C */ + +static const double +C1 = 4.16666666666666019037e-02, /* 0x3FA55555, 0x5555554C */ +C2 = -1.38888888888741095749e-03, /* 0xBF56C16C, 0x16C15177 */ +C3 = 2.48015872894767294178e-05, /* 0x3EFA01A0, 0x19CB1590 */ +C4 = -2.75573143513906633035e-07, /* 0xBE927E4F, 0x809C52AD */ +C5 = 2.08757232129817482790e-09, /* 0x3E21EE9E, 0xBDB4B1C4 */ +C6 = -1.13596475577881948265e-11; /* 0xBDA8FAE9, 0xBE8838D4 */ + +static inline void +__kernel_sincos(double x, double y, int iy, double *sn, double *cs) +{ + double hz, r, v, w, z; + + z = x * x; + w = z * z; + r = S2 + z * (S3 + z * S4) + z * w * (S5 + z * S6); + v = z * x; + + if (iy == 0) + *sn = x + v * (S1 + z * r); + else + *sn = x - ((z * (y / 2 - v * r) - y) - v * S1); + + r = z * (C1 + z * (C2 + z * C3)) + w * w * (C4 + z * (C5 + z * C6)); + hz = z / 2; + w = 1 - hz; + *cs = w + (((1 - w) - hz) + (z * r - x * y)); +} diff --git a/lib/msun/src/k_sincosf.h b/lib/msun/src/k_sincosf.h new file mode 100644 index 0000000..073986d --- /dev/null +++ b/lib/msun/src/k_sincosf.h @@ -0,0 +1,43 @@ +/*- + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + * + * k_sinf.c and k_cosf.c merged by Steven G. Kargl. + */ + +#include <sys/cdefs.h> +__FBSDID("$FreeBSD$"); + +/* |sin(x)/x - s(x)| < 2**-37.5 (~[-4.89e-12, 4.824e-12]). */ +static const double +S1 = -0x15555554cbac77.0p-55, /* -0.166666666416265235595 */ +S2 = 0x111110896efbb2.0p-59, /* 0.0083333293858894631756 */ +S3 = -0x1a00f9e2cae774.0p-65, /* -0.000198393348360966317347 */ +S4 = 0x16cd878c3b46a7.0p-71; /* 0.0000027183114939898219064 */ + +/* |cos(x) - c(x)| < 2**-34.1 (~[-5.37e-11, 5.295e-11]). */ +static const double +C0 = -0x1ffffffd0c5e81.0p-54, /* -0.499999997251031003120 */ +C1 = 0x155553e1053a42.0p-57, /* 0.0416666233237390631894 */ +C2 = -0x16c087e80f1e27.0p-62, /* -0.00138867637746099294692 */ +C3 = 0x199342e0ee5069.0p-68; /* 0.0000243904487962774090654 */ + +static inline void +__kernel_sincosdf(double x, float *sn, float *cs) +{ + double r, s, w, z; + + z = x * x; + w = z * z; + r = S3 + z * S4; + s = z * x; + *sn = (x + s * (S1 + z * S2)) + s * w * r; + r = C2 + z * C3; + *cs = ((1 + z * C0) + w * C1) + (w * z) * r; +} diff --git a/lib/msun/src/k_sincosl.h b/lib/msun/src/k_sincosl.h new file mode 100644 index 0000000..3a609bc --- /dev/null +++ b/lib/msun/src/k_sincosl.h @@ -0,0 +1,134 @@ +/*- + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * Copyright (c) 2008 Steven G. Kargl, David Schultz, Bruce D. Evans. + * + * Developed at SunSoft, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + * + * k_sinl.c and k_cosl.c merged by Steven G. Kargl + */ + +#include <sys/cdefs.h> +__FBSDID("$FreeBSD$"); + +#if LDBL_MANT_DIG == 64 /* ld80 version of k_sincosl.c. */ + +#if defined(__amd64__) || defined(__i386__) +/* Long double constants are slow on these arches, and broken on i386. */ +static const volatile double +C1hi = 0.041666666666666664, /* 0x15555555555555.0p-57 */ +C1lo = 2.2598839032744733e-18, /* 0x14d80000000000.0p-111 */ +S1hi = -0.16666666666666666, /* -0x15555555555555.0p-55 */ +S1lo = -9.2563760475949941e-18; /* -0x15580000000000.0p-109 */ +#define S1 ((long double)S1hi + S1lo) +#define C1 ((long double)C1hi + C1lo) +#else +static const long double +C1 = 0.0416666666666666666136L; /* 0xaaaaaaaaaaaaaa9b.0p-68 */ +S1 = -0.166666666666666666671L, /* -0xaaaaaaaaaaaaaaab.0p-66 */ +#endif + +static const double +C2 = -0.0013888888888888874, /* -0x16c16c16c16c10.0p-62 */ +C3 = 0.000024801587301571716, /* 0x1a01a01a018e22.0p-68 */ +C4 = -0.00000027557319215507120, /* -0x127e4fb7602f22.0p-74 */ +C5 = 0.0000000020876754400407278, /* 0x11eed8caaeccf1.0p-81 */ +C6 = -1.1470297442401303e-11, /* -0x19393412bd1529.0p-89 */ +C7 = 4.7383039476436467e-14, /* 0x1aac9d9af5c43e.0p-97 */ +S2 = 0.0083333333333333332, /* 0x11111111111111.0p-59 */ +S3 = -0.00019841269841269427, /* -0x1a01a01a019f81.0p-65 */ +S4 = 0.0000027557319223597490, /* 0x171de3a55560f7.0p-71 */ +S5 = -0.000000025052108218074604, /* -0x1ae64564f16cad.0p-78 */ +S6 = 1.6059006598854211e-10, /* 0x161242b90243b5.0p-85 */ +S7 = -7.6429779983024564e-13, /* -0x1ae42ebd1b2e00.0p-93 */ +S8 = 2.6174587166648325e-15; /* 0x179372ea0b3f64.0p-101 */ + +static inline void +__kernel_sincosl(long double x, long double y, int iy, long double *sn, + long double *cs) +{ + long double hz, r, v, w, z; + + z = x * x; + v = z * x; + /* + * XXX Replace Horner scheme with an algorithm suitable for CPUs + * with more complex pipelines. + */ + r = S2 + z * (S3 + z * (S4 + z * (S5 + z * (S6 + z * (S7 + z * S8))))); + + if (iy == 0) + *sn = x + v * (S1 + z * r); + else + *sn = x - ((z * (y / 2 - v * r) - y) - v * S1); + + hz = z / 2; + w = 1 - hz; + r = z * (C1 + z * (C2 + z * (C3 + z * (C4 + z * (C5 + z * (C6 + + z * C7)))))); + *cs = w + (((1 - w) - hz) + (z * r - x * y)); +} + +#elif LDBL_MANT_DIG == 113 /* ld128 version of k_sincosl.c. */ + +static const long double +C1 = 0.04166666666666666666666666666666658424671L, +C2 = -0.001388888888888888888888888888863490893732L, +C3 = 0.00002480158730158730158730158600795304914210L, +C4 = -0.2755731922398589065255474947078934284324e-6L, +C5 = 0.2087675698786809897659225313136400793948e-8L, +C6 = -0.1147074559772972315817149986812031204775e-10L, +C7 = 0.4779477332386808976875457937252120293400e-13L, +S1 = -0.16666666666666666666666666666666666606732416116558L, +S2 = 0.0083333333333333333333333333333331135404851288270047L, +S3 = -0.00019841269841269841269841269839935785325638310428717L, +S4 = 0.27557319223985890652557316053039946268333231205686e-5L, +S5 = -0.25052108385441718775048214826384312253862930064745e-7L, +S6 = 0.16059043836821614596571832194524392581082444805729e-9L, +S7 = -0.76471637318198151807063387954939213287488216303768e-12L, +S8 = 0.28114572543451292625024967174638477283187397621303e-14L; + +static const double +C8 = -0.1561920696721507929516718307820958119868e-15, +C9 = 0.4110317413744594971475941557607804508039e-18, +C10 = -0.8896592467191938803288521958313920156409e-21, +C11 = 0.1601061435794535138244346256065192782581e-23, +S9 = -0.82206352458348947812512122163446202498005154296863e-17, +S10 = 0.19572940011906109418080609928334380560135358385256e-19, +S11 = -0.38680813379701966970673724299207480965452616911420e-22, +S12 = 0.64038150078671872796678569586315881020659912139412e-25; + +static inline void +__kernel_sincosl(long double x, long double y, int iy, long double *sn, + long double *cs) +{ + long double hz, r, v, w, z; + + z = x * x; + v = z * x; + /* + * XXX Replace Horner scheme with an algorithm suitable for CPUs + * with more complex pipelines. + */ + r = S2 + z * (S3 + z * (S4 + z * (S5 + z * (S6 + z * (S7 + z * (S8 + + z * (S9 + z * (S10 + z * (S11 + z * S12))))))))); + + if (iy == 0) + *sn = x + v * (S1 + z * r); + else + *cs = x - ((z * (y / 2 - v * r) - y) - v * S1); + + hz = z / 2; + w = 1 - hz; + r = z * (C1 + z * (C2 + z * (C3 + z * (C4 + z * (C5 + z * (C6 + + z * (C7 + z * (C8 + z * (C9 + z * (C10 + z * C11)))))))))); + + *cs = w + (((1 - w) - hz) + (z * r - x * y)); +} +#else +#error "Unsupported long double format" +#endif diff --git a/lib/msun/src/math.h b/lib/msun/src/math.h index 2214c07..a8f4554 100644 --- a/lib/msun/src/math.h +++ b/lib/msun/src/math.h @@ -500,6 +500,9 @@ long double truncl(long double); #if __BSD_VISIBLE long double lgammal_r(long double, int *); +void sincos(double, double *, double *); +void sincosf(float, float *, float *); +void sincosl(long double, long double *, long double *); #endif __END_DECLS diff --git a/lib/msun/src/math_private.h b/lib/msun/src/math_private.h index afaf201..7487902 100644 --- a/lib/msun/src/math_private.h +++ b/lib/msun/src/math_private.h @@ -306,9 +306,21 @@ do { \ fpsetprec(__oprec); \ RETURNF(__retval); \ } while (0) +#define ENTERV() \ + fp_prec_t __oprec; \ + \ + if ((__oprec = fpgetprec()) != FP_PE) \ + fpsetprec(FP_PE) +#define RETURNV() do { \ + if (__oprec != FP_PE) \ + fpsetprec(__oprec); \ + return; \ +} while (0) #else -#define ENTERI(x) +#define ENTERI() #define RETURNI(x) RETURNF(x) +#define ENTERV() +#define RETURNV() return #endif /* Default return statement if hack*_t() is not used. */ diff --git a/lib/msun/src/s_fabs.c b/lib/msun/src/s_fabs.c index 15529e5..48aab25 100644 --- a/lib/msun/src/s_fabs.c +++ b/lib/msun/src/s_fabs.c @@ -10,9 +10,8 @@ * ==================================================== */ -#ifndef lint -static char rcsid[] = "$FreeBSD$"; -#endif +#include <sys/cdefs.h> +__FBSDID("$FreeBSD$"); /* * fabs(x) returns the absolute value of x. diff --git a/lib/msun/src/s_fmax.c b/lib/msun/src/s_fmax.c index b51b865..77e27be 100644 --- a/lib/msun/src/s_fmax.c +++ b/lib/msun/src/s_fmax.c @@ -27,6 +27,7 @@ #include <sys/cdefs.h> __FBSDID("$FreeBSD$"); +#include <float.h> #include <math.h> #include "fpmath.h" @@ -51,3 +52,7 @@ fmax(double x, double y) return (x > y ? x : y); } + +#if (LDBL_MANT_DIG == 53) +__weak_reference(fmax, fmaxl); +#endif diff --git a/lib/msun/src/s_fmin.c b/lib/msun/src/s_fmin.c index 3500c84..5545ff2 100644 --- a/lib/msun/src/s_fmin.c +++ b/lib/msun/src/s_fmin.c @@ -27,6 +27,7 @@ #include <sys/cdefs.h> __FBSDID("$FreeBSD$"); +#include <float.h> #include <math.h> #include "fpmath.h" @@ -51,3 +52,7 @@ fmin(double x, double y) return (x < y ? x : y); } + +#if (LDBL_MANT_DIG == 53) +__weak_reference(fmin, fminl); +#endif diff --git a/lib/msun/src/s_logbl.c b/lib/msun/src/s_logbl.c index 7e88e36..ee1a91f 100644 --- a/lib/msun/src/s_logbl.c +++ b/lib/msun/src/s_logbl.c @@ -10,9 +10,8 @@ * ==================================================== */ -#ifndef lint -static char rcsid[] = "$FreeBSD$"; -#endif +#include <sys/cdefs.h> +__FBSDID("$FreeBSD$"); #include <float.h> #include <limits.h> diff --git a/lib/msun/src/s_scalbn.c b/lib/msun/src/s_scalbn.c index e7efaab..b048b05 100644 --- a/lib/msun/src/s_scalbn.c +++ b/lib/msun/src/s_scalbn.c @@ -10,9 +10,8 @@ * ==================================================== */ -#ifndef lint -static char rcsid[] = "$FreeBSD$"; -#endif +#include <sys/cdefs.h> +__FBSDID("$FreeBSD$"); /* * scalbn (double x, int n) @@ -21,7 +20,6 @@ static char rcsid[] = "$FreeBSD$"; * exponentiation or a multiplication. */ -#include <sys/cdefs.h> #include <float.h> #include "math.h" @@ -51,10 +49,12 @@ scalbn (double x, int n) if (k > 0x7fe) return huge*copysign(huge,x); /* overflow */ if (k > 0) /* normal result */ {SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20)); return x;} - if (k <= -54) + if (k <= -54) { if (n > 50000) /* in case integer overflow in n+k */ return huge*copysign(huge,x); /*overflow*/ - else return tiny*copysign(tiny,x); /*underflow*/ + else + return tiny*copysign(tiny,x); /*underflow*/ + } k += 54; /* subnormal result */ SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20)); return x*twom54; diff --git a/lib/msun/src/s_scalbnf.c b/lib/msun/src/s_scalbnf.c index 7666c74..21d001c 100644 --- a/lib/msun/src/s_scalbnf.c +++ b/lib/msun/src/s_scalbnf.c @@ -13,11 +13,8 @@ * ==================================================== */ -#ifndef lint -static char rcsid[] = "$FreeBSD$"; -#endif - #include <sys/cdefs.h> +__FBSDID("$FreeBSD$"); #include "math.h" #include "math_private.h" @@ -46,10 +43,12 @@ scalbnf (float x, int n) if (k > 0xfe) return huge*copysignf(huge,x); /* overflow */ if (k > 0) /* normal result */ {SET_FLOAT_WORD(x,(ix&0x807fffff)|(k<<23)); return x;} - if (k <= -25) + if (k <= -25) { if (n > 50000) /* in case integer overflow in n+k */ return huge*copysignf(huge,x); /*overflow*/ - else return tiny*copysignf(tiny,x); /*underflow*/ + else + return tiny*copysignf(tiny,x); /*underflow*/ + } k += 25; /* subnormal result */ SET_FLOAT_WORD(x,(ix&0x807fffff)|(k<<23)); return x*twom25; diff --git a/lib/msun/src/s_scalbnl.c b/lib/msun/src/s_scalbnl.c index fc89f8d..28b0cf9 100644 --- a/lib/msun/src/s_scalbnl.c +++ b/lib/msun/src/s_scalbnl.c @@ -10,9 +10,8 @@ * ==================================================== */ -#ifndef lint -static char rcsid[] = "$FreeBSD$"; -#endif +#include <sys/cdefs.h> +__FBSDID("$FreeBSD$"); /* * scalbnl (long double x, int n) @@ -27,7 +26,6 @@ static char rcsid[] = "$FreeBSD$"; * for scalbn(), so we don't use this routine. */ -#include <sys/cdefs.h> #include <float.h> #include <math.h> @@ -59,10 +57,12 @@ scalbnl (long double x, int n) if (k >= 0x7fff) return huge*copysignl(huge,x); /* overflow */ if (k > 0) /* normal result */ {u.bits.exp = k; return u.e;} - if (k <= -128) + if (k <= -128) { if (n > 50000) /* in case integer overflow in n+k */ return huge*copysign(huge,x); /*overflow*/ - else return tiny*copysign(tiny,x); /*underflow*/ + else + return tiny*copysign(tiny,x); /*underflow*/ + } k += 128; /* subnormal result */ u.bits.exp = k; return u.e*0x1p-128; diff --git a/lib/msun/src/s_sincos.c b/lib/msun/src/s_sincos.c new file mode 100644 index 0000000..85e8d74 --- /dev/null +++ b/lib/msun/src/s_sincos.c @@ -0,0 +1,80 @@ +/*- + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + * + * s_sin.c and s_cos.c merged by Steven G. Kargl. Descriptions of the + * algorithms are contained in the original files. + */ + +#include <sys/cdefs.h> +__FBSDID("$FreeBSD$"); + +#include <float.h> + +#include "math.h" +#define INLINE_REM_PIO2 +#include "math_private.h" +#include "e_rem_pio2.c" +#include "k_sincos.h" + +void +sincos(double x, double *sn, double *cs) +{ + double y[2]; + int32_t n, ix; + + /* High word of x. */ + GET_HIGH_WORD(ix, x); + + /* |x| ~< pi/4 */ + ix &= 0x7fffffff; + if (ix <= 0x3fe921fb) { + if (ix < 0x3e400000) { /* |x| < 2**-27 */ + if ((int)x == 0) { /* Generate inexact. */ + *sn = x; + *cs = 1; + return; + } + } + __kernel_sincos(x, 0, 0, sn, cs); + return; + } + + /* If x = Inf or NaN, then sin(x) = NaN and cos(x) = NaN. */ + if (ix >= 0x7ff00000) { + *sn = x - x; + *cs = x - x; + return; + } + + /* Argument reduction. */ + n = __ieee754_rem_pio2(x, y); + + switch(n & 3) { + case 0: + __kernel_sincos(y[0], y[1], 1, sn, cs); + break; + case 1: + __kernel_sincos(y[0], y[1], 1, cs, sn); + *cs = -*cs; + break; + case 2: + __kernel_sincos(y[0], y[1], 1, sn, cs); + *sn = -*sn; + *cs = -*cs; + break; + default: + __kernel_sincos(y[0], y[1], 1, cs, sn); + *sn = -*sn; + } +} + +#if (LDBL_MANT_DIG == 53) +__weak_reference(sincos, sincosl); +#endif diff --git a/lib/msun/src/s_sincosf.c b/lib/msun/src/s_sincosf.c new file mode 100644 index 0000000..755ff05 --- /dev/null +++ b/lib/msun/src/s_sincosf.c @@ -0,0 +1,126 @@ +/*- + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* s_sincosf.c -- float version of s_sincos.c. + * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. + * Optimized by Bruce D. Evans. + * Merged s_sinf.c and s_cosf.c by Steven G. Kargl. + */ + +#include <sys/cdefs.h> +__FBSDID("$FreeBSD$"); + +#include <float.h> + +#include "math.h" +#define INLINE_REM_PIO2F +#include "math_private.h" +#include "e_rem_pio2f.c" +#include "k_sincosf.h" + +/* Small multiples of pi/2 rounded to double precision. */ +static const double +p1pio2 = 1*M_PI_2, /* 0x3FF921FB, 0x54442D18 */ +p2pio2 = 2*M_PI_2, /* 0x400921FB, 0x54442D18 */ +p3pio2 = 3*M_PI_2, /* 0x4012D97C, 0x7F3321D2 */ +p4pio2 = 4*M_PI_2; /* 0x401921FB, 0x54442D18 */ + +void +sincosf(float x, float *sn, float *cs) +{ + float c, s; + double y; + int32_t n, hx, ix; + + GET_FLOAT_WORD(hx, x); + ix = hx & 0x7fffffff; + + if (ix <= 0x3f490fda) { /* |x| ~<= pi/4 */ + if (ix < 0x39800000) { /* |x| < 2**-12 */ + if ((int)x == 0) { + *sn = x; /* x with inexact if x != 0 */ + *cs = 1; + return; + } + } + __kernel_sincosdf(x, sn, cs); + return; + } + + if (ix <= 0x407b53d1) { /* |x| ~<= 5*pi/4 */ + if (ix <= 0x4016cbe3) { /* |x| ~<= 3pi/4 */ + if (hx > 0) { + __kernel_sincosdf(x - p1pio2, cs, sn); + *cs = -*cs; + } else { + __kernel_sincosdf(x + p1pio2, cs, sn); + *sn = -*sn; + } + } else { + if (hx > 0) + __kernel_sincosdf(x - p2pio2, sn, cs); + else + __kernel_sincosdf(x + p2pio2, sn, cs); + *sn = -*sn; + *cs = -*cs; + } + return; + } + + if (ix <= 0x40e231d5) { /* |x| ~<= 9*pi/4 */ + if (ix <= 0x40afeddf) { /* |x| ~<= 7*pi/4 */ + if (hx > 0) { + __kernel_sincosdf(x - p3pio2, cs, sn); + *sn = -*sn; + } else { + __kernel_sincosdf(x + p3pio2, cs, sn); + *cs = -*cs; + } + } else { + if (hx > 0) + __kernel_sincosdf(x - p4pio2, sn, cs); + else + __kernel_sincosdf(x + p4pio2, sn, cs); + } + return; + } + + /* If x = Inf or NaN, then sin(x) = NaN and cos(x) = NaN. */ + if (ix >= 0x7f800000) { + *sn = x - x; + *cs = x - x; + return; + } + + /* Argument reduction. */ + n = __ieee754_rem_pio2f(x, &y); + __kernel_sincosdf(y, &s, &c); + + switch(n & 3) { + case 0: + *sn = s; + *cs = c; + break; + case 1: + *sn = c; + *cs = -s; + break; + case 2: + *sn = -s; + *cs = -c; + break; + default: + *sn = -c; + *cs = s; + } +} + + diff --git a/lib/msun/src/s_sincosl.c b/lib/msun/src/s_sincosl.c new file mode 100644 index 0000000..aef36c2 --- /dev/null +++ b/lib/msun/src/s_sincosl.c @@ -0,0 +1,105 @@ +/*- + * Copyright (c) 2007, 2010-2013 Steven G. Kargl + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * 1. Redistributions of source code must retain the above copyright + * notice unmodified, this list of conditions, and the following + * disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR + * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES + * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. + * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, + * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT + * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, + * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY + * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF + * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + * + * s_sinl.c and s_cosl.c merged by Steven G. Kargl. + */ + +#include <sys/cdefs.h> +__FBSDID("$FreeBSD$"); + +#include <float.h> +#ifdef __i386__ +#include <ieeefp.h> +#endif + +#include "math.h" +#include "math_private.h" +#include "k_sincosl.h" + +#if LDBL_MANT_DIG == 64 +#include "../ld80/e_rem_pio2l.h" +#elif LDBL_MANT_DIG == 113 +#include "../ld128/e_rem_pio2l.h" +#else +#error "Unsupported long double format" +#endif + +void +sincosl(long double x, long double *sn, long double *cs) +{ + union IEEEl2bits z; + int e0, sgn; + long double y[2]; + + z.e = x; + sgn = z.bits.sign; + z.bits.sign = 0; + + ENTERV(); + + /* Optimize the case where x is already within range. */ + if (z.e < M_PI_4) { + /* + * If x = +-0 or x is a subnormal number, then sin(x) = x and + * cos(x) = 1. + */ + if (z.bits.exp == 0) { + *sn = x; + *cs = 1; + } else + __kernel_sincosl(x, 0, 0, sn, cs); + RETURNV(); + } + + /* If x = NaN or Inf, then sin(x) and cos(x) are NaN. */ + if (z.bits.exp == 32767) { + *sn = x - x; + *cs = x - x; + RETURNV(); + } + + /* Range reduction. */ + e0 = __ieee754_rem_pio2l(x, y); + + switch (e0 & 3) { + case 0: + __kernel_sincosl(y[0], y[1], 1, sn, cs); + break; + case 1: + __kernel_sincosl(y[0], y[1], 1, cs, sn); + *cs = -*cs; + break; + case 2: + __kernel_sincosl(y[0], y[1], 1, sn, cs); + *sn = -*sn; + *cs = -*cs; + break; + default: + __kernel_sincosl(y[0], y[1], 1, cs, sn); + *sn = -*sn; + } + + RETURNV(); +} diff --git a/lib/msun/src/s_tanhl.c b/lib/msun/src/s_tanhl.c index 886158b..b753186 100644 --- a/lib/msun/src/s_tanhl.c +++ b/lib/msun/src/s_tanhl.c @@ -113,7 +113,9 @@ long double tanhl(long double x) { long double hi,lo,s,x2,x4,z; +#if LDBL_MANT_DIG == 113 double dx2; +#endif int16_t jx,ix; GET_LDBL_EXPSIGN(jx,x); |