1 /* This Source Code Form is subject to the terms of the Mozilla Public 2 * License, v. 2.0. If a copy of the MPL was not distributed with this 3 * file, You can obtain one at http://mozilla.org/MPL/2.0/. */ 4 5 /** 6 * This package provides mathematical functions. 7 * 8 * The $(D_PSYMBOL tanya.math) package itself provides only representation 9 * functions for built-in types, such as functions that provide information 10 * about internal representation of floating-point numbers and low-level 11 * operatons on these. Actual mathematical functions and additional types can 12 * be found in its submodules. $(D_PSYMBOL tanya.math) doesn't import any 13 * submodules publically, they should be imported explicitly. 14 * 15 * Copyright: Eugene Wissner 2016-2020. 16 * License: $(LINK2 https://www.mozilla.org/en-US/MPL/2.0/, 17 * Mozilla Public License, v. 2.0). 18 * Authors: $(LINK2 mailto:info@caraus.de, Eugene Wissner) 19 * Source: $(LINK2 https://github.com/caraus-ecms/tanya/blob/master/source/tanya/math/package.d, 20 * tanya/math/package.d) 21 */ 22 module tanya.math; 23 24 import tanya.math.nbtheory; 25 import tanya.meta.trait; 26 import tanya.meta.transform; 27 28 /// Floating-point number precisions according to IEEE-754. 29 enum IEEEPrecision : ubyte 30 { 31 single = 4, /// Single precision: 64-bit. 32 double_ = 8, /// Single precision: 64-bit. 33 doubleExtended = 10, /// Double extended precision: 80-bit. 34 } 35 36 /** 37 * Tests the precision of floating-point type $(D_PARAM F). 38 * 39 * For $(D_KEYWORD float) $(D_PSYMBOL ieeePrecision) always evaluates to 40 * $(D_INLINECODE IEEEPrecision.single); for $(D_KEYWORD double) - to 41 * $(D_INLINECODE IEEEPrecision.double). It returns different values only 42 * for $(D_KEYWORD real), since $(D_KEYWORD real) is a platform-dependent type. 43 * 44 * If $(D_PARAM F) is a $(D_KEYWORD real) and the target platform isn't 45 * currently supported, static assertion error will be raised (you can use 46 * $(D_INLINECODE is(typeof(ieeePrecision!F))) for testing the platform support 47 * without a compilation error). 48 * 49 * Params: 50 * F = Type to be tested. 51 * 52 * Returns: Precision according to IEEE-754. 53 * 54 * See_Also: $(D_PSYMBOL IEEEPrecision). 55 */ 56 template ieeePrecision(F) 57 if (isFloatingPoint!F) 58 { 59 static if (F.sizeof == float.sizeof) 60 { 61 enum IEEEPrecision ieeePrecision = IEEEPrecision.single; 62 } 63 else static if (F.sizeof == double.sizeof) 64 { 65 enum IEEEPrecision ieeePrecision = IEEEPrecision.double_; 66 } 67 else version (X86) 68 { 69 enum IEEEPrecision ieeePrecision = IEEEPrecision.doubleExtended; 70 } 71 else version (X86_64) 72 { 73 enum IEEEPrecision ieeePrecision = IEEEPrecision.doubleExtended; 74 } 75 else 76 { 77 static assert(false, "Unsupported IEEE 754 floating point precision"); 78 } 79 } 80 81 /// 82 @nogc nothrow pure @safe unittest 83 { 84 static assert(ieeePrecision!float == IEEEPrecision.single); 85 static assert(ieeePrecision!double == IEEEPrecision.double_); 86 } 87 88 package(tanya) union FloatBits(F) 89 { 90 Unqual!F floating; 91 static if (ieeePrecision!F == IEEEPrecision.single) 92 { 93 uint integral; 94 enum uint expMask = 0x7f800000; 95 } 96 else static if (ieeePrecision!F == IEEEPrecision.double_) 97 { 98 ulong integral; 99 enum ulong expMask = 0x7ff0000000000000; 100 } 101 else static if (ieeePrecision!F == IEEEPrecision.doubleExtended) 102 { 103 struct // Little-endian. 104 { 105 ulong mantissa; 106 ushort exp; 107 } 108 enum ulong mantissaMask = 0x7fffffffffffffff; 109 enum uint expMask = 0x7fff; 110 } 111 else 112 { 113 static assert(false, "Unsupported IEEE 754 floating point precision"); 114 } 115 } 116 117 /** 118 * Floating-point number classifications. 119 */ 120 enum FloatingPointClass : ubyte 121 { 122 /** 123 * Not a Number. 124 * 125 * See_Also: $(D_PSYMBOL isNaN). 126 */ 127 nan, 128 129 /// Zero. 130 zero, 131 132 /** 133 * Infinity. 134 * 135 * See_Also: $(D_PSYMBOL isInfinity). 136 */ 137 infinite, 138 139 /** 140 * Denormalized number. 141 * 142 * See_Also: $(D_PSYMBOL isSubnormal). 143 */ 144 subnormal, 145 146 /** 147 * Normalized number. 148 * 149 * See_Also: $(D_PSYMBOL isNormal). 150 */ 151 normal, 152 } 153 154 /** 155 * Returns whether $(D_PARAM x) is a NaN, zero, infinity, subnormal or 156 * normalized number. 157 * 158 * This function doesn't distinguish between negative and positive infinity, 159 * negative and positive NaN or negative and positive zero. 160 * 161 * Params: 162 * F = Type of the floating point number. 163 * x = Floating point number. 164 * 165 * Returns: Classification of $(D_PARAM x). 166 */ 167 FloatingPointClass classify(F)(F x) 168 if (isFloatingPoint!F) 169 { 170 if (x == 0) 171 { 172 return FloatingPointClass.zero; 173 } 174 FloatBits!F bits; 175 bits.floating = abs(x); 176 177 static if (ieeePrecision!F == IEEEPrecision.single) 178 { 179 if (bits.integral > bits.expMask) 180 { 181 return FloatingPointClass.nan; 182 } 183 else if (bits.integral == bits.expMask) 184 { 185 return FloatingPointClass.infinite; 186 } 187 else if (bits.integral < (1 << 23)) 188 { 189 return FloatingPointClass.subnormal; 190 } 191 } 192 else static if (ieeePrecision!F == IEEEPrecision.double_) 193 { 194 if (bits.integral > bits.expMask) 195 { 196 return FloatingPointClass.nan; 197 } 198 else if (bits.integral == bits.expMask) 199 { 200 return FloatingPointClass.infinite; 201 } 202 else if (bits.integral < (1L << 52)) 203 { 204 return FloatingPointClass.subnormal; 205 } 206 } 207 else static if (ieeePrecision!F == IEEEPrecision.doubleExtended) 208 { 209 if (bits.exp == bits.expMask) 210 { 211 if ((bits.mantissa & bits.mantissaMask) == 0) 212 { 213 return FloatingPointClass.infinite; 214 } 215 else 216 { 217 return FloatingPointClass.nan; 218 } 219 } 220 else if (bits.exp == 0) 221 { 222 return FloatingPointClass.subnormal; 223 } 224 else if (bits.mantissa < (1L << 63)) // "Unnormal". 225 { 226 return FloatingPointClass.nan; 227 } 228 } 229 230 return FloatingPointClass.normal; 231 } 232 233 /// 234 @nogc nothrow pure @safe unittest 235 { 236 assert(classify(0.0) == FloatingPointClass.zero); 237 assert(classify(double.nan) == FloatingPointClass.nan); 238 assert(classify(double.infinity) == FloatingPointClass.infinite); 239 assert(classify(-double.infinity) == FloatingPointClass.infinite); 240 assert(classify(1.4) == FloatingPointClass.normal); 241 assert(classify(1.11254e-307 / 10) == FloatingPointClass.subnormal); 242 243 assert(classify(0.0f) == FloatingPointClass.zero); 244 assert(classify(float.nan) == FloatingPointClass.nan); 245 assert(classify(float.infinity) == FloatingPointClass.infinite); 246 assert(classify(-float.infinity) == FloatingPointClass.infinite); 247 assert(classify(0.3) == FloatingPointClass.normal); 248 assert(classify(5.87747e-38f / 10) == FloatingPointClass.subnormal); 249 250 assert(classify(0.0L) == FloatingPointClass.zero); 251 assert(classify(real.nan) == FloatingPointClass.nan); 252 assert(classify(real.infinity) == FloatingPointClass.infinite); 253 assert(classify(-real.infinity) == FloatingPointClass.infinite); 254 } 255 256 /** 257 * Determines whether $(D_PARAM x) is a finite number. 258 * 259 * Params: 260 * F = Type of the floating point number. 261 * x = Floating point number. 262 * 263 * Returns: $(D_KEYWORD true) if $(D_PARAM x) is a finite number, 264 * $(D_KEYWORD false) otherwise. 265 * 266 * See_Also: $(D_PSYMBOL isInfinity). 267 */ 268 bool isFinite(F)(F x) 269 if (isFloatingPoint!F) 270 { 271 FloatBits!F bits; 272 static if (ieeePrecision!F == IEEEPrecision.single 273 || ieeePrecision!F == IEEEPrecision.double_) 274 { 275 bits.floating = x; 276 bits.integral &= bits.expMask; 277 return bits.integral != bits.expMask; 278 } 279 else static if (ieeePrecision!F == IEEEPrecision.doubleExtended) 280 { 281 bits.floating = abs(x); 282 return (bits.exp != bits.expMask) 283 && (bits.exp == 0 || bits.mantissa >= (1L << 63)); 284 } 285 } 286 287 /// 288 @nogc nothrow pure @safe unittest 289 { 290 assert(!isFinite(float.infinity)); 291 assert(!isFinite(-double.infinity)); 292 assert(isFinite(0.0)); 293 assert(!isFinite(float.nan)); 294 assert(isFinite(5.87747e-38f / 10)); 295 assert(isFinite(1.11254e-307 / 10)); 296 assert(isFinite(0.5)); 297 } 298 299 /** 300 * Determines whether $(D_PARAM x) is $(B n)ot $(B a) $(B n)umber (NaN). 301 * 302 * Params: 303 * F = Type of the floating point number. 304 * x = Floating point number. 305 * 306 * Returns: $(D_KEYWORD true) if $(D_PARAM x) is not a number, 307 * $(D_KEYWORD false) otherwise. 308 */ 309 bool isNaN(F)(F x) 310 if (isFloatingPoint!F) 311 { 312 FloatBits!F bits; 313 bits.floating = abs(x); 314 315 static if (ieeePrecision!F == IEEEPrecision.single 316 || ieeePrecision!F == IEEEPrecision.double_) 317 { 318 return bits.integral > bits.expMask; 319 } 320 else static if (ieeePrecision!F == IEEEPrecision.doubleExtended) 321 { 322 const maskedMantissa = bits.mantissa & bits.mantissaMask; 323 if ((bits.exp == bits.expMask && maskedMantissa != 0) 324 || ((bits.exp != 0) && (bits.mantissa < (1L << 63)))) 325 { 326 return true; 327 } 328 return false; 329 } 330 } 331 332 /// 333 @nogc nothrow pure @safe unittest 334 { 335 assert(isNaN(float.init)); 336 assert(isNaN(double.init)); 337 assert(isNaN(real.init)); 338 } 339 340 /** 341 * Determines whether $(D_PARAM x) is a positive or negative infinity. 342 * 343 * Params: 344 * F = Type of the floating point number. 345 * x = Floating point number. 346 * 347 * Returns: $(D_KEYWORD true) if $(D_PARAM x) is infinity, $(D_KEYWORD false) 348 * otherwise. 349 * 350 * See_Also: $(D_PSYMBOL isFinite). 351 */ 352 bool isInfinity(F)(F x) 353 if (isFloatingPoint!F) 354 { 355 FloatBits!F bits; 356 bits.floating = abs(x); 357 static if (ieeePrecision!F == IEEEPrecision.single 358 || ieeePrecision!F == IEEEPrecision.double_) 359 { 360 return bits.integral == bits.expMask; 361 } 362 else static if (ieeePrecision!F == IEEEPrecision.doubleExtended) 363 { 364 return (bits.exp == bits.expMask) 365 && ((bits.mantissa & bits.mantissaMask) == 0); 366 } 367 } 368 369 /// 370 @nogc nothrow pure @safe unittest 371 { 372 assert(isInfinity(float.infinity)); 373 assert(isInfinity(-float.infinity)); 374 assert(isInfinity(double.infinity)); 375 assert(isInfinity(-double.infinity)); 376 assert(isInfinity(real.infinity)); 377 assert(isInfinity(-real.infinity)); 378 } 379 380 /** 381 * Determines whether $(D_PARAM x) is a denormilized number or not. 382 * 383 * Denormalized number is a number between `0` and `1` that cannot be 384 * represented as 385 * 386 * <pre> 387 * m*2<sup>e</sup> 388 * </pre> 389 * 390 * where $(I m) is the mantissa and $(I e) is an exponent that fits into the 391 * exponent field of the type $(D_PARAM F). 392 * 393 * `0` is neither normalized nor denormalized. 394 * 395 * Params: 396 * F = Type of the floating point number. 397 * x = Floating point number. 398 * 399 * Returns: $(D_KEYWORD true) if $(D_PARAM x) is a denormilized number, 400 * $(D_KEYWORD false) otherwise. 401 * 402 * See_Also: $(D_PSYMBOL isNormal). 403 */ 404 bool isSubnormal(F)(F x) 405 if (isFloatingPoint!F) 406 { 407 FloatBits!F bits; 408 bits.floating = abs(x); 409 static if (ieeePrecision!F == IEEEPrecision.single) 410 { 411 return bits.integral < (1 << 23) && bits.integral > 0; 412 } 413 else static if (ieeePrecision!F == IEEEPrecision.double_) 414 { 415 return bits.integral < (1L << 52) && bits.integral > 0; 416 } 417 else static if (ieeePrecision!F == IEEEPrecision.doubleExtended) 418 { 419 return bits.exp == 0 && bits.mantissa != 0; 420 } 421 } 422 423 /// 424 @nogc nothrow pure @safe unittest 425 { 426 assert(!isSubnormal(0.0f)); 427 assert(!isSubnormal(float.nan)); 428 assert(!isSubnormal(float.infinity)); 429 assert(!isSubnormal(0.3f)); 430 assert(isSubnormal(5.87747e-38f / 10)); 431 432 assert(!isSubnormal(0.0)); 433 assert(!isSubnormal(double.nan)); 434 assert(!isSubnormal(double.infinity)); 435 assert(!isSubnormal(1.4)); 436 assert(isSubnormal(1.11254e-307 / 10)); 437 438 assert(!isSubnormal(0.0L)); 439 assert(!isSubnormal(real.nan)); 440 assert(!isSubnormal(real.infinity)); 441 } 442 443 /** 444 * Determines whether $(D_PARAM x) is a normilized number or not. 445 * 446 * Normalized number is a number that can be represented as 447 * 448 * <pre> 449 * m*2<sup>e</sup> 450 * </pre> 451 * 452 * where $(I m) is the mantissa and $(I e) is an exponent that fits into the 453 * exponent field of the type $(D_PARAM F). 454 * 455 * `0` is neither normalized nor denormalized. 456 * 457 * Params: 458 * F = Type of the floating point number. 459 * x = Floating point number. 460 * 461 * Returns: $(D_KEYWORD true) if $(D_PARAM x) is a normilized number, 462 * $(D_KEYWORD false) otherwise. 463 * 464 * See_Also: $(D_PSYMBOL isSubnormal). 465 */ 466 bool isNormal(F)(F x) 467 if (isFloatingPoint!F) 468 { 469 static if (ieeePrecision!F == IEEEPrecision.single 470 || ieeePrecision!F == IEEEPrecision.double_) 471 { 472 FloatBits!F bits; 473 bits.floating = x; 474 bits.integral &= bits.expMask; 475 return bits.integral != 0 && bits.integral != bits.expMask; 476 } 477 else static if (ieeePrecision!F == IEEEPrecision.doubleExtended) 478 { 479 return classify(x) == FloatingPointClass.normal; 480 } 481 } 482 483 /// 484 @nogc nothrow pure @safe unittest 485 { 486 assert(!isNormal(0.0f)); 487 assert(!isNormal(float.nan)); 488 assert(!isNormal(float.infinity)); 489 assert(isNormal(0.3f)); 490 assert(!isNormal(5.87747e-38f / 10)); 491 492 assert(!isNormal(0.0)); 493 assert(!isNormal(double.nan)); 494 assert(!isNormal(double.infinity)); 495 assert(isNormal(1.4)); 496 assert(!isNormal(1.11254e-307 / 10)); 497 498 assert(!isNormal(0.0L)); 499 assert(!isNormal(real.nan)); 500 assert(!isNormal(real.infinity)); 501 } 502 503 /** 504 * Determines whether the sign bit of $(D_PARAM x) is set or not. 505 * 506 * If the sign bit, $(D_PARAM x) is a negative number, otherwise positive. 507 * 508 * Params: 509 * F = Type of the floating point number. 510 * x = Floating point number. 511 * 512 * Returns: $(D_KEYWORD true) if the sign bit of $(D_PARAM x) is set, 513 * $(D_KEYWORD false) otherwise. 514 */ 515 bool signBit(F)(F x) 516 if (isFloatingPoint!F) 517 { 518 FloatBits!F bits; 519 bits.floating = x; 520 static if (ieeePrecision!F == IEEEPrecision.single) 521 { 522 return (bits.integral & (1 << 31)) != 0; 523 } 524 else static if (ieeePrecision!F == IEEEPrecision.double_) 525 { 526 return (bits.integral & (1L << 63)) != 0; 527 } 528 else static if (ieeePrecision!F == IEEEPrecision.doubleExtended) 529 { 530 return (bits.exp & (1 << 15)) != 0; 531 } 532 } 533 534 /// 535 @nogc nothrow pure @safe unittest 536 { 537 assert(signBit(-1.0f)); 538 assert(!signBit(1.0f)); 539 540 assert(signBit(-1.0)); 541 assert(!signBit(1.0)); 542 543 assert(signBit(-1.0L)); 544 assert(!signBit(1.0L)); 545 } 546 547 /** 548 * Computes $(D_PARAM x) to the power $(D_PARAM y) modulo $(D_PARAM z). 549 * 550 * Params: 551 * I = Base type. 552 * G = Exponent type. 553 * H = Divisor type: 554 * x = Base. 555 * y = Exponent. 556 * z = Divisor. 557 * 558 * Returns: Reminder of the division of $(D_PARAM x) to the power $(D_PARAM y) 559 * by $(D_PARAM z). 560 * 561 * Precondition: $(D_INLINECODE z > 0) 562 */ 563 H pow(I, G, H)(in auto ref I x, in auto ref G y, in auto ref H z) 564 if (isIntegral!I && isIntegral!G && isIntegral!H) 565 in 566 { 567 assert(z > 0, "Division by zero"); 568 } 569 do 570 { 571 G mask = G.max / 2 + 1; 572 H result; 573 574 if (y == 0) 575 { 576 return 1 % z; 577 } 578 else if (y == 1) 579 { 580 return x % z; 581 } 582 do 583 { 584 immutable bit = y & mask; 585 if (!result && bit) 586 { 587 result = x; 588 continue; 589 } 590 591 result *= result; 592 if (bit) 593 { 594 result *= x; 595 } 596 result %= z; 597 } 598 while (mask >>= 1); 599 600 return result; 601 } 602 603 /// 604 @nogc nothrow pure @safe unittest 605 { 606 assert(pow(3, 5, 7) == 5); 607 assert(pow(2, 2, 1) == 0); 608 assert(pow(3, 3, 3) == 0); 609 assert(pow(7, 4, 2) == 1); 610 assert(pow(53, 0, 2) == 1); 611 assert(pow(53, 1, 3) == 2); 612 assert(pow(53, 2, 5) == 4); 613 assert(pow(0, 0, 5) == 1); 614 assert(pow(0, 5, 5) == 0); 615 }