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