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