Commit 4ccac1de416bff8846352ef3bf2fb592c5a419a9 | zig

author	mihael <hi@mihaelm.com> 2026-04-02 23:53:13+02:00
committer	mihael <hi@mihaelm.com> 2026-04-02 23:54:20+02:00
log	4ccac1de416bff8846352ef3bf2fb592c5a419a9
tree	989d270f4df185757eea5cd3c31cf6479c55cc0c
parent	ad10c7600793547393848528e7d4a608d9fde59e
signature	Signed by SSH key SHA256:aoFoShdYLdrqMichqKXSSieTKUfACUIDJHsKc4V2tQg

`compiler_rt`: Make `long double` trig impls non-generic

5 files changed, 336 insertions(+), 236 deletions(-)

lib
- compiler_rt

lib/compiler_rt/cos.zig+30-14

 @@ -123,36 +123,52 @@ pub fn cos(x: f64) callconv(.c) f64 {
     };
 }
-fn coslGeneric(comptime T: type, x: T) T {
+pub fn cosx(x: f80) callconv(.c) f80 {
     const se = ld.signExponent(x) & 0x7fff;
     if (se == 0x7fff) {
         return x - x;
     }
     if (@abs(x) < trig.pi_4) {
-        if (se < 0x3fff - math.floatMantissaBits(T)) {
+        if (se < 0x3fff - math.floatMantissaBits(f80)) {
             // raise inexact if x!=0
             return 1.0 + x;
         }
-        return trig.cosl(T, x, 0.0);
+        return trig.cosx(x, 0.0);
     }
-    var y: [2]T = undefined;
+    var y: [2]f80 = undefined;
-    const n = rem_pio2l(T, x, &y);
+    const n = rem_pio2l(f80, x, &y);
     return switch (n & 3) {
-        0 => trig.cosl(T, y[0], y[1]),
+        0 => trig.cosx(y[0], y[1]),
-        1 => -trig.sinl(T, y[0], y[1], 1),
+        1 => -trig.sinx(y[0], y[1], 1),
-        2 => -trig.cosl(T, y[0], y[1]),
+        2 => -trig.cosx(y[0], y[1]),
-        else => trig.sinl(T, y[0], y[1], 1),
+        else => trig.sinx(y[0], y[1], 1),
     };
 }
-pub fn cosx(x: f80) callconv(.c) f80 {
-    return coslGeneric(f80, x);
-}
 pub fn cosq(x: f128) callconv(.c) f128 {
-    return coslGeneric(f128, x);
+    const se = ld.signExponent(x) & 0x7fff;
+    if (se == 0x7fff) {
+        return x - x;
+    }
+    if (@abs(x) < trig.pi_4) {
+        if (se < 0x3fff - math.floatMantissaBits(f128)) {
+            // raise inexact if x!=0
+            return 1.0 + x;
+        }
+        return trig.cosq(x, 0.0);
+    }
+    var y: [2]f128 = undefined;
+    const n = rem_pio2l(f128, x, &y);
+    return switch (n & 3) {
+        0 => trig.cosq(y[0], y[1]),
+        1 => -trig.sinq(y[0], y[1], 1),
+        2 => -trig.cosq(y[0], y[1]),
+        else => trig.sinq(y[0], y[1], 1),
+    };
 }
 pub fn cosl(x: c_longdouble) callconv(.c) c_longdouble {

lib/compiler_rt/sin.zig+33-14

 @@ -133,39 +133,58 @@ pub fn sin(x: f64) callconv(.c) f64 {
     };
 }
-fn sinlGeneric(comptime T: type, x: T) T {
+fn sinx(x: f80) callconv(.c) f80 {
     const se = ld.signExponent(x) & 0x7fff;
     if (se == 0x7fff) {
         return x - x;
     }
     if (@abs(x) < trig.pi_4) {
-        if (se < 0x3fff - (math.floatMantissaBits(T) / 2)) {
+        if (se < 0x3fff - (math.floatMantissaBits(f80) / 2)) {
             // raise inexact if x!=0 and underflow if subnormal
             if (compiler_rt.want_float_exceptions) {
                 mem.doNotOptimizeAway(if (se == 0) x * 0x1p-120 else x + 0x1p120);
             }
             return x;
         }
-        return trig.sinl(T, x, 0.0, 0);
+        return trig.sinx(x, 0.0, 0);
     }
-    var y: [2]T = undefined;
+    var y: [2]f80 = undefined;
-    const n = rem_pio2l(T, x, &y);
+    const n = rem_pio2l(f80, x, &y);
     return switch (n & 3) {
-        0 => trig.sinl(T, y[0], y[1], 1),
+        0 => trig.sinx(y[0], y[1], 1),
-        1 => trig.cosl(T, y[0], y[1]),
+        1 => trig.cosx(y[0], y[1]),
-        2 => -trig.sinl(T, y[0], y[1], 1),
+        2 => -trig.sinx(y[0], y[1], 1),
-        else => -trig.cosl(T, y[0], y[1]),
+        else => -trig.cosx(y[0], y[1]),
     };
 }
-pub fn sinx(x: f80) callconv(.c) f80 {
-    return sinlGeneric(f80, x);
-}
 pub fn sinq(x: f128) callconv(.c) f128 {
-    return sinlGeneric(f128, x);
+    const se = ld.signExponent(x) & 0x7fff;
+    if (se == 0x7fff) {
+        return x - x;
+    }
+    if (@abs(x) < trig.pi_4) {
+        if (se < 0x3fff - (math.floatMantissaBits(f128) / 2)) {
+            // raise inexact if x!=0 and underflow if subnormal
+            if (compiler_rt.want_float_exceptions) {
+                mem.doNotOptimizeAway(if (se == 0) x * 0x1p-120 else x + 0x1p120);
+            }
+            return x;
+        }
+        return trig.sinq(x, 0.0, 0);
+    }
+    var y: [2]f128 = undefined;
+    const n = rem_pio2l(f128, x, &y);
+    return switch (n & 3) {
+        0 => trig.sinq(y[0], y[1], 1),
+        1 => trig.cosq(y[0], y[1]),
+        2 => -trig.sinq(y[0], y[1], 1),
+        else => -trig.cosq(y[0], y[1]),
+    };
 }
 pub fn sinl(x: c_longdouble) callconv(.c) c_longdouble {

lib/compiler_rt/sincos.zig+54-17

...	@@ -192,11 +192,7 @@ pub fn sincos(x: f64, r_sin: f64, r_cos: f64) callconv(.c) void {	...	@@ -192,11 +192,7 @@ pub fn sincos(x: f64, r_sin: f64, r_cos: f64) callconv(.c) void {
192	}	192	}
193	}	193	}
194		194
195	fn sincoslGeneric(comptime T: type, x: T, r_sin: T, r_cos: T) void {	195	pub fn sincosx(x: f80, r_sin: f80, r_cos: f80) callconv(.c) void {
196	if (T != f80 and T != f128) {
197	@compileError("`sincoslGeneric` implemented only for `f80` and `f128`, got: " ++ @typeName(T));
198	}
199
200	const se = ld.signExponent(x) & 0x7fff;	196	const se = ld.signExponent(x) & 0x7fff;
201	if (se == 0x7fff) {	197	if (se == 0x7fff) {
202	const result = x - x;	198	const result = x - x;
...	@@ -206,7 +202,7 @@ fn sincoslGeneric(comptime T: type, x: T, r_sin: T, r_cos: T) void {	...	@@ -206,7 +202,7 @@ fn sincoslGeneric(comptime T: type, x: T, r_sin: T, r_cos: T) void {
206	}	202	}
207		203
208	if (@abs(x) < trig.pi_4) {	204	if (@abs(x) < trig.pi_4) {
209	if (se < 0x3fff - math.floatMantissaBits(T)) {	205	if (se < 0x3fff - math.floatMantissaBits(f80)) {
210	// raise underflow if subnormal	206	// raise underflow if subnormal
211	if (compiler_rt.want_float_exceptions and se == 0) {	207	if (compiler_rt.want_float_exceptions and se == 0) {
212	mem.doNotOptimizeAway(x * 0x1p-120);	208	mem.doNotOptimizeAway(x * 0x1p-120);
...	@@ -216,15 +212,15 @@ fn sincoslGeneric(comptime T: type, x: T, r_sin: T, r_cos: T) void {	...	@@ -216,15 +212,15 @@ fn sincoslGeneric(comptime T: type, x: T, r_sin: T, r_cos: T) void {
216	r_cos.* = 1.0 + x;	212	r_cos.* = 1.0 + x;
217	return;	213	return;
218	}	214	}
219	r_sin.* = trig.sinl(T, x, 0.0, 0);	215	r_sin.* = trig.sinx(x, 0.0, 0);
220	r_cos.* = trig.cosl(T, x, 0.0);	216	r_cos.* = trig.cosx(x, 0.0);
221	return;	217	return;
222	}	218	}
223		219
224	var y: [2]T = undefined;	220	var y: [2]f80 = undefined;
225	const n = rem_pio2l(T, x, &y);	221	const n = rem_pio2l(f80, x, &y);
226	const s = trig.sinl(T, y[0], y[1], 1);	222	const s = trig.sinx(y[0], y[1], 1);
227	const c = trig.cosl(T, y[0], y[1]);	223	const c = trig.cosx(y[0], y[1]);
228	switch (n & 3) {	224	switch (n & 3) {
229	0 => {	225	0 => {
230	r_sin.* = s;	226	r_sin.* = s;
...	@@ -245,12 +241,53 @@ fn sincoslGeneric(comptime T: type, x: T, r_sin: T, r_cos: T) void {	...	@@ -245,12 +241,53 @@ fn sincoslGeneric(comptime T: type, x: T, r_sin: T, r_cos: T) void {
245	}	241	}
246	}	242	}
247		243
248	pub fn sincosx(x: f80, r_sin: f80, r_cos: f80) callconv(.c) void {
249	return sincoslGeneric(f80, x, r_sin, r_cos);
250	}
251
252	pub fn sincosq(x: f128, r_sin: f128, r_cos: f128) callconv(.c) void {	244	pub fn sincosq(x: f128, r_sin: f128, r_cos: f128) callconv(.c) void {
253	return sincoslGeneric(f128, x, r_sin, r_cos);	245	const se = ld.signExponent(x) & 0x7fff;
		246	if (se == 0x7fff) {
		247	const result = x - x;
		248	r_sin.* = result;
		249	r_cos.* = result;
		250	return;
		251	}
		252
		253	if (@abs(x) < trig.pi_4) {
		254	if (se < 0x3fff - math.floatMantissaBits(f128)) {
		255	// raise underflow if subnormal
		256	if (compiler_rt.want_float_exceptions and se == 0) {
		257	mem.doNotOptimizeAway(x * 0x1p-120);
		258	}
		259	r_sin.* = x;
		260	// raise inexact if x!=0
		261	r_cos.* = 1.0 + x;
		262	return;
		263	}
		264	r_sin.* = trig.sinq(x, 0.0, 0);
		265	r_cos.* = trig.cosq(x, 0.0);
		266	return;
		267	}
		268
		269	var y: [2]f128 = undefined;
		270	const n = rem_pio2l(f128, x, &y);
		271	const s = trig.sinq(y[0], y[1], 1);
		272	const c = trig.cosq(y[0], y[1]);
		273	switch (n & 3) {
		274	0 => {
		275	r_sin.* = s;
		276	r_cos.* = c;
		277	},
		278	1 => {
		279	r_sin.* = c;
		280	r_cos.* = -s;
		281	},
		282	2 => {
		283	r_sin.* = -s;
		284	r_cos.* = -c;
		285	},
		286	else => {
		287	r_sin.* = -c;
		288	r_cos.* = s;
		289	},
		290	}
254	}	291	}
255		292
256	pub fn sincosl(x: c_longdouble, r_sin: c_longdouble, r_cos: c_longdouble) callconv(.c) void {	293	pub fn sincosl(x: c_longdouble, r_sin: c_longdouble, r_cos: c_longdouble) callconv(.c) void {

lib/compiler_rt/tan.zig+24-15

 @@ -120,37 +120,46 @@ pub fn tan(x: f64) callconv(.c) f64 {
     return kernel.tan(y[0], y[1], n & 1 != 0);
 }
-fn tanlGeneric(comptime T: type, x: T) T {
+pub fn tanx(x: f80) callconv(.c) f80 {
-    if (!(T == f80 or T == f128)) {
-        @compileError("`tanlGeneric` implemented only for `f80` and `f128`, got: " ++ T);
-    }
     const se = ld.signExponent(x) & 0x7fff;
     if (se == 0x7fff) {
         return x - x;
     }
     if (@abs(x) < kernel.pi_4) {
-        if (se < 0x3fff - math.floatMantissaBits(T) / 2) {
+        if (se < 0x3fff - math.floatMantissaBits(f80) / 2) {
             if (compiler_rt.want_float_exceptions) {
                 mem.doNotOptimizeAway(if (se == 0) x * 0x1p-120 else x + 0x1p120);
             }
             return x;
         }
-        return kernel.tanl(T, x, 0.0, 0);
+        return kernel.tanx(x, 0.0, 0);
     }
-    var y: [2]T = undefined;
+    var y: [2]f80 = undefined;
-    const n = rem_pio2l(T, x, &y);
+    const n = rem_pio2l(f80, x, &y);
-    return kernel.tanl(T, y[0], y[1], n & 1);
+    return kernel.tanx(y[0], y[1], n & 1);
-}
-pub fn tanx(x: f80) callconv(.c) f80 {
-    return tanlGeneric(f80, x);
 }
 pub fn tanq(x: f128) callconv(.c) f128 {
-    return tanlGeneric(f128, x);
+    const se = ld.signExponent(x) & 0x7fff;
+    if (se == 0x7fff) {
+        return x - x;
+    }
+    if (@abs(x) < kernel.pi_4) {
+        if (se < 0x3fff - math.floatMantissaBits(f128) / 2) {
+            if (compiler_rt.want_float_exceptions) {
+                mem.doNotOptimizeAway(if (se == 0) x * 0x1p-120 else x + 0x1p120);
+            }
+            return x;
+        }
+        return kernel.tanq(x, 0.0, 0);
+    }
+    var y: [2]f128 = undefined;
+    const n = rem_pio2l(f128, x, &y);
+    return kernel.tanq(y[0], y[1], n & 1);
 }
 pub fn tanl(x: c_longdouble) callconv(.c) c_longdouble {

lib/compiler_rt/trig.zig+195-176

...	@@ -80,49 +80,43 @@ pub fn cosdf(x: f64) f32 {	...	@@ -80,49 +80,43 @@ pub fn cosdf(x: f64) f32 {
80	return @floatCast(((1.0 + z * C0) + w * C1) + (w * z) * r);	80	return @floatCast(((1.0 + z * C0) + w * C1) + (w * z) * r);
81	}	81	}
82		82
83	pub fn cosl(comptime T: type, x: T, y: T) T {	83	pub fn cosx(x: f80, y: f80) f80 {
84	const impl = switch (T) {	84	const C1: f80 = 0.0416666666666666666136;
85	f80 => struct {	85	const C2: f64 = -0.0013888888888888874;
86	const C1: T = 0.0416666666666666666136;	86	const C3: f64 = 0.000024801587301571716;
87		87	const C4: f64 = -0.00000027557319215507120;
88	const C2: f64 = -0.0013888888888888874;	88	const C5: f64 = 0.0000000020876754400407278;
89	const C3: f64 = 0.000024801587301571716;	89	const C6: f64 = -1.1470297442401303e-11;
90	const C4: f64 = -0.00000027557319215507120;	90	const C7: f64 = 4.7383039476436467e-14;
91	const C5: f64 = 0.0000000020876754400407278;
92	const C6: f64 = -1.1470297442401303e-11;
93	const C7: f64 = 4.7383039476436467e-14;
94
95	inline fn poly(z: T) T {
96	return z * (C1 + z * (C2 + z * (C3 + z * (C4 +
97	z * (C5 + z * (C6 + z * C7))))));
98	}
99	},
100	f128 => struct {
101	const C1: T = 0.04166666666666666666666666666666658424671;
102	const C2: T = -0.001388888888888888888888888888863490893732;
103	const C3: T = 0.00002480158730158730158730158600795304914210;
104	const C4: T = -0.2755731922398589065255474947078934284324e-6;
105	const C5: T = 0.2087675698786809897659225313136400793948e-8;
106	const C6: T = -0.1147074559772972315817149986812031204775e-10;
107	const C7: T = 0.4779477332386808976875457937252120293400e-13;
108
109	const C8: f64 = -0.1561920696721507929516718307820958119868e-15;
110	const C9: f64 = 0.4110317413744594971475941557607804508039e-18;
111	const C10: f64 = -0.8896592467191938803288521958313920156409e-21;
112	const C11: f64 = 0.1601061435794535138244346256065192782581e-23;
113
114	inline fn poly(z: T) T {
115	return z * (C1 + z * (C2 + z * (C3 + z * (C4 + z * (C5 + z * (C6 +
116	z * (C7 + z * (C8 + z * (C9 + z * (C10 + z * C11))))))))));
117	}
118	},
119	else => @compileError("cosl supports only f80 and f128, got: " ++ @typeName(T)),
120	};
121		91
122	const z = x * x;	92	const z = x * x;
123	const r = impl.poly(z);	93	const r = z * (C1 + z * (C2 + z * (C3 + z * (C4 +
		94	z * (C5 + z * (C6 + z * C7))))));
124	const hz = 0.5 * z;	95	const hz = 0.5 * z;
125	const w = 1.0 - hz;	96	const w = 1.0 - hz;
		97
		98	return w + (((1.0 - w) - hz) + (z * r - x * y));
		99	}
		100
		101	pub fn cosq(x: f128, y: f128) f128 {
		102	const C1: f128 = 0.04166666666666666666666666666666658424671;
		103	const C2: f128 = -0.001388888888888888888888888888863490893732;
		104	const C3: f128 = 0.00002480158730158730158730158600795304914210;
		105	const C4: f128 = -0.2755731922398589065255474947078934284324e-6;
		106	const C5: f128 = 0.2087675698786809897659225313136400793948e-8;
		107	const C6: f128 = -0.1147074559772972315817149986812031204775e-10;
		108	const C7: f128 = 0.4779477332386808976875457937252120293400e-13;
		109	const C8: f64 = -0.1561920696721507929516718307820958119868e-15;
		110	const C9: f64 = 0.4110317413744594971475941557607804508039e-18;
		111	const C10: f64 = -0.8896592467191938803288521958313920156409e-21;
		112	const C11: f64 = 0.1601061435794535138244346256065192782581e-23;
		113
		114	const z = x * x;
		115	const r = z * (C1 + z * (C2 + z * (C3 + z * (C4 + z * (C5 + z * (C6 +
		116	z * (C7 + z * (C8 + z * (C9 + z * (C10 + z * C11))))))))));
		117	const hz = 0.5 * z;
		118	const w = 1.0 - hz;
		119
126	return w + (((1.0 - w) - hz) + (z * r - x * y));	120	return w + (((1.0 - w) - hz) + (z * r - x * y));
127	}	121	}
128		122
...	@@ -172,58 +166,6 @@ pub fn sin(x: f64, y: f64, iy: i32) f64 {	...	@@ -172,58 +166,6 @@ pub fn sin(x: f64, y: f64, iy: i32) f64 {
172	}	166	}
173	}	167	}
174		168
175	pub fn sinl(comptime T: type, x: T, y: T, iy: i32) T {
176	const impl = switch (T) {
177	f80 => struct {
178	const S1: T = -0.166666666666666666671;
179
180	const S2: f64 = 0.0083333333333333332;
181	const S3: f64 = -0.00019841269841269427;
182	const S4: f64 = 0.0000027557319223597490;
183	const S5: f64 = -0.000000025052108218074604;
184	const S6: f64 = 1.6059006598854211e-10;
185	const S7: f64 = -7.6429779983024564e-13;
186	const S8: f64 = 2.6174587166648325e-15;
187
188	inline fn poly(z: T) T {
189	return S2 + z * (S3 + z * (S4 + z * (S5 +
190	z * (S6 + z * (S7 + z * S8)))));
191	}
192	},
193	f128 => struct {
194	const S1: T = -0.16666666666666666666666666666666666606732416116558;
195	const S2: T = 0.0083333333333333333333333333333331135404851288270047;
196	const S3: T = -0.00019841269841269841269841269839935785325638310428717;
197	const S4: T = 0.27557319223985890652557316053039946268333231205686e-5;
198	const S5: T = -0.25052108385441718775048214826384312253862930064745e-7;
199	const S6: T = 0.16059043836821614596571832194524392581082444805729e-9;
200	const S7: T = -0.76471637318198151807063387954939213287488216303768e-12;
201	const S8: T = 0.28114572543451292625024967174638477283187397621303e-14;
202
203	const S9: f64 = -0.82206352458348947812512122163446202498005154296863e-17;
204	const S10: f64 = 0.19572940011906109418080609928334380560135358385256e-19;
205	const S11: f64 = -0.38680813379701966970673724299207480965452616911420e-22;
206	const S12: f64 = 0.64038150078671872796678569586315881020659912139412e-25;
207
208	inline fn poly(z: T) T {
209	return S2 + z * (S3 + z * (S4 + z * (S5 + z * (S6 + z * (S7 + z * (S8 +
210	z * (S9 + z * (S10 + z * (S11 + z * S12)))))))));
211	}
212	},
213	else => @compileError("sinl supports only f80 and f128, got: " ++ @typeName(T)),
214	};
215
216	const z = x * x;
217	const v = z * x;
218	const r = impl.poly(z);
219
220	if (iy == 0) {
221	return x + v * (impl.S1 + z * r);
222	}
223
224	return x - ((z * (0.5 * y - v * r) - y) - v * impl.S1);
225	}
226
227	pub fn sindf(x: f64) f32 {	169	pub fn sindf(x: f64) f32 {
228	// \|sin(x)/x - s(x)\| < 2**-37.5 (~[-4.89e-12, 4.824e-12]).	170	// \|sin(x)/x - s(x)\| < 2**-37.5 (~[-4.89e-12, 4.824e-12]).
229	const S1 = -0x15555554cbac77.0p-55; // -0.166666666416265235595	171	const S1 = -0x15555554cbac77.0p-55; // -0.166666666416265235595
...	@@ -239,6 +181,52 @@ pub fn sindf(x: f64) f32 {	...	@@ -239,6 +181,52 @@ pub fn sindf(x: f64) f32 {
239	return @floatCast((x + s * (S1 + z * S2)) + s * w * r);	181	return @floatCast((x + s * (S1 + z * S2)) + s * w * r);
240	}	182	}
241		183
		184	pub fn sinx(x: f80, y: f80, iy: i32) f80 {
		185	const S1: f80 = -0.166666666666666666671;
		186	const S2: f64 = 0.0083333333333333332;
		187	const S3: f64 = -0.00019841269841269427;
		188	const S4: f64 = 0.0000027557319223597490;
		189	const S5: f64 = -0.000000025052108218074604;
		190	const S6: f64 = 1.6059006598854211e-10;
		191	const S7: f64 = -7.6429779983024564e-13;
		192	const S8: f64 = 2.6174587166648325e-15;
		193
		194	const z = x * x;
		195	const v = z * x;
		196	const r = S2 + z * (S3 + z * (S4 + z * (S5 +
		197	z * (S6 + z * (S7 + z * S8)))));
		198
		199	if (iy == 0)
		200	return x + v * (S1 + z * r);
		201
		202	return x - ((z * (0.5 * y - v * r) - y) - v * S1);
		203	}
		204
		205	pub fn sinq(x: f128, y: f128, iy: i32) f128 {
		206	const S1: f128 = -0.16666666666666666666666666666666666606732416116558;
		207	const S2: f128 = 0.0083333333333333333333333333333331135404851288270047;
		208	const S3: f128 = -0.00019841269841269841269841269839935785325638310428717;
		209	const S4: f128 = 0.27557319223985890652557316053039946268333231205686e-5;
		210	const S5: f128 = -0.25052108385441718775048214826384312253862930064745e-7;
		211	const S6: f128 = 0.16059043836821614596571832194524392581082444805729e-9;
		212	const S7: f128 = -0.76471637318198151807063387954939213287488216303768e-12;
		213	const S8: f128 = 0.28114572543451292625024967174638477283187397621303e-14;
		214	const S9: f64 = -0.82206352458348947812512122163446202498005154296863e-17;
		215	const S10: f64 = 0.19572940011906109418080609928334380560135358385256e-19;
		216	const S11: f64 = -0.38680813379701966970673724299207480965452616911420e-22;
		217	const S12: f64 = 0.64038150078671872796678569586315881020659912139412e-25;
		218
		219	const z = x * x;
		220	const v = z * x;
		221	const r = S2 + z * (S3 + z * (S4 + z * (S5 + z * (S6 + z * (S7 + z * (S8 +
		222	z * (S9 + z * (S10 + z * (S11 + z * S12)))))))));
		223
		224	if (iy == 0)
		225	return x + v * (S1 + z * r);
		226
		227	return x - ((z * (0.5 * y - v * r) - y) - v * S1);
		228	}
		229
242	/// kernel tan function on ~[-pi/4, pi/4] (except on -0), pi/4 ~ 0.7854	230	/// kernel tan function on ~[-pi/4, pi/4] (except on -0), pi/4 ~ 0.7854
243	/// Input x is assumed to be bounded by ~pi/4 in magnitude.	231	/// Input x is assumed to be bounded by ~pi/4 in magnitude.
244	/// Input y is the tail of x.	232	/// Input y is the tail of x.
...	@@ -377,88 +365,112 @@ pub fn tandf(x: f64, odd: bool) f32 {	...	@@ -377,88 +365,112 @@ pub fn tandf(x: f64, odd: bool) f32 {
377	return @floatCast(if (odd) -1.0 / r0 else r0);	365	return @floatCast(if (odd) -1.0 / r0 else r0);
378	}	366	}
379		367
380	pub fn tanl(comptime T: type, x_: T, y_: T, odd: i32) T {	368	pub fn tanx(x_: f80, y_: f80, odd: i32) f80 {
		369	const pio4: f80 = 0.785398163397448309628;
		370	const pio4lo: f80 = -1.25413940316708300586e-20;
		371
		372	const T3: f80 = 0.333333333333333333180;
		373	const T5: f80 = 0.133333333333333372290;
		374	const T7: f80 = 0.0539682539682504975744;
		375	const T9: f64 = 0.021869488536312216;
		376	const T11: f64 = 0.0088632355256619590;
		377	const T13: f64 = 0.0035921281113786528;
		378	const T15: f64 = 0.0014558334756312418;
		379	const T17: f64 = 0.00059003538700862256;
		380	const T19: f64 = 0.00023907843576635544;
		381	const T21: f64 = 0.000097154625656538905;
		382	const T23: f64 = 0.000038440165747303162;
		383	const T25: f64 = 0.000018082171885432524;
		384	const T27: f64 = 0.0000024196006108814377;
		385	const T29: f64 = 0.0000078293456938132840;
		386	const T31: f64 = -0.0000032609076735050182;
		387	const T33: f64 = 0.0000023261313142559411;
		388
		389	var x = x_;
		390	var y = y_;
		391	const big = @abs(x) >= 0.67434;
		392	var sign: i8 = 0;
		393
		394	if (big) {
		395	if (x < 0) {
		396	sign = -1;
		397	x = -x;
		398	y = -y;
		399	}
		400	x = (pio4 - x) + (pio4lo - y);
		401	y = 0.0;
		402	}
		403
		404	var z = x * x;
		405	var w = z * z;
		406
		407	var r = T5 + w * (T9 + w * (T13 + w * (T17 + w * (T21 +
		408	w * (T25 + w * (T29 + w * T33))))));
		409
		410	var v = z * (T7 + w * (T11 + w * (T15 + w * (T19 + w * (T23 +
		411	w * (T27 + w * T31))))));
		412
		413	var s = z * x;
		414	r = y + z * (s * (r + v) + y) + T3 * s;
		415	w = x + r;
		416
		417	if (big) {
		418	s = @as(f80, @floatFromInt(1 - 2 * odd));
		419	v = s - 2.0 * (x + (r - w * w / (w + s)));
		420	return if (sign == -1) -v else v;
		421	}
		422
		423	if (odd == 0) {
		424	return w;
		425	}
		426
		427	// if allow error up to 2 ulp, simply return
		428	// -1.0 / (x+r) here
		429	//
		430	// compute -1.0 / (x+r) accurately
		431	z = w + 0x1p32 - 0x1p32;
		432	v = r - (z - x);
		433	const a = -1.0 / w;
		434	const t = a + 0x1p32 - 0x1p32;
		435	s = 1.0 + t * z;
		436	return t + a * (s + t * v);
		437	}
		438
		439	pub fn tanq(x_: f128, y_: f128, odd: i32) f128 {
		440	const pio4: f128 = 0x1.921fb54442d18469898cc51701b8p-1;
		441	const pio4lo: f128 = 0x1.cd129024e088a67cc74020bbea60p-116;
		442
		443	const T3: f128 = 0x1.5555555555555555555555555553p-2;
		444	const T5: f128 = 0x1.1111111111111111111111111eb5p-3;
		445	const T7: f128 = 0x1.ba1ba1ba1ba1ba1ba1ba1b694cd6p-5;
		446	const T9: f128 = 0x1.664f4882c10f9f32d6bbe09d8bcdp-6;
		447	const T11: f128 = 0x1.226e355e6c23c8f5b4f5762322eep-7;
		448	const T13: f128 = 0x1.d6d3d0e157ddfb5fed8e84e27b37p-9;
		449	const T15: f128 = 0x1.7da36452b75e2b5fce9ee7c2c92ep-10;
		450	const T17: f128 = 0x1.355824803674477dfcf726649efep-11;
		451	const T19: f128 = 0x1.f57d7734d1656e0aceb716f614c2p-13;
		452	const T21: f128 = 0x1.967e18afcb180ed942dfdc518d6cp-14;
		453	const T23: f128 = 0x1.497d8eea21e95bc7e2aa79b9f2cdp-15;
		454	const T25: f128 = 0x1.0b132d39f055c81be49eff7afd50p-16;
		455	const T27: f128 = 0x1.b0f72d33eff7bfa2fbc1059d90b6p-18;
		456	const T29: f128 = 0x1.5ef2daf21d1113df38d0fbc00267p-19;
		457	const T31: f128 = 0x1.1c77d6eac0234988cdaa04c96626p-20;
		458	const T33: f128 = 0x1.cd2a5a292b180e0bdd701057dfe3p-22;
		459	const T35: f128 = 0x1.75c7357d0298c01a31d0a6f7d518p-23;
		460	const T37: f128 = 0x1.2f3190f4718a9a520f98f50081fcp-24;
		461	const T39: f64 = 0.000000028443389121318352;
		462	const T41: f64 = 0.000000011981013102001973;
		463	const T43: f64 = 0.0000000038303578044958070;
		464	const T45: f64 = 0.0000000034664378216909893;
		465	const T47: f64 = -0.0000000015090641701997785;
		466	const T49: f64 = 0.0000000029449552300483952;
		467	const T51: f64 = -0.0000000022006995706097711;
		468	const T53: f64 = 0.0000000015468200913196612;
		469	const T55: f64 = -0.00000000061311613386849674;
		470	const T57: f64 = 1.4912469681508012e-10;
		471
381	var x = x_;	472	var x = x_;
382	var y = y_;	473	var y = y_;
383	const impl = switch (T) {
384	f80 => struct {
385	const pio4: T = 0.785398163397448309628;
386	const pio4lo: T = -1.25413940316708300586e-20;
387
388	const T3: T = 0.333333333333333333180;
389	const T5: T = 0.133333333333333372290;
390	const T7: T = 0.0539682539682504975744;
391	const T9: f64 = 0.021869488536312216;
392	const T11: f64 = 0.0088632355256619590;
393	const T13: f64 = 0.0035921281113786528;
394	const T15: f64 = 0.0014558334756312418;
395	const T17: f64 = 0.00059003538700862256;
396	const T19: f64 = 0.00023907843576635544;
397	const T21: f64 = 0.000097154625656538905;
398	const T23: f64 = 0.000038440165747303162;
399	const T25: f64 = 0.000018082171885432524;
400	const T27: f64 = 0.0000024196006108814377;
401	const T29: f64 = 0.0000078293456938132840;
402	const T31: f64 = -0.0000032609076735050182;
403	const T33: f64 = 0.0000023261313142559411;
404
405	inline fn rpoly(w: T) T {
406	return T5 + w * (T9 + w * (T13 + w * (T17 + w * (T21 +
407	w * (T25 + w * (T29 + w * T33))))));
408	}
409
410	inline fn vpoly(w: T) T {
411	return T7 + w * (T11 + w * (T15 + w * (T19 + w * (T23 +
412	w * (T27 + w * T31)))));
413	}
414	},
415	f128 => struct {
416	const pio4: T = 0x1.921fb54442d18469898cc51701b8p-1;
417	const pio4lo: T = 0x1.cd129024e088a67cc74020bbea60p-116;
418
419	const T3: T = 0x1.5555555555555555555555555553p-2;
420	const T5: T = 0x1.1111111111111111111111111eb5p-3;
421	const T7: T = 0x1.ba1ba1ba1ba1ba1ba1ba1b694cd6p-5;
422	const T9: T = 0x1.664f4882c10f9f32d6bbe09d8bcdp-6;
423	const T11: T = 0x1.226e355e6c23c8f5b4f5762322eep-7;
424	const T13: T = 0x1.d6d3d0e157ddfb5fed8e84e27b37p-9;
425	const T15: T = 0x1.7da36452b75e2b5fce9ee7c2c92ep-10;
426	const T17: T = 0x1.355824803674477dfcf726649efep-11;
427	const T19: T = 0x1.f57d7734d1656e0aceb716f614c2p-13;
428	const T21: T = 0x1.967e18afcb180ed942dfdc518d6cp-14;
429	const T23: T = 0x1.497d8eea21e95bc7e2aa79b9f2cdp-15;
430	const T25: T = 0x1.0b132d39f055c81be49eff7afd50p-16;
431	const T27: T = 0x1.b0f72d33eff7bfa2fbc1059d90b6p-18;
432	const T29: T = 0x1.5ef2daf21d1113df38d0fbc00267p-19;
433	const T31: T = 0x1.1c77d6eac0234988cdaa04c96626p-20;
434	const T33: T = 0x1.cd2a5a292b180e0bdd701057dfe3p-22;
435	const T35: T = 0x1.75c7357d0298c01a31d0a6f7d518p-23;
436	const T37: T = 0x1.2f3190f4718a9a520f98f50081fcp-24;
437	const T39: f64 = 0.000000028443389121318352;
438	const T41: f64 = 0.000000011981013102001973;
439	const T43: f64 = 0.0000000038303578044958070;
440	const T45: f64 = 0.0000000034664378216909893;
441	const T47: f64 = -0.0000000015090641701997785;
442	const T49: f64 = 0.0000000029449552300483952;
443	const T51: f64 = -0.0000000022006995706097711;
444	const T53: f64 = 0.0000000015468200913196612;
445	const T55: f64 = -0.00000000061311613386849674;
446	const T57: f64 = 1.4912469681508012e-10;
447
448	inline fn rpoly(w: T) T {
449	return T5 + w * (T9 + w * (T13 + w * (T17 + w * (T21 +
450	w * (T25 + w * (T29 + w * (T33 + w * (T37 + w * (T41 +
451	w * (T45 + w * (T49 + w * (T53 + w * T57))))))))))));
452	}
453
454	inline fn vpoly(w: T) T {
455	return T7 + w * (T11 + w * (T15 + w * (T19 + w * (T23 +
456	w * (T27 + w * (T31 + w * (T35 + w * (T39 + w * (T43 +
457	w * (T47 + w * (T51 + w * T55)))))))))));
458	}
459	},
460	else => @compileError("tanl supports only f80 and f128, got: " ++ @typeName(T)),
461	};
462		474
463	const big = @abs(x) >= 0.67434;	475	const big = @abs(x) >= 0.67434;
464	var sign: i8 = 0;	476	var sign: i8 = 0;
...	@@ -469,20 +481,27 @@ pub fn tanl(comptime T: type, x_: T, y_: T, odd: i32) T {	...	@@ -469,20 +481,27 @@ pub fn tanl(comptime T: type, x_: T, y_: T, odd: i32) T {
469	x = -x;	481	x = -x;
470	y = -y;	482	y = -y;
471	}	483	}
472	x = (impl.pio4 - x) + (impl.pio4lo - y);	484	x = (pio4 - x) + (pio4lo - y);
473	y = 0.0;	485	y = 0.0;
474	}	486	}
475		487
476	var z = x * x;	488	var z = x * x;
477	var w = z * z;	489	var w = z * z;
478	var r = impl.rpoly(w);	490
479	var v = z * impl.vpoly(w);	491	var r = T5 + w * (T9 + w * (T13 + w * (T17 + w * (T21 +
		492	w * (T25 + w * (T29 + w * (T33 + w * (T37 + w * (T41 +
		493	w * (T45 + w * (T49 + w * (T53 + w * T57))))))))))));
		494
		495	var v = z * (T7 + w * (T11 + w * (T15 + w * (T19 + w * (T23 +
		496	w * (T27 + w * (T31 + w * (T35 + w * (T39 + w * (T43 +
		497	w * (T47 + w * (T51 + w * T55))))))))))));
		498
480	var s = z * x;	499	var s = z * x;
481	r = y + z * (s * (r + v) + y) + impl.T3 * s;	500	r = y + z * (s * (r + v) + y) + T3 * s;
482	w = x + r;	501	w = x + r;
483		502
484	if (big) {	503	if (big) {
485	s = @as(T, @floatFromInt(1 - 2 * odd));	504	s = @as(f128, @floatFromInt(1 - 2 * odd));
486	v = s - 2.0 * (x + (r - w * w / (w + s)));	505	v = s - 2.0 * (x + (r - w * w / (w + s)));
487	return if (sign == -1) -v else v;	506	return if (sign == -1) -v else v;
488	}	507	}