authorgravatar for andrew@ziglang.orgAndrew Kelley <andrew@ziglang.org> 2026-02-13 12:26:11-08:00
committergravatar for andrew@ziglang.orgAndrew Kelley <andrew@ziglang.org> 2026-02-13 12:26:11-08:00
loge262a32ad1d039ee575009c8ea770743b72e8160
treeb8b0268088b056e295422c115fdd5f36df253653
parent07b08b96389a1dd60bc150a5ed72586ee57b2303
parent12e6ac2e8fb034f823390e5d7f1680411c903e7c

Merge branch 'std.math.atan: Add @Vector(?,f32) and @Vector(?,f64) support'

Add SIMD support for atan (f32 and f64), based on the [ARM impl](https://github.com/ARM-software/optimized-routines/blob/master/math/aarch64/advsimd/atanf.c). To reduce branching, more polynomial approximation is used. Reviewed-on: https://codeberg.org/ziglang/zig/pulls/31195

1 files changed, 205 insertions(+), 8 deletions(-)

lib/std/math/atan.zig+205-8
...@@ -4,6 +4,12 @@...@@ -4,6 +4,12 @@
4// https://git.musl-libc.org/cgit/musl/tree/src/math/atanf.c4// https://git.musl-libc.org/cgit/musl/tree/src/math/atanf.c
5// https://git.musl-libc.org/cgit/musl/tree/src/math/atan.c5// https://git.musl-libc.org/cgit/musl/tree/src/math/atan.c
6// https://git.musl-libc.org/cgit/musl/tree/src/math/atanl.c6// https://git.musl-libc.org/cgit/musl/tree/src/math/atanl.c
7//
8// Ported from ARM-software, which is licensed under the MIT license:
9// https://github.com/ARM-software/optimized-routines/blob/master/LICENSE
10//
11// https://github.com/ARM-software/optimized-routines/blob/master/math/aarch64/advsimd/atanf.c
12// https://github.com/ARM-software/optimized-routines/blob/master/math/aarch64/advsimd/atan.c
713
8const std = @import("../std.zig");14const std = @import("../std.zig");
9const math = std.math;15const math = std.math;
...@@ -17,14 +23,22 @@ const testing = std.testing;...@@ -17,14 +23,22 @@ const testing = std.testing;
17/// - atan(+-inf) = +-pi/223/// - atan(+-inf) = +-pi/2
18pub fn atan(x: anytype) @TypeOf(x) {24pub fn atan(x: anytype) @TypeOf(x) {
19 const T = @TypeOf(x);25 const T = @TypeOf(x);
20 return switch (T) {26 switch (@typeInfo(T)) {
21 f16 => atanBinary16(x),27 .float => |info| switch (info.bits) {
22 f32 => atanBinary32(x),28 16 => return atanBinary16(x),
23 f64 => atanBinary64(x),29 32 => return atanBinary32(x),
24 f80 => atanExtended80(x),30 64 => return atanBinary64(x),
25 f128 => atanBinary128(x),31 80 => return atanExtended80(x),
26 else => @compileError("atan not implemented for " ++ @typeName(T)),32 128 => return atanBinary128(x),
27 };33 else => comptime unreachable,
34 },
35 .vector => |info| switch (info.child) {
36 f32 => return atanBinary32Vec(info.len, x),
37 f64 => return atanBinary64Vec(info.len, x),
38 else => @compileError("unimplemented"),
39 },
40 else => comptime unreachable,
41 }
28}42}
2943
30fn atanBinary16(x: f16) f16 {44fn atanBinary16(x: f16) f16 {
...@@ -579,3 +593,186 @@ test "atanBinary128" {...@@ -579,3 +593,186 @@ test "atanBinary128" {
579 try testing.expectApproxEqAbs(atanBinary128(0x1.299d54ac7d6afc5154643b601519p1), 0x1.2a24e22d861debfd6f974500567fp0, math.floatEpsAt(f128, 0x1.2a24e22d861debfd6f974500567fp0));593 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));594 try testing.expectApproxEqAbs(atanBinary128(-0x1.0264fb9f3d50e4f0f966f0686064p1), -0x1.1c617825f97512b7f38656ab12cdp0, math.floatEpsAt(f128, -0x1.1c617825f97512b7f38656ab12cdp0));
581}595}
596
597fn atanBinary32Vec(comptime vec_len: comptime_int, x: @Vector(vec_len, f32)) @TypeOf(x) {
598 const sign_mask: @Vector(vec_len, u32) = @splat(0x80000000);
599 const neg_one: @Vector(vec_len, f32) = @splat(-1.0);
600 const pi_over_2: @Vector(vec_len, u32) = @splat(0x3fc90fdb);
601 const zero: @Vector(vec_len, u32) = @splat(0);
602 const c0: @Vector(vec_len, f32) = @splat(-0x1.5554dcp-2);
603 const c1: @Vector(vec_len, f32) = @splat(0x1.9978ecp-3);
604 const c2: @Vector(vec_len, f32) = @splat(-0x1.230a94p-3);
605 const c3: @Vector(vec_len, f32) = @splat(0x1.b4debp-4);
606 const c4: @Vector(vec_len, f32) = @splat(-0x1.3550dap-4);
607 const c5: @Vector(vec_len, f32) = @splat(0x1.61eebp-5);
608 const c6: @Vector(vec_len, f32) = @splat(-0x1.0c17d4p-6);
609 const c7: @Vector(vec_len, f32) = @splat(0x1.7ea694p-9);
610
611 const ix: @Vector(vec_len, u32) = @bitCast(x);
612 const sign = ix & sign_mask;
613 const pred = @abs(x) > @abs(neg_one);
614 const z = @select(f32, pred, neg_one / x, x);
615 const shift: @Vector(vec_len, f32) = @bitCast(@select(u32, pred, pi_over_2 ^ sign, zero));
616 const z2 = z * z;
617 const z3 = z * z2;
618 const z4 = z2 * z2;
619 const z8 = z4 * z4;
620 const p0_1 = @mulAdd(@Vector(vec_len, f32), z2, c1, c0);
621 const p2_3 = @mulAdd(@Vector(vec_len, f32), z2, c3, c2);
622 const p4_5 = @mulAdd(@Vector(vec_len, f32), z2, c5, c4);
623 const p6_7 = @mulAdd(@Vector(vec_len, f32), z2, c7, c6);
624 const p0_3 = @mulAdd(@Vector(vec_len, f32), z4, p2_3, p0_1);
625 const p4_7 = @mulAdd(@Vector(vec_len, f32), z4, p6_7, p4_5);
626 const p0_7 = @mulAdd(@Vector(vec_len, f32), z8, p4_7, p0_3);
627 return @mulAdd(@Vector(vec_len, f32), z3, p0_7, shift + z);
628}
629
630fn atanBinary64Vec(comptime vec_len: comptime_int, x: @Vector(vec_len, f64)) @TypeOf(x) {
631 const sign_mask: @Vector(vec_len, u64) = @splat(0x8000000000000000);
632 const neg_one: @Vector(vec_len, f64) = @splat(-1.0);
633 const pi_over_2: @Vector(vec_len, u64) = @splat(0x3ff921fb54442d18);
634 const zero: @Vector(vec_len, u64) = @splat(0);
635 const c0: @Vector(vec_len, f64) = @splat(-0x1.555555555552ap-2);
636 const c1: @Vector(vec_len, f64) = @splat(0x1.9999999995aebp-3);
637 const c2: @Vector(vec_len, f64) = @splat(-0x1.24924923923f6p-3);
638 const c3: @Vector(vec_len, f64) = @splat(0x1.c71c7184288a2p-4);
639 const c4: @Vector(vec_len, f64) = @splat(-0x1.745d11fb3d32bp-4);
640 const c5: @Vector(vec_len, f64) = @splat(0x1.3b136a18051b9p-4);
641 const c6: @Vector(vec_len, f64) = @splat(-0x1.110e6d985f496p-4);
642 const c7: @Vector(vec_len, f64) = @splat(0x1.e1bcf7f08801dp-5);
643 const c8: @Vector(vec_len, f64) = @splat(-0x1.ae644e28058c3p-5);
644 const c9: @Vector(vec_len, f64) = @splat(0x1.82eeb1fed85c6p-5);
645 const c10: @Vector(vec_len, f64) = @splat(-0x1.59d7f901566cbp-5);
646 const c11: @Vector(vec_len, f64) = @splat(0x1.2c982855ab069p-5);
647 const c12: @Vector(vec_len, f64) = @splat(-0x1.eb49592998177p-6);
648 const c13: @Vector(vec_len, f64) = @splat(0x1.69d8b396e3d38p-6);
649 const c14: @Vector(vec_len, f64) = @splat(-0x1.ca980345c4204p-7);
650 const c15: @Vector(vec_len, f64) = @splat(0x1.dc050eafde0b3p-8);
651 const c16: @Vector(vec_len, f64) = @splat(-0x1.7ea70755b8eccp-9);
652 const c17: @Vector(vec_len, f64) = @splat(0x1.ba3da3de903e8p-11);
653 const c18: @Vector(vec_len, f64) = @splat(-0x1.44a4b059b6f67p-13);
654 const c19: @Vector(vec_len, f64) = @splat(0x1.c4a45029e5a91p-17);
655
656 const ix: @Vector(vec_len, u64) = @bitCast(x);
657 const sign = ix & sign_mask;
658 const pred = @abs(x) > @abs(neg_one);
659 const shift: @Vector(vec_len, f64) = @bitCast(@select(u64, pred, pi_over_2 ^ sign, zero));
660 const z = @select(f64, pred, neg_one / x, x);
661 const z2 = z * z;
662 const z3 = z * z2;
663 const z4 = z2 * z2;
664 const z8 = z4 * z4;
665 const z16 = z8 * z8;
666 const p0_1 = @mulAdd(@Vector(vec_len, f64), z2, c1, c0);
667 const p2_3 = @mulAdd(@Vector(vec_len, f64), z2, c3, c2);
668 const p0_3 = @mulAdd(@Vector(vec_len, f64), z4, p2_3, p0_1);
669 const p4_5 = @mulAdd(@Vector(vec_len, f64), z2, c5, c4);
670 const p6_7 = @mulAdd(@Vector(vec_len, f64), z2, c7, c6);
671 const p4_7 = @mulAdd(@Vector(vec_len, f64), z4, p6_7, p4_5);
672 const p0_7 = @mulAdd(@Vector(vec_len, f64), z8, p4_7, p0_3);
673 const p8_9 = @mulAdd(@Vector(vec_len, f64), z2, c9, c8);
674 const p10_11 = @mulAdd(@Vector(vec_len, f64), z2, c11, c10);
675 const p8_11 = @mulAdd(@Vector(vec_len, f64), z4, p10_11, p8_9);
676 const p12_13 = @mulAdd(@Vector(vec_len, f64), z2, c13, c12);
677 const p14_15 = @mulAdd(@Vector(vec_len, f64), z2, c15, c14);
678 const p12_15 = @mulAdd(@Vector(vec_len, f64), z4, p14_15, p12_13);
679 const p16_17 = @mulAdd(@Vector(vec_len, f64), z2, c17, c16);
680 const p18_19 = @mulAdd(@Vector(vec_len, f64), z2, c19, c18);
681 const p16_19 = @mulAdd(@Vector(vec_len, f64), z4, p18_19, p16_17);
682 const p8_15 = @mulAdd(@Vector(vec_len, f64), z8, p12_15, p8_11);
683 const p8_19 = @mulAdd(@Vector(vec_len, f64), z16, p16_19, p8_15);
684 const p0_19 = @mulAdd(@Vector(vec_len, f64), p8_19, z16, p0_7);
685 return @mulAdd(@Vector(vec_len, f64), z3, p0_19, shift + z);
686}
687
688test "atanBinary32Vec.special" {
689 const input: @Vector(7, f32) = .{
690 0x0p+0,
691 -0x0p+0,
692 0x1p+0,
693 -0x1p+0,
694 math.inf(f32),
695 -math.inf(f32),
696 math.nan(f32),
697 };
698 const output = atanBinary32Vec(7, input);
699 try testing.expectEqual(output[0], 0x0p+0);
700 try testing.expectEqual(output[1], -0x0p+0);
701 try testing.expectApproxEqAbs(output[2], 0x1.921fb6p-1, math.floatEpsAt(f32, 0x1.921fb6p-1));
702 try testing.expectApproxEqAbs(output[3], -0x1.921fb6p-1, math.floatEpsAt(f32, -0x1.921fb6p-1));
703 try testing.expectApproxEqAbs(output[4], 0x1.921fb6p+0, math.floatEpsAt(f32, 0x1.921fb6p+0));
704 try testing.expectApproxEqAbs(output[5], -0x1.921fb6p+0, math.floatEpsAt(f32, -0x1.921fb6p+0));
705 try testing.expect(math.isNan(output[6]));
706}
707
708test "atanBinary32Vec" {
709 const input: @Vector(10, f32) = .{
710 -0x1.8629dp-2,
711 -0x1.59d42ep1,
712 -0x1.d2dbe2p0,
713 -0x1.5f314ep-1,
714 0x1.5869bp1,
715 -0x1.b13a06p-2,
716 0x1.3cb0f2p1,
717 -0x1.0ed746p-2,
718 0x1.299d54p1,
719 -0x1.0264fcp1,
720 };
721 const output = atanBinary32Vec(10, input);
722 try testing.expectApproxEqAbs(output[0], -0x1.74c62p-2, math.floatEpsAt(f32, -0x1.74c62p-2));
723 try testing.expectApproxEqAbs(output[1], -0x1.375fd8p0, math.floatEpsAt(f32, -0x1.375fd8p0));
724 try testing.expectApproxEqAbs(output[2], -0x1.11b8aep0, math.floatEpsAt(f32, -0x1.11b8aep0));
725 try testing.expectApproxEqAbs(output[3], -0x1.33d28cp-1, math.floatEpsAt(f32, -0x1.33d28cp-1));
726 try testing.expectApproxEqAbs(output[4], 0x1.37082ep0, math.floatEpsAt(f32, 0x1.37082ep0));
727 try testing.expectApproxEqAbs(output[5], -0x1.99d7cap-2, math.floatEpsAt(f32, -0x1.99d7cap-2));
728 try testing.expectApproxEqAbs(output[6], 0x1.2fcb12p0, math.floatEpsAt(f32, 0x1.2fcb12p0));
729 try testing.expectApproxEqAbs(output[7], -0x1.08c71ap-2, math.floatEpsAt(f32, -0x1.08c71ap-2));
730 try testing.expectApproxEqAbs(output[8], 0x1.2a24e2p0, math.floatEpsAt(f32, 0x1.2a24e2p0));
731 try testing.expectApproxEqAbs(output[9], -0x1.1c6178p0, math.floatEpsAt(f32, -0x1.1c6178p0));
732}
733
734test "atanBinary64Vec.special" {
735 const input: @Vector(7, f64) = .{
736 0x0p+0,
737 -0x0p+0,
738 0x1p+0,
739 -0x1p+0,
740 math.inf(f64),
741 -math.inf(f64),
742 math.nan(f64),
743 };
744 const output = atanBinary64Vec(7, input);
745 try testing.expectEqual(output[0], 0x0p+0);
746 try testing.expectEqual(output[1], -0x0p+0);
747 try testing.expectApproxEqAbs(output[2], 0x1.921fb54442d18p-1, math.floatEpsAt(f64, 0x1.921fb54442d18p-1));
748 try testing.expectApproxEqAbs(output[3], -0x1.921fb54442d18p-1, math.floatEpsAt(f64, -0x1.921fb54442d18p-1));
749 try testing.expectApproxEqAbs(output[4], 0x1.921fb54442d18p+0, math.floatEpsAt(f64, 0x1.921fb54442d18p+0));
750 try testing.expectApproxEqAbs(output[5], -0x1.921fb54442d18p+0, math.floatEpsAt(f64, -0x1.921fb54442d18p+0));
751 try testing.expect(math.isNan(output[6]));
752}
753
754test "atanBinary64Vec" {
755 const input: @Vector(10, f64) = .{
756 -0x1.8629d0244cdccp-2,
757 -0x1.59d42d4659937p1,
758 -0x1.d2dbe23d04f06p0,
759 -0x1.5f314e72398e8p-1,
760 0x1.5869af37b7d08p1,
761 -0x1.b13a05a662618p-2,
762 0x1.3cb0f12f39d8ap1,
763 -0x1.0ed746b39cbb7p-2,
764 0x1.299d54ac7d6bp1,
765 -0x1.0264fb9f3d50ep1,
766 };
767 const output = atanBinary64Vec(10, input);
768 try testing.expectApproxEqAbs(output[0], -0x1.74c61f4377016p-2, math.floatEpsAt(f64, -0x1.74c61f4377016p-2));
769 try testing.expectApproxEqAbs(output[1], -0x1.375fd7987cc2p0, math.floatEpsAt(f64, -0x1.375fd7987cc2p0));
770 try testing.expectApproxEqAbs(output[2], -0x1.11b8adeba5616p0, math.floatEpsAt(f64, -0x1.11b8adeba5616p0));
771 try testing.expectApproxEqAbs(output[3], -0x1.33d28ca762539p-1, math.floatEpsAt(f64, -0x1.33d28ca762539p-1));
772 try testing.expectApproxEqAbs(output[4], 0x1.37082ce2dd03p0, math.floatEpsAt(f64, 0x1.37082ce2dd03p0));
773 try testing.expectApproxEqAbs(output[5], -0x1.99d7cac66dd44p-2, math.floatEpsAt(f64, -0x1.99d7cac66dd44p-2));
774 try testing.expectApproxEqAbs(output[6], 0x1.2fcb120468e8ep0, math.floatEpsAt(f64, 0x1.2fcb120468e8ep0));
775 try testing.expectApproxEqAbs(output[7], -0x1.08c71aa0e509p-2, math.floatEpsAt(f64, -0x1.08c71aa0e509p-2));
776 try testing.expectApproxEqAbs(output[8], 0x1.2a24e22d861dfp0, math.floatEpsAt(f64, 0x1.2a24e22d861dfp0));
777 try testing.expectApproxEqAbs(output[9], -0x1.1c617825f9751p0, math.floatEpsAt(f64, -0x1.1c617825f9751p0));
778}