authorgravatar for 3575188313@qq.comZhenming-Lin <3575188313@qq.com> 2026-01-30 07:18:13+08:00
committergravatar for andrew@ziglang.orgAndrew Kelley <andrew@ziglang.org> 2026-01-30 06:20:14+01:00
logccd82ae7cc5ce3f708a919a070bffc7572416c62
tree997e5a08c1e6637a8a6e77bbcfab40d57395a8f5
parentad0458f5826a9283df3e152fea24869d0914634d

Add `f16`, `f80` and `f128` support for `atan`

12 files changed, 540 insertions(+), 652 deletions(-)

lib/c/math.zig+45-22
......@@ -1,4 +1,5 @@
11const std = @import("std");
2const math = std.math;
23const common = @import("common.zig");
34const builtin = @import("builtin");
45
......@@ -11,20 +12,20 @@ comptime {
1112 @export(&isnanl, .{ .name = "isnanl", .linkage = common.linkage, .visibility = common.visibility });
1213 @export(&isnanl, .{ .name = "__isnanl", .linkage = common.linkage, .visibility = common.visibility });
1314
14 @export(&std.math.nan(f64), .{ .name = "__QNAN", .linkage = common.linkage, .visibility = common.visibility });
15 @export(&std.math.snan(f64), .{ .name = "__SNAN", .linkage = common.linkage, .visibility = common.visibility });
16 @export(&std.math.inf(f64), .{ .name = "__INF", .linkage = common.linkage, .visibility = common.visibility });
17 @export(&std.math.floatTrueMin(f64), .{ .name = "__DENORM", .linkage = common.linkage, .visibility = common.visibility });
15 @export(&math.nan(f64), .{ .name = "__QNAN", .linkage = common.linkage, .visibility = common.visibility });
16 @export(&math.snan(f64), .{ .name = "__SNAN", .linkage = common.linkage, .visibility = common.visibility });
17 @export(&math.inf(f64), .{ .name = "__INF", .linkage = common.linkage, .visibility = common.visibility });
18 @export(&math.floatTrueMin(f64), .{ .name = "__DENORM", .linkage = common.linkage, .visibility = common.visibility });
1819
19 @export(&std.math.nan(f32), .{ .name = "__QNANF", .linkage = common.linkage, .visibility = common.visibility });
20 @export(&std.math.snan(f32), .{ .name = "__SNANF", .linkage = common.linkage, .visibility = common.visibility });
21 @export(&std.math.inf(f32), .{ .name = "__INFF", .linkage = common.linkage, .visibility = common.visibility });
22 @export(&std.math.floatTrueMin(f32), .{ .name = "__DENORMF", .linkage = common.linkage, .visibility = common.visibility });
20 @export(&math.nan(f32), .{ .name = "__QNANF", .linkage = common.linkage, .visibility = common.visibility });
21 @export(&math.snan(f32), .{ .name = "__SNANF", .linkage = common.linkage, .visibility = common.visibility });
22 @export(&math.inf(f32), .{ .name = "__INFF", .linkage = common.linkage, .visibility = common.visibility });
23 @export(&math.floatTrueMin(f32), .{ .name = "__DENORMF", .linkage = common.linkage, .visibility = common.visibility });
2324
24 @export(&std.math.nan(c_longdouble), .{ .name = "__QNANL", .linkage = common.linkage, .visibility = common.visibility });
25 @export(&std.math.snan(c_longdouble), .{ .name = "__SNANL", .linkage = common.linkage, .visibility = common.visibility });
26 @export(&std.math.inf(c_longdouble), .{ .name = "__INFL", .linkage = common.linkage, .visibility = common.visibility });
27 @export(&std.math.floatTrueMin(c_longdouble), .{ .name = "__DENORML", .linkage = common.linkage, .visibility = common.visibility });
25 @export(&math.nan(c_longdouble), .{ .name = "__QNANL", .linkage = common.linkage, .visibility = common.visibility });
26 @export(&math.snan(c_longdouble), .{ .name = "__SNANL", .linkage = common.linkage, .visibility = common.visibility });
27 @export(&math.inf(c_longdouble), .{ .name = "__INFL", .linkage = common.linkage, .visibility = common.visibility });
28 @export(&math.floatTrueMin(c_longdouble), .{ .name = "__DENORML", .linkage = common.linkage, .visibility = common.visibility });
2829 }
2930
3031 if (builtin.target.isMinGW() or builtin.target.isMuslLibC() or builtin.target.isWasiLibC()) {
......@@ -35,6 +36,9 @@ comptime {
3536
3637 if (builtin.target.isMuslLibC() or builtin.target.isWasiLibC()) {
3738 @export(&acos, .{ .name = "acos", .linkage = common.linkage, .visibility = common.visibility });
39 @export(&atanf, .{ .name = "atanf", .linkage = common.linkage, .visibility = common.visibility });
40 @export(&atan, .{ .name = "atan", .linkage = common.linkage, .visibility = common.visibility });
41 @export(&atanl, .{ .name = "atanl", .linkage = common.linkage, .visibility = common.visibility });
3842 }
3943
4044 if (builtin.target.isMuslLibC()) {
......@@ -45,41 +49,60 @@ comptime {
4549}
4650
4751fn acos(x: f64) callconv(.c) f64 {
48 return std.math.acos(x);
52 return math.acos(x);
53}
54
55fn atanf(x: f32) callconv(.c) f32 {
56 return math.atan(x);
57}
58
59fn atan(x: f64) callconv(.c) f64 {
60 return math.atan(x);
61}
62
63fn atanl(x: c_longdouble) callconv(.c) c_longdouble {
64 return switch (@typeInfo(@TypeOf(x)).float.bits) {
65 16 => math.atan(@as(f16, @floatCast(x))),
66 32 => math.atan(@as(f32, @floatCast(x))),
67 64 => math.atan(@as(f64, @floatCast(x))),
68 80 => math.atan(@as(f80, @floatCast(x))),
69 128 => math.atan(@as(f128, @floatCast(x))),
70 else => unreachable,
71 };
4972}
5073
5174fn isnan(x: f64) callconv(.c) c_int {
52 return if (std.math.isNan(x)) 1 else 0;
75 return if (math.isNan(x)) 1 else 0;
5376}
5477
5578fn isnanf(x: f32) callconv(.c) c_int {
56 return if (std.math.isNan(x)) 1 else 0;
79 return if (math.isNan(x)) 1 else 0;
5780}
5881
5982fn isnanl(x: c_longdouble) callconv(.c) c_int {
60 return if (std.math.isNan(x)) 1 else 0;
83 return if (math.isNan(x)) 1 else 0;
6184}
6285
6386fn nan(_: [*:0]const c_char) callconv(.c) f64 {
64 return std.math.nan(f64);
87 return math.nan(f64);
6588}
6689
6790fn nanf(_: [*:0]const c_char) callconv(.c) f32 {
68 return std.math.nan(f32);
91 return math.nan(f32);
6992}
7093
7194fn nanl(_: [*:0]const c_char) callconv(.c) c_longdouble {
72 return std.math.nan(c_longdouble);
95 return math.nan(c_longdouble);
7396}
7497
7598fn copysignf(x: f32, y: f32) callconv(.c) f32 {
76 return std.math.copysign(x, y);
99 return math.copysign(x, y);
77100}
78101
79102fn copysign(x: f64, y: f64) callconv(.c) f64 {
80 return std.math.copysign(x, y);
103 return math.copysign(x, y);
81104}
82105
83106fn copysignl(x: c_longdouble, y: c_longdouble) callconv(.c) c_longdouble {
84 return std.math.copysign(x, y);
107 return math.copysign(x, y);
85108}
lib/libc/musl/src/math/atan.c deleted-116
......@@ -1,116 +0,0 @@
1/* origin: FreeBSD /usr/src/lib/msun/src/s_atan.c */
2/*
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5 *
6 * Developed at SunPro, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice
9 * is preserved.
10 * ====================================================
11 */
12/* atan(x)
13 * Method
14 * 1. Reduce x to positive by atan(x) = -atan(-x).
15 * 2. According to the integer k=4t+0.25 chopped, t=x, the argument
16 * is further reduced to one of the following intervals and the
17 * arctangent of t is evaluated by the corresponding formula:
18 *
19 * [0,7/16] atan(x) = t-t^3*(a1+t^2*(a2+...(a10+t^2*a11)...)
20 * [7/16,11/16] atan(x) = atan(1/2) + atan( (t-0.5)/(1+t/2) )
21 * [11/16.19/16] atan(x) = atan( 1 ) + atan( (t-1)/(1+t) )
22 * [19/16,39/16] atan(x) = atan(3/2) + atan( (t-1.5)/(1+1.5t) )
23 * [39/16,INF] atan(x) = atan(INF) + atan( -1/t )
24 *
25 * Constants:
26 * The hexadecimal values are the intended ones for the following
27 * constants. The decimal values may be used, provided that the
28 * compiler will convert from decimal to binary accurately enough
29 * to produce the hexadecimal values shown.
30 */
31
32
33#include "libm.h"
34
35static const double atanhi[] = {
36 4.63647609000806093515e-01, /* atan(0.5)hi 0x3FDDAC67, 0x0561BB4F */
37 7.85398163397448278999e-01, /* atan(1.0)hi 0x3FE921FB, 0x54442D18 */
38 9.82793723247329054082e-01, /* atan(1.5)hi 0x3FEF730B, 0xD281F69B */
39 1.57079632679489655800e+00, /* atan(inf)hi 0x3FF921FB, 0x54442D18 */
40};
41
42static const double atanlo[] = {
43 2.26987774529616870924e-17, /* atan(0.5)lo 0x3C7A2B7F, 0x222F65E2 */
44 3.06161699786838301793e-17, /* atan(1.0)lo 0x3C81A626, 0x33145C07 */
45 1.39033110312309984516e-17, /* atan(1.5)lo 0x3C700788, 0x7AF0CBBD */
46 6.12323399573676603587e-17, /* atan(inf)lo 0x3C91A626, 0x33145C07 */
47};
48
49static const double aT[] = {
50 3.33333333333329318027e-01, /* 0x3FD55555, 0x5555550D */
51 -1.99999999998764832476e-01, /* 0xBFC99999, 0x9998EBC4 */
52 1.42857142725034663711e-01, /* 0x3FC24924, 0x920083FF */
53 -1.11111104054623557880e-01, /* 0xBFBC71C6, 0xFE231671 */
54 9.09088713343650656196e-02, /* 0x3FB745CD, 0xC54C206E */
55 -7.69187620504482999495e-02, /* 0xBFB3B0F2, 0xAF749A6D */
56 6.66107313738753120669e-02, /* 0x3FB10D66, 0xA0D03D51 */
57 -5.83357013379057348645e-02, /* 0xBFADDE2D, 0x52DEFD9A */
58 4.97687799461593236017e-02, /* 0x3FA97B4B, 0x24760DEB */
59 -3.65315727442169155270e-02, /* 0xBFA2B444, 0x2C6A6C2F */
60 1.62858201153657823623e-02, /* 0x3F90AD3A, 0xE322DA11 */
61};
62
63double atan(double x)
64{
65 double_t w,s1,s2,z;
66 uint32_t ix,sign;
67 int id;
68
69 GET_HIGH_WORD(ix, x);
70 sign = ix >> 31;
71 ix &= 0x7fffffff;
72 if (ix >= 0x44100000) { /* if |x| >= 2^66 */
73 if (isnan(x))
74 return x;
75 z = atanhi[3] + 0x1p-120f;
76 return sign ? -z : z;
77 }
78 if (ix < 0x3fdc0000) { /* |x| < 0.4375 */
79 if (ix < 0x3e400000) { /* |x| < 2^-27 */
80 if (ix < 0x00100000)
81 /* raise underflow for subnormal x */
82 FORCE_EVAL((float)x);
83 return x;
84 }
85 id = -1;
86 } else {
87 x = fabs(x);
88 if (ix < 0x3ff30000) { /* |x| < 1.1875 */
89 if (ix < 0x3fe60000) { /* 7/16 <= |x| < 11/16 */
90 id = 0;
91 x = (2.0*x-1.0)/(2.0+x);
92 } else { /* 11/16 <= |x| < 19/16 */
93 id = 1;
94 x = (x-1.0)/(x+1.0);
95 }
96 } else {
97 if (ix < 0x40038000) { /* |x| < 2.4375 */
98 id = 2;
99 x = (x-1.5)/(1.0+1.5*x);
100 } else { /* 2.4375 <= |x| < 2^66 */
101 id = 3;
102 x = -1.0/x;
103 }
104 }
105 }
106 /* end of argument reduction */
107 z = x*x;
108 w = z*z;
109 /* break sum from i=0 to 10 aT[i]z**(i+1) into odd and even poly */
110 s1 = z*(aT[0]+w*(aT[2]+w*(aT[4]+w*(aT[6]+w*(aT[8]+w*aT[10])))));
111 s2 = w*(aT[1]+w*(aT[3]+w*(aT[5]+w*(aT[7]+w*aT[9]))));
112 if (id < 0)
113 return x - x*(s1+s2);
114 z = atanhi[id] - (x*(s1+s2) - atanlo[id] - x);
115 return sign ? -z : z;
116}
lib/libc/musl/src/math/atanf.c deleted-94
......@@ -1,94 +0,0 @@
1/* origin: FreeBSD /usr/src/lib/msun/src/s_atanf.c */
2/*
3 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
4 */
5/*
6 * ====================================================
7 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
8 *
9 * Developed at SunPro, a Sun Microsystems, Inc. business.
10 * Permission to use, copy, modify, and distribute this
11 * software is freely granted, provided that this notice
12 * is preserved.
13 * ====================================================
14 */
15
16
17#include "libm.h"
18
19static const float atanhi[] = {
20 4.6364760399e-01, /* atan(0.5)hi 0x3eed6338 */
21 7.8539812565e-01, /* atan(1.0)hi 0x3f490fda */
22 9.8279368877e-01, /* atan(1.5)hi 0x3f7b985e */
23 1.5707962513e+00, /* atan(inf)hi 0x3fc90fda */
24};
25
26static const float atanlo[] = {
27 5.0121582440e-09, /* atan(0.5)lo 0x31ac3769 */
28 3.7748947079e-08, /* atan(1.0)lo 0x33222168 */
29 3.4473217170e-08, /* atan(1.5)lo 0x33140fb4 */
30 7.5497894159e-08, /* atan(inf)lo 0x33a22168 */
31};
32
33static const float aT[] = {
34 3.3333328366e-01,
35 -1.9999158382e-01,
36 1.4253635705e-01,
37 -1.0648017377e-01,
38 6.1687607318e-02,
39};
40
41float atanf(float x)
42{
43 float_t w,s1,s2,z;
44 uint32_t ix,sign;
45 int id;
46
47 GET_FLOAT_WORD(ix, x);
48 sign = ix>>31;
49 ix &= 0x7fffffff;
50 if (ix >= 0x4c800000) { /* if |x| >= 2**26 */
51 if (isnan(x))
52 return x;
53 z = atanhi[3] + 0x1p-120f;
54 return sign ? -z : z;
55 }
56 if (ix < 0x3ee00000) { /* |x| < 0.4375 */
57 if (ix < 0x39800000) { /* |x| < 2**-12 */
58 if (ix < 0x00800000)
59 /* raise underflow for subnormal x */
60 FORCE_EVAL(x*x);
61 return x;
62 }
63 id = -1;
64 } else {
65 x = fabsf(x);
66 if (ix < 0x3f980000) { /* |x| < 1.1875 */
67 if (ix < 0x3f300000) { /* 7/16 <= |x| < 11/16 */
68 id = 0;
69 x = (2.0f*x - 1.0f)/(2.0f + x);
70 } else { /* 11/16 <= |x| < 19/16 */
71 id = 1;
72 x = (x - 1.0f)/(x + 1.0f);
73 }
74 } else {
75 if (ix < 0x401c0000) { /* |x| < 2.4375 */
76 id = 2;
77 x = (x - 1.5f)/(1.0f + 1.5f*x);
78 } else { /* 2.4375 <= |x| < 2**26 */
79 id = 3;
80 x = -1.0f/x;
81 }
82 }
83 }
84 /* end of argument reduction */
85 z = x*x;
86 w = z*z;
87 /* break sum from i=0 to 10 aT[i]z**(i+1) into odd and even poly */
88 s1 = z*(aT[0]+w*(aT[2]+w*aT[4]));
89 s2 = w*(aT[1]+w*aT[3]);
90 if (id < 0)
91 return x - x*(s1+s2);
92 z = atanhi[id] - ((x*(s1+s2) - atanlo[id]) - x);
93 return sign ? -z : z;
94}
lib/libc/musl/src/math/atanl.c deleted-184
......@@ -1,184 +0,0 @@
1/* origin: FreeBSD /usr/src/lib/msun/src/s_atanl.c */
2/*
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5 *
6 * Developed at SunPro, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice
9 * is preserved.
10 * ====================================================
11 */
12/*
13 * See comments in atan.c.
14 * Converted to long double by David Schultz <das@FreeBSD.ORG>.
15 */
16
17#include "libm.h"
18
19#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
20long double atanl(long double x)
21{
22 return atan(x);
23}
24#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
25
26#if LDBL_MANT_DIG == 64
27#define EXPMAN(u) ((u.i.se & 0x7fff)<<8 | (u.i.m>>55 & 0xff))
28
29static const long double atanhi[] = {
30 4.63647609000806116202e-01L,
31 7.85398163397448309628e-01L,
32 9.82793723247329067960e-01L,
33 1.57079632679489661926e+00L,
34};
35
36static const long double atanlo[] = {
37 1.18469937025062860669e-20L,
38 -1.25413940316708300586e-20L,
39 2.55232234165405176172e-20L,
40 -2.50827880633416601173e-20L,
41};
42
43static const long double aT[] = {
44 3.33333333333333333017e-01L,
45 -1.99999999999999632011e-01L,
46 1.42857142857046531280e-01L,
47 -1.11111111100562372733e-01L,
48 9.09090902935647302252e-02L,
49 -7.69230552476207730353e-02L,
50 6.66661718042406260546e-02L,
51 -5.88158892835030888692e-02L,
52 5.25499891539726639379e-02L,
53 -4.70119845393155721494e-02L,
54 4.03539201366454414072e-02L,
55 -2.91303858419364158725e-02L,
56 1.24822046299269234080e-02L,
57};
58
59static long double T_even(long double x)
60{
61 return aT[0] + x * (aT[2] + x * (aT[4] + x * (aT[6] +
62 x * (aT[8] + x * (aT[10] + x * aT[12])))));
63}
64
65static long double T_odd(long double x)
66{
67 return aT[1] + x * (aT[3] + x * (aT[5] + x * (aT[7] +
68 x * (aT[9] + x * aT[11]))));
69}
70#elif LDBL_MANT_DIG == 113
71#define EXPMAN(u) ((u.i.se & 0x7fff)<<8 | u.i.top>>8)
72
73static const long double atanhi[] = {
74 4.63647609000806116214256231461214397e-01L,
75 7.85398163397448309615660845819875699e-01L,
76 9.82793723247329067985710611014666038e-01L,
77 1.57079632679489661923132169163975140e+00L,
78};
79
80static const long double atanlo[] = {
81 4.89509642257333492668618435220297706e-36L,
82 2.16795253253094525619926100651083806e-35L,
83 -2.31288434538183565909319952098066272e-35L,
84 4.33590506506189051239852201302167613e-35L,
85};
86
87static const long double aT[] = {
88 3.33333333333333333333333333333333125e-01L,
89 -1.99999999999999999999999999999180430e-01L,
90 1.42857142857142857142857142125269827e-01L,
91 -1.11111111111111111111110834490810169e-01L,
92 9.09090909090909090908522355708623681e-02L,
93 -7.69230769230769230696553844935357021e-02L,
94 6.66666666666666660390096773046256096e-02L,
95 -5.88235294117646671706582985209643694e-02L,
96 5.26315789473666478515847092020327506e-02L,
97 -4.76190476189855517021024424991436144e-02L,
98 4.34782608678695085948531993458097026e-02L,
99 -3.99999999632663469330634215991142368e-02L,
100 3.70370363987423702891250829918659723e-02L,
101 -3.44827496515048090726669907612335954e-02L,
102 3.22579620681420149871973710852268528e-02L,
103 -3.03020767654269261041647570626778067e-02L,
104 2.85641979882534783223403715930946138e-02L,
105 -2.69824879726738568189929461383741323e-02L,
106 2.54194698498808542954187110873675769e-02L,
107 -2.35083879708189059926183138130183215e-02L,
108 2.04832358998165364349957325067131428e-02L,
109 -1.54489555488544397858507248612362957e-02L,
110 8.64492360989278761493037861575248038e-03L,
111 -2.58521121597609872727919154569765469e-03L,
112};
113
114static long double T_even(long double x)
115{
116 return (aT[0] + x * (aT[2] + x * (aT[4] + x * (aT[6] + x * (aT[8] +
117 x * (aT[10] + x * (aT[12] + x * (aT[14] + x * (aT[16] +
118 x * (aT[18] + x * (aT[20] + x * aT[22])))))))))));
119}
120
121static long double T_odd(long double x)
122{
123 return (aT[1] + x * (aT[3] + x * (aT[5] + x * (aT[7] + x * (aT[9] +
124 x * (aT[11] + x * (aT[13] + x * (aT[15] + x * (aT[17] +
125 x * (aT[19] + x * (aT[21] + x * aT[23])))))))))));
126}
127#endif
128
129long double atanl(long double x)
130{
131 union ldshape u = {x};
132 long double w, s1, s2, z;
133 int id;
134 unsigned e = u.i.se & 0x7fff;
135 unsigned sign = u.i.se >> 15;
136 unsigned expman;
137
138 if (e >= 0x3fff + LDBL_MANT_DIG + 1) { /* if |x| is large, atan(x)~=pi/2 */
139 if (isnan(x))
140 return x;
141 return sign ? -atanhi[3] : atanhi[3];
142 }
143 /* Extract the exponent and the first few bits of the mantissa. */
144 expman = EXPMAN(u);
145 if (expman < ((0x3fff - 2) << 8) + 0xc0) { /* |x| < 0.4375 */
146 if (e < 0x3fff - (LDBL_MANT_DIG+1)/2) { /* if |x| is small, atanl(x)~=x */
147 /* raise underflow if subnormal */
148 if (e == 0)
149 FORCE_EVAL((float)x);
150 return x;
151 }
152 id = -1;
153 } else {
154 x = fabsl(x);
155 if (expman < (0x3fff << 8) + 0x30) { /* |x| < 1.1875 */
156 if (expman < ((0x3fff - 1) << 8) + 0x60) { /* 7/16 <= |x| < 11/16 */
157 id = 0;
158 x = (2.0*x-1.0)/(2.0+x);
159 } else { /* 11/16 <= |x| < 19/16 */
160 id = 1;
161 x = (x-1.0)/(x+1.0);
162 }
163 } else {
164 if (expman < ((0x3fff + 1) << 8) + 0x38) { /* |x| < 2.4375 */
165 id = 2;
166 x = (x-1.5)/(1.0+1.5*x);
167 } else { /* 2.4375 <= |x| */
168 id = 3;
169 x = -1.0/x;
170 }
171 }
172 }
173 /* end of argument reduction */
174 z = x*x;
175 w = z*z;
176 /* break sum aT[i]z**(i+1) into odd and even poly */
177 s1 = z*T_even(w);
178 s2 = w*T_odd(w);
179 if (id < 0)
180 return x - x*(s1+s2);
181 z = atanhi[id] - ((x*(s1+s2) - atanlo[id]) - x);
182 return sign ? -z : z;
183}
184#endif
lib/libc/musl/src/math/i386/atan.s deleted-16
......@@ -1,16 +0,0 @@
1.global atan
2.type atan,@function
3atan:
4 fldl 4(%esp)
5 mov 8(%esp),%eax
6 add %eax,%eax
7 cmp $0x00200000,%eax
8 jb 1f
9 fld1
10 fpatan
11 fstpl 4(%esp)
12 fldl 4(%esp)
13 ret
14 # subnormal x, return x with underflow
151: fsts 4(%esp)
16 ret
lib/libc/musl/src/math/i386/atanf.s deleted-18
......@@ -1,18 +0,0 @@
1.global atanf
2.type atanf,@function
3atanf:
4 flds 4(%esp)
5 mov 4(%esp),%eax
6 add %eax,%eax
7 cmp $0x01000000,%eax
8 jb 1f
9 fld1
10 fpatan
11 fstps 4(%esp)
12 flds 4(%esp)
13 ret
14 # subnormal x, return x with underflow
151: fld %st(0)
16 fmul %st(1)
17 fstps 4(%esp)
18 ret
lib/libc/musl/src/math/i386/atanl.s deleted-7
......@@ -1,7 +0,0 @@
1.global atanl
2.type atanl,@function
3atanl:
4 fldt 4(%esp)
5 fld1
6 fpatan
7 ret
lib/libc/musl/src/math/x32/atanl.s deleted-7
......@@ -1,7 +0,0 @@
1.global atanl
2.type atanl,@function
3atanl:
4 fldt 8(%esp)
5 fld1
6 fpatan
7 ret
lib/libc/musl/src/math/x86_64/atanl.s deleted-7
......@@ -1,7 +0,0 @@
1.global atanl
2.type atanl,@function
3atanl:
4 fldt 8(%rsp)
5 fld1
6 fpatan
7 ret
lib/std/math/atan.zig+495-170
......@@ -3,11 +3,12 @@
33//
44// https://git.musl-libc.org/cgit/musl/tree/src/math/atanf.c
55// https://git.musl-libc.org/cgit/musl/tree/src/math/atan.c
6// https://git.musl-libc.org/cgit/musl/tree/src/math/atanl.c
67
78const std = @import("../std.zig");
89const math = std.math;
910const mem = std.mem;
10const expect = std.testing.expect;
11const testing = std.testing;
1112
1213/// Returns the arc-tangent of x.
1314///
......@@ -17,28 +18,100 @@ const expect = std.testing.expect;
1718pub fn atan(x: anytype) @TypeOf(x) {
1819 const T = @TypeOf(x);
1920 return switch (T) {
20 f32 => atan32(x),
21 f64 => atan64(x),
21 f16 => atanBinary16(x),
22 f32 => atanBinary32(x),
23 f64 => atanBinary64(x),
24 f80 => atanExtended80(x),
25 f128 => atanBinary128(x),
2226 else => @compileError("atan not implemented for " ++ @typeName(T)),
2327 };
2428}
2529
26fn atan32(x_: f32) f32 {
27 const atanhi = [_]f32{
28 4.6364760399e-01, // atan(0.5)hi
29 7.8539812565e-01, // atan(1.0)hi
30 9.8279368877e-01, // atan(1.5)hi
31 1.5707962513e+00, // atan(inf)hi
30fn atanBinary16(x: f16) f16 {
31 const atanhi: []const f32 = &.{
32 4.6364760399e-01, // atan(0.5)hi 0x3eed6338
33 7.8539812565e-01, // atan(1.0)hi 0x3f490fda
34 9.8279368877e-01, // atan(1.5)hi 0x3f7b985e
35 1.5707962513e+00, // atan(inf)hi 0x3fc90fda
36 };
37 const aT: []const f32 = &.{
38 0x1.fffcccp-1,
39 -0x1.52e8ccp-2,
40 0x1.522336p-3,
3241 };
3342
34 const atanlo = [_]f32{
35 5.0121582440e-09, // atan(0.5)lo
36 3.7748947079e-08, // atan(1.0)lo
37 3.4473217170e-08, // atan(1.5)lo
38 7.5497894159e-08, // atan(inf)lo
43 const hx: u16 = @bitCast(x);
44 const ix = hx & 0x7fff;
45 const sign = (hx >> 15) != 0;
46 // if |x| >= 2^11
47 if (ix >= 0x6800) {
48 if (math.isNan(x)) {
49 return x;
50 }
51 const z = atanhi[3] + 0x1p-120;
52 return @floatCast(if (sign) -z else z);
53 }
54 const x_: f32, const id: ?usize = blk: {
55 // |x| < 0.4375
56 if (ix < 0x3700) {
57 // |x| < 2^(-6)
58 if (ix < 0x2400) {
59 if (ix < 0x400) {
60 // raise underflow for subnormal x
61 mem.doNotOptimizeAway(x * x);
62 }
63 return x;
64 }
65 break :blk .{ @floatCast(x), null };
66 } else {
67 const x_: f32 = @floatCast(@abs(x));
68 // |x| < 1.1875
69 if (ix < 0x3cc0) {
70 // 7/16 <= |x| < 11/16
71 if (ix < 0x3980) {
72 break :blk .{ (2.0 * x_ - 1.0) / (2.0 + x_), 0 };
73 }
74 // 11/16 <= |x| < 19/16
75 else {
76 break :blk .{ (x_ - 1.0) / (x_ + 1.0), 1 };
77 }
78 } else {
79 // |x| < 2.4375
80 if (ix < 0x40e0) {
81 break :blk .{ (x_ - 1.5) / (1.0 + 1.5 * x_), 2 };
82 }
83 // 2.4375 <= |x| < 2^11
84 else {
85 break :blk .{ -1.0 / x_, 3 };
86 }
87 }
88 }
3989 };
90 // end of argument reduction
91 const z = x_ * x_;
92 const s = aT[0] + z * (aT[1] + z * aT[2]);
93 if (id) |id_| {
94 const z_ = atanhi[id_] + x_ * s;
95 return @floatCast(if (sign) -z_ else z_);
96 } else {
97 return @floatCast(x_ * s);
98 }
99}
40100
41 const aT = [_]f32{
101fn atanBinary32(x: f32) f32 {
102 const atanhi: []const f32 = &.{
103 4.6364760399e-01, // atan(0.5)hi 0x3eed6338
104 7.8539812565e-01, // atan(1.0)hi 0x3f490fda
105 9.8279368877e-01, // atan(1.5)hi 0x3f7b985e
106 1.5707962513e+00, // atan(inf)hi 0x3fc90fda
107 };
108 const atanlo: []const f32 = &.{
109 5.0121582440e-09, // atan(0.5)lo 0x31ac3769
110 3.7748947079e-08, // atan(1.0)lo 0x33222168
111 3.4473217170e-08, // atan(1.5)lo 0x33140fb4
112 7.5497894159e-08, // atan(inf)lo 0x33a22168
113 };
114 const aT: []const f32 = &.{
42115 3.3333328366e-01,
43116 -1.9999158382e-01,
44117 1.4253635705e-01,
......@@ -46,211 +119,463 @@ fn atan32(x_: f32) f32 {
46119 6.1687607318e-02,
47120 };
48121
49 var x = x_;
50 var ix: u32 = @as(u32, @bitCast(x));
51 const sign = ix >> 31;
52 ix &= 0x7FFFFFFF;
53
54 // |x| >= 2^26
55 if (ix >= 0x4C800000) {
122 const hx: u32 = @bitCast(x);
123 const ix = hx & 0x7fff_ffff;
124 const sign = (hx >> 31) != 0;
125 // if |x| >= 2^26
126 if (ix >= 0x4c80_0000) {
56127 if (math.isNan(x)) {
57128 return x;
58 } else {
59 const z = atanhi[3] + 0x1.0p-120;
60 return if (sign != 0) -z else z;
61129 }
130 const z = atanhi[3] + 0x1p-120;
131 return if (sign) -z else z;
62132 }
63
64 var id: ?usize = undefined;
65
66 // |x| < 0.4375
67 if (ix < 0x3EE00000) {
68 // |x| < 2^(-12)
69 if (ix < 0x39800000) {
70 if (ix < 0x00800000) {
71 mem.doNotOptimizeAway(x * x);
72 }
73 return x;
74 }
75 id = null;
76 } else {
77 x = @abs(x);
78 // |x| < 1.1875
79 if (ix < 0x3F980000) {
80 // 7/16 <= |x| < 11/16
81 if (ix < 0x3F300000) {
82 id = 0;
83 x = (2.0 * x - 1.0) / (2.0 + x);
84 }
85 // 11/16 <= |x| < 19/16
86 else {
87 id = 1;
88 x = (x - 1.0) / (x + 1.0);
133 const x_, const id: ?usize = blk: {
134 // |x| < 0.4375
135 if (ix < 0x3ee00000) {
136 // |x| < 2^(-12)
137 if (ix < 0x39800000) {
138 if (ix < 0x00800000) {
139 // raise underflow for subnormal x
140 mem.doNotOptimizeAway(x * x);
141 }
142 return x;
89143 }
144 break :blk .{ x, null };
90145 } else {
91 // |x| < 2.4375
92 if (ix < 0x401C0000) {
93 id = 2;
94 x = (x - 1.5) / (1.0 + 1.5 * x);
95 }
96 // 2.4375 <= |x| < 2^26
97 else {
98 id = 3;
99 x = -1.0 / x;
146 const x_ = @abs(x);
147 // |x| < 1.1875
148 if (ix < 0x3f98_0000) {
149 // 7/16 <= |x| < 11/16
150 if (ix < 0x3f30_0000) {
151 break :blk .{ (2.0 * x_ - 1.0) / (2.0 + x_), 0 };
152 }
153 // 11/16 <= |x| < 19/16
154 else {
155 break :blk .{ (x_ - 1.0) / (x_ + 1.0), 1 };
156 }
157 } else {
158 // |x| < 2.4375
159 if (ix < 0x401c_0000) {
160 break :blk .{ (x_ - 1.5) / (1.0 + 1.5 * x_), 2 };
161 }
162 // 2.4375 <= |x| < 2^26
163 else {
164 break :blk .{ -1.0 / x_, 3 };
165 }
100166 }
101167 }
102 }
103
104 const z = x * x;
168 };
169 // end of argument reduction
170 const z = x_ * x_;
105171 const w = z * z;
172 // break sum from i=0 to 10 aT[i]z^(i+1) into odd and even poly
106173 const s1 = z * (aT[0] + w * (aT[2] + w * aT[4]));
107174 const s2 = w * (aT[1] + w * aT[3]);
108
109 if (id) |id_value| {
110 const zz = atanhi[id_value] - ((x * (s1 + s2) - atanlo[id_value]) - x);
111 return if (sign != 0) -zz else zz;
175 if (id) |id_| {
176 const z_ = atanhi[id_] - ((x_ * (s1 + s2) - atanlo[id_]) - x_);
177 return if (sign) -z_ else z_;
112178 } else {
113 return x - x * (s1 + s2);
179 return x_ - x_ * (s1 + s2);
114180 }
115181}
116182
117fn atan64(x_: f64) f64 {
118 const atanhi = [_]f64{
119 4.63647609000806093515e-01, // atan(0.5)hi
120 7.85398163397448278999e-01, // atan(1.0)hi
121 9.82793723247329054082e-01, // atan(1.5)hi
122 1.57079632679489655800e+00, // atan(inf)hi
183fn atanBinary64(x: f64) f64 {
184 const atanhi: []const f64 = &.{
185 4.63647609000806093515e-01, // atan(0.5)hi 0x3FDDAC67, 0x0561BB4F
186 7.85398163397448278999e-01, // atan(1.0)hi 0x3FE921FB, 0x54442D18
187 9.82793723247329054082e-01, // atan(1.5)hi 0x3FEF730B, 0xD281F69B
188 1.57079632679489655800e+00, // atan(inf)hi 0x3FF921FB, 0x54442D18
123189 };
124
125 const atanlo = [_]f64{
126 2.26987774529616870924e-17, // atan(0.5)lo
127 3.06161699786838301793e-17, // atan(1.0)lo
128 1.39033110312309984516e-17, // atan(1.5)lo
129 6.12323399573676603587e-17, // atan(inf)lo
190 const atanlo: []const f64 = &.{
191 2.26987774529616870924e-17, // atan(0.5)lo 0x3C7A2B7F, 0x222F65E2
192 3.06161699786838301793e-17, // atan(1.0)lo 0x3C81A626, 0x33145C07
193 1.39033110312309984516e-17, // atan(1.5)lo 0x3C700788, 0x7AF0CBBD
194 6.12323399573676603587e-17, // atan(inf)lo 0x3C91A626, 0x33145C07
130195 };
131
132 const aT = [_]f64{
133 3.33333333333329318027e-01,
134 -1.99999999998764832476e-01,
135 1.42857142725034663711e-01,
136 -1.11111104054623557880e-01,
137 9.09088713343650656196e-02,
138 -7.69187620504482999495e-02,
139 6.66107313738753120669e-02,
140 -5.83357013379057348645e-02,
141 4.97687799461593236017e-02,
142 -3.65315727442169155270e-02,
143 1.62858201153657823623e-02,
196 const aT: []const f64 = &.{
197 3.33333333333329318027e-01, // 0x3FD55555, 0x5555550D
198 -1.99999999998764832476e-01, // 0xBFC99999, 0x9998EBC4
199 1.42857142725034663711e-01, // 0x3FC24924, 0x920083FF
200 -1.11111104054623557880e-01, // 0xBFBC71C6, 0xFE231671
201 9.09088713343650656196e-02, // 0x3FB745CD, 0xC54C206E
202 -7.69187620504482999495e-02, // 0xBFB3B0F2, 0xAF749A6D
203 6.66107313738753120669e-02, // 0x3FB10D66, 0xA0D03D51
204 -5.83357013379057348645e-02, // 0xBFADDE2D, 0x52DEFD9A
205 4.97687799461593236017e-02, // 0x3FA97B4B, 0x24760DEB
206 -3.65315727442169155270e-02, // 0xBFA2B444, 0x2C6A6C2F
207 1.62858201153657823623e-02, // 0x3F90AD3A, 0xE322DA11
144208 };
145209
146 var x = x_;
147 const ux: u64 = @bitCast(x);
148 var ix: u32 = @intCast(ux >> 32);
149 const sign = ix >> 31;
150 ix &= 0x7FFFFFFF;
151
152 // |x| >= 2^66
210 const hx: u64 = @bitCast(x);
211 const ix: u32 = @truncate((hx >> 32) & 0x7fffffff);
212 const sign = (hx >> 63) != 0;
213 // if |x| >= 2^66
153214 if (ix >= 0x44100000) {
154215 if (math.isNan(x)) {
155216 return x;
217 }
218 const z = atanhi[3] + 0x1p-120;
219 return if (sign) -z else z;
220 }
221 const x_, const id: ?usize = blk: {
222 // |x| < 0.4375
223 if (ix < 0x3fdc_0000) {
224 // |x| < 2^(-27)
225 if (ix < 0x3e40_0000) {
226 if (ix < 0x0010_0000) {
227 // raise underflow for subnormal x
228 mem.doNotOptimizeAway(@as(f32, @floatCast(x)));
229 }
230 return x;
231 }
232 break :blk .{ x, null };
156233 } else {
157 const z = atanhi[3] + 0x1.0p-120;
158 return if (sign != 0) -z else z;
234 const x_ = @abs(x);
235 // |x| < 1.1875
236 if (ix < 0x3ff3_0000) {
237 // 7/16 <= |x| < 11/16
238 if (ix < 0x3fe6_0000) {
239 break :blk .{ (2.0 * x_ - 1.0) / (2.0 + x_), 0 };
240 }
241 // 11/16 <= |x| < 19/16
242 else {
243 break :blk .{ (x_ - 1.0) / (x_ + 1.0), 1 };
244 }
245 } else {
246 // |x| < 2.4375
247 if (ix < 0x4003_8000) {
248 break :blk .{ (x_ - 1.5) / (1.0 + 1.5 * x_), 2 };
249 }
250 // 2.4375 <= |x| < 2^66
251 else {
252 break :blk .{ -1.0 / x_, 3 };
253 }
254 }
159255 }
256 };
257 // end of argument reduction
258 const z = x_ * x_;
259 const w = z * z;
260 // break sum from i=0 to 10 aT[i]z^(i+1) into odd and even poly
261 const s1 = z * (aT[0] + w * (aT[2] + w * (aT[4] + w * (aT[6] + w * (aT[8] + w * aT[10])))));
262 const s2 = w * (aT[1] + w * (aT[3] + w * (aT[5] + w * (aT[7] + w * aT[9]))));
263 if (id) |id_| {
264 const z_ = atanhi[id_] - (x_ * (s1 + s2) - atanlo[id_] - x_);
265 return if (sign) -z_ else z_;
266 } else {
267 return x_ - x_ * (s1 + s2);
160268 }
269}
161270
162 var id: ?usize = undefined;
271fn atanExtended80(x: f80) f80 {
272 const atanhi: []const f80 = &.{
273 4.63647609000806116202e-01,
274 7.85398163397448309628e-01,
275 9.82793723247329067960e-01,
276 1.57079632679489661926e+00,
277 };
278 const atanlo: []const f80 = &.{
279 1.18469937025062860669e-20,
280 -1.25413940316708300586e-20,
281 2.55232234165405176172e-20,
282 -2.50827880633416601173e-20,
283 };
284 const aT: []const f80 = &.{
285 3.33333333333333333017e-01,
286 -1.99999999999999632011e-01,
287 1.42857142857046531280e-01,
288 -1.11111111100562372733e-01,
289 9.09090902935647302252e-02,
290 -7.69230552476207730353e-02,
291 6.66661718042406260546e-02,
292 -5.88158892835030888692e-02,
293 5.25499891539726639379e-02,
294 -4.70119845393155721494e-02,
295 4.03539201366454414072e-02,
296 -2.91303858419364158725e-02,
297 1.24822046299269234080e-02,
298 };
163299
164 // |x| < 0.4375
165 if (ix < 0x3FDC0000) {
166 // |x| < 2^(-27)
167 if (ix < 0x3E400000) {
168 if (ix < 0x00100000) {
169 mem.doNotOptimizeAway(@as(f32, @floatCast(x)));
170 }
300 const hx: u80 = @bitCast(x);
301 const se: u16 = @truncate(hx >> 64);
302 const e = se & 0x7fff;
303 const sign = se >> 15 != 0;
304 // if |x| is large, atan(x)~=pi/2
305 if (e >= 0x3fff + math.floatMantissaBits(f80) + 1) {
306 if (math.isNan(x)) {
171307 return x;
172308 }
173 id = null;
174 } else {
175 x = @abs(x);
176 // |x| < 1.1875
177 if (ix < 0x3FF30000) {
178 // 7/16 <= |x| < 11/16
179 if (ix < 0x3FE60000) {
180 id = 0;
181 x = (2.0 * x - 1.0) / (2.0 + x);
182 }
183 // 11/16 <= |x| < 19/16
184 else {
185 id = 1;
186 x = (x - 1.0) / (x + 1.0);
309 return if (sign) -atanhi[3] else atanhi[3];
310 }
311 // Extract the exponent and the first few bits of the mantissa.
312 const m: u64 = @truncate(hx & 0x0000_ffff_ffff_ffff_ffff);
313 const expman = ((@as(u32, @intCast(se)) & 0x7fff) << 8) | (@as(u32, @truncate(m >> 55)) & 0xff);
314 const x_, const id: ?usize = blk: {
315 // |x| < 0.4375
316 if (expman < ((0x3fff - 2) << 8) + 0xc0) {
317 // if |x| is small, atanl(x)~=x
318 if (e < 0x3fff - (math.floatMantissaBits(f80) + 1) / 2) {
319 // raise underflow if subnormal
320 if (e == 0) {
321 std.mem.doNotOptimizeAway(@as(f32, @floatCast(x)));
322 }
323 return x;
187324 }
325 break :blk .{ x, null };
188326 } else {
189 // |x| < 2.4375
190 if (ix < 0x40038000) {
191 id = 2;
192 x = (x - 1.5) / (1.0 + 1.5 * x);
193 }
194 // 2.4375 <= |x| < 2^66
195 else {
196 id = 3;
197 x = -1.0 / x;
327 const x_ = @abs(x);
328 // |x| < 1.1875
329 if (expman < (0x3fff << 8) + 0x30) {
330 // 7/16 <= |x| < 11/16
331 if (expman < ((0x3fff - 1) << 8) + 0x60) {
332 break :blk .{ (2.0 * x_ - 1.0) / (2.0 + x_), 0 };
333 }
334 // 11/16 <= |x| < 19/16
335 else {
336 break :blk .{ (x_ - 1.0) / (x_ + 1.0), 1 };
337 }
338 } else {
339 // |x| < 2.4375
340 if (expman < ((0x3fff + 1) << 8) + 0x38) {
341 break :blk .{ (x_ - 1.5) / (1.0 + 1.5 * x_), 2 };
342 }
343 // 2.4375 <= |x|
344 else {
345 break :blk .{ -1.0 / x_, 3 };
346 }
198347 }
199348 }
349 };
350 // end of argument reduction
351 const z = x_ * x_;
352 const w = z * z;
353 // break sum aT[i]z^(i+1) into odd and even poly
354 const s1 = z * (aT[0] + w * (aT[2] + w * (aT[4] + w * (aT[6] + w * (aT[8] + w * (aT[10] + w * aT[12]))))));
355 const s2 = w * (aT[1] + w * (aT[3] + w * (aT[5] + w * (aT[7] + w * (aT[9] + w * aT[11])))));
356 if (id) |id_| {
357 const z_ = atanhi[id_] - ((x_ * (s1 + s2) - atanlo[id_]) - x_);
358 return if (sign) -z_ else z_;
359 } else {
360 return x_ - x_ * (s1 + s2);
200361 }
362}
201363
202 const z = x * x;
203 const w = z * z;
204 const s1 = z * (aT[0] + w * (aT[2] + w * (aT[4] + w * (aT[6] + w * (aT[8] + w * aT[10])))));
205 const s2 = w * (aT[1] + w * (aT[3] + w * (aT[5] + w * (aT[7] + w * aT[9]))));
364fn atanBinary128(x: f128) f128 {
365 const atanhi: []const f128 = &.{
366 4.63647609000806116214256231461214397e-01,
367 7.85398163397448309615660845819875699e-01,
368 9.82793723247329067985710611014666038e-01,
369 1.57079632679489661923132169163975140e+00,
370 };
371 const atanlo: []const f128 = &.{
372 4.89509642257333492668618435220297706e-36,
373 2.16795253253094525619926100651083806e-35,
374 -2.31288434538183565909319952098066272e-35,
375 4.33590506506189051239852201302167613e-35,
376 };
377 const aT: []const f128 = &.{
378 3.33333333333333333333333333333333125e-01,
379 -1.99999999999999999999999999999180430e-01,
380 1.42857142857142857142857142125269827e-01,
381 -1.11111111111111111111110834490810169e-01,
382 9.09090909090909090908522355708623681e-02,
383 -7.69230769230769230696553844935357021e-02,
384 6.66666666666666660390096773046256096e-02,
385 -5.88235294117646671706582985209643694e-02,
386 5.26315789473666478515847092020327506e-02,
387 -4.76190476189855517021024424991436144e-02,
388 4.34782608678695085948531993458097026e-02,
389 -3.99999999632663469330634215991142368e-02,
390 3.70370363987423702891250829918659723e-02,
391 -3.44827496515048090726669907612335954e-02,
392 3.22579620681420149871973710852268528e-02,
393 -3.03020767654269261041647570626778067e-02,
394 2.85641979882534783223403715930946138e-02,
395 -2.69824879726738568189929461383741323e-02,
396 2.54194698498808542954187110873675769e-02,
397 -2.35083879708189059926183138130183215e-02,
398 2.04832358998165364349957325067131428e-02,
399 -1.54489555488544397858507248612362957e-02,
400 8.64492360989278761493037861575248038e-03,
401 -2.58521121597609872727919154569765469e-03,
402 };
206403
207 if (id) |id_value| {
208 const zz = atanhi[id_value] - ((x * (s1 + s2) - atanlo[id_value]) - x);
209 return if (sign != 0) -zz else zz;
404 const hx: u128 = @bitCast(x);
405 const se: u16 = @truncate(hx >> 112);
406 const e = se & 0x7fff;
407 const sign = se >> 15 != 0;
408 // if |x| is large, atan(x)~=pi/2
409 if (e >= 0x3fff + math.floatMantissaBits(f128) + 2) {
410 if (math.isNan(x)) {
411 return x;
412 }
413 return if (sign) -atanhi[3] else atanhi[3];
414 }
415 // Extract the exponent and the first few bits of the mantissa.
416 const top: u16 = @truncate((hx >> 96) & 0x0000_ffff);
417 const expman = ((@as(u32, @intCast(se)) & 0x7fff) << 8) | (@as(u32, @intCast(top)) >> 8);
418 const x_, const id: ?usize = blk: {
419 // |x| < 0.4375
420 if (expman < ((0x3fff - 2) << 8) + 0xc0) {
421 // if |x| is small, atanl(x)~=x
422 if (e < 0x3fff - (math.floatMantissaBits(f128) + 2) / 2) {
423 // raise underflow if subnormal
424 if (e == 0) {
425 mem.doNotOptimizeAway(@as(f32, @floatCast(x)));
426 }
427 return x;
428 }
429 break :blk .{ x, null };
430 } else {
431 const x_ = @abs(x);
432 // |x| < 1.1875
433 if (expman < (0x3fff << 8) + 0x30) {
434 // 7/16 <= |x| < 11/16
435 if (expman < ((0x3fff - 1) << 8) + 0x60) {
436 break :blk .{ (2.0 * x_ - 1.0) / (2.0 + x_), 0 };
437 }
438 // 11/16 <= |x| < 19/16
439 else {
440 break :blk .{ (x_ - 1.0) / (x_ + 1.0), 1 };
441 }
442 } else {
443 // |x| < 2.4375
444 if (expman < ((0x3fff + 1) << 8) + 0x38) {
445 break :blk .{ (x_ - 1.5) / (1.0 + 1.5 * x_), 2 };
446 }
447 // 2.4375 <= |x|
448 else {
449 break :blk .{ -1.0 / x_, 3 };
450 }
451 }
452 }
453 };
454 // end of argument reduction
455 const z = x_ * x_;
456 const w = z * z;
457 // break sum aT[i]z^(i+1) into odd and even poly
458 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])))))))))));
459 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])))))))))));
460 if (id) |id_| {
461 const z_ = atanhi[id_] - ((x_ * (s1 + s2) - atanlo[id_]) - x_);
462 return if (sign) -z_ else z_;
210463 } else {
211 return x - x * (s1 + s2);
464 return x_ - x_ * (s1 + s2);
212465 }
213466}
214467
215test atan {
216 try expect(@as(u32, @bitCast(atan(@as(f32, 0.2)))) == @as(u32, @bitCast(atan32(0.2))));
217 try expect(atan(@as(f64, 0.2)) == atan64(0.2));
468test "atanBinary16.special" {
469 try testing.expectEqual(atanBinary16(0x0p+0), 0x0p+0);
470 try testing.expectEqual(atanBinary16(-0x0p+0), -0x0p+0);
471 try testing.expectApproxEqAbs(atanBinary16(0x1p+0), 0x1.92p-1, math.floatEpsAt(f16, 0x1.92p-1));
472 try testing.expectApproxEqAbs(atanBinary16(-0x1p+0), -0x1.92p-1, math.floatEpsAt(f16, -0x1.92p-1));
473 try testing.expectApproxEqAbs(atanBinary16(math.inf(f16)), 0x1.92p0, math.floatEpsAt(f16, 0x1.92p0));
474 try testing.expectApproxEqAbs(atanBinary16(-math.inf(f16)), -0x1.92p0, math.floatEpsAt(f16, -0x1.92p0));
475 try testing.expect(math.isNan(atanBinary16(math.nan(f16))));
218476}
219477
220test atan32 {
221 const epsilon = 0.000001;
478test "atanBinary16" {
479 try testing.expectApproxEqAbs(atanBinary16(-0x1.864p-2), -0x1.74cp-2, math.floatEpsAt(f16, -0x1.74cp-2));
480 try testing.expectApproxEqAbs(atanBinary16(-0x1.59cp1), -0x1.374p0, math.floatEpsAt(f16, -0x1.374p0));
481 try testing.expectApproxEqAbs(atanBinary16(-0x1.d2cp0), -0x1.11cp0, math.floatEpsAt(f16, -0x1.11cp0));
482 try testing.expectApproxEqAbs(atanBinary16(-0x1.5f4p-1), -0x1.33cp-1, math.floatEpsAt(f16, -0x1.33cp-1));
483 try testing.expectApproxEqAbs(atanBinary16(0x1.588p1), 0x1.37p0, math.floatEpsAt(f16, 0x1.37p0));
484 try testing.expectApproxEqAbs(atanBinary16(-0x1.b14p-2), -0x1.99cp-2, math.floatEpsAt(f16, -0x1.99cp-2));
485 try testing.expectApproxEqAbs(atanBinary16(0x1.3ccp1), 0x1.2fcp0, math.floatEpsAt(f16, 0x1.2fcp0));
486 try testing.expectApproxEqAbs(atanBinary16(-0x1.0ecp-2), -0x1.08cp-2, math.floatEpsAt(f16, -0x1.08cp-2));
487 try testing.expectApproxEqAbs(atanBinary16(0x1.298p1), 0x1.2ap0, math.floatEpsAt(f16, 0x1.2ap0));
488 try testing.expectApproxEqAbs(atanBinary16(-0x1.028p1), -0x1.1c8p0, math.floatEpsAt(f16, -0x1.1c8p0));
489}
222490
223 try expect(math.approxEqAbs(f32, atan32(0.2), 0.197396, epsilon));
224 try expect(math.approxEqAbs(f32, atan32(-0.2), -0.197396, epsilon));
225 try expect(math.approxEqAbs(f32, atan32(0.3434), 0.330783, epsilon));
226 try expect(math.approxEqAbs(f32, atan32(0.8923), 0.728545, epsilon));
227 try expect(math.approxEqAbs(f32, atan32(1.5), 0.982794, epsilon));
491test "atanBinary32.special" {
492 try testing.expectEqual(atanBinary32(0x0p+0), 0x0p+0);
493 try testing.expectEqual(atanBinary32(-0x0p+0), -0x0p+0);
494 try testing.expectApproxEqAbs(atanBinary32(0x1p+0), 0x1.921fb6p-1, math.floatEpsAt(f32, 0x1.921fb6p-1));
495 try testing.expectApproxEqAbs(atanBinary32(-0x1p+0), -0x1.921fb6p-1, math.floatEpsAt(f32, -0x1.921fb6p-1));
496 try testing.expectApproxEqAbs(atanBinary32(math.inf(f32)), 0x1.921fb6p+0, math.floatEpsAt(f32, 0x1.921fb6p+0));
497 try testing.expectApproxEqAbs(atanBinary32(-math.inf(f32)), -0x1.921fb6p+0, math.floatEpsAt(f32, -0x1.921fb6p+0));
498 try testing.expect(math.isNan(atanBinary32(math.nan(f32))));
228499}
229500
230test atan64 {
231 const epsilon = 0.000001;
501test "atanBinary32" {
502 try testing.expectApproxEqAbs(atanBinary32(-0x1.8629dp-2), -0x1.74c62p-2, math.floatEpsAt(f32, -0x1.74c62p-2));
503 try testing.expectApproxEqAbs(atanBinary32(-0x1.59d42ep1), -0x1.375fd8p0, math.floatEpsAt(f32, -0x1.375fd8p0));
504 try testing.expectApproxEqAbs(atanBinary32(-0x1.d2dbe2p0), -0x1.11b8aep0, math.floatEpsAt(f32, -0x1.11b8aep0));
505 try testing.expectApproxEqAbs(atanBinary32(-0x1.5f314ep-1), -0x1.33d28cp-1, math.floatEpsAt(f32, -0x1.33d28cp-1));
506 try testing.expectApproxEqAbs(atanBinary32(0x1.5869bp1), 0x1.37082ep0, math.floatEpsAt(f32, 0x1.37082ep0));
507 try testing.expectApproxEqAbs(atanBinary32(-0x1.b13a06p-2), -0x1.99d7cap-2, math.floatEpsAt(f32, -0x1.99d7cap-2));
508 try testing.expectApproxEqAbs(atanBinary32(0x1.3cb0f2p1), 0x1.2fcb12p0, math.floatEpsAt(f32, 0x1.2fcb12p0));
509 try testing.expectApproxEqAbs(atanBinary32(-0x1.0ed746p-2), -0x1.08c71ap-2, math.floatEpsAt(f32, -0x1.08c71ap-2));
510 try testing.expectApproxEqAbs(atanBinary32(0x1.299d54p1), 0x1.2a24e2p0, math.floatEpsAt(f32, 0x1.2a24e2p0));
511 try testing.expectApproxEqAbs(atanBinary32(-0x1.0264fcp1), -0x1.1c6178p0, math.floatEpsAt(f32, -0x1.1c6178p0));
512}
232513
233 try expect(math.approxEqAbs(f64, atan64(0.2), 0.197396, epsilon));
234 try expect(math.approxEqAbs(f64, atan64(-0.2), -0.197396, epsilon));
235 try expect(math.approxEqAbs(f64, atan64(0.3434), 0.330783, epsilon));
236 try expect(math.approxEqAbs(f64, atan64(0.8923), 0.728545, epsilon));
237 try expect(math.approxEqAbs(f64, atan64(1.5), 0.982794, epsilon));
514test "atanBinary64.special" {
515 try testing.expectEqual(atanBinary64(0x0p+0), 0x0p+0);
516 try testing.expectEqual(atanBinary64(-0x0p+0), -0x0p+0);
517 try testing.expectApproxEqAbs(atanBinary64(0x1p+0), 0x1.921fb54442d18p-1, math.floatEpsAt(f64, 0x1.921fb54442d18p-1));
518 try testing.expectApproxEqAbs(atanBinary64(-0x1p+0), -0x1.921fb54442d18p-1, math.floatEpsAt(f64, -0x1.921fb54442d18p-1));
519 try testing.expectApproxEqAbs(atanBinary64(math.inf(f64)), 0x1.921fb54442d18p+0, math.floatEpsAt(f64, 0x1.921fb54442d18p+0));
520 try testing.expectApproxEqAbs(atanBinary64(-math.inf(f64)), -0x1.921fb54442d18p+0, math.floatEpsAt(f64, -0x1.921fb54442d18p+0));
521 try testing.expect(math.isNan(atanBinary64(math.nan(f64))));
238522}
239523
240test "atan32.special" {
241 const epsilon = 0.000001;
524test "atanBinary64" {
525 try testing.expectApproxEqAbs(atanBinary64(-0x1.8629d0244cdccp-2), -0x1.74c61f4377016p-2, math.floatEpsAt(f64, -0x1.74c61f4377016p-2));
526 try testing.expectApproxEqAbs(atanBinary64(-0x1.59d42d4659937p1), -0x1.375fd7987cc2p0, math.floatEpsAt(f64, -0x1.375fd7987cc2p0));
527 try testing.expectApproxEqAbs(atanBinary64(-0x1.d2dbe23d04f06p0), -0x1.11b8adeba5616p0, math.floatEpsAt(f64, -0x1.11b8adeba5616p0));
528 try testing.expectApproxEqAbs(atanBinary64(-0x1.5f314e72398e8p-1), -0x1.33d28ca762539p-1, math.floatEpsAt(f64, -0x1.33d28ca762539p-1));
529 try testing.expectApproxEqAbs(atanBinary64(0x1.5869af37b7d08p1), 0x1.37082ce2dd03p0, math.floatEpsAt(f64, 0x1.37082ce2dd03p0));
530 try testing.expectApproxEqAbs(atanBinary64(-0x1.b13a05a662618p-2), -0x1.99d7cac66dd44p-2, math.floatEpsAt(f64, -0x1.99d7cac66dd44p-2));
531 try testing.expectApproxEqAbs(atanBinary64(0x1.3cb0f12f39d8ap1), 0x1.2fcb120468e8ep0, math.floatEpsAt(f64, 0x1.2fcb120468e8ep0));
532 try testing.expectApproxEqAbs(atanBinary64(-0x1.0ed746b39cbb7p-2), -0x1.08c71aa0e509p-2, math.floatEpsAt(f64, -0x1.08c71aa0e509p-2));
533 try testing.expectApproxEqAbs(atanBinary64(0x1.299d54ac7d6bp1), 0x1.2a24e22d861dfp0, math.floatEpsAt(f64, 0x1.2a24e22d861dfp0));
534 try testing.expectApproxEqAbs(atanBinary64(-0x1.0264fb9f3d50ep1), -0x1.1c617825f9751p0, math.floatEpsAt(f64, -0x1.1c617825f9751p0));
535}
536
537test "atanExtended80.special" {
538 try testing.expectEqual(atanExtended80(0x0p+0), 0x0p+0);
539 try testing.expectEqual(atanExtended80(-0x0p+0), -0x0p+0);
540 try testing.expectApproxEqAbs(atanExtended80(0x1p+0), 0x1.921fb54442d1846ap-1, math.floatEpsAt(f80, 0x1.921fb54442d1846ap-1));
541 try testing.expectApproxEqAbs(atanExtended80(-0x1p+0), -0x1.921fb54442d1846ap-1, math.floatEpsAt(f80, -0x1.921fb54442d1846ap-1));
542 try testing.expectApproxEqAbs(atanExtended80(math.inf(f80)), 0x1.921fb54442d1846ap0, math.floatEpsAt(f80, 0x1.921fb54442d1846ap0));
543 try testing.expectApproxEqAbs(atanExtended80(-math.inf(f80)), -0x1.921fb54442d1846ap0, math.floatEpsAt(f80, -0x1.921fb54442d1846ap0));
544 try testing.expect(math.isNan(atanExtended80(math.nan(f80))));
545}
242546
243 try expect(math.isPositiveZero(atan32(0.0)));
244 try expect(math.isNegativeZero(atan32(-0.0)));
245 try expect(math.approxEqAbs(f32, atan32(math.inf(f32)), math.pi / 2.0, epsilon));
246 try expect(math.approxEqAbs(f32, atan32(-math.inf(f32)), -math.pi / 2.0, epsilon));
547test "atanExtended80" {
548 try testing.expectApproxEqAbs(atanExtended80(-0x1.8629d0244cdcbed8p-2), -0x1.74c61f437701661p-2, math.floatEpsAt(f80, -0x1.74c61f437701661p-2));
549 try testing.expectApproxEqAbs(atanExtended80(-0x1.59d42d4659936d9ep1), -0x1.375fd7987cc1fd02p0, math.floatEpsAt(f80, -0x1.375fd7987cc1fd02p0));
550 try testing.expectApproxEqAbs(atanExtended80(-0x1.d2dbe23d04f067b4p0), -0x1.11b8adeba5615e04p0, math.floatEpsAt(f80, -0x1.11b8adeba5615e04p0));
551 try testing.expectApproxEqAbs(atanExtended80(-0x1.5f314e72398e7dbcp-1), -0x1.33d28ca76253964cp-1, math.floatEpsAt(f80, -0x1.33d28ca76253964cp-1));
552 try testing.expectApproxEqAbs(atanExtended80(0x1.5869af37b7d078cap1), 0x1.37082ce2dd03010cp0, math.floatEpsAt(f80, 0x1.37082ce2dd03010cp0));
553 try testing.expectApproxEqAbs(atanExtended80(-0x1.b13a05a66261821ap-2), -0x1.99d7cac66dd4438p-2, math.floatEpsAt(f80, -0x1.99d7cac66dd4438p-2));
554 try testing.expectApproxEqAbs(atanExtended80(0x1.3cb0f12f39d899cp1), 0x1.2fcb120468e8d9ecp0, math.floatEpsAt(f80, 0x1.2fcb120468e8d9ecp0));
555 try testing.expectApproxEqAbs(atanExtended80(-0x1.0ed746b39cbb7614p-2), -0x1.08c71aa0e5090998p-2, math.floatEpsAt(f80, -0x1.08c71aa0e5090998p-2));
556 try testing.expectApproxEqAbs(atanExtended80(0x1.299d54ac7d6afc52p1), 0x1.2a24e22d861debfep0, math.floatEpsAt(f80, 0x1.2a24e22d861debfep0));
557 try testing.expectApproxEqAbs(atanExtended80(-0x1.0264fb9f3d50e4fp1), -0x1.1c617825f97512b8p0, math.floatEpsAt(f80, -0x1.1c617825f97512b8p0));
247558}
248559
249test "atan64.special" {
250 const epsilon = 0.000001;
560test "atanBinary128.special" {
561 try testing.expectEqual(atanBinary128(0x0p+0), 0x0p+0);
562 try testing.expectEqual(atanBinary128(-0x0p+0), -0x0p+0);
563 try testing.expectApproxEqAbs(atanBinary128(0x1p+0), 0x1.921fb54442d18469898cc51701b8p-1, math.floatEpsAt(f128, 0x1.921fb54442d18469898cc51701b8p-1));
564 try testing.expectApproxEqAbs(atanBinary128(-0x1p+0), -0x1.921fb54442d18469898cc51701b8p-1, math.floatEpsAt(f128, -0x1.921fb54442d18469898cc51701b8p-1));
565 try testing.expectApproxEqAbs(atanBinary128(math.inf(f128)), 0x1.921fb54442d18469898cc51701b8p0, math.floatEpsAt(f128, 0x1.921fb54442d18469898cc51701b8p0));
566 try testing.expectApproxEqAbs(atanBinary128(-math.inf(f128)), -0x1.921fb54442d18469898cc51701b8p0, math.floatEpsAt(f128, -0x1.921fb54442d18469898cc51701b8p0));
567 try testing.expect(math.isNan(atanBinary128(math.nan(f128))));
568}
251569
252 try expect(math.isPositiveZero(atan64(0.0)));
253 try expect(math.isNegativeZero(atan64(-0.0)));
254 try expect(math.approxEqAbs(f64, atan64(math.inf(f64)), math.pi / 2.0, epsilon));
255 try expect(math.approxEqAbs(f64, atan64(-math.inf(f64)), -math.pi / 2.0, epsilon));
570test "atanBinary128" {
571 try testing.expectApproxEqAbs(atanBinary128(-0x1.8629d0244cdcbed71792ccdec26dp-2), -0x1.74c61f437701660ff76989d23707p-2, math.floatEpsAt(f128, -0x1.74c61f437701660ff76989d23707p-2));
572 try testing.expectApproxEqAbs(atanBinary128(-0x1.59d42d4659936d9e22b5dea4faefp1), -0x1.375fd7987cc1fd0119cf0cc5b708p0, math.floatEpsAt(f128, -0x1.375fd7987cc1fd0119cf0cc5b708p0));
573 try testing.expectApproxEqAbs(atanBinary128(-0x1.d2dbe23d04f067b42da3f8efdf57p0), -0x1.11b8adeba5615e0370722b511231p0, math.floatEpsAt(f128, -0x1.11b8adeba5615e0370722b511231p0));
574 try testing.expectApproxEqAbs(atanBinary128(-0x1.5f314e72398e7dbbe70fb072983ep-1), -0x1.33d28ca76253964cb5d3581cdd88p-1, math.floatEpsAt(f128, -0x1.33d28ca76253964cb5d3581cdd88p-1));
575 try testing.expectApproxEqAbs(atanBinary128(0x1.5869af37b7d078caa3456c44aecep1), 0x1.37082ce2dd03010bbea814dc5882p0, math.floatEpsAt(f128, 0x1.37082ce2dd03010bbea814dc5882p0));
576 try testing.expectApproxEqAbs(atanBinary128(-0x1.b13a05a66261821a364ad8c6c999p-2), -0x1.99d7cac66dd4438077284b491a91p-2, math.floatEpsAt(f128, -0x1.99d7cac66dd4438077284b491a91p-2));
577 try testing.expectApproxEqAbs(atanBinary128(0x1.3cb0f12f39d899c0d963ac413297p1), 0x1.2fcb120468e8d9ebdb74702314c8p0, math.floatEpsAt(f128, 0x1.2fcb120468e8d9ebdb74702314c8p0));
578 try testing.expectApproxEqAbs(atanBinary128(-0x1.0ed746b39cbb7614d8735e8315a8p-2), -0x1.08c71aa0e5090998206fbbe2090fp-2, math.floatEpsAt(f128, -0x1.08c71aa0e5090998206fbbe2090fp-2));
579 try testing.expectApproxEqAbs(atanBinary128(0x1.299d54ac7d6afc5154643b601519p1), 0x1.2a24e22d861debfd6f974500567fp0, math.floatEpsAt(f128, 0x1.2a24e22d861debfd6f974500567fp0));
580 try testing.expectApproxEqAbs(atanBinary128(-0x1.0264fb9f3d50e4f0f966f0686064p1), -0x1.1c617825f97512b7f38656ab12cdp0, math.floatEpsAt(f128, -0x1.1c617825f97512b7f38656ab12cdp0));
256581}
src/libs/musl.zig-8
......@@ -836,12 +836,9 @@ const src_files = [_][]const u8{
836836 "musl/src/math/atan2.c",
837837 "musl/src/math/atan2f.c",
838838 "musl/src/math/atan2l.c",
839 "musl/src/math/atan.c",
840 "musl/src/math/atanf.c",
841839 "musl/src/math/atanh.c",
842840 "musl/src/math/atanhf.c",
843841 "musl/src/math/atanhl.c",
844 "musl/src/math/atanl.c",
845842 "musl/src/math/cbrt.c",
846843 "musl/src/math/cbrtf.c",
847844 "musl/src/math/cbrtl.c",
......@@ -892,9 +889,6 @@ const src_files = [_][]const u8{
892889 "musl/src/math/i386/atan2f.s",
893890 "musl/src/math/i386/atan2l.s",
894891 "musl/src/math/i386/atan2.s",
895 "musl/src/math/i386/atanf.s",
896 "musl/src/math/i386/atanl.s",
897 "musl/src/math/i386/atan.s",
898892 "musl/src/math/i386/exp2l.s",
899893 "musl/src/math/i386/exp_ld.s",
900894 "musl/src/math/i386/expl.s",
......@@ -1076,7 +1070,6 @@ const src_files = [_][]const u8{
10761070 "musl/src/math/x32/acosl.s",
10771071 "musl/src/math/x32/asinl.s",
10781072 "musl/src/math/x32/atan2l.s",
1079 "musl/src/math/x32/atanl.s",
10801073 "musl/src/math/x32/exp2l.s",
10811074 "musl/src/math/x32/expl.s",
10821075 "musl/src/math/x32/expm1l.s",
......@@ -1098,7 +1091,6 @@ const src_files = [_][]const u8{
10981091 "musl/src/math/x86_64/acosl.s",
10991092 "musl/src/math/x86_64/asinl.s",
11001093 "musl/src/math/x86_64/atan2l.s",
1101 "musl/src/math/x86_64/atanl.s",
11021094 "musl/src/math/x86_64/exp2l.s",
11031095 "musl/src/math/x86_64/expl.s",
11041096 "musl/src/math/x86_64/expm1l.s",
src/libs/wasi_libc.zig-3
......@@ -698,12 +698,9 @@ const libc_top_half_src_files = [_][]const u8{
698698 "musl/src/math/atan2.c",
699699 "musl/src/math/atan2f.c",
700700 "musl/src/math/atan2l.c",
701 "musl/src/math/atan.c",
702 "musl/src/math/atanf.c",
703701 "musl/src/math/atanh.c",
704702 "musl/src/math/atanhf.c",
705703 "musl/src/math/atanhl.c",
706 "musl/src/math/atanl.c",
707704 "musl/src/math/cbrt.c",
708705 "musl/src/math/cbrtf.c",
709706 "musl/src/math/cbrtl.c",