sincos.zig (8779B) - Raw
1 const std = @import("std"); 2 const builtin = @import("builtin"); 3 const arch = builtin.cpu.arch; 4 const math = std.math; 5 const mem = std.mem; 6 const trig = @import("trig.zig"); 7 const rem_pio2 = @import("rem_pio2.zig").rem_pio2; 8 const rem_pio2f = @import("rem_pio2f.zig").rem_pio2f; 9 const common = @import("common.zig"); 10 11 pub const panic = common.panic; 12 13 comptime { 14 @export(&__sincosh, .{ .name = "__sincosh", .linkage = common.linkage, .visibility = common.visibility }); 15 @export(&sincosf, .{ .name = "sincosf", .linkage = common.linkage, .visibility = common.visibility }); 16 @export(&sincos, .{ .name = "sincos", .linkage = common.linkage, .visibility = common.visibility }); 17 @export(&__sincosx, .{ .name = "__sincosx", .linkage = common.linkage, .visibility = common.visibility }); 18 if (common.want_ppc_abi) { 19 @export(&sincosq, .{ .name = "sincosf128", .linkage = common.linkage, .visibility = common.visibility }); 20 } 21 @export(&sincosq, .{ .name = "sincosq", .linkage = common.linkage, .visibility = common.visibility }); 22 @export(&sincosl, .{ .name = "sincosl", .linkage = common.linkage, .visibility = common.visibility }); 23 } 24 25 pub fn __sincosh(x: f16, r_sin: *f16, r_cos: *f16) callconv(.c) void { 26 // TODO: more efficient implementation 27 var big_sin: f32 = undefined; 28 var big_cos: f32 = undefined; 29 sincosf(x, &big_sin, &big_cos); 30 r_sin.* = @as(f16, @floatCast(big_sin)); 31 r_cos.* = @as(f16, @floatCast(big_cos)); 32 } 33 34 pub fn sincosf(x: f32, r_sin: *f32, r_cos: *f32) callconv(.c) void { 35 const sc1pio2: f64 = 1.0 * math.pi / 2.0; // 0x3FF921FB, 0x54442D18 36 const sc2pio2: f64 = 2.0 * math.pi / 2.0; // 0x400921FB, 0x54442D18 37 const sc3pio2: f64 = 3.0 * math.pi / 2.0; // 0x4012D97C, 0x7F3321D2 38 const sc4pio2: f64 = 4.0 * math.pi / 2.0; // 0x401921FB, 0x54442D18 39 40 const pre_ix = @as(u32, @bitCast(x)); 41 const sign = pre_ix >> 31 != 0; 42 const ix = pre_ix & 0x7fffffff; 43 44 // |x| ~<= pi/4 45 if (ix <= 0x3f490fda) { 46 // |x| < 2**-12 47 if (ix < 0x39800000) { 48 // raise inexact if x!=0 and underflow if subnormal 49 if (common.want_float_exceptions) { 50 if (ix < 0x00100000) { 51 mem.doNotOptimizeAway(x / 0x1p120); 52 } else { 53 mem.doNotOptimizeAway(x + 0x1p120); 54 } 55 } 56 r_sin.* = x; 57 r_cos.* = 1.0; 58 return; 59 } 60 r_sin.* = trig.__sindf(x); 61 r_cos.* = trig.__cosdf(x); 62 return; 63 } 64 65 // |x| ~<= 5*pi/4 66 if (ix <= 0x407b53d1) { 67 // |x| ~<= 3pi/4 68 if (ix <= 0x4016cbe3) { 69 if (sign) { 70 r_sin.* = -trig.__cosdf(x + sc1pio2); 71 r_cos.* = trig.__sindf(x + sc1pio2); 72 } else { 73 r_sin.* = trig.__cosdf(sc1pio2 - x); 74 r_cos.* = trig.__sindf(sc1pio2 - x); 75 } 76 return; 77 } 78 // -sin(x+c) is not correct if x+c could be 0: -0 vs +0 79 r_sin.* = -trig.__sindf(if (sign) x + sc2pio2 else x - sc2pio2); 80 r_cos.* = -trig.__cosdf(if (sign) x + sc2pio2 else x - sc2pio2); 81 return; 82 } 83 84 // |x| ~<= 9*pi/4 85 if (ix <= 0x40e231d5) { 86 // |x| ~<= 7*pi/4 87 if (ix <= 0x40afeddf) { 88 if (sign) { 89 r_sin.* = trig.__cosdf(x + sc3pio2); 90 r_cos.* = -trig.__sindf(x + sc3pio2); 91 } else { 92 r_sin.* = -trig.__cosdf(x - sc3pio2); 93 r_cos.* = trig.__sindf(x - sc3pio2); 94 } 95 return; 96 } 97 r_sin.* = trig.__sindf(if (sign) x + sc4pio2 else x - sc4pio2); 98 r_cos.* = trig.__cosdf(if (sign) x + sc4pio2 else x - sc4pio2); 99 return; 100 } 101 102 // sin(Inf or NaN) is NaN 103 if (ix >= 0x7f800000) { 104 const result = x - x; 105 r_sin.* = result; 106 r_cos.* = result; 107 return; 108 } 109 110 // general argument reduction needed 111 var y: f64 = undefined; 112 const n = rem_pio2f(x, &y); 113 const s = trig.__sindf(y); 114 const c = trig.__cosdf(y); 115 switch (n & 3) { 116 0 => { 117 r_sin.* = s; 118 r_cos.* = c; 119 }, 120 1 => { 121 r_sin.* = c; 122 r_cos.* = -s; 123 }, 124 2 => { 125 r_sin.* = -s; 126 r_cos.* = -c; 127 }, 128 else => { 129 r_sin.* = -c; 130 r_cos.* = s; 131 }, 132 } 133 } 134 135 pub fn sincos(x: f64, r_sin: *f64, r_cos: *f64) callconv(.c) void { 136 const ix = @as(u32, @truncate(@as(u64, @bitCast(x)) >> 32)) & 0x7fffffff; 137 138 // |x| ~< pi/4 139 if (ix <= 0x3fe921fb) { 140 // if |x| < 2**-27 * sqrt(2) 141 if (ix < 0x3e46a09e) { 142 // raise inexact if x != 0 and underflow if subnormal 143 if (common.want_float_exceptions) { 144 if (ix < 0x00100000) { 145 mem.doNotOptimizeAway(x / 0x1p120); 146 } else { 147 mem.doNotOptimizeAway(x + 0x1p120); 148 } 149 } 150 r_sin.* = x; 151 r_cos.* = 1.0; 152 return; 153 } 154 r_sin.* = trig.__sin(x, 0.0, 0); 155 r_cos.* = trig.__cos(x, 0.0); 156 return; 157 } 158 159 // sincos(Inf or NaN) is NaN 160 if (ix >= 0x7ff00000) { 161 const result = x - x; 162 r_sin.* = result; 163 r_cos.* = result; 164 return; 165 } 166 167 // argument reduction needed 168 var y: [2]f64 = undefined; 169 const n = rem_pio2(x, &y); 170 const s = trig.__sin(y[0], y[1], 1); 171 const c = trig.__cos(y[0], y[1]); 172 switch (n & 3) { 173 0 => { 174 r_sin.* = s; 175 r_cos.* = c; 176 }, 177 1 => { 178 r_sin.* = c; 179 r_cos.* = -s; 180 }, 181 2 => { 182 r_sin.* = -s; 183 r_cos.* = -c; 184 }, 185 else => { 186 r_sin.* = -c; 187 r_cos.* = s; 188 }, 189 } 190 } 191 192 pub fn __sincosx(x: f80, r_sin: *f80, r_cos: *f80) callconv(.c) void { 193 // TODO: more efficient implementation 194 //return sincos_generic(f80, x, r_sin, r_cos); 195 var big_sin: f128 = undefined; 196 var big_cos: f128 = undefined; 197 sincosq(x, &big_sin, &big_cos); 198 r_sin.* = @as(f80, @floatCast(big_sin)); 199 r_cos.* = @as(f80, @floatCast(big_cos)); 200 } 201 202 pub fn sincosq(x: f128, r_sin: *f128, r_cos: *f128) callconv(.c) void { 203 // TODO: more correct implementation 204 //return sincos_generic(f128, x, r_sin, r_cos); 205 var small_sin: f64 = undefined; 206 var small_cos: f64 = undefined; 207 sincos(@as(f64, @floatCast(x)), &small_sin, &small_cos); 208 r_sin.* = small_sin; 209 r_cos.* = small_cos; 210 } 211 212 pub fn sincosl(x: c_longdouble, r_sin: *c_longdouble, r_cos: *c_longdouble) callconv(.c) void { 213 switch (@typeInfo(c_longdouble).float.bits) { 214 16 => return __sincosh(x, r_sin, r_cos), 215 32 => return sincosf(x, r_sin, r_cos), 216 64 => return sincos(x, r_sin, r_cos), 217 80 => return __sincosx(x, r_sin, r_cos), 218 128 => return sincosq(x, r_sin, r_cos), 219 else => @compileError("unreachable"), 220 } 221 } 222 223 pub const rem_pio2_generic = @compileError("TODO"); 224 225 /// Ported from musl sincosl.c. Needs the following dependencies to be complete: 226 /// * rem_pio2_generic ported from __rem_pio2l.c 227 /// * trig.sin_generic ported from __sinl.c 228 /// * trig.cos_generic ported from __cosl.c 229 inline fn sincos_generic(comptime F: type, x: F, r_sin: *F, r_cos: *F) void { 230 const sc1pio4: F = 1.0 * math.pi / 4.0; 231 const bits = @typeInfo(F).float.bits; 232 const I = std.meta.Int(.unsigned, bits); 233 const ix = @as(I, @bitCast(x)) & (math.maxInt(I) >> 1); 234 const se: u16 = @truncate(ix >> (bits - 16)); 235 236 if (se == 0x7fff) { 237 const result = x - x; 238 r_sin.* = result; 239 r_cos.* = result; 240 return; 241 } 242 243 if (@as(F, @bitCast(ix)) < sc1pio4) { 244 if (se < 0x3fff - math.floatFractionalBits(F) - 1) { 245 // raise underflow if subnormal 246 if (se == 0) { 247 if (common.want_float_exceptions) mem.doNotOptimizeAway(x * 0x1p-120); 248 } 249 r_sin.* = x; 250 // raise inexact if x!=0 251 r_cos.* = 1.0 + x; 252 return; 253 } 254 r_sin.* = trig.sin_generic(F, x, 0, 0); 255 r_cos.* = trig.cos_generic(F, x, 0); 256 return; 257 } 258 259 var y: [2]F = undefined; 260 const n = rem_pio2_generic(F, x, &y); 261 const s = trig.sin_generic(F, y[0], y[1], 1); 262 const c = trig.cos_generic(F, y[0], y[1]); 263 switch (n & 3) { 264 0 => { 265 r_sin.* = s; 266 r_cos.* = c; 267 }, 268 1 => { 269 r_sin.* = c; 270 r_cos.* = -s; 271 }, 272 2 => { 273 r_sin.* = -s; 274 r_cos.* = -c; 275 }, 276 else => { 277 r_sin.* = -c; 278 r_cos.* = s; 279 }, 280 } 281 }