sin.zig (6783B) - Raw
1 //! Ported from musl, which is licensed under the MIT license: 2 //! https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT 3 //! 4 //! https://git.musl-libc.org/cgit/musl/tree/src/math/sinf.c 5 //! https://git.musl-libc.org/cgit/musl/tree/src/math/sin.c 6 7 const std = @import("std"); 8 const builtin = @import("builtin"); 9 const arch = builtin.cpu.arch; 10 const math = std.math; 11 const mem = std.mem; 12 const expect = std.testing.expect; 13 const common = @import("common.zig"); 14 15 const trig = @import("trig.zig"); 16 const rem_pio2 = @import("rem_pio2.zig").rem_pio2; 17 const rem_pio2f = @import("rem_pio2f.zig").rem_pio2f; 18 19 pub const panic = common.panic; 20 21 comptime { 22 @export(&__sinh, .{ .name = "__sinh", .linkage = common.linkage, .visibility = common.visibility }); 23 @export(&sinf, .{ .name = "sinf", .linkage = common.linkage, .visibility = common.visibility }); 24 @export(&sin, .{ .name = "sin", .linkage = common.linkage, .visibility = common.visibility }); 25 @export(&__sinx, .{ .name = "__sinx", .linkage = common.linkage, .visibility = common.visibility }); 26 if (common.want_ppc_abi) { 27 @export(&sinq, .{ .name = "sinf128", .linkage = common.linkage, .visibility = common.visibility }); 28 } 29 @export(&sinq, .{ .name = "sinq", .linkage = common.linkage, .visibility = common.visibility }); 30 @export(&sinl, .{ .name = "sinl", .linkage = common.linkage, .visibility = common.visibility }); 31 } 32 33 pub fn __sinh(x: f16) callconv(.c) f16 { 34 // TODO: more efficient implementation 35 return @floatCast(sinf(x)); 36 } 37 38 pub fn sinf(x: f32) callconv(.c) f32 { 39 // Small multiples of pi/2 rounded to double precision. 40 const s1pio2: f64 = 1.0 * math.pi / 2.0; // 0x3FF921FB, 0x54442D18 41 const s2pio2: f64 = 2.0 * math.pi / 2.0; // 0x400921FB, 0x54442D18 42 const s3pio2: f64 = 3.0 * math.pi / 2.0; // 0x4012D97C, 0x7F3321D2 43 const s4pio2: f64 = 4.0 * math.pi / 2.0; // 0x401921FB, 0x54442D18 44 45 var ix: u32 = @bitCast(x); 46 const sign = ix >> 31 != 0; 47 ix &= 0x7fffffff; 48 49 if (ix <= 0x3f490fda) { // |x| ~<= pi/4 50 if (ix < 0x39800000) { // |x| < 2**-12 51 // raise inexact if x!=0 and underflow if subnormal 52 if (common.want_float_exceptions) { 53 if (ix < 0x00800000) { 54 mem.doNotOptimizeAway(x / 0x1p120); 55 } else { 56 mem.doNotOptimizeAway(x + 0x1p120); 57 } 58 } 59 return x; 60 } 61 return trig.__sindf(x); 62 } 63 if (ix <= 0x407b53d1) { // |x| ~<= 5*pi/4 64 if (ix <= 0x4016cbe3) { // |x| ~<= 3pi/4 65 if (sign) { 66 return -trig.__cosdf(x + s1pio2); 67 } else { 68 return trig.__cosdf(x - s1pio2); 69 } 70 } 71 return trig.__sindf(if (sign) -(x + s2pio2) else -(x - s2pio2)); 72 } 73 if (ix <= 0x40e231d5) { // |x| ~<= 9*pi/4 74 if (ix <= 0x40afeddf) { // |x| ~<= 7*pi/4 75 if (sign) { 76 return trig.__cosdf(x + s3pio2); 77 } else { 78 return -trig.__cosdf(x - s3pio2); 79 } 80 } 81 return trig.__sindf(if (sign) x + s4pio2 else x - s4pio2); 82 } 83 84 // sin(Inf or NaN) is NaN 85 if (ix >= 0x7f800000) { 86 return x - x; 87 } 88 89 var y: f64 = undefined; 90 const n = rem_pio2f(x, &y); 91 return switch (n & 3) { 92 0 => trig.__sindf(y), 93 1 => trig.__cosdf(y), 94 2 => trig.__sindf(-y), 95 else => -trig.__cosdf(y), 96 }; 97 } 98 99 pub fn sin(x: f64) callconv(.c) f64 { 100 var ix = @as(u64, @bitCast(x)) >> 32; 101 ix &= 0x7fffffff; 102 103 // |x| ~< pi/4 104 if (ix <= 0x3fe921fb) { 105 if (ix < 0x3e500000) { // |x| < 2**-26 106 // raise inexact if x != 0 and underflow if subnormal 107 if (common.want_float_exceptions) { 108 if (ix < 0x00100000) { 109 mem.doNotOptimizeAway(x / 0x1p120); 110 } else { 111 mem.doNotOptimizeAway(x + 0x1p120); 112 } 113 } 114 return x; 115 } 116 return trig.__sin(x, 0.0, 0); 117 } 118 119 // sin(Inf or NaN) is NaN 120 if (ix >= 0x7ff00000) { 121 return x - x; 122 } 123 124 var y: [2]f64 = undefined; 125 const n = rem_pio2(x, &y); 126 return switch (n & 3) { 127 0 => trig.__sin(y[0], y[1], 1), 128 1 => trig.__cos(y[0], y[1]), 129 2 => -trig.__sin(y[0], y[1], 1), 130 else => -trig.__cos(y[0], y[1]), 131 }; 132 } 133 134 pub fn __sinx(x: f80) callconv(.c) f80 { 135 // TODO: more efficient implementation 136 return @floatCast(sinq(x)); 137 } 138 139 pub fn sinq(x: f128) callconv(.c) f128 { 140 // TODO: more correct implementation 141 return sin(@floatCast(x)); 142 } 143 144 pub fn sinl(x: c_longdouble) callconv(.c) c_longdouble { 145 switch (@typeInfo(c_longdouble).float.bits) { 146 16 => return __sinh(x), 147 32 => return sinf(x), 148 64 => return sin(x), 149 80 => return __sinx(x), 150 128 => return sinq(x), 151 else => @compileError("unreachable"), 152 } 153 } 154 155 test "sin32" { 156 const epsilon = 0.00001; 157 158 try expect(math.approxEqAbs(f32, sinf(0.0), 0.0, epsilon)); 159 try expect(math.approxEqAbs(f32, sinf(0.2), 0.198669, epsilon)); 160 try expect(math.approxEqAbs(f32, sinf(0.8923), 0.778517, epsilon)); 161 try expect(math.approxEqAbs(f32, sinf(1.5), 0.997495, epsilon)); 162 try expect(math.approxEqAbs(f32, sinf(-1.5), -0.997495, epsilon)); 163 try expect(math.approxEqAbs(f32, sinf(37.45), -0.246544, epsilon)); 164 try expect(math.approxEqAbs(f32, sinf(89.123), 0.916166, epsilon)); 165 } 166 167 test "sin64" { 168 const epsilon = 0.000001; 169 170 try expect(math.approxEqAbs(f64, sin(0.0), 0.0, epsilon)); 171 try expect(math.approxEqAbs(f64, sin(0.2), 0.198669, epsilon)); 172 try expect(math.approxEqAbs(f64, sin(0.8923), 0.778517, epsilon)); 173 try expect(math.approxEqAbs(f64, sin(1.5), 0.997495, epsilon)); 174 try expect(math.approxEqAbs(f64, sin(-1.5), -0.997495, epsilon)); 175 try expect(math.approxEqAbs(f64, sin(37.45), -0.246543, epsilon)); 176 try expect(math.approxEqAbs(f64, sin(89.123), 0.916166, epsilon)); 177 } 178 179 test "sin32.special" { 180 try expect(sinf(0.0) == 0.0); 181 try expect(sinf(-0.0) == -0.0); 182 try expect(math.isNan(sinf(math.inf(f32)))); 183 try expect(math.isNan(sinf(-math.inf(f32)))); 184 try expect(math.isNan(sinf(math.nan(f32)))); 185 } 186 187 test "sin64.special" { 188 try expect(sin(0.0) == 0.0); 189 try expect(sin(-0.0) == -0.0); 190 try expect(math.isNan(sin(math.inf(f64)))); 191 try expect(math.isNan(sin(-math.inf(f64)))); 192 try expect(math.isNan(sin(math.nan(f64)))); 193 } 194 195 test "sin32 #9901" { 196 const float: f32 = @bitCast(@as(u32, 0b11100011111111110000000000000000)); 197 _ = sinf(float); 198 } 199 200 test "sin64 #9901" { 201 const float: f64 = @bitCast(@as(u64, 0b1111111101000001000000001111110111111111100000000000000000000001)); 202 _ = sin(float); 203 }