master
1const std = @import("std");
2const math = std.math;
3const mem = std.mem;
4const expect = std.testing.expect;
5const common = @import("common.zig");
6
7pub const panic = common.panic;
8
9const trig = @import("trig.zig");
10const rem_pio2 = @import("rem_pio2.zig").rem_pio2;
11const rem_pio2f = @import("rem_pio2f.zig").rem_pio2f;
12
13comptime {
14 @export(&__cosh, .{ .name = "__cosh", .linkage = common.linkage, .visibility = common.visibility });
15 @export(&cosf, .{ .name = "cosf", .linkage = common.linkage, .visibility = common.visibility });
16 @export(&cos, .{ .name = "cos", .linkage = common.linkage, .visibility = common.visibility });
17 @export(&__cosx, .{ .name = "__cosx", .linkage = common.linkage, .visibility = common.visibility });
18 if (common.want_ppc_abi) {
19 @export(&cosq, .{ .name = "cosf128", .linkage = common.linkage, .visibility = common.visibility });
20 }
21 @export(&cosq, .{ .name = "cosq", .linkage = common.linkage, .visibility = common.visibility });
22 @export(&cosl, .{ .name = "cosl", .linkage = common.linkage, .visibility = common.visibility });
23}
24
25pub fn __cosh(a: f16) callconv(.c) f16 {
26 // TODO: more efficient implementation
27 return @floatCast(cosf(a));
28}
29
30pub fn cosf(x: f32) callconv(.c) f32 {
31 // Small multiples of pi/2 rounded to double precision.
32 const c1pio2: f64 = 1.0 * math.pi / 2.0; // 0x3FF921FB, 0x54442D18
33 const c2pio2: f64 = 2.0 * math.pi / 2.0; // 0x400921FB, 0x54442D18
34 const c3pio2: f64 = 3.0 * math.pi / 2.0; // 0x4012D97C, 0x7F3321D2
35 const c4pio2: f64 = 4.0 * math.pi / 2.0; // 0x401921FB, 0x54442D18
36
37 var ix: u32 = @bitCast(x);
38 const sign = ix >> 31 != 0;
39 ix &= 0x7fffffff;
40
41 if (ix <= 0x3f490fda) { // |x| ~<= pi/4
42 if (ix < 0x39800000) { // |x| < 2**-12
43 // raise inexact if x != 0
44 if (common.want_float_exceptions) mem.doNotOptimizeAway(x + 0x1p120);
45 return 1.0;
46 }
47 return trig.__cosdf(x);
48 }
49 if (ix <= 0x407b53d1) { // |x| ~<= 5*pi/4
50 if (ix > 0x4016cbe3) { // |x| ~> 3*pi/4
51 return -trig.__cosdf(if (sign) x + c2pio2 else x - c2pio2);
52 } else {
53 if (sign) {
54 return trig.__sindf(x + c1pio2);
55 } else {
56 return trig.__sindf(c1pio2 - x);
57 }
58 }
59 }
60 if (ix <= 0x40e231d5) { // |x| ~<= 9*pi/4
61 if (ix > 0x40afeddf) { // |x| ~> 7*pi/4
62 return trig.__cosdf(if (sign) x + c4pio2 else x - c4pio2);
63 } else {
64 if (sign) {
65 return trig.__sindf(-x - c3pio2);
66 } else {
67 return trig.__sindf(x - c3pio2);
68 }
69 }
70 }
71
72 // cos(Inf or NaN) is NaN
73 if (ix >= 0x7f800000) {
74 return x - x;
75 }
76
77 var y: f64 = undefined;
78 const n = rem_pio2f(x, &y);
79 return switch (n & 3) {
80 0 => trig.__cosdf(y),
81 1 => trig.__sindf(-y),
82 2 => -trig.__cosdf(y),
83 else => trig.__sindf(y),
84 };
85}
86
87pub fn cos(x: f64) callconv(.c) f64 {
88 var ix = @as(u64, @bitCast(x)) >> 32;
89 ix &= 0x7fffffff;
90
91 // |x| ~< pi/4
92 if (ix <= 0x3fe921fb) {
93 if (ix < 0x3e46a09e) { // |x| < 2**-27 * sqrt(2)
94 // raise inexact if x!=0
95 if (common.want_float_exceptions) mem.doNotOptimizeAway(x + 0x1p120);
96 return 1.0;
97 }
98 return trig.__cos(x, 0);
99 }
100
101 // cos(Inf or NaN) is NaN
102 if (ix >= 0x7ff00000) {
103 return x - x;
104 }
105
106 var y: [2]f64 = undefined;
107 const n = rem_pio2(x, &y);
108 return switch (n & 3) {
109 0 => trig.__cos(y[0], y[1]),
110 1 => -trig.__sin(y[0], y[1], 1),
111 2 => -trig.__cos(y[0], y[1]),
112 else => trig.__sin(y[0], y[1], 1),
113 };
114}
115
116pub fn __cosx(a: f80) callconv(.c) f80 {
117 // TODO: more efficient implementation
118 return @floatCast(cosq(a));
119}
120
121pub fn cosq(a: f128) callconv(.c) f128 {
122 // TODO: more correct implementation
123 return cos(@floatCast(a));
124}
125
126pub fn cosl(x: c_longdouble) callconv(.c) c_longdouble {
127 switch (@typeInfo(c_longdouble).float.bits) {
128 16 => return __cosh(x),
129 32 => return cosf(x),
130 64 => return cos(x),
131 80 => return __cosx(x),
132 128 => return cosq(x),
133 else => @compileError("unreachable"),
134 }
135}
136
137test "cos32" {
138 const epsilon = 0.00001;
139
140 try expect(math.approxEqAbs(f32, cosf(0.0), 1.0, epsilon));
141 try expect(math.approxEqAbs(f32, cosf(0.2), 0.980067, epsilon));
142 try expect(math.approxEqAbs(f32, cosf(0.8923), 0.627623, epsilon));
143 try expect(math.approxEqAbs(f32, cosf(1.5), 0.070737, epsilon));
144 try expect(math.approxEqAbs(f32, cosf(-1.5), 0.070737, epsilon));
145 try expect(math.approxEqAbs(f32, cosf(37.45), 0.969132, epsilon));
146 try expect(math.approxEqAbs(f32, cosf(89.123), 0.400798, epsilon));
147}
148
149test "cos64" {
150 const epsilon = 0.000001;
151
152 try expect(math.approxEqAbs(f64, cos(0.0), 1.0, epsilon));
153 try expect(math.approxEqAbs(f64, cos(0.2), 0.980067, epsilon));
154 try expect(math.approxEqAbs(f64, cos(0.8923), 0.627623, epsilon));
155 try expect(math.approxEqAbs(f64, cos(1.5), 0.070737, epsilon));
156 try expect(math.approxEqAbs(f64, cos(-1.5), 0.070737, epsilon));
157 try expect(math.approxEqAbs(f64, cos(37.45), 0.969132, epsilon));
158 try expect(math.approxEqAbs(f64, cos(89.123), 0.40080, epsilon));
159}
160
161test "cos32.special" {
162 try expect(math.isNan(cosf(math.inf(f32))));
163 try expect(math.isNan(cosf(-math.inf(f32))));
164 try expect(math.isNan(cosf(math.nan(f32))));
165}
166
167test "cos64.special" {
168 try expect(math.isNan(cos(math.inf(f64))));
169 try expect(math.isNan(cos(-math.inf(f64))));
170 try expect(math.isNan(cos(math.nan(f64))));
171}