Author: mmokhi (ports committer) Date: Wed Feb 15 07:59:54 2017 New Revision: 313761 URL: https://svnweb.freebsd.org/changeset/base/313761
Log: 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: s...@troutmask.apl.washington.edu Reviewed by: bde, mat, theraven Approved by: bde (src committer), mat (mentor) Differential Revision: https://reviews.freebsd.org/D9491 Added: head/lib/msun/src/catrigl.c - copied, changed from r313732, head/lib/msun/src/catrigf.c Modified: head/lib/msun/Makefile head/lib/msun/Symbol.map Modified: head/lib/msun/Makefile ============================================================================== --- head/lib/msun/Makefile Wed Feb 15 07:22:47 2017 (r313760) +++ head/lib/msun/Makefile Wed Feb 15 07:59:54 2017 (r313761) @@ -107,7 +107,7 @@ COMMON_SRCS+= e_acoshl.c e_acosl.c e_asi .endif # C99 complex functions -COMMON_SRCS+= catrig.c catrigf.c \ +COMMON_SRCS+= catrig.c catrigf.c catrigl.c \ s_ccosh.c s_ccoshf.c s_cexp.c s_cexpf.c \ s_cimag.c s_cimagf.c s_cimagl.c \ s_conj.c s_conjf.c s_conjl.c \ Modified: head/lib/msun/Symbol.map ============================================================================== --- head/lib/msun/Symbol.map Wed Feb 15 07:22:47 2017 (r313760) +++ head/lib/msun/Symbol.map Wed Feb 15 07:59:54 2017 (r313761) @@ -285,3 +285,13 @@ FBSD_1.3 { FBSD_1.4 { lgammal_r; }; + +/* First added in 12.0-CURRENT */ +FBSD_1.5 { + cacoshl; + cacosl; + casinhl; + casinl; + catanl; + catanhl; +}; Copied and modified: head/lib/msun/src/catrigl.c (from r313732, head/lib/msun/src/catrigf.c) ============================================================================== --- head/lib/msun/src/catrigf.c Tue Feb 14 16:49:32 2017 (r313732, copy source) +++ head/lib/msun/src/catrigl.c Wed Feb 15 07:59:54 2017 (r313761) @@ -1,5 +1,6 @@ /*- * Copyright (c) 2012 Stephen Montgomery-Smith <step...@freebsd.org> + * Copyright (c) 2017 Mahdi Mokhtari <mmo...@freebsd.org> * All rights reserved. * * Redistribution and use in source and binary forms, with or without @@ -44,38 +45,53 @@ __FBSDID("$FreeBSD$"); #include <complex.h> #include <float.h> +#include "invtrig.h" #include "math.h" #include "math_private.h" #undef isinf -#define isinf(x) (fabsf(x) == INFINITY) +#define isinf(x) (fabsl(x) == INFINITY) #undef isnan #define isnan(x) ((x) != (x)) #define raise_inexact() do { volatile float junk = 1 + tiny; } while(0) #undef signbit -#define signbit(x) (__builtin_signbitf(x)) +#define signbit(x) (__builtin_signbitl(x)) -static const float +static const long double A_crossover = 10, B_crossover = 0.6417, -FOUR_SQRT_MIN = 0x1p-61, -QUARTER_SQRT_MAX = 0x1p61, -m_e = 2.7182818285e0, /* 0xadf854.0p-22 */ -m_ln2 = 6.9314718056e-1, /* 0xb17218.0p-24 */ -pio2_hi = 1.5707962513e0, /* 0xc90fda.0p-23 */ -RECIP_EPSILON = 1 / FLT_EPSILON, -SQRT_3_EPSILON = 5.9801995673e-4, /* 0x9cc471.0p-34 */ -SQRT_6_EPSILON = 8.4572793338e-4, /* 0xddb3d7.0p-34 */ -SQRT_MIN = 0x1p-63; +FOUR_SQRT_MIN = 0x1p-8189L, +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 -pio2_lo = 7.5497899549e-8, /* 0xa22169.0p-47 */ tiny = 0x1p-100; -static float complex clog_for_large_values(float complex z); +static long double complex clog_for_large_values(long double complex z); -static inline float -f(float a, float b, float hypot_a_b) +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); @@ -85,40 +101,40 @@ f(float a, float b, float hypot_a_b) } static inline void -do_hard_work(float x, float y, float *rx, int *B_is_usable, float *B, - float *sqrt_A2my2, float *new_y) +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) { - float R, S, A; - float Am1, Amy; + long double R, S, A; + long double Am1, Amy; - R = hypotf(x, y + 1); - S = hypotf(x, y - 1); + 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 < FLT_EPSILON * FLT_EPSILON / 128) { - *rx = sqrtf(x); - } else if (x >= FLT_EPSILON * fabsf(y - 1)) { + 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 = log1pf(Am1 + sqrtf(Am1 * (A + 1))); + *rx = log1pl(Am1 + sqrtl(Am1 * (A + 1))); } else if (y < 1) { - *rx = x / sqrtf((1 - y) * (1 + y)); + *rx = x / sqrtl((1 - y) * (1 + y)); } else { - *rx = log1pf((y - 1) + sqrtf((y - 1) * (y + 1))); + *rx = log1pl((y - 1) + sqrtl((y - 1) * (y + 1))); } } else { - *rx = logf(A + sqrtf(A * A - 1)); + *rx = logl(A + sqrtl(A * A - 1)); } *new_y = y; if (y < FOUR_SQRT_MIN) { *B_is_usable = 0; - *sqrt_A2my2 = A * (2 / FLT_EPSILON); - *new_y = y * (2 / FLT_EPSILON); + *sqrt_A2my2 = A * (2 / LDBL_EPSILON); + *new_y = y * (2 / LDBL_EPSILON); return; } @@ -127,41 +143,41 @@ do_hard_work(float x, float y, float *rx if (*B > B_crossover) { *B_is_usable = 0; - if (y == 1 && x < FLT_EPSILON / 128) { - *sqrt_A2my2 = sqrtf(x) * sqrtf((A + y) / 2); - } else if (x >= FLT_EPSILON * fabsf(y - 1)) { + 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 = sqrtf(Amy * (A + y)); + *sqrt_A2my2 = sqrtl(Amy * (A + y)); } else if (y > 1) { - *sqrt_A2my2 = x * (4 / FLT_EPSILON / FLT_EPSILON) * y / - sqrtf((y + 1) * (y - 1)); - *new_y = y * (4 / FLT_EPSILON / FLT_EPSILON); + *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 = sqrtf((1 - y) * (1 + y)); + *sqrt_A2my2 = sqrtl((1 - y) * (1 + y)); } } } -float complex -casinhf(float complex z) +long double complex +casinhl(long double complex z) { - float x, y, ax, ay, rx, ry, B, sqrt_A2my2, new_y; + long double x, y, ax, ay, rx, ry, B, sqrt_A2my2, new_y; int B_is_usable; - float complex w; + long double complex w; - x = crealf(z); - y = cimagf(z); - ax = fabsf(x); - ay = fabsf(y); + x = creall(z); + y = cimagl(z); + ax = fabsl(x); + ay = fabsl(y); if (isnan(x) || isnan(y)) { if (isinf(x)) - return (CMPLXF(x, y + y)); + return (CMPLXL(x, y + y)); if (isinf(y)) - return (CMPLXF(y, x + x)); + return (CMPLXL(y, x + x)); if (y == 0) - return (CMPLXF(x + x, y)); - return (CMPLXF(x + 0.0L + (y + 0), x + 0.0L + (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) { @@ -169,8 +185,8 @@ casinhf(float complex z) w = clog_for_large_values(z) + m_ln2; else w = clog_for_large_values(-z) + m_ln2; - return (CMPLXF(copysignf(crealf(w), x), - copysignf(cimagf(w), y))); + return (CMPLXL(copysignl(creall(w), x), + copysignl(cimagl(w), y))); } if (x == 0 && y == 0) @@ -183,125 +199,126 @@ casinhf(float complex z) do_hard_work(ax, ay, &rx, &B_is_usable, &B, &sqrt_A2my2, &new_y); if (B_is_usable) - ry = asinf(B); + ry = asinl(B); else - ry = atan2f(new_y, sqrt_A2my2); - return (CMPLXF(copysignf(rx, x), copysignf(ry, y))); + ry = atan2l(new_y, sqrt_A2my2); + return (CMPLXL(copysignl(rx, x), copysignl(ry, y))); } -float complex -casinf(float complex z) +long double complex +casinl(long double complex z) { - float complex w = casinhf(CMPLXF(cimagf(z), crealf(z))); + long double complex w; - return (CMPLXF(cimagf(w), crealf(w))); + w = casinhl(CMPLXL(cimagl(z), creall(z))); + return (CMPLXL(cimagl(w), creall(w))); } -float complex -cacosf(float complex z) +long double complex +cacosl(long double complex z) { - float x, y, ax, ay, rx, ry, B, sqrt_A2mx2, new_x; + long double x, y, ax, ay, rx, ry, B, sqrt_A2mx2, new_x; int sx, sy; int B_is_usable; - float complex w; + long double complex w; - x = crealf(z); - y = cimagf(z); + x = creall(z); + y = cimagl(z); sx = signbit(x); sy = signbit(y); - ax = fabsf(x); - ay = fabsf(y); + ax = fabsl(x); + ay = fabsl(y); if (isnan(x) || isnan(y)) { if (isinf(x)) - return (CMPLXF(y + y, -INFINITY)); + return (CMPLXL(y + y, -INFINITY)); if (isinf(y)) - return (CMPLXF(x + x, -y)); + return (CMPLXL(x + x, -y)); if (x == 0) - return (CMPLXF(pio2_hi + pio2_lo, y + y)); - return (CMPLXF(x + 0.0L + (y + 0), x + 0.0L + (y + 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 = fabsf(cimagf(w)); - ry = crealf(w) + m_ln2; + rx = fabsl(cimagl(w)); + ry = creall(w) + m_ln2; if (sy == 0) ry = -ry; - return (CMPLXF(rx, ry)); + return (CMPLXL(rx, ry)); } if (x == 1 && y == 0) - return (CMPLXF(0, -y)); + return (CMPLXL(0, -y)); raise_inexact(); if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4) - return (CMPLXF(pio2_hi - (x - pio2_lo), -y)); + 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 = acosf(B); + rx = acosl(B); else - rx = acosf(-B); + rx = acosl(-B); } else { if (sx == 0) - rx = atan2f(sqrt_A2mx2, new_x); + rx = atan2l(sqrt_A2mx2, new_x); else - rx = atan2f(sqrt_A2mx2, -new_x); + rx = atan2l(sqrt_A2mx2, -new_x); } if (sy == 0) ry = -ry; - return (CMPLXF(rx, ry)); + return (CMPLXL(rx, ry)); } -float complex -cacoshf(float complex z) +long double complex +cacoshl(long double complex z) { - float complex w; - float rx, ry; + long double complex w; + long double rx, ry; - w = cacosf(z); - rx = crealf(w); - ry = cimagf(w); + w = cacosl(z); + rx = creall(w); + ry = cimagl(w); if (isnan(rx) && isnan(ry)) - return (CMPLXF(ry, rx)); + return (CMPLXL(ry, rx)); if (isnan(rx)) - return (CMPLXF(fabsf(ry), rx)); + return (CMPLXL(fabsl(ry), rx)); if (isnan(ry)) - return (CMPLXF(ry, ry)); - return (CMPLXF(fabsf(ry), copysignf(rx, cimagf(z)))); + return (CMPLXL(ry, ry)); + return (CMPLXL(fabsl(ry), copysignl(rx, cimagl(z)))); } -static float complex -clog_for_large_values(float complex z) +static long double complex +clog_for_large_values(long double complex z) { - float x, y; - float ax, ay, t; + long double x, y; + long double ax, ay, t; - x = crealf(z); - y = cimagf(z); - ax = fabsf(x); - ay = fabsf(y); + x = creall(z); + y = cimagl(z); + ax = fabsl(x); + ay = fabsl(y); if (ax < ay) { t = ax; ax = ay; ay = t; } - if (ax > FLT_MAX / 2) - return (CMPLXF(logf(hypotf(x / m_e, y / m_e)) + 1, - atan2f(y, x))); + if (ax > LDBL_MAX / 2) + return (CMPLXL(logl(hypotl(x / m_e, y / m_e)) + 1, + atan2l(y, x))); if (ax > QUARTER_SQRT_MAX || ay < SQRT_MIN) - return (CMPLXF(logf(hypotf(x, y)), atan2f(y, x))); + return (CMPLXL(logl(hypotl(x, y)), atan2l(y, x))); - return (CMPLXF(logf(ax * ax + ay * ay) / 2, atan2f(y, x))); + return (CMPLXL(logl(ax * ax + ay * ay) / 2, atan2l(y, x))); } -static inline float -sum_squares(float x, float y) +static inline long double +sum_squares(long double x, long double y) { if (y < SQRT_MIN) @@ -310,84 +327,86 @@ sum_squares(float x, float y) return (x * x + y * y); } -static inline float -real_part_reciprocal(float x, float y) +static inline long double +real_part_reciprocal(long double x, long double y) { - float scale; - uint32_t hx, hy; - int32_t ix, iy; - - GET_FLOAT_WORD(hx, x); - ix = hx & 0x7f800000; - GET_FLOAT_WORD(hy, y); - iy = hy & 0x7f800000; -#define BIAS (FLT_MAX_EXP - 1) -#define CUTOFF (FLT_MANT_DIG / 2 + 1) - if (ix - iy >= CUTOFF << 23 || isinf(x)) + 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 << 23) + if (iy - ix >= CUTOFF) return (x / y / y); - if (ix <= (BIAS + FLT_MAX_EXP / 2 - CUTOFF) << 23) + if (ix <= BIAS + LDBL_MAX_EXP / 2 - CUTOFF) return (x / (x * x + y * y)); - SET_FLOAT_WORD(scale, 0x7f800000 - ix); + scale = 1; + SET_LDBL_EXPSIGN(scale, 0x7fff - ix); x *= scale; y *= scale; return (x / (x * x + y * y) * scale); } -float complex -catanhf(float complex z) +long double complex +catanhl(long double complex z) { - float x, y, ax, ay, rx, ry; + long double x, y, ax, ay, rx, ry; - x = crealf(z); - y = cimagf(z); - ax = fabsf(x); - ay = fabsf(y); + x = creall(z); + y = cimagl(z); + ax = fabsl(x); + ay = fabsl(y); if (y == 0 && ax <= 1) - return (CMPLXF(atanhf(x), y)); + return (CMPLXL(atanhl(x), y)); if (x == 0) - return (CMPLXF(x, atanf(y))); + return (CMPLXL(x, atanl(y))); if (isnan(x) || isnan(y)) { if (isinf(x)) - return (CMPLXF(copysignf(0, x), y + y)); + return (CMPLXL(copysignl(0, x), y + y)); if (isinf(y)) - return (CMPLXF(copysignf(0, x), - copysignf(pio2_hi + pio2_lo, y))); - return (CMPLXF(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); + 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 (CMPLXF(real_part_reciprocal(x, y), - copysignf(pio2_hi + pio2_lo, y))); + 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 < FLT_EPSILON) - rx = (m_ln2 - logf(ay)) / 2; + if (ax == 1 && ay < LDBL_EPSILON) + rx = (m_ln2 - logl(ay)) / 2; else - rx = log1pf(4 * ax / sum_squares(ax - 1, ay)) / 4; + rx = log1pl(4 * ax / sum_squares(ax - 1, ay)) / 4; if (ax == 1) - ry = atan2f(2, -ay) / 2; - else if (ay < FLT_EPSILON) - ry = atan2f(2 * ay, (1 - ax) * (1 + ax)) / 2; + ry = atan2l(2, -ay) / 2; + else if (ay < LDBL_EPSILON) + ry = atan2l(2 * ay, (1 - ax) * (1 + ax)) / 2; else - ry = atan2f(2 * ay, (1 - ax) * (1 + ax) - ay * ay) / 2; + ry = atan2l(2 * ay, (1 - ax) * (1 + ax) - ay * ay) / 2; - return (CMPLXF(copysignf(rx, x), copysignf(ry, y))); + return (CMPLXL(copysignl(rx, x), copysignl(ry, y))); } -float complex -catanf(float complex z) +long double complex +catanl(long double complex z) { - float complex w = catanhf(CMPLXF(cimagf(z), crealf(z))); + long double complex w; - return (CMPLXF(cimagf(w), crealf(w))); + w = catanhl(CMPLXL(cimagl(z), creall(z))); + return (CMPLXL(cimagl(w), creall(w))); } _______________________________________________ svn-src-head@freebsd.org mailing list https://lists.freebsd.org/mailman/listinfo/svn-src-head To unsubscribe, send any mail to "svn-src-head-unsubscr...@freebsd.org"
