blob c315adc4 (6812B) - Raw
1 // Special Cases: 2 // 3 // - atan(+-0) = +-0 4 // - atan(+-inf) = +-pi/2 5 6 const std = @import("../index.zig"); 7 const math = std.math; 8 const assert = std.debug.assert; 9 10 pub fn atan(x: var) @typeOf(x) { 11 const T = @typeOf(x); 12 return switch (T) { 13 f32 => atan32(x), 14 f64 => atan64(x), 15 else => @compileError("atan not implemented for " ++ @typeName(T)), 16 }; 17 } 18 19 fn atan32(x_: f32) f32 { 20 const atanhi = []const f32 { 21 4.6364760399e-01, // atan(0.5)hi 22 7.8539812565e-01, // atan(1.0)hi 23 9.8279368877e-01, // atan(1.5)hi 24 1.5707962513e+00, // atan(inf)hi 25 }; 26 27 const atanlo = []const f32 { 28 5.0121582440e-09, // atan(0.5)lo 29 3.7748947079e-08, // atan(1.0)lo 30 3.4473217170e-08, // atan(1.5)lo 31 7.5497894159e-08, // atan(inf)lo 32 }; 33 34 const aT = []const f32 { 35 3.3333328366e-01, 36 -1.9999158382e-01, 37 1.4253635705e-01, 38 -1.0648017377e-01, 39 6.1687607318e-02, 40 }; 41 42 var x = x_; 43 var ix: u32 = @bitCast(u32, x); 44 const sign = ix >> 31; 45 ix &= 0x7FFFFFFF; 46 47 // |x| >= 2^26 48 if (ix >= 0x4C800000) { 49 if (math.isNan(x)) { 50 return x; 51 } else { 52 const z = atanhi[3] + 0x1.0p-120; 53 return if (sign != 0) -z else z; 54 } 55 } 56 57 var id: ?usize = undefined; 58 59 // |x| < 0.4375 60 if (ix < 0x3EE00000) { 61 // |x| < 2^(-12) 62 if (ix < 0x39800000) { 63 if (ix < 0x00800000) { 64 math.forceEval(x * x); 65 } 66 return x; 67 } 68 id = null; 69 } else { 70 x = math.fabs(x); 71 // |x| < 1.1875 72 if (ix < 0x3F980000) { 73 // 7/16 <= |x| < 11/16 74 if (ix < 0x3F300000) { 75 id = 0; 76 x = (2.0 * x - 1.0) / (2.0 + x); 77 } 78 // 11/16 <= |x| < 19/16 79 else { 80 id = 1; 81 x = (x - 1.0) / (x + 1.0); 82 } 83 } 84 else { 85 // |x| < 2.4375 86 if (ix < 0x401C0000) { 87 id = 2; 88 x = (x - 1.5) / (1.0 + 1.5 * x); 89 } 90 // 2.4375 <= |x| < 2^26 91 else { 92 id = 3; 93 x = -1.0 / x; 94 } 95 } 96 } 97 98 const z = x * x; 99 const w = z * z; 100 const s1 = z * (aT[0] + w * (aT[2] + w * aT[4])); 101 const s2 = w * (aT[1] + w * aT[3]); 102 103 if (id) |id_value| { 104 const zz = atanhi[id_value] - ((x * (s1 + s2) - atanlo[id_value]) - x); 105 return if (sign != 0) -zz else zz; 106 } else { 107 return x - x * (s1 + s2); 108 } 109 } 110 111 fn atan64(x_: f64) f64 { 112 const atanhi = []const f64 { 113 4.63647609000806093515e-01, // atan(0.5)hi 114 7.85398163397448278999e-01, // atan(1.0)hi 115 9.82793723247329054082e-01, // atan(1.5)hi 116 1.57079632679489655800e+00, // atan(inf)hi 117 }; 118 119 const atanlo = []const f64 { 120 2.26987774529616870924e-17, // atan(0.5)lo 121 3.06161699786838301793e-17, // atan(1.0)lo 122 1.39033110312309984516e-17, // atan(1.5)lo 123 6.12323399573676603587e-17, // atan(inf)lo 124 }; 125 126 const aT = []const f64 { 127 3.33333333333329318027e-01, 128 -1.99999999998764832476e-01, 129 1.42857142725034663711e-01, 130 -1.11111104054623557880e-01, 131 9.09088713343650656196e-02, 132 -7.69187620504482999495e-02, 133 6.66107313738753120669e-02, 134 -5.83357013379057348645e-02, 135 4.97687799461593236017e-02, 136 -3.65315727442169155270e-02, 137 1.62858201153657823623e-02, 138 }; 139 140 var x = x_; 141 var ux = @bitCast(u64, x); 142 var ix = u32(ux >> 32); 143 const sign = ix >> 31; 144 ix &= 0x7FFFFFFF; 145 146 // |x| >= 2^66 147 if (ix >= 0x44100000) { 148 if (math.isNan(x)) { 149 return x; 150 } else { 151 const z = atanhi[3] + 0x1.0p-120; 152 return if (sign != 0) -z else z; 153 } 154 } 155 156 var id: ?usize = undefined; 157 158 // |x| < 0.4375 159 if (ix < 0x3DFC0000) { 160 // |x| < 2^(-27) 161 if (ix < 0x3E400000) { 162 if (ix < 0x00100000) { 163 math.forceEval(f32(x)); 164 } 165 return x; 166 } 167 id = null; 168 } else { 169 x = math.fabs(x); 170 // |x| < 1.1875 171 if (ix < 0x3FF30000) { 172 // 7/16 <= |x| < 11/16 173 if (ix < 0x3FE60000) { 174 id = 0; 175 x = (2.0 * x - 1.0) / (2.0 + x); 176 } 177 // 11/16 <= |x| < 19/16 178 else { 179 id = 1; 180 x = (x - 1.0) / (x + 1.0); 181 } 182 } 183 else { 184 // |x| < 2.4375 185 if (ix < 0x40038000) { 186 id = 2; 187 x = (x - 1.5) / (1.0 + 1.5 * x); 188 } 189 // 2.4375 <= |x| < 2^66 190 else { 191 id = 3; 192 x = -1.0 / x; 193 } 194 } 195 } 196 197 const z = x * x; 198 const w = z * z; 199 const s1 = z * (aT[0] + w * (aT[2] + w * (aT[4] + w * (aT[6] + w * (aT[8] + w * aT[10]))))); 200 const s2 = w * (aT[1] + w * (aT[3] + w * (aT[5] + w * (aT[7] + w * aT[9])))); 201 202 if (id) |id_value| { 203 const zz = atanhi[id_value] - ((x * (s1 + s2) - atanlo[id_value]) - x); 204 return if (sign != 0) -zz else zz; 205 } else { 206 return x - x * (s1 + s2); 207 } 208 } 209 210 test "math.atan" { 211 assert(@bitCast(u32, atan(f32(0.2))) == @bitCast(u32, atan32(0.2))); 212 assert(atan(f64(0.2)) == atan64(0.2)); 213 } 214 215 test "math.atan32" { 216 const epsilon = 0.000001; 217 218 assert(math.approxEq(f32, atan32(0.2), 0.197396, epsilon)); 219 assert(math.approxEq(f32, atan32(-0.2), -0.197396, epsilon)); 220 assert(math.approxEq(f32, atan32(0.3434), 0.330783, epsilon)); 221 assert(math.approxEq(f32, atan32(0.8923), 0.728545, epsilon)); 222 assert(math.approxEq(f32, atan32(1.5), 0.982794, epsilon)); 223 } 224 225 test "math.atan64" { 226 const epsilon = 0.000001; 227 228 assert(math.approxEq(f64, atan64(0.2), 0.197396, epsilon)); 229 assert(math.approxEq(f64, atan64(-0.2), -0.197396, epsilon)); 230 assert(math.approxEq(f64, atan64(0.3434), 0.330783, epsilon)); 231 assert(math.approxEq(f64, atan64(0.8923), 0.728545, epsilon)); 232 assert(math.approxEq(f64, atan64(1.5), 0.982794, epsilon)); 233 } 234 235 test "math.atan32.special" { 236 const epsilon = 0.000001; 237 238 assert(atan32(0.0) == 0.0); 239 assert(atan32(-0.0) == -0.0); 240 assert(math.approxEq(f32, atan32(math.inf(f32)), math.pi / 2.0, epsilon)); 241 assert(math.approxEq(f32, atan32(-math.inf(f32)), -math.pi / 2.0, epsilon)); 242 } 243 244 test "math.atan64.special" { 245 const epsilon = 0.000001; 246 247 assert(atan64(0.0) == 0.0); 248 assert(atan64(-0.0) == -0.0); 249 assert(math.approxEq(f64, atan64(math.inf(f64)), math.pi / 2.0, epsilon)); 250 assert(math.approxEq(f64, atan64(-math.inf(f64)), -math.pi / 2.0, epsilon)); 251 }