cos.zig (5691B) - Raw
1 const std = @import("std"); 2 const math = std.math; 3 const mem = std.mem; 4 const expect = std.testing.expect; 5 const common = @import("common.zig"); 6 7 pub const panic = common.panic; 8 9 const trig = @import("trig.zig"); 10 const rem_pio2 = @import("rem_pio2.zig").rem_pio2; 11 const rem_pio2f = @import("rem_pio2f.zig").rem_pio2f; 12 13 comptime { 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 25 pub fn __cosh(a: f16) callconv(.c) f16 { 26 // TODO: more efficient implementation 27 return @floatCast(cosf(a)); 28 } 29 30 pub 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 87 pub 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 116 pub fn __cosx(a: f80) callconv(.c) f80 { 117 // TODO: more efficient implementation 118 return @floatCast(cosq(a)); 119 } 120 121 pub fn cosq(a: f128) callconv(.c) f128 { 122 // TODO: more correct implementation 123 return cos(@floatCast(a)); 124 } 125 126 pub 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 137 test "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 149 test "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 161 test "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 167 test "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 }