1 module dcnum; 2 3 import std.traits : isIntegral, isUnsigned; 4 import std.typecons : tuple, Tuple; 5 import std.conv : to, ConvOverflowException; 6 import std.math : log10; 7 import std.format : format; 8 import std.array : join; 9 import std.algorithm : reverse, min, max, all, map; 10 import std.exception; 11 import std.random : rndGen, Random, uniform; 12 13 const BASE = 10; 14 15 class DCNumException : Exception 16 { 17 this(string msg, string file = __FILE__, size_t line = __LINE__) pure 18 { 19 super(msg, file, line); 20 } 21 } 22 23 struct DCNum 24 { 25 private: 26 bool sign; // when true this value is negative; 27 uint scale; // number of digits after the decimal point. 28 ubyte[] value; // array of digits. from higher digit to lower digit. 29 30 long len() const pure 31 { 32 return this.value.length - this.scale; 33 } 34 35 this(T)(bool sign, T scale, in ubyte[] value) pure if (isIntegral!T) 36 { 37 this.sign = sign; 38 this.scale = scale.to!uint; 39 this.value = value.dup; 40 } 41 42 this(in ubyte[] value) pure 43 { 44 this.sign = false; 45 this.scale = 0; 46 this.value = value.dup; 47 } 48 49 void rescale(uint new_scale) pure 50 { 51 if (new_scale > this.scale) 52 { 53 this.value ~= new ubyte[](new_scale - this.scale); 54 } 55 else if (new_scale < this.scale) 56 { 57 this.value.length = this.value.length - (this.scale - new_scale); 58 } 59 this.scale = new_scale; 60 } 61 62 DCNum rescaled(uint new_scale) pure const 63 { 64 DCNum copy = DCNum(this); 65 copy.rescale(new_scale); 66 return copy; 67 } 68 69 unittest 70 { 71 assert(DCNum(1).rescaled(5) == DCNum("1.00000")); 72 assert(DCNum("1.234567").rescaled(5) == DCNum("1.23456")); 73 assert(DCNum("625.0000").rescaled(1) == DCNum("625.0")); 74 } 75 76 public: 77 /// copy constructor 78 this(in DCNum num) pure 79 { 80 this.sign = num.sign; 81 this.scale = num.scale; 82 this.value = num.value.dup; 83 } 84 /// constructor from Integral values 85 this(INT)(INT val) pure if (isIntegral!INT) 86 { 87 88 // zero is special case 89 if (val == 0) 90 { 91 this.sign = false; 92 this.value = [0]; 93 this.scale = 0; 94 return; 95 } 96 97 this.sign = false; 98 static if (!isUnsigned!INT) 99 { 100 if (val < 0) 101 { 102 this.sign = true; 103 val = cast(INT)(0 - val); 104 } 105 } 106 107 uint index = 0; 108 ubyte[] buf = new ubyte[](log10(INT.max).to!uint + 1); 109 while (val != 0) 110 { 111 buf[index++] = cast(ubyte)(val % BASE); 112 val = val / BASE; 113 } 114 115 this.value = reverse(buf[0 .. index]); 116 this.scale = 0; 117 } 118 119 unittest 120 { 121 const zero = DCNum(0); 122 assert(zero.sign == false); 123 assert(zero.len == 1); 124 assert(zero.scale == 0); 125 assert(zero.value == [0]); 126 127 const one = DCNum(1); 128 assert(one.sign == false); 129 assert(one.len == 1); 130 assert(one.scale == 0); 131 assert(one.value[0] == 1); 132 133 const minus = DCNum(-1); 134 assert(minus.sign == true); 135 assert(minus.len == 1); 136 assert(minus.scale == 0); 137 assert(minus.value[0] == 1); 138 139 const large = DCNum(9876543210987654321u); 140 assert(large.sign == false); 141 assert(large.len == 19); 142 assert(large.scale == 0); 143 assert(large.value == [ 144 9, 8, 7, 6, 5, 4, 3, 2, 1, 0, 9, 8, 7, 6, 5, 4, 3, 2, 1 145 ]); 146 } 147 148 /// constructor from string. 149 this(string value) pure 150 { 151 // sign 152 ulong p = 0; 153 this.sign = false; 154 155 // "+" and "-" are invalid. so when leading + or - exists, buffer size must be longer than 1 156 if (value.length > p + 1 && value[p] == '+') 157 { 158 p++; 159 } 160 else if (value.length > p + 1 && value[p] == '-') 161 { 162 this.sign = true; 163 p++; 164 } 165 166 // skip leading zeroes 167 while (value.length > p && value[p] == '0') 168 { 169 p++; 170 } 171 172 // read values 173 this.value = new ubyte[](value.length - p); 174 long index = 0; 175 int len = 0; 176 this.scale = 0; 177 while (value.length > p) 178 { 179 if ('0' <= value[p] && value[p] <= '9') 180 { 181 this.value[index++] = cast(ubyte)(value[p] - '0'); 182 p++; 183 len++; 184 } 185 else if (value[p] == '.' && value.length > p + 1) 186 { 187 p++; 188 break; 189 } 190 else 191 { 192 throw new DCNumException("illegal character: %c".format(value[p])); 193 } 194 } 195 while (value.length > p) 196 { 197 if ('0' <= value[p] && value[p] <= '9') 198 { 199 this.value[index++] = cast(ubyte)(value[p] - '0'); 200 p++; 201 this.scale++; 202 } 203 else 204 { 205 throw new DCNumException("illegal character: %c".format(value[p])); 206 } 207 } 208 209 // shrink 210 this.value.length = (len + this.scale); 211 212 // zero 213 if (this.len == 0 && this.scale == 0) 214 { 215 this.sign = false; 216 this.scale = 0; 217 this.value = [0]; 218 } 219 } 220 221 unittest 222 { 223 const zero = DCNum("0"); 224 assert(zero.sign == false); 225 assert(zero.len == 1); 226 assert(zero.scale == 0); 227 assert(zero.value == [0]); 228 229 const minus_zero = DCNum("-0"); 230 assert(minus_zero.sign == false); 231 assert(minus_zero.len == 1); 232 assert(minus_zero.scale == 0); 233 assert(minus_zero.value == [0]); 234 235 const floating = DCNum("-1234.56789"); 236 assert(floating.sign == true); 237 assert(floating.len == 4); 238 assert(floating.scale == 5); 239 assert(floating.value == [1, 2, 3, 4, 5, 6, 7, 8, 9]); 240 241 const floating_zero = DCNum("0.56789"); 242 assert(floating_zero.sign == false); 243 assert(floating_zero.len == 0); 244 assert(floating_zero.scale == 5); 245 assert(floating_zero.value == [5, 6, 7, 8, 9]); 246 247 const trailing_zero = DCNum("-1234.567890"); 248 assert(trailing_zero.sign == true); 249 assert(trailing_zero.len == 4); 250 assert(trailing_zero.scale == 6); 251 assert(trailing_zero.value == [1, 2, 3, 4, 5, 6, 7, 8, 9, 0]); 252 253 const large = DCNum("98765432109876543210"); 254 assert(large.sign == false); 255 assert(large.len == 20); 256 assert(large.scale == 0); 257 assert(large.value == [ 258 9, 8, 7, 6, 5, 4, 3, 2, 1, 0, 9, 8, 7, 6, 5, 4, 3, 2, 1, 0 259 ]); 260 261 assertThrown!DCNumException(DCNum("+")); 262 assertThrown!DCNumException(DCNum("-")); 263 assertThrown!DCNumException(DCNum(".")); 264 assertThrown!DCNumException(DCNum("0.")); 265 assertThrown!DCNumException(DCNum("1.")); 266 assertThrown!DCNumException(DCNum("1..23")); 267 assertNotThrown!DCNumException(DCNum(".0")); 268 } 269 270 string toString() const pure 271 { 272 string s = ""; 273 if (this.sign) 274 { 275 s ~= "-"; 276 } 277 278 if (this.len == 0) 279 { 280 s ~= "0"; 281 } 282 else 283 { 284 s ~= this.value[0 .. this.len].map!(x => x.to!string).join(""); 285 } 286 287 if (this.scale > 0) 288 { 289 s ~= "."; 290 s ~= this.value[this.len .. this.len + this.scale].map!(x => x.to!string).join(""); 291 } 292 293 return s; 294 } 295 296 unittest 297 { 298 assert(DCNum("0").to!string == "0"); 299 assert(DCNum("-0").to!string == "0"); 300 assert(DCNum("100").to!string == "100"); 301 assert(DCNum("0.100").to!string == "0.100"); 302 assert(DCNum("-10.1").to!string == "-10.1"); 303 } 304 305 T to(T)() const pure if (is(T : string)) 306 { 307 return this.toString(); 308 } 309 310 T to(T)() const pure if (isIntegral!(T)) 311 { 312 if (DCNum.cmp(this, DCNum(T.max)) > 0) 313 { 314 throw new ConvOverflowException("Conversion positive overflow"); 315 } 316 static if (isUnsigned!(T)) 317 { 318 if (this.sign) 319 { 320 throw new ConvOverflowException("Conversion negative overflow"); 321 } 322 } 323 324 T v = cast(T)(0); 325 foreach (x; this.value[0 .. this.len]) 326 { 327 v = cast(T)(v * BASE + x); 328 } 329 if (this.sign) 330 { 331 v = cast(T)(0 - v); 332 } 333 return v; 334 } 335 336 unittest 337 { 338 assert(DCNum("10000").to!int == 10000); 339 assert(DCNum("10000.2345").to!int == 10000); 340 assert(DCNum("10000.9").to!int == 10000); 341 assert(DCNum("0.9").to!int == 0); 342 assert(DCNum(-1).to!byte == -1); 343 assertThrown!ConvOverflowException(DCNum(10000).to!ubyte); 344 assertThrown!ConvOverflowException(DCNum(-1).to!ubyte); 345 assertThrown!ConvOverflowException(DCNum(10000).to!ubyte); 346 } 347 348 bool isZero() pure const 349 { 350 return this.value.all!"a == 0"; 351 } 352 353 unittest 354 { 355 assert(DCNum("0").isZero); 356 assert(DCNum("-0").isZero); 357 assert((DCNum("10") - DCNum("10")).isZero); 358 assert(!(DCNum("10") - DCNum("9.99999999")).isZero); 359 } 360 361 /// if lhs is larger than rhs, return 1. if lhs is smaller than rhs, then return -1. if lhs and rhs are same, return 0; 362 private static int cmp(in DCNum lhs, in DCNum rhs, bool ignore_sign = false) pure 363 { 364 // compare sign 365 if (lhs.sign == false && rhs.sign == true && ignore_sign == false) 366 { 367 return 1; 368 } 369 else if (lhs.sign == true && rhs.sign == false && ignore_sign == false) 370 { 371 return -1; 372 } 373 int sign = (lhs.sign && ignore_sign == false) ? -1 : 1; 374 375 // compare size of integer part 376 if (lhs.len > rhs.len) 377 { 378 return sign; 379 } 380 else if (lhs.len < rhs.len) 381 { 382 return -sign; 383 } 384 385 // compare integer parts 386 foreach (i; 0 .. lhs.len) 387 { 388 if (lhs.value[i] > rhs.value[i]) 389 { 390 return sign; 391 } 392 else if (lhs.value[i] < rhs.value[i]) 393 { 394 return -sign; 395 } 396 } 397 398 // compare fraction parts 399 foreach (i; 0 .. min(lhs.scale, rhs.scale)) 400 { 401 if (lhs.value[lhs.len + i] > rhs.value[rhs.len + i]) 402 { 403 return sign; 404 } 405 else if (lhs.value[lhs.len + i] < rhs.value[rhs.len + i]) 406 { 407 return -sign; 408 } 409 } 410 411 if (lhs.scale > rhs.scale) 412 { 413 if (lhs.value[(lhs.len + rhs.scale) .. $].all!"a == 0") 414 { 415 return 0; 416 } 417 return sign; 418 } 419 else 420 { 421 if (rhs.value[(rhs.len + lhs.scale) .. $].all!"a == 0") 422 { 423 return 0; 424 } 425 return -sign; 426 } 427 } 428 429 bool opEquals(in DCNum rhs) const pure 430 { 431 return DCNum.cmp(this, rhs) == 0; 432 } 433 434 bool opEquals(T)(T rhs) const pure if (isIntegral!(T)) 435 { 436 if (this.scale != 0) 437 { 438 return false; 439 } 440 if (rhs == 0 && this.len == 1 && this.value[0] == 0) 441 { 442 return true; 443 } 444 if (rhs < 0) 445 { 446 if (this.sign == false) 447 { 448 return false; 449 } 450 rhs = -rhs; 451 } 452 453 long i = this.len - 1; 454 while (true) 455 { 456 const bool x = rhs == 0; 457 const bool y = i < 0; 458 if (x && y) 459 { 460 return true; 461 } 462 if (x || y) 463 { 464 return false; 465 } 466 467 if (rhs % BASE != this.value[i]) 468 { 469 return false; 470 } 471 rhs = rhs / BASE; 472 i--; 473 } 474 assert(false); // unreachable 475 } 476 477 int opCmp(in DCNum rhs) const pure 478 { 479 return DCNum.cmp(this, rhs); 480 } 481 482 unittest 483 { 484 assert(DCNum(0) == DCNum(0)); 485 assert(DCNum(1) > DCNum(0)); 486 assert(DCNum(1) < DCNum(2)); 487 assert(DCNum(1) > DCNum(-2)); 488 assert(DCNum(-1) > DCNum(-2)); 489 assert(DCNum("0.1") == DCNum("0.1")); 490 assert(DCNum("0.1") == DCNum("0.10")); 491 492 assert(DCNum(0) == 0); 493 assert(DCNum(1) == 1); 494 assert(DCNum(-1) == -1); 495 } 496 497 /// add lhs and rhs 498 /// this function ignore sign 499 private static DCNum add(in DCNum lhs, in DCNum rhs) pure 500 { 501 ubyte[] lhs_value = lhs.value.dup; 502 ubyte[] rhs_value = rhs.value.dup; 503 504 // align digits before the decimal point 505 if (lhs.len > rhs.len) 506 { 507 rhs_value = new ubyte[](lhs.len - rhs.len) ~ rhs_value; 508 } 509 else if (lhs.len < rhs.len) 510 { 511 lhs_value = new ubyte[](rhs.len - lhs.len) ~ lhs_value; 512 } 513 514 // align digits after the decimal point 515 if (lhs.scale > rhs.scale) 516 { 517 rhs_value = rhs_value ~ new ubyte[](lhs.scale - rhs.scale); 518 } 519 else if (lhs.scale < rhs.scale) 520 { 521 lhs_value = lhs_value ~ new ubyte[](rhs.scale - lhs.scale); 522 } 523 524 // addition digit by digit 525 ubyte carry = 0; 526 ubyte[] buf = new ubyte[](lhs_value.length); 527 foreach_reverse (i; 0 .. lhs_value.length) 528 { 529 buf[i] = (lhs_value[i] + rhs_value[i] + carry).to!(ubyte); 530 if (buf[i] >= BASE) 531 { 532 carry = 1; 533 buf[i] = buf[i] % BASE; 534 } 535 else 536 { 537 carry = 0; 538 } 539 } 540 541 const uint new_scale = max(lhs.scale, rhs.scale); 542 if (carry) 543 { 544 buf = cast(ubyte[])[1] ~ buf; 545 } 546 return DCNum(false, new_scale, buf); 547 } 548 549 /// subtract rhs from lhs. lhs must be larger than rhs on magnitude 550 private static DCNum sub(in DCNum lhs, in DCNum rhs) pure 551 { 552 ubyte[] lhs_value = lhs.value.dup; 553 554 // align digits after the decimal point 555 if (lhs.scale < rhs.scale) 556 { 557 lhs_value = lhs_value ~ new ubyte[](rhs.scale - lhs.scale); 558 } 559 560 ubyte borrow = 0; 561 const pad = lhs.len - rhs.len; 562 foreach_reverse (i; 0 .. (rhs.len + rhs.scale)) 563 { 564 if (lhs_value[pad + i] < rhs.value[i] + borrow) 565 { 566 lhs_value[pad + i] = cast(ubyte)((10 + lhs_value[pad + i]) - (rhs.value[i] + borrow)); 567 borrow = 1; 568 } 569 else 570 { 571 lhs_value[pad + i] = cast(ubyte)(lhs_value[pad + i] - (rhs.value[i] + borrow)); 572 borrow = 0; 573 } 574 } 575 if (borrow) 576 { 577 lhs_value[pad - 1]--; 578 } 579 580 uint new_scale = max(lhs.scale, rhs.scale); 581 582 // shrink value 583 uint p = 0; 584 while (p < lhs.len && lhs_value[p] == 0) 585 { 586 p++; 587 } 588 lhs_value = lhs_value[p .. $]; 589 590 return DCNum(false, new_scale, lhs_value); 591 } 592 593 DCNum opBinary(string op : "+")(in DCNum rhs) const pure 594 { 595 if (this.sign == false && rhs.sign == false) 596 { 597 return DCNum.add(this, rhs); 598 } 599 else if (this.sign == true && rhs.sign == true) 600 { 601 auto v = DCNum.add(this, rhs); 602 v.sign = true; 603 return v; 604 } 605 else 606 { 607 switch (DCNum.cmp(this, rhs, true)) 608 { 609 case 0: 610 return DCNum(0); 611 case 1: 612 auto v = DCNum.sub(this, rhs); 613 v.sign = this.sign; 614 return v; 615 case -1: 616 auto v = DCNum.sub(rhs, this); 617 v.sign = rhs.sign; 618 return v; 619 default: 620 assert(0); 621 } 622 } 623 } 624 625 unittest 626 { 627 assert((DCNum(1) + DCNum(1)).to!string == "2"); 628 assert((DCNum(1) + DCNum(2)).to!string == "3"); 629 assert((DCNum(-2) + DCNum(1)).to!string == "-1"); 630 assert((DCNum(2) + DCNum(-1)).to!string == "1"); 631 assert((DCNum(2) + DCNum(-4)).to!string == "-2"); 632 assert((DCNum(-2) + DCNum(-4)).to!string == "-6"); 633 634 assert((DCNum("0.1") + DCNum("-0.05")).to!string == "0.05"); 635 assert((DCNum("0.1") + DCNum("-0.050")).to!string == "0.050"); 636 assert((DCNum("10.1") + DCNum("-5")).to!string == "5.1"); 637 } 638 639 DCNum opBinary(string op : "-")(in DCNum rhs) pure const 640 { 641 if (this.sign == false && rhs.sign == true) 642 { 643 return DCNum.add(this, rhs); 644 } 645 else if (this.sign == true && rhs.sign == false) 646 { 647 auto v = DCNum.add(this, rhs); 648 v.sign = true; 649 return v; 650 } 651 else 652 { 653 switch (DCNum.cmp(this, rhs, true)) 654 { 655 case 0: 656 return DCNum(0); 657 case 1: 658 auto v = DCNum.sub(this, rhs); 659 v.sign = this.sign; 660 return v; 661 case -1: 662 auto v = DCNum.sub(rhs, this); 663 v.sign = !rhs.sign; 664 return v; 665 default: 666 assert(0); 667 } 668 } 669 } 670 671 unittest 672 { 673 assert((DCNum(1) - DCNum(1)).to!string == "0"); 674 assert((DCNum(1) - DCNum(2)).to!string == "-1"); 675 assert((DCNum(-2) - DCNum(1)).to!string == "-3"); 676 assert((DCNum(2) - DCNum(-1)).to!string == "3"); 677 assert((DCNum(2) - DCNum(-4)).to!string == "6"); 678 assert((DCNum(-2) - DCNum(-4)).to!string == "2"); 679 680 assert((DCNum("0.1") - DCNum("-0.05")).to!string == "0.15"); 681 assert((DCNum("0.1") - DCNum("-0.050")).to!string == "0.150"); 682 assert((DCNum("10.1") - DCNum("-5")).to!string == "15.1"); 683 } 684 685 /// simple multiply by karatsuba method. 686 /// this function ignores sign and scale, caller should set len and scale ownself 687 private static DCNum mul(in DCNum lhs, in DCNum rhs) pure 688 { 689 // utilities 690 const bytes_to_long = (in ubyte[] xs) { 691 long y = 0; 692 foreach (x; xs) 693 { 694 y = y * 10 + x; 695 } 696 return y; 697 }; 698 const long_to_bytes = (long x) { 699 ubyte[] xs = []; 700 while (x != 0) 701 { 702 xs ~= cast(ubyte)(x % 10); 703 x /= 10; 704 } 705 return reverse(xs); 706 }; 707 const max_base_size = cast(int)(log10(long.max) / 2); 708 709 // decide base size 710 const lhs_size = (lhs.len + lhs.scale + 1) / 2; 711 const rhs_size = (rhs.len + rhs.scale + 1) / 2; 712 713 if (lhs_size > max_base_size) 714 { 715 // split lhs = high || low 716 const half_len = (lhs.len + lhs.scale) / 2; 717 const lhs_low = DCNum(lhs.value[half_len .. $]); 718 const lhs_high = DCNum(lhs.value[0 .. half_len]); 719 720 // high * rhs || low * rhs 721 const low_result = DCNum.mul(lhs_low, rhs); 722 const high_result = DCNum.mul(lhs_high, rhs); 723 const high_value = high_result.value ~ new ubyte[](lhs.value.length - half_len); 724 return DCNum(high_value) + low_result; 725 } 726 if (rhs_size > max_base_size) 727 { 728 // split rhs = high || low 729 const half_len = (rhs.len + rhs.scale) / 2; 730 const rhs_low = DCNum(rhs.value[half_len .. $]); 731 const rhs_high = DCNum(rhs.value[0 .. half_len]); 732 733 // high * lhs || low * lhs 734 const low_result = DCNum.mul(rhs_low, lhs); 735 const high_result = DCNum.mul(rhs_high, lhs); 736 const high_value = high_result.value ~ new ubyte[](rhs.value.length - half_len); 737 return DCNum(high_value) + low_result; 738 } 739 740 const base_size = max(lhs_size, rhs_size); 741 742 // if both are small, calculate as long 743 if (lhs.len + lhs.scale + rhs.len + rhs.scale < max_base_size) 744 { 745 const long x = bytes_to_long(lhs.value); 746 const long y = bytes_to_long(rhs.value); 747 const long z = x * y; 748 ubyte[] buf = long_to_bytes(z); 749 const uint new_scale = lhs.scale + rhs.scale; 750 return DCNum(buf); 751 } 752 753 assert(base_size <= max_base_size); // TODO 754 755 // split x -> x1 || x0 756 const long x0 = bytes_to_long(lhs.value[max(cast(int)($ - base_size), 0) .. $]); 757 const long x1 = bytes_to_long(lhs.value[0 .. max(cast(int)($ - base_size), 0)]); 758 const long y0 = bytes_to_long(rhs.value[max(cast(int)($ - base_size), 0) .. $]); 759 const long y1 = bytes_to_long(rhs.value[0 .. max(cast(int)($ - base_size), 0)]); 760 761 // karatsuba algorithm 762 const long z0 = x0 * y0; 763 const long z2 = x1 * y1; 764 const long z1 = z2 + z0 - (x1 - x0) * (y1 - y0); // buf0 = z0, buf1 = z1 * base, buf2 = z2 * base^2 765 const buf0 = long_to_bytes(z0); 766 const buf1 = long_to_bytes(z1) ~ new ubyte[](base_size); 767 const buf2 = long_to_bytes(z2) ~ new ubyte[](base_size * 2); // convert to DCNum 768 const v0 = DCNum(buf0); 769 const v1 = DCNum(buf1); 770 const v2 = DCNum(buf2); // summing 771 auto v = v0 + v1 + v2; 772 773 // remove reading zeroes 774 long p = 0; 775 while (p < v.len && v.value[p] == 0) 776 { 777 p++; 778 } 779 v.value = v.value[p .. $]; 780 return v; 781 } 782 783 /// multiply lhs by integer rhs 784 /// this function ignores len and scale properties 785 private static DCNum mul_by_const(T)(in DCNum lhs, T rhs) pure 786 if (isIntegral!(T)) 787 in 788 { 789 assert((BASE - 1) * rhs <= ubyte.max); 790 } 791 do 792 { 793 if (rhs == 0) 794 { 795 return DCNum(0); 796 } 797 if (rhs == 1) 798 { 799 return DCNum(lhs); 800 } 801 802 ubyte[] buf = new ubyte[](lhs.value.length + 1); 803 int carry = 0; 804 long i = lhs.value.length; 805 foreach_reverse (x; lhs.value) 806 { 807 const y = x * rhs + carry; 808 buf[i--] = (y % BASE).to!ubyte; 809 carry = cast(int)(y / BASE); 810 } 811 if (carry != 0) 812 { 813 buf[i--] = carry.to!ubyte; 814 } 815 return DCNum(buf[i + 1 .. $]); 816 } 817 818 unittest 819 { 820 alias TestCase = Tuple!(DCNum, int, string); 821 auto testcases = [ 822 TestCase(DCNum(500), 2, "1000"), 823 TestCase(DCNum("33333333"), 3, "99999999"), 824 TestCase(DCNum("123456789"), 5, "617283945") 825 ]; 826 foreach (t; testcases) 827 { 828 auto r = DCNum.mul_by_const(t[0], t[1]).to!string; 829 assert(r == t[2], "Case: %s, Got: %s".format(t, r)); 830 } 831 } 832 833 DCNum mul(in DCNum rhs, uint scale) const pure 834 { 835 if (this.isZero || rhs.isZero) 836 { 837 return DCNum(0); 838 } 839 840 auto v = DCNum.mul(this, rhs); 841 int v_scale = this.scale + rhs.scale; 842 if (v.value.length < v_scale) 843 { 844 v.value = new ubyte[](v_scale - v.value.length) ~ v.value; 845 } 846 v.scale = max(scale, min(this.scale, rhs.scale)); 847 if (v.scale < v_scale) 848 { 849 v.value.length = v.value.length - (v_scale - v.scale); 850 } 851 852 if (this.sign != rhs.sign) 853 { 854 v.sign = true; 855 } 856 return v; 857 } 858 859 unittest 860 { 861 import std.stdio; 862 863 assert(DCNum("1.00").mul(DCNum("0.5"), 2) == DCNum("0.50")); 864 assert(DCNum("1.00").mul(DCNum("0.5"), 3) == DCNum("0.50")); 865 assert(DCNum("1.000").mul(DCNum("0.5"), 2) == DCNum("0.50")); 866 assert(DCNum("3.00").mul(DCNum("0.5"), 2) == DCNum("1.50")); 867 } 868 869 DCNum opBinary(string op : "*")(in DCNum rhs) const pure 870 { 871 return this.mul(rhs, this.scale + rhs.scale); 872 } 873 874 unittest 875 { 876 // normal values 877 assert((DCNum(0) * DCNum(1)).to!string == "0"); 878 assert((DCNum(1) * DCNum(1)).to!string == "1"); 879 assert((DCNum(-1) * DCNum(1)).to!string == "-1"); 880 assert((DCNum(-1) * DCNum(-1)).to!string == "1"); 881 assert((DCNum(100000) * DCNum(70)).to!string == "7000000"); 882 assert((DCNum("987654321123456789") * DCNum("100000000")) 883 .to!string == "98765432112345678900000000"); 884 885 // decimal values 886 assert((DCNum("2") * DCNum("0.5")) == DCNum("1.0")); 887 assert((DCNum("100000.123") * DCNum("100")).to!string == "10000012.300"); 888 assert((DCNum("0.123") * DCNum("0.01")).to!string == "0.00123"); // very large number 889 assert((DCNum("100000.123") * DCNum("0.01")).to!string == "1000.00123"); // very large number 890 assert((DCNum("9876543211234567899") * DCNum("100000000")) 891 .to!string == "987654321123456789900000000"); 892 893 // rsa 894 const p = DCNum("7857843807357322428021380248993576655206988614176418792379176652835565059295420706559476442876718401226381634767797717201944708260927696952220575206571167"); 895 const q = DCNum("11022865926124182806180388767382016652191532553299862348953322076764410256860215835703669245291707730752129977734682225290368939156485722324614823488258901"); 896 assert((p * q).to!string == "86615958756924946592957282448568720038805999499540908216698775245619824596674378512195525165203154029569489225605263626685364659699870945114711447932248705556536031296400659122760841627071717950914771235328300476962435317906251410048014717963467669606882231758796075711787284426301244369129372556726977707467"); 897 } 898 899 /// divide lhs by rhs. this function ignores its scale and sign 900 /// this function using Knuth's Algorithm D 901 private static DCNum div(in DCNum lhs, in DCNum rhs) pure 902 { 903 if (rhs.isZero) 904 { 905 throw new DCNumException("Divide by zero"); 906 } 907 if (lhs.isZero) 908 { 909 return DCNum(0); 910 } 911 912 DCNum l = DCNum(lhs), r = DCNum(rhs); 913 l.sign = false; 914 r.sign = false; 915 l.scale = 0; 916 r.scale = 0; 917 918 if (DCNum.cmp(l, r) < 0) 919 { 920 return DCNum(0); 921 } 922 while (l.value.length <= 2 || r.value.length <= 2) 923 { 924 l.value ~= [0]; 925 r.value ~= [0]; 926 } 927 928 // normalize 929 auto d = ((BASE - 1) / r.value[0]).to!long; 930 auto rr = DCNum.mul_by_const(r, d); 931 if (rr.value[0] < BASE / 2) 932 { 933 d = (BASE / (r.value[0] + 1)).to!long; 934 } 935 if (d != 1) 936 { 937 l = DCNum.mul_by_const(l, d); 938 r = DCNum.mul_by_const(r, d); 939 } 940 941 ubyte[] q = []; 942 const n = l.len - r.len; 943 for (long j = n; j >= 0; j--) 944 { 945 // guess q 946 auto qguess = (l.value[0] * BASE + l.value[1]) / r.value[0]; 947 auto rguess = (l.value[0] * BASE + l.value[1]) - qguess * r.value[0]; 948 if (qguess == BASE && rguess == 0) 949 { 950 auto x = DCNum(r); 951 x.value ~= new ubyte[](j); 952 if (DCNum.cmp(l, x) == 0) 953 { 954 q ~= cast(ubyte[])[1] ~ new ubyte[](j); 955 break; 956 } 957 } 958 while (rguess < BASE && (qguess >= BASE || qguess * r.value[1] 959 > BASE * rguess + l.value[2])) 960 { 961 qguess--; 962 rguess += r.value[0]; 963 } 964 965 // multiple and substract 966 auto x = DCNum.mul_by_const(r, qguess); 967 if (j > 0) 968 { 969 x.value ~= new ubyte[](j - 1); 970 } 971 972 // fix guess 973 while (DCNum.cmp(l, x) < 0) 974 { 975 qguess--; 976 x = DCNum.mul_by_const(r, qguess); 977 if (j > 0) 978 { 979 x.value ~= new ubyte[](j - 1); 980 } 981 } 982 983 l = DCNum.sub(l, x); 984 q ~= qguess.to!ubyte; 985 if (DCNum.cmp(l, r) <= 0) 986 { 987 if (j > 0) 988 { 989 q ~= new ubyte[](j - 1); 990 } 991 break; 992 } 993 } 994 995 return DCNum(q); 996 } 997 998 unittest 999 { 1000 assert(DCNum.div(DCNum(10), DCNum(5)) == 2); 1001 assert(DCNum.div(DCNum(1000), DCNum(5)) == 200); 1002 assert(DCNum.div(DCNum(1000), DCNum(50)) == 20); 1003 assert(DCNum.div(DCNum("100000000000"), DCNum(5)).to!string == "20000000000"); 1004 assert(DCNum.div(DCNum("11044102452163169934"), 1005 DCNum("10319522097943752571")).to!string == "1"); 1006 1007 // random case 1008 foreach (_; 0 .. 100) 1009 { 1010 auto rnd = rndGen(); 1011 auto x = uniform!ulong(rnd); 1012 auto y = uniform!ulong(rnd); 1013 auto r = DCNum.div(DCNum(x), DCNum(y)).to!string; 1014 assert(r == (x / y).to!string, "Case: %d / %d, Got: %s".format(x, y, r)); 1015 } 1016 } 1017 1018 DCNum div(in DCNum rhs, uint scale) const pure 1019 { 1020 if (this.isZero) 1021 { 1022 return DCNum(0); 1023 } 1024 DCNum lhs = DCNum(this); 1025 lhs.value ~= new ubyte[](scale + rhs.scale); 1026 lhs.scale = 0; 1027 1028 auto v = DCNum.div(lhs, rhs); 1029 v.value.length = v.value.length - this.scale; 1030 if (v.value.length < scale) 1031 { 1032 v.value = new ubyte[](scale - v.value.length) ~ v.value; 1033 } 1034 v.scale = scale; 1035 if (this.sign != rhs.sign) 1036 { 1037 v.sign = true; 1038 } 1039 return v; 1040 } 1041 1042 unittest 1043 { 1044 assert(DCNum(2).div(DCNum(2), 10).to!string == "1.0000000000"); 1045 assert(DCNum(10).div(DCNum(3), 10).to!string == "3.3333333333"); 1046 assert(DCNum(10).div(DCNum("3.0"), 10).to!string == "3.3333333333"); 1047 assert(DCNum(10).div(DCNum("3.3"), 10).to!string == "3.0303030303"); 1048 assert(DCNum(10).div(DCNum("3.333333333333333333"), 10).to!string == "3.0000000000"); 1049 assert(DCNum(5).div(DCNum(10), 10).to!string == "0.5000000000"); 1050 assert(DCNum("123456789.0987654321").div(DCNum("555"), 10).to!string == "222444.6650428205"); 1051 assert(DCNum("10.000").div(DCNum("2"), 1).to!string == "5.0"); 1052 assert(DCNum("10.000").div(DCNum("-2"), 1).to!string == "-5.0"); 1053 assert(DCNum("-10.000").div(DCNum("2"), 1).to!string == "-5.0"); 1054 assert(DCNum("138458412558.000000").div(DCNum("74.4200006"), 5) == DCNum("1860500019.37248")); 1055 assert(DCNum("33333333333333333333333333") 1056 .div(DCNum("248352686608866080.9427714159"), 11) == DCNum("134217727.97582189752")); 1057 assert(DCNum(2).div(DCNum("1.5"), 1) == DCNum("1.3")); 1058 assert(DCNum(1).div(DCNum("625.000000"), 4) == DCNum("0.0016")); 1059 } 1060 1061 DCNum opBinary(string op : "/")(in DCNum rhs) const pure 1062 { 1063 return this.div(rhs, max(this.scale, rhs.scale)); 1064 } 1065 1066 unittest 1067 { 1068 assert((DCNum(5) / DCNum(3)).to!string == "1"); 1069 assert((DCNum(10) / DCNum(3)).to!string == "3"); 1070 assert((DCNum(10) / DCNum("3.0")).to!string == "3.3"); 1071 assert((DCNum(-10) / DCNum("0.50")).to!string == "-20.00"); 1072 const n = DCNum("86615958756924946592957282448568720038805999499540908216698775245619824596674378512195525165203154029569489225605263626685364659699870945114711447932248705556536031296400659122760841627071717950914771235328300476962435317906251410048014717963467669606882231758796075711787284426301244369129372556726977707467"); 1073 const q = DCNum("11022865926124182806180388767382016652191532553299862348953322076764410256860215835703669245291707730752129977734682225290368939156485722324614823488258901"); 1074 assert((n / q).to!string == "7857843807357322428021380248993576655206988614176418792379176652835565059295420706559476442876718401226381634767797717201944708260927696952220575206571167"); 1075 } 1076 1077 DCNum mod(in DCNum rhs, uint scale) const pure 1078 { 1079 auto r = this.div(rhs, scale); 1080 auto m = this - r * rhs; 1081 return m; 1082 } 1083 1084 unittest 1085 { 1086 assert(DCNum(10).mod(DCNum(3), 0) == 1); 1087 assert(DCNum(10).mod(DCNum(3), 2) == DCNum("0.01")); 1088 assert(DCNum(10).mod(DCNum("3.0"), 0) == DCNum("1.0")); 1089 assert(DCNum(10).mod(DCNum("3.0"), 2) == DCNum("0.010")); 1090 assert(DCNum(10).mod(DCNum("5.0"), 1000) == DCNum(0)); 1091 assert(DCNum(-2).mod(DCNum("1.60"), 0) == DCNum("-0.40")); 1092 } 1093 1094 DCNum opBinary(string op : "%")(in DCNum rhs) const pure 1095 { 1096 return this.mod(rhs, max(this.scale, rhs.scale)); 1097 } 1098 1099 unittest 1100 { 1101 assert(DCNum(10) % DCNum(3) == 1); 1102 assert(DCNum(-10) % DCNum(3) == -1); 1103 assert(DCNum(10) % DCNum(-3) == 1); 1104 assert(DCNum(10) % DCNum("3.0") == DCNum("0.1")); 1105 assert(DCNum(10) % DCNum("3.000") == DCNum("0.001")); 1106 } 1107 1108 /// Find square root by Newton's algorithm 1109 /// this function assumes this >= 0 1110 /// the returned value has specified scale 1111 DCNum sqrt(uint scale) const pure 1112 { 1113 // check this 1114 if (DCNum.cmp(this, DCNum(0)) <= 0) 1115 { 1116 throw new DCNumException("negative or zero value given for sqrt"); 1117 } 1118 1119 // guess the start value 1120 DCNum guess; 1121 final switch (DCNum.cmp(this, DCNum(1))) 1122 { 1123 case 0: 1124 return DCNum(1).rescaled(scale); 1125 case -1: 1126 guess = DCNum(1); 1127 break; 1128 case 1: 1129 guess = DCNum(cast(ubyte[])([1]) ~ new ubyte[](this.len / 2)); 1130 break; 1131 } 1132 1133 // newton's algorithm 1134 while (true) 1135 { 1136 // update guess 1137 DCNum new_guess = (this.div(guess, scale + 1) + guess).mul(DCNum("0.5"), scale); 1138 DCNum diff = new_guess - guess; 1139 guess = DCNum(new_guess); 1140 // check diff is near the zero 1141 if (diff.value[0 .. $ - 1].all!"a == 0" && diff.value[$ - 1] <= 0) 1142 { 1143 break; 1144 } 1145 } 1146 1147 guess.rescale(scale); 1148 return guess; 1149 } 1150 1151 unittest 1152 { 1153 assertThrown!DCNumException(DCNum(0).sqrt(1)); 1154 assertThrown!DCNumException(DCNum(-4).sqrt(1)); 1155 1156 assert(DCNum(2).sqrt(0) == DCNum("1")); 1157 assert(DCNum(2).sqrt(1) == DCNum("1.4")); 1158 assert(DCNum(2).sqrt(5) == DCNum("1.41421")); 1159 assert(DCNum(2).sqrt(10) == DCNum("1.4142135623")); 1160 1161 assert(DCNum(4).sqrt(0) == DCNum("2")); 1162 assert(DCNum(16).sqrt(1) == DCNum("4.0")); 1163 assert(DCNum("9.99").sqrt(5) == DCNum("3.16069")); 1164 assert(DCNum("33333333333333333333333333").sqrt(10) == DCNum("5773502691896.2576450914")); 1165 } 1166 1167 /// raise this to the exponent power 1168 DCNum raise(long exponent, uint scale) pure const 1169 { 1170 if (exponent == 0) 1171 { 1172 return DCNum(1); 1173 } 1174 1175 bool neg; 1176 if (exponent < 0) 1177 { 1178 neg = true; 1179 exponent = -exponent; 1180 } 1181 else 1182 { 1183 neg = false; 1184 scale = min(this.scale * exponent, max(this.scale, scale)).to!uint; 1185 } 1186 1187 DCNum pow = DCNum(this); 1188 uint pow_scale = this.scale; 1189 while ((exponent & 1) == 0) 1190 { 1191 pow_scale <<= 1; 1192 pow = pow.mul(pow, pow_scale); 1193 exponent >>= 1; 1194 } 1195 1196 DCNum result = DCNum(pow); 1197 uint calc_scale = pow_scale; 1198 exponent >>= 1; 1199 while (exponent > 0) 1200 { 1201 pow_scale <<= 1; 1202 pow = pow.mul(pow, pow_scale); 1203 if (exponent & 1) 1204 { 1205 calc_scale += pow_scale; 1206 result = result.mul(pow, calc_scale); 1207 } 1208 exponent >>= 1; 1209 } 1210 1211 if (neg) 1212 { 1213 return DCNum(1).div(result, scale); 1214 } 1215 else 1216 { 1217 return result.rescaled(scale); 1218 } 1219 } 1220 1221 unittest 1222 { 1223 import std.stdio; 1224 1225 assert(DCNum(2).raise(3, 0) == DCNum(8)); 1226 assert(DCNum("5.0").raise(4, 0) == DCNum("625.0")); 1227 assert(DCNum("1.234").raise(5, 0) == DCNum("2.861")); 1228 assert(DCNum("1.234").raise(5, 10) == DCNum("2.8613817210")); 1229 assert(DCNum("1.234").raise(5, 20) == DCNum("2.861381721051424")); 1230 1231 assert(DCNum(2).raise(-3, 0) == DCNum("0")); 1232 assert(DCNum(2).raise(-3, 1) == DCNum("0.1")); 1233 assert(DCNum(2).raise(-3, 5) == DCNum("0.12500")); 1234 assert(DCNum("5.0000000").raise(-4, 0) == DCNum("0")); 1235 assert(DCNum("5.0000000").raise(-4, 10) == DCNum("0.0016000000")); 1236 } 1237 }