1 module math; 2 3 import options; 4 5 import std.traits: Select, isFloatingPoint, isIntegral; 6 7 import std.math: 8 math_pi = PI, 9 math_sin = sin, 10 math_rint = rint, 11 math_floor = floor, 12 math_cos = cos, 13 math_acos = acos, 14 math_sqrt = sqrt, 15 math_atan2 = atan2, 16 math_tan = tan, 17 math_asin = asin, 18 math_abs = abs, 19 math_ceil = ceil, 20 math_round = round, 21 math_exp = exp, 22 math_fma = fma, 23 math_isInfinite = isInfinity; 24 25 // These aren't math library but oh well 26 import std.algorithm: 27 math_min = min, 28 math_max = max; 29 30 // I'm not indent matching this 31 import std.random: 32 math_random = Random, 33 math_unpredictableSeed = unpredictableSeed, 34 math_uniform = uniform; 35 36 import rounding_mode; 37 38 // All the java comments are now wrong oh nooo 39 40 41 /* 42 * The MIT License 43 * 44 * Copyright (c) 2015-2021 DOML 45 +*%#$%# Translated by jordan4ibanez 46 * 47 * Permission is hereby granted, free of charge, to any person obtaining a copy 48 * of this software and associated documentation files (the "Software"), to deal 49 * in the Software without restriction, including without limitation the rights 50 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell 51 * copies of the Software, and to permit persons to whom the Software is 52 * furnished to do so, subject to the following conditions: 53 * 54 * The above copyright notice and this permission notice shall be included in 55 * all copies or substantial portions of the Software. 56 * 57 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR 58 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 59 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE 60 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER 61 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, 62 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN 63 * THE SOFTWARE. 64 */ 65 66 /** 67 * Contains fast approximations of some {@link java.lang.Math} operations. 68 * <p> 69 * By default, {@link java.lang.Math} methods will be used by all other DOML classes. In order to use the approximations in this class, start the JVM with the parameter <code>-DDOML.fastmath</code>. 70 * <p> 71 * There are two algorithms for approximating sin/cos: 72 * <ol> 73 * <li>arithmetic <a href="http://www.java-gaming.org/topics/DOML-1-8-0-release/37491/msg/361815/view.html#msg361815">polynomial approximation</a> contributed by roquendm 74 * <li>theagentd's <a href="http://www.java-gaming.org/topics/extremely-fast-sine-cosine/36469/msg/346213/view.html#msg346213">linear interpolation</a> variant of Riven's algorithm from 75 * <a href="http://www.java-gaming.org/topics/extremely-fast-sine-cosine/36469/view.html">http://www.java-gaming.org/</a> 76 * </ol> 77 * By default, the first algorithm is being used. In order to use the second one, start the JVM with <code>-DDOML.sinLookup</code>. The lookup table bit length of the second algorithm can also be adjusted 78 * for improved accuracy via <code>-DDOML.sinLookup.bits=<n></code>, where <n> is the number of bits of the lookup table. 79 * 80 * @author Kai Burjack 81 */ 82 83 /* 84 * The following implementation of an approximation of sine and cosine was 85 * thankfully donated by Riven from http://java-gaming.org/. 86 * 87 * The code for linear interpolation was gratefully donated by theagentd 88 * from the same site. 89 */ 90 91 immutable double PI = math_pi; 92 immutable double PI2 = PI * 2.0; 93 immutable double PIHalf = PI * 0.5; 94 immutable double PI_4 = PI * 0.25; 95 immutable double PI_INV = 1.0 / PI; 96 97 98 // Credit: https://forum.dlang.org/post/bug-5900-3@http.d.puremagic.com%2Fissues%2F 99 /// Converts from degrees to radians. 100 @safe pure nothrow Select!(isFloatingPoint!T || isIntegral!T, T, double) 101 radians(T)(in T x) { 102 return x * (PI / 180); 103 } 104 105 /// Converts from radians to degrees. 106 @safe pure nothrow Select!(isFloatingPoint!T || isIntegral!T, T, double) 107 degrees(T)(in T x) { 108 return x / (PI / 180); 109 } 110 111 @safe pure nothrow Select!(isFloatingPoint!T || isIntegral!T, T, double) 112 math_signum(T)(in T x) { 113 return (T(0) < x) - (x < T(0)); 114 } 115 116 // Thanks to adr, a Java function in D 117 double longbits(long a) { 118 return *cast(double*) &a; 119 } 120 int floatToRawIntBits(float a) { 121 return *cast(int*) &a; 122 } 123 long doubleToLongBits(double a) { 124 return *cast(long*) &a; 125 } 126 127 bool isInfinite(double value) { 128 return math_isInfinite(value); 129 } 130 131 132 private const double c1 = longbits(-4_628_199_217_061_079_772L); 133 private const double c2 = longbits(4_575_957_461_383_582_011L); 134 private const double c3 = longbits(-4_671_919_876_300_759_001L); 135 private const double c4 = longbits(4_523_617_214_285_661_942L); 136 private const double c5 = longbits(-4_730_215_272_828_025_532L); 137 private const double c6 = longbits(4_460_272_573_143_870_633L); 138 private const double c7 = longbits(-4_797_767_418_267_846_529L); 139 140 /** 141 * @author theagentd 142 */ 143 double sin_theagentd_arith(double x){ 144 double xi = math_floor((x + PI_4) * PI_INV); 145 double x_ = x - xi * PI; 146 double sign = (cast(int)xi & 1) * -2 + 1; 147 double x2 = x_ * x_; 148 double sin = x_; 149 double tx = x_ * x2; 150 sin += tx * c1; tx *= x2; 151 sin += tx * c2; tx *= x2; 152 sin += tx * c3; tx *= x2; 153 sin += tx * c4; tx *= x2; 154 sin += tx * c5; tx *= x2; 155 sin += tx * c6; tx *= x2; 156 sin += tx * c7; 157 return sign * sin; 158 } 159 160 /** 161 * Reference: <a href="http://www.java-gaming.org/topics/DOML-1-8-0-release/37491/msg/361718/view.html#msg361718">http://www.java-gaming.org/</a> 162 */ 163 double sin_roquen_arith(double x) { 164 double xi = math_floor((x + PI_4) * PI_INV); 165 double x_ = x - xi * PI; 166 double sign = (cast(int)xi & 1) * -2 + 1; 167 double x2 = x_ * x_; 168 169 // code from sin_theagentd_arith: 170 // double sin = x_; 171 // double tx = x_ * x2; 172 // sin += tx * c1; tx *= x2; 173 // sin += tx * c2; tx *= x2; 174 // sin += tx * c3; tx *= x2; 175 // sin += tx * c4; tx *= x2; 176 // sin += tx * c5; tx *= x2; 177 // sin += tx * c6; tx *= x2; 178 // sin += tx * c7; 179 // return sign * sin; 180 181 double sin; 182 x_ = sign*x_; 183 sin = c7; 184 sin = sin*x2 + c6; 185 sin = sin*x2 + c5; 186 sin = sin*x2 + c4; 187 sin = sin*x2 + c3; 188 sin = sin*x2 + c2; 189 sin = sin*x2 + c1; 190 return x_ + x_*x2*sin; 191 } 192 193 private const double s5 = longbits(4_523_227_044_276_562_163L); 194 private const double s4 = longbits(-4_671_934_770_969_572_232L); 195 private const double s3 = longbits(4_575_957_211_482_072_852L); 196 private const double s2 = longbits(-4_628_199_223_918_090_387L); 197 private const double s1 = longbits(4_607_182_418_589_157_889L); 198 199 /** 200 * Reference: <a href="http://www.java-gaming.org/topics/DOML-1-8-0-release/37491/msg/361815/view.html#msg361815">http://www.java-gaming.org/</a> 201 */ 202 double sin_roquen_9(double v) { 203 double i = math_rint(v*PI_INV); 204 double x = v - i * PI; 205 double qs = 1-2*(cast(int)i & 1); 206 double x2 = x*x; 207 double r; 208 x = qs*x; 209 r = s5; 210 r = r*x2 + s4; 211 r = r*x2 + s3; 212 r = r*x2 + s2; 213 r = r*x2 + s1; 214 return x*r; 215 } 216 217 private const double k1 = longbits(-4_628_199_217_061_079_959L); 218 private const double k2 = longbits(4_575_957_461_383_549_981L); 219 private const double k3 = longbits(-4_671_919_876_307_284_301L); 220 private const double k4 = longbits(4_523_617_213_632_129_738L); 221 private const double k5 = longbits(-4_730_215_344_060_517_252L); 222 private const double k6 = longbits(4_460_268_259_291_226_124L); 223 private const double k7 = longbits(-4_798_040_743_777_455_072L); 224 225 /** 226 * Reference: <a href="http://www.java-gaming.org/topics/DOML-1-8-0-release/37491/msg/361815/view.html#msg361815">http://www.java-gaming.org/</a> 227 */ 228 double sin_roquen_newk(double v) { 229 double i = math_rint(v*PI_INV); 230 double x = v - i * PI; 231 double qs = 1-2*(cast(int)i & 1); 232 double x2 = x*x; 233 double r; 234 x = qs*x; 235 r = k7; 236 r = r*x2 + k6; 237 r = r*x2 + k5; 238 r = r*x2 + k4; 239 r = r*x2 + k3; 240 r = r*x2 + k2; 241 r = r*x2 + k1; 242 return x + x*x2*r; 243 } 244 245 double sin(double rad) { 246 return math_sin(rad); 247 } 248 249 double cos(double rad) { 250 return math_cos(rad); 251 } 252 253 double cosFromSin(double sin, double angle) { 254 //if (Options.FASTMATH){ 255 // return math_sin(angle + PIHalf); 256 // } 257 // sin(x)^2 + cos(x)^2 = 1 258 double cos = math_sqrt(1.0 - sin * sin); 259 double a = angle + PIHalf; 260 double b = a - cast(int)(a / PI2) * PI2; 261 if (b < 0.0) 262 b = PI2 + b; 263 if (b >= PI) 264 return -cos; 265 return cos; 266 } 267 268 /* Other math functions not yet approximated */ 269 270 double sqrt(double r) { 271 return math_sqrt(r); 272 } 273 274 double invsqrt(double r) { 275 return 1.0 / math_sqrt(r); 276 } 277 278 double tan(double r) { 279 return math_tan(r); 280 } 281 282 double acos(double r) { 283 return math_acos(r); 284 } 285 286 double safeAcos(double v) { 287 if (v < -1.0) 288 return PI; 289 else if (v > +1.0) 290 return 0.0; 291 else 292 return math_acos(v); 293 } 294 295 296 double atan2(double y, double x) { 297 return math_atan2(y, x); 298 } 299 300 double asin(double r) { 301 return math_asin(r); 302 } 303 304 double abs(double r) { 305 return math_abs(r); 306 } 307 308 // Taken from runtime 309 bool equals(double a, double b, double delta) { 310 return doubleToLongBits(a) == doubleToLongBits(b) || abs(a - b) <= delta; 311 } 312 313 bool absEqualsOne(double r) { 314 return (doubleToLongBits(r) & 0x7FFFFFFFFFFFFFFFL) == 0x3FF0000000000000L; 315 } 316 317 int abs(int r) { 318 return math_abs(r); 319 } 320 321 int max(int x, int y) { 322 return math_max(x, y); 323 } 324 325 int min(int x, int y) { 326 return math_min(x, y); 327 } 328 329 double min(double a, double b) { 330 return a < b ? a : b; 331 } 332 333 double max(double a, double b) { 334 return a > b ? a : b; 335 } 336 337 ulong max(ulong a, ulong b) { 338 return a > b ? a : b; 339 } 340 341 double clamp(double a, double b, double val) { 342 return max(a,math_min(b,val)); 343 } 344 int clamp(int a, int b, int val) { 345 return max(a, math_min(b, val)); 346 } 347 348 349 double toRadians(double angles) { 350 return radians(angles); 351 } 352 353 double toDegrees(double angles) { 354 return degrees(angles); 355 } 356 357 double floor(double v) { 358 return math_floor(v); 359 } 360 361 double ceil(double v) { 362 return math_ceil(v); 363 } 364 365 long round(double v) { 366 return cast(long) math_round(v); 367 } 368 369 double exp(double a) { 370 return math_exp(a); 371 } 372 373 bool isFinite(double d) { 374 return math_abs(d) <= double.max; 375 } 376 377 double fma(double a, double b, double c) { 378 return math_fma(a, b, c); 379 } 380 381 int roundUsing(double v, int mode) { 382 switch (mode) { 383 case RoundingMode.TRUNCATE: 384 return cast(int) v; 385 case RoundingMode.CEILING: 386 return cast(int) math_ceil(v); 387 case RoundingMode.FLOOR: 388 return cast(int) math_floor(v); 389 case RoundingMode.HALF_DOWN: 390 return roundHalfDown(v); 391 case RoundingMode.HALF_UP: 392 return roundHalfUp(v); 393 case RoundingMode.HALF_EVEN: 394 return roundHalfEven(v); 395 default: { 396 return 0; 397 } 398 } 399 } 400 401 double lerp(double a, double b, double t) { 402 return math_fma(b - a, t, a); 403 } 404 405 double biLerp(double q00, double q10, double q01, double q11, double tx, double ty) { 406 double lerpX1 = lerp(q00, q10, tx); 407 double lerpX2 = lerp(q01, q11, tx); 408 return lerp(lerpX1, lerpX2, ty); 409 } 410 411 double triLerp(double q000, double q100, double q010, double q110, double q001, double q101, double q011, double q111, 412 double tx, double ty, double tz) { 413 double x00 = lerp(q000, q100, tx); 414 double x10 = lerp(q010, q110, tx); 415 double x01 = lerp(q001, q101, tx); 416 double x11 = lerp(q011, q111, tx); 417 double y0 = lerp(x00, x10, ty); 418 double y1 = lerp(x01, x11, ty); 419 return lerp(y0, y1, tz); 420 } 421 422 int roundHalfEven(double v) { 423 return cast(int) math_rint(v); 424 } 425 int roundHalfDown(double v) { 426 return (v > 0) ? cast(int) math_ceil(v - 0.5) : cast(int) math_floor(v + 0.5); 427 } 428 int roundHalfUp(double v) { 429 return (v > 0) ? cast(int) math_floor(v + 0.5) : cast(int) math_ceil(v - 0.5); 430 } 431 432 double random() { 433 math_random random = math_random(math_unpredictableSeed()); 434 return math_uniform(0.0, 1.0, random); 435 } 436 437 double safeAsin(double r) { 438 return r <= -1.0 ? -PIHalf : r >= 1.0 ? PIHalf : math_asin(r); 439 } 440 441 double signum(double v) { 442 return math_signum(v); 443 } 444 int signum(int v) { 445 int r; 446 r = math_signum(v); 447 r = (v >> 31) | (-v >>> 31); 448 return r; 449 } 450 int signum(long v) { 451 int r; 452 r = cast(int)math_signum(v); 453 r = cast(int) ((v >> 63) | (-v >>> 63)); 454 return r; 455 } 456