1 /++ 2 Note: 3 The module doesn't provide full arithmetic API for now. 4 +/ 5 module mir.bignum.fp; 6 7 import mir.bitop; 8 import mir.utility; 9 import std.traits; 10 11 package enum half(uint hs) = (){ 12 import mir.bignum.fixed: UInt; 13 UInt!hs ret; ret.signBit = true; return ret; 14 }(); 15 16 /++ 17 Software floating point number. 18 19 Params: 20 size = coefficient size in bits 21 +/ 22 struct Fp(uint size) 23 if (size % (uint.sizeof * 8) == 0 && size >= (uint.sizeof * 8)) 24 { 25 import mir.bignum.fixed: UInt; 26 27 bool sign; 28 long exponent; 29 UInt!size coefficient; 30 31 /++ 32 +/ 33 nothrow 34 this(bool sign, long exponent, UInt!size normalizedCoefficient) 35 { 36 this.coefficient = normalizedCoefficient; 37 this.exponent = exponent; 38 this.sign = sign; 39 } 40 41 /++ 42 Constructs $(LREF Fp) from hardaware floating point number. 43 Params: 44 value = Hardware floating point number. Special values `nan` and `inf` aren't allowed. 45 normalize = flag to indicate if the normalization should be performed. 46 +/ 47 this(T)(const T value, bool normalize = true) 48 @safe pure nothrow @nogc 49 if (isFloatingPoint!T && T.mant_dig <= size) 50 { 51 import mir.math.common : fabs; 52 import mir.math.ieee : frexp, signbit, ldexp; 53 this.sign = value.signbit != 0; 54 if (value == 0) 55 return; 56 T x = value.fabs; 57 if (_expect(!(x < T.infinity), false)) 58 { 59 this.exponent = this.exponent.max; 60 this.coefficient = x != T.infinity; 61 return; 62 } 63 int exp; 64 { 65 enum scale = T(2) ^^ T.mant_dig; 66 x = frexp(x, exp) * scale; 67 } 68 69 static if (T.mant_dig < 64) 70 { 71 auto xx = cast(ulong)cast(long)x; 72 if (normalize) 73 { 74 auto shift = ctlz(xx); 75 exp -= shift + T.mant_dig + size - 64; 76 xx <<= shift; 77 this.coefficient = UInt!64(xx).rightExtend!(size - 64); 78 } 79 else 80 { 81 this.coefficient = xx; 82 } 83 } 84 else 85 static if (T.mant_dig == 64) 86 { 87 auto xx = cast(ulong)x; 88 if (normalize) 89 { 90 auto shift = ctlz(xx); 91 exp -= shift + T.mant_dig + size - 64; 92 xx <<= shift; 93 this.coefficient = UInt!64(xx).rightExtend!(size - 64); 94 } 95 else 96 { 97 this.coefficient = xx; 98 } 99 } 100 else 101 { 102 enum scale = T(2) ^^ 64; 103 enum scaleInv = 1 / scale; 104 x *= scaleInv; 105 long high = cast(long) x; 106 if (high > x) 107 --high; 108 x -= high; 109 x *= scale; 110 auto most = ulong(high); 111 auto least = cast(ulong)x; 112 version(LittleEndian) 113 ulong[2] pair = [least, most]; 114 else 115 ulong[2] pair = [most, least]; 116 117 if (normalize) 118 { 119 this.coefficient = UInt!128(pair).rightExtend!(size - 128); 120 auto shift = most ? ctlz(most) : ctlz(least) + 64; 121 exp -= shift + T.mant_dig + size - 64 * (1 + (T.mant_dig > 64)); 122 this.coefficient <<= shift; 123 } 124 else 125 { 126 this.coefficient = pair; 127 } 128 } 129 if (!normalize) 130 { 131 exp -= T.mant_dig; 132 int shift = T.min_exp - T.mant_dig - exp; 133 if (shift > 0) 134 { 135 this.coefficient >>= shift; 136 exp = T.min_exp - T.mant_dig; 137 } 138 } 139 this.exponent = exp; 140 } 141 142 static if (size == 128) 143 /// 144 version(mir_bignum_test) 145 @safe pure @nogc nothrow 146 unittest 147 { 148 enum h = -33.0 * 2.0 ^^ -10; 149 auto f = Fp!64(h); 150 assert(f.sign); 151 assert(f.exponent == -10 - (64 - 6)); 152 assert(f.coefficient == 33UL << (64 - 6)); 153 assert(cast(double) f == h); 154 155 // CTFE 156 static assert(cast(double) Fp!64(h) == h); 157 158 f = Fp!64(-0.0); 159 assert(f.sign); 160 assert(f.exponent == 0); 161 assert(f.coefficient == 0); 162 163 // subnormals 164 static assert(cast(float) Fp!64(float.min_normal / 2) == float.min_normal / 2); 165 static assert(cast(float) Fp!64(float.min_normal * float.epsilon) == float.min_normal * float.epsilon); 166 // subnormals 167 static assert(cast(double) Fp!64(double.min_normal / 2) == double.min_normal / 2); 168 static assert(cast(double) Fp!64(double.min_normal * double.epsilon) == double.min_normal * double.epsilon); 169 // subnormals 170 static if (real.mant_dig <= 64) 171 { 172 static assert(cast(real) Fp!128(real.min_normal / 2) == real.min_normal / 2); 173 static assert(cast(real) Fp!128(real.min_normal * real.epsilon) == real.min_normal * real.epsilon); 174 } 175 176 enum d = cast(float) Fp!64(float.min_normal / 2, false); 177 178 // subnormals 179 static assert(cast(float) Fp!64(float.min_normal / 2, false) == float.min_normal / 2, d.stringof); 180 static assert(cast(float) Fp!64(float.min_normal * float.epsilon, false) == float.min_normal * float.epsilon); 181 // subnormals 182 static assert(cast(double) Fp!64(double.min_normal / 2, false) == double.min_normal / 2); 183 static assert(cast(double) Fp!64(double.min_normal * double.epsilon, false) == double.min_normal * double.epsilon); 184 // subnormals 185 static if (real.mant_dig <= 64) 186 { 187 static assert(cast(real) Fp!64(real.min_normal / 2, false) == real.min_normal / 2); 188 static assert(cast(real) Fp!64(real.min_normal * real.epsilon, false) == real.min_normal * real.epsilon); 189 } 190 191 import mir.bignum.fixed: UInt; 192 193 assert(cast(double)Fp!128(+double.infinity) == +double.infinity); 194 assert(cast(double)Fp!128(-double.infinity) == -double.infinity); 195 196 import mir.math.ieee : signbit; 197 auto r = cast(double)Fp!128(-double.nan); 198 assert(r != r && r.signbit); 199 } 200 201 // static if (size == 128) 202 // /// Without normalization 203 // version(mir_bignum_test) 204 // @safe pure @nogc nothrow 205 // unittest 206 // { 207 // auto f = Fp!64(-33.0 * 2.0 ^^ -10, false); 208 // assert(f.sign); 209 // assert(f.exponent == -10 - (double.mant_dig - 6)); 210 // assert(f.coefficient == 33UL << (double.mant_dig - 6)); 211 // } 212 213 /++ 214 +/ 215 this(uint isize)(UInt!isize integer, bool normalizedInteger = false) 216 nothrow 217 { 218 import mir.bignum.fixed: UInt; 219 static if (isize < size) 220 { 221 if (normalizedInteger) 222 { 223 this(false, long(isize) - size, integer.rightExtend!(size - isize)); 224 } 225 else 226 { 227 this(integer.toSize!size, false); 228 } 229 } 230 else 231 { 232 if (!integer) 233 return; 234 this.exponent = isize - size; 235 if (!normalizedInteger) 236 { 237 auto c = integer.ctlz; 238 integer <<= c; 239 this.exponent -= c; 240 } 241 static if (isize == size) 242 { 243 coefficient = integer; 244 } 245 else 246 { 247 enum N = coefficient.data.length; 248 version (LittleEndian) 249 coefficient.data = integer.data[$ - N .. $]; 250 else 251 coefficient.data = integer.data[0 .. N]; 252 enum tailSize = isize - size; 253 auto cr = integer.toSize!tailSize.opCmp(half!tailSize); 254 if (cr > 0 || cr == 0 && coefficient.bt(0)) 255 { 256 if (auto overflow = coefficient += 1) 257 { 258 coefficient = half!size; 259 exponent++; 260 } 261 } 262 } 263 } 264 } 265 266 static if (size == 128) 267 /// 268 version(mir_bignum_test) 269 @safe pure @nogc 270 unittest 271 { 272 import mir.bignum.fixed: UInt; 273 274 auto fp = Fp!128(UInt!128.fromHexString("afbbfae3cd0aff2714a1de7022b0029d")); 275 assert(fp.exponent == 0); 276 assert(fp.coefficient == UInt!128.fromHexString("afbbfae3cd0aff2714a1de7022b0029d")); 277 278 fp = Fp!128(UInt!128.fromHexString("afbbfae3cd0aff2714a1de7022b0029d"), true); 279 assert(fp.exponent == 0); 280 assert(fp.coefficient == UInt!128.fromHexString("afbbfae3cd0aff2714a1de7022b0029d")); 281 282 fp = Fp!128(UInt!128.fromHexString("ae3cd0aff2714a1de7022b0029d")); 283 assert(fp.exponent == -20); 284 assert(fp.coefficient == UInt!128.fromHexString("ae3cd0aff2714a1de7022b0029d00000")); 285 286 fp = Fp!128(UInt!128.fromHexString("e7022b0029d")); 287 assert(fp.exponent == -84); 288 assert(fp.coefficient == UInt!128.fromHexString("e7022b0029d000000000000000000000")); 289 290 fp = Fp!128(UInt!64.fromHexString("e7022b0029d")); 291 assert(fp.exponent == -84); 292 assert(fp.coefficient == UInt!128.fromHexString("e7022b0029d000000000000000000000")); 293 294 fp = Fp!128(UInt!64.fromHexString("e7022b0029dd0aff"), true); 295 assert(fp.exponent == -64); 296 assert(fp.coefficient == UInt!128.fromHexString("e7022b0029dd0aff0000000000000000")); 297 298 fp = Fp!128(UInt!64.fromHexString("e7022b0029d")); 299 assert(fp.exponent == -84); 300 assert(fp.coefficient == UInt!128.fromHexString("e7022b0029d000000000000000000000")); 301 302 fp = Fp!128(UInt!192.fromHexString("ffffffffffffffffffffffffffffffff1000000000000000")); 303 assert(fp.exponent == 64); 304 assert(fp.coefficient == UInt!128.fromHexString("ffffffffffffffffffffffffffffffff")); 305 306 fp = Fp!128(UInt!192.fromHexString("ffffffffffffffffffffffffffffffff8000000000000000")); 307 assert(fp.exponent == 65); 308 assert(fp.coefficient == UInt!128.fromHexString("80000000000000000000000000000000")); 309 310 fp = Fp!128(UInt!192.fromHexString("fffffffffffffffffffffffffffffffe8000000000000000")); 311 assert(fp.exponent == 64); 312 assert(fp.coefficient == UInt!128.fromHexString("fffffffffffffffffffffffffffffffe")); 313 314 fp = Fp!128(UInt!192.fromHexString("fffffffffffffffffffffffffffffffe8000000000000001")); 315 assert(fp.exponent == 64); 316 assert(fp.coefficient == UInt!128.fromHexString("ffffffffffffffffffffffffffffffff")); 317 } 318 319 /// ditto 320 this(ulong value) 321 { 322 this(UInt!64(value)); 323 } 324 325 /// 326 this(long value) 327 { 328 this(ulong(value >= 0 ? value : -value)); 329 this.sign = !(value >= 0); 330 } 331 332 /// 333 this(int value) 334 { 335 this(long(value)); 336 } 337 338 /// 339 this(uint value) 340 { 341 this(ulong(value)); 342 } 343 344 /// 345 bool isNaN() scope const @property 346 { 347 return this.exponent == this.exponent.max && this.coefficient != this.coefficient.init; 348 } 349 350 /// 351 bool isInfinity() scope const @property 352 { 353 return this.exponent == this.exponent.max && this.coefficient == coefficient.init; 354 } 355 356 /// 357 bool isSpecial() scope const @property 358 { 359 return this.exponent == this.exponent.max; 360 } 361 362 /// 363 bool opEquals(const Fp rhs) scope const 364 { 365 if (this.exponent != rhs.exponent) 366 return false; 367 if (this.coefficient != rhs.coefficient) 368 return false; 369 if (this.coefficient == 0) 370 return !this.isSpecial || this.sign == rhs.sign; 371 if (this.sign != rhs.sign) 372 return false; 373 return !this.isSpecial; 374 } 375 376 /// 377 ref Fp opOpAssign(string op)(Fp rhs) nothrow scope return 378 if (op == "*" || op == "/") 379 { 380 this = this.opBinary!op(rhs); 381 return this; 382 } 383 384 /// 385 Fp!(max(size, rhsSize)) opBinary(string op : "*", uint rhsSize)(Fp!rhsSize rhs) nothrow const 386 { 387 return cast(Fp) .extendedMul(cast()this, rhs); 388 } 389 390 static if (size == 128) 391 /// 392 version(mir_bignum_test) 393 @safe pure @nogc 394 unittest 395 { 396 import mir.bignum.fixed: UInt; 397 398 auto a = Fp!128(0, -13, UInt!128.fromHexString("dfbbfae3cd0aff2714a1de7022b0029d")); 399 auto b = Fp!128(1, 100, UInt!128.fromHexString("e3251bacb112c88b71ad3f85a970a314")); 400 auto fp = a.opBinary!"*"(b); 401 assert(fp.sign); 402 assert(fp.exponent == 100 - 13 + 128); 403 assert(fp.coefficient == UInt!128.fromHexString("c6841dd302415d785373ab6d93712988")); 404 } 405 406 /// Uses approximate division for now 407 /// TODO: use full precision division for void when Fp division is ready 408 Fp!(max(size, rhsSize)) opBinary(string op : "/", uint rhsSize)(Fp!rhsSize rhs) nothrow const 409 { 410 Fp a = this; 411 alias b = rhs; 412 auto exponent = a.exponent - b.exponent; 413 a.exponent = b.exponent = -long(size); 414 auto ret = typeof(return)(cast(real) a / cast(real) b); 415 ret.exponent += exponent; 416 return ret; 417 } 418 419 /// 420 T opCast(T)() nothrow const 421 if (is(Unqual!T == bool)) 422 { 423 return exponent || coefficient; 424 } 425 426 /// 427 T opCast(T, bool noSpecial = false, bool noHalf = false)() nothrow const 428 if (is(T == float) || is(T == double) || is(T == real)) 429 { 430 import mir.math.ieee: ldexp; 431 static if (!noSpecial) 432 { 433 if (_expect(this.isSpecial, false)) 434 { 435 T ret = this.coefficient ? T.nan : T.infinity; 436 if (this.sign) 437 ret = -ret; 438 return ret; 439 } 440 } 441 auto exp = cast()this.exponent; 442 static if (size == 32) 443 { 444 T c = cast(uint) coefficient; 445 } 446 else 447 static if (size == 64) 448 { 449 T c = cast(ulong) coefficient; 450 } 451 else 452 { 453 enum shift = size - T.mant_dig; 454 enum rMask = (UInt!size(1) << shift) - UInt!size(1); 455 enum rHalf = UInt!size(1) << (shift - 1); 456 enum rInc = UInt!size(1) << shift; 457 UInt!size adC = this.coefficient; 458 static if (!noHalf) 459 { 460 auto cr = (this.coefficient & rMask).opCmp(rHalf); 461 if ((cr > 0) | (cr == 0) & this.coefficient.bt(shift)) 462 { 463 if (auto overflow = adC += rInc) 464 { 465 adC = half!size; 466 exp++; 467 } 468 } 469 } 470 adC >>= shift; 471 exp += shift; 472 T c = cast(ulong) adC; 473 static if (T.mant_dig > 64) // 474 { 475 static assert (T.mant_dig <= 128); 476 c += ldexp(cast(T) cast(ulong) (adC >> 64), 64); 477 } 478 } 479 if (this.sign) 480 c = -c; 481 static if (exp.sizeof > int.sizeof) 482 { 483 import mir.utility: min, max; 484 exp = exp.max(int.min).min(int.max); 485 } 486 return ldexp(c, cast(int)exp); 487 } 488 489 static if (size == 128) 490 /// 491 version(mir_bignum_test) 492 @safe pure @nogc 493 unittest 494 { 495 import mir.bignum.fixed: UInt; 496 auto fp = Fp!128(1, 100, UInt!128.fromHexString("e3251bacb112cb8b71ad3f85a970a314")); 497 assert(cast(double)fp == -0xE3251BACB112C8p+172); 498 499 fp = Fp!128(1, long.max, UInt!128.init); 500 assert(cast(double)fp == -double.infinity); 501 502 import mir.math.ieee : signbit; 503 fp = Fp!128(1, long.max, UInt!128(123)); 504 auto r = cast(double)fp; 505 assert(r != r && r.signbit); 506 } 507 508 static if (size == 128) 509 /// 510 version(mir_bignum_test) 511 @safe pure @nogc 512 unittest 513 { 514 import mir.bignum.fixed: UInt; 515 auto fp = Fp!128(1, 100, UInt!128.fromHexString("e3251bacb112cb8b71ad3f85a970a314")); 516 static if (real.mant_dig == 64) 517 assert(cast(real)fp == -0xe3251bacb112cb8bp+164L); 518 } 519 520 static if (size == 128) 521 /// 522 version(mir_bignum_test) 523 @safe pure @nogc 524 unittest 525 { 526 import mir.bignum.fixed: UInt; 527 auto fp = Fp!64(1, 100, UInt!64(0xe3251bacb112cb8b)); 528 version (DigitalMars) 529 { 530 // https://issues.dlang.org/show_bug.cgi?id=20963 531 assert(cast(double)fp == -0xE3251BACB112C8p+108 532 || cast(double)fp == -0xE3251BACB112D0p+108); 533 } 534 else 535 { 536 assert(cast(double)fp == -0xE3251BACB112C8p+108); 537 } 538 } 539 // -0x1.c64a375962259p+163 = 540 // -0xe.3251bacb112cb8bp+160 = 541 // -0x1.c64a37596225ap+163 = 542 // -0xe.3251bacb112cb8bp+160 = 543 static if (size == 128) 544 /// 545 version(mir_bignum_test) 546 @safe pure @nogc 547 unittest 548 { 549 import mir.bignum.fixed: UInt; 550 auto fp = Fp!64(1, 100, UInt!64(0xe3251bacb112cb8b)); 551 static if (real.mant_dig == 64) 552 assert(cast(real)fp == -0xe3251bacb112cb8bp+100L); 553 } 554 555 /// 556 T opCast(T : Fp!newSize, bool noSpecial = false, size_t newSize)() nothrow const 557 if (newSize != size) 558 { 559 Fp!newSize ret; 560 ret.sign = this.sign; 561 562 static if (!noSpecial) 563 { 564 if (_expect(this.isSpecial, false)) 565 { 566 ret.exponent = ret.exponent.max; 567 ret.coefficient = !!this.coefficient; 568 return ret; 569 } 570 if (!this) 571 { 572 return ret; 573 } 574 } 575 576 UInt!size coefficient = this.coefficient; 577 int shift; 578 // subnormal 579 580 static if (!noSpecial) 581 { 582 if (this.exponent == this.exponent.min) 583 { 584 shift = cast(int)coefficient.ctlz; 585 coefficient <<= shift; 586 } 587 } 588 589 ret = Fp!newSize(coefficient, true); 590 ret.exponent -= shift; 591 ret.sign = this.sign; 592 593 import mir.checkedint: adds; 594 /// overflow 595 596 static if (!noSpecial) 597 { 598 bool overflow; 599 ret.exponent = adds(ret.exponent, this.exponent, overflow); 600 if (_expect(overflow, false)) 601 { 602 // overflow 603 if (this.exponent > 0) 604 { 605 ret.exponent = ret.exponent.max; 606 ret.coefficient = 0u; 607 } 608 // underflow 609 else 610 { 611 ret.coefficient >>= cast(uint)(ret.exponent - exponent.min); 612 ret.exponent = ret.coefficient ? ret.exponent.min : 0; 613 } 614 } 615 } 616 else 617 { 618 ret.exponent += this.exponent; 619 } 620 return ret; 621 } 622 623 static if (size == 128) 624 /// 625 version(mir_bignum_test) 626 @safe pure @nogc 627 unittest 628 { 629 import mir.bignum.fixed: UInt; 630 auto fp = cast(Fp!64) Fp!128(UInt!128.fromHexString("afbbfae3cd0aff2784a1de7022b0029d")); 631 assert(fp.exponent == 64); 632 assert(fp.coefficient == UInt!64.fromHexString("afbbfae3cd0aff28")); 633 634 assert(Fp!128(-double.infinity) * Fp!128(1) == Fp!128(-double.infinity)); 635 } 636 } 637 638 /// 639 Fp!(coefficientizeA + coefficientizeB) extendedMul(bool noSpecial = false, uint coefficientizeA, uint coefficientizeB)(Fp!coefficientizeA a, Fp!coefficientizeB b) 640 @safe pure nothrow @nogc 641 { 642 import mir.bignum.fixed: extendedMul; 643 import mir.checkedint: adds; 644 645 typeof(return) ret = void; 646 ret.coefficient = extendedMul(a.coefficient, b.coefficient); 647 static if (noSpecial) 648 { 649 ret.exponent = a.exponent + b.exponent; 650 if (!ret.coefficient.signBit) 651 { 652 ret.exponent -= 1; // check overflow 653 ret.coefficient = ret.coefficient.smallLeftShift(1); 654 } 655 } 656 else 657 { 658 // nan * any -> nan 659 // inf * fin -> inf 660 if (_expect(a.isSpecial | b.isSpecial, false)) 661 { // set nan 662 ret.exponent = ret.exponent.max; 663 // nan inf case 664 if (a.isSpecial & b.isSpecial) 665 ret.coefficient = a.coefficient | b.coefficient; 666 } 667 else 668 { 669 bool overflow; 670 ret.exponent = adds(a.exponent, b.exponent, overflow); 671 // exponent underflow -> 0 or subnormal 672 // overflow -> inf 673 if (_expect(overflow, false)) 674 { 675 // overflow 676 if (a.exponent > 0) // && b.exponent > 0 is always true 677 { 678 ret.exponent = ret.exponent.max; 679 ret.coefficient = 0; 680 } 681 // underflow 682 else // a.exponent < 0 and b.exponent < 0 683 { 684 // TODO: subnormal 685 ret.exponent = 0; 686 ret.coefficient = 0; 687 } 688 } 689 else 690 if (!ret.coefficient.signBit) 691 { 692 auto normal = ret.exponent != ret.exponent.min; 693 ret.exponent -= normal; // check overflow 694 ret.coefficient = ret.coefficient.smallLeftShift(normal); 695 } 696 } 697 } 698 ret.sign = a.sign ^ b.sign; 699 return ret; 700 } 701 702 /// 703 template fp_log2(T) 704 if (is(T == float) || is(T == double) || is(T == real)) 705 { 706 /// 707 T fp_log2(uint size)(Fp!size x) 708 { 709 import mir.math.common: log2; 710 auto exponent = x.exponent + size; 711 if (!x.isSpecial) 712 x.exponent = -long(size); 713 return log2(cast(T)x) + exponent; 714 } 715 } 716 717 /// 718 version(mir_bignum_test) 719 @safe pure nothrow @nogc 720 unittest 721 { 722 import mir.math.common: log2, approxEqual; 723 import mir.bignum.fp: fp_log2; 724 725 double x = 123456789.0e+123; 726 assert(approxEqual(x.Fp!128.fp_log2!double, x.log2)); 727 } 728 729 /// 730 template fp_log(T) 731 if (is(T == float) || is(T == double) || is(T == real)) 732 { 733 /// 734 T fp_log(uint size)(Fp!size x) 735 { 736 import mir.math.constant: LN2; 737 return T(LN2) * fp_log2!T(x); 738 } 739 } 740 741 /// 742 version(mir_bignum_test) 743 @safe pure nothrow @nogc 744 unittest 745 { 746 import mir.math.common: log, approxEqual; 747 import mir.bignum.fp: fp_log; 748 749 double x = 123456789.0e+123; 750 assert(approxEqual(x.Fp!128.fp_log!double, x.log)); 751 }