1/* 2 * Copyright (C) 1987-2008 Sun Microsystems, Inc. All Rights Reserved. 3 * Copyright (C) 2008-2012 Robert Ancell 4 * 5 * This program is free software: you can redistribute it and/or modify it under 6 * the terms of the GNU General Public License as published by the Free Software 7 * Foundation, either version 3 of the License, or (at your option) any later 8 * version. See http://www.gnu.org/copyleft/gpl.html the full text of the 9 * license. 10 */ 11 12/* This maths library is based on the MP multi-precision floating-point 13 * arithmetic package originally written in FORTRAN by Richard Brent, 14 * Computer Centre, Australian National University in the 1970's. 15 * 16 * It has been converted from FORTRAN into C using the freely available 17 * f2c translator, available via netlib on research.att.com. 18 * 19 * The subsequently converted C code has then been tidied up, mainly to 20 * remove any dependencies on the libI77 and libF77 support libraries. 21 * 22 * FOR A GENERAL DESCRIPTION OF THE PHILOSOPHY AND DESIGN OF MP, 23 * SEE - R. P. BRENT, A FORTRAN MULTIPLE-PRECISION ARITHMETIC 24 * SOME ADDITIONAL DETAILS ARE GIVEN IN THE SAME ISSUE, 71-81. 25 * FOR DETAILS OF THE IMPLEMENTATION, CALLING SEQUENCES ETC. SEE 26 * THE MP USERS GUIDE. 27 */ 28 29using MPC; 30 31private delegate int BitwiseFunc (int v1, int v2); 32 33public enum AngleUnit 34{ 35 RADIANS, 36 DEGREES, 37 GRADIANS 38} 39 40/* Object for a high precision floating point number representation */ 41public class Number : Object 42{ 43 /* real and imaginary part of a Number */ 44 private Complex num = Complex (precision); 45 46 public static MPFR.Precision precision { get; set; default = 1000; } 47 48 /* Stores the error msg if an error occurs during calculation. Otherwise should be null */ 49 public static string? error { get; set; default = null; } 50 51 public Number.integer (int64 real, int64 imag = 0) 52 { 53 num.set_signed_integer ((long) real, (long) imag); 54 } 55 56 public Number.unsigned_integer (uint64 real, uint64 imag = 0) 57 { 58 num.set_unsigned_integer ((ulong) real, (ulong) imag); 59 } 60 61 public Number.fraction (int64 numerator, int64 denominator) 62 { 63 if (denominator < 0) 64 { 65 numerator = -numerator; 66 denominator = -denominator; 67 } 68 69 this.integer (numerator); 70 if (denominator != 1) 71 { 72 num.divide_unsigned_integer (num, (long) denominator); 73 } 74 } 75 76 /* Helper constructor. Creates new Number from already existing MPFR.Real. */ 77 public Number.mpreal (MPFR.Real real, MPFR.Real? imag = null) 78 { 79 num.set_mpreal (real, imag); 80 } 81 82 public Number.double (double real, double imag = 0) 83 { 84 num.set_double (real, imag); 85 } 86 87 public Number.complex (Number r, Number i) 88 { 89 num.set_mpreal (r.num.get_real ().val, i.num.get_real ().val); 90 } 91 92 public Number.polar (Number r, Number theta, AngleUnit unit = AngleUnit.RADIANS) 93 { 94 var x = theta.cos (unit); 95 var y = theta.sin (unit); 96 this.complex (x.multiply (r), y.multiply (r)); 97 } 98 99 public Number.eulers () 100 { 101 num.get_real ().val.set_unsigned_integer (1); 102 /* e^1, since mpfr doesn't have a function to return e */ 103 num.get_real ().val.exp (num.get_real ().val); 104 num.get_imag ().val.set_zero (); 105 } 106 107 public Number.i () 108 { 109 num.set_signed_integer (0, 1); 110 } 111 112 public Number.pi () 113 { 114 num.get_real ().val.const_pi (); 115 num.get_imag ().val.set_zero (); 116 } 117 118 public Number.tau () 119 { 120 num.get_real ().val.const_tau (); 121 num.get_imag ().val.set_zero (); 122 } 123 124 /* Sets z to be a uniform random number in the range [0, 1] */ 125 public Number.random () 126 { 127 this.double (Random.next_double ()); 128 } 129 130 public int64 to_integer () 131 { 132 return num.get_real ().val.get_signed_integer (); 133 } 134 135 public uint64 to_unsigned_integer () 136 { 137 return num.get_real ().val.get_unsigned_integer (); 138 } 139 140 public float to_float () 141 { 142 return num.get_real ().val.get_float (MPFR.Round.NEAREST); 143 } 144 145 public double to_double () 146 { 147 return num.get_real ().val.get_double (MPFR.Round.NEAREST); 148 } 149 150 /* Return true if the value is x == 0 */ 151 public bool is_zero () 152 { 153 return num.is_zero (); 154 } 155 156 /* Return true if x < 0 */ 157 public bool is_negative () 158 { 159 return num.get_real ().val.sgn () < 0; 160 } 161 162 /* Return true if x is integer */ 163 public bool is_integer () 164 { 165 if (is_complex ()) 166 return false; 167 168 return num.get_real ().val.is_integer () != 0; 169 } 170 171 /* Return true if x is a positive integer */ 172 public bool is_positive_integer () 173 { 174 if (is_complex ()) 175 return false; 176 else 177 return num.get_real ().val.sgn () >= 0 && is_integer (); 178 } 179 180 /* Return true if x is a natural number (an integer ≥ 0) */ 181 public bool is_natural () 182 { 183 if (is_complex ()) 184 return false; 185 else 186 return num.get_real ().val.sgn () > 0 && is_integer (); 187 } 188 189 /* Return true if x has an imaginary component */ 190 public bool is_complex () 191 { 192 return !num.get_imag ().val.is_zero (); 193 } 194 195 /* Return error if overflow or underflow */ 196 public static void check_flags () 197 { 198 if (MPFR.mpfr_is_underflow () != 0) 199 { 200 /* Translators: Error displayed when underflow error occured */ 201 error = _("Underflow error"); 202 } 203 else if (MPFR.mpfr_is_overflow () != 0) 204 { 205 /* Translators: Error displayed when overflow error occured */ 206 error = _("Overflow error"); 207 } 208 } 209 210 /* Return true if x == y */ 211 public bool equals (Number y) 212 { 213 return num.is_equal (y.num); 214 } 215 216 /* Returns: 217 * 0 if x == y 218 * <0 if x < y 219 * >0 if x > y 220 */ 221 public int compare (Number y) 222 { 223 return num.get_real ().val.cmp (y.num.get_real ().val); 224 } 225 226 /* Sets z = sgn (x) */ 227 public Number sgn () 228 { 229 var z = new Number.integer (num.get_real ().val.sgn ()); 230 return z; 231 } 232 233 /* Sets z = −x */ 234 public Number invert_sign () 235 { 236 var z = new Number (); 237 z.num.neg (num); 238 return z; 239 } 240 241 /* Sets z = |x| */ 242 public Number abs () 243 { 244 var z = new Number (); 245 z.num.get_imag ().val.set_zero (); 246 MPC.abs (z.num.get_real ().val, num); 247 return z; 248 } 249 250 /* Sets z = Arg (x) */ 251 public Number arg (AngleUnit unit = AngleUnit.RADIANS) 252 { 253 if (is_zero ()) 254 { 255 /* Translators: Error display when attempting to take argument of zero */ 256 error = _("Argument not defined for zero"); 257 return new Number.integer (0); 258 } 259 var z = new Number (); 260 z.num.get_imag ().val.set_zero (); 261 MPC.arg (z.num.get_real ().val, num); 262 mpc_from_radians (z.num, z.num, unit); 263 // MPC returns -π for the argument of negative real numbers if 264 // their imaginary part is -0 (which it is in the numbers 265 // created by test-equation), we want +π for all real negative 266 // numbers 267 if (!is_complex () && is_negative ()) 268 z.num.get_real ().val.abs (z.num.get_real ().val); 269 270 return z; 271 } 272 273 /* Sets z = ‾̅x */ 274 public Number conjugate () 275 { 276 var z = new Number (); 277 z.num.conj (num); 278 return z; 279 } 280 281 /* Sets z = Re (x) */ 282 public Number real_component () 283 { 284 var z = new Number (); 285 z.num.set_mpreal (num.get_real ().val); 286 return z; 287 } 288 289 /* Sets z = Im (x) */ 290 public Number imaginary_component () 291 { 292 /* Copy imaginary component to real component */ 293 var z = new Number (); 294 z.num.set_mpreal (num.get_imag ().val); 295 return z; 296 } 297 298 public Number integer_component () 299 { 300 var z = new Number (); 301 z.num.get_imag ().val.set_zero (); 302 z.num.get_real ().val.trunc (num.get_real ().val); 303 return z; 304 } 305 306 /* Sets z = x mod 1 */ 307 public Number fractional_component () 308 { 309 var z = new Number (); 310 z.num.get_imag ().val.set_zero (); 311 z.num.get_real ().val.frac (num.get_real ().val); 312 return z; 313 } 314 315 /* Sets z = {x} */ 316 public Number fractional_part () 317 { 318 return subtract (floor ()); 319 } 320 321 /* Sets z = ⌊x⌋ */ 322 public Number floor () 323 { 324 var z = new Number (); 325 z.num.get_imag ().val.set_zero (); 326 z.num.get_real ().val.floor (num.get_real ().val); 327 return z; 328 } 329 330 /* Sets z = ⌈x⌉ */ 331 public Number ceiling () 332 { 333 var z = new Number (); 334 z.num.get_imag ().val.set_zero (); 335 z.num.get_real ().val.ceil (num.get_real ().val); 336 return z; 337 } 338 339 /* Sets z = [x] */ 340 public Number round () 341 { 342 var z = new Number (); 343 z.num.get_imag ().val.set_zero (); 344 z.num.get_real ().val.round (num.get_real ().val); 345 return z; 346 } 347 348 /* Sets z = 1 ÷ x */ 349 public Number reciprocal () 350 { 351 var z = new Number (); 352 z.num.set_signed_integer (1); 353 z.num.mpreal_divide (z.num.get_real ().val, num); 354 return z; 355 } 356 357 /* Sets z = e^x */ 358 public Number epowy () 359 { 360 var z = new Number (); 361 z.num.exp (num); 362 return z; 363 } 364 365 /* Sets z = x^y */ 366 public Number xpowy (Number y) 367 { 368 /* 0^-n invalid */ 369 if (is_zero () && y.is_negative ()) 370 { 371 /* Translators: Error displayed when attempted to raise 0 to a negative re_exponent */ 372 error = _("The power of zero is undefined for a negative exponent"); 373 return new Number.integer (0); 374 } 375 376 /* 0^0 is indeterminate */ 377 if (is_zero () && y.is_zero ()) 378 { 379 /* Translators: Error displayed when attempted to raise 0 to power of zero */ 380 error = _("Zero raised to zero is undefined"); 381 return new Number.integer (0); 382 } 383 if (!is_complex () && !y.is_complex () && !y.is_integer ()) 384 { 385 var reciprocal = y.reciprocal (); 386 if (reciprocal.is_integer ()) 387 return root (reciprocal.to_integer ()); 388 } 389 390 var z = new Number (); 391 z.num.power (num, y.num); 392 return z; 393 } 394 395 /* Sets z = x^y */ 396 public Number xpowy_integer (int64 n) 397 { 398 /* 0^-n invalid */ 399 if (is_zero () && n < 0) 400 { 401 /* Translators: Error displayed when attempted to raise 0 to a negative re_exponent */ 402 error = _("The power of zero is undefined for a negative exponent"); 403 return new Number.integer (0); 404 } 405 406 /* 0^0 is indeterminate */ 407 if (is_zero () && n == 0) 408 { 409 /* Translators: Error displayed when attempted to raise 0 to power of zero */ 410 error = _("Zero raised to zero is undefined"); 411 return new Number.integer (0); 412 } 413 var z = new Number (); 414 z.num.power_integer (num, (long) n); 415 return z; 416 } 417 418 /* Sets z = n√x */ 419 public Number root (int64 n) 420 { 421 uint64 p; 422 var z = new Number (); 423 if (n < 0) 424 { 425 z.num.unsigned_integer_divide (1, num); 426 if (n == int64.MIN) 427 p = (uint64) int64.MAX + 1; 428 else 429 p = -n; 430 } else if (n > 0) { 431 z.num.@set (num); 432 p = n; 433 } else { 434 error = _("The zeroth root of a number is undefined"); 435 return new Number.integer (0); 436 } 437 438 if (!is_complex () && (!is_negative () || (p & 1) == 1)) 439 { 440 z.num.get_real ().val.root (z.num.get_real ().val, (ulong) p); 441 z.num.get_imag().val.set_zero(); 442 } else { 443 var tmp = MPFR.Real (precision); 444 tmp.set_unsigned_integer ((ulong) p); 445 tmp.unsigned_integer_divide (1, tmp); 446 z.num.power_mpreal (z.num, tmp); 447 } 448 return z; 449 } 450 451 /* Sets z = √x */ 452 public Number sqrt () 453 { 454 return root(2); 455 } 456 457 /* Sets z = ln x */ 458 public Number ln () 459 { 460 /* ln (0) undefined */ 461 if (is_zero ()) 462 { 463 /* Translators: Error displayed when attempting to take logarithm of zero */ 464 error = _("Logarithm of zero is undefined"); 465 return new Number.integer (0); 466 } 467 468 /* ln (-x) complex */ 469 /* FIXME: Make complex numbers optional */ 470 /*if (is_negative ()) 471 { 472 // Translators: Error displayed attempted to take logarithm of negative value 473 mperr (_("Logarithm of negative values is undefined")); 474 return new Number.integer (0); 475 }*/ 476 477 var z = new Number (); 478 z.num.log (num); 479 // MPC returns -π for the imaginary part of the log of 480 // negative real numbers if their imaginary part is -0 (which 481 // it is in the numbers created by test-equation), we want +π 482 if (!is_complex () && is_negative ()) 483 z.num.get_imag ().val.abs (z.num.get_imag ().val); 484 485 return z; 486 } 487 488 /* Sets z = log_n x */ 489 public Number logarithm (int64 n) 490 { 491 /* log (0) undefined */ 492 if (is_zero ()) 493 { 494 /* Translators: Error displayed when attempting to take logarithm of zero */ 495 error = _("Logarithm of zero is undefined"); 496 return new Number.integer (0); 497 } 498 499 /* logn (x) = ln (x) / ln (n) */ 500 var t1 = new Number.integer (n); 501 return ln ().divide (t1.ln ()); 502 } 503 504 /* Sets z = x! */ 505 public Number factorial () 506 { 507 /* 0! == 1 */ 508 if (is_zero ()) 509 return new Number.integer (1); 510 if (!is_natural ()) 511 { 512 513 /* Factorial Not defined for Complex or for negative numbers */ 514 if (is_negative () || is_complex ()) 515 { 516 /* Translators: Error displayed when attempted take the factorial of a negative or complex number */ 517 error = _("Factorial is only defined for non-negative real numbers"); 518 return new Number.integer (0); 519 } 520 521 var tmp = add (new Number.integer (1)); 522 var tmp2 = MPFR.Real (precision); 523 524 /* Factorial(x) = Gamma(x+1) - This is the formula used to calculate Factorial.*/ 525 tmp2.gamma (tmp.num.get_real ().val); 526 527 return new Number.mpreal (tmp2); 528 } 529 530 /* Convert to integer - if couldn't be converted then the factorial would be too big anyway */ 531 var value = to_integer (); 532 var z = this; 533 for (var i = 2; i < value; i++) 534 z = z.multiply_integer (i); 535 536 return z; 537 } 538 539 /* Sets z = x + y */ 540 public Number add (Number y) 541 { 542 var z = new Number (); 543 z.num.add (num, y.num); 544 return z; 545 } 546 547 /* Sets z = x − y */ 548 public Number subtract (Number y) 549 { 550 var z = new Number (); 551 z.num.subtract (num, y.num); 552 return z; 553 } 554 555 /* Sets z = x × y */ 556 public Number multiply (Number y) 557 { 558 var z = new Number (); 559 z.num.multiply (num, y.num); 560 return z; 561 } 562 563 /* Sets z = x × y */ 564 public Number multiply_integer (int64 y) 565 { 566 var z = new Number (); 567 z.num.multiply_signed_integer (num, (long) y); 568 return z; 569 } 570 571 /* Sets z = x ÷ y */ 572 public Number divide (Number y) 573 { 574 if (y.is_zero ()) 575 { 576 /* Translators: Error displayed attempted to divide by zero */ 577 error = _("Division by zero is undefined"); 578 return new Number.integer (0); 579 } 580 581 var z = new Number (); 582 z.num.divide (num, y.num); 583 return z; 584 } 585 586 /* Sets z = x ÷ y */ 587 public Number divide_integer (int64 y) 588 { 589 return divide (new Number.integer (y)); 590 } 591 592 /* Sets z = x mod y */ 593 public Number modulus_divide (Number y) 594 { 595 if (!is_integer () || !y.is_integer ()) 596 { 597 /* Translators: Error displayed when attemping to do a modulus division on non-integer numbers */ 598 error = _("Modulus division is only defined for integers"); 599 return new Number.integer (0); 600 } 601 602 var t1 = divide (y).floor (); 603 var t2 = t1.multiply (y); 604 var z = subtract (t2); 605 606 t1 = new Number.integer (0); 607 if ((y.compare (t1) < 0 && z.compare (t1) > 0) || (y.compare (t1) > 0 && z.compare (t1) < 0)) 608 z = z.add (y); 609 610 return z; 611 } 612 613 /* Sets z = x ^ y mod p */ 614 public Number modular_exponentiation (Number exp, Number mod) 615 { 616 var base_value = copy (); 617 if (exp.is_negative ()) 618 base_value = base_value.reciprocal (); 619 var exp_value = exp.abs (); 620 var ans = new Number.integer (1); 621 var two = new Number.integer (2); 622 while (!exp_value.is_zero ()) 623 { 624 bool is_even = exp_value.modulus_divide (two).is_zero (); 625 if (!is_even) 626 { 627 ans = ans.multiply (base_value); 628 ans = ans.modulus_divide (mod); 629 } 630 base_value = base_value.multiply (base_value); 631 base_value = base_value.modulus_divide (mod); 632 exp_value = exp_value.divide_integer (2).floor (); 633 } 634 return ans.modulus_divide (mod); 635 } 636 637 /* Sets z = sin x */ 638 public Number sin (AngleUnit unit = AngleUnit.RADIANS) 639 { 640 var z = new Number (); 641 if (is_complex ()) 642 z.num.@set (num); 643 else 644 mpc_to_radians (z.num, num, unit); 645 z.num.sin (z.num); 646 return z; 647 } 648 649 /* Sets z = cos x */ 650 public Number cos (AngleUnit unit = AngleUnit.RADIANS) 651 { 652 var z = new Number (); 653 if (is_complex ()) 654 z.num.@set (num); 655 else 656 mpc_to_radians (z.num, num, unit); 657 z.num.cos (z.num); 658 return z; 659 } 660 661 /* Sets z = tan x */ 662 public Number tan (AngleUnit unit = AngleUnit.RADIANS) 663 { 664 /* Check for undefined values */ 665 var x_radians = to_radians (unit); 666 var check = x_radians.subtract (new Number.pi ().divide_integer (2)).divide (new Number.pi ()); 667 668 if (check.is_integer ()) 669 { 670 /* Translators: Error displayed when tangent value is undefined */ 671 error = _("Tangent is undefined for angles that are multiples of π (180°) from π∕2 (90°)"); 672 return new Number.integer (0); 673 } 674 675 var z = new Number (); 676 if (is_complex ()) 677 z.num.@set (num); 678 else 679 mpc_to_radians (z.num, num, unit); 680 z.num.tan (z.num); 681 return z; 682 } 683 684 /* Sets z = sin⁻¹ x */ 685 public Number asin (AngleUnit unit = AngleUnit.RADIANS) 686 { 687 if (compare (new Number.integer (1)) > 0 || compare (new Number.integer (-1)) < 0) 688 { 689 /* Translators: Error displayed when inverse sine value is undefined */ 690 error = _("Inverse sine is undefined for values outside [-1, 1]"); 691 return new Number.integer (0); 692 } 693 694 var z = new Number (); 695 z.num.asin (num); 696 if (!z.is_complex ()) 697 mpc_from_radians (z.num, z.num, unit); 698 return z; 699 } 700 701 /* Sets z = cos⁻¹ x */ 702 public Number acos (AngleUnit unit = AngleUnit.RADIANS) 703 { 704 if (compare (new Number.integer (1)) > 0 || compare (new Number.integer (-1)) < 0) 705 { 706 /* Translators: Error displayed when inverse cosine value is undefined */ 707 error = _("Inverse cosine is undefined for values outside [-1, 1]"); 708 return new Number.integer (0); 709 } 710 711 var z = new Number (); 712 z.num.acos (num); 713 if (!z.is_complex ()) 714 mpc_from_radians (z.num, z.num, unit); 715 return z; 716 } 717 718 /* Sets z = tan⁻¹ x */ 719 public Number atan (AngleUnit unit = AngleUnit.RADIANS) 720 { 721 /* Check x != i and x != -i */ 722 if (equals (new Number.integer (0,1)) || equals (new Number.integer(0,-1))) 723 { 724 /* Translators: Error displayed when trying to calculate undefined atan(i) or atan (-i) */ 725 error = _("Arctangent function is undefined for values i and -i"); 726 return new Number.integer (0); 727 } 728 var z = new Number (); 729 z.num.atan (num); 730 if (!z.is_complex ()) 731 mpc_from_radians (z.num, z.num, unit); 732 return z; 733 } 734 735 /* Sets z = sinh x */ 736 public Number sinh () 737 { 738 var z = new Number (); 739 z.num.sinh (num); 740 return z; 741 } 742 743 /* Sets z = cosh x */ 744 public Number cosh () 745 { 746 var z = new Number (); 747 z.num.cosh (num); 748 return z; 749 } 750 751 /* Sets z = tanh x */ 752 public Number tanh () 753 { 754 var z = new Number (); 755 z.num.tanh (num); 756 return z; 757 } 758 759 /* Sets z = sinh⁻¹ x */ 760 public Number asinh () 761 { 762 var z = new Number (); 763 z.num.asinh (num); 764 return z; 765 } 766 767 /* Sets z = cosh⁻¹ x */ 768 public Number acosh () 769 { 770 /* Check x >= 1 */ 771 var t = new Number.integer (1); 772 if (compare (t) < 0) 773 { 774 /* Translators: Error displayed when inverse hyperbolic cosine value is undefined */ 775 error = _("Inverse hyperbolic cosine is undefined for values less than one"); 776 return new Number.integer (0); 777 } 778 779 var z = new Number (); 780 z.num.acosh (num); 781 return z; 782 } 783 784 /* Sets z = tanh⁻¹ x */ 785 public Number atanh () 786 { 787 /* Check -1 <= x <= 1 */ 788 if (compare (new Number.integer (1)) >= 0 || compare (new Number.integer (-1)) <= 0) 789 { 790 /* Translators: Error displayed when inverse hyperbolic tangent value is undefined */ 791 error = _("Inverse hyperbolic tangent is undefined for values outside [-1, 1]"); 792 return new Number.integer (0); 793 } 794 795 var z = new Number (); 796 z.num.atanh (num); 797 return z; 798 } 799 800 /* Sets z = boolean AND for each bit in x and z */ 801 public Number and (Number y) 802 { 803 if (! 804 is_positive_integer () || !y.is_positive_integer ()) 805 { 806 /* Translators: Error displayed when boolean AND attempted on non-integer values */ 807 error = _("Boolean AND is only defined for positive integers"); 808 } 809 810 return bitwise (y, (v1, v2) => { return v1 & v2; }, 0); 811 } 812 813 /* Sets z = boolean OR for each bit in x and z */ 814 public Number or (Number y) 815 { 816 if (!is_positive_integer () || !y.is_positive_integer ()) 817 { 818 /* Translators: Error displayed when boolean OR attempted on non-integer values */ 819 error = _("Boolean OR is only defined for positive integers"); 820 } 821 822 return bitwise (y, (v1, v2) => { return v1 | v2; }, 0); 823 } 824 825 /* Sets z = boolean XOR for each bit in x and z */ 826 public Number xor (Number y) 827 { 828 if (!is_positive_integer () || !y.is_positive_integer ()) 829 { 830 /* Translators: Error displayed when boolean XOR attempted on non-integer values */ 831 error = _("Boolean XOR is only defined for positive integers"); 832 } 833 834 return bitwise (y, (v1, v2) => { return v1 ^ v2; }, 0); 835 } 836 837 /* Sets z = boolean NOT for each bit in x and z for word of length 'wordlen' */ 838 public Number not (int wordlen) 839 { 840 if (!is_positive_integer ()) 841 { 842 /* Translators: Error displayed when boolean XOR attempted on non-integer values */ 843 error = _("Boolean NOT is only defined for positive integers"); 844 } 845 846 return bitwise (new Number.integer (0), (v1, v2) => { return v1 ^ 0xF; }, wordlen); 847 } 848 849 /* Sets z = x masked to 'wordlen' bits */ 850 public Number mask (Number x, int wordlen) 851 { 852 /* Convert to a hexadecimal string and use last characters */ 853 var text = x.to_hex_string (); 854 var len = text.length; 855 var offset = wordlen / 4; 856 offset = len > offset ? (int) len - offset: 0; 857 return mp_set_from_string (text.substring (offset), 16); 858 } 859 860 /* Sets z = x shifted by 'count' bits. Positive shift increases the value, negative decreases */ 861 public Number shift (int64 count) 862 { 863 if (!is_integer ()) 864 { 865 /* Translators: Error displayed when bit shift attempted on non-integer values */ 866 error = _("Shift is only possible on integer values"); 867 return new Number.integer (0); 868 } 869 870 if (count >= 0) 871 { 872 var multiplier = 1; 873 for (var i = 0; i < count; i++) 874 multiplier *= 2; 875 return multiply_integer (multiplier); 876 } 877 else 878 { 879 var multiplier = 1; 880 for (var i = 0; i < -count; i++) 881 multiplier *= 2; 882 return divide_integer (multiplier).floor (); 883 } 884 } 885 886 /* Sets z to be the ones complement of x for word of length 'wordlen' */ 887 public Number ones_complement (int wordlen) 888 { 889 return bitwise (new Number.integer (0), (v1, v2) => { return v1 ^ v2; }, wordlen).not (wordlen); 890 } 891 892 /* Sets z to be the twos complement of x for word of length 'wordlen' */ 893 public Number twos_complement (int wordlen) 894 { 895 return ones_complement (wordlen).add (new Number.integer (1)); 896 } 897 898 /* Returns a list of all prime factors in x as Numbers */ 899 public List<Number?> factorize () 900 { 901 var factors = new List<Number?> (); 902 903 var value = abs (); 904 905 if (value.is_zero ()) 906 { 907 factors.append (value); 908 return factors; 909 } 910 911 if (value.equals (new Number.integer (1))) 912 { 913 factors.append (this); 914 return factors; 915 } 916 917 // if value < 2^64-1, call for factorize_uint64 function which deals in integers 918 919 uint64 num = 1; 920 num = num << 63; 921 num += (num - 1); 922 var int_max = new Number.unsigned_integer (num); 923 924 if (value.compare (int_max) <= 0) 925 { 926 var factors_int64 = factorize_uint64 (value.to_unsigned_integer ()); 927 if (is_negative ()) 928 factors_int64.data = factors_int64.data.invert_sign (); 929 return factors_int64; 930 } 931 932 var divisor = new Number.integer (2); 933 while (true) 934 { 935 var tmp = value.divide (divisor); 936 if (tmp.is_integer ()) 937 { 938 value = tmp; 939 factors.append (divisor); 940 } 941 else 942 break; 943 } 944 945 divisor = new Number.integer (3); 946 var root = value.sqrt (); 947 while (divisor.compare (root) <= 0) 948 { 949 var tmp = value.divide (divisor); 950 if (tmp.is_integer ()) 951 { 952 value = tmp; 953 root = value.sqrt (); 954 factors.append (divisor); 955 } 956 else 957 { 958 tmp = divisor.add (new Number.integer (2)); 959 divisor = tmp; 960 } 961 } 962 963 if (value.compare (new Number.integer (1)) > 0) 964 factors.append (value); 965 966 if (is_negative ()) 967 factors.data = factors.data.invert_sign (); 968 969 return factors; 970 } 971 972 public List<Number?> factorize_uint64 (uint64 n) 973 { 974 var factors = new List<Number?> (); 975 while (n % 2 == 0) 976 { 977 n /= 2; 978 factors.append (new Number.unsigned_integer (2)); 979 } 980 981 for (uint64 divisor = 3; divisor <= n / divisor; divisor += 2) 982 { 983 while (n % divisor == 0) 984 { 985 n /= divisor; 986 factors.append (new Number.unsigned_integer (divisor)); 987 } 988 } 989 990 if (n > 1) 991 factors.append (new Number.unsigned_integer (n)); 992 return factors; 993 } 994 995 private Number copy () 996 { 997 var z = new Number (); 998 z.num.@set (num); 999 return z; 1000 } 1001 1002 private static void mpc_from_radians (Complex res, Complex op, AngleUnit unit) 1003 { 1004 int i; 1005 1006 switch (unit) 1007 { 1008 default: 1009 case AngleUnit.RADIANS: 1010 if (res != op) 1011 res.@set (op); 1012 return; 1013 1014 case AngleUnit.DEGREES: 1015 i = 180; 1016 break; 1017 1018 case AngleUnit.GRADIANS: 1019 i=200; 1020 break; 1021 1022 } 1023 var scale = MPFR.Real (precision); 1024 scale.const_pi (); 1025 scale.signed_integer_divide (i, scale); 1026 res.multiply_mpreal (op, scale); 1027 } 1028 1029 private static void mpc_to_radians (Complex res, Complex op, AngleUnit unit) 1030 { 1031 int i; 1032 1033 switch (unit) 1034 { 1035 default: 1036 case AngleUnit.RADIANS: 1037 if (res != op) 1038 res.@set (op); 1039 return; 1040 1041 case AngleUnit.DEGREES: 1042 i = 180; 1043 break; 1044 1045 case AngleUnit.GRADIANS: 1046 i=200; 1047 break; 1048 } 1049 var scale = MPFR.Real (precision); 1050 scale.const_pi (); 1051 scale.divide_signed_integer (scale, i); 1052 res.multiply_mpreal (op, scale); 1053 } 1054 1055 /* Convert x to radians */ 1056 private Number to_radians (AngleUnit unit) 1057 { 1058 var z = new Number (); 1059 mpc_to_radians (z.num, num, unit); 1060 return z; 1061 } 1062 1063 private Number bitwise (Number y, BitwiseFunc bitwise_operator, int wordlen) 1064 { 1065 var text1 = to_hex_string (); 1066 var text2 = y.to_hex_string (); 1067 var offset1 = text1.length - 1; 1068 var offset2 = text2.length - 1; 1069 var offset_out = wordlen / 4 - 1; 1070 if (offset_out <= 0) 1071 offset_out = offset1 > offset2 ? offset1 : offset2; 1072 if (offset_out > 0 && (offset_out < offset1 || offset_out < offset2)) 1073 { 1074 error = ("Overflow. Try a bigger word size"); 1075 return new Number.integer (0); 1076 } 1077 1078 var text_out = new char[offset_out + 2]; 1079 1080 /* Perform bitwise operator on each character from right to left */ 1081 for (text_out[offset_out+1] = '\0'; offset_out >= 0; offset_out--) 1082 { 1083 int v1 = 0, v2 = 0; 1084 const char digits[] = { '0', '1', '2', '3', '4', '5', '6', '7', '8', '9', 'A', 'B', 'C', 'D', 'E', 'F' }; 1085 1086 if (offset1 >= 0) 1087 { 1088 v1 = hex_to_int (text1[offset1]); 1089 offset1--; 1090 } 1091 if (offset2 >= 0) 1092 { 1093 v2 = hex_to_int (text2[offset2]); 1094 offset2--; 1095 } 1096 text_out[offset_out] = digits[bitwise_operator (v1, v2)]; 1097 } 1098 1099 return mp_set_from_string ((string) text_out, 16); 1100 } 1101 1102 private int hex_to_int (char digit) 1103 { 1104 if (digit >= '0' && digit <= '9') 1105 return digit - '0'; 1106 if (digit >= 'A' && digit <= 'F') 1107 return digit - 'A' + 10; 1108 if (digit >= 'a' && digit <= 'f') 1109 return digit - 'a' + 10; 1110 return 0; 1111 } 1112 1113 private string to_hex_string () 1114 { 1115 var serializer = new Serializer (DisplayFormat.FIXED, 16, 0); 1116 return serializer.to_string (this); 1117 } 1118} 1119 1120private static int parse_literal_prefix (string str, ref int prefix_len) 1121{ 1122 var new_base = 0; 1123 1124 if (str.length < 3 || str[0] != '0') 1125 return new_base; 1126 1127 var prefix = str[1].tolower (); 1128 1129 if (prefix == 'b') 1130 new_base = 2; 1131 else if (prefix == 'o') 1132 new_base = 8; 1133 else if (prefix == 'x') 1134 new_base = 16; 1135 1136 if (new_base != 0) 1137 prefix_len = 2; 1138 1139 return new_base; 1140} 1141 1142// FIXME: Should all be in the class 1143 1144// FIXME: Re-add overflow and underflow detection 1145 1146/* Sets z from a string representation in 'text'. */ 1147public Number? mp_set_from_string (string str, int default_base = 10) 1148{ 1149 if (str.index_of_char ('°') >= 0) 1150 return set_from_sexagesimal (str); 1151 1152 /* Find the base */ 1153 const unichar base_digits[] = {'₀', '₁', '₂', '₃', '₄', '₅', '₆', '₇', '₈', '₉'}; 1154 var index = 0; 1155 var base_prefix = 0; 1156 unichar c; 1157 while (str.get_next_char (ref index, out c)); 1158 var end = index; 1159 var number_base = 0; 1160 var literal_base = 0; 1161 var base_multiplier = 1; 1162 while (str.get_prev_char (ref index, out c)) 1163 { 1164 var value = -1; 1165 for (var i = 0; i < base_digits.length; i++) 1166 { 1167 if (c == base_digits[i]) 1168 { 1169 value = i; 1170 break; 1171 } 1172 } 1173 if (value < 0) 1174 break; 1175 1176 end = index; 1177 number_base += value * base_multiplier; 1178 base_multiplier *= 10; 1179 } 1180 1181 literal_base = parse_literal_prefix (str, ref base_prefix); 1182 1183 if (number_base != 0 && literal_base != 0 && literal_base != number_base) 1184 return null; 1185 1186 if (number_base == 0) 1187 number_base = (literal_base != 0) ? literal_base : default_base; 1188 1189 /* Check if this has a sign */ 1190 var negate = false; 1191 index = base_prefix; 1192 str.get_next_char (ref index, out c); 1193 if (c == '+') 1194 negate = false; 1195 else if (c == '-' || c == '−') 1196 negate = true; 1197 else 1198 str.get_prev_char (ref index, out c); 1199 1200 /* Convert integer part */ 1201 var z = new Number.integer (0); 1202 1203 while (str.get_next_char (ref index, out c)) 1204 { 1205 var i = char_val (c, number_base); 1206 if (i > number_base) 1207 return null; 1208 if (i < 0) 1209 { 1210 str.get_prev_char (ref index, out c); 1211 break; 1212 } 1213 1214 z = z.multiply_integer (number_base).add (new Number.integer (i)); 1215 } 1216 1217 /* Look for fraction characters, e.g. ⅚ */ 1218 const unichar fractions[] = {'½', '⅓', '⅔', '¼', '¾', '⅕', '⅖', '⅗', '⅘', '⅙', '⅚', '⅛', '⅜', '⅝', '⅞'}; 1219 const int numerators[] = { 1, 1, 2, 1, 3, 1, 2, 3, 4, 1, 5, 1, 3, 5, 7}; 1220 const int denominators[] = { 2, 3, 3, 4, 4, 5, 5, 5, 5, 6, 6, 8, 8, 8, 8}; 1221 var has_fraction = false; 1222 if (str.get_next_char (ref index, out c)) 1223 { 1224 for (var i = 0; i < fractions.length; i++) 1225 { 1226 if (c == fractions[i]) 1227 { 1228 var fraction = new Number.fraction (numerators[i], denominators[i]); 1229 z = z.add (fraction); 1230 1231 /* Must end with fraction */ 1232 if (!str.get_next_char (ref index, out c)) 1233 return z; 1234 else 1235 return null; 1236 } 1237 } 1238 1239 /* Check for decimal point */ 1240 if (c == '.') 1241 has_fraction = true; 1242 else 1243 str.get_prev_char (ref index, out c); 1244 } 1245 1246 /* Convert fractional part */ 1247 if (has_fraction) 1248 { 1249 var numerator = new Number.integer (0); 1250 var denominator = new Number.integer (1); 1251 1252 while (str.get_next_char (ref index, out c)) 1253 { 1254 var i = char_val (c, number_base); 1255 if (i < 0) 1256 { 1257 str.get_prev_char (ref index, out c); 1258 break; 1259 } 1260 1261 denominator = denominator.multiply_integer (number_base); 1262 numerator = numerator.multiply_integer (number_base); 1263 numerator = numerator.add (new Number.integer (i)); 1264 } 1265 1266 numerator = numerator.divide (denominator); 1267 z = z.add (numerator); 1268 } 1269 1270 if (index != end) 1271 return null; 1272 1273 if (negate) 1274 z = z.invert_sign (); 1275 1276 return z; 1277} 1278 1279private int char_val (unichar c, int number_base) 1280{ 1281 if (!c.isxdigit ()) 1282 return -1; 1283 1284 var value = c.xdigit_value (); 1285 1286 if (value >= number_base) 1287 return -1; 1288 1289 return value; 1290} 1291 1292private Number? set_from_sexagesimal (string str) 1293{ 1294 var degree_index = str.index_of_char ('°'); 1295 if (degree_index < 0) 1296 return null; 1297 var degrees = mp_set_from_string (str.substring (0, degree_index)); 1298 if (degrees == null) 1299 return null; 1300 var minute_start = degree_index; 1301 unichar c; 1302 str.get_next_char (ref minute_start, out c); 1303 1304 if (str[minute_start] == '\0') 1305 return degrees; 1306 var minute_index = str.index_of_char ('\'', minute_start); 1307 if (minute_index < 0) 1308 return null; 1309 var minutes = mp_set_from_string (str.substring (minute_start, minute_index - minute_start)); 1310 if (minutes == null) 1311 return null; 1312 degrees = degrees.add (minutes.divide_integer (60)); 1313 var second_start = minute_index; 1314 str.get_next_char (ref second_start, out c); 1315 1316 if (str[second_start] == '\0') 1317 return degrees; 1318 var second_index = str.index_of_char ('"', second_start); 1319 if (second_index < 0) 1320 return null; 1321 var seconds = mp_set_from_string (str.substring (second_start, second_index - second_start)); 1322 if (seconds == null) 1323 return null; 1324 degrees = degrees.add (seconds.divide_integer (3600)); 1325 str.get_next_char (ref second_index, out c); 1326 1327 /* Skip over second marker and expect no more characters */ 1328 if (str[second_index] == '\0') 1329 return degrees; 1330 else 1331 return null; 1332} 1333 1334/* Returns true if x is cannot be represented in a binary word of length 'wordlen' */ 1335public bool mp_is_overflow (Number x, int wordlen) 1336{ 1337 var t2 = new Number.integer (2).xpowy_integer (wordlen); 1338 return t2.compare (x) > 0; 1339} 1340