From ccd82ae7cc5ce3f708a919a070bffc7572416c62 Mon Sep 17 00:00:00 2001 From: lzm-build <3575188313@qq.com> Date: Fri, 30 Jan 2026 07:18:13 +0800 Subject: [PATCH] Add `f16`, `f80` and `f128` support for `atan` --- lib/c/math.zig | 67 ++- lib/libc/musl/src/math/atan.c | 116 ----- lib/libc/musl/src/math/atanf.c | 94 ---- lib/libc/musl/src/math/atanl.c | 184 ------- lib/libc/musl/src/math/i386/atan.s | 16 - lib/libc/musl/src/math/i386/atanf.s | 18 - lib/libc/musl/src/math/i386/atanl.s | 7 - lib/libc/musl/src/math/x32/atanl.s | 7 - lib/libc/musl/src/math/x86_64/atanl.s | 7 - lib/std/math/atan.zig | 673 +++++++++++++++++++------- src/libs/musl.zig | 8 - src/libs/wasi_libc.zig | 3 - 12 files changed, 544 insertions(+), 656 deletions(-) delete mode 100644 lib/libc/musl/src/math/atan.c delete mode 100644 lib/libc/musl/src/math/atanf.c delete mode 100644 lib/libc/musl/src/math/atanl.c delete mode 100644 lib/libc/musl/src/math/i386/atan.s delete mode 100644 lib/libc/musl/src/math/i386/atanf.s delete mode 100644 lib/libc/musl/src/math/i386/atanl.s delete mode 100644 lib/libc/musl/src/math/x32/atanl.s delete mode 100644 lib/libc/musl/src/math/x86_64/atanl.s diff --git a/lib/c/math.zig b/lib/c/math.zig index 8a6c91be2efc735326a5a56d1fe3452f8837bb2b..a9e7c808a80adab4758c439209d9c8470319d442 100644 --- a/lib/c/math.zig +++ b/lib/c/math.zig @@ -1,4 +1,5 @@ const std = @import("std"); +const math = std.math; const common = @import("common.zig"); const builtin = @import("builtin"); @@ -11,20 +12,20 @@ comptime { @export(&isnanl, .{ .name = "isnanl", .linkage = common.linkage, .visibility = common.visibility }); @export(&isnanl, .{ .name = "__isnanl", .linkage = common.linkage, .visibility = common.visibility }); - @export(&std.math.nan(f64), .{ .name = "__QNAN", .linkage = common.linkage, .visibility = common.visibility }); - @export(&std.math.snan(f64), .{ .name = "__SNAN", .linkage = common.linkage, .visibility = common.visibility }); - @export(&std.math.inf(f64), .{ .name = "__INF", .linkage = common.linkage, .visibility = common.visibility }); - @export(&std.math.floatTrueMin(f64), .{ .name = "__DENORM", .linkage = common.linkage, .visibility = common.visibility }); + @export(&math.nan(f64), .{ .name = "__QNAN", .linkage = common.linkage, .visibility = common.visibility }); + @export(&math.snan(f64), .{ .name = "__SNAN", .linkage = common.linkage, .visibility = common.visibility }); + @export(&math.inf(f64), .{ .name = "__INF", .linkage = common.linkage, .visibility = common.visibility }); + @export(&math.floatTrueMin(f64), .{ .name = "__DENORM", .linkage = common.linkage, .visibility = common.visibility }); - @export(&std.math.nan(f32), .{ .name = "__QNANF", .linkage = common.linkage, .visibility = common.visibility }); - @export(&std.math.snan(f32), .{ .name = "__SNANF", .linkage = common.linkage, .visibility = common.visibility }); - @export(&std.math.inf(f32), .{ .name = "__INFF", .linkage = common.linkage, .visibility = common.visibility }); - @export(&std.math.floatTrueMin(f32), .{ .name = "__DENORMF", .linkage = common.linkage, .visibility = common.visibility }); + @export(&math.nan(f32), .{ .name = "__QNANF", .linkage = common.linkage, .visibility = common.visibility }); + @export(&math.snan(f32), .{ .name = "__SNANF", .linkage = common.linkage, .visibility = common.visibility }); + @export(&math.inf(f32), .{ .name = "__INFF", .linkage = common.linkage, .visibility = common.visibility }); + @export(&math.floatTrueMin(f32), .{ .name = "__DENORMF", .linkage = common.linkage, .visibility = common.visibility }); - @export(&std.math.nan(c_longdouble), .{ .name = "__QNANL", .linkage = common.linkage, .visibility = common.visibility }); - @export(&std.math.snan(c_longdouble), .{ .name = "__SNANL", .linkage = common.linkage, .visibility = common.visibility }); - @export(&std.math.inf(c_longdouble), .{ .name = "__INFL", .linkage = common.linkage, .visibility = common.visibility }); - @export(&std.math.floatTrueMin(c_longdouble), .{ .name = "__DENORML", .linkage = common.linkage, .visibility = common.visibility }); + @export(&math.nan(c_longdouble), .{ .name = "__QNANL", .linkage = common.linkage, .visibility = common.visibility }); + @export(&math.snan(c_longdouble), .{ .name = "__SNANL", .linkage = common.linkage, .visibility = common.visibility }); + @export(&math.inf(c_longdouble), .{ .name = "__INFL", .linkage = common.linkage, .visibility = common.visibility }); + @export(&math.floatTrueMin(c_longdouble), .{ .name = "__DENORML", .linkage = common.linkage, .visibility = common.visibility }); } if (builtin.target.isMinGW() or builtin.target.isMuslLibC() or builtin.target.isWasiLibC()) { @@ -35,6 +36,9 @@ comptime { if (builtin.target.isMuslLibC() or builtin.target.isWasiLibC()) { @export(&acos, .{ .name = "acos", .linkage = common.linkage, .visibility = common.visibility }); + @export(&atanf, .{ .name = "atanf", .linkage = common.linkage, .visibility = common.visibility }); + @export(&atan, .{ .name = "atan", .linkage = common.linkage, .visibility = common.visibility }); + @export(&atanl, .{ .name = "atanl", .linkage = common.linkage, .visibility = common.visibility }); } if (builtin.target.isMuslLibC()) { @@ -45,41 +49,60 @@ comptime { } fn acos(x: f64) callconv(.c) f64 { - return std.math.acos(x); + return math.acos(x); +} + +fn atanf(x: f32) callconv(.c) f32 { + return math.atan(x); +} + +fn atan(x: f64) callconv(.c) f64 { + return math.atan(x); +} + +fn atanl(x: c_longdouble) callconv(.c) c_longdouble { + return switch (@typeInfo(@TypeOf(x)).float.bits) { + 16 => math.atan(@as(f16, @floatCast(x))), + 32 => math.atan(@as(f32, @floatCast(x))), + 64 => math.atan(@as(f64, @floatCast(x))), + 80 => math.atan(@as(f80, @floatCast(x))), + 128 => math.atan(@as(f128, @floatCast(x))), + else => unreachable, + }; } fn isnan(x: f64) callconv(.c) c_int { - return if (std.math.isNan(x)) 1 else 0; + return if (math.isNan(x)) 1 else 0; } fn isnanf(x: f32) callconv(.c) c_int { - return if (std.math.isNan(x)) 1 else 0; + return if (math.isNan(x)) 1 else 0; } fn isnanl(x: c_longdouble) callconv(.c) c_int { - return if (std.math.isNan(x)) 1 else 0; + return if (math.isNan(x)) 1 else 0; } fn nan(_: [*:0]const c_char) callconv(.c) f64 { - return std.math.nan(f64); + return math.nan(f64); } fn nanf(_: [*:0]const c_char) callconv(.c) f32 { - return std.math.nan(f32); + return math.nan(f32); } fn nanl(_: [*:0]const c_char) callconv(.c) c_longdouble { - return std.math.nan(c_longdouble); + return math.nan(c_longdouble); } fn copysignf(x: f32, y: f32) callconv(.c) f32 { - return std.math.copysign(x, y); + return math.copysign(x, y); } fn copysign(x: f64, y: f64) callconv(.c) f64 { - return std.math.copysign(x, y); + return math.copysign(x, y); } fn copysignl(x: c_longdouble, y: c_longdouble) callconv(.c) c_longdouble { - return std.math.copysign(x, y); + return math.copysign(x, y); } diff --git a/lib/libc/musl/src/math/atan.c b/lib/libc/musl/src/math/atan.c deleted file mode 100644 index 63b0ab25e3cf02ea81bab5a9ee4d99d6c40bb582..0000000000000000000000000000000000000000 --- a/lib/libc/musl/src/math/atan.c +++ /dev/null @@ -1,116 +0,0 @@ -/* origin: FreeBSD /usr/src/lib/msun/src/s_atan.c */ -/* - * ==================================================== - * 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. - * ==================================================== - */ -/* atan(x) - * Method - * 1. Reduce x to positive by atan(x) = -atan(-x). - * 2. According to the integer k=4t+0.25 chopped, t=x, the argument - * is further reduced to one of the following intervals and the - * arctangent of t is evaluated by the corresponding formula: - * - * [0,7/16] atan(x) = t-t^3*(a1+t^2*(a2+...(a10+t^2*a11)...) - * [7/16,11/16] atan(x) = atan(1/2) + atan( (t-0.5)/(1+t/2) ) - * [11/16.19/16] atan(x) = atan( 1 ) + atan( (t-1)/(1+t) ) - * [19/16,39/16] atan(x) = atan(3/2) + atan( (t-1.5)/(1+1.5t) ) - * [39/16,INF] atan(x) = atan(INF) + atan( -1/t ) - * - * Constants: - * The hexadecimal values are the intended ones for the following - * constants. The decimal values may be used, provided that the - * compiler will convert from decimal to binary accurately enough - * to produce the hexadecimal values shown. - */ - - -#include "libm.h" - -static const double atanhi[] = { - 4.63647609000806093515e-01, /* atan(0.5)hi 0x3FDDAC67, 0x0561BB4F */ - 7.85398163397448278999e-01, /* atan(1.0)hi 0x3FE921FB, 0x54442D18 */ - 9.82793723247329054082e-01, /* atan(1.5)hi 0x3FEF730B, 0xD281F69B */ - 1.57079632679489655800e+00, /* atan(inf)hi 0x3FF921FB, 0x54442D18 */ -}; - -static const double atanlo[] = { - 2.26987774529616870924e-17, /* atan(0.5)lo 0x3C7A2B7F, 0x222F65E2 */ - 3.06161699786838301793e-17, /* atan(1.0)lo 0x3C81A626, 0x33145C07 */ - 1.39033110312309984516e-17, /* atan(1.5)lo 0x3C700788, 0x7AF0CBBD */ - 6.12323399573676603587e-17, /* atan(inf)lo 0x3C91A626, 0x33145C07 */ -}; - -static const double aT[] = { - 3.33333333333329318027e-01, /* 0x3FD55555, 0x5555550D */ - -1.99999999998764832476e-01, /* 0xBFC99999, 0x9998EBC4 */ - 1.42857142725034663711e-01, /* 0x3FC24924, 0x920083FF */ - -1.11111104054623557880e-01, /* 0xBFBC71C6, 0xFE231671 */ - 9.09088713343650656196e-02, /* 0x3FB745CD, 0xC54C206E */ - -7.69187620504482999495e-02, /* 0xBFB3B0F2, 0xAF749A6D */ - 6.66107313738753120669e-02, /* 0x3FB10D66, 0xA0D03D51 */ - -5.83357013379057348645e-02, /* 0xBFADDE2D, 0x52DEFD9A */ - 4.97687799461593236017e-02, /* 0x3FA97B4B, 0x24760DEB */ - -3.65315727442169155270e-02, /* 0xBFA2B444, 0x2C6A6C2F */ - 1.62858201153657823623e-02, /* 0x3F90AD3A, 0xE322DA11 */ -}; - -double atan(double x) -{ - double_t w,s1,s2,z; - uint32_t ix,sign; - int id; - - GET_HIGH_WORD(ix, x); - sign = ix >> 31; - ix &= 0x7fffffff; - if (ix >= 0x44100000) { /* if |x| >= 2^66 */ - if (isnan(x)) - return x; - z = atanhi[3] + 0x1p-120f; - return sign ? -z : z; - } - if (ix < 0x3fdc0000) { /* |x| < 0.4375 */ - if (ix < 0x3e400000) { /* |x| < 2^-27 */ - if (ix < 0x00100000) - /* raise underflow for subnormal x */ - FORCE_EVAL((float)x); - return x; - } - id = -1; - } else { - x = fabs(x); - if (ix < 0x3ff30000) { /* |x| < 1.1875 */ - if (ix < 0x3fe60000) { /* 7/16 <= |x| < 11/16 */ - id = 0; - x = (2.0*x-1.0)/(2.0+x); - } else { /* 11/16 <= |x| < 19/16 */ - id = 1; - x = (x-1.0)/(x+1.0); - } - } else { - if (ix < 0x40038000) { /* |x| < 2.4375 */ - id = 2; - x = (x-1.5)/(1.0+1.5*x); - } else { /* 2.4375 <= |x| < 2^66 */ - id = 3; - x = -1.0/x; - } - } - } - /* end of argument reduction */ - z = x*x; - w = z*z; - /* break sum from i=0 to 10 aT[i]z**(i+1) into odd and even poly */ - s1 = z*(aT[0]+w*(aT[2]+w*(aT[4]+w*(aT[6]+w*(aT[8]+w*aT[10]))))); - s2 = w*(aT[1]+w*(aT[3]+w*(aT[5]+w*(aT[7]+w*aT[9])))); - if (id < 0) - return x - x*(s1+s2); - z = atanhi[id] - (x*(s1+s2) - atanlo[id] - x); - return sign ? -z : z; -} diff --git a/lib/libc/musl/src/math/atanf.c b/lib/libc/musl/src/math/atanf.c deleted file mode 100644 index 178341b670fa249fa50157d878ac2a66bd7f1843..0000000000000000000000000000000000000000 --- a/lib/libc/musl/src/math/atanf.c +++ /dev/null @@ -1,94 +0,0 @@ -/* origin: FreeBSD /usr/src/lib/msun/src/s_atanf.c */ -/* - * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. - */ -/* - * ==================================================== - * 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. - * ==================================================== - */ - - -#include "libm.h" - -static const float atanhi[] = { - 4.6364760399e-01, /* atan(0.5)hi 0x3eed6338 */ - 7.8539812565e-01, /* atan(1.0)hi 0x3f490fda */ - 9.8279368877e-01, /* atan(1.5)hi 0x3f7b985e */ - 1.5707962513e+00, /* atan(inf)hi 0x3fc90fda */ -}; - -static const float atanlo[] = { - 5.0121582440e-09, /* atan(0.5)lo 0x31ac3769 */ - 3.7748947079e-08, /* atan(1.0)lo 0x33222168 */ - 3.4473217170e-08, /* atan(1.5)lo 0x33140fb4 */ - 7.5497894159e-08, /* atan(inf)lo 0x33a22168 */ -}; - -static const float aT[] = { - 3.3333328366e-01, - -1.9999158382e-01, - 1.4253635705e-01, - -1.0648017377e-01, - 6.1687607318e-02, -}; - -float atanf(float x) -{ - float_t w,s1,s2,z; - uint32_t ix,sign; - int id; - - GET_FLOAT_WORD(ix, x); - sign = ix>>31; - ix &= 0x7fffffff; - if (ix >= 0x4c800000) { /* if |x| >= 2**26 */ - if (isnan(x)) - return x; - z = atanhi[3] + 0x1p-120f; - return sign ? -z : z; - } - if (ix < 0x3ee00000) { /* |x| < 0.4375 */ - if (ix < 0x39800000) { /* |x| < 2**-12 */ - if (ix < 0x00800000) - /* raise underflow for subnormal x */ - FORCE_EVAL(x*x); - return x; - } - id = -1; - } else { - x = fabsf(x); - if (ix < 0x3f980000) { /* |x| < 1.1875 */ - if (ix < 0x3f300000) { /* 7/16 <= |x| < 11/16 */ - id = 0; - x = (2.0f*x - 1.0f)/(2.0f + x); - } else { /* 11/16 <= |x| < 19/16 */ - id = 1; - x = (x - 1.0f)/(x + 1.0f); - } - } else { - if (ix < 0x401c0000) { /* |x| < 2.4375 */ - id = 2; - x = (x - 1.5f)/(1.0f + 1.5f*x); - } else { /* 2.4375 <= |x| < 2**26 */ - id = 3; - x = -1.0f/x; - } - } - } - /* end of argument reduction */ - z = x*x; - w = z*z; - /* break sum from i=0 to 10 aT[i]z**(i+1) into odd and even poly */ - s1 = z*(aT[0]+w*(aT[2]+w*aT[4])); - s2 = w*(aT[1]+w*aT[3]); - if (id < 0) - return x - x*(s1+s2); - z = atanhi[id] - ((x*(s1+s2) - atanlo[id]) - x); - return sign ? -z : z; -} diff --git a/lib/libc/musl/src/math/atanl.c b/lib/libc/musl/src/math/atanl.c deleted file mode 100644 index c3b0c9268db309060614a94e11d3755e2d26110a..0000000000000000000000000000000000000000 --- a/lib/libc/musl/src/math/atanl.c +++ /dev/null @@ -1,184 +0,0 @@ -/* origin: FreeBSD /usr/src/lib/msun/src/s_atanl.c */ -/* - * ==================================================== - * 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. - * ==================================================== - */ -/* - * See comments in atan.c. - * Converted to long double by David Schultz . - */ - -#include "libm.h" - -#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 -long double atanl(long double x) -{ - return atan(x); -} -#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 - -#if LDBL_MANT_DIG == 64 -#define EXPMAN(u) ((u.i.se & 0x7fff)<<8 | (u.i.m>>55 & 0xff)) - -static const long double atanhi[] = { - 4.63647609000806116202e-01L, - 7.85398163397448309628e-01L, - 9.82793723247329067960e-01L, - 1.57079632679489661926e+00L, -}; - -static const long double atanlo[] = { - 1.18469937025062860669e-20L, - -1.25413940316708300586e-20L, - 2.55232234165405176172e-20L, - -2.50827880633416601173e-20L, -}; - -static const long double aT[] = { - 3.33333333333333333017e-01L, - -1.99999999999999632011e-01L, - 1.42857142857046531280e-01L, - -1.11111111100562372733e-01L, - 9.09090902935647302252e-02L, - -7.69230552476207730353e-02L, - 6.66661718042406260546e-02L, - -5.88158892835030888692e-02L, - 5.25499891539726639379e-02L, - -4.70119845393155721494e-02L, - 4.03539201366454414072e-02L, - -2.91303858419364158725e-02L, - 1.24822046299269234080e-02L, -}; - -static long double T_even(long double x) -{ - return aT[0] + x * (aT[2] + x * (aT[4] + x * (aT[6] + - x * (aT[8] + x * (aT[10] + x * aT[12]))))); -} - -static long double T_odd(long double x) -{ - return aT[1] + x * (aT[3] + x * (aT[5] + x * (aT[7] + - x * (aT[9] + x * aT[11])))); -} -#elif LDBL_MANT_DIG == 113 -#define EXPMAN(u) ((u.i.se & 0x7fff)<<8 | u.i.top>>8) - -static const long double atanhi[] = { - 4.63647609000806116214256231461214397e-01L, - 7.85398163397448309615660845819875699e-01L, - 9.82793723247329067985710611014666038e-01L, - 1.57079632679489661923132169163975140e+00L, -}; - -static const long double atanlo[] = { - 4.89509642257333492668618435220297706e-36L, - 2.16795253253094525619926100651083806e-35L, - -2.31288434538183565909319952098066272e-35L, - 4.33590506506189051239852201302167613e-35L, -}; - -static const long double aT[] = { - 3.33333333333333333333333333333333125e-01L, - -1.99999999999999999999999999999180430e-01L, - 1.42857142857142857142857142125269827e-01L, - -1.11111111111111111111110834490810169e-01L, - 9.09090909090909090908522355708623681e-02L, - -7.69230769230769230696553844935357021e-02L, - 6.66666666666666660390096773046256096e-02L, - -5.88235294117646671706582985209643694e-02L, - 5.26315789473666478515847092020327506e-02L, - -4.76190476189855517021024424991436144e-02L, - 4.34782608678695085948531993458097026e-02L, - -3.99999999632663469330634215991142368e-02L, - 3.70370363987423702891250829918659723e-02L, - -3.44827496515048090726669907612335954e-02L, - 3.22579620681420149871973710852268528e-02L, - -3.03020767654269261041647570626778067e-02L, - 2.85641979882534783223403715930946138e-02L, - -2.69824879726738568189929461383741323e-02L, - 2.54194698498808542954187110873675769e-02L, - -2.35083879708189059926183138130183215e-02L, - 2.04832358998165364349957325067131428e-02L, - -1.54489555488544397858507248612362957e-02L, - 8.64492360989278761493037861575248038e-03L, - -2.58521121597609872727919154569765469e-03L, -}; - -static long double T_even(long double x) -{ - return (aT[0] + x * (aT[2] + x * (aT[4] + x * (aT[6] + x * (aT[8] + - x * (aT[10] + x * (aT[12] + x * (aT[14] + x * (aT[16] + - x * (aT[18] + x * (aT[20] + x * aT[22]))))))))))); -} - -static long double T_odd(long double x) -{ - return (aT[1] + x * (aT[3] + x * (aT[5] + x * (aT[7] + x * (aT[9] + - x * (aT[11] + x * (aT[13] + x * (aT[15] + x * (aT[17] + - x * (aT[19] + x * (aT[21] + x * aT[23]))))))))))); -} -#endif - -long double atanl(long double x) -{ - union ldshape u = {x}; - long double w, s1, s2, z; - int id; - unsigned e = u.i.se & 0x7fff; - unsigned sign = u.i.se >> 15; - unsigned expman; - - if (e >= 0x3fff + LDBL_MANT_DIG + 1) { /* if |x| is large, atan(x)~=pi/2 */ - if (isnan(x)) - return x; - return sign ? -atanhi[3] : atanhi[3]; - } - /* Extract the exponent and the first few bits of the mantissa. */ - expman = EXPMAN(u); - if (expman < ((0x3fff - 2) << 8) + 0xc0) { /* |x| < 0.4375 */ - if (e < 0x3fff - (LDBL_MANT_DIG+1)/2) { /* if |x| is small, atanl(x)~=x */ - /* raise underflow if subnormal */ - if (e == 0) - FORCE_EVAL((float)x); - return x; - } - id = -1; - } else { - x = fabsl(x); - if (expman < (0x3fff << 8) + 0x30) { /* |x| < 1.1875 */ - if (expman < ((0x3fff - 1) << 8) + 0x60) { /* 7/16 <= |x| < 11/16 */ - id = 0; - x = (2.0*x-1.0)/(2.0+x); - } else { /* 11/16 <= |x| < 19/16 */ - id = 1; - x = (x-1.0)/(x+1.0); - } - } else { - if (expman < ((0x3fff + 1) << 8) + 0x38) { /* |x| < 2.4375 */ - id = 2; - x = (x-1.5)/(1.0+1.5*x); - } else { /* 2.4375 <= |x| */ - id = 3; - x = -1.0/x; - } - } - } - /* end of argument reduction */ - z = x*x; - w = z*z; - /* break sum aT[i]z**(i+1) into odd and even poly */ - s1 = z*T_even(w); - s2 = w*T_odd(w); - if (id < 0) - return x - x*(s1+s2); - z = atanhi[id] - ((x*(s1+s2) - atanlo[id]) - x); - return sign ? -z : z; -} -#endif diff --git a/lib/libc/musl/src/math/i386/atan.s b/lib/libc/musl/src/math/i386/atan.s deleted file mode 100644 index 2c57f6b309d3d5e54430f17db69171e0ce5bc833..0000000000000000000000000000000000000000 --- a/lib/libc/musl/src/math/i386/atan.s +++ /dev/null @@ -1,16 +0,0 @@ -.global atan -.type atan,@function -atan: - fldl 4(%esp) - mov 8(%esp),%eax - add %eax,%eax - cmp $0x00200000,%eax - jb 1f - fld1 - fpatan - fstpl 4(%esp) - fldl 4(%esp) - ret - # subnormal x, return x with underflow -1: fsts 4(%esp) - ret diff --git a/lib/libc/musl/src/math/i386/atanf.s b/lib/libc/musl/src/math/i386/atanf.s deleted file mode 100644 index c2cbe2e0267f93272a195725489afcbbc24c8876..0000000000000000000000000000000000000000 --- a/lib/libc/musl/src/math/i386/atanf.s +++ /dev/null @@ -1,18 +0,0 @@ -.global atanf -.type atanf,@function -atanf: - flds 4(%esp) - mov 4(%esp),%eax - add %eax,%eax - cmp $0x01000000,%eax - jb 1f - fld1 - fpatan - fstps 4(%esp) - flds 4(%esp) - ret - # subnormal x, return x with underflow -1: fld %st(0) - fmul %st(1) - fstps 4(%esp) - ret diff --git a/lib/libc/musl/src/math/i386/atanl.s b/lib/libc/musl/src/math/i386/atanl.s deleted file mode 100644 index c508bc465b1a510efca440808bc88fe6c2a7d426..0000000000000000000000000000000000000000 --- a/lib/libc/musl/src/math/i386/atanl.s +++ /dev/null @@ -1,7 +0,0 @@ -.global atanl -.type atanl,@function -atanl: - fldt 4(%esp) - fld1 - fpatan - ret diff --git a/lib/libc/musl/src/math/x32/atanl.s b/lib/libc/musl/src/math/x32/atanl.s deleted file mode 100644 index f475fe0e9ee73ca2940fd47c9326323769207953..0000000000000000000000000000000000000000 --- a/lib/libc/musl/src/math/x32/atanl.s +++ /dev/null @@ -1,7 +0,0 @@ -.global atanl -.type atanl,@function -atanl: - fldt 8(%esp) - fld1 - fpatan - ret diff --git a/lib/libc/musl/src/math/x86_64/atanl.s b/lib/libc/musl/src/math/x86_64/atanl.s deleted file mode 100644 index df76de5de4f12834a4a3710a2e7370d10245e862..0000000000000000000000000000000000000000 --- a/lib/libc/musl/src/math/x86_64/atanl.s +++ /dev/null @@ -1,7 +0,0 @@ -.global atanl -.type atanl,@function -atanl: - fldt 8(%rsp) - fld1 - fpatan - ret diff --git a/lib/std/math/atan.zig b/lib/std/math/atan.zig index 377d96897ff7acad0e76aece1866fdbc4fd31d4e..8b783bf3474b7cf260c9ca4930e35200dd6704bd 100644 --- a/lib/std/math/atan.zig +++ b/lib/std/math/atan.zig @@ -3,11 +3,12 @@ // // https://git.musl-libc.org/cgit/musl/tree/src/math/atanf.c // https://git.musl-libc.org/cgit/musl/tree/src/math/atan.c +// https://git.musl-libc.org/cgit/musl/tree/src/math/atanl.c const std = @import("../std.zig"); const math = std.math; const mem = std.mem; -const expect = std.testing.expect; +const testing = std.testing; /// Returns the arc-tangent of x. /// @@ -17,28 +18,100 @@ const expect = std.testing.expect; pub fn atan(x: anytype) @TypeOf(x) { const T = @TypeOf(x); return switch (T) { - f32 => atan32(x), - f64 => atan64(x), + f16 => atanBinary16(x), + f32 => atanBinary32(x), + f64 => atanBinary64(x), + f80 => atanExtended80(x), + f128 => atanBinary128(x), else => @compileError("atan not implemented for " ++ @typeName(T)), }; } -fn atan32(x_: f32) f32 { - const atanhi = [_]f32{ - 4.6364760399e-01, // atan(0.5)hi - 7.8539812565e-01, // atan(1.0)hi - 9.8279368877e-01, // atan(1.5)hi - 1.5707962513e+00, // atan(inf)hi +fn atanBinary16(x: f16) f16 { + const atanhi: []const f32 = &.{ + 4.6364760399e-01, // atan(0.5)hi 0x3eed6338 + 7.8539812565e-01, // atan(1.0)hi 0x3f490fda + 9.8279368877e-01, // atan(1.5)hi 0x3f7b985e + 1.5707962513e+00, // atan(inf)hi 0x3fc90fda + }; + const aT: []const f32 = &.{ + 0x1.fffcccp-1, + -0x1.52e8ccp-2, + 0x1.522336p-3, }; - const atanlo = [_]f32{ - 5.0121582440e-09, // atan(0.5)lo - 3.7748947079e-08, // atan(1.0)lo - 3.4473217170e-08, // atan(1.5)lo - 7.5497894159e-08, // atan(inf)lo + const hx: u16 = @bitCast(x); + const ix = hx & 0x7fff; + const sign = (hx >> 15) != 0; + // if |x| >= 2^11 + if (ix >= 0x6800) { + if (math.isNan(x)) { + return x; + } + const z = atanhi[3] + 0x1p-120; + return @floatCast(if (sign) -z else z); + } + const x_: f32, const id: ?usize = blk: { + // |x| < 0.4375 + if (ix < 0x3700) { + // |x| < 2^(-6) + if (ix < 0x2400) { + if (ix < 0x400) { + // raise underflow for subnormal x + mem.doNotOptimizeAway(x * x); + } + return x; + } + break :blk .{ @floatCast(x), null }; + } else { + const x_: f32 = @floatCast(@abs(x)); + // |x| < 1.1875 + if (ix < 0x3cc0) { + // 7/16 <= |x| < 11/16 + if (ix < 0x3980) { + break :blk .{ (2.0 * x_ - 1.0) / (2.0 + x_), 0 }; + } + // 11/16 <= |x| < 19/16 + else { + break :blk .{ (x_ - 1.0) / (x_ + 1.0), 1 }; + } + } else { + // |x| < 2.4375 + if (ix < 0x40e0) { + break :blk .{ (x_ - 1.5) / (1.0 + 1.5 * x_), 2 }; + } + // 2.4375 <= |x| < 2^11 + else { + break :blk .{ -1.0 / x_, 3 }; + } + } + } }; + // end of argument reduction + const z = x_ * x_; + const s = aT[0] + z * (aT[1] + z * aT[2]); + if (id) |id_| { + const z_ = atanhi[id_] + x_ * s; + return @floatCast(if (sign) -z_ else z_); + } else { + return @floatCast(x_ * s); + } +} - const aT = [_]f32{ +fn atanBinary32(x: f32) f32 { + const atanhi: []const f32 = &.{ + 4.6364760399e-01, // atan(0.5)hi 0x3eed6338 + 7.8539812565e-01, // atan(1.0)hi 0x3f490fda + 9.8279368877e-01, // atan(1.5)hi 0x3f7b985e + 1.5707962513e+00, // atan(inf)hi 0x3fc90fda + }; + const atanlo: []const f32 = &.{ + 5.0121582440e-09, // atan(0.5)lo 0x31ac3769 + 3.7748947079e-08, // atan(1.0)lo 0x33222168 + 3.4473217170e-08, // atan(1.5)lo 0x33140fb4 + 7.5497894159e-08, // atan(inf)lo 0x33a22168 + }; + const aT: []const f32 = &.{ 3.3333328366e-01, -1.9999158382e-01, 1.4253635705e-01, @@ -46,211 +119,463 @@ fn atan32(x_: f32) f32 { 6.1687607318e-02, }; - var x = x_; - var ix: u32 = @as(u32, @bitCast(x)); - const sign = ix >> 31; - ix &= 0x7FFFFFFF; - - // |x| >= 2^26 - if (ix >= 0x4C800000) { + const hx: u32 = @bitCast(x); + const ix = hx & 0x7fff_ffff; + const sign = (hx >> 31) != 0; + // if |x| >= 2^26 + if (ix >= 0x4c80_0000) { if (math.isNan(x)) { return x; - } else { - const z = atanhi[3] + 0x1.0p-120; - return if (sign != 0) -z else z; } + const z = atanhi[3] + 0x1p-120; + return if (sign) -z else z; } - - var id: ?usize = undefined; - - // |x| < 0.4375 - if (ix < 0x3EE00000) { - // |x| < 2^(-12) - if (ix < 0x39800000) { - if (ix < 0x00800000) { - mem.doNotOptimizeAway(x * x); - } - return x; - } - id = null; - } else { - x = @abs(x); - // |x| < 1.1875 - if (ix < 0x3F980000) { - // 7/16 <= |x| < 11/16 - if (ix < 0x3F300000) { - id = 0; - x = (2.0 * x - 1.0) / (2.0 + x); - } - // 11/16 <= |x| < 19/16 - else { - id = 1; - x = (x - 1.0) / (x + 1.0); + const x_, const id: ?usize = blk: { + // |x| < 0.4375 + if (ix < 0x3ee00000) { + // |x| < 2^(-12) + if (ix < 0x39800000) { + if (ix < 0x00800000) { + // raise underflow for subnormal x + mem.doNotOptimizeAway(x * x); + } + return x; } + break :blk .{ x, null }; } else { - // |x| < 2.4375 - if (ix < 0x401C0000) { - id = 2; - x = (x - 1.5) / (1.0 + 1.5 * x); - } - // 2.4375 <= |x| < 2^26 - else { - id = 3; - x = -1.0 / x; + const x_ = @abs(x); + // |x| < 1.1875 + if (ix < 0x3f98_0000) { + // 7/16 <= |x| < 11/16 + if (ix < 0x3f30_0000) { + break :blk .{ (2.0 * x_ - 1.0) / (2.0 + x_), 0 }; + } + // 11/16 <= |x| < 19/16 + else { + break :blk .{ (x_ - 1.0) / (x_ + 1.0), 1 }; + } + } else { + // |x| < 2.4375 + if (ix < 0x401c_0000) { + break :blk .{ (x_ - 1.5) / (1.0 + 1.5 * x_), 2 }; + } + // 2.4375 <= |x| < 2^26 + else { + break :blk .{ -1.0 / x_, 3 }; + } } } - } - - const z = x * x; + }; + // end of argument reduction + const z = x_ * x_; const w = z * z; + // break sum from i=0 to 10 aT[i]z^(i+1) into odd and even poly const s1 = z * (aT[0] + w * (aT[2] + w * aT[4])); const s2 = w * (aT[1] + w * aT[3]); - - if (id) |id_value| { - const zz = atanhi[id_value] - ((x * (s1 + s2) - atanlo[id_value]) - x); - return if (sign != 0) -zz else zz; + if (id) |id_| { + const z_ = atanhi[id_] - ((x_ * (s1 + s2) - atanlo[id_]) - x_); + return if (sign) -z_ else z_; } else { - return x - x * (s1 + s2); + return x_ - x_ * (s1 + s2); } } -fn atan64(x_: f64) f64 { - const atanhi = [_]f64{ - 4.63647609000806093515e-01, // atan(0.5)hi - 7.85398163397448278999e-01, // atan(1.0)hi - 9.82793723247329054082e-01, // atan(1.5)hi - 1.57079632679489655800e+00, // atan(inf)hi +fn atanBinary64(x: f64) f64 { + const atanhi: []const f64 = &.{ + 4.63647609000806093515e-01, // atan(0.5)hi 0x3FDDAC67, 0x0561BB4F + 7.85398163397448278999e-01, // atan(1.0)hi 0x3FE921FB, 0x54442D18 + 9.82793723247329054082e-01, // atan(1.5)hi 0x3FEF730B, 0xD281F69B + 1.57079632679489655800e+00, // atan(inf)hi 0x3FF921FB, 0x54442D18 }; - - const atanlo = [_]f64{ - 2.26987774529616870924e-17, // atan(0.5)lo - 3.06161699786838301793e-17, // atan(1.0)lo - 1.39033110312309984516e-17, // atan(1.5)lo - 6.12323399573676603587e-17, // atan(inf)lo + const atanlo: []const f64 = &.{ + 2.26987774529616870924e-17, // atan(0.5)lo 0x3C7A2B7F, 0x222F65E2 + 3.06161699786838301793e-17, // atan(1.0)lo 0x3C81A626, 0x33145C07 + 1.39033110312309984516e-17, // atan(1.5)lo 0x3C700788, 0x7AF0CBBD + 6.12323399573676603587e-17, // atan(inf)lo 0x3C91A626, 0x33145C07 }; - - const aT = [_]f64{ - 3.33333333333329318027e-01, - -1.99999999998764832476e-01, - 1.42857142725034663711e-01, - -1.11111104054623557880e-01, - 9.09088713343650656196e-02, - -7.69187620504482999495e-02, - 6.66107313738753120669e-02, - -5.83357013379057348645e-02, - 4.97687799461593236017e-02, - -3.65315727442169155270e-02, - 1.62858201153657823623e-02, + const aT: []const f64 = &.{ + 3.33333333333329318027e-01, // 0x3FD55555, 0x5555550D + -1.99999999998764832476e-01, // 0xBFC99999, 0x9998EBC4 + 1.42857142725034663711e-01, // 0x3FC24924, 0x920083FF + -1.11111104054623557880e-01, // 0xBFBC71C6, 0xFE231671 + 9.09088713343650656196e-02, // 0x3FB745CD, 0xC54C206E + -7.69187620504482999495e-02, // 0xBFB3B0F2, 0xAF749A6D + 6.66107313738753120669e-02, // 0x3FB10D66, 0xA0D03D51 + -5.83357013379057348645e-02, // 0xBFADDE2D, 0x52DEFD9A + 4.97687799461593236017e-02, // 0x3FA97B4B, 0x24760DEB + -3.65315727442169155270e-02, // 0xBFA2B444, 0x2C6A6C2F + 1.62858201153657823623e-02, // 0x3F90AD3A, 0xE322DA11 }; - var x = x_; - const ux: u64 = @bitCast(x); - var ix: u32 = @intCast(ux >> 32); - const sign = ix >> 31; - ix &= 0x7FFFFFFF; - - // |x| >= 2^66 + const hx: u64 = @bitCast(x); + const ix: u32 = @truncate((hx >> 32) & 0x7fffffff); + const sign = (hx >> 63) != 0; + // if |x| >= 2^66 if (ix >= 0x44100000) { if (math.isNan(x)) { return x; - } else { - const z = atanhi[3] + 0x1.0p-120; - return if (sign != 0) -z else z; } + const z = atanhi[3] + 0x1p-120; + return if (sign) -z else z; } - - var id: ?usize = undefined; - - // |x| < 0.4375 - if (ix < 0x3FDC0000) { - // |x| < 2^(-27) - if (ix < 0x3E400000) { - if (ix < 0x00100000) { - mem.doNotOptimizeAway(@as(f32, @floatCast(x))); - } - return x; - } - id = null; - } else { - x = @abs(x); - // |x| < 1.1875 - if (ix < 0x3FF30000) { - // 7/16 <= |x| < 11/16 - if (ix < 0x3FE60000) { - id = 0; - x = (2.0 * x - 1.0) / (2.0 + x); - } - // 11/16 <= |x| < 19/16 - else { - id = 1; - x = (x - 1.0) / (x + 1.0); + const x_, const id: ?usize = blk: { + // |x| < 0.4375 + if (ix < 0x3fdc_0000) { + // |x| < 2^(-27) + if (ix < 0x3e40_0000) { + if (ix < 0x0010_0000) { + // raise underflow for subnormal x + mem.doNotOptimizeAway(@as(f32, @floatCast(x))); + } + return x; } + break :blk .{ x, null }; } else { - // |x| < 2.4375 - if (ix < 0x40038000) { - id = 2; - x = (x - 1.5) / (1.0 + 1.5 * x); - } - // 2.4375 <= |x| < 2^66 - else { - id = 3; - x = -1.0 / x; + const x_ = @abs(x); + // |x| < 1.1875 + if (ix < 0x3ff3_0000) { + // 7/16 <= |x| < 11/16 + if (ix < 0x3fe6_0000) { + break :blk .{ (2.0 * x_ - 1.0) / (2.0 + x_), 0 }; + } + // 11/16 <= |x| < 19/16 + else { + break :blk .{ (x_ - 1.0) / (x_ + 1.0), 1 }; + } + } else { + // |x| < 2.4375 + if (ix < 0x4003_8000) { + break :blk .{ (x_ - 1.5) / (1.0 + 1.5 * x_), 2 }; + } + // 2.4375 <= |x| < 2^66 + else { + break :blk .{ -1.0 / x_, 3 }; + } } } - } - - const z = x * x; + }; + // end of argument reduction + const z = x_ * x_; const w = z * z; + // break sum from i=0 to 10 aT[i]z^(i+1) into odd and even poly const s1 = z * (aT[0] + w * (aT[2] + w * (aT[4] + w * (aT[6] + w * (aT[8] + w * aT[10]))))); const s2 = w * (aT[1] + w * (aT[3] + w * (aT[5] + w * (aT[7] + w * aT[9])))); + if (id) |id_| { + const z_ = atanhi[id_] - (x_ * (s1 + s2) - atanlo[id_] - x_); + return if (sign) -z_ else z_; + } else { + return x_ - x_ * (s1 + s2); + } +} - if (id) |id_value| { - const zz = atanhi[id_value] - ((x * (s1 + s2) - atanlo[id_value]) - x); - return if (sign != 0) -zz else zz; +fn atanExtended80(x: f80) f80 { + const atanhi: []const f80 = &.{ + 4.63647609000806116202e-01, + 7.85398163397448309628e-01, + 9.82793723247329067960e-01, + 1.57079632679489661926e+00, + }; + const atanlo: []const f80 = &.{ + 1.18469937025062860669e-20, + -1.25413940316708300586e-20, + 2.55232234165405176172e-20, + -2.50827880633416601173e-20, + }; + const aT: []const f80 = &.{ + 3.33333333333333333017e-01, + -1.99999999999999632011e-01, + 1.42857142857046531280e-01, + -1.11111111100562372733e-01, + 9.09090902935647302252e-02, + -7.69230552476207730353e-02, + 6.66661718042406260546e-02, + -5.88158892835030888692e-02, + 5.25499891539726639379e-02, + -4.70119845393155721494e-02, + 4.03539201366454414072e-02, + -2.91303858419364158725e-02, + 1.24822046299269234080e-02, + }; + + const hx: u80 = @bitCast(x); + const se: u16 = @truncate(hx >> 64); + const e = se & 0x7fff; + const sign = se >> 15 != 0; + // if |x| is large, atan(x)~=pi/2 + if (e >= 0x3fff + math.floatMantissaBits(f80) + 1) { + if (math.isNan(x)) { + return x; + } + return if (sign) -atanhi[3] else atanhi[3]; + } + // Extract the exponent and the first few bits of the mantissa. + const m: u64 = @truncate(hx & 0x0000_ffff_ffff_ffff_ffff); + const expman = ((@as(u32, @intCast(se)) & 0x7fff) << 8) | (@as(u32, @truncate(m >> 55)) & 0xff); + const x_, const id: ?usize = blk: { + // |x| < 0.4375 + if (expman < ((0x3fff - 2) << 8) + 0xc0) { + // if |x| is small, atanl(x)~=x + if (e < 0x3fff - (math.floatMantissaBits(f80) + 1) / 2) { + // raise underflow if subnormal + if (e == 0) { + std.mem.doNotOptimizeAway(@as(f32, @floatCast(x))); + } + return x; + } + break :blk .{ x, null }; + } else { + const x_ = @abs(x); + // |x| < 1.1875 + if (expman < (0x3fff << 8) + 0x30) { + // 7/16 <= |x| < 11/16 + if (expman < ((0x3fff - 1) << 8) + 0x60) { + break :blk .{ (2.0 * x_ - 1.0) / (2.0 + x_), 0 }; + } + // 11/16 <= |x| < 19/16 + else { + break :blk .{ (x_ - 1.0) / (x_ + 1.0), 1 }; + } + } else { + // |x| < 2.4375 + if (expman < ((0x3fff + 1) << 8) + 0x38) { + break :blk .{ (x_ - 1.5) / (1.0 + 1.5 * x_), 2 }; + } + // 2.4375 <= |x| + else { + break :blk .{ -1.0 / x_, 3 }; + } + } + } + }; + // end of argument reduction + const z = x_ * x_; + const w = z * z; + // break sum aT[i]z^(i+1) into odd and even poly + const s1 = z * (aT[0] + w * (aT[2] + w * (aT[4] + w * (aT[6] + w * (aT[8] + w * (aT[10] + w * aT[12])))))); + const s2 = w * (aT[1] + w * (aT[3] + w * (aT[5] + w * (aT[7] + w * (aT[9] + w * aT[11]))))); + if (id) |id_| { + const z_ = atanhi[id_] - ((x_ * (s1 + s2) - atanlo[id_]) - x_); + return if (sign) -z_ else z_; } else { - return x - x * (s1 + s2); + return x_ - x_ * (s1 + s2); } } -test atan { - try expect(@as(u32, @bitCast(atan(@as(f32, 0.2)))) == @as(u32, @bitCast(atan32(0.2)))); - try expect(atan(@as(f64, 0.2)) == atan64(0.2)); +fn atanBinary128(x: f128) f128 { + const atanhi: []const f128 = &.{ + 4.63647609000806116214256231461214397e-01, + 7.85398163397448309615660845819875699e-01, + 9.82793723247329067985710611014666038e-01, + 1.57079632679489661923132169163975140e+00, + }; + const atanlo: []const f128 = &.{ + 4.89509642257333492668618435220297706e-36, + 2.16795253253094525619926100651083806e-35, + -2.31288434538183565909319952098066272e-35, + 4.33590506506189051239852201302167613e-35, + }; + const aT: []const f128 = &.{ + 3.33333333333333333333333333333333125e-01, + -1.99999999999999999999999999999180430e-01, + 1.42857142857142857142857142125269827e-01, + -1.11111111111111111111110834490810169e-01, + 9.09090909090909090908522355708623681e-02, + -7.69230769230769230696553844935357021e-02, + 6.66666666666666660390096773046256096e-02, + -5.88235294117646671706582985209643694e-02, + 5.26315789473666478515847092020327506e-02, + -4.76190476189855517021024424991436144e-02, + 4.34782608678695085948531993458097026e-02, + -3.99999999632663469330634215991142368e-02, + 3.70370363987423702891250829918659723e-02, + -3.44827496515048090726669907612335954e-02, + 3.22579620681420149871973710852268528e-02, + -3.03020767654269261041647570626778067e-02, + 2.85641979882534783223403715930946138e-02, + -2.69824879726738568189929461383741323e-02, + 2.54194698498808542954187110873675769e-02, + -2.35083879708189059926183138130183215e-02, + 2.04832358998165364349957325067131428e-02, + -1.54489555488544397858507248612362957e-02, + 8.64492360989278761493037861575248038e-03, + -2.58521121597609872727919154569765469e-03, + }; + + const hx: u128 = @bitCast(x); + const se: u16 = @truncate(hx >> 112); + const e = se & 0x7fff; + const sign = se >> 15 != 0; + // if |x| is large, atan(x)~=pi/2 + if (e >= 0x3fff + math.floatMantissaBits(f128) + 2) { + if (math.isNan(x)) { + return x; + } + return if (sign) -atanhi[3] else atanhi[3]; + } + // Extract the exponent and the first few bits of the mantissa. + const top: u16 = @truncate((hx >> 96) & 0x0000_ffff); + const expman = ((@as(u32, @intCast(se)) & 0x7fff) << 8) | (@as(u32, @intCast(top)) >> 8); + const x_, const id: ?usize = blk: { + // |x| < 0.4375 + if (expman < ((0x3fff - 2) << 8) + 0xc0) { + // if |x| is small, atanl(x)~=x + if (e < 0x3fff - (math.floatMantissaBits(f128) + 2) / 2) { + // raise underflow if subnormal + if (e == 0) { + mem.doNotOptimizeAway(@as(f32, @floatCast(x))); + } + return x; + } + break :blk .{ x, null }; + } else { + const x_ = @abs(x); + // |x| < 1.1875 + if (expman < (0x3fff << 8) + 0x30) { + // 7/16 <= |x| < 11/16 + if (expman < ((0x3fff - 1) << 8) + 0x60) { + break :blk .{ (2.0 * x_ - 1.0) / (2.0 + x_), 0 }; + } + // 11/16 <= |x| < 19/16 + else { + break :blk .{ (x_ - 1.0) / (x_ + 1.0), 1 }; + } + } else { + // |x| < 2.4375 + if (expman < ((0x3fff + 1) << 8) + 0x38) { + break :blk .{ (x_ - 1.5) / (1.0 + 1.5 * x_), 2 }; + } + // 2.4375 <= |x| + else { + break :blk .{ -1.0 / x_, 3 }; + } + } + } + }; + // end of argument reduction + const z = x_ * x_; + const w = z * z; + // break sum aT[i]z^(i+1) into odd and even poly + const s1 = z * (aT[0] + w * (aT[2] + w * (aT[4] + w * (aT[6] + w * (aT[8] + w * (aT[10] + w * (aT[12] + w * (aT[14] + w * (aT[16] + w * (aT[18] + w * (aT[20] + w * aT[22]))))))))))); + const s2 = w * (aT[1] + w * (aT[3] + w * (aT[5] + w * (aT[7] + w * (aT[9] + w * (aT[11] + w * (aT[13] + w * (aT[15] + w * (aT[17] + w * (aT[19] + w * (aT[21] + w * aT[23]))))))))))); + if (id) |id_| { + const z_ = atanhi[id_] - ((x_ * (s1 + s2) - atanlo[id_]) - x_); + return if (sign) -z_ else z_; + } else { + return x_ - x_ * (s1 + s2); + } } -test atan32 { - const epsilon = 0.000001; +test "atanBinary16.special" { + try testing.expectEqual(atanBinary16(0x0p+0), 0x0p+0); + try testing.expectEqual(atanBinary16(-0x0p+0), -0x0p+0); + try testing.expectApproxEqAbs(atanBinary16(0x1p+0), 0x1.92p-1, math.floatEpsAt(f16, 0x1.92p-1)); + try testing.expectApproxEqAbs(atanBinary16(-0x1p+0), -0x1.92p-1, math.floatEpsAt(f16, -0x1.92p-1)); + try testing.expectApproxEqAbs(atanBinary16(math.inf(f16)), 0x1.92p0, math.floatEpsAt(f16, 0x1.92p0)); + try testing.expectApproxEqAbs(atanBinary16(-math.inf(f16)), -0x1.92p0, math.floatEpsAt(f16, -0x1.92p0)); + try testing.expect(math.isNan(atanBinary16(math.nan(f16)))); +} + +test "atanBinary16" { + try testing.expectApproxEqAbs(atanBinary16(-0x1.864p-2), -0x1.74cp-2, math.floatEpsAt(f16, -0x1.74cp-2)); + try testing.expectApproxEqAbs(atanBinary16(-0x1.59cp1), -0x1.374p0, math.floatEpsAt(f16, -0x1.374p0)); + try testing.expectApproxEqAbs(atanBinary16(-0x1.d2cp0), -0x1.11cp0, math.floatEpsAt(f16, -0x1.11cp0)); + try testing.expectApproxEqAbs(atanBinary16(-0x1.5f4p-1), -0x1.33cp-1, math.floatEpsAt(f16, -0x1.33cp-1)); + try testing.expectApproxEqAbs(atanBinary16(0x1.588p1), 0x1.37p0, math.floatEpsAt(f16, 0x1.37p0)); + try testing.expectApproxEqAbs(atanBinary16(-0x1.b14p-2), -0x1.99cp-2, math.floatEpsAt(f16, -0x1.99cp-2)); + try testing.expectApproxEqAbs(atanBinary16(0x1.3ccp1), 0x1.2fcp0, math.floatEpsAt(f16, 0x1.2fcp0)); + try testing.expectApproxEqAbs(atanBinary16(-0x1.0ecp-2), -0x1.08cp-2, math.floatEpsAt(f16, -0x1.08cp-2)); + try testing.expectApproxEqAbs(atanBinary16(0x1.298p1), 0x1.2ap0, math.floatEpsAt(f16, 0x1.2ap0)); + try testing.expectApproxEqAbs(atanBinary16(-0x1.028p1), -0x1.1c8p0, math.floatEpsAt(f16, -0x1.1c8p0)); +} + +test "atanBinary32.special" { + try testing.expectEqual(atanBinary32(0x0p+0), 0x0p+0); + try testing.expectEqual(atanBinary32(-0x0p+0), -0x0p+0); + try testing.expectApproxEqAbs(atanBinary32(0x1p+0), 0x1.921fb6p-1, math.floatEpsAt(f32, 0x1.921fb6p-1)); + try testing.expectApproxEqAbs(atanBinary32(-0x1p+0), -0x1.921fb6p-1, math.floatEpsAt(f32, -0x1.921fb6p-1)); + try testing.expectApproxEqAbs(atanBinary32(math.inf(f32)), 0x1.921fb6p+0, math.floatEpsAt(f32, 0x1.921fb6p+0)); + try testing.expectApproxEqAbs(atanBinary32(-math.inf(f32)), -0x1.921fb6p+0, math.floatEpsAt(f32, -0x1.921fb6p+0)); + try testing.expect(math.isNan(atanBinary32(math.nan(f32)))); +} - try expect(math.approxEqAbs(f32, atan32(0.2), 0.197396, epsilon)); - try expect(math.approxEqAbs(f32, atan32(-0.2), -0.197396, epsilon)); - try expect(math.approxEqAbs(f32, atan32(0.3434), 0.330783, epsilon)); - try expect(math.approxEqAbs(f32, atan32(0.8923), 0.728545, epsilon)); - try expect(math.approxEqAbs(f32, atan32(1.5), 0.982794, epsilon)); +test "atanBinary32" { + try testing.expectApproxEqAbs(atanBinary32(-0x1.8629dp-2), -0x1.74c62p-2, math.floatEpsAt(f32, -0x1.74c62p-2)); + try testing.expectApproxEqAbs(atanBinary32(-0x1.59d42ep1), -0x1.375fd8p0, math.floatEpsAt(f32, -0x1.375fd8p0)); + try testing.expectApproxEqAbs(atanBinary32(-0x1.d2dbe2p0), -0x1.11b8aep0, math.floatEpsAt(f32, -0x1.11b8aep0)); + try testing.expectApproxEqAbs(atanBinary32(-0x1.5f314ep-1), -0x1.33d28cp-1, math.floatEpsAt(f32, -0x1.33d28cp-1)); + try testing.expectApproxEqAbs(atanBinary32(0x1.5869bp1), 0x1.37082ep0, math.floatEpsAt(f32, 0x1.37082ep0)); + try testing.expectApproxEqAbs(atanBinary32(-0x1.b13a06p-2), -0x1.99d7cap-2, math.floatEpsAt(f32, -0x1.99d7cap-2)); + try testing.expectApproxEqAbs(atanBinary32(0x1.3cb0f2p1), 0x1.2fcb12p0, math.floatEpsAt(f32, 0x1.2fcb12p0)); + try testing.expectApproxEqAbs(atanBinary32(-0x1.0ed746p-2), -0x1.08c71ap-2, math.floatEpsAt(f32, -0x1.08c71ap-2)); + try testing.expectApproxEqAbs(atanBinary32(0x1.299d54p1), 0x1.2a24e2p0, math.floatEpsAt(f32, 0x1.2a24e2p0)); + try testing.expectApproxEqAbs(atanBinary32(-0x1.0264fcp1), -0x1.1c6178p0, math.floatEpsAt(f32, -0x1.1c6178p0)); } -test atan64 { - const epsilon = 0.000001; +test "atanBinary64.special" { + try testing.expectEqual(atanBinary64(0x0p+0), 0x0p+0); + try testing.expectEqual(atanBinary64(-0x0p+0), -0x0p+0); + try testing.expectApproxEqAbs(atanBinary64(0x1p+0), 0x1.921fb54442d18p-1, math.floatEpsAt(f64, 0x1.921fb54442d18p-1)); + try testing.expectApproxEqAbs(atanBinary64(-0x1p+0), -0x1.921fb54442d18p-1, math.floatEpsAt(f64, -0x1.921fb54442d18p-1)); + try testing.expectApproxEqAbs(atanBinary64(math.inf(f64)), 0x1.921fb54442d18p+0, math.floatEpsAt(f64, 0x1.921fb54442d18p+0)); + try testing.expectApproxEqAbs(atanBinary64(-math.inf(f64)), -0x1.921fb54442d18p+0, math.floatEpsAt(f64, -0x1.921fb54442d18p+0)); + try testing.expect(math.isNan(atanBinary64(math.nan(f64)))); +} - try expect(math.approxEqAbs(f64, atan64(0.2), 0.197396, epsilon)); - try expect(math.approxEqAbs(f64, atan64(-0.2), -0.197396, epsilon)); - try expect(math.approxEqAbs(f64, atan64(0.3434), 0.330783, epsilon)); - try expect(math.approxEqAbs(f64, atan64(0.8923), 0.728545, epsilon)); - try expect(math.approxEqAbs(f64, atan64(1.5), 0.982794, epsilon)); +test "atanBinary64" { + try testing.expectApproxEqAbs(atanBinary64(-0x1.8629d0244cdccp-2), -0x1.74c61f4377016p-2, math.floatEpsAt(f64, -0x1.74c61f4377016p-2)); + try testing.expectApproxEqAbs(atanBinary64(-0x1.59d42d4659937p1), -0x1.375fd7987cc2p0, math.floatEpsAt(f64, -0x1.375fd7987cc2p0)); + try testing.expectApproxEqAbs(atanBinary64(-0x1.d2dbe23d04f06p0), -0x1.11b8adeba5616p0, math.floatEpsAt(f64, -0x1.11b8adeba5616p0)); + try testing.expectApproxEqAbs(atanBinary64(-0x1.5f314e72398e8p-1), -0x1.33d28ca762539p-1, math.floatEpsAt(f64, -0x1.33d28ca762539p-1)); + try testing.expectApproxEqAbs(atanBinary64(0x1.5869af37b7d08p1), 0x1.37082ce2dd03p0, math.floatEpsAt(f64, 0x1.37082ce2dd03p0)); + try testing.expectApproxEqAbs(atanBinary64(-0x1.b13a05a662618p-2), -0x1.99d7cac66dd44p-2, math.floatEpsAt(f64, -0x1.99d7cac66dd44p-2)); + try testing.expectApproxEqAbs(atanBinary64(0x1.3cb0f12f39d8ap1), 0x1.2fcb120468e8ep0, math.floatEpsAt(f64, 0x1.2fcb120468e8ep0)); + try testing.expectApproxEqAbs(atanBinary64(-0x1.0ed746b39cbb7p-2), -0x1.08c71aa0e509p-2, math.floatEpsAt(f64, -0x1.08c71aa0e509p-2)); + try testing.expectApproxEqAbs(atanBinary64(0x1.299d54ac7d6bp1), 0x1.2a24e22d861dfp0, math.floatEpsAt(f64, 0x1.2a24e22d861dfp0)); + try testing.expectApproxEqAbs(atanBinary64(-0x1.0264fb9f3d50ep1), -0x1.1c617825f9751p0, math.floatEpsAt(f64, -0x1.1c617825f9751p0)); } -test "atan32.special" { - const epsilon = 0.000001; +test "atanExtended80.special" { + try testing.expectEqual(atanExtended80(0x0p+0), 0x0p+0); + try testing.expectEqual(atanExtended80(-0x0p+0), -0x0p+0); + try testing.expectApproxEqAbs(atanExtended80(0x1p+0), 0x1.921fb54442d1846ap-1, math.floatEpsAt(f80, 0x1.921fb54442d1846ap-1)); + try testing.expectApproxEqAbs(atanExtended80(-0x1p+0), -0x1.921fb54442d1846ap-1, math.floatEpsAt(f80, -0x1.921fb54442d1846ap-1)); + try testing.expectApproxEqAbs(atanExtended80(math.inf(f80)), 0x1.921fb54442d1846ap0, math.floatEpsAt(f80, 0x1.921fb54442d1846ap0)); + try testing.expectApproxEqAbs(atanExtended80(-math.inf(f80)), -0x1.921fb54442d1846ap0, math.floatEpsAt(f80, -0x1.921fb54442d1846ap0)); + try testing.expect(math.isNan(atanExtended80(math.nan(f80)))); +} - try expect(math.isPositiveZero(atan32(0.0))); - try expect(math.isNegativeZero(atan32(-0.0))); - try expect(math.approxEqAbs(f32, atan32(math.inf(f32)), math.pi / 2.0, epsilon)); - try expect(math.approxEqAbs(f32, atan32(-math.inf(f32)), -math.pi / 2.0, epsilon)); +test "atanExtended80" { + try testing.expectApproxEqAbs(atanExtended80(-0x1.8629d0244cdcbed8p-2), -0x1.74c61f437701661p-2, math.floatEpsAt(f80, -0x1.74c61f437701661p-2)); + try testing.expectApproxEqAbs(atanExtended80(-0x1.59d42d4659936d9ep1), -0x1.375fd7987cc1fd02p0, math.floatEpsAt(f80, -0x1.375fd7987cc1fd02p0)); + try testing.expectApproxEqAbs(atanExtended80(-0x1.d2dbe23d04f067b4p0), -0x1.11b8adeba5615e04p0, math.floatEpsAt(f80, -0x1.11b8adeba5615e04p0)); + try testing.expectApproxEqAbs(atanExtended80(-0x1.5f314e72398e7dbcp-1), -0x1.33d28ca76253964cp-1, math.floatEpsAt(f80, -0x1.33d28ca76253964cp-1)); + try testing.expectApproxEqAbs(atanExtended80(0x1.5869af37b7d078cap1), 0x1.37082ce2dd03010cp0, math.floatEpsAt(f80, 0x1.37082ce2dd03010cp0)); + try testing.expectApproxEqAbs(atanExtended80(-0x1.b13a05a66261821ap-2), -0x1.99d7cac66dd4438p-2, math.floatEpsAt(f80, -0x1.99d7cac66dd4438p-2)); + try testing.expectApproxEqAbs(atanExtended80(0x1.3cb0f12f39d899cp1), 0x1.2fcb120468e8d9ecp0, math.floatEpsAt(f80, 0x1.2fcb120468e8d9ecp0)); + try testing.expectApproxEqAbs(atanExtended80(-0x1.0ed746b39cbb7614p-2), -0x1.08c71aa0e5090998p-2, math.floatEpsAt(f80, -0x1.08c71aa0e5090998p-2)); + try testing.expectApproxEqAbs(atanExtended80(0x1.299d54ac7d6afc52p1), 0x1.2a24e22d861debfep0, math.floatEpsAt(f80, 0x1.2a24e22d861debfep0)); + try testing.expectApproxEqAbs(atanExtended80(-0x1.0264fb9f3d50e4fp1), -0x1.1c617825f97512b8p0, math.floatEpsAt(f80, -0x1.1c617825f97512b8p0)); } -test "atan64.special" { - const epsilon = 0.000001; +test "atanBinary128.special" { + try testing.expectEqual(atanBinary128(0x0p+0), 0x0p+0); + try testing.expectEqual(atanBinary128(-0x0p+0), -0x0p+0); + try testing.expectApproxEqAbs(atanBinary128(0x1p+0), 0x1.921fb54442d18469898cc51701b8p-1, math.floatEpsAt(f128, 0x1.921fb54442d18469898cc51701b8p-1)); + try testing.expectApproxEqAbs(atanBinary128(-0x1p+0), -0x1.921fb54442d18469898cc51701b8p-1, math.floatEpsAt(f128, -0x1.921fb54442d18469898cc51701b8p-1)); + try testing.expectApproxEqAbs(atanBinary128(math.inf(f128)), 0x1.921fb54442d18469898cc51701b8p0, math.floatEpsAt(f128, 0x1.921fb54442d18469898cc51701b8p0)); + try testing.expectApproxEqAbs(atanBinary128(-math.inf(f128)), -0x1.921fb54442d18469898cc51701b8p0, math.floatEpsAt(f128, -0x1.921fb54442d18469898cc51701b8p0)); + try testing.expect(math.isNan(atanBinary128(math.nan(f128)))); +} - try expect(math.isPositiveZero(atan64(0.0))); - try expect(math.isNegativeZero(atan64(-0.0))); - try expect(math.approxEqAbs(f64, atan64(math.inf(f64)), math.pi / 2.0, epsilon)); - try expect(math.approxEqAbs(f64, atan64(-math.inf(f64)), -math.pi / 2.0, epsilon)); +test "atanBinary128" { + try testing.expectApproxEqAbs(atanBinary128(-0x1.8629d0244cdcbed71792ccdec26dp-2), -0x1.74c61f437701660ff76989d23707p-2, math.floatEpsAt(f128, -0x1.74c61f437701660ff76989d23707p-2)); + try testing.expectApproxEqAbs(atanBinary128(-0x1.59d42d4659936d9e22b5dea4faefp1), -0x1.375fd7987cc1fd0119cf0cc5b708p0, math.floatEpsAt(f128, -0x1.375fd7987cc1fd0119cf0cc5b708p0)); + try testing.expectApproxEqAbs(atanBinary128(-0x1.d2dbe23d04f067b42da3f8efdf57p0), -0x1.11b8adeba5615e0370722b511231p0, math.floatEpsAt(f128, -0x1.11b8adeba5615e0370722b511231p0)); + try testing.expectApproxEqAbs(atanBinary128(-0x1.5f314e72398e7dbbe70fb072983ep-1), -0x1.33d28ca76253964cb5d3581cdd88p-1, math.floatEpsAt(f128, -0x1.33d28ca76253964cb5d3581cdd88p-1)); + try testing.expectApproxEqAbs(atanBinary128(0x1.5869af37b7d078caa3456c44aecep1), 0x1.37082ce2dd03010bbea814dc5882p0, math.floatEpsAt(f128, 0x1.37082ce2dd03010bbea814dc5882p0)); + try testing.expectApproxEqAbs(atanBinary128(-0x1.b13a05a66261821a364ad8c6c999p-2), -0x1.99d7cac66dd4438077284b491a91p-2, math.floatEpsAt(f128, -0x1.99d7cac66dd4438077284b491a91p-2)); + try testing.expectApproxEqAbs(atanBinary128(0x1.3cb0f12f39d899c0d963ac413297p1), 0x1.2fcb120468e8d9ebdb74702314c8p0, math.floatEpsAt(f128, 0x1.2fcb120468e8d9ebdb74702314c8p0)); + try testing.expectApproxEqAbs(atanBinary128(-0x1.0ed746b39cbb7614d8735e8315a8p-2), -0x1.08c71aa0e5090998206fbbe2090fp-2, math.floatEpsAt(f128, -0x1.08c71aa0e5090998206fbbe2090fp-2)); + try testing.expectApproxEqAbs(atanBinary128(0x1.299d54ac7d6afc5154643b601519p1), 0x1.2a24e22d861debfd6f974500567fp0, math.floatEpsAt(f128, 0x1.2a24e22d861debfd6f974500567fp0)); + try testing.expectApproxEqAbs(atanBinary128(-0x1.0264fb9f3d50e4f0f966f0686064p1), -0x1.1c617825f97512b7f38656ab12cdp0, math.floatEpsAt(f128, -0x1.1c617825f97512b7f38656ab12cdp0)); } diff --git a/src/libs/musl.zig b/src/libs/musl.zig index 3c2c4dd25fed0d205bbbf9d9be0b9f3fb78467a3..28602cd6f6ec14f15f586c01bf857234bf1ce8b7 100644 --- a/src/libs/musl.zig +++ b/src/libs/musl.zig @@ -836,12 +836,9 @@ const src_files = [_][]const u8{ "musl/src/math/atan2.c", "musl/src/math/atan2f.c", "musl/src/math/atan2l.c", - "musl/src/math/atan.c", - "musl/src/math/atanf.c", "musl/src/math/atanh.c", "musl/src/math/atanhf.c", "musl/src/math/atanhl.c", - "musl/src/math/atanl.c", "musl/src/math/cbrt.c", "musl/src/math/cbrtf.c", "musl/src/math/cbrtl.c", @@ -892,9 +889,6 @@ const src_files = [_][]const u8{ "musl/src/math/i386/atan2f.s", "musl/src/math/i386/atan2l.s", "musl/src/math/i386/atan2.s", - "musl/src/math/i386/atanf.s", - "musl/src/math/i386/atanl.s", - "musl/src/math/i386/atan.s", "musl/src/math/i386/exp2l.s", "musl/src/math/i386/exp_ld.s", "musl/src/math/i386/expl.s", @@ -1076,7 +1070,6 @@ const src_files = [_][]const u8{ "musl/src/math/x32/acosl.s", "musl/src/math/x32/asinl.s", "musl/src/math/x32/atan2l.s", - "musl/src/math/x32/atanl.s", "musl/src/math/x32/exp2l.s", "musl/src/math/x32/expl.s", "musl/src/math/x32/expm1l.s", @@ -1098,7 +1091,6 @@ const src_files = [_][]const u8{ "musl/src/math/x86_64/acosl.s", "musl/src/math/x86_64/asinl.s", "musl/src/math/x86_64/atan2l.s", - "musl/src/math/x86_64/atanl.s", "musl/src/math/x86_64/exp2l.s", "musl/src/math/x86_64/expl.s", "musl/src/math/x86_64/expm1l.s", diff --git a/src/libs/wasi_libc.zig b/src/libs/wasi_libc.zig index 61d0f2cd0cb2d432d130aecd4b5f2b97fef1ad81..b1a1159778257d4b16297f1ce87f5431635dad23 100644 --- a/src/libs/wasi_libc.zig +++ b/src/libs/wasi_libc.zig @@ -698,12 +698,9 @@ const libc_top_half_src_files = [_][]const u8{ "musl/src/math/atan2.c", "musl/src/math/atan2f.c", "musl/src/math/atan2l.c", - "musl/src/math/atan.c", - "musl/src/math/atanf.c", "musl/src/math/atanh.c", "musl/src/math/atanhf.c", "musl/src/math/atanhl.c", - "musl/src/math/atanl.c", "musl/src/math/cbrt.c", "musl/src/math/cbrtf.c", "musl/src/math/cbrtl.c", -- 2.54.0