1 /++ 2 Common floating point math functions. 3 4 This module has generic LLVM-oriented API compatible with all D compilers. 5 6 License: $(HTTP www.apache.org/licenses/LICENSE-2.0, Apache-2.0) 7 Authors: Ilya Yaroshenko, Phobos Team 8 +/ 9 module mir.math.common; 10 11 import mir.internal.utility: isComplex, isComplexOf, isFloatingPoint; 12 13 version(LDC) 14 { 15 static import ldc.attributes; 16 17 private alias AliasSeq(T...) = T; 18 19 /++ 20 Functions attribute, an alias for `AliasSeq!(llvmFastMathFlag("contract"));`. 21 22 $(UL 23 $(LI 1. Allow floating-point contraction (e.g. fusing a multiply followed by an addition into a fused multiply-and-add). ) 24 ) 25 26 Note: Can be used with all compilers. 27 +/ 28 alias fmamath = AliasSeq!(ldc.attributes.llvmFastMathFlag("contract")); 29 30 /++ 31 Functions attribute, an alias for `AliasSeq!(llvmFastMathFlag("fast"))`. 32 33 It is similar to $(LREF fastmath), but does not allow unsafe-fp-math. 34 This flag does NOT force LDC to use the reciprocal of an argument rather than perform division. 35 36 This flag is default for string lambdas. 37 38 Note: Can be used with all compilers. 39 +/ 40 alias optmath = AliasSeq!(ldc.attributes.llvmFastMathFlag("fast")); 41 42 /++ 43 Functions attribute, an alias for `ldc.attributes.fastmath` . 44 45 $(UL 46 47 $(LI 1. Enable optimizations that make unsafe assumptions about IEEE math (e.g. that addition is associative) or may not work for all input ranges. 48 These optimizations allow the code generator to make use of some instructions which would otherwise not be usable (such as fsin on X86). ) 49 50 $(LI 2. Allow optimizations to assume the arguments and result are not NaN. 51 Such optimizations are required to retain defined behavior over NaNs, 52 but the value of the result is undefined. ) 53 54 $(LI 3. Allow optimizations to assume the arguments and result are not +$(BACKTICK)-inf. 55 Such optimizations are required to retain defined behavior over +$(BACKTICK)-Inf, 56 but the value of the result is undefined. ) 57 58 $(LI 4. Allow optimizations to treat the sign of a zero argument or result as insignificant. ) 59 60 $(LI 5. Allow optimizations to use the reciprocal of an argument rather than perform division. ) 61 62 $(LI 6. Allow floating-point contraction (e.g. fusing a multiply followed by an addition into a fused multiply-and-add). ) 63 64 $(LI 7. Allow algebraically equivalent transformations that may dramatically change results in floating point (e.g. reassociate). ) 65 ) 66 67 Note: Can be used with all compilers. 68 +/ 69 alias fastmath = ldc.attributes.fastmath; 70 } 71 else 72 enum 73 { 74 /++ 75 Functions attribute, an alias for `AliasSeq!(llvmFastMathFlag("contract"));`. 76 77 $(UL 78 $(LI Allow floating-point contraction (e.g. fusing a multiply followed by an addition into a fused multiply-and-add). ) 79 ) 80 81 Note: Can be used with all compilers. 82 +/ 83 fmamath, 84 85 /++ 86 Functions attribute, an alias for `AliasSeq!(llvmAttr("unsafe-fp-math", "false"), llvmFastMathFlag("fast"))`. 87 88 It is similar to $(LREF fastmath), but does not allow unsafe-fp-math. 89 This flag does NOT force LDC to use the reciprocal of an argument rather than perform division. 90 91 This flag is default for string lambdas. 92 93 Note: Can be used with all compilers. 94 +/ 95 optmath, 96 97 /++ 98 Functions attribute, an alias for `ldc.attributes.fastmath = AliasSeq!(llvmAttr("unsafe-fp-math", "true"), llvmFastMathFlag("fast"))` . 99 100 $(UL 101 102 $(LI Enable optimizations that make unsafe assumptions about IEEE math (e.g. that addition is associative) or may not work for all input ranges. 103 These optimizations allow the code generator to make use of some instructions which would otherwise not be usable (such as fsin on X86). ) 104 105 $(LI Allow optimizations to assume the arguments and result are not NaN. 106 Such optimizations are required to retain defined behavior over NaNs, 107 but the value of the result is undefined. ) 108 109 $(LI Allow optimizations to assume the arguments and result are not +$(BACKTICK)-inf. 110 Such optimizations are required to retain defined behavior over +$(BACKTICK)-Inf, 111 but the value of the result is undefined. ) 112 113 $(LI Allow optimizations to treat the sign of a zero argument or result as insignificant. ) 114 115 $(LI Allow optimizations to use the reciprocal of an argument rather than perform division. ) 116 117 $(LI Allow floating-point contraction (e.g. fusing a multiply followed by an addition into a fused multiply-and-add). ) 118 119 $(LI Allow algebraically equivalent transformations that may dramatically change results in floating point (e.g. reassociate). ) 120 ) 121 122 Note: Can be used with all compilers. 123 +/ 124 fastmath 125 } 126 127 version(LDC) 128 { 129 nothrow @nogc pure @safe: 130 131 pragma(LDC_intrinsic, "llvm.sqrt.f#") 132 /// 133 T sqrt(T)(in T val) if (isFloatingPoint!T); 134 135 pragma(LDC_intrinsic, "llvm.sin.f#") 136 /// 137 T sin(T)(in T val) if (isFloatingPoint!T); 138 139 pragma(LDC_intrinsic, "llvm.cos.f#") 140 /// 141 T cos(T)(in T val) if (isFloatingPoint!T); 142 143 pragma(LDC_intrinsic, "llvm.powi.f#") 144 /// 145 T powi(T)(in T val, int power) if (isFloatingPoint!T); 146 147 pragma(LDC_intrinsic, "llvm.pow.f#") 148 /// 149 T pow(T)(in T val, in T power) if (isFloatingPoint!T); 150 151 pragma(LDC_intrinsic, "llvm.exp.f#") 152 /// 153 T exp(T)(in T val) if (isFloatingPoint!T); 154 155 pragma(LDC_intrinsic, "llvm.log.f#") 156 /// 157 T log(T)(in T val) if (isFloatingPoint!T); 158 159 pragma(LDC_intrinsic, "llvm.fma.f#") 160 /// 161 T fma(T)(T vala, T valb, T valc) if (isFloatingPoint!T); 162 163 pragma(LDC_intrinsic, "llvm.fabs.f#") 164 /// 165 T fabs(T)(in T val) if (isFloatingPoint!T); 166 167 pragma(LDC_intrinsic, "llvm.floor.f#") 168 /// 169 T floor(T)(in T val) if (isFloatingPoint!T); 170 171 pragma(LDC_intrinsic, "llvm.exp2.f#") 172 /// 173 T exp2(T)(in T val) if (isFloatingPoint!T); 174 175 pragma(LDC_intrinsic, "llvm.log10.f#") 176 /// 177 T log10(T)(in T val) if (isFloatingPoint!T); 178 179 pragma(LDC_intrinsic, "llvm.log2.f#") 180 /// 181 T log2(T)(in T val) if (isFloatingPoint!T); 182 183 pragma(LDC_intrinsic, "llvm.ceil.f#") 184 /// 185 T ceil(T)(in T val) if (isFloatingPoint!T); 186 187 pragma(LDC_intrinsic, "llvm.trunc.f#") 188 /// 189 T trunc(T)(in T val) if (isFloatingPoint!T); 190 191 pragma(LDC_intrinsic, "llvm.rint.f#") 192 /// 193 T rint(T)(in T val) if (isFloatingPoint!T); 194 195 pragma(LDC_intrinsic, "llvm.nearbyint.f#") 196 /// 197 T nearbyint(T)(in T val) if (isFloatingPoint!T); 198 199 pragma(LDC_intrinsic, "llvm.copysign.f#") 200 /// 201 T copysign(T)(in T mag, in T sgn) if (isFloatingPoint!T); 202 203 pragma(LDC_intrinsic, "llvm.round.f#") 204 /// 205 T round(T)(in T val) if (isFloatingPoint!T); 206 207 pragma(LDC_intrinsic, "llvm.fmuladd.f#") 208 /// 209 T fmuladd(T)(in T vala, in T valb, in T valc) if (isFloatingPoint!T); 210 211 pragma(LDC_intrinsic, "llvm.minnum.f#") 212 /// 213 T fmin(T)(in T vala, in T valb) if (isFloatingPoint!T); 214 215 pragma(LDC_intrinsic, "llvm.maxnum.f#") 216 /// 217 T fmax(T)(in T vala, in T valb) if (isFloatingPoint!T); 218 } 219 else version(GNU) 220 { 221 static import gcc.builtins; 222 223 // Calls GCC builtin for either float (suffix "f"), double (no suffix), or real (suffix "l"). 224 private enum mixinGCCBuiltin(string fun) = 225 `static if (T.mant_dig == float.mant_dig) return gcc.builtins.__builtin_`~fun~`f(x);`~ 226 ` else static if (T.mant_dig == double.mant_dig) return gcc.builtins.__builtin_`~fun~`(x);`~ 227 ` else static if (T.mant_dig == real.mant_dig) return gcc.builtins.__builtin_`~fun~`l(x);`~ 228 ` else static assert(0);`; 229 230 // As above but for two-argument function. 231 private enum mixinGCCBuiltin2(string fun) = 232 `static if (T.mant_dig == float.mant_dig) return gcc.builtins.__builtin_`~fun~`f(x, y);`~ 233 ` else static if (T.mant_dig == double.mant_dig) return gcc.builtins.__builtin_`~fun~`(x, y);`~ 234 ` else static if (T.mant_dig == real.mant_dig) return gcc.builtins.__builtin_`~fun~`l(x, y);`~ 235 ` else static assert(0);`; 236 237 /// 238 T sqrt(T)(in T x) if (isFloatingPoint!T) { mixin(mixinGCCBuiltin!`sqrt`); } 239 /// 240 T sin(T)(in T x) if (isFloatingPoint!T) { mixin(mixinGCCBuiltin!`sin`); } 241 /// 242 T cos(T)(in T x) if (isFloatingPoint!T) { mixin(mixinGCCBuiltin!`cos`); } 243 /// 244 T pow(T)(in T x, in T power) if (isFloatingPoint!T) { alias y = power; mixin(mixinGCCBuiltin2!`pow`); } 245 /// 246 T powi(T)(in T x, int power) if (isFloatingPoint!T) { alias y = power; mixin(mixinGCCBuiltin2!`powi`); } 247 /// 248 T exp(T)(in T x) if (isFloatingPoint!T) { mixin(mixinGCCBuiltin!`exp`); } 249 /// 250 T log(T)(in T x) if (isFloatingPoint!T) { mixin(mixinGCCBuiltin!`log`); } 251 /// 252 T fabs(T)(in T x) if (isFloatingPoint!T) { mixin(mixinGCCBuiltin!`fabs`); } 253 /// 254 T floor(T)(in T x) if (isFloatingPoint!T) { mixin(mixinGCCBuiltin!`floor`); } 255 /// 256 T exp2(T)(in T x) if (isFloatingPoint!T) { mixin(mixinGCCBuiltin!`exp2`); } 257 /// 258 T log10(T)(in T x) if (isFloatingPoint!T) { mixin(mixinGCCBuiltin!`log10`); } 259 /// 260 T log2(T)(in T x) if (isFloatingPoint!T) { mixin(mixinGCCBuiltin!`log2`); } 261 /// 262 T ceil(T)(in T x) if (isFloatingPoint!T) { mixin(mixinGCCBuiltin!`ceil`); } 263 /// 264 T trunc(T)(in T x) if (isFloatingPoint!T) { mixin(mixinGCCBuiltin!`trunc`); } 265 /// 266 T rint(T)(in T x) if (isFloatingPoint!T) { mixin(mixinGCCBuiltin!`rint`); } 267 /// 268 T nearbyint(T)(in T x) if (isFloatingPoint!T) { mixin(mixinGCCBuiltin!`nearbyint`); } 269 /// 270 T copysign(T)(in T mag, in T sgn) if (isFloatingPoint!T) { alias y = sgn; mixin(mixinGCCBuiltin2!`copysign`); } 271 /// 272 T round(T)(in T x) if (isFloatingPoint!T) { mixin(mixinGCCBuiltin!`round`); } 273 /// 274 T fmuladd(T)(in T a, in T b, in T c) if (isFloatingPoint!T) 275 { 276 static if (T.mant_dig == float.mant_dig) 277 return gcc.builtins.__builtin_fmaf(a, b, c); 278 else static if (T.mant_dig == double.mant_dig) 279 return gcc.builtins.__builtin_fma(a, b, c); 280 else static if (T.mant_dig == real.mant_dig) 281 return gcc.builtins.__builtin_fmal(a, b, c); 282 else 283 static assert(0); 284 } 285 version(mir_core_test) 286 unittest { assert(fmuladd!double(2, 3, 4) == 2 * 3 + 4); } 287 /// 288 T fmin(T)(in T x, in T y) if (isFloatingPoint!T) { mixin(mixinGCCBuiltin2!`fmin`); } 289 /// 290 T fmax(T)(in T x, in T y) if (isFloatingPoint!T) { mixin(mixinGCCBuiltin2!`fmax`); } 291 } 292 else static if (__VERSION__ >= 2082) // DMD 2.082 onward. 293 { 294 static import std.math; 295 static import core.stdc.math; 296 297 // Calls either std.math or cmath function for either float (suffix "f") 298 // or double (no suffix). std.math will always be used during CTFE or for 299 // arguments with greater than double precision or if the cmath function 300 // is impure. 301 private enum mixinCMath(string fun) = 302 `pragma(inline, true); 303 static if (!is(typeof(std.math.`~fun~`(0.5f)) == float) 304 && is(typeof(() pure => core.stdc.math.`~fun~`f(0.5f)))) 305 if (!__ctfe) 306 { 307 static if (T.mant_dig == float.mant_dig) return core.stdc.math.`~fun~`f(x); 308 else static if (T.mant_dig == double.mant_dig) return core.stdc.math.`~fun~`(x); 309 } 310 return std.math.`~fun~`(x);`; 311 312 // As above but for two-argument function (both arguments must be floating point). 313 private enum mixinCMath2(string fun) = 314 `pragma(inline, true); 315 static if (!is(typeof(std.math.`~fun~`(0.5f, 0.5f)) == float) 316 && is(typeof(() pure => core.stdc.math.`~fun~`f(0.5f, 0.5f)))) 317 if (!__ctfe) 318 { 319 static if (T.mant_dig == float.mant_dig) return core.stdc.math.`~fun~`f(x, y); 320 else static if (T.mant_dig == double.mant_dig) return core.stdc.math.`~fun~`(x, y); 321 } 322 return std.math.`~fun~`(x, y);`; 323 324 // Some std.math functions have appropriate return types (float, 325 // double, real) without need for a wrapper. We can alias them 326 // directly but we leave the templates afterwards for documentation 327 // purposes and so explicit template instantiation still works. 328 // The aliases will always match before the templates. 329 // Note that you cannot put any "static if" around the aliases or 330 // compilation will fail due to conflict with the templates! 331 alias sqrt = std.math.sqrt; 332 alias sin = std.math.sin; 333 alias cos = std.math.cos; 334 alias exp = std.math.exp; 335 //alias fabs = std.math.fabs; 336 alias floor = std.math.floor; 337 alias exp2 = std.math.exp2; 338 alias ceil = std.math.ceil; 339 alias rint = std.math.rint; 340 341 /// 342 T sqrt(T)(in T x) if (isFloatingPoint!T) { return std.math.sqrt(x); } 343 /// 344 T sin(T)(in T x) if (isFloatingPoint!T) { return std.math.sin(x); } 345 /// 346 T cos(T)(in T x) if (isFloatingPoint!T) { return std.math.cos(x); } 347 /// 348 T pow(T)(in T x, in T power) if (isFloatingPoint!T) { alias y = power; mixin(mixinCMath2!`pow`); } 349 /// 350 T powi(T)(in T x, int power) if (isFloatingPoint!T) { alias y = power; mixin(mixinCMath2!`pow`); } 351 /// 352 T exp(T)(in T x) if (isFloatingPoint!T) { return std.math.exp(x); } 353 /// 354 T log(T)(in T x) if (isFloatingPoint!T) { mixin(mixinCMath!`log`); } 355 /// 356 T fabs(T)(in T x) if (isFloatingPoint!T) { return std.math.fabs(x); } 357 /// 358 T floor(T)(in T x) if (isFloatingPoint!T) { return std.math.floor(x); } 359 /// 360 T exp2(T)(in T x) if (isFloatingPoint!T) { return std.math.exp2(x); } 361 /// 362 T log10(T)(in T x) if (isFloatingPoint!T) { mixin(mixinCMath!`log10`); } 363 /// 364 T log2(T)(in T x) if (isFloatingPoint!T) { mixin(mixinCMath!`log2`); } 365 /// 366 T ceil(T)(in T x) if (isFloatingPoint!T) { return std.math.ceil(x); } 367 /// 368 T trunc(T)(in T x) if (isFloatingPoint!T) { mixin(mixinCMath!`trunc`); } 369 /// 370 T rint(T)(in T x) if (isFloatingPoint!T) { return std.math.rint(x); } 371 /// 372 T nearbyint(T)(in T x) if (isFloatingPoint!T) { mixin(mixinCMath!`nearbyint`); } 373 /// 374 T copysign(T)(in T mag, in T sgn) if (isFloatingPoint!T) 375 { 376 alias x = mag; 377 alias y = sgn; 378 mixin(mixinCMath2!`copysign`); 379 } 380 /// 381 T round(T)(in T x) if (isFloatingPoint!T) { mixin(mixinCMath!`round`); } 382 /// 383 T fmuladd(T)(in T a, in T b, in T c) if (isFloatingPoint!T) { return a * b + c; } 384 version(mir_core_test) 385 unittest { assert(fmuladd!double(2, 3, 4) == 2 * 3 + 4); } 386 /// 387 T fmin(T)(in T x, in T y) if (isFloatingPoint!T) 388 { 389 version (Windows) // https://issues.dlang.org/show_bug.cgi?id=19798 390 { 391 version (CRuntime_Microsoft) 392 mixin(mixinCMath2!`fmin`); 393 else 394 return std.math.fmin(x, y); 395 } 396 else 397 mixin(mixinCMath2!`fmin`); 398 } 399 /// 400 T fmax(T)(in T x, in T y) if (isFloatingPoint!T) 401 { 402 version (Windows) // https://issues.dlang.org/show_bug.cgi?id=19798 403 { 404 version (CRuntime_Microsoft) 405 mixin(mixinCMath2!`fmax`); 406 else 407 return std.math.fmax(x, y); 408 } 409 else 410 mixin(mixinCMath2!`fmax`); 411 } 412 413 version (mir_core_test) @nogc nothrow pure @safe unittest 414 { 415 // Check the aliases are correct. 416 static assert(is(typeof(sqrt(1.0f)) == float)); 417 static assert(is(typeof(sin(1.0f)) == float)); 418 static assert(is(typeof(cos(1.0f)) == float)); 419 static assert(is(typeof(exp(1.0f)) == float)); 420 static assert(is(typeof(fabs(1.0f)) == float)); 421 static assert(is(typeof(floor(1.0f)) == float)); 422 static assert(is(typeof(exp2(1.0f)) == float)); 423 static assert(is(typeof(ceil(1.0f)) == float)); 424 static assert(is(typeof(rint(1.0f)) == float)); 425 426 auto x = sqrt!float(2.0f); // Explicit template instantiation still works. 427 auto fp = &sqrt!float; // Can still take function address. 428 429 // Test for DMD linker problem with fmin on Windows. 430 static assert(is(typeof(fmin!float(1.0f, 1.0f)))); 431 static assert(is(typeof(fmax!float(1.0f, 1.0f)))); 432 } 433 } 434 else // DMD version prior to 2.082 435 { 436 static import std.math; 437 static import core.stdc.math; 438 439 // Calls either std.math or cmath function for either float (suffix "f") 440 // or double (no suffix). std.math will always be used during CTFE or for 441 // arguments with greater than double precision or if the cmath function 442 // is impure. 443 private enum mixinCMath(string fun) = 444 `pragma(inline, true); 445 static if (!is(typeof(std.math.`~fun~`(0.5f)) == float) 446 && is(typeof(() pure => core.stdc.math.`~fun~`f(0.5f)))) 447 if (!__ctfe) 448 { 449 static if (T.mant_dig == float.mant_dig) return core.stdc.math.`~fun~`f(x); 450 else static if (T.mant_dig == double.mant_dig) return core.stdc.math.`~fun~`(x); 451 } 452 return std.math.`~fun~`(x);`; 453 454 // As above but for two-argument function (both arguments must be floating point). 455 private enum mixinCMath2(string fun) = 456 `pragma(inline, true); 457 static if (!is(typeof(std.math.`~fun~`(0.5f, 0.5f)) == float) 458 && is(typeof(() pure => core.stdc.math.`~fun~`f(0.5f, 0.5f)))) 459 if (!__ctfe) 460 { 461 static if (T.mant_dig == float.mant_dig) return core.stdc.math.`~fun~`f(x, y); 462 else static if (T.mant_dig == double.mant_dig) return core.stdc.math.`~fun~`(x, y); 463 } 464 return std.math.`~fun~`(x, y);`; 465 466 // Some std.math functions have appropriate return types (float, 467 // double, real) without need for a wrapper. 468 alias sqrt = std.math.sqrt; 469 470 /// 471 T sqrt(T)(in T x) if (isFloatingPoint!T) { return std.math.sqrt(x); } 472 /// 473 T sin(T)(in T x) if (isFloatingPoint!T) { mixin(mixinCMath!`sin`); } 474 /// 475 T cos(T)(in T x) if (isFloatingPoint!T) { mixin(mixinCMath!`cos`); } 476 /// 477 T pow(T)(in T x, in T power) if (isFloatingPoint!T) { alias y = power; mixin(mixinCMath2!`pow`); } 478 /// 479 T powi(T)(in T x, int power) if (isFloatingPoint!T) { alias y = power; mixin(mixinCMath2!`pow`); } 480 /// 481 T exp(T)(in T x) if (isFloatingPoint!T) { mixin(mixinCMath!`exp`); } 482 /// 483 T log(T)(in T x) if (isFloatingPoint!T) { mixin(mixinCMath!`log`); } 484 /// 485 T fabs(T)(in T x) if (isFloatingPoint!T) { mixin(mixinCMath!`fabs`); } 486 /// 487 T floor(T)(in T x) if (isFloatingPoint!T) { mixin(mixinCMath!`floor`); } 488 /// 489 T exp2(T)(in T x) if (isFloatingPoint!T) { mixin(mixinCMath!`exp2`); } 490 /// 491 T log10(T)(in T x) if (isFloatingPoint!T) { mixin(mixinCMath!`log10`); } 492 /// 493 T log2(T)(in T x) if (isFloatingPoint!T) { mixin(mixinCMath!`log2`); } 494 /// 495 T ceil(T)(in T x) if (isFloatingPoint!T) { mixin(mixinCMath!`ceil`); } 496 /// 497 T trunc(T)(in T x) if (isFloatingPoint!T) { mixin(mixinCMath!`trunc`); } 498 /// 499 T rint(T)(in T x) if (isFloatingPoint!T) { mixin(mixinCMath!`rint`); } 500 /// 501 T nearbyint(T)(in T x) if (isFloatingPoint!T) { mixin(mixinCMath!`nearbyint`); } 502 /// 503 T copysign(T)(in T mag, in T sgn) if (isFloatingPoint!T) 504 { 505 alias x = mag; 506 alias y = sgn; 507 mixin(mixinCMath2!`copysign`); 508 } 509 /// 510 T round(T)(in T x) if (isFloatingPoint!T) { mixin(mixinCMath!`round`); } 511 /// 512 T fmuladd(T)(in T a, in T b, in T c) if (isFloatingPoint!T) { return a * b + c; } 513 version(mir_core_test) 514 unittest { assert(fmuladd!double(2, 3, 4) == 2 * 3 + 4); } 515 /// 516 T fmin(T)(in T x, in T y) if (isFloatingPoint!T) 517 { 518 version (Windows) // https://issues.dlang.org/show_bug.cgi?id=19798 519 { 520 version (CRuntime_Microsoft) 521 mixin(mixinCMath2!`fmin`); 522 else 523 return std.math.fmin(x, y); 524 } 525 else 526 mixin(mixinCMath2!`fmin`); 527 } 528 /// 529 T fmax(T)(in T x, in T y) if (isFloatingPoint!T) 530 { 531 version (Windows) // https://issues.dlang.org/show_bug.cgi?id=19798 532 { 533 version (CRuntime_Microsoft) 534 mixin(mixinCMath2!`fmax`); 535 else 536 return std.math.fmax(x, y); 537 } 538 else 539 mixin(mixinCMath2!`fmax`); 540 } 541 542 version (mir_core_test) @nogc nothrow pure @safe unittest 543 { 544 // Check the aliases are correct. 545 static assert(is(typeof(sqrt(1.0f)) == float)); 546 auto x = sqrt!float(2.0f); // Explicit template instantiation still works. 547 auto fp = &sqrt!float; // Can still take function address. 548 549 // Test for DMD linker problem with fmin on Windows. 550 static assert(is(typeof(fmin!float(1.0f, 1.0f)))); 551 static assert(is(typeof(fmax!float(1.0f, 1.0f)))); 552 } 553 } 554 555 version (mir_core_test) 556 @nogc nothrow pure @safe unittest 557 { 558 import mir.math: PI, feqrel; 559 assert(feqrel(pow(2.0L, -0.5L), cos(PI / 4)) >= real.mant_dig - 1); 560 } 561 562 /// Overload for cdouble, cfloat and creal 563 @optmath auto fabs(T)(in T x) 564 if (isComplex!T) 565 { 566 return x.re * x.re + x.im * x.im; 567 } 568 569 /// 570 version(mir_core_test) unittest 571 { 572 import mir.complex; 573 assert(fabs(Complex!double(3, 4)) == 25); 574 } 575 576 /++ 577 Computes whether two values are approximately equal, admitting a maximum 578 relative difference, and a maximum absolute difference. 579 Params: 580 lhs = First item to compare. 581 rhs = Second item to compare. 582 maxRelDiff = Maximum allowable difference relative to `rhs`. Defaults to `0.5 ^^ 20`. 583 maxAbsDiff = Maximum absolute difference. Defaults to `0.5 ^^ 20`. 584 585 Returns: 586 `true` if the two items are equal or approximately equal under either criterium. 587 +/ 588 bool approxEqual(T)(const T lhs, const T rhs, const T maxRelDiff = T(0x1p-20f), const T maxAbsDiff = T(0x1p-20f)) 589 if (isFloatingPoint!T) 590 { 591 if (rhs == lhs) // infs 592 return true; 593 auto diff = fabs(lhs - rhs); 594 if (diff <= maxAbsDiff) 595 return true; 596 diff /= fabs(rhs); 597 return diff <= maxRelDiff; 598 } 599 600 /// 601 @safe pure nothrow @nogc version(mir_core_test) unittest 602 { 603 assert(approxEqual(1.0, 1.0000001)); 604 assert(approxEqual(1.0f, 1.0000001f)); 605 assert(approxEqual(1.0L, 1.0000001L)); 606 607 assert(approxEqual(10000000.0, 10000001)); 608 assert(approxEqual(10000000f, 10000001f)); 609 assert(!approxEqual(100000.0L, 100001L)); 610 } 611 612 /// ditto 613 bool approxEqual(T)(const T lhs, const T rhs, float maxRelDiff = 0x1p-20f, float maxAbsDiff = 0x1p-20f) 614 if (isComplexOf!(T, float)) 615 { 616 return approxEqual(lhs.re, rhs.re, maxRelDiff, maxAbsDiff) 617 && approxEqual(lhs.im, rhs.im, maxRelDiff, maxAbsDiff); 618 } 619 620 /// ditto 621 bool approxEqual(T)(const T lhs, const T rhs, double maxRelDiff = 0x1p-20f, double maxAbsDiff = 0x1p-20f) 622 if (isComplexOf!(T, double)) 623 { 624 return approxEqual(lhs.re, rhs.re, maxRelDiff, maxAbsDiff) 625 && approxEqual(lhs.im, rhs.im, maxRelDiff, maxAbsDiff); 626 } 627 628 /// ditto 629 bool approxEqual(T)(const T lhs, const T rhs, real maxRelDiff = 0x1p-20f, real maxAbsDiff = 0x1p-20f) 630 if (isComplexOf!(T, real)) 631 { 632 return approxEqual(lhs.re, rhs.re, maxRelDiff, maxAbsDiff) 633 && approxEqual(lhs.im, rhs.im, maxRelDiff, maxAbsDiff); 634 } 635 636 /// Complex types works as `approxEqual(l.re, r.re) && approxEqual(l.im, r.im)` 637 @safe pure nothrow @nogc version(mir_core_test) unittest 638 { 639 import mir.internal.utility: isComplexOf; 640 static struct UserComplex(T) { T re, im; } 641 alias _cdouble = UserComplex!double; 642 643 static assert(isComplexOf!(_cdouble, double)); 644 645 assert(approxEqual(_cdouble(1.0, 1), _cdouble(1.0000001, 1), 1.0000001)); 646 assert(!approxEqual(_cdouble(100000.0L, 0), _cdouble(100001L, 0))); 647 }