authorgravatar for andrew@ziglang.orgAndrew Kelley <andrew@ziglang.org> 2021-12-05 14:49:04-08:00
committergravatar for noreply@github.comGitHub <noreply@github.com> 2021-12-05 14:49:04-08:00
loga7828c261a0909a1ed672761adf332c978d492c7
tree101e26c5c26cdef3f689a3866ed8deaca6597336
parent9e03cf948948ea25feb3f413f234e98e71d55786
parentb02384e03d1e52d72ccf7dd7d3cec30f05b51082
signaturebadge-question-mark Signed by PGP key 4AEE18F83AFDEB23

Merge pull request #10276 from tiehuis/musl-trig-reimpl

sin/cos/tan musl reimplementation

7 files changed, 1376 insertions(+), 232 deletions(-)

lib/std/math/__rem_pio2.zig created+198
...@@ -0,0 +1,198 @@
1// Ported from musl, which is licensed under the MIT license:
2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
3//
4// https://git.musl-libc.org/cgit/musl/tree/src/math/__rem_pio2.c
5
6const std = @import("../std.zig");
7const __rem_pio2_large = @import("__rem_pio2_large.zig").__rem_pio2_large;
8const math = std.math;
9
10const toint = 1.5 / math.epsilon(f64);
11// pi/4
12const pio4 = 0x1.921fb54442d18p-1;
13// invpio2: 53 bits of 2/pi
14const invpio2 = 6.36619772367581382433e-01; // 0x3FE45F30, 0x6DC9C883
15// pio2_1: first 33 bit of pi/2
16const pio2_1 = 1.57079632673412561417e+00; // 0x3FF921FB, 0x54400000
17// pio2_1t: pi/2 - pio2_1
18const pio2_1t = 6.07710050650619224932e-11; // 0x3DD0B461, 0x1A626331
19// pio2_2: second 33 bit of pi/2
20const pio2_2 = 6.07710050630396597660e-11; // 0x3DD0B461, 0x1A600000
21// pio2_2t: pi/2 - (pio2_1+pio2_2)
22const pio2_2t = 2.02226624879595063154e-21; // 0x3BA3198A, 0x2E037073
23// pio2_3: third 33 bit of pi/2
24const pio2_3 = 2.02226624871116645580e-21; // 0x3BA3198A, 0x2E000000
25// pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3)
26const pio2_3t = 8.47842766036889956997e-32; // 0x397B839A, 0x252049C1
27
28fn U(x: anytype) usize {
29 return @intCast(usize, x);
30}
31
32fn medium(ix: u32, x: f64, y: *[2]f64) i32 {
33 var w: f64 = undefined;
34 var t: f64 = undefined;
35 var r: f64 = undefined;
36 var @"fn": f64 = undefined;
37 var n: i32 = undefined;
38 var ex: i32 = undefined;
39 var ey: i32 = undefined;
40 var ui: u64 = undefined;
41
42 // rint(x/(pi/2))
43 @"fn" = x * invpio2 + toint - toint;
44 n = @floatToInt(i32, @"fn");
45 r = x - @"fn" * pio2_1;
46 w = @"fn" * pio2_1t; // 1st round, good to 85 bits
47 // Matters with directed rounding.
48 if (r - w < -pio4) {
49 n -= 1;
50 @"fn" -= 1;
51 r = x - @"fn" * pio2_1;
52 w = @"fn" * pio2_1t;
53 } else if (r - w > pio4) {
54 n += 1;
55 @"fn" += 1;
56 r = x - @"fn" * pio2_1;
57 w = @"fn" * pio2_1t;
58 }
59 y[0] = r - w;
60 ui = @bitCast(u64, y[0]);
61 ey = @intCast(i32, (ui >> 52) & 0x7ff);
62 ex = @intCast(i32, ix >> 20);
63 if (ex - ey > 16) { // 2nd round, good to 118 bits
64 t = r;
65 w = @"fn" * pio2_2;
66 r = t - w;
67 w = @"fn" * pio2_2t - ((t - r) - w);
68 y[0] = r - w;
69 ui = @bitCast(u64, y[0]);
70 ey = @intCast(i32, (ui >> 52) & 0x7ff);
71 if (ex - ey > 49) { // 3rd round, good to 151 bits, covers all cases
72 t = r;
73 w = @"fn" * pio2_3;
74 r = t - w;
75 w = @"fn" * pio2_3t - ((t - r) - w);
76 y[0] = r - w;
77 }
78 }
79 y[1] = (r - y[0]) - w;
80 return n;
81}
82
83// Returns the remainder of x rem pi/2 in y[0]+y[1]
84//
85// use __rem_pio2_large() for large x
86//
87// caller must handle the case when reduction is not needed: |x| ~<= pi/4 */
88pub fn __rem_pio2(x: f64, y: *[2]f64) i32 {
89 var z: f64 = undefined;
90 var tx: [3]f64 = undefined;
91 var ty: [2]f64 = undefined;
92 var n: i32 = undefined;
93 var ix: u32 = undefined;
94 var sign: bool = undefined;
95 var i: i32 = undefined;
96 var ui: u64 = undefined;
97
98 ui = @bitCast(u64, x);
99 sign = ui >> 63 != 0;
100 ix = @truncate(u32, (ui >> 32) & 0x7fffffff);
101 if (ix <= 0x400f6a7a) { // |x| ~<= 5pi/4
102 if ((ix & 0xfffff) == 0x921fb) { // |x| ~= pi/2 or 2pi/2
103 return medium(ix, x, y);
104 }
105 if (ix <= 0x4002d97c) { // |x| ~<= 3pi/4
106 if (!sign) {
107 z = x - pio2_1; // one round good to 85 bits
108 y[0] = z - pio2_1t;
109 y[1] = (z - y[0]) - pio2_1t;
110 return 1;
111 } else {
112 z = x + pio2_1;
113 y[0] = z + pio2_1t;
114 y[1] = (z - y[0]) + pio2_1t;
115 return -1;
116 }
117 } else {
118 if (!sign) {
119 z = x - 2 * pio2_1;
120 y[0] = z - 2 * pio2_1t;
121 y[1] = (z - y[0]) - 2 * pio2_1t;
122 return 2;
123 } else {
124 z = x + 2 * pio2_1;
125 y[0] = z + 2 * pio2_1t;
126 y[1] = (z - y[0]) + 2 * pio2_1t;
127 return -2;
128 }
129 }
130 }
131 if (ix <= 0x401c463b) { // |x| ~<= 9pi/4
132 if (ix <= 0x4015fdbc) { // |x| ~<= 7pi/4
133 if (ix == 0x4012d97c) { // |x| ~= 3pi/2
134 return medium(ix, x, y);
135 }
136 if (!sign) {
137 z = x - 3 * pio2_1;
138 y[0] = z - 3 * pio2_1t;
139 y[1] = (z - y[0]) - 3 * pio2_1t;
140 return 3;
141 } else {
142 z = x + 3 * pio2_1;
143 y[0] = z + 3 * pio2_1t;
144 y[1] = (z - y[0]) + 3 * pio2_1t;
145 return -3;
146 }
147 } else {
148 if (ix == 0x401921fb) { // |x| ~= 4pi/2 */
149 return medium(ix, x, y);
150 }
151 if (!sign) {
152 z = x - 4 * pio2_1;
153 y[0] = z - 4 * pio2_1t;
154 y[1] = (z - y[0]) - 4 * pio2_1t;
155 return 4;
156 } else {
157 z = x + 4 * pio2_1;
158 y[0] = z + 4 * pio2_1t;
159 y[1] = (z - y[0]) + 4 * pio2_1t;
160 return -4;
161 }
162 }
163 }
164 if (ix < 0x413921fb) { // |x| ~< 2^20*(pi/2), medium size
165 return medium(ix, x, y);
166 }
167 // all other (large) arguments
168 if (ix >= 0x7ff00000) { // x is inf or NaN
169 y[0] = x - x;
170 y[1] = y[0];
171 return 0;
172 }
173 // set z = scalbn(|x|,-ilogb(x)+23)
174 ui = @bitCast(u64, x);
175 ui &= std.math.maxInt(u64) >> 12;
176 ui |= @as(u64, 0x3ff + 23) << 52;
177 z = @bitCast(f64, ui);
178
179 i = 0;
180 while (i < 2) : (i += 1) {
181 tx[U(i)] = @intToFloat(f64, @floatToInt(i32, z));
182 z = (z - tx[U(i)]) * 0x1p24;
183 }
184 tx[U(i)] = z;
185 // skip zero terms, first term is non-zero
186 while (tx[U(i)] == 0.0) {
187 i -= 1;
188 }
189 n = __rem_pio2_large(tx[0..], ty[0..], @intCast(i32, (ix >> 20)) - (0x3ff + 23), i + 1, 1);
190 if (sign) {
191 y[0] = -ty[0];
192 y[1] = -ty[1];
193 return -n;
194 }
195 y[0] = ty[0];
196 y[1] = ty[1];
197 return n;
198}
lib/std/math/__rem_pio2_large.zig created+510
...@@ -0,0 +1,510 @@
1// Ported from musl, which is licensed under the MIT license:
2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
3//
4// https://git.musl-libc.org/cgit/musl/tree/src/math/__rem_pio2_large.c
5
6const std = @import("../std.zig");
7const math = std.math;
8
9const init_jk = [_]i32{ 3, 4, 4, 6 }; // initial value for jk
10
11//
12// Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi
13//
14// integer array, contains the (24*i)-th to (24*i+23)-th
15// bit of 2/pi after binary point. The corresponding
16// floating value is
17//
18// ipio2[i] * 2^(-24(i+1)).
19//
20// NB: This table must have at least (e0-3)/24 + jk terms.
21// For quad precision (e0 <= 16360, jk = 6), this is 686.
22///
23const ipio2 = [_]i32{
24 0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62,
25 0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A,
26 0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129,
27 0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41,
28 0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8,
29 0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF,
30 0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5,
31 0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08,
32 0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3,
33 0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880,
34 0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B,
35
36 //#if LDBL_MAX_EXP > 1024
37 0x47C419, 0xC367CD, 0xDCE809, 0x2A8359, 0xC4768B, 0x961CA6,
38 0xDDAF44, 0xD15719, 0x053EA5, 0xFF0705, 0x3F7E33, 0xE832C2,
39 0xDE4F98, 0x327DBB, 0xC33D26, 0xEF6B1E, 0x5EF89F, 0x3A1F35,
40 0xCAF27F, 0x1D87F1, 0x21907C, 0x7C246A, 0xFA6ED5, 0x772D30,
41 0x433B15, 0xC614B5, 0x9D19C3, 0xC2C4AD, 0x414D2C, 0x5D000C,
42 0x467D86, 0x2D71E3, 0x9AC69B, 0x006233, 0x7CD2B4, 0x97A7B4,
43 0xD55537, 0xF63ED7, 0x1810A3, 0xFC764D, 0x2A9D64, 0xABD770,
44 0xF87C63, 0x57B07A, 0xE71517, 0x5649C0, 0xD9D63B, 0x3884A7,
45 0xCB2324, 0x778AD6, 0x23545A, 0xB91F00, 0x1B0AF1, 0xDFCE19,
46 0xFF319F, 0x6A1E66, 0x615799, 0x47FBAC, 0xD87F7E, 0xB76522,
47 0x89E832, 0x60BFE6, 0xCDC4EF, 0x09366C, 0xD43F5D, 0xD7DE16,
48 0xDE3B58, 0x929BDE, 0x2822D2, 0xE88628, 0x4D58E2, 0x32CAC6,
49 0x16E308, 0xCB7DE0, 0x50C017, 0xA71DF3, 0x5BE018, 0x34132E,
50 0x621283, 0x014883, 0x5B8EF5, 0x7FB0AD, 0xF2E91E, 0x434A48,
51 0xD36710, 0xD8DDAA, 0x425FAE, 0xCE616A, 0xA4280A, 0xB499D3,
52 0xF2A606, 0x7F775C, 0x83C2A3, 0x883C61, 0x78738A, 0x5A8CAF,
53 0xBDD76F, 0x63A62D, 0xCBBFF4, 0xEF818D, 0x67C126, 0x45CA55,
54 0x36D9CA, 0xD2A828, 0x8D61C2, 0x77C912, 0x142604, 0x9B4612,
55 0xC459C4, 0x44C5C8, 0x91B24D, 0xF31700, 0xAD43D4, 0xE54929,
56 0x10D5FD, 0xFCBE00, 0xCC941E, 0xEECE70, 0xF53E13, 0x80F1EC,
57 0xC3E7B3, 0x28F8C7, 0x940593, 0x3E71C1, 0xB3092E, 0xF3450B,
58 0x9C1288, 0x7B20AB, 0x9FB52E, 0xC29247, 0x2F327B, 0x6D550C,
59 0x90A772, 0x1FE76B, 0x96CB31, 0x4A1679, 0xE27941, 0x89DFF4,
60 0x9794E8, 0x84E6E2, 0x973199, 0x6BED88, 0x365F5F, 0x0EFDBB,
61 0xB49A48, 0x6CA467, 0x427271, 0x325D8D, 0xB8159F, 0x09E5BC,
62 0x25318D, 0x3974F7, 0x1C0530, 0x010C0D, 0x68084B, 0x58EE2C,
63 0x90AA47, 0x02E774, 0x24D6BD, 0xA67DF7, 0x72486E, 0xEF169F,
64 0xA6948E, 0xF691B4, 0x5153D1, 0xF20ACF, 0x339820, 0x7E4BF5,
65 0x6863B2, 0x5F3EDD, 0x035D40, 0x7F8985, 0x295255, 0xC06437,
66 0x10D86D, 0x324832, 0x754C5B, 0xD4714E, 0x6E5445, 0xC1090B,
67 0x69F52A, 0xD56614, 0x9D0727, 0x50045D, 0xDB3BB4, 0xC576EA,
68 0x17F987, 0x7D6B49, 0xBA271D, 0x296996, 0xACCCC6, 0x5414AD,
69 0x6AE290, 0x89D988, 0x50722C, 0xBEA404, 0x940777, 0x7030F3,
70 0x27FC00, 0xA871EA, 0x49C266, 0x3DE064, 0x83DD97, 0x973FA3,
71 0xFD9443, 0x8C860D, 0xDE4131, 0x9D3992, 0x8C70DD, 0xE7B717,
72 0x3BDF08, 0x2B3715, 0xA0805C, 0x93805A, 0x921110, 0xD8E80F,
73 0xAF806C, 0x4BFFDB, 0x0F9038, 0x761859, 0x15A562, 0xBBCB61,
74 0xB989C7, 0xBD4010, 0x04F2D2, 0x277549, 0xF6B6EB, 0xBB22DB,
75 0xAA140A, 0x2F2689, 0x768364, 0x333B09, 0x1A940E, 0xAA3A51,
76 0xC2A31D, 0xAEEDAF, 0x12265C, 0x4DC26D, 0x9C7A2D, 0x9756C0,
77 0x833F03, 0xF6F009, 0x8C402B, 0x99316D, 0x07B439, 0x15200C,
78 0x5BC3D8, 0xC492F5, 0x4BADC6, 0xA5CA4E, 0xCD37A7, 0x36A9E6,
79 0x9492AB, 0x6842DD, 0xDE6319, 0xEF8C76, 0x528B68, 0x37DBFC,
80 0xABA1AE, 0x3115DF, 0xA1AE00, 0xDAFB0C, 0x664D64, 0xB705ED,
81 0x306529, 0xBF5657, 0x3AFF47, 0xB9F96A, 0xF3BE75, 0xDF9328,
82 0x3080AB, 0xF68C66, 0x15CB04, 0x0622FA, 0x1DE4D9, 0xA4B33D,
83 0x8F1B57, 0x09CD36, 0xE9424E, 0xA4BE13, 0xB52333, 0x1AAAF0,
84 0xA8654F, 0xA5C1D2, 0x0F3F0B, 0xCD785B, 0x76F923, 0x048B7B,
85 0x721789, 0x53A6C6, 0xE26E6F, 0x00EBEF, 0x584A9B, 0xB7DAC4,
86 0xBA66AA, 0xCFCF76, 0x1D02D1, 0x2DF1B1, 0xC1998C, 0x77ADC3,
87 0xDA4886, 0xA05DF7, 0xF480C6, 0x2FF0AC, 0x9AECDD, 0xBC5C3F,
88 0x6DDED0, 0x1FC790, 0xB6DB2A, 0x3A25A3, 0x9AAF00, 0x9353AD,
89 0x0457B6, 0xB42D29, 0x7E804B, 0xA707DA, 0x0EAA76, 0xA1597B,
90 0x2A1216, 0x2DB7DC, 0xFDE5FA, 0xFEDB89, 0xFDBE89, 0x6C76E4,
91 0xFCA906, 0x70803E, 0x156E85, 0xFF87FD, 0x073E28, 0x336761,
92 0x86182A, 0xEABD4D, 0xAFE7B3, 0x6E6D8F, 0x396795, 0x5BBF31,
93 0x48D784, 0x16DF30, 0x432DC7, 0x356125, 0xCE70C9, 0xB8CB30,
94 0xFD6CBF, 0xA200A4, 0xE46C05, 0xA0DD5A, 0x476F21, 0xD21262,
95 0x845CB9, 0x496170, 0xE0566B, 0x015299, 0x375550, 0xB7D51E,
96 0xC4F133, 0x5F6E13, 0xE4305D, 0xA92E85, 0xC3B21D, 0x3632A1,
97 0xA4B708, 0xD4B1EA, 0x21F716, 0xE4698F, 0x77FF27, 0x80030C,
98 0x2D408D, 0xA0CD4F, 0x99A520, 0xD3A2B3, 0x0A5D2F, 0x42F9B4,
99 0xCBDA11, 0xD0BE7D, 0xC1DB9B, 0xBD17AB, 0x81A2CA, 0x5C6A08,
100 0x17552E, 0x550027, 0xF0147F, 0x8607E1, 0x640B14, 0x8D4196,
101 0xDEBE87, 0x2AFDDA, 0xB6256B, 0x34897B, 0xFEF305, 0x9EBFB9,
102 0x4F6A68, 0xA82A4A, 0x5AC44F, 0xBCF82D, 0x985AD7, 0x95C7F4,
103 0x8D4D0D, 0xA63A20, 0x5F57A4, 0xB13F14, 0x953880, 0x0120CC,
104 0x86DD71, 0xB6DEC9, 0xF560BF, 0x11654D, 0x6B0701, 0xACB08C,
105 0xD0C0B2, 0x485551, 0x0EFB1E, 0xC37295, 0x3B06A3, 0x3540C0,
106 0x7BDC06, 0xCC45E0, 0xFA294E, 0xC8CAD6, 0x41F3E8, 0xDE647C,
107 0xD8649B, 0x31BED9, 0xC397A4, 0xD45877, 0xC5E369, 0x13DAF0,
108 0x3C3ABA, 0x461846, 0x5F7555, 0xF5BDD2, 0xC6926E, 0x5D2EAC,
109 0xED440E, 0x423E1C, 0x87C461, 0xE9FD29, 0xF3D6E7, 0xCA7C22,
110 0x35916F, 0xC5E008, 0x8DD7FF, 0xE26A6E, 0xC6FDB0, 0xC10893,
111 0x745D7C, 0xB2AD6B, 0x9D6ECD, 0x7B723E, 0x6A11C6, 0xA9CFF7,
112 0xDF7329, 0xBAC9B5, 0x5100B7, 0x0DB2E2, 0x24BA74, 0x607DE5,
113 0x8AD874, 0x2C150D, 0x0C1881, 0x94667E, 0x162901, 0x767A9F,
114 0xBEFDFD, 0xEF4556, 0x367ED9, 0x13D9EC, 0xB9BA8B, 0xFC97C4,
115 0x27A831, 0xC36EF1, 0x36C594, 0x56A8D8, 0xB5A8B4, 0x0ECCCF,
116 0x2D8912, 0x34576F, 0x89562C, 0xE3CE99, 0xB920D6, 0xAA5E6B,
117 0x9C2A3E, 0xCC5F11, 0x4A0BFD, 0xFBF4E1, 0x6D3B8E, 0x2C86E2,
118 0x84D4E9, 0xA9B4FC, 0xD1EEEF, 0xC9352E, 0x61392F, 0x442138,
119 0xC8D91B, 0x0AFC81, 0x6A4AFB, 0xD81C2F, 0x84B453, 0x8C994E,
120 0xCC2254, 0xDC552A, 0xD6C6C0, 0x96190B, 0xB8701A, 0x649569,
121 0x605A26, 0xEE523F, 0x0F117F, 0x11B5F4, 0xF5CBFC, 0x2DBC34,
122 0xEEBC34, 0xCC5DE8, 0x605EDD, 0x9B8E67, 0xEF3392, 0xB817C9,
123 0x9B5861, 0xBC57E1, 0xC68351, 0x103ED8, 0x4871DD, 0xDD1C2D,
124 0xA118AF, 0x462C21, 0xD7F359, 0x987AD9, 0xC0549E, 0xFA864F,
125 0xFC0656, 0xAE79E5, 0x362289, 0x22AD38, 0xDC9367, 0xAAE855,
126 0x382682, 0x9BE7CA, 0xA40D51, 0xB13399, 0x0ED7A9, 0x480569,
127 0xF0B265, 0xA7887F, 0x974C88, 0x36D1F9, 0xB39221, 0x4A827B,
128 0x21CF98, 0xDC9F40, 0x5547DC, 0x3A74E1, 0x42EB67, 0xDF9DFE,
129 0x5FD45E, 0xA4677B, 0x7AACBA, 0xA2F655, 0x23882B, 0x55BA41,
130 0x086E59, 0x862A21, 0x834739, 0xE6E389, 0xD49EE5, 0x40FB49,
131 0xE956FF, 0xCA0F1C, 0x8A59C5, 0x2BFA94, 0xC5C1D3, 0xCFC50F,
132 0xAE5ADB, 0x86C547, 0x624385, 0x3B8621, 0x94792C, 0x876110,
133 0x7B4C2A, 0x1A2C80, 0x12BF43, 0x902688, 0x893C78, 0xE4C4A8,
134 0x7BDBE5, 0xC23AC4, 0xEAF426, 0x8A67F7, 0xBF920D, 0x2BA365,
135 0xB1933D, 0x0B7CBD, 0xDC51A4, 0x63DD27, 0xDDE169, 0x19949A,
136 0x9529A8, 0x28CE68, 0xB4ED09, 0x209F44, 0xCA984E, 0x638270,
137 0x237C7E, 0x32B90F, 0x8EF5A7, 0xE75614, 0x08F121, 0x2A9DB5,
138 0x4D7E6F, 0x5119A5, 0xABF9B5, 0xD6DF82, 0x61DD96, 0x023616,
139 0x9F3AC4, 0xA1A283, 0x6DED72, 0x7A8D39, 0xA9B882, 0x5C326B,
140 0x5B2746, 0xED3400, 0x7700D2, 0x55F4FC, 0x4D5901,
141 0x8071E0,
142 //#endif
143};
144
145const PIo2 = [_]f64{
146 1.57079625129699707031e+00, // 0x3FF921FB, 0x40000000
147 7.54978941586159635335e-08, // 0x3E74442D, 0x00000000
148 5.39030252995776476554e-15, // 0x3CF84698, 0x80000000
149 3.28200341580791294123e-22, // 0x3B78CC51, 0x60000000
150 1.27065575308067607349e-29, // 0x39F01B83, 0x80000000
151 1.22933308981111328932e-36, // 0x387A2520, 0x40000000
152 2.73370053816464559624e-44, // 0x36E38222, 0x80000000
153 2.16741683877804819444e-51, // 0x3569F31D, 0x00000000
154};
155
156fn U(x: anytype) usize {
157 return @intCast(usize, x);
158}
159
160// Returns the last three digits of N with y = x - N*pi/2 so that |y| < pi/2.
161//
162// The method is to compute the integer (mod 8) and fraction parts of
163// (2/pi)*x without doing the full multiplication. In general we
164// skip the part of the product that are known to be a huge integer (
165// more accurately, = 0 mod 8 ). Thus the number of operations are
166// independent of the exponent of the input.
167//
168// (2/pi) is represented by an array of 24-bit integers in ipio2[].
169//
170// Input parameters:
171// x[] The input value (must be positive) is broken into nx
172// pieces of 24-bit integers in double precision format.
173// x[i] will be the i-th 24 bit of x. The scaled exponent
174// of x[0] is given in input parameter e0 (i.e., x[0]*2^e0
175// match x's up to 24 bits.
176//
177// Example of breaking a double positive z into x[0]+x[1]+x[2]:
178// e0 = ilogb(z)-23
179// z = scalbn(z,-e0)
180// for i = 0,1,2
181// x[i] = floor(z)
182// z = (z-x[i])*2**24
183//
184//
185// y[] ouput result in an array of double precision numbers.
186// The dimension of y[] is:
187// 24-bit precision 1
188// 53-bit precision 2
189// 64-bit precision 2
190// 113-bit precision 3
191// The actual value is the sum of them. Thus for 113-bit
192// precison, one may have to do something like:
193//
194// long double t,w,r_head, r_tail;
195// t = (long double)y[2] + (long double)y[1];
196// w = (long double)y[0];
197// r_head = t+w;
198// r_tail = w - (r_head - t);
199//
200// e0 The exponent of x[0]. Must be <= 16360 or you need to
201// expand the ipio2 table.
202//
203// nx dimension of x[]
204//
205// prec an integer indicating the precision:
206// 0 24 bits (single)
207// 1 53 bits (double)
208// 2 64 bits (extended)
209// 3 113 bits (quad)
210//
211// Here is the description of some local variables:
212//
213// jk jk+1 is the initial number of terms of ipio2[] needed
214// in the computation. The minimum and recommended value
215// for jk is 3,4,4,6 for single, double, extended, and quad.
216// jk+1 must be 2 larger than you might expect so that our
217// recomputation test works. (Up to 24 bits in the integer
218// part (the 24 bits of it that we compute) and 23 bits in
219// the fraction part may be lost to cancelation before we
220// recompute.)
221//
222// jz local integer variable indicating the number of
223// terms of ipio2[] used.
224//
225// jx nx - 1
226//
227// jv index for pointing to the suitable ipio2[] for the
228// computation. In general, we want
229// ( 2^e0*x[0] * ipio2[jv-1]*2^(-24jv) )/8
230// is an integer. Thus
231// e0-3-24*jv >= 0 or (e0-3)/24 >= jv
232// Hence jv = max(0,(e0-3)/24).
233//
234// jp jp+1 is the number of terms in PIo2[] needed, jp = jk.
235//
236// q[] double array with integral value, representing the
237// 24-bits chunk of the product of x and 2/pi.
238//
239// q0 the corresponding exponent of q[0]. Note that the
240// exponent for q[i] would be q0-24*i.
241//
242// PIo2[] double precision array, obtained by cutting pi/2
243// into 24 bits chunks.
244//
245// f[] ipio2[] in floating point
246//
247// iq[] integer array by breaking up q[] in 24-bits chunk.
248//
249// fq[] final product of x*(2/pi) in fq[0],..,fq[jk]
250//
251// ih integer. If >0 it indicates q[] is >= 0.5, hence
252// it also indicates the *sign* of the result.
253//
254///
255//
256// Constants:
257// The hexadecimal values are the intended ones for the following
258// constants. The decimal values may be used, provided that the
259// compiler will convert from decimal to binary accurately enough
260// to produce the hexadecimal values shown.
261///
262pub fn __rem_pio2_large(x: []f64, y: []f64, e0: i32, nx: i32, prec: usize) i32 {
263 var jz: i32 = undefined;
264 var jx: i32 = undefined;
265 var jv: i32 = undefined;
266 var jp: i32 = undefined;
267 var jk: i32 = undefined;
268 var carry: i32 = undefined;
269 var n: i32 = undefined;
270 var iq: [20]i32 = undefined;
271 var i: i32 = undefined;
272 var j: i32 = undefined;
273 var k: i32 = undefined;
274 var m: i32 = undefined;
275 var q0: i32 = undefined;
276 var ih: i32 = undefined;
277
278 var z: f64 = undefined;
279 var fw: f64 = undefined;
280 var f: [20]f64 = undefined;
281 var fq: [20]f64 = undefined;
282 var q: [20]f64 = undefined;
283
284 // initialize jk
285 jk = init_jk[prec];
286 jp = jk;
287
288 // determine jx,jv,q0, note that 3>q0
289 jx = nx - 1;
290 jv = @divFloor(e0 - 3, 24);
291 if (jv < 0) jv = 0;
292 q0 = e0 - 24 * (jv + 1);
293
294 // set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk]
295 j = jv - jx;
296 m = jx + jk;
297 i = 0;
298 while (i <= m) : ({
299 i += 1;
300 j += 1;
301 }) {
302 f[U(i)] = if (j < 0) 0.0 else @intToFloat(f64, ipio2[U(j)]);
303 }
304
305 // compute q[0],q[1],...q[jk]
306 i = 0;
307 while (i <= jk) : (i += 1) {
308 j = 0;
309 fw = 0;
310 while (j <= jx) : (j += 1) {
311 fw += x[U(j)] * f[U(jx + i - j)];
312 }
313 q[U(i)] = fw;
314 }
315
316 jz = jk;
317
318 // This is to handle a non-trivial goto translation from C.
319 // An unconditional return statement is found at the end of this loop.
320 recompute: while (true) {
321 // distill q[] into iq[] reversingly
322 i = 0;
323 j = jz;
324 z = q[U(jz)];
325 while (j > 0) : ({
326 i += 1;
327 j -= 1;
328 }) {
329 fw = @intToFloat(f64, @floatToInt(i32, 0x1p-24 * z));
330 iq[U(i)] = @floatToInt(i32, z - 0x1p24 * fw);
331 z = q[U(j - 1)] + fw;
332 }
333
334 // compute n
335 z = math.scalbn(z, q0); // actual value of z
336 z -= 8.0 * math.floor(z * 0.125); // trim off integer >= 8
337 n = @floatToInt(i32, z);
338 z -= @intToFloat(f64, n);
339 ih = 0;
340 if (q0 > 0) { // need iq[jz-1] to determine n
341 i = iq[U(jz - 1)] >> @intCast(u5, 24 - q0);
342 n += i;
343 iq[U(jz - 1)] -= i << @intCast(u5, 24 - q0);
344 ih = iq[U(jz - 1)] >> @intCast(u5, 23 - q0);
345 } else if (q0 == 0) {
346 ih = iq[U(jz - 1)] >> 23;
347 } else if (z >= 0.5) {
348 ih = 2;
349 }
350
351 if (ih > 0) { // q > 0.5
352 n += 1;
353 carry = 0;
354 i = 0;
355 while (i < jz) : (i += 1) { // compute 1-q
356 j = iq[U(i)];
357 if (carry == 0) {
358 if (j != 0) {
359 carry = 1;
360 iq[U(i)] = 0x1000000 - j;
361 }
362 } else {
363 iq[U(i)] = 0xffffff - j;
364 }
365 }
366 if (q0 > 0) { // rare case: chance is 1 in 12
367 switch (q0) {
368 1 => iq[U(jz - 1)] &= 0x7fffff,
369 2 => iq[U(jz - 1)] &= 0x3fffff,
370 else => unreachable,
371 }
372 }
373 if (ih == 2) {
374 z = 1.0 - z;
375 if (carry != 0) {
376 z -= math.scalbn(@as(f64, 1.0), q0);
377 }
378 }
379 }
380
381 // check if recomputation is needed
382 if (z == 0.0) {
383 j = 0;
384 i = jz - 1;
385 while (i >= jk) : (i -= 1) {
386 j |= iq[U(i)];
387 }
388
389 if (j == 0) { // need recomputation
390 k = 1;
391 while (iq[U(jk - k)] == 0) : (k += 1) {
392 // k = no. of terms needed
393 }
394
395 i = jz + 1;
396 while (i <= jz + k) : (i += 1) { // add q[jz+1] to q[jz+k]
397 f[U(jx + i)] = @intToFloat(f64, ipio2[U(jv + i)]);
398 j = 0;
399 fw = 0;
400 while (j <= jx) : (j += 1) {
401 fw += x[U(j)] * f[U(jx + i - j)];
402 }
403 q[U(i)] = fw;
404 }
405 jz += k;
406 continue :recompute; // mimic goto recompute
407 }
408 }
409
410 // chop off zero terms
411 if (z == 0.0) {
412 jz -= 1;
413 q0 -= 24;
414 while (iq[U(jz)] == 0) {
415 jz -= 1;
416 q0 -= 24;
417 }
418 } else { // break z into 24-bit if necessary
419 z = math.scalbn(z, -q0);
420 if (z >= 0x1p24) {
421 fw = @intToFloat(f64, @floatToInt(i32, 0x1p-24 * z));
422 iq[U(jz)] = @floatToInt(i32, z - 0x1p24 * fw);
423 jz += 1;
424 q0 += 24;
425 iq[U(jz)] = @floatToInt(i32, fw);
426 } else {
427 iq[U(jz)] = @floatToInt(i32, z);
428 }
429 }
430
431 // convert integer "bit" chunk to floating-point value
432 fw = math.scalbn(@as(f64, 1.0), q0);
433 i = jz;
434 while (i >= 0) : (i -= 1) {
435 q[U(i)] = fw * @intToFloat(f64, iq[U(i)]);
436 fw *= 0x1p-24;
437 }
438
439 // compute PIo2[0,...,jp]*q[jz,...,0]
440 i = jz;
441 while (i >= 0) : (i -= 1) {
442 fw = 0;
443 k = 0;
444 while (k <= jp and k <= jz - i) : (k += 1) {
445 fw += PIo2[U(k)] * q[U(i + k)];
446 }
447 fq[U(jz - i)] = fw;
448 }
449
450 // compress fq[] into y[]
451 switch (prec) {
452 0 => {
453 fw = 0.0;
454 i = jz;
455 while (i >= 0) : (i -= 1) {
456 fw += fq[U(i)];
457 }
458 y[0] = if (ih == 0) fw else -fw;
459 },
460
461 1, 2 => {
462 fw = 0.0;
463 i = jz;
464 while (i >= 0) : (i -= 1) {
465 fw += fq[U(i)];
466 }
467 // TODO: drop excess precision here once double_t is used
468 fw = fw;
469 y[0] = if (ih == 0) fw else -fw;
470 fw = fq[0] - fw;
471 i = 1;
472 while (i <= jz) : (i += 1) {
473 fw += fq[U(i)];
474 }
475 y[1] = if (ih == 0) fw else -fw;
476 },
477 3 => { // painful
478 i = jz;
479 while (i > 0) : (i -= 1) {
480 fw = fq[U(i - 1)] + fq[U(i)];
481 fq[U(i)] += fq[U(i - 1)] - fw;
482 fq[U(i - 1)] = fw;
483 }
484 i = jz;
485 while (i > 1) : (i -= 1) {
486 fw = fq[U(i - 1)] + fq[U(i)];
487 fq[U(i)] += fq[U(i - 1)] - fw;
488 fq[U(i - 1)] = fw;
489 }
490 fw = 0;
491 i = jz;
492 while (i >= 2) : (i -= 1) {
493 fw += fq[U(i)];
494 }
495 if (ih == 0) {
496 y[0] = fq[0];
497 y[1] = fq[1];
498 y[2] = fw;
499 } else {
500 y[0] = -fq[0];
501 y[1] = -fq[1];
502 y[2] = -fw;
503 }
504 },
505 else => unreachable,
506 }
507
508 return n & 7;
509 }
510}
lib/std/math/__rem_pio2f.zig created+70
...@@ -0,0 +1,70 @@
1// Ported from musl, which is licensed under the MIT license:
2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
3//
4// https://git.musl-libc.org/cgit/musl/tree/src/math/__rem_pio2f.c
5
6const std = @import("../std.zig");
7const __rem_pio2_large = @import("__rem_pio2_large.zig").__rem_pio2_large;
8const math = std.math;
9
10const toint = 1.5 / math.epsilon(f64);
11// pi/4
12const pio4 = 0x1.921fb6p-1;
13// invpio2: 53 bits of 2/pi
14const invpio2 = 6.36619772367581382433e-01; // 0x3FE45F30, 0x6DC9C883
15// pio2_1: first 25 bits of pi/2
16const pio2_1 = 1.57079631090164184570e+00; // 0x3FF921FB, 0x50000000
17// pio2_1t: pi/2 - pio2_1
18const pio2_1t = 1.58932547735281966916e-08; // 0x3E5110b4, 0x611A6263
19
20// Returns the remainder of x rem pi/2 in *y
21// use double precision for everything except passing x
22// use __rem_pio2_large() for large x
23pub fn __rem_pio2f(x: f32, y: *f64) i32 {
24 var tx: [1]f64 = undefined;
25 var ty: [1]f64 = undefined;
26 var @"fn": f64 = undefined;
27 var ix: u32 = undefined;
28 var n: i32 = undefined;
29 var sign: bool = undefined;
30 var e0: u32 = undefined;
31 var ui: u32 = undefined;
32
33 ui = @bitCast(u32, x);
34 ix = ui & 0x7fffffff;
35
36 // 25+53 bit pi is good enough for medium size
37 if (ix < 0x4dc90fdb) { // |x| ~< 2^28*(pi/2), medium size
38 // Use a specialized rint() to get fn.
39 @"fn" = @floatCast(f64, x) * invpio2 + toint - toint;
40 n = @floatToInt(i32, @"fn");
41 y.* = x - @"fn" * pio2_1 - @"fn" * pio2_1t;
42 // Matters with directed rounding.
43 if (y.* < -pio4) {
44 n -= 1;
45 @"fn" -= 1;
46 y.* = x - @"fn" * pio2_1 - @"fn" * pio2_1t;
47 } else if (y.* > pio4) {
48 n += 1;
49 @"fn" += 1;
50 y.* = x - @"fn" * pio2_1 - @"fn" * pio2_1t;
51 }
52 return n;
53 }
54 if (ix >= 0x7f800000) { // x is inf or NaN
55 y.* = x - x;
56 return 0;
57 }
58 // scale x into [2^23, 2^24-1]
59 sign = ui >> 31 != 0;
60 e0 = (ix >> 23) - (0x7f + 23); // e0 = ilogb(|x|)-23, positive
61 ui = ix - (e0 << 23);
62 tx[0] = @bitCast(f32, ui);
63 n = __rem_pio2_large(&tx, &ty, @intCast(i32, e0), 1, 0);
64 if (sign) {
65 y.* = -ty[0];
66 return -n;
67 }
68 y.* = ty[0];
69 return n;
70}
lib/std/math/__trig.zig created+273
...@@ -0,0 +1,273 @@
1// Ported from musl, which is licensed under the MIT license:
2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
3//
4// https://git.musl-libc.org/cgit/musl/tree/src/math/__cos.c
5// https://git.musl-libc.org/cgit/musl/tree/src/math/__cosdf.c
6// https://git.musl-libc.org/cgit/musl/tree/src/math/__sin.c
7// https://git.musl-libc.org/cgit/musl/tree/src/math/__sindf.c
8// https://git.musl-libc.org/cgit/musl/tree/src/math/__tand.c
9// https://git.musl-libc.org/cgit/musl/tree/src/math/__tandf.c
10
11// kernel cos function on [-pi/4, pi/4], pi/4 ~ 0.785398164
12// Input x is assumed to be bounded by ~pi/4 in magnitude.
13// Input y is the tail of x.
14//
15// Algorithm
16// 1. Since cos(-x) = cos(x), we need only to consider positive x.
17// 2. if x < 2^-27 (hx<0x3e400000 0), return 1 with inexact if x!=0.
18// 3. cos(x) is approximated by a polynomial of degree 14 on
19// [0,pi/4]
20// 4 14
21// cos(x) ~ 1 - x*x/2 + C1*x + ... + C6*x
22// where the remez error is
23//
24// | 2 4 6 8 10 12 14 | -58
25// |cos(x)-(1-.5*x +C1*x +C2*x +C3*x +C4*x +C5*x +C6*x )| <= 2
26// | |
27//
28// 4 6 8 10 12 14
29// 4. let r = C1*x +C2*x +C3*x +C4*x +C5*x +C6*x , then
30// cos(x) ~ 1 - x*x/2 + r
31// since cos(x+y) ~ cos(x) - sin(x)*y
32// ~ cos(x) - x*y,
33// a correction term is necessary in cos(x) and hence
34// cos(x+y) = 1 - (x*x/2 - (r - x*y))
35// For better accuracy, rearrange to
36// cos(x+y) ~ w + (tmp + (r-x*y))
37// where w = 1 - x*x/2 and tmp is a tiny correction term
38// (1 - x*x/2 == w + tmp exactly in infinite precision).
39// The exactness of w + tmp in infinite precision depends on w
40// and tmp having the same precision as x. If they have extra
41// precision due to compiler bugs, then the extra precision is
42// only good provided it is retained in all terms of the final
43// expression for cos(). Retention happens in all cases tested
44// under FreeBSD, so don't pessimize things by forcibly clipping
45// any extra precision in w.
46pub fn __cos(x: f64, y: f64) f64 {
47 const C1 = 4.16666666666666019037e-02; // 0x3FA55555, 0x5555554C
48 const C2 = -1.38888888888741095749e-03; // 0xBF56C16C, 0x16C15177
49 const C3 = 2.48015872894767294178e-05; // 0x3EFA01A0, 0x19CB1590
50 const C4 = -2.75573143513906633035e-07; // 0xBE927E4F, 0x809C52AD
51 const C5 = 2.08757232129817482790e-09; // 0x3E21EE9E, 0xBDB4B1C4
52 const C6 = -1.13596475577881948265e-11; // 0xBDA8FAE9, 0xBE8838D4
53
54 const z = x * x;
55 const zs = z * z;
56 const r = z * (C1 + z * (C2 + z * C3)) + zs * zs * (C4 + z * (C5 + z * C6));
57 const hz = 0.5 * z;
58 const w = 1.0 - hz;
59 return w + (((1.0 - w) - hz) + (z * r - x * y));
60}
61
62pub fn __cosdf(x: f64) f32 {
63 // |cos(x) - c(x)| < 2**-34.1 (~[-5.37e-11, 5.295e-11]).
64 const C0 = -0x1ffffffd0c5e81.0p-54; // -0.499999997251031003120
65 const C1 = 0x155553e1053a42.0p-57; // 0.0416666233237390631894
66 const C2 = -0x16c087e80f1e27.0p-62; // -0.00138867637746099294692
67 const C3 = 0x199342e0ee5069.0p-68; // 0.0000243904487962774090654
68
69 // Try to optimize for parallel evaluation as in __tandf.c.
70 const z = x * x;
71 const w = z * z;
72 const r = C2 + z * C3;
73 return @floatCast(f32, ((1.0 + z * C0) + w * C1) + (w * z) * r);
74}
75
76// kernel sin function on ~[-pi/4, pi/4] (except on -0), pi/4 ~ 0.7854
77// Input x is assumed to be bounded by ~pi/4 in magnitude.
78// Input y is the tail of x.
79// Input iy indicates whether y is 0. (if iy=0, y assume to be 0).
80//
81// Algorithm
82// 1. Since sin(-x) = -sin(x), we need only to consider positive x.
83// 2. Callers must return sin(-0) = -0 without calling here since our
84// odd polynomial is not evaluated in a way that preserves -0.
85// Callers may do the optimization sin(x) ~ x for tiny x.
86// 3. sin(x) is approximated by a polynomial of degree 13 on
87// [0,pi/4]
88// 3 13
89// sin(x) ~ x + S1*x + ... + S6*x
90// where
91//
92// |sin(x) 2 4 6 8 10 12 | -58
93// |----- - (1+S1*x +S2*x +S3*x +S4*x +S5*x +S6*x )| <= 2
94// | x |
95//
96// 4. sin(x+y) = sin(x) + sin'(x')*y
97// ~ sin(x) + (1-x*x/2)*y
98// For better accuracy, let
99// 3 2 2 2 2
100// r = x *(S2+x *(S3+x *(S4+x *(S5+x *S6))))
101// then 3 2
102// sin(x) = x + (S1*x + (x *(r-y/2)+y))
103pub fn __sin(x: f64, y: f64, iy: i32) f64 {
104 const S1 = -1.66666666666666324348e-01; // 0xBFC55555, 0x55555549
105 const S2 = 8.33333333332248946124e-03; // 0x3F811111, 0x1110F8A6
106 const S3 = -1.98412698298579493134e-04; // 0xBF2A01A0, 0x19C161D5
107 const S4 = 2.75573137070700676789e-06; // 0x3EC71DE3, 0x57B1FE7D
108 const S5 = -2.50507602534068634195e-08; // 0xBE5AE5E6, 0x8A2B9CEB
109 const S6 = 1.58969099521155010221e-10; // 0x3DE5D93A, 0x5ACFD57C
110
111 const z = x * x;
112 const w = z * z;
113 const r = S2 + z * (S3 + z * S4) + z * w * (S5 + z * S6);
114 const v = z * x;
115 if (iy == 0) {
116 return x + v * (S1 + z * r);
117 } else {
118 return x - ((z * (0.5 * y - v * r) - y) - v * S1);
119 }
120}
121
122pub fn __sindf(x: f64) f32 {
123 // |sin(x)/x - s(x)| < 2**-37.5 (~[-4.89e-12, 4.824e-12]).
124 const S1 = -0x15555554cbac77.0p-55; // -0.166666666416265235595
125 const S2 = 0x111110896efbb2.0p-59; // 0.0083333293858894631756
126 const S3 = -0x1a00f9e2cae774.0p-65; // -0.000198393348360966317347
127 const S4 = 0x16cd878c3b46a7.0p-71; // 0.0000027183114939898219064
128
129 // Try to optimize for parallel evaluation as in __tandf.c.
130 const z = x * x;
131 const w = z * z;
132 const r = S3 + z * S4;
133 const s = z * x;
134 return @floatCast(f32, (x + s * (S1 + z * S2)) + s * w * r);
135}
136
137// kernel tan function on ~[-pi/4, pi/4] (except on -0), pi/4 ~ 0.7854
138// Input x is assumed to be bounded by ~pi/4 in magnitude.
139// Input y is the tail of x.
140// Input odd indicates whether tan (if odd = 0) or -1/tan (if odd = 1) is returned.
141//
142// Algorithm
143// 1. Since tan(-x) = -tan(x), we need only to consider positive x.
144// 2. Callers must return tan(-0) = -0 without calling here since our
145// odd polynomial is not evaluated in a way that preserves -0.
146// Callers may do the optimization tan(x) ~ x for tiny x.
147// 3. tan(x) is approximated by a odd polynomial of degree 27 on
148// [0,0.67434]
149// 3 27
150// tan(x) ~ x + T1*x + ... + T13*x
151// where
152//
153// |tan(x) 2 4 26 | -59.2
154// |----- - (1+T1*x +T2*x +.... +T13*x )| <= 2
155// | x |
156//
157// Note: tan(x+y) = tan(x) + tan'(x)*y
158// ~ tan(x) + (1+x*x)*y
159// Therefore, for better accuracy in computing tan(x+y), let
160// 3 2 2 2 2
161// r = x *(T2+x *(T3+x *(...+x *(T12+x *T13))))
162// then
163// 3 2
164// tan(x+y) = x + (T1*x + (x *(r+y)+y))
165//
166// 4. For x in [0.67434,pi/4], let y = pi/4 - x, then
167// tan(x) = tan(pi/4-y) = (1-tan(y))/(1+tan(y))
168// = 1 - 2*(tan(y) - (tan(y)^2)/(1+tan(y)))
169pub fn __tan(x_: f64, y_: f64, odd: bool) f64 {
170 var x = x_;
171 var y = y_;
172
173 const T = [_]f64{
174 3.33333333333334091986e-01, // 3FD55555, 55555563
175 1.33333333333201242699e-01, // 3FC11111, 1110FE7A
176 5.39682539762260521377e-02, // 3FABA1BA, 1BB341FE
177 2.18694882948595424599e-02, // 3F9664F4, 8406D637
178 8.86323982359930005737e-03, // 3F8226E3, E96E8493
179 3.59207910759131235356e-03, // 3F6D6D22, C9560328
180 1.45620945432529025516e-03, // 3F57DBC8, FEE08315
181 5.88041240820264096874e-04, // 3F4344D8, F2F26501
182 2.46463134818469906812e-04, // 3F3026F7, 1A8D1068
183 7.81794442939557092300e-05, // 3F147E88, A03792A6
184 7.14072491382608190305e-05, // 3F12B80F, 32F0A7E9
185 -1.85586374855275456654e-05, // BEF375CB, DB605373
186 2.59073051863633712884e-05, // 3EFB2A70, 74BF7AD4
187 };
188 const pio4 = 7.85398163397448278999e-01; // 3FE921FB, 54442D18
189 const pio4lo = 3.06161699786838301793e-17; // 3C81A626, 33145C07
190
191 var z: f64 = undefined;
192 var r: f64 = undefined;
193 var v: f64 = undefined;
194 var w: f64 = undefined;
195 var s: f64 = undefined;
196 var a: f64 = undefined;
197 var w0: f64 = undefined;
198 var a0: f64 = undefined;
199 var hx: u32 = undefined;
200 var sign: bool = undefined;
201
202 hx = @intCast(u32, @bitCast(u64, x) >> 32);
203 const big = (hx & 0x7fffffff) >= 0x3FE59428; // |x| >= 0.6744
204 if (big) {
205 sign = hx >> 31 != 0;
206 if (sign) {
207 x = -x;
208 y = -y;
209 }
210 x = (pio4 - x) + (pio4lo - y);
211 y = 0.0;
212 }
213 z = x * x;
214 w = z * z;
215
216 // Break x^5*(T[1]+x^2*T[2]+...) into
217 // x^5(T[1]+x^4*T[3]+...+x^20*T[11]) +
218 // x^5(x^2*(T[2]+x^4*T[4]+...+x^22*[T12]))
219 r = T[1] + w * (T[3] + w * (T[5] + w * (T[7] + w * (T[9] + w * T[11]))));
220 v = z * (T[2] + w * (T[4] + w * (T[6] + w * (T[8] + w * (T[10] + w * T[12])))));
221 s = z * x;
222 r = y + z * (s * (r + v) + y) + s * T[0];
223 w = x + r;
224 if (big) {
225 s = 1 - 2 * @intToFloat(f64, @boolToInt(odd));
226 v = s - 2.0 * (x + (r - w * w / (w + s)));
227 return if (sign) -v else v;
228 }
229 if (!odd) {
230 return w;
231 }
232 // -1.0/(x+r) has up to 2ulp error, so compute it accurately
233 w0 = w;
234 w0 = @bitCast(f64, @bitCast(u64, w0) & 0xffffffff00000000);
235 v = r - (w0 - x); // w0+v = r+x
236 a = -1.0 / w;
237 a0 = a;
238 a0 = @bitCast(f64, @bitCast(u64, a0) & 0xffffffff00000000);
239 return a0 + a * (1.0 + a0 * w0 + a0 * v);
240}
241
242pub fn __tandf(x: f64, odd: bool) f32 {
243 // |tan(x)/x - t(x)| < 2**-25.5 (~[-2e-08, 2e-08]).
244 const T = [_]f64{
245 0x15554d3418c99f.0p-54, // 0.333331395030791399758
246 0x1112fd38999f72.0p-55, // 0.133392002712976742718
247 0x1b54c91d865afe.0p-57, // 0.0533812378445670393523
248 0x191df3908c33ce.0p-58, // 0.0245283181166547278873
249 0x185dadfcecf44e.0p-61, // 0.00297435743359967304927
250 0x1362b9bf971bcd.0p-59, // 0.00946564784943673166728
251 };
252
253 const z = x * x;
254 // Split up the polynomial into small independent terms to give
255 // opportunities for parallel evaluation. The chosen splitting is
256 // micro-optimized for Athlons (XP, X64). It costs 2 multiplications
257 // relative to Horner's method on sequential machines.
258 //
259 // We add the small terms from lowest degree up for efficiency on
260 // non-sequential machines (the lowest degree terms tend to be ready
261 // earlier). Apart from this, we don't care about order of
262 // operations, and don't need to to care since we have precision to
263 // spare. However, the chosen splitting is good for accuracy too,
264 // and would give results as accurate as Horner's method if the
265 // small terms were added from highest degree down.
266 const r = T[4] + z * T[5];
267 const t = T[2] + z * T[3];
268 const w = z * z;
269 const s = z * x;
270 const u = T[0] + z * T[1];
271 const r0 = (x + s * u) + (s * w) * (t + w * r);
272 return @floatCast(f32, if (odd) -1.0 / r0 else r0);
273}
lib/std/math/cos.zig+110-79
...@@ -1,12 +1,17 @@...@@ -1,12 +1,17 @@
1// Ported from go, which is licensed under a BSD-3 license.1// Ported from musl, which is licensed under the MIT license:
2// https://golang.org/LICENSE2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
3//3//
4// https://golang.org/src/math/sin.go4// https://git.musl-libc.org/cgit/musl/tree/src/math/cosf.c
5// https://git.musl-libc.org/cgit/musl/tree/src/math/cos.c
56
6const std = @import("../std.zig");7const std = @import("../std.zig");
7const math = std.math;8const math = std.math;
8const expect = std.testing.expect;9const expect = std.testing.expect;
910
11const kernel = @import("__trig.zig");
12const __rem_pio2 = @import("__rem_pio2.zig").__rem_pio2;
13const __rem_pio2f = @import("__rem_pio2f.zig").__rem_pio2f;
14
10/// Returns the cosine of the radian value x.15/// Returns the cosine of the radian value x.
11///16///
12/// Special Cases:17/// Special Cases:
...@@ -15,109 +20,135 @@ const expect = std.testing.expect;...@@ -15,109 +20,135 @@ const expect = std.testing.expect;
15pub fn cos(x: anytype) @TypeOf(x) {20pub fn cos(x: anytype) @TypeOf(x) {
16 const T = @TypeOf(x);21 const T = @TypeOf(x);
17 return switch (T) {22 return switch (T) {
18 f32 => cos_(f32, x),23 f32 => cos32(x),
19 f64 => cos_(f64, x),24 f64 => cos64(x),
20 else => @compileError("cos not implemented for " ++ @typeName(T)),25 else => @compileError("cos not implemented for " ++ @typeName(T)),
21 };26 };
22}27}
2328
24// sin polynomial coefficients29fn cos32(x: f32) f32 {
25const S0 = 1.58962301576546568060E-10;30 // Small multiples of pi/2 rounded to double precision.
26const S1 = -2.50507477628578072866E-8;31 const c1pio2: f64 = 1.0 * math.pi / 2.0; // 0x3FF921FB, 0x54442D18
27const S2 = 2.75573136213857245213E-6;32 const c2pio2: f64 = 2.0 * math.pi / 2.0; // 0x400921FB, 0x54442D18
28const S3 = -1.98412698295895385996E-4;33 const c3pio2: f64 = 3.0 * math.pi / 2.0; // 0x4012D97C, 0x7F3321D2
29const S4 = 8.33333333332211858878E-3;34 const c4pio2: f64 = 4.0 * math.pi / 2.0; // 0x401921FB, 0x54442D18
30const S5 = -1.66666666666666307295E-1;35
3136 var ix = @bitCast(u32, x);
32// cos polynomial coeffiecients37 const sign = ix >> 31 != 0;
33const C0 = -1.13585365213876817300E-11;38 ix &= 0x7fffffff;
34const C1 = 2.08757008419747316778E-9;39
35const C2 = -2.75573141792967388112E-7;40 if (ix <= 0x3f490fda) { // |x| ~<= pi/4
36const C3 = 2.48015872888517045348E-5;41 if (ix < 0x39800000) { // |x| < 2**-12
37const C4 = -1.38888888888730564116E-3;42 // raise inexact if x != 0
38const C5 = 4.16666666666665929218E-2;43 math.doNotOptimizeAway(x + 0x1p120);
3944 return 1.0;
40const pi4a = 7.85398125648498535156e-1;45 }
41const pi4b = 3.77489470793079817668E-8;46 return kernel.__cosdf(x);
42const pi4c = 2.69515142907905952645E-15;47 }
43const m4pi = 1.273239544735162542821171882678754627704620361328125;48 if (ix <= 0x407b53d1) { // |x| ~<= 5*pi/4
4449 if (ix > 0x4016cbe3) { // |x| ~> 3*pi/4
45fn cos_(comptime T: type, x_: T) T {50 return -kernel.__cosdf(if (sign) x + c2pio2 else x - c2pio2);
46 const I = std.meta.Int(.signed, @typeInfo(T).Float.bits);51 } else {
4752 if (sign) {
48 var x = x_;53 return kernel.__sindf(x + c1pio2);
49 if (math.isNan(x) or math.isInf(x)) {54 } else {
50 return math.nan(T);55 return kernel.__sindf(c1pio2 - x);
56 }
57 }
58 }
59 if (ix <= 0x40e231d5) { // |x| ~<= 9*pi/4
60 if (ix > 0x40afeddf) { // |x| ~> 7*pi/4
61 return kernel.__cosdf(if (sign) x + c4pio2 else x - c4pio2);
62 } else {
63 if (sign) {
64 return kernel.__sindf(-x - c3pio2);
65 } else {
66 return kernel.__sindf(x - c3pio2);
67 }
68 }
51 }69 }
5270
53 var sign = false;71 // cos(Inf or NaN) is NaN
54 x = math.fabs(x);72 if (ix >= 0x7f800000) {
73 return x - x;
74 }
5575
56 var y = math.floor(x * m4pi);76 var y: f64 = undefined;
57 var j = @floatToInt(I, y);77 const n = __rem_pio2f(x, &y);
78 return switch (n & 3) {
79 0 => kernel.__cosdf(y),
80 1 => kernel.__sindf(-y),
81 2 => -kernel.__cosdf(y),
82 else => kernel.__sindf(y),
83 };
84}
5885
59 if (j & 1 == 1) {86fn cos64(x: f64) f64 {
60 j += 1;87 var ix = @bitCast(u64, x) >> 32;
61 y += 1;88 ix &= 0x7fffffff;
89
90 // |x| ~< pi/4
91 if (ix <= 0x3fe921fb) {
92 if (ix < 0x3e46a09e) { // |x| < 2**-27 * sqrt(2)
93 // raise inexact if x!=0
94 math.doNotOptimizeAway(x + 0x1p120);
95 return 1.0;
96 }
97 return kernel.__cos(x, 0);
62 }98 }
6399
64 j &= 7;100 // cos(Inf or NaN) is NaN
65 if (j > 3) {101 if (ix >= 0x7ff00000) {
66 j -= 4;102 return x - x;
67 sign = !sign;
68 }103 }
69 if (j > 1) {
70 sign = !sign;
71 }
72
73 const z = ((x - y * pi4a) - y * pi4b) - y * pi4c;
74 const w = z * z;
75
76 const r = if (j == 1 or j == 2)
77 z + z * w * (S5 + w * (S4 + w * (S3 + w * (S2 + w * (S1 + w * S0)))))
78 else
79 1.0 - 0.5 * w + w * w * (C5 + w * (C4 + w * (C3 + w * (C2 + w * (C1 + w * C0)))));
80104
81 return if (sign) -r else r;105 var y: [2]f64 = undefined;
106 const n = __rem_pio2(x, &y);
107 return switch (n & 3) {
108 0 => kernel.__cos(y[0], y[1]),
109 1 => -kernel.__sin(y[0], y[1], 1),
110 2 => -kernel.__cos(y[0], y[1]),
111 else => kernel.__sin(y[0], y[1], 1),
112 };
82}113}
83114
84test "math.cos" {115test "math.cos" {
85 try expect(cos(@as(f32, 0.0)) == cos_(f32, 0.0));116 try expect(cos(@as(f32, 0.0)) == cos32(0.0));
86 try expect(cos(@as(f64, 0.0)) == cos_(f64, 0.0));117 try expect(cos(@as(f64, 0.0)) == cos64(0.0));
87}118}
88119
89test "math.cos32" {120test "math.cos32" {
90 const epsilon = 0.000001;121 const epsilon = 0.00001;
91122
92 try expect(math.approxEqAbs(f32, cos_(f32, 0.0), 1.0, epsilon));123 try expect(math.approxEqAbs(f32, cos32(0.0), 1.0, epsilon));
93 try expect(math.approxEqAbs(f32, cos_(f32, 0.2), 0.980067, epsilon));124 try expect(math.approxEqAbs(f32, cos32(0.2), 0.980067, epsilon));
94 try expect(math.approxEqAbs(f32, cos_(f32, 0.8923), 0.627623, epsilon));125 try expect(math.approxEqAbs(f32, cos32(0.8923), 0.627623, epsilon));
95 try expect(math.approxEqAbs(f32, cos_(f32, 1.5), 0.070737, epsilon));126 try expect(math.approxEqAbs(f32, cos32(1.5), 0.070737, epsilon));
96 try expect(math.approxEqAbs(f32, cos_(f32, -1.5), 0.070737, epsilon));127 try expect(math.approxEqAbs(f32, cos32(-1.5), 0.070737, epsilon));
97 try expect(math.approxEqAbs(f32, cos_(f32, 37.45), 0.969132, epsilon));128 try expect(math.approxEqAbs(f32, cos32(37.45), 0.969132, epsilon));
98 try expect(math.approxEqAbs(f32, cos_(f32, 89.123), 0.400798, epsilon));129 try expect(math.approxEqAbs(f32, cos32(89.123), 0.400798, epsilon));
99}130}
100131
101test "math.cos64" {132test "math.cos64" {
102 const epsilon = 0.000001;133 const epsilon = 0.000001;
103134
104 try expect(math.approxEqAbs(f64, cos_(f64, 0.0), 1.0, epsilon));135 try expect(math.approxEqAbs(f64, cos64(0.0), 1.0, epsilon));
105 try expect(math.approxEqAbs(f64, cos_(f64, 0.2), 0.980067, epsilon));136 try expect(math.approxEqAbs(f64, cos64(0.2), 0.980067, epsilon));
106 try expect(math.approxEqAbs(f64, cos_(f64, 0.8923), 0.627623, epsilon));137 try expect(math.approxEqAbs(f64, cos64(0.8923), 0.627623, epsilon));
107 try expect(math.approxEqAbs(f64, cos_(f64, 1.5), 0.070737, epsilon));138 try expect(math.approxEqAbs(f64, cos64(1.5), 0.070737, epsilon));
108 try expect(math.approxEqAbs(f64, cos_(f64, -1.5), 0.070737, epsilon));139 try expect(math.approxEqAbs(f64, cos64(-1.5), 0.070737, epsilon));
109 try expect(math.approxEqAbs(f64, cos_(f64, 37.45), 0.969132, epsilon));140 try expect(math.approxEqAbs(f64, cos64(37.45), 0.969132, epsilon));
110 try expect(math.approxEqAbs(f64, cos_(f64, 89.123), 0.40080, epsilon));141 try expect(math.approxEqAbs(f64, cos64(89.123), 0.40080, epsilon));
111}142}
112143
113test "math.cos32.special" {144test "math.cos32.special" {
114 try expect(math.isNan(cos_(f32, math.inf(f32))));145 try expect(math.isNan(cos32(math.inf(f32))));
115 try expect(math.isNan(cos_(f32, -math.inf(f32))));146 try expect(math.isNan(cos32(-math.inf(f32))));
116 try expect(math.isNan(cos_(f32, math.nan(f32))));147 try expect(math.isNan(cos32(math.nan(f32))));
117}148}
118149
119test "math.cos64.special" {150test "math.cos64.special" {
120 try expect(math.isNan(cos_(f64, math.inf(f64))));151 try expect(math.isNan(cos64(math.inf(f64))));
121 try expect(math.isNan(cos_(f64, -math.inf(f64))));152 try expect(math.isNan(cos64(-math.inf(f64))));
122 try expect(math.isNan(cos_(f64, math.nan(f64))));153 try expect(math.isNan(cos64(math.nan(f64))));
123}154}
lib/std/math/sin.zig+122-83
...@@ -1,12 +1,17 @@...@@ -1,12 +1,17 @@
1// Ported from go, which is licensed under a BSD-3 license.1// Ported from musl, which is licensed under the MIT license:
2// https://golang.org/LICENSE2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
3//
4// https://git.musl-libc.org/cgit/musl/tree/src/math/sinf.c
5// https://git.musl-libc.org/cgit/musl/tree/src/math/sin.c
3//6//
4// https://golang.org/src/math/sin.go
5
6const std = @import("../std.zig");7const std = @import("../std.zig");
7const math = std.math;8const math = std.math;
8const expect = std.testing.expect;9const expect = std.testing.expect;
910
11const kernel = @import("__trig.zig");
12const __rem_pio2 = @import("__rem_pio2.zig").__rem_pio2;
13const __rem_pio2f = @import("__rem_pio2f.zig").__rem_pio2f;
14
10/// Returns the sine of the radian value x.15/// Returns the sine of the radian value x.
11///16///
12/// Special Cases:17/// Special Cases:
...@@ -16,114 +21,148 @@ const expect = std.testing.expect;...@@ -16,114 +21,148 @@ const expect = std.testing.expect;
16pub fn sin(x: anytype) @TypeOf(x) {21pub fn sin(x: anytype) @TypeOf(x) {
17 const T = @TypeOf(x);22 const T = @TypeOf(x);
18 return switch (T) {23 return switch (T) {
19 f32 => sin_(T, x),24 f32 => sin32(x),
20 f64 => sin_(T, x),25 f64 => sin64(x),
21 else => @compileError("sin not implemented for " ++ @typeName(T)),26 else => @compileError("sin not implemented for " ++ @typeName(T)),
22 };27 };
23}28}
2429
25// sin polynomial coefficients30fn sin32(x: f32) f32 {
26const S0 = 1.58962301576546568060E-10;31 // Small multiples of pi/2 rounded to double precision.
27const S1 = -2.50507477628578072866E-8;32 const s1pio2: f64 = 1.0 * math.pi / 2.0; // 0x3FF921FB, 0x54442D18
28const S2 = 2.75573136213857245213E-6;33 const s2pio2: f64 = 2.0 * math.pi / 2.0; // 0x400921FB, 0x54442D18
29const S3 = -1.98412698295895385996E-4;34 const s3pio2: f64 = 3.0 * math.pi / 2.0; // 0x4012D97C, 0x7F3321D2
30const S4 = 8.33333333332211858878E-3;35 const s4pio2: f64 = 4.0 * math.pi / 2.0; // 0x401921FB, 0x54442D18
31const S5 = -1.66666666666666307295E-1;36
3237 var ix = @bitCast(u32, x);
33// cos polynomial coeffiecients38 const sign = ix >> 31 != 0;
34const C0 = -1.13585365213876817300E-11;39 ix &= 0x7fffffff;
35const C1 = 2.08757008419747316778E-9;40
36const C2 = -2.75573141792967388112E-7;41 if (ix <= 0x3f490fda) { // |x| ~<= pi/4
37const C3 = 2.48015872888517045348E-5;42 if (ix < 0x39800000) { // |x| < 2**-12
38const C4 = -1.38888888888730564116E-3;43 // raise inexact if x!=0 and underflow if subnormal
39const C5 = 4.16666666666665929218E-2;44 math.doNotOptimizeAway(if (ix < 0x00800000) x / 0x1p120 else x + 0x1p120);
4045 return x;
41const pi4a = 7.85398125648498535156e-1;46 }
42const pi4b = 3.77489470793079817668E-8;47 return kernel.__sindf(x);
43const pi4c = 2.69515142907905952645E-15;48 }
44const m4pi = 1.273239544735162542821171882678754627704620361328125;49 if (ix <= 0x407b53d1) { // |x| ~<= 5*pi/4
4550 if (ix <= 0x4016cbe3) { // |x| ~<= 3pi/4
46fn sin_(comptime T: type, x_: T) T {51 if (sign) {
47 const I = std.meta.Int(.signed, @typeInfo(T).Float.bits);52 return -kernel.__cosdf(x + s1pio2);
4853 } else {
49 var x = x_;54 return kernel.__cosdf(x - s1pio2);
50 if (x == 0 or math.isNan(x)) {55 }
51 return x;56 }
57 return kernel.__sindf(if (sign) -(x + s2pio2) else -(x - s2pio2));
52 }58 }
53 if (math.isInf(x)) {59 if (ix <= 0x40e231d5) { // |x| ~<= 9*pi/4
54 return math.nan(T);60 if (ix <= 0x40afeddf) { // |x| ~<= 7*pi/4
61 if (sign) {
62 return kernel.__cosdf(x + s3pio2);
63 } else {
64 return -kernel.__cosdf(x - s3pio2);
65 }
66 }
67 return kernel.__sindf(if (sign) x + s4pio2 else x - s4pio2);
55 }68 }
5669
57 var sign = x < 0;70 // sin(Inf or NaN) is NaN
58 x = math.fabs(x);71 if (ix >= 0x7f800000) {
72 return x - x;
73 }
5974
60 var y = math.floor(x * m4pi);75 var y: f64 = undefined;
61 var j = @floatToInt(I, y);76 const n = __rem_pio2f(x, &y);
77 return switch (n & 3) {
78 0 => kernel.__sindf(y),
79 1 => kernel.__cosdf(y),
80 2 => kernel.__sindf(-y),
81 else => -kernel.__cosdf(y),
82 };
83}
6284
63 if (j & 1 == 1) {85fn sin64(x: f64) f64 {
64 j += 1;86 var ix = @bitCast(u64, x) >> 32;
65 y += 1;87 ix &= 0x7fffffff;
88
89 // |x| ~< pi/4
90 if (ix <= 0x3fe921fb) {
91 if (ix < 0x3e500000) { // |x| < 2**-26
92 // raise inexact if x != 0 and underflow if subnormal
93 math.doNotOptimizeAway(if (ix < 0x00100000) x / 0x1p120 else x + 0x1p120);
94 return x;
95 }
96 return kernel.__sin(x, 0.0, 0);
66 }97 }
6798
68 j &= 7;99 // sin(Inf or NaN) is NaN
69 if (j > 3) {100 if (ix >= 0x7ff00000) {
70 j -= 4;101 return x - x;
71 sign = !sign;
72 }102 }
73103
74 const z = ((x - y * pi4a) - y * pi4b) - y * pi4c;104 var y: [2]f64 = undefined;
75 const w = z * z;105 const n = __rem_pio2(x, &y);
76106 return switch (n & 3) {
77 const r = if (j == 1 or j == 2)107 0 => kernel.__sin(y[0], y[1], 1),
78 1.0 - 0.5 * w + w * w * (C5 + w * (C4 + w * (C3 + w * (C2 + w * (C1 + w * C0)))))108 1 => kernel.__cos(y[0], y[1]),
79 else109 2 => -kernel.__sin(y[0], y[1], 1),
80 z + z * w * (S5 + w * (S4 + w * (S3 + w * (S2 + w * (S1 + w * S0)))));110 else => -kernel.__cos(y[0], y[1]),
81111 };
82 return if (sign) -r else r;
83}112}
84113
85test "math.sin" {114test "math.sin" {
86 try expect(sin(@as(f32, 0.0)) == sin_(f32, 0.0));115 try expect(sin(@as(f32, 0.0)) == sin32(0.0));
87 try expect(sin(@as(f64, 0.0)) == sin_(f64, 0.0));116 try expect(sin(@as(f64, 0.0)) == sin64(0.0));
88 try expect(comptime (math.sin(@as(f64, 2))) == math.sin(@as(f64, 2)));117 try expect(comptime (math.sin(@as(f64, 2))) == math.sin(@as(f64, 2)));
89}118}
90119
91test "math.sin32" {120test "math.sin32" {
92 const epsilon = 0.000001;121 const epsilon = 0.00001;
93122
94 try expect(math.approxEqAbs(f32, sin_(f32, 0.0), 0.0, epsilon));123 try expect(math.approxEqAbs(f32, sin32(0.0), 0.0, epsilon));
95 try expect(math.approxEqAbs(f32, sin_(f32, 0.2), 0.198669, epsilon));124 try expect(math.approxEqAbs(f32, sin32(0.2), 0.198669, epsilon));
96 try expect(math.approxEqAbs(f32, sin_(f32, 0.8923), 0.778517, epsilon));125 try expect(math.approxEqAbs(f32, sin32(0.8923), 0.778517, epsilon));
97 try expect(math.approxEqAbs(f32, sin_(f32, 1.5), 0.997495, epsilon));126 try expect(math.approxEqAbs(f32, sin32(1.5), 0.997495, epsilon));
98 try expect(math.approxEqAbs(f32, sin_(f32, -1.5), -0.997495, epsilon));127 try expect(math.approxEqAbs(f32, sin32(-1.5), -0.997495, epsilon));
99 try expect(math.approxEqAbs(f32, sin_(f32, 37.45), -0.246544, epsilon));128 try expect(math.approxEqAbs(f32, sin32(37.45), -0.246544, epsilon));
100 try expect(math.approxEqAbs(f32, sin_(f32, 89.123), 0.916166, epsilon));129 try expect(math.approxEqAbs(f32, sin32(89.123), 0.916166, epsilon));
101}130}
102131
103test "math.sin64" {132test "math.sin64" {
104 const epsilon = 0.000001;133 const epsilon = 0.000001;
105134
106 try expect(math.approxEqAbs(f64, sin_(f64, 0.0), 0.0, epsilon));135 try expect(math.approxEqAbs(f64, sin64(0.0), 0.0, epsilon));
107 try expect(math.approxEqAbs(f64, sin_(f64, 0.2), 0.198669, epsilon));136 try expect(math.approxEqAbs(f64, sin64(0.2), 0.198669, epsilon));
108 try expect(math.approxEqAbs(f64, sin_(f64, 0.8923), 0.778517, epsilon));137 try expect(math.approxEqAbs(f64, sin64(0.8923), 0.778517, epsilon));
109 try expect(math.approxEqAbs(f64, sin_(f64, 1.5), 0.997495, epsilon));138 try expect(math.approxEqAbs(f64, sin64(1.5), 0.997495, epsilon));
110 try expect(math.approxEqAbs(f64, sin_(f64, -1.5), -0.997495, epsilon));139 try expect(math.approxEqAbs(f64, sin64(-1.5), -0.997495, epsilon));
111 try expect(math.approxEqAbs(f64, sin_(f64, 37.45), -0.246543, epsilon));140 try expect(math.approxEqAbs(f64, sin64(37.45), -0.246543, epsilon));
112 try expect(math.approxEqAbs(f64, sin_(f64, 89.123), 0.916166, epsilon));141 try expect(math.approxEqAbs(f64, sin64(89.123), 0.916166, epsilon));
113}142}
114143
115test "math.sin32.special" {144test "math.sin32.special" {
116 try expect(sin_(f32, 0.0) == 0.0);145 try expect(sin32(0.0) == 0.0);
117 try expect(sin_(f32, -0.0) == -0.0);146 try expect(sin32(-0.0) == -0.0);
118 try expect(math.isNan(sin_(f32, math.inf(f32))));147 try expect(math.isNan(sin32(math.inf(f32))));
119 try expect(math.isNan(sin_(f32, -math.inf(f32))));148 try expect(math.isNan(sin32(-math.inf(f32))));
120 try expect(math.isNan(sin_(f32, math.nan(f32))));149 try expect(math.isNan(sin32(math.nan(f32))));
121}150}
122151
123test "math.sin64.special" {152test "math.sin64.special" {
124 try expect(sin_(f64, 0.0) == 0.0);153 try expect(sin64(0.0) == 0.0);
125 try expect(sin_(f64, -0.0) == -0.0);154 try expect(sin64(-0.0) == -0.0);
126 try expect(math.isNan(sin_(f64, math.inf(f64))));155 try expect(math.isNan(sin64(math.inf(f64))));
127 try expect(math.isNan(sin_(f64, -math.inf(f64))));156 try expect(math.isNan(sin64(-math.inf(f64))));
128 try expect(math.isNan(sin_(f64, math.nan(f64))));157 try expect(math.isNan(sin64(math.nan(f64))));
158}
159
160test "math.sin32 #9901" {
161 const float = @bitCast(f32, @as(u32, 0b11100011111111110000000000000000));
162 _ = std.math.sin(float);
163}
164
165test "math.sin64 #9901" {
166 const float = @bitCast(f64, @as(u64, 0b1111111101000001000000001111110111111111100000000000000000000001));
167 _ = std.math.sin(float);
129}168}
lib/std/math/tan.zig+93-70
...@@ -1,12 +1,18 @@...@@ -1,12 +1,18 @@
1// Ported from go, which is licensed under a BSD-3 license.1// Ported from musl, which is licensed under the MIT license:
2// https://golang.org/LICENSE2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
3//3//
4// https://git.musl-libc.org/cgit/musl/tree/src/math/tanf.c
5// https://git.musl-libc.org/cgit/musl/tree/src/math/tan.c
4// https://golang.org/src/math/tan.go6// https://golang.org/src/math/tan.go
57
6const std = @import("../std.zig");8const std = @import("../std.zig");
7const math = std.math;9const math = std.math;
8const expect = std.testing.expect;10const expect = std.testing.expect;
911
12const kernel = @import("__trig.zig");
13const __rem_pio2 = @import("__rem_pio2.zig").__rem_pio2;
14const __rem_pio2f = @import("__rem_pio2f.zig").__rem_pio2f;
15
10/// Returns the tangent of the radian value x.16/// Returns the tangent of the radian value x.
11///17///
12/// Special Cases:18/// Special Cases:
...@@ -16,102 +22,119 @@ const expect = std.testing.expect;...@@ -16,102 +22,119 @@ const expect = std.testing.expect;
16pub fn tan(x: anytype) @TypeOf(x) {22pub fn tan(x: anytype) @TypeOf(x) {
17 const T = @TypeOf(x);23 const T = @TypeOf(x);
18 return switch (T) {24 return switch (T) {
19 f32 => tan_(f32, x),25 f32 => tan32(x),
20 f64 => tan_(f64, x),26 f64 => tan64(x),
21 else => @compileError("tan not implemented for " ++ @typeName(T)),27 else => @compileError("tan not implemented for " ++ @typeName(T)),
22 };28 };
23}29}
2430
25const Tp0 = -1.30936939181383777646E4;31fn tan32(x: f32) f32 {
26const Tp1 = 1.15351664838587416140E6;32 // Small multiples of pi/2 rounded to double precision.
27const Tp2 = -1.79565251976484877988E7;33 const t1pio2: f64 = 1.0 * math.pi / 2.0; // 0x3FF921FB, 0x54442D18
2834 const t2pio2: f64 = 2.0 * math.pi / 2.0; // 0x400921FB, 0x54442D18
29const Tq1 = 1.36812963470692954678E4;35 const t3pio2: f64 = 3.0 * math.pi / 2.0; // 0x4012D97C, 0x7F3321D2
30const Tq2 = -1.32089234440210967447E6;36 const t4pio2: f64 = 4.0 * math.pi / 2.0; // 0x401921FB, 0x54442D18
31const Tq3 = 2.50083801823357915839E7;37
32const Tq4 = -5.38695755929454629881E7;38 var ix = @bitCast(u32, x);
3339 const sign = ix >> 31 != 0;
34const pi4a = 7.85398125648498535156e-1;40 ix &= 0x7fffffff;
35const pi4b = 3.77489470793079817668E-8;41
36const pi4c = 2.69515142907905952645E-15;42 if (ix <= 0x3f490fda) { // |x| ~<= pi/4
37const m4pi = 1.273239544735162542821171882678754627704620361328125;43 if (ix < 0x39800000) { // |x| < 2**-12
3844 // raise inexact if x!=0 and underflow if subnormal
39fn tan_(comptime T: type, x_: T) T {45 math.doNotOptimizeAway(if (ix < 0x00800000) x / 0x1p120 else x + 0x1p120);
40 const I = std.meta.Int(.signed, @typeInfo(T).Float.bits);46 return x;
4147 }
42 var x = x_;48 return kernel.__tandf(x, false);
43 if (x == 0 or math.isNan(x)) {
44 return x;
45 }49 }
46 if (math.isInf(x)) {50 if (ix <= 0x407b53d1) { // |x| ~<= 5*pi/4
47 return math.nan(T);51 if (ix <= 0x4016cbe3) { // |x| ~<= 3pi/4
52 return kernel.__tandf((if (sign) x + t1pio2 else x - t1pio2), true);
53 } else {
54 return kernel.__tandf((if (sign) x + t2pio2 else x - t2pio2), false);
55 }
56 }
57 if (ix <= 0x40e231d5) { // |x| ~<= 9*pi/4
58 if (ix <= 0x40afeddf) { // |x| ~<= 7*pi/4
59 return kernel.__tandf((if (sign) x + t3pio2 else x - t3pio2), true);
60 } else {
61 return kernel.__tandf((if (sign) x + t4pio2 else x - t4pio2), false);
62 }
48 }63 }
4964
50 var sign = x < 0;65 // tan(Inf or NaN) is NaN
51 x = math.fabs(x);66 if (ix >= 0x7f800000) {
5267 return x - x;
53 var y = math.floor(x * m4pi);
54 var j = @floatToInt(I, y);
55
56 if (j & 1 == 1) {
57 j += 1;
58 y += 1;
59 }68 }
6069
61 const z = ((x - y * pi4a) - y * pi4b) - y * pi4c;70 var y: f64 = undefined;
62 const w = z * z;71 const n = __rem_pio2f(x, &y);
72 return kernel.__tandf(y, n & 1 != 0);
73}
6374
64 var r = if (w > 1e-14)75fn tan64(x: f64) f64 {
65 z + z * (w * ((Tp0 * w + Tp1) * w + Tp2) / ((((w + Tq1) * w + Tq2) * w + Tq3) * w + Tq4))76 var ix = @bitCast(u64, x) >> 32;
66 else77 ix &= 0x7fffffff;
67 z;78
79 // |x| ~< pi/4
80 if (ix <= 0x3fe921fb) {
81 if (ix < 0x3e400000) { // |x| < 2**-27
82 // raise inexact if x!=0 and underflow if subnormal
83 math.doNotOptimizeAway(if (ix < 0x00100000) x / 0x1p120 else x + 0x1p120);
84 return x;
85 }
86 return kernel.__tan(x, 0.0, false);
87 }
6888
69 if (j & 2 == 2) {89 // tan(Inf or NaN) is NaN
70 r = -1 / r;90 if (ix >= 0x7ff00000) {
91 return x - x;
71 }92 }
7293
73 return if (sign) -r else r;94 var y: [2]f64 = undefined;
95 const n = __rem_pio2(x, &y);
96 return kernel.__tan(y[0], y[1], n & 1 != 0);
74}97}
7598
76test "math.tan" {99test "math.tan" {
77 try expect(tan(@as(f32, 0.0)) == tan_(f32, 0.0));100 try expect(tan(@as(f32, 0.0)) == tan32(0.0));
78 try expect(tan(@as(f64, 0.0)) == tan_(f64, 0.0));101 try expect(tan(@as(f64, 0.0)) == tan64(0.0));
79}102}
80103
81test "math.tan32" {104test "math.tan32" {
82 const epsilon = 0.000001;105 const epsilon = 0.00001;
83106
84 try expect(math.approxEqAbs(f32, tan_(f32, 0.0), 0.0, epsilon));107 try expect(math.approxEqAbs(f32, tan32(0.0), 0.0, epsilon));
85 try expect(math.approxEqAbs(f32, tan_(f32, 0.2), 0.202710, epsilon));108 try expect(math.approxEqAbs(f32, tan32(0.2), 0.202710, epsilon));
86 try expect(math.approxEqAbs(f32, tan_(f32, 0.8923), 1.240422, epsilon));109 try expect(math.approxEqAbs(f32, tan32(0.8923), 1.240422, epsilon));
87 try expect(math.approxEqAbs(f32, tan_(f32, 1.5), 14.101420, epsilon));110 try expect(math.approxEqAbs(f32, tan32(1.5), 14.101420, epsilon));
88 try expect(math.approxEqAbs(f32, tan_(f32, 37.45), -0.254397, epsilon));111 try expect(math.approxEqAbs(f32, tan32(37.45), -0.254397, epsilon));
89 try expect(math.approxEqAbs(f32, tan_(f32, 89.123), 2.285852, epsilon));112 try expect(math.approxEqAbs(f32, tan32(89.123), 2.285852, epsilon));
90}113}
91114
92test "math.tan64" {115test "math.tan64" {
93 const epsilon = 0.000001;116 const epsilon = 0.000001;
94117
95 try expect(math.approxEqAbs(f64, tan_(f64, 0.0), 0.0, epsilon));118 try expect(math.approxEqAbs(f64, tan64(0.0), 0.0, epsilon));
96 try expect(math.approxEqAbs(f64, tan_(f64, 0.2), 0.202710, epsilon));119 try expect(math.approxEqAbs(f64, tan64(0.2), 0.202710, epsilon));
97 try expect(math.approxEqAbs(f64, tan_(f64, 0.8923), 1.240422, epsilon));120 try expect(math.approxEqAbs(f64, tan64(0.8923), 1.240422, epsilon));
98 try expect(math.approxEqAbs(f64, tan_(f64, 1.5), 14.101420, epsilon));121 try expect(math.approxEqAbs(f64, tan64(1.5), 14.101420, epsilon));
99 try expect(math.approxEqAbs(f64, tan_(f64, 37.45), -0.254397, epsilon));122 try expect(math.approxEqAbs(f64, tan64(37.45), -0.254397, epsilon));
100 try expect(math.approxEqAbs(f64, tan_(f64, 89.123), 2.2858376, epsilon));123 try expect(math.approxEqAbs(f64, tan64(89.123), 2.2858376, epsilon));
101}124}
102125
103test "math.tan32.special" {126test "math.tan32.special" {
104 try expect(tan_(f32, 0.0) == 0.0);127 try expect(tan32(0.0) == 0.0);
105 try expect(tan_(f32, -0.0) == -0.0);128 try expect(tan32(-0.0) == -0.0);
106 try expect(math.isNan(tan_(f32, math.inf(f32))));129 try expect(math.isNan(tan32(math.inf(f32))));
107 try expect(math.isNan(tan_(f32, -math.inf(f32))));130 try expect(math.isNan(tan32(-math.inf(f32))));
108 try expect(math.isNan(tan_(f32, math.nan(f32))));131 try expect(math.isNan(tan32(math.nan(f32))));
109}132}
110133
111test "math.tan64.special" {134test "math.tan64.special" {
112 try expect(tan_(f64, 0.0) == 0.0);135 try expect(tan64(0.0) == 0.0);
113 try expect(tan_(f64, -0.0) == -0.0);136 try expect(tan64(-0.0) == -0.0);
114 try expect(math.isNan(tan_(f64, math.inf(f64))));137 try expect(math.isNan(tan64(math.inf(f64))));
115 try expect(math.isNan(tan_(f64, -math.inf(f64))));138 try expect(math.isNan(tan64(-math.inf(f64))));
116 try expect(math.isNan(tan_(f64, math.nan(f64))));139 try expect(math.isNan(tan64(math.nan(f64))));
117}140}