diff --git a/app/src/main/java/midpcalc/Real.java b/app/src/main/java/midpcalc/Real.java
new file mode 100755
index 00000000..d6a6e0d0
--- /dev/null
+++ b/app/src/main/java/midpcalc/Real.java
@@ -0,0 +1,6021 @@
+
+
+
+package midpcalc;
+
+
+
+/**
+ * Java integer implementation of 63-bit precision floating point.
+ *
Version 1.13
+ *
+ * Copyright 2003-2009 Roar Lauritzsen
+ *
+ *
+ *
+ * This library is free software; you can redistribute it and/or modify it
+ * under the terms of the GNU General Public License as published by the Free
+ * Software Foundation; either version 2 of the License, or (at your option)
+ * any later version.
+ *
+ *
This library is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+ * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ * for more details.
+ *
+ *
The following link provides a copy of the GNU General Public License:
+ *
http://www.gnu.org/licenses/gpl.txt
+ *
If you are unable to obtain the copy from this address, write to the
+ * Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
+ * 02111-1307 USA
+ *
+ *
+ *
+ * General notes
+ *
+ *
+ * Real objects are not immutable, like Java
+ * Double or BigDecimal. This means that you
+ * should not think of a Real object as a "number", but rather
+ * as a "register holding a number". This design choice is done to encourage
+ * object reuse and limit garbage production for more efficient execution on
+ * e.g. a limited MIDP device. The design choice is reflected in the API,
+ * where an operation like {@link #add(Real) add} does not return a new
+ * object containing the result (as with {@link
+ * java.math.BigDecimal#add(java.math.BigDecimal) BigDecimal}), but rather
+ * adds the argument to the object itself, and returns nothing.
+ *
+ * - This library implements infinities and NaN (Not-a-Number) following
+ * the IEEE 754 logic. If an operation produces a result larger (in
+ * magnitude) than the largest representable number, a value representing
+ * positive or negative infinity is generated. If an operation produces a
+ * result smaller than the smallest representable number, a positive or
+ * negative zero is generated. If an operation is undefined, a NaN value is
+ * produced. Abnormal numbers are often fine to use in further
+ * calculations. In most cases where the final result would be meaningful,
+ * abnormal numbers accomplish this, e.g. atan(1/0)=π/2. In most cases
+ * where the final result is not meaningful, a NaN will be produced.
+ * No exception is ever (deliberately) thrown.
+ *
+ *
- Error bounds listed under Method Detail
+ * are calculated using William Rossi's rossi.dfp.dfp at 40 decimal digits
+ * accuracy. Error bounds are for "typical arguments" and may increase when
+ * results approach zero or
+ * infinity. The abbreviation {@link Math#ulp(double) ULP} means Unit in the
+ * Last Place. An error bound of ½ ULP means that the result is correctly
+ * rounded. The relative execution time listed under each method is the
+ * average from running on SonyEricsson T610 (R3C), K700i, and Nokia 6230i.
+ *
+ *
- The library is not thread-safe. Static
Real objects are
+ * used extensively as temporary values to avoid garbage production and the
+ * overhead of new. To make the library thread-safe, references
+ * to all these static objects must be replaced with code that instead
+ * allocates new Real objects in their place.
+ *
+ * - There is one bug that occurs again and again and is really difficult
+ * to debug. Although the pre-calculated constants are declared
static
+ * final, Java cannot really protect the contents of the objects in
+ * the same way as consts are protected in C/C++. Consequently,
+ * you can accidentally change these values if you send them into a function
+ * that modifies its arguments. If you were to modify {@link #ONE Real.ONE}
+ * for instance, many of the succeeding calculations would be wrong because
+ * the same variable is used extensively in the internal calculations of
+ * Real.java.
+ *
+ *
+ */
+public final class Real
+{
+ /**
+ * The mantissa of a Real. To maintain numbers in a
+ * normalized state and to preserve the integrity of abnormal numbers, it
+ * is discouraged to modify the inner representation of a
+ * Real directly.
+ *
+ * The number represented by a Real equals:
+ * -1sign · mantissa · 2-62 · 2exponent-0x40000000
+ *
+ *
The normalized mantissa of a finite Real must be
+ * between 0x4000000000000000L and
+ * 0x7fffffffffffffffL. Using a denormalized
+ * Real in any operation other than {@link
+ * #normalize()} may produce undefined results. The mantissa of zero and
+ * of an infinite value is 0x0000000000000000L.
+ *
+ *
The mantissa of a NaN is any nonzero value. However, it is
+ * recommended to use the value 0x4000000000000000L. Any
+ * other values are reserved for future extensions.
+ */
+ public long mantissa;
+ /**
+ * The exponent of a Real. To maintain numbers in a
+ * normalized state and to preserve the integrity of abnormal numbers, it
+ * is discouraged to modify the inner representation of a
+ * Real directly.
+ *
+ *
The exponent of a finite Real must be between
+ * 0x00000000 and 0x7fffffff. The exponent of
+ * zero 0x00000000.
+ *
+ *
The exponent of an infinite value and of a NaN is any negative
+ * value. However, it is recommended to use the value
+ * 0x80000000. Any other values are reserved for future
+ * extensions.
+ */
+ public int exponent;
+ /**
+ * The sign of a Real. To maintain numbers in a normalized
+ * state and to preserve the integrity of abnormal numbers, it is
+ * discouraged to modify the inner representation of a Real
+ * directly.
+ *
+ *
The sign of a finite, zero or infinite Real is 0 for
+ * positive values and 1 for negative values. Any other values may produce
+ * undefined results.
+ *
+ *
The sign of a NaN is ignored. However, it is recommended to use the
+ * value 0. Any other values are reserved for future
+ * extensions.
+ */
+ public byte sign;
+ /**
+ * Set to false during numerical algorithms to favor accuracy
+ * over prettyness. This flag is initially set to true.
+ *
+ *
The flag controls the operation of a subtraction of two
+ * almost-identical numbers that differ only in the last three bits of the
+ * mantissa. With this flag enabled, the result of such a subtraction is
+ * rounded down to zero. Probabilistically, this is the correct course of
+ * action in an overwhelmingly large percentage of calculations.
+ * However, certain numerical algorithms such as differentiation depend
+ * on keeping maximum accuracy during subtraction.
+ *
+ *
Note, that because of magicRounding,
+ * a.sub(b) may produce zero even though
+ * a.equalTo(b) returns false. This must be
+ * considered e.g. when trying to avoid division by zero.
+ */
+ public static boolean magicRounding = true;
+ /**
+ * A Real constant holding the exact value of 0. Among other
+ * uses, this value is used as a result when a positive underflow occurs.
+ */
+ public static final Real ZERO = new Real(0,0x00000000,0x0000000000000000L);
+ /**
+ * A Real constant holding the exact value of 1.
+ */
+ public static final Real ONE = new Real(0,0x40000000,0x4000000000000000L);
+ /**
+ * A Real constant holding the exact value of 2.
+ */
+ public static final Real TWO = new Real(0,0x40000001,0x4000000000000000L);
+ /**
+ * A Real constant holding the exact value of 3.
+ */
+ public static final Real THREE= new Real(0,0x40000001,0x6000000000000000L);
+ /**
+ * A Real constant holding the exact value of 5.
+ */
+ public static final Real FIVE = new Real(0,0x40000002,0x5000000000000000L);
+ /**
+ * A Real constant holding the exact value of 10.
+ */
+ public static final Real TEN = new Real(0,0x40000003,0x5000000000000000L);
+ /**
+ * A Real constant holding the exact value of 100.
+ */
+ public static final Real HUNDRED=new Real(0,0x40000006,0x6400000000000000L);
+ /**
+ * A Real constant holding the exact value of 1/2.
+ */
+ public static final Real HALF = new Real(0,0x3fffffff,0x4000000000000000L);
+ /**
+ * A Real constant that is closer than any other to 1/3.
+ */
+ public static final Real THIRD= new Real(0,0x3ffffffe,0x5555555555555555L);
+ /**
+ * A Real constant that is closer than any other to 1/10.
+ */
+ public static final Real TENTH= new Real(0,0x3ffffffc,0x6666666666666666L);
+ /**
+ * A Real constant that is closer than any other to 1/100.
+ */
+ public static final Real PERCENT=new Real(0,0x3ffffff9,0x51eb851eb851eb85L);
+ /**
+ * A Real constant that is closer than any other to the
+ * square root of 2.
+ */
+ public static final Real SQRT2= new Real(0,0x40000000,0x5a827999fcef3242L);
+ /**
+ * A Real constant that is closer than any other to the
+ * square root of 1/2.
+ */
+ public static final Real SQRT1_2=new Real(0,0x3fffffff,0x5a827999fcef3242L);
+ /**
+ * A Real constant that is closer than any other to 2π.
+ */
+ public static final Real PI2 = new Real(0,0x40000002,0x6487ed5110b4611aL);
+ /**
+ * A Real constant that is closer than any other to π, the
+ * ratio of the circumference of a circle to its diameter.
+ */
+ public static final Real PI = new Real(0,0x40000001,0x6487ed5110b4611aL);
+ /**
+ * A Real constant that is closer than any other to π/2.
+ */
+ public static final Real PI_2 = new Real(0,0x40000000,0x6487ed5110b4611aL);
+ /**
+ * A Real constant that is closer than any other to π/4.
+ */
+ public static final Real PI_4 = new Real(0,0x3fffffff,0x6487ed5110b4611aL);
+ /**
+ * A Real constant that is closer than any other to π/8.
+ */
+ public static final Real PI_8 = new Real(0,0x3ffffffe,0x6487ed5110b4611aL);
+ /**
+ * A Real constant that is closer than any other to e,
+ * the base of the natural logarithms.
+ */
+ public static final Real E = new Real(0,0x40000001,0x56fc2a2c515da54dL);
+ /**
+ * A Real constant that is closer than any other to the
+ * natural logarithm of 2.
+ */
+ public static final Real LN2 = new Real(0,0x3fffffff,0x58b90bfbe8e7bcd6L);
+ /**
+ * A Real constant that is closer than any other to the
+ * natural logarithm of 10.
+ */
+ public static final Real LN10 = new Real(0,0x40000001,0x49aec6eed554560bL);
+ /**
+ * A Real constant that is closer than any other to the
+ * base-2 logarithm of e.
+ */
+ public static final Real LOG2E= new Real(0,0x40000000,0x5c551d94ae0bf85eL);
+ /**
+ * A Real constant that is closer than any other to the
+ * base-10 logarithm of e.
+ */
+ public static final Real LOG10E=new Real(0,0x3ffffffe,0x6f2dec549b9438cbL);
+ /**
+ * A Real constant holding the maximum non-infinite positive
+ * number = 4.197e323228496.
+ */
+ public static final Real MAX = new Real(0,0x7fffffff,0x7fffffffffffffffL);
+ /**
+ * A Real constant holding the minimum non-zero positive
+ * number = 2.383e-323228497.
+ */
+ public static final Real MIN = new Real(0,0x00000000,0x4000000000000000L);
+ /**
+ * A Real constant holding the value of NaN (not-a-number).
+ * This value is always used as a result to signal an invalid operation.
+ */
+ public static final Real NAN = new Real(0,0x80000000,0x4000000000000000L);
+ /**
+ * A Real constant holding the value of positive infinity.
+ * This value is always used as a result to signal a positive overflow.
+ */
+ public static final Real INF = new Real(0,0x80000000,0x0000000000000000L);
+ /**
+ * A Real constant holding the value of negative infinity.
+ * This value is always used as a result to signal a negative overflow.
+ */
+ public static final Real INF_N= new Real(1,0x80000000,0x0000000000000000L);
+ /**
+ * A Real constant holding the value of negative zero. This
+ * value is used as a result e.g. when a negative underflow occurs.
+ */
+ public static final Real ZERO_N=new Real(1,0x00000000,0x0000000000000000L);
+ /**
+ * A Real constant holding the exact value of -1.
+ */
+ public static final Real ONE_N= new Real(1,0x40000000,0x4000000000000000L);
+ private static final int clz_magic = 0x7c4acdd;
+ private static final byte[] clz_tab =
+ { 31,22,30,21,18,10,29, 2,20,17,15,13, 9, 6,28, 1,
+ 23,19,11, 3,16,14, 7,24,12, 4, 8,25, 5,26,27, 0 };
+ /**
+ * Creates a new Real with a value of zero.
+ */
+ public Real() {
+ }
+ /**
+ * Creates a new Real, assigning the value of another
+ * Real. See {@link #assign(Real)}.
+ *
+ * @param a the Real to assign.
+ */
+ public Real(Real a) {
+ { this.mantissa = a.mantissa; this.exponent = a.exponent; this.sign = a.sign; };
+ }
+ /**
+ * Creates a new Real, assigning the value of an integer. See
+ * {@link #assign(int)}.
+ *
+ * @param a the int to assign.
+ */
+ public Real(int a) {
+ assign(a);
+ }
+ /**
+ * Creates a new Real, assigning the value of a long
+ * integer. See {@link #assign(long)}.
+ *
+ * @param a the long to assign.
+ */
+ public Real(long a) {
+ assign(a);
+ }
+ /**
+ * Creates a new Real, assigning the value encoded in a
+ * String using base-10. See {@link #assign(String)}.
+ *
+ * @param a the String to assign.
+ */
+ public Real(String a) {
+ assign(a,10);
+ }
+ /**
+ * Creates a new Real, assigning the value encoded in a
+ * String using the specified number base. See {@link
+ * #assign(String,int)}.
+ *
+ * @param a the String to assign.
+ * @param base the number base of a. Valid base values are 2,
+ * 8, 10 and 16.
+ */
+ public Real(String a, int base) {
+ assign(a,base);
+ }
+ /**
+ * Creates a new Real, assigning a value by directly setting
+ * the fields of the internal representation. The arguments must represent
+ * a valid, normalized Real. This is the fastest way of
+ * creating a constant value. See {@link #assign(int,int,long)}.
+ *
+ * @param s {@link #sign} bit, 0 for positive sign, 1 for negative sign
+ * @param e {@link #exponent}
+ * @param m {@link #mantissa}
+ */
+ public Real(int s, int e, long m) {
+ { this.sign=(byte)s; this.exponent=e; this.mantissa=m; };
+ }
+ /**
+ * Creates a new Real, assigning the value previously encoded
+ * into twelve consecutive bytes in a byte array using {@link
+ * #toBytes(byte[],int) toBytes}. See {@link #assign(byte[],int)}.
+ *
+ * @param data byte array to decode into this Real.
+ * @param offset offset to start encoding from. The bytes
+ * data[offset]...data[offset+11] will be
+ * read.
+ */
+ public Real(byte [] data, int offset) {
+ assign(data,offset);
+ }
+ /**
+ * Assigns this Real the value of another Real.
+ *
+ *
+ * Equivalent double code: |
+ * this = a;
+ * |
| Error bound: |
+ * 0 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 0.3
+ * |
+ *
+ * @param a the Real to assign.
+ */
+ public void assign(Real a) {
+ if (a == null) {
+ makeZero();
+ return;
+ }
+ sign = a.sign;
+ exponent = a.exponent;
+ mantissa = a.mantissa;
+ }
+ /**
+ * Assigns this Real the value of an integer.
+ * All integer values can be represented without loss of accuracy.
+ *
+ *
+ * Equivalent double code: |
+ * this = (double)a;
+ * |
| Error bound: |
+ * 0 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 0.6
+ * |
+ *
+ * @param a the int to assign.
+ */
+ public void assign(int a) {
+ if (a==0) {
+ makeZero();
+ return;
+ }
+ sign = 0;
+ if (a<0) {
+ sign = 1;
+ a = -a; // Also works for 0x80000000
+ }
+ // Normalize int
+ int t=a; t|=t>>1; t|=t>>2; t|=t>>4; t|=t>>8; t|=t>>16;
+ t = clz_tab[(t*clz_magic)>>>27]-1;
+ exponent = 0x4000001E-t;
+ mantissa = ((long)a)<<(32+t);
+ }
+ /**
+ * Assigns this Real the value of a signed long integer.
+ * All long values can be represented without loss of accuracy.
+ *
+ *
+ * Equivalent double code: |
+ * this = (double)a;
+ * |
| Error bound: |
+ * 0 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 1.0
+ * |
+ *
+ * @param a the long to assign.
+ */
+ public void assign(long a) {
+ sign = 0;
+ if (a<0) {
+ sign = 1;
+ a = -a; // Also works for 0x8000000000000000
+ }
+ exponent = 0x4000003E;
+ mantissa = a;
+ normalize();
+ }
+ /**
+ * Assigns this Real a value encoded in a String
+ * using base-10, as specified in {@link #assign(String,int)}.
+ *
+ *
+ * Equivalent double code: |
+ * this = Double.{@link Double#valueOf(String) valueOf}(a);
+ * |
| Approximate error bound: |
+ * ½-1 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 80
+ * |
+ *
+ * @param a the String to assign.
+ */
+ public void assign(String a) {
+ assign(a,10);
+ }
+ /**
+ * Assigns this Real a value encoded in a String
+ * using the specified number base. The string is parsed as follows:
+ *
+ *
+ * - If the string is
null or an empty string, zero is
+ * assigned.
+ * - Leading spaces are ignored.
+ *
- An optional sign, '+', '-' or '/', where '/' precedes a negative
+ * two's-complement number, reading: "an infinite number of 1-bits
+ * preceding the number".
+ *
- Optional digits preceding the radix, in the specified base.
+ *
+ * - In base-2, allowed digits are '01'.
+ *
- In base-8, allowed digits are '01234567'.
+ *
- In base-10, allowed digits are '0123456789'.
+ *
- In base-16, allowed digits are '0123456789ABCDEF'.
+
+ * - An optional radix character, '.' or ','.
+ *
- Optional digits following the radix.
+ *
- The following spaces are ignored.
+ *
- An optional exponent indicator, 'e'. If not base-16, or after a
+ * space, 'E' is also accepted.
+ *
- An optional sign, '+' or '-'.
+ *
- Optional exponent digits in base-10.
+ *
+ *
+ * Valid examples:
+ * base-2: "-.110010101e+5"
+ * base-8: "+5462E-99"
+ * base-10: " 3,1415927"
+ * base-16: "/FFF800C.CCCE e64"
+ *
+ *
The number is parsed until the end of the string or an unknown
+ * character is encountered, then silently returns even if the whole
+ * string has not been parsed. Please note that specifying an
+ * excessive number of digits in base-10 may in fact decrease the
+ * accuracy of the result because of the extra multiplications performed.
+ *
+ *
+ * Equivalent double code: |
+ *
+ * this = Double.{@link Double#valueOf(String) valueOf}(a);
+ * // Works only for base-10
+ * |
|
+ * Approximate error bound:
+ * | base-10 |
+ * ½-1 ULPs
+ * |
| 2/8/16 |
+ * ½ ULPs
+ * |
|
+ * Execution time relative to add:
+ * | base-2 |
+ * 54
+ * |
| base-8 |
+ * 60
+ * |
| base-10 |
+ * 80
+ * |
| base-16 |
+ * 60
+ * |
+ *
+ * @param a the String to assign.
+ * @param base the number base of a. Valid base values are
+ * 2, 8, 10 and 16.
+ */
+ public void assign(String a, int base) {
+ if (a==null || a.length()==0) {
+ assign(ZERO);
+ return;
+ }
+ atof(a,base);
+ }
+ /**
+ * Assigns this Real a value by directly setting the fields
+ * of the internal representation. The arguments must represent a valid,
+ * normalized Real. This is the fastest way of assigning a
+ * constant value.
+ *
+ *
+ * Equivalent double code: |
+ * this = (1-2*s) * m *
+ * Math.{@link Math#pow(double,double)
+ * pow}(2.0,e-0x400000e3);
+ * |
| Error bound: |
+ * 0 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 0.3
+ * |
+ *
+ * @param s {@link #sign} bit, 0 for positive sign, 1 for negative sign
+ * @param e {@link #exponent}
+ * @param m {@link #mantissa}
+ */
+ public void assign(int s, int e, long m) {
+ sign = (byte)s;
+ exponent = e;
+ mantissa = m;
+ }
+ /**
+ * Assigns this Real a value previously encoded into into
+ * twelve consecutive bytes in a byte array using {@link
+ * #toBytes(byte[],int) toBytes}.
+ *
+ * |
+ * Error bound: |
+ * 0 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 1.2
+ * |
+ *
+ * @param data byte array to decode into this Real.
+ * @param offset offset to start encoding from. The bytes
+ * data[offset]...data[offset+11] will be
+ * read.
+ */
+ public void assign(byte [] data, int offset) {
+ sign = (byte)((data[offset+4]>>7)&1);
+ exponent = (((data[offset ]&0xff)<<24)+
+ ((data[offset +1]&0xff)<<16)+
+ ((data[offset +2]&0xff)<<8)+
+ ((data[offset +3]&0xff)));
+ mantissa = (((long)(data[offset+ 4]&0x7f)<<56)+
+ ((long)(data[offset+ 5]&0xff)<<48)+
+ ((long)(data[offset+ 6]&0xff)<<40)+
+ ((long)(data[offset+ 7]&0xff)<<32)+
+ ((long)(data[offset+ 8]&0xff)<<24)+
+ ((long)(data[offset+ 9]&0xff)<<16)+
+ ((long)(data[offset+10]&0xff)<< 8)+
+ ( (data[offset+11]&0xff)));
+ }
+ /**
+ * Encodes an accurate representation of this Real value into
+ * twelve consecutive bytes in a byte array. Can be decoded using {@link
+ * #assign(byte[],int)}.
+ *
+ * |
+ * Execution time relative to add:
+ * |
+ * 1.2
+ * |
+ *
+ * @param data byte array to save this Real in.
+ * @param offset offset to start encoding to. The bytes
+ * data[offset]...data[offset+11] will be
+ * written.
+ */
+ public void toBytes(byte [] data, int offset) {
+ data[offset ] = (byte)(exponent>>24);
+ data[offset+ 1] = (byte)(exponent>>16);
+ data[offset+ 2] = (byte)(exponent>>8);
+ data[offset+ 3] = (byte)(exponent);
+ data[offset+ 4] = (byte)((sign<<7)+(mantissa>>56));
+ data[offset+ 5] = (byte)(mantissa>>48);
+ data[offset+ 6] = (byte)(mantissa>>40);
+ data[offset+ 7] = (byte)(mantissa>>32);
+ data[offset+ 8] = (byte)(mantissa>>24);
+ data[offset+ 9] = (byte)(mantissa>>16);
+ data[offset+10] = (byte)(mantissa>>8);
+ data[offset+11] = (byte)(mantissa);
+ }
+ /**
+ * Assigns this Real the value corresponding to a given bit
+ * representation. The argument is considered to be a representation of a
+ * floating-point value according to the IEEE 754 floating-point "single
+ * format" bit layout.
+ *
+ *
+ * Equivalent float code: |
+ * this = Float.{@link Float#intBitsToFloat(int)
+ * intBitsToFloat}(bits);
+ * |
| Error bound: |
+ * 0 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 0.6
+ * |
+ *
+ * @param bits a float value encoded in an int.
+ */
+ public void assignFloatBits(int bits) {
+ sign = (byte)(bits>>>31);
+ exponent = (bits>>23)&0xff;
+ mantissa = (long)(bits&0x007fffff)<<39;
+ if (exponent == 0 && mantissa == 0)
+ return; // Valid zero
+ if (exponent == 0 && mantissa != 0) {
+ // Degenerate small float
+ exponent = 0x40000000-126;
+ normalize();
+ return;
+ }
+ if (exponent <= 254) {
+ // Normal IEEE 754 float
+ exponent += 0x40000000-127;
+ mantissa |= 1L<<62;
+ return;
+ }
+ if (mantissa == 0)
+ makeInfinity(sign);
+ else
+ makeNan();
+ }
+ /**
+ * Assigns this Real the value corresponding to a given bit
+ * representation. The argument is considered to be a representation of a
+ * floating-point value according to the IEEE 754 floating-point "double
+ * format" bit layout.
+ *
+ *
+ * Equivalent double code: |
+ * this = Double.{@link Double#longBitsToDouble(long)
+ * longBitsToDouble}(bits);
+ * |
| Error bound: |
+ * 0 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 0.6
+ * |
+ *
+ * @param bits a double value encoded in a long.
+ */
+ public void assignDoubleBits(long bits) {
+ sign = (byte)((bits>>63)&1);
+ exponent = (int)((bits>>52)&0x7ff);
+ mantissa = (bits&0x000fffffffffffffL)<<10;
+ if (exponent == 0 && mantissa == 0)
+ return; // Valid zero
+ if (exponent == 0 && mantissa != 0) {
+ // Degenerate small float
+ exponent = 0x40000000-1022;
+ normalize();
+ return;
+ }
+ if (exponent <= 2046) {
+ // Normal IEEE 754 float
+ exponent += 0x40000000-1023;
+ mantissa |= 1L<<62;
+ return;
+ }
+ if (mantissa == 0)
+ makeInfinity(sign);
+ else
+ makeNan();
+ }
+ /**
+ * Returns a representation of this Real according to the
+ * IEEE 754 floating-point "single format" bit layout.
+ *
+ *
+ * Equivalent float code: |
+ * Float.{@link Float#floatToIntBits(float)
+ * floatToIntBits}(this)
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 0.7
+ * |
+ *
+ * @return the bits that represent the floating-point number.
+ */
+ public int toFloatBits() {
+ if ((this.exponent < 0 && this.mantissa != 0))
+ return 0x7fffffff; // nan
+ int e = exponent-0x40000000+127;
+ long m = mantissa;
+ // Round properly!
+ m += 1L<<38;
+ if (m<0) {
+ m >>>= 1;
+ e++;
+ if (exponent < 0) // Overflow
+ return (sign<<31)|0x7f800000; // inf
+ }
+ if ((this.exponent < 0 && this.mantissa == 0) || e > 254)
+ return (sign<<31)|0x7f800000; // inf
+ if ((this.exponent == 0 && this.mantissa == 0) || e < -22)
+ return (sign<<31); // zero
+ if (e <= 0) // Degenerate small float
+ return (sign<<31)|((int)(m>>>(40-e))&0x007fffff);
+ // Normal IEEE 754 float
+ return (sign<<31)|(e<<23)|((int)(m>>>39)&0x007fffff);
+ }
+ /**
+ * Returns a representation of this Real according to the
+ * IEEE 754 floating-point "double format" bit layout.
+ *
+ *
+ * Equivalent double code: |
+ * Double.{@link Double#doubleToLongBits(double)
+ * doubleToLongBits}(this)
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 0.7
+ * |
+ *
+ * @return the bits that represent the floating-point number.
+ */
+ public long toDoubleBits() {
+ if ((this.exponent < 0 && this.mantissa != 0))
+ return 0x7fffffffffffffffL; // nan
+ int e = exponent-0x40000000+1023;
+ long m = mantissa;
+ // Round properly!
+ m += 1L<<9;
+ if (m<0) {
+ m >>>= 1;
+ e++;
+ if (exponent < 0)
+ return ((long)sign<<63)|0x7ff0000000000000L; // inf
+ }
+ if ((this.exponent < 0 && this.mantissa == 0) || e > 2046)
+ return ((long)sign<<63)|0x7ff0000000000000L; // inf
+ if ((this.exponent == 0 && this.mantissa == 0) || e < -51)
+ return ((long)sign<<63); // zero
+ if (e <= 0) // Degenerate small double
+ return ((long)sign<<63)|((m>>>(11-e))&0x000fffffffffffffL);
+ // Normal IEEE 754 double
+ return ((long)sign<<63)|((long)e<<52)|((m>>>10)&0x000fffffffffffffL);
+ }
+ /**
+ * Makes this Real the value of positive zero.
+ *
+ *
+ * Equivalent double code: |
+ * this = 0;
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 0.2
+ * |
+ */
+ public void makeZero() {
+ sign = 0;
+ mantissa = 0;
+ exponent = 0;
+ }
+ /**
+ * Makes this Real the value of zero with the specified sign.
+ *
+ *
+ * Equivalent double code: |
+ * this = 0.0 * (1-2*s);
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 0.2
+ * |
+ *
+ * @param s sign bit, 0 to make a positive zero, 1 to make a negative zero
+ */
+ public void makeZero(int s) {
+ sign = (byte)s;
+ mantissa = 0;
+ exponent = 0;
+ }
+ /**
+ * Makes this Real the value of infinity with the specified
+ * sign.
+ *
+ *
+ * Equivalent double code: |
+ * this = Double.{@link Double#POSITIVE_INFINITY POSITIVE_INFINITY}
+ * * (1-2*s);
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 0.3
+ * |
+ *
+ * @param s sign bit, 0 to make positive infinity, 1 to make negative
+ * infinity
+ */
+ public void makeInfinity(int s) {
+ sign = (byte)s;
+ mantissa = 0;
+ exponent = 0x80000000;
+ }
+ /**
+ * Makes this Real the value of Not-a-Number (NaN).
+ *
+ *
+ * Equivalent double code: |
+ * this = Double.{@link Double#NaN NaN};
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 0.3
+ * |
+ */
+ public void makeNan() {
+ sign = 0;
+ mantissa = 0x4000000000000000L;
+ exponent = 0x80000000;
+ }
+ /**
+ * Returns true if the value of this Real is
+ * zero, false otherwise.
+ *
+ *
+ * Equivalent double code: |
+ * (this == 0)
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 0.3
+ * |
+ *
+ * @return true if the value represented by this object is
+ * zero, false otherwise.
+ */
+ public boolean isZero() {
+ return (exponent == 0 && mantissa == 0);
+ }
+ /**
+ * Returns true if the value of this Real is
+ * infinite, false otherwise.
+ *
+ *
+ * Equivalent double code: |
+ * Double.{@link Double#isInfinite(double) isInfinite}(this)
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 0.3
+ * |
+ *
+ * @return true if the value represented by this object is
+ * infinite, false if it is finite or NaN.
+ */
+ public boolean isInfinity() {
+ return (exponent < 0 && mantissa == 0);
+ }
+ /**
+ * Returns true if the value of this Real is
+ * Not-a-Number (NaN), false otherwise.
+ *
+ *
+ * Equivalent double code: |
+ * Double.{@link Double#isNaN(double) isNaN}(this)
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 0.3
+ * |
+ *
+ * @return true if the value represented by this object is
+ * NaN, false otherwise.
+ */
+ public boolean isNan() {
+ return (exponent < 0 && mantissa != 0);
+ }
+ /**
+ * Returns true if the value of this Real is
+ * finite, false otherwise.
+ *
+ *
+ * Equivalent double code: |
+ * (!Double.{@link Double#isNaN(double) isNaN}(this) &&
+ * !Double.{@link Double#isInfinite(double)
+ * isInfinite}(this))
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 0.3
+ * |
+ *
+ * @return true if the value represented by this object is
+ * finite, false if it is infinite or NaN.
+ */
+ public boolean isFinite() {
+ // That is, non-infinite and non-nan
+ return (exponent >= 0);
+ }
+ /**
+ * Returns true if the value of this Real is
+ * finite and nonzero, false otherwise.
+ *
+ *
+ * Equivalent double code: |
+ * (!Double.{@link Double#isNaN(double) isNaN}(this) &&
+ * !Double.{@link Double#isInfinite(double) isInfinite}(this) &&
+ * (this!=0))
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 0.3
+ * |
+ *
+ * @return true if the value represented by this object is
+ * finite and nonzero, false if it is infinite, NaN or
+ * zero.
+ */
+ public boolean isFiniteNonZero() {
+ // That is, non-infinite and non-nan and non-zero
+ return (exponent >= 0 && mantissa != 0);
+ }
+ /**
+ * Returns true if the value of this Real is
+ * negative, false otherwise.
+ *
+ *
+ * Equivalent double code: |
+ * (this < 0)
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 0.3
+ * |
+ *
+ * @return true if the value represented by this object
+ * is negative, false if it is positive or NaN.
+ */
+ public boolean isNegative() {
+ return sign!=0;
+ }
+ /**
+ * Calculates the absolute value.
+ * Replaces the contents of this Real with the result.
+ *
+ *
+ * Equivalent double code: |
+ * this = Math.{@link Math#abs(double) abs}(this);
+ * |
| Error bound: |
+ * 0 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 0.2
+ * |
+ */
+ public void abs() {
+ sign = 0;
+ }
+ /**
+ * Negates this Real.
+ * Replaces the contents of this Real with the result.
+ *
+ *
+ * Equivalent double code: |
+ * this = -this;
+ * |
| Error bound: |
+ * 0 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 0.2
+ * |
+ */
+ public void neg() {
+ if (!(this.exponent < 0 && this.mantissa != 0))
+ sign ^= 1;
+ }
+ /**
+ * Copies the sign from a.
+ * Replaces the contents of this Real with the result.
+ *
+ *
+ * Equivalent double code: |
+ * this = Math.{@link Math#abs(double)
+ * abs}(this)*Math.{@link Math#signum(double) signum}(a);
+ * |
| Error bound: |
+ * 0 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 0.2
+ * |
+ *
+ * @param a the Real to copy the sign from.
+ */
+ public void copysign(Real a) {
+ if ((this.exponent < 0 && this.mantissa != 0) || (a.exponent < 0 && a.mantissa != 0)) {
+ makeNan();
+ return;
+ }
+ sign = a.sign;
+ }
+ /**
+ * Readjusts the mantissa of this Real. The exponent is
+ * adjusted accordingly. This is necessary when the mantissa has been
+ * {@link #mantissa modified directly} for some purpose and may be
+ * denormalized. The normalized mantissa of a finite Real
+ * must have bit 63 cleared and bit 62 set. Using a denormalized
+ * Real in any other operation may produce undefined
+ * results.
+ *
+ * |
+ * Error bound: |
+ * ½ ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 0.7
+ * |
+ */
+ public void normalize() {
+ if ((this.exponent >= 0)) {
+ if (mantissa > 0)
+ {
+ int clz = 0;
+ int t = (int)(mantissa>>>32);
+ if (t == 0) { clz = 32; t = (int)mantissa; }
+ t|=t>>1; t|=t>>2; t|=t>>4; t|=t>>8; t|=t>>16;
+ clz += clz_tab[(t*clz_magic)>>>27]-1;
+ mantissa <<= clz;
+ exponent -= clz;
+ if (exponent < 0) // Underflow
+ makeZero(sign);
+ }
+ else if (mantissa < 0)
+ {
+ mantissa = (mantissa+1)>>>1;
+ exponent ++;
+ if (mantissa == 0) { // Ooops, it was 0xffffffffffffffffL
+ mantissa = 0x4000000000000000L;
+ exponent ++;
+ }
+ if (exponent < 0) // Overflow
+ makeInfinity(sign);
+ }
+ else // mantissa == 0
+ {
+ exponent = 0;
+ }
+ }
+ }
+ /**
+ * Readjusts the mantissa of a Real with extended
+ * precision. The exponent is adjusted accordingly. This is necessary when
+ * the mantissa has been {@link #mantissa modified directly} for some
+ * purpose and may be denormalized. The normalized mantissa of a finite
+ * Real must have bit 63 cleared and bit 62 set. Using a
+ * denormalized Real in any other operation may
+ * produce undefined results.
+ *
+ * |
+ * Approximate error bound: |
+ * 2-64 ULPs (i.e. of a normal precision Real)
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 0.7
+ * |
+ *
+ * @param extra the extra 64 bits of mantissa of this extended precision
+ * Real.
+ * @return the extra 64 bits of mantissa of the resulting extended
+ * precision Real.
+ */
+ public long normalize128(long extra) {
+ if (!(this.exponent >= 0))
+ return 0;
+ if (mantissa == 0) {
+ if (extra == 0) {
+ exponent = 0;
+ return 0;
+ }
+ mantissa = extra;
+ extra = 0;
+ exponent -= 64;
+ if (exponent < 0) { // Underflow
+ makeZero(sign);
+ return 0;
+ }
+ }
+ if (mantissa < 0) {
+ extra = (mantissa<<63)+(extra>>>1);
+ mantissa >>>= 1;
+ exponent ++;
+ if (exponent < 0) { // Overflow
+ makeInfinity(sign);
+ return 0;
+ }
+ return extra;
+ }
+ int clz = 0;
+ int t = (int)(mantissa>>>32);
+ if (t == 0) { clz = 32; t = (int)mantissa; }
+ t|=t>>1; t|=t>>2; t|=t>>4; t|=t>>8; t|=t>>16;
+ clz += clz_tab[(t*clz_magic)>>>27]-1;
+ if (clz == 0)
+ return extra;
+ mantissa = (mantissa<>>(64-clz));
+ extra <<= clz;
+ exponent -= clz;
+ if (exponent < 0) { // Underflow
+ makeZero(sign);
+ return 0;
+ }
+ return extra;
+ }
+ /**
+ * Rounds an extended precision Real to the nearest
+ * Real of normal precision. Replaces the contents of this
+ * Real with the result.
+ *
+ * |
+ * Error bound: |
+ * ½ ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 1.0
+ * |
+ *
+ * @param extra the extra 64 bits of mantissa of this extended precision
+ * Real.
+ */
+ public void roundFrom128(long extra) {
+ mantissa += (extra>>63)&1;
+ normalize();
+ }
+ /**
+ * Returns true if this Java object is the same
+ * object as a. Since a Real should be
+ * thought of as a "register holding a number", this method compares the
+ * object references, not the contents of the two objects.
+ * This is very different from {@link #equalTo(Real)}.
+ *
+ * @param a the object to compare to this.
+ * @return true if this object is the same as a.
+ */
+ public boolean equals(Object a) {
+ return this==a;
+ }
+ private int compare(Real a) {
+ // Compare of normal floats, zeros, but not nan or equal-signed inf
+ if ((this.exponent == 0 && this.mantissa == 0) && (a.exponent == 0 && a.mantissa == 0))
+ return 0;
+ if (sign != a.sign)
+ return a.sign-sign;
+ int s = (this.sign==0) ? 1 : -1;
+ if ((this.exponent < 0 && this.mantissa == 0))
+ return s;
+ if ((a.exponent < 0 && a.mantissa == 0))
+ return -s;
+ if (exponent != a.exponent)
+ return exponenttrue if this Real is equal to
+ * a.
+ * If the numbers are incomparable, i.e. the values are infinities of
+ * the same sign or any of them is NaN, false is always
+ * returned. This method must not be confused with {@link #equals(Object)}.
+ *
+ *
+ * Equivalent double code: |
+ * (this == a)
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 1.0
+ * |
+ *
+ * @param a the Real to compare to this.
+ * @return true if the value represented by this object is
+ * equal to the value represented by a. false
+ * otherwise, or if the numbers are incomparable.
+ */
+ public boolean equalTo(Real a) {
+ if (invalidCompare(a))
+ return false;
+ return compare(a) == 0;
+ }
+ /**
+ * Returns true if this Real is equal to
+ * the integer a.
+ *
+ *
+ * Equivalent double code: |
+ * (this == a)
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 1.7
+ * |
+ *
+ * @param a the int to compare to this.
+ * @return true if the value represented by this object is
+ * equal to the integer a. false
+ * otherwise.
+ */
+ public boolean equalTo(int a) {
+ tmp0.assign(a);
+ return equalTo(tmp0);
+ }
+ /**
+ * Returns true if this Real is not equal to
+ * a.
+ * If the numbers are incomparable, i.e. the values are infinities of
+ * the same sign or any of them is NaN, false is always
+ * returned.
+ * This distinguishes notEqualTo(a) from the expression
+ * !equalTo(a).
+ *
+ *
+ * Equivalent double code: |
+ * (this != a)
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 1.0
+ * |
+ *
+ * @param a the Real to compare to this.
+ * @return true if the value represented by this object is not
+ * equal to the value represented by a. false
+ * otherwise, or if the numbers are incomparable.
+ */
+ public boolean notEqualTo(Real a) {
+ if (invalidCompare(a))
+ return false;
+ return compare(a) != 0;
+ }
+ /**
+ * Returns true if this Real is not equal to
+ * the integer a.
+ * If this Real is NaN, false is always
+ * returned.
+ * This distinguishes notEqualTo(a) from the expression
+ * !equalTo(a).
+ *
+ *
+ * Equivalent double code: |
+ * (this != a)
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 1.7
+ * |
+ *
+ * @param a the int to compare to this.
+ * @return true if the value represented by this object is not
+ * equal to the integer a. false
+ * otherwise, or if this Real is NaN.
+ */
+ public boolean notEqualTo(int a) {
+ tmp0.assign(a);
+ return notEqualTo(tmp0);
+ }
+ /**
+ * Returns true if this Real is less than
+ * a.
+ * If the numbers are incomparable, i.e. the values are infinities of
+ * the same sign or any of them is NaN, false is always
+ * returned.
+ *
+ *
+ * Equivalent double code: |
+ * (this < a)
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 1.0
+ * |
+ *
+ * @param a the Real to compare to this.
+ * @return true if the value represented by this object is
+ * less than the value represented by a.
+ * false otherwise, or if the numbers are incomparable.
+ */
+ public boolean lessThan(Real a) {
+ if (invalidCompare(a))
+ return false;
+ return compare(a) < 0;
+ }
+ /**
+ * Returns true if this Real is less than
+ * the integer a.
+ * If this Real is NaN, false is always
+ * returned.
+ *
+ *
+ * Equivalent double code: |
+ * (this < a)
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 1.7
+ * |
+ *
+ * @param a the int to compare to this.
+ * @return true if the value represented by this object is
+ * less than the integer a. false otherwise,
+ * or if this Real is NaN.
+ */
+ public boolean lessThan(int a) {
+ tmp0.assign(a);
+ return lessThan(tmp0);
+ }
+ /**
+ * Returns true if this Real is less than or
+ * equal to a.
+ * If the numbers are incomparable, i.e. the values are infinities of
+ * the same sign or any of them is NaN, false is always
+ * returned.
+ *
+ *
+ * Equivalent double code: |
+ * (this <= a)
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 1.0
+ * |
+ *
+ * @param a the Real to compare to this.
+ * @return true if the value represented by this object is
+ * less than or equal to the value represented by a.
+ * false otherwise, or if the numbers are incomparable.
+ */
+ public boolean lessEqual(Real a) {
+ if (invalidCompare(a))
+ return false;
+ return compare(a) <= 0;
+ }
+ /**
+ * Returns true if this Real is less than or
+ * equal to the integer a.
+ * If this Real is NaN, false is always
+ * returned.
+ *
+ *
+ * Equivalent double code: |
+ * (this <= a)
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 1.7
+ * |
+ *
+ * @param a the int to compare to this.
+ * @return true if the value represented by this object is
+ * less than or equal to the integer a. false
+ * otherwise, or if this Real is NaN.
+ */
+ public boolean lessEqual(int a) {
+ tmp0.assign(a);
+ return lessEqual(tmp0);
+ }
+ /**
+ * Returns true if this Real is greater than
+ * a.
+ * If the numbers are incomparable, i.e. the values are infinities of
+ * the same sign or any of them is NaN, false is always
+ * returned.
+ *
+ *
+ * Equivalent double code: |
+ * (this > a)
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 1.0
+ * |
+ *
+ * @param a the Real to compare to this.
+ * @return true if the value represented by this object is
+ * greater than the value represented by a.
+ * false otherwise, or if the numbers are incomparable.
+ */
+ public boolean greaterThan(Real a) {
+ if (invalidCompare(a))
+ return false;
+ return compare(a) > 0;
+ }
+ /**
+ * Returns true if this Real is greater than
+ * the integer a.
+ * If this Real is NaN, false is always
+ * returned.
+ *
+ *
+ * Equivalent double code: |
+ * (this > a)
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 1.7
+ * |
+ *
+ * @param a the int to compare to this.
+ * @return true if the value represented by this object is
+ * greater than the integer a.
+ * false otherwise, or if this Real is NaN.
+ */
+ public boolean greaterThan(int a) {
+ tmp0.assign(a);
+ return greaterThan(tmp0);
+ }
+ /**
+ * Returns true if this Real is greater than
+ * or equal to a.
+ * If the numbers are incomparable, i.e. the values are infinities of
+ * the same sign or any of them is NaN, false is always
+ * returned.
+ *
+ *
+ * Equivalent double code: |
+ * (this >= a)
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 1.0
+ * |
+ *
+ * @param a the Real to compare to this.
+ * @return true if the value represented by this object is
+ * greater than or equal to the value represented by a.
+ * false otherwise, or if the numbers are incomparable.
+ */
+ public boolean greaterEqual(Real a) {
+ if (invalidCompare(a))
+ return false;
+ return compare(a) >= 0;
+ }
+ /**
+ * Returns true if this Real is greater than
+ * or equal to the integer a.
+ * If this Real is NaN, false is always
+ * returned.
+ *
+ *
+ * Equivalent double code: |
+ * (this >= a)
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 1.7
+ * |
+ *
+ * @param a the int to compare to this.
+ * @return true if the value represented by this object is
+ * greater than or equal to the integer a.
+ * false otherwise, or if this Real is NaN.
+ */
+ public boolean greaterEqual(int a) {
+ tmp0.assign(a);
+ return greaterEqual(tmp0);
+ }
+ /**
+ * Returns true if the absolute value of this
+ * Real is less than the absolute value of
+ * a.
+ * If the numbers are incomparable, i.e. the values are both infinite
+ * or any of them is NaN, false is always returned.
+ *
+ *
+ * Equivalent double code: |
+ * (Math.{@link Math#abs(double) abs}(this) <
+ * Math.{@link Math#abs(double) abs}(a))
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 0.5
+ * |
+ *
+ * @param a the Real to compare to this.
+ * @return true if the absolute of the value represented by
+ * this object is less than the absolute of the value represented by
+ * a.
+ * false otherwise, or if the numbers are incomparable.
+ */
+ public boolean absLessThan(Real a) {
+ if ((this.exponent < 0 && this.mantissa != 0) || (a.exponent < 0 && a.mantissa != 0) || (this.exponent < 0 && this.mantissa == 0))
+ return false;
+ if ((a.exponent < 0 && a.mantissa == 0))
+ return true;
+ if (exponent != a.exponent)
+ return exponentReal by 2 to the power of n.
+ * Replaces the contents of this Real with the result.
+ * This operation is faster than normal multiplication since it only
+ * involves adding to the exponent.
+ *
+ *
+ * Equivalent double code: |
+ * this *= Math.{@link Math#pow(double,double) pow}(2.0,n);
+ * |
| Error bound: |
+ * 0 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 0.3
+ * |
+ *
+ * @param n the integer argument.
+ */
+ public void scalbn(int n) {
+ if (!(this.exponent >= 0 && this.mantissa != 0))
+ return;
+ exponent += n;
+ if (exponent < 0) {
+ if (n<0)
+ makeZero(sign); // Underflow
+ else
+ makeInfinity(sign); // Overflow
+ }
+ }
+ /**
+ * Calculates the next representable neighbour of this Real
+ * in the direction towards a.
+ * Replaces the contents of this Real with the result.
+ * If the two values are equal, nothing happens.
+ *
+ *
+ * Equivalent double code: |
+ * this += Math.{@link Math#ulp(double) ulp}(this)*Math.{@link
+ * Math#signum(double) signum}(a-this);
+ * |
| Error bound: |
+ * 0 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 0.8
+ * |
+ *
+ * @param a the Real argument.
+ */
+ public void nextafter(Real a) {
+ if ((this.exponent < 0 && this.mantissa != 0) || (a.exponent < 0 && a.mantissa != 0)) {
+ makeNan();
+ return;
+ }
+ if ((this.exponent < 0 && this.mantissa == 0) && (a.exponent < 0 && a.mantissa == 0) && sign == a.sign)
+ return;
+ int dir = -compare(a);
+ if (dir == 0)
+ return;
+ if ((this.exponent == 0 && this.mantissa == 0)) {
+ { this.mantissa = MIN.mantissa; this.exponent = MIN.exponent; this.sign = MIN.sign; };
+ sign = (byte)(dir<0 ? 1 : 0);
+ return;
+ }
+ if ((this.exponent < 0 && this.mantissa == 0)) {
+ { this.mantissa = MAX.mantissa; this.exponent = MAX.exponent; this.sign = MAX.sign; };
+ sign = (byte)(dir<0 ? 0 : 1);
+ return;
+ }
+ if ((this.sign==0) ^ dir<0) {
+ mantissa ++;
+ } else {
+ if (mantissa == 0x4000000000000000L) {
+ mantissa <<= 1;
+ exponent--;
+ }
+ mantissa --;
+ }
+ normalize();
+ }
+ /**
+ * Calculates the largest (closest to positive infinity)
+ * Real value that is less than or equal to this
+ * Real and is equal to a mathematical integer.
+ * Replaces the contents of this Real with the result.
+ *
+ *
+ * Equivalent double code: |
+ * this = Math.{@link Math#floor(double) floor}(this);
+ * |
| Approximate error bound: |
+ * 0 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 0.5
+ * |
+ */
+ public void floor() {
+ if (!(this.exponent >= 0 && this.mantissa != 0))
+ return;
+ if (exponent < 0x40000000) {
+ if ((this.sign==0))
+ makeZero(sign);
+ else {
+ exponent = ONE.exponent;
+ mantissa = ONE.mantissa;
+ // sign unchanged!
+ }
+ return;
+ }
+ int shift = 0x4000003e-exponent;
+ if (shift <= 0)
+ return;
+ if ((this.sign!=0))
+ mantissa += ((1L<Real value that is greater than or equal to this
+ * Real and is equal to a mathematical integer.
+ * Replaces the contents of this Real with the result.
+ *
+ *
+ * Equivalent double code: |
+ * this = Math.{@link Math#ceil(double) ceil}(this);
+ * |
| Approximate error bound: |
+ * 0 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 1.8
+ * |
+ */
+ public void ceil() {
+ neg();
+ floor();
+ neg();
+ }
+ /**
+ * Rounds this Real value to the closest value that is equal
+ * to a mathematical integer. If two Real values that are
+ * mathematical integers are equally close, the result is the integer
+ * value with the largest magnitude (positive or negative). Replaces the
+ * contents of this Real with the result.
+ *
+ *
+ * Equivalent double code: |
+ * this = Math.{@link Math#rint(double) rint}(this);
+ * |
| Approximate error bound: |
+ * 0 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 1.3
+ * |
+ */
+ public void round() {
+ if (!(this.exponent >= 0 && this.mantissa != 0))
+ return;
+ if (exponent < 0x3fffffff) {
+ makeZero(sign);
+ return;
+ }
+ int shift = 0x4000003e-exponent;
+ if (shift <= 0)
+ return;
+ mantissa += 1L<<(shift-1); // Bla-bla, this works almost
+ mantissa &= ~((1L<Real value to the closest value towards
+ * zero that is equal to a mathematical integer.
+ * Replaces the contents of this Real with the result.
+ *
+ *
+ * Equivalent double code: |
+ * this = (double)((long)this);
+ * |
| Approximate error bound: |
+ * 0 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 1.2
+ * |
+ */
+ public void trunc() {
+ if (!(this.exponent >= 0 && this.mantissa != 0))
+ return;
+ if (exponent < 0x40000000) {
+ makeZero(sign);
+ return;
+ }
+ int shift = 0x4000003e-exponent;
+ if (shift <= 0)
+ return;
+ mantissa &= ~((1L<Real by subtracting
+ * the closest value towards zero that is equal to a mathematical integer.
+ * Replaces the contents of this Real with the result.
+ *
+ *
+ * Equivalent double code: |
+ * this -= (double)((long)this);
+ * |
| Approximate error bound: |
+ * 0 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 1.2
+ * |
+ */
+ public void frac() {
+ if (!(this.exponent >= 0 && this.mantissa != 0) || exponent < 0x40000000)
+ return;
+ int shift = 0x4000003e-exponent;
+ if (shift <= 0) {
+ makeZero(sign);
+ return;
+ }
+ mantissa &= ((1L<Real value to the closest int
+ * value towards zero.
+ *
+ * If the value of this Real is too large, {@link
+ * Integer#MAX_VALUE} is returned. However, if the value of this
+ * Real is too small, -Integer.MAX_VALUE is
+ * returned, not {@link Integer#MIN_VALUE}. This is done to ensure that
+ * the sign will be correct if you calculate
+ * -this.toInteger(). A NaN is converted to 0.
+ *
+ *
+ * Equivalent double code: |
+ * (int)this
+ * |
| Approximate error bound: |
+ * 0 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 0.6
+ * |
+ *
+ * @return an int representation of this Real.
+ */
+ public int toInteger() {
+ if ((this.exponent == 0 && this.mantissa == 0) || (this.exponent < 0 && this.mantissa != 0))
+ return 0;
+ if ((this.exponent < 0 && this.mantissa == 0)) {
+ return ((this.sign==0)) ? 0x7fffffff : 0x80000001;
+ // 0x80000001, so that you can take -x.toInteger()
+ }
+ if (exponent < 0x40000000)
+ return 0;
+ int shift = 0x4000003e-exponent;
+ if (shift < 32) {
+ return ((this.sign==0)) ? 0x7fffffff : 0x80000001;
+ // 0x80000001, so that you can take -x.toInteger()
+ }
+ return (this.sign==0) ?
+ (int)(mantissa>>>shift) : -(int)(mantissa>>>shift);
+ }
+ /**
+ * Converts this Real value to the closest long
+ * value towards zero.
+ *
+ * If the value of this Real is too large, {@link
+ * Long#MAX_VALUE} is returned. However, if the value of this
+ * Real is too small, -Long.MAX_VALUE is
+ * returned, not {@link Long#MIN_VALUE}. This is done to ensure that the
+ * sign will be correct if you calculate -this.toLong().
+ * A NaN is converted to 0.
+ *
+ *
+ * Equivalent double code: |
+ * (long)this
+ * |
| Approximate error bound: |
+ * 0 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 0.5
+ * |
+ *
+ * @return a long representation of this Real.
+ */
+ public long toLong() {
+ if ((this.exponent == 0 && this.mantissa == 0) || (this.exponent < 0 && this.mantissa != 0))
+ return 0;
+ if ((this.exponent < 0 && this.mantissa == 0)) {
+ return ((this.sign==0))? 0x7fffffffffffffffL:0x8000000000000001L;
+ // 0x8000000000000001L, so that you can take -x.toLong()
+ }
+ if (exponent < 0x40000000)
+ return 0;
+ int shift = 0x4000003e-exponent;
+ if (shift < 0) {
+ return ((this.sign==0))? 0x7fffffffffffffffL:0x8000000000000001L;
+ // 0x8000000000000001L, so that you can take -x.toLong()
+ }
+ return (this.sign==0) ? (mantissa>>>shift) : -(mantissa>>>shift);
+ }
+ /**
+ * Returns true if the value of this Real
+ * represents a mathematical integer. If the value is too large to
+ * determine if it is an integer, true is returned.
+ *
+ *
+ * Equivalent double code: |
+ * (this == (long)this)
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 0.6
+ * |
+ *
+ * @return true if the value represented by this object
+ * represents a mathematical integer, false otherwise.
+ */
+ public boolean isIntegral() {
+ if ((this.exponent < 0 && this.mantissa != 0))
+ return false;
+ if ((this.exponent == 0 && this.mantissa == 0) || (this.exponent < 0 && this.mantissa == 0))
+ return true;
+ if (exponent < 0x40000000)
+ return false;
+ int shift = 0x4000003e-exponent;
+ if (shift <= 0)
+ return true;
+ return (mantissa&((1L<true if the mathematical integer represented
+ * by this Real is odd. You must first determine
+ * that the value is actually an integer using {@link
+ * #isIntegral()}. If the value is too large to determine if the
+ * integer is odd, false is returned.
+ *
+ *
+ * Equivalent double code: |
+ * ((((long)this)&1) == 1)
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 0.6
+ * |
+ *
+ * @return true if the mathematical integer represented by
+ * this Real is odd, false otherwise.
+ */
+ public boolean isOdd() {
+ if (!(this.exponent >= 0 && this.mantissa != 0) ||
+ exponent < 0x40000000 || exponent > 0x4000003e)
+ return false;
+ int shift = 0x4000003e-exponent;
+ return ((mantissa>>>shift)&1) != 0;
+ }
+ /**
+ * Exchanges the contents of this Real and a.
+ *
+ *
+ * Equivalent double code: |
+ * tmp=this; this=a; a=tmp;
+ * |
| Error bound: |
+ * 0 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 0.5
+ * |
+ *
+ * @param a the Real to exchange with this.
+ */
+ public void swap(Real a) {
+ long tmpMantissa=mantissa; mantissa=a.mantissa; a.mantissa=tmpMantissa;
+ int tmpExponent=exponent; exponent=a.exponent; a.exponent=tmpExponent;
+ byte tmpSign =sign; sign =a.sign; a.sign =tmpSign;
+ }
+ // Temporary values used by functions (to avoid "new" inside functions)
+ private static Real tmp0 = new Real(); // tmp for basic functions
+ private static Real recipTmp = new Real();
+ private static Real recipTmp2 = new Real();
+ private static Real sqrtTmp = new Real();
+ private static Real expTmp = new Real();
+ private static Real expTmp2 = new Real();
+ private static Real expTmp3 = new Real();
+ private static Real tmp1 = new Real();
+ private static Real tmp2 = new Real();
+ private static Real tmp3 = new Real();
+ private static Real tmp4 = new Real();
+ private static Real tmp5 = new Real();
+ /**
+ * Calculates the sum of this Real and a.
+ * Replaces the contents of this Real with the result.
+ *
+ *
+ * Equivalent double code: |
+ * this += a;
+ * |
| Error bound: |
+ * ½ ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * «« 1.0 »»
+ * |
+ *
+ * @param a the Real to add to this.
+ */
+ public void add(Real a) {
+ if ((this.exponent < 0 && this.mantissa != 0) || (a.exponent < 0 && a.mantissa != 0)) {
+ makeNan();
+ return;
+ }
+ if ((this.exponent < 0 && this.mantissa == 0) || (a.exponent < 0 && a.mantissa == 0)) {
+ if ((this.exponent < 0 && this.mantissa == 0) && (a.exponent < 0 && a.mantissa == 0) && sign != a.sign)
+ makeNan();
+ else
+ makeInfinity((this.exponent < 0 && this.mantissa == 0) ? sign : a.sign);
+ return;
+ }
+ if ((this.exponent == 0 && this.mantissa == 0) || (a.exponent == 0 && a.mantissa == 0)) {
+ if ((this.exponent == 0 && this.mantissa == 0))
+ { this.mantissa = a.mantissa; this.exponent = a.exponent; this.sign = a.sign; };
+ if ((this.exponent == 0 && this.mantissa == 0))
+ sign=0;
+ return;
+ }
+ byte s;
+ int e;
+ long m;
+ if (exponent > a.exponent ||
+ (exponent == a.exponent && mantissa>=a.mantissa))
+ {
+ s = a.sign;
+ e = a.exponent;
+ m = a.mantissa;
+ } else {
+ s = sign;
+ e = exponent;
+ m = mantissa;
+ sign = a.sign;
+ exponent = a.exponent;
+ mantissa = a.mantissa;
+ }
+ int shift = exponent-e;
+ if (shift>=64)
+ return;
+ if (sign == s) {
+ mantissa += m>>>shift;
+ if (mantissa >= 0 && shift>0 && ((m>>>(shift-1))&1) != 0)
+ mantissa ++; // We don't need normalization, so round now
+ if (mantissa < 0) {
+ // Simplified normalize()
+ mantissa = (mantissa+1)>>>1;
+ exponent ++;
+ if (exponent < 0) { // Overflow
+ makeInfinity(sign);
+ return;
+ }
+ }
+ } else {
+ if (shift>0) {
+ // Shift mantissa up to increase accuracy
+ mantissa <<= 1;
+ exponent --;
+ shift --;
+ }
+ m = -m;
+ mantissa += m>>shift;
+ if (mantissa >= 0 && shift>0 && ((m>>>(shift-1))&1) != 0)
+ mantissa ++; // We don't need to shift down, so round now
+ if (mantissa < 0) {
+ // Simplified normalize()
+ mantissa = (mantissa+1)>>>1;
+ exponent ++; // Can't overflow
+ } else if (shift==0) {
+ // Operands have equal exponents => many bits may be cancelled
+ // Magic rounding: if result of subtract leaves only a few bits
+ // standing, the result should most likely be 0...
+ if (magicRounding && mantissa > 0 && mantissa <= 7) {
+ // If arguments were integers <= 2^63-1, then don't
+ // do the magic rounding anyway.
+ // This is a bit "post mortem" investigation but it happens
+ // so seldom that it's no problem to spend the extra time.
+ m = -m;
+ if (exponent == 0x4000003c || exponent == 0x4000003d ||
+ (exponent == 0x4000003e && mantissa+m > 0)) {
+ long mask = (1<<(0x4000003e-exponent))-1;
+ if ((mantissa & mask) != 0 || (m & mask) != 0)
+ mantissa = 0;
+ } else
+ mantissa = 0;
+ }
+ normalize();
+ } // else... if (shift>=1 && mantissa>=0) it should be a-ok
+ }
+ if ((this.exponent == 0 && this.mantissa == 0))
+ sign=0;
+ }
+ /**
+ * Calculates the sum of this Real and the integer
+ * a.
+ * Replaces the contents of this Real with the result.
+ *
+ *
+ * Equivalent double code: |
+ * this += a;
+ * |
| Error bound: |
+ * ½ ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 1.8
+ * |
+ *
+ * @param a the int to add to this.
+ */
+ public void add(int a) {
+ tmp0.assign(a);
+ add(tmp0);
+ }
+ /**
+ * Calculates the sum of this Real and a with
+ * extended precision. Replaces the contents of this Real
+ * with the result. Returns the extra mantissa of the extended precision
+ * result.
+ *
+ * An extra 64 bits of mantissa is added to both arguments for extended
+ * precision. If any of the arguments are not of extended precision, use
+ * 0 for the extra mantissa.
+ *
+ *
Extended prevision can be useful in many situations. For instance,
+ * when accumulating a lot of very small values it is advantageous for the
+ * accumulator to have extended precision. To convert the extended
+ * precision value back to a normal Real for further
+ * processing, use {@link #roundFrom128(long)}.
+ *
+ *
+ * Equivalent double code: |
+ * this += a;
+ * |
| Approximate error bound: |
+ * 2-62 ULPs (i.e. of a normal precision Real)
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 2.0
+ * |
+ *
+ * @param extra the extra 64 bits of mantissa of this extended precision
+ * Real.
+ * @param a the Real to add to this.
+ * @param aExtra the extra 64 bits of mantissa of the extended precision
+ * value a.
+ * @return the extra 64 bits of mantissa of the resulting extended
+ * precision Real.
+ */
+ public long add128(long extra, Real a, long aExtra) {
+ if ((this.exponent < 0 && this.mantissa != 0) || (a.exponent < 0 && a.mantissa != 0)) {
+ makeNan();
+ return 0;
+ }
+ if ((this.exponent < 0 && this.mantissa == 0) || (a.exponent < 0 && a.mantissa == 0)) {
+ if ((this.exponent < 0 && this.mantissa == 0) && (a.exponent < 0 && a.mantissa == 0) && sign != a.sign)
+ makeNan();
+ else
+ makeInfinity((this.exponent < 0 && this.mantissa == 0) ? sign : a.sign);
+ return 0;
+ }
+ if ((this.exponent == 0 && this.mantissa == 0) || (a.exponent == 0 && a.mantissa == 0)) {
+ if ((this.exponent == 0 && this.mantissa == 0)) {
+ { this.mantissa = a.mantissa; this.exponent = a.exponent; this.sign = a.sign; };
+ extra = aExtra;
+ }
+ if ((this.exponent == 0 && this.mantissa == 0))
+ sign=0;
+ return extra;
+ }
+ byte s;
+ int e;
+ long m;
+ long x;
+ if (exponent > a.exponent ||
+ (exponent == a.exponent && mantissa>a.mantissa) ||
+ (exponent == a.exponent && mantissa==a.mantissa &&
+ (extra>>>1)>=(aExtra>>>1)))
+ {
+ s = a.sign;
+ e = a.exponent;
+ m = a.mantissa;
+ x = aExtra;
+ } else {
+ s = sign;
+ e = exponent;
+ m = mantissa;
+ x = extra;
+ sign = a.sign;
+ exponent = a.exponent;
+ mantissa = a.mantissa;
+ extra = aExtra;
+ }
+ int shift = exponent-e;
+ if (shift>=127)
+ return extra;
+ if (shift>=64) {
+ x = m>>>(shift-64);
+ m = 0;
+ } else if (shift>0) {
+ x = (x>>>shift)+(m<<(64-shift));
+ m >>>= shift;
+ }
+ extra >>>= 1;
+ x >>>= 1;
+ if (sign == s) {
+ extra += x;
+ mantissa += (extra>>63)&1;
+ mantissa += m;
+ } else {
+ extra -= x;
+ mantissa -= (extra>>63)&1;
+ mantissa -= m;
+ // Magic rounding: if result of subtract leaves only a few bits
+ // standing, the result should most likely be 0...
+ if (mantissa == 0 && extra > 0 && extra <= 0x1f)
+ extra = 0;
+ }
+ extra <<= 1;
+ extra = normalize128(extra);
+ if ((this.exponent == 0 && this.mantissa == 0))
+ sign=0;
+ return extra;
+ }
+ /**
+ * Calculates the difference between this Real and
+ * a.
+ * Replaces the contents of this Real with the result.
+ *
+ * (To achieve extended precision subtraction, it is enough to call
+ * a.{@link #neg() neg}() before calling {@link
+ * #add128(long,Real,long) add128}(extra,a,aExtra), since only
+ * the sign bit of a need to be changed.)
+ *
+ *
+ * Equivalent double code: |
+ * this -= a;
+ * |
| Error bound: |
+ * ½ ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 2.0
+ * |
+ *
+ * @param a the Real to subtract from this.
+ */
+ public void sub(Real a) {
+ tmp0.mantissa = a.mantissa;
+ tmp0.exponent = a.exponent;
+ tmp0.sign = (byte)(a.sign^1);
+ add(tmp0);
+ }
+ /**
+ * Calculates the difference between this Real and the
+ * integer a.
+ * Replaces the contents of this Real with the result.
+ *
+ *
+ * Equivalent double code: |
+ * this -= a;
+ * |
| Error bound: |
+ * ½ ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 2.4
+ * |
+ *
+ * @param a the int to subtract from this.
+ */
+ public void sub(int a) {
+ tmp0.assign(a);
+ sub(tmp0);
+ }
+ /**
+ * Calculates the product of this Real and a.
+ * Replaces the contents of this Real with the result.
+ *
+ *
+ * Equivalent double code: |
+ * this *= a;
+ * |
| Error bound: |
+ * ½ ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 1.3
+ * |
+ *
+ * @param a the Real to multiply to this.
+ */
+ public void mul(Real a) {
+ if ((this.exponent < 0 && this.mantissa != 0) || (a.exponent < 0 && a.mantissa != 0)) {
+ makeNan();
+ return;
+ }
+ sign ^= a.sign;
+ if ((this.exponent == 0 && this.mantissa == 0) || (a.exponent == 0 && a.mantissa == 0)) {
+ if ((this.exponent < 0 && this.mantissa == 0) || (a.exponent < 0 && a.mantissa == 0))
+ makeNan();
+ else
+ makeZero(sign);
+ return;
+ }
+ if ((this.exponent < 0 && this.mantissa == 0) || (a.exponent < 0 && a.mantissa == 0)) {
+ makeInfinity(sign);
+ return;
+ }
+ long a0 = mantissa & 0x7fffffff;
+ long a1 = mantissa >>> 31;
+ long b0 = a.mantissa & 0x7fffffff;
+ long b1 = a.mantissa >>> 31;
+ mantissa = a1*b1;
+ // If we're going to need normalization, we don't want to round twice
+ int round = (mantissa<0) ? 0 : 0x40000000;
+ mantissa += ((a0*b1 + a1*b0 + ((a0*b0)>>>31) + round)>>>31);
+ int aExp = a.exponent;
+ exponent += aExp-0x40000000;
+ if (exponent < 0) {
+ if (exponent == -1 && aExp < 0x40000000 && mantissa < 0) {
+ // Not underflow after all, it will be corrected in the
+ // normalization below
+ } else {
+ if (aExp < 0x40000000)
+ makeZero(sign); // Underflow
+ else
+ makeInfinity(sign); // Overflow
+ return;
+ }
+ }
+ // Simplified normalize()
+ if (mantissa < 0) {
+ mantissa = (mantissa+1)>>>1;
+ exponent ++;
+ if (exponent < 0) // Overflow
+ makeInfinity(sign);
+ }
+ }
+ /**
+ * Calculates the product of this Real and the integer
+ * a.
+ * Replaces the contents of this Real with the result.
+ *
+ *
+ * Equivalent double code: |
+ * this *= a;
+ * |
| Error bound: |
+ * ½ ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 1.3
+ * |
+ *
+ * @param a the int to multiply to this.
+ */
+ public void mul(int a) {
+ if ((this.exponent < 0 && this.mantissa != 0))
+ return;
+ if (a<0) {
+ sign ^= 1;
+ a = -a;
+ }
+ if ((this.exponent == 0 && this.mantissa == 0) || a==0) {
+ if ((this.exponent < 0 && this.mantissa == 0))
+ makeNan();
+ else
+ makeZero(sign);
+ return;
+ }
+ if ((this.exponent < 0 && this.mantissa == 0))
+ return;
+ // Normalize int
+ int t=a; t|=t>>1; t|=t>>2; t|=t>>4; t|=t>>8; t|=t>>16;
+ t = clz_tab[(t*clz_magic)>>>27];
+ exponent += 0x1F-t;
+ a <<= t;
+ if (exponent < 0) {
+ makeInfinity(sign); // Overflow
+ return;
+ }
+ long a0 = mantissa & 0x7fffffff;
+ long a1 = mantissa >>> 31;
+ long b0 = a & 0xffffffffL;
+ mantissa = a1*b0;
+ // If we're going to need normalization, we don't want to round twice
+ int round = (mantissa<0) ? 0 : 0x40000000;
+ mantissa += ((a0*b0 + round)>>>31);
+ // Simplified normalize()
+ if (mantissa < 0) {
+ mantissa = (mantissa+1)>>>1;
+ exponent ++;
+ if (exponent < 0) // Overflow
+ makeInfinity(sign);
+ }
+ }
+ /**
+ * Calculates the product of this Real and a with
+ * extended precision.
+ * Replaces the contents of this Real with the result.
+ * Returns the extra mantissa of the extended precision result.
+ *
+ * An extra 64 bits of mantissa is added to both arguments for
+ * extended precision. If any of the arguments are not of extended
+ * precision, use 0 for the extra mantissa. See also {@link
+ * #add128(long,Real,long)}.
+ *
+ *
+ * Equivalent double code: |
+ * this *= a;
+ * |
| Approximate error bound: |
+ * 2-60 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 3.1
+ * |
+ *
+ * @param extra the extra 64 bits of mantissa of this extended precision
+ * Real.
+ * @param a the Real to multiply to this.
+ * @param aExtra the extra 64 bits of mantissa of the extended precision
+ * value a.
+ * @return the extra 64 bits of mantissa of the resulting extended
+ * precision Real.
+ */
+ public long mul128(long extra, Real a, long aExtra) {
+ if ((this.exponent < 0 && this.mantissa != 0) || (a.exponent < 0 && a.mantissa != 0)) {
+ makeNan();
+ return 0;
+ }
+ sign ^= a.sign;
+ if ((this.exponent == 0 && this.mantissa == 0) || (a.exponent == 0 && a.mantissa == 0)) {
+ if ((this.exponent < 0 && this.mantissa == 0) || (a.exponent < 0 && a.mantissa == 0))
+ makeNan();
+ else
+ makeZero(sign);
+ return 0;
+ }
+ if ((this.exponent < 0 && this.mantissa == 0) || (a.exponent < 0 && a.mantissa == 0)) {
+ makeInfinity(sign);
+ return 0;
+ }
+ int aExp = a.exponent;
+ exponent += aExp-0x40000000;
+ if (exponent < 0) {
+ if (aExp < 0x40000000)
+ makeZero(sign); // Underflow
+ else
+ makeInfinity(sign); // Overflow
+ return 0;
+ }
+ long ffffffffL = 0xffffffffL;
+ long a0 = extra & ffffffffL;
+ long a1 = extra >>> 32;
+ long a2 = mantissa & ffffffffL;
+ long a3 = mantissa >>> 32;
+ long b0 = aExtra & ffffffffL;
+ long b1 = aExtra >>> 32;
+ long b2 = a.mantissa & ffffffffL;
+ long b3 = a.mantissa >>> 32;
+ a0 = ((a3*b0>>>2)+
+ (a2*b1>>>2)+
+ (a1*b2>>>2)+
+ (a0*b3>>>2)+
+ 0x60000000)>>>28;
+ //(a2*b0>>>34)+(a1*b1>>>34)+(a0*b2>>>34)+0x08000000)>>>28;
+ a1 *= b3;
+ b0 = a2*b2;
+ b1 *= a3;
+ a0 += ((a1<<2)&ffffffffL) + ((b0<<2)&ffffffffL) + ((b1<<2)&ffffffffL);
+ a1 = (a0>>>32) + (a1>>>30) + (b0>>>30) + (b1>>>30);
+ a0 &= ffffffffL;
+ a2 *= b3;
+ b2 *= a3;
+ a1 += ((a2<<2)&ffffffffL) + ((b2<<2)&ffffffffL);
+ extra = (a1<<32) + a0;
+ mantissa = ((a3*b3)<<2) + (a1>>>32) + (a2>>>30) + (b2>>>30);
+ extra = normalize128(extra);
+ return extra;
+ }
+ private void mul10() {
+ if (!(this.exponent >= 0 && this.mantissa != 0))
+ return;
+ mantissa += (mantissa+2)>>>2;
+ exponent += 3;
+ if (mantissa < 0) {
+ mantissa = (mantissa+1)>>>1;
+ exponent++;
+ }
+ if (exponent < 0)
+ makeInfinity(sign); // Overflow
+ }
+ /**
+ * Calculates the square of this Real.
+ * Replaces the contents of this Real with the result.
+ *
+ *
+ * Equivalent double code: |
+ * this = this*this;
+ * |
| Error bound: |
+ * ½ ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 1.1
+ * |
+ */
+ public void sqr() {
+ sign = 0;
+ if (!(this.exponent >= 0 && this.mantissa != 0))
+ return;
+ int e = exponent;
+ exponent += exponent-0x40000000;
+ if (exponent < 0) {
+ if (e < 0x40000000)
+ makeZero(sign); // Underflow
+ else
+ makeInfinity(sign); // Overflow
+ return;
+ }
+ long a0 = mantissa&0x7fffffff;
+ long a1 = mantissa>>>31;
+ mantissa = a1*a1;
+ // If we're going to need normalization, we don't want to round twice
+ int round = (mantissa<0) ? 0 : 0x40000000;
+ mantissa += ((((a0*a1)<<1) + ((a0*a0)>>>31) + round)>>>31);
+ // Simplified normalize()
+ if (mantissa < 0) {
+ mantissa = (mantissa+1)>>>1;
+ exponent ++;
+ if (exponent < 0) // Overflow
+ makeInfinity(sign);
+ }
+ }
+ private static long ldiv(long a, long b) {
+ // Calculate (a<<63)/b, where a<2**64, b<2**63, b<=a and a<2*b The
+ // result will always be 63 bits, leading to a 3-stage radix-2**21
+ // (very high radix) algorithm, as described here:
+ // S.F. Oberman and M.J. Flynn, "Division Algorithms and
+ // Implementations," IEEE Trans. Computers, vol. 46, no. 8,
+ // pp. 833-854, Aug 1997 Section 4: "Very High Radix Algorithms"
+ int bInv24; // Approximate 1/b, never more than 24 bits
+ int aHi24; // High 24 bits of a (sometimes 25 bits)
+ int next21; // The next 21 bits of result, possibly 1 less
+ long q; // Resulting quotient: round((a<<63)/b)
+ // Preparations
+ bInv24 = (int)(0x400000000000L/((b>>>40)+1));
+ aHi24 = (int)(a>>32)>>>8;
+ a <<= 20; // aHi24 and a overlap by 4 bits
+ // Now perform the division
+ next21 = (int)(((long)aHi24*(long)bInv24)>>>26);
+ a -= next21*b; // Bits above 2**64 will always be cancelled
+ // No need to remove remainder, this will be cared for in next block
+ q = next21;
+ aHi24 = (int)(a>>32)>>>7;
+ a <<= 21;
+ // Two more almost identical blocks...
+ next21 = (int)(((long)aHi24*(long)bInv24)>>>26);
+ a -= next21*b;
+ q = (q<<21)+next21;
+ aHi24 = (int)(a>>32)>>>7;
+ a <<= 21;
+ next21 = (int)(((long)aHi24*(long)bInv24)>>>26);
+ a -= next21*b;
+ q = (q<<21)+next21;
+ // Remove final remainder
+ if (a<0 || a>=b) { q++; a -= b; }
+ a <<= 1;
+ // Round correctly
+ if (a<0 || a>=b) q++;
+ return q;
+ }
+ /**
+ * Calculates the quotient of this Real and a.
+ * Replaces the contents of this Real with the result.
+ *
+ * (To achieve extended precision division, call
+ * aExtra=a.{@link #recip128(long) recip128}(aExtra) before
+ * calling {@link #mul128(long,Real,long)
+ * mul128}(extra,a,aExtra).)
+ *
+ *
+ * Equivalent double code: |
+ * this /= a;
+ * |
| Error bound: |
+ * ½ ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 2.6
+ * |
+ *
+ * @param a the Real to divide this with.
+ */
+ public void div(Real a) {
+ if ((this.exponent < 0 && this.mantissa != 0) || (a.exponent < 0 && a.mantissa != 0)) {
+ makeNan();
+ return;
+ }
+ sign ^= a.sign;
+ if ((this.exponent < 0 && this.mantissa == 0)) {
+ if ((a.exponent < 0 && a.mantissa == 0))
+ makeNan();
+ return;
+ }
+ if ((a.exponent < 0 && a.mantissa == 0)) {
+ makeZero(sign);
+ return;
+ }
+ if ((this.exponent == 0 && this.mantissa == 0)) {
+ if ((a.exponent == 0 && a.mantissa == 0))
+ makeNan();
+ return;
+ }
+ if ((a.exponent == 0 && a.mantissa == 0)) {
+ makeInfinity(sign);
+ return;
+ }
+ exponent += 0x40000000-a.exponent;
+ if (mantissa < a.mantissa) {
+ mantissa <<= 1;
+ exponent--;
+ }
+ if (exponent < 0) {
+ if (a.exponent >= 0x40000000)
+ makeZero(sign); // Underflow
+ else
+ makeInfinity(sign); // Overflow
+ return;
+ }
+ if (a.mantissa == 0x4000000000000000L)
+ return;
+ mantissa = ldiv(mantissa,a.mantissa);
+ }
+ /**
+ * Calculates the quotient of this Real and the integer
+ * a.
+ * Replaces the contents of this Real with the result.
+ *
+ *
+ * Equivalent double code: |
+ * this /= a;
+ * |
| Error bound: |
+ * ½ ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 2.6
+ * |
+ *
+ * @param a the int to divide this with.
+ */
+ public void div(int a) {
+ if ((this.exponent < 0 && this.mantissa != 0))
+ return;
+ if (a<0) {
+ sign ^= 1;
+ a = -a;
+ }
+ if ((this.exponent < 0 && this.mantissa == 0))
+ return;
+ if ((this.exponent == 0 && this.mantissa == 0)) {
+ if (a==0)
+ makeNan();
+ return;
+ }
+ if (a==0) {
+ makeInfinity(sign);
+ return;
+ }
+ long denom = a & 0xffffffffL;
+ long remainder = mantissa%denom;
+ mantissa /= denom;
+ // Normalizing mantissa and scaling remainder accordingly
+ int clz = 0;
+ int t = (int)(mantissa>>>32);
+ if (t == 0) { clz = 32; t = (int)mantissa; }
+ t|=t>>1; t|=t>>2; t|=t>>4; t|=t>>8; t|=t>>16;
+ clz += clz_tab[(t*clz_magic)>>>27]-1;
+ mantissa <<= clz;
+ remainder <<= clz;
+ exponent -= clz;
+ // Final division, correctly rounded
+ remainder = (remainder+denom/2)/denom;
+ mantissa += remainder;
+ if (exponent < 0) // Underflow
+ makeZero(sign);
+ }
+ /**
+ * Calculates the quotient of a and this Real.
+ * Replaces the contents of this Real with the result.
+ *
+ *
+ * Equivalent double code: |
+ * this = a/this;
+ * |
| Error bound: |
+ * ½ ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 3.1
+ * |
+ *
+ * @param a the Real to be divided by this.
+ */
+ public void rdiv(Real a) {
+ { recipTmp.mantissa = a.mantissa; recipTmp.exponent = a.exponent; recipTmp.sign = a.sign; };
+ recipTmp.div(this);
+ { this.mantissa = recipTmp.mantissa; this.exponent = recipTmp.exponent; this.sign = recipTmp.sign; };
+ }
+ /**
+ * Calculates the quotient of the integer a and this
+ * Real.
+ * Replaces the contents of this Real with the result.
+ *
+ *
+ * Equivalent double code: |
+ * this = a/this;
+ * |
| Error bound: |
+ * ½ ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 3.9
+ * |
+ *
+ * @param a the int to be divided by this.
+ */
+ public void rdiv(int a) {
+ tmp0.assign(a);
+ rdiv(tmp0);
+ }
+ /**
+ * Calculates the reciprocal of this Real.
+ * Replaces the contents of this Real with the result.
+ *
+ *
+ * Equivalent double code: |
+ * this = 1/this;
+ * |
| Error bound: |
+ * ½ ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 2.3
+ * |
+ */
+ public void recip() {
+ if ((this.exponent < 0 && this.mantissa != 0))
+ return;
+ if ((this.exponent < 0 && this.mantissa == 0)) {
+ makeZero(sign);
+ return;
+ }
+ if ((this.exponent == 0 && this.mantissa == 0)) {
+ makeInfinity(sign);
+ return;
+ }
+ exponent = 0x80000000-exponent;
+ if (mantissa == 0x4000000000000000L) {
+ if (exponent < 0)
+ makeInfinity(sign); // Overflow
+ return;
+ }
+ exponent--;
+ mantissa = ldiv(0x8000000000000000L,mantissa);
+ }
+ /**
+ * Calculates the reciprocal of this Real with
+ * extended precision.
+ * Replaces the contents of this Real with the result.
+ * Returns the extra mantissa of the extended precision result.
+ *
+ * An extra 64 bits of mantissa is added for extended precision.
+ * If the argument is not of extended precision, use 0
+ * for the extra mantissa. See also {@link #add128(long,Real,long)}.
+ *
+ *
+ * Equivalent double code: |
+ * this = 1/this;
+ * |
| Approximate error bound: |
+ * 2-60 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 17
+ * |
+ *
+ * @param extra the extra 64 bits of mantissa of this extended precision
+ * Real.
+ * @return the extra 64 bits of mantissa of the resulting extended
+ * precision Real.
+ */
+ public long recip128(long extra) {
+ if ((this.exponent < 0 && this.mantissa != 0))
+ return 0;
+ if ((this.exponent < 0 && this.mantissa == 0)) {
+ makeZero(sign);
+ return 0;
+ }
+ if ((this.exponent == 0 && this.mantissa == 0)) {
+ makeInfinity(sign);
+ return 0;
+ }
+ byte s = sign;
+ sign = 0;
+ // Special case, simple power of 2
+ if (mantissa == 0x4000000000000000L && extra == 0) {
+ exponent = 0x80000000-exponent;
+ if (exponent<0) // Overflow
+ makeInfinity(s);
+ return 0;
+ }
+ // Normalize exponent
+ int exp = 0x40000000-exponent;
+ exponent = 0x40000000;
+ // Save -A
+ { recipTmp.mantissa = this.mantissa; recipTmp.exponent = this.exponent; recipTmp.sign = this.sign; };
+ long recipTmpExtra = extra;
+ recipTmp.neg();
+ // First establish approximate result (actually 63 bit accurate)
+ recip();
+ // Perform one Newton-Raphson iteration
+ // Xn+1 = Xn + Xn*(1-A*Xn)
+ { recipTmp2.mantissa = this.mantissa; recipTmp2.exponent = this.exponent; recipTmp2.sign = this.sign; };
+ extra = mul128(0,recipTmp,recipTmpExtra);
+ extra = add128(extra,ONE,0);
+ extra = mul128(extra,recipTmp2,0);
+ extra = add128(extra,recipTmp2,0);
+ // Fix exponent
+ scalbn(exp);
+ // Fix sign
+ if (!isNan())
+ sign = s;
+ return extra;
+ }
+ /**
+ * Calculates the mathematical integer that is less than or equal to
+ * this Real divided by a.
+ * Replaces the contents of this Real with the result.
+ *
+ *
+ * Equivalent double code: |
+ * this = Math.{@link Math#floor(double) floor}(this/a);
+ * |
| Error bound: |
+ * 0 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 22
+ * |
+ *
+ * @param a the Real argument.
+ */
+ public void divf(Real a) {
+ if ((this.exponent < 0 && this.mantissa != 0) || (a.exponent < 0 && a.mantissa != 0)) {
+ makeNan();
+ return;
+ }
+ if ((this.exponent < 0 && this.mantissa == 0)) {
+ if ((a.exponent < 0 && a.mantissa == 0))
+ makeNan();
+ return;
+ }
+ if ((a.exponent < 0 && a.mantissa == 0)) {
+ makeZero(sign);
+ return;
+ }
+ if ((this.exponent == 0 && this.mantissa == 0)) {
+ if ((a.exponent == 0 && a.mantissa == 0))
+ makeNan();
+ return;
+ }
+ if ((a.exponent == 0 && a.mantissa == 0)) {
+ makeInfinity(sign);
+ return;
+ }
+ { tmp0.mantissa = a.mantissa; tmp0.exponent = a.exponent; tmp0.sign = a.sign; }; // tmp0 should be free
+ // Perform same division as with mod, and don't round up
+ long extra = tmp0.recip128(0);
+ extra = mul128(0,tmp0,extra);
+ if (((tmp0.sign!=0) && (extra < 0 || extra > 0x1f)) ||
+ (!(tmp0.sign!=0) && extra < 0 && extra > 0xffffffe0))
+ {
+ // For accurate floor()
+ mantissa++;
+ normalize();
+ }
+ floor();
+ }
+ private void modInternal(/*long thisExtra,*/ Real a, long aExtra) {
+ { tmp0.mantissa = a.mantissa; tmp0.exponent = a.exponent; tmp0.sign = a.sign; }; // tmp0 should be free
+ long extra = tmp0.recip128(aExtra);
+ extra = tmp0.mul128(extra,this,0/*thisExtra*/); // tmp0 == this/a
+ if (tmp0.exponent > 0x4000003e) {
+ // floor() will be inaccurate
+ makeZero(a.sign); // What else can be done? makeNan?
+ return;
+ }
+ if (((tmp0.sign!=0) && (extra < 0 || extra > 0x1f)) ||
+ (!(tmp0.sign!=0) && extra < 0 && extra > 0xffffffe0))
+ {
+ // For accurate floor() with a bit of "magical rounding"
+ tmp0.mantissa++;
+ tmp0.normalize();
+ }
+ tmp0.floor();
+ tmp0.neg(); // tmp0 == -floor(this/a)
+ extra = tmp0.mul128(0,a,aExtra);
+ extra = add128(0/*thisExtra*/,tmp0,extra);
+ roundFrom128(extra);
+ }
+ /**
+ * Calculates the value of this Real modulo a.
+ * Replaces the contents of this Real with the result.
+ * The modulo in this case is defined as the remainder after subtracting
+ * a multiplied by the mathematical integer that is less than
+ * or equal to this Real divided by a.
+ *
+ *
+ * Equivalent double code: |
+ * this = this -
+ * a*Math.{@link Math#floor(double) floor}(this/a);
+ * |
| Error bound: |
+ * 0 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 27
+ * |
+ *
+ * @param a the Real argument.
+ */
+ public void mod(Real a) {
+ if ((this.exponent < 0 && this.mantissa != 0) || (a.exponent < 0 && a.mantissa != 0)) {
+ makeNan();
+ return;
+ }
+ if ((this.exponent < 0 && this.mantissa == 0)) {
+ makeNan();
+ return;
+ }
+ if ((this.exponent == 0 && this.mantissa == 0)) {
+ if ((a.exponent == 0 && a.mantissa == 0))
+ makeNan();
+ else
+ sign = a.sign;
+ return;
+ }
+ if ((a.exponent < 0 && a.mantissa == 0)) {
+ if (sign != a.sign)
+ makeInfinity(a.sign);
+ return;
+ }
+ if ((a.exponent == 0 && a.mantissa == 0)) {
+ makeZero(a.sign);
+ return;
+ }
+ modInternal(a,0);
+ }
+ /**
+ * Calculates the logical AND of this Real and
+ * a.
+ * Replaces the contents of this Real with the result.
+ *
+ * Semantics of bitwise logical operations exactly mimic those of
+ * Java's bitwise integer operators. In these operations, the
+ * internal binary representation of the numbers are used. If the
+ * values represented by the operands are not mathematical
+ * integers, the fractional bits are also included in the operation.
+ *
+ *
Negative numbers are interpreted as two's-complement,
+ * generalized to real numbers: Negating the number inverts all
+ * bits, including an infinite number of 1-bits before the radix
+ * point and an infinite number of 1-bits after the radix point. The
+ * infinite number of 1-bits after the radix is rounded upwards
+ * producing an infinite number of 0-bits, until the first 0-bit is
+ * encountered which will be switched to a 1 (rounded or not, these
+ * two forms are mathematically equivalent). For example, the number
+ * "1" negated, becomes (in binary form)
+ * ...1111110.111111.... Rounding of the infinite
+ * number of 1's after the radix gives the number
+ * ...1111111.000000..., which is exactly the way we
+ * usually see "-1" as two's-complement.
+ *
+ *
This method calculates a negative value if and only
+ * if this and a are both negative.
+ *
+ *
+ * Equivalent int code: |
+ * this &= a;
+ * |
| Error bound: |
+ * 0 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 1.5
+ * |
+ *
+ * @param a the Real argument
+ */
+ public void and(Real a) {
+ if ((this.exponent < 0 && this.mantissa != 0) || (a.exponent < 0 && a.mantissa != 0)) {
+ makeNan();
+ return;
+ }
+ if ((this.exponent == 0 && this.mantissa == 0) || (a.exponent == 0 && a.mantissa == 0)) {
+ makeZero();
+ return;
+ }
+ if ((this.exponent < 0 && this.mantissa == 0) || (a.exponent < 0 && a.mantissa == 0)) {
+ if (!(this.exponent < 0 && this.mantissa == 0) && (this.sign!=0)) {
+ { this.mantissa = a.mantissa; this.exponent = a.exponent; this.sign = a.sign; };
+ } else if (!(a.exponent < 0 && a.mantissa == 0) && (a.sign!=0))
+ ; // ASSIGN(this,this)
+ else if ((this.exponent < 0 && this.mantissa == 0) && (a.exponent < 0 && a.mantissa == 0) &&
+ (this.sign!=0) && (a.sign!=0))
+ ; // makeInfinity(1)
+ else
+ makeZero();
+ return;
+ }
+ byte s;
+ int e;
+ long m;
+ if (exponent >= a.exponent) {
+ s = a.sign;
+ e = a.exponent;
+ m = a.mantissa;
+ } else {
+ s = sign;
+ e = exponent;
+ m = mantissa;
+ sign = a.sign;
+ exponent = a.exponent;
+ mantissa = a.mantissa;
+ }
+ int shift = exponent-e;
+ if (shift>=64) {
+ if (s == 0)
+ makeZero(sign);
+ return;
+ }
+ if (s != 0)
+ m = -m;
+ if ((this.sign!=0))
+ mantissa = -mantissa;
+ mantissa &= m>>shift;
+ sign = 0;
+ if (mantissa < 0) {
+ mantissa = -mantissa;
+ sign = 1;
+ }
+ normalize();
+ }
+ /**
+ * Calculates the logical OR of this Real and
+ * a.
+ * Replaces the contents of this Real with the result.
+ *
+ * See {@link #and(Real)} for an explanation of the
+ * interpretation of a Real in bitwise operations.
+ * This method calculates a negative value if and only
+ * if either this or a is negative.
+ *
+ *
+ * Equivalent int code: |
+ * this |= a;
+ * |
| Error bound: |
+ * 0 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 1.6
+ * |
+ *
+ * @param a the Real argument
+ */
+ public void or(Real a) {
+ if ((this.exponent < 0 && this.mantissa != 0) || (a.exponent < 0 && a.mantissa != 0)) {
+ makeNan();
+ return;
+ }
+ if ((this.exponent == 0 && this.mantissa == 0) || (a.exponent == 0 && a.mantissa == 0)) {
+ if ((this.exponent == 0 && this.mantissa == 0))
+ { this.mantissa = a.mantissa; this.exponent = a.exponent; this.sign = a.sign; };
+ return;
+ }
+ if ((this.exponent < 0 && this.mantissa == 0) || (a.exponent < 0 && a.mantissa == 0)) {
+ if (!(this.exponent < 0 && this.mantissa == 0) && (this.sign!=0))
+ ; // ASSIGN(this,this);
+ else if (!(a.exponent < 0 && a.mantissa == 0) && (a.sign!=0)) {
+ { this.mantissa = a.mantissa; this.exponent = a.exponent; this.sign = a.sign; };
+ } else
+ makeInfinity(sign | a.sign);
+ return;
+ }
+ byte s;
+ int e;
+ long m;
+ if (((this.sign!=0) && exponent <= a.exponent) ||
+ ((a.sign==0) && exponent >= a.exponent))
+ {
+ s = a.sign;
+ e = a.exponent;
+ m = a.mantissa;
+ } else {
+ s = sign;
+ e = exponent;
+ m = mantissa;
+ sign = a.sign;
+ exponent = a.exponent;
+ mantissa = a.mantissa;
+ }
+ int shift = exponent-e;
+ if (shift>=64 || shift<=-64)
+ return;
+ if (s != 0)
+ m = -m;
+ if ((this.sign!=0))
+ mantissa = -mantissa;
+ if (shift>=0)
+ mantissa |= m>>shift;
+ else
+ mantissa |= m<<(-shift);
+ sign = 0;
+ if (mantissa < 0) {
+ mantissa = -mantissa;
+ sign = 1;
+ }
+ normalize();
+ }
+ /**
+ * Calculates the logical XOR of this Real and
+ * a.
+ * Replaces the contents of this Real with the result.
+ *
+ * See {@link #and(Real)} for an explanation of the
+ * interpretation of a Real in bitwise operations.
+ * This method calculates a negative value if and only
+ * if exactly one of this and a is negative.
+ *
+ *
The operation NOT has been omitted in this library
+ * because it cannot be generalized to fractional numbers. If this
+ * Real represents a mathematical integer, the
+ * operation NOT can be calculated as "this XOR -1",
+ * which is equivalent to "this XOR
+ * /FFFFFFFF.0000".
+ *
+ *
+ * Equivalent int code: |
+ * this ^= a;
+ * |
| Error bound: |
+ * 0 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 1.5
+ * |
+ *
+ * @param a the Real argument
+ */
+ public void xor(Real a) {
+ if ((this.exponent < 0 && this.mantissa != 0) || (a.exponent < 0 && a.mantissa != 0)) {
+ makeNan();
+ return;
+ }
+ if ((this.exponent == 0 && this.mantissa == 0) || (a.exponent == 0 && a.mantissa == 0)) {
+ if ((this.exponent == 0 && this.mantissa == 0))
+ { this.mantissa = a.mantissa; this.exponent = a.exponent; this.sign = a.sign; };
+ return;
+ }
+ if ((this.exponent < 0 && this.mantissa == 0) || (a.exponent < 0 && a.mantissa == 0)) {
+ makeInfinity(sign ^ a.sign);
+ return;
+ }
+ byte s;
+ int e;
+ long m;
+ if (exponent >= a.exponent) {
+ s = a.sign;
+ e = a.exponent;
+ m = a.mantissa;
+ } else {
+ s = sign;
+ e = exponent;
+ m = mantissa;
+ sign = a.sign;
+ exponent = a.exponent;
+ mantissa = a.mantissa;
+ }
+ int shift = exponent-e;
+ if (shift>=64)
+ return;
+ if (s != 0)
+ m = -m;
+ if ((this.sign!=0))
+ mantissa = -mantissa;
+ mantissa ^= m>>shift;
+ sign = 0;
+ if (mantissa < 0) {
+ mantissa = -mantissa;
+ sign = 1;
+ }
+ normalize();
+ }
+ /**
+ * Calculates the value of this Real AND NOT
+ * a. The opeation is read as "bit clear".
+ * Replaces the contents of this Real with the result.
+ *
+ * See {@link #and(Real)} for an explanation of the
+ * interpretation of a Real in bitwise operations.
+ * This method calculates a negative value if and only
+ * if this is negative and not a is negative.
+ *
+ *
+ * Equivalent int code: |
+ * this &= ~a;
+ * |
| Error bound: |
+ * 0 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 1.5
+ * |
+ *
+ * @param a the Real argument
+ */
+ public void bic(Real a) {
+ if ((this.exponent < 0 && this.mantissa != 0) || (a.exponent < 0 && a.mantissa != 0)) {
+ makeNan();
+ return;
+ }
+ if ((this.exponent == 0 && this.mantissa == 0) || (a.exponent == 0 && a.mantissa == 0))
+ return;
+ if ((this.exponent < 0 && this.mantissa == 0) || (a.exponent < 0 && a.mantissa == 0)) {
+ if (!(this.exponent < 0 && this.mantissa == 0)) {
+ if ((this.sign!=0))
+ if ((a.sign!=0))
+ makeInfinity(0);
+ else
+ makeInfinity(1);
+ } else if ((a.sign!=0)) {
+ if ((a.exponent < 0 && a.mantissa == 0))
+ makeInfinity(0);
+ else
+ makeZero();
+ }
+ return;
+ }
+ int shift = exponent-a.exponent;
+ if (shift>=64 || (shift<=-64 && (this.sign==0)))
+ return;
+ long m = a.mantissa;
+ if ((a.sign!=0))
+ m = -m;
+ if ((this.sign!=0))
+ mantissa = -mantissa;
+ if (shift<0) {
+ if ((this.sign!=0)) {
+ if (shift<=-64)
+ mantissa = ~m;
+ else
+ mantissa = (mantissa>>(-shift)) & ~m;
+ exponent = a.exponent;
+ } else
+ mantissa &= ~(m<<(-shift));
+ } else
+ mantissa &= ~(m>>shift);
+ sign = 0;
+ if (mantissa < 0) {
+ mantissa = -mantissa;
+ sign = 1;
+ }
+ normalize();
+ }
+ private int compare(int a) {
+ tmp0.assign(a);
+ return compare(tmp0);
+ }
+ /**
+ * Calculates the square root of this Real.
+ * Replaces the contents of this Real with the result.
+ *
+ *
+ * Equivalent double code: |
+ * this = Math.{@link Math#sqrt(double) sqrt}(this);
+ * |
| Approximate error bound: |
+ * 1 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 19
+ * |
+ */
+ public void sqrt() {
+ /*
+ * Adapted from:
+ * Cephes Math Library Release 2.2: December, 1990
+ * Copyright 1984, 1990 by Stephen L. Moshier
+ *
+ * sqrtl.c
+ *
+ * long double sqrtl(long double x);
+ */
+ if ((this.exponent < 0 && this.mantissa != 0))
+ return;
+ if ((this.exponent == 0 && this.mantissa == 0)) {
+ sign=0;
+ return;
+ }
+ if ((this.sign!=0)) {
+ makeNan();
+ return;
+ }
+ if ((this.exponent < 0 && this.mantissa == 0))
+ return;
+ // Save X
+ { recipTmp.mantissa = this.mantissa; recipTmp.exponent = this.exponent; recipTmp.sign = this.sign; };
+ // normalize to range [0.5, 1)
+ int e = exponent-0x3fffffff;
+ exponent = 0x3fffffff;
+ // quadratic approximation, relative error 6.45e-4
+ { recipTmp2.mantissa = this.mantissa; recipTmp2.exponent = this.exponent; recipTmp2.sign = this.sign; };
+ { sqrtTmp.sign=(byte)1; sqrtTmp.exponent=0x3ffffffd; sqrtTmp.mantissa=0x68a7e193370ff21bL; };//-0.2044058315473477195990
+ mul(sqrtTmp);
+ { sqrtTmp.sign=(byte)0; sqrtTmp.exponent=0x3fffffff; sqrtTmp.mantissa=0x71f1e120690deae8L; };//0.89019407351052789754347
+ add(sqrtTmp);
+ mul(recipTmp2);
+ { sqrtTmp.sign=(byte)0; sqrtTmp.exponent=0x3ffffffe; sqrtTmp.mantissa=0x5045ee6baf28677aL; };//0.31356706742295303132394
+ add(sqrtTmp);
+ // adjust for odd powers of 2
+ if ((e&1) != 0)
+ mul(SQRT2);
+ // calculate exponent
+ exponent += e>>1;
+ // Newton iteratios:
+ // Yn+1 = (Yn + X/Yn)/2
+ for (int i=0; i<3; i++) {
+ { recipTmp2.mantissa = recipTmp.mantissa; recipTmp2.exponent = recipTmp.exponent; recipTmp2.sign = recipTmp.sign; };
+ recipTmp2.div(this);
+ add(recipTmp2);
+ scalbn(-1);
+ }
+ }
+ /**
+ * Calculates the reciprocal square root of this Real.
+ * Replaces the contents of this Real with the result.
+ *
+ *
+ * Equivalent double code: |
+ * this = 1/Math.{@link Math#sqrt(double) sqrt}(this);
+ * |
| Approximate error bound: |
+ * 1 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 21
+ * |
+ */
+ public void rsqrt() {
+ sqrt();
+ recip();
+ }
+ /**
+ * Calculates the cube root of this Real.
+ * Replaces the contents of this Real with the result.
+ * The cube root of a negative value is the negative of the cube
+ * root of that value's magnitude.
+ *
+ *
+ * Equivalent double code: |
+ * this = Math.{@link Math#cbrt(double) cbrt}(this);
+ * |
| Approximate error bound: |
+ * 2 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 32
+ * |
+ */
+ public void cbrt() {
+ if (!(this.exponent >= 0 && this.mantissa != 0))
+ return;
+ byte s = sign;
+ sign = 0;
+ // Calculates recipocal cube root of normalized Real,
+ // not zero, nan or infinity
+ final long start = 0x5120000000000000L;
+ // Save -A
+ { recipTmp.mantissa = this.mantissa; recipTmp.exponent = this.exponent; recipTmp.sign = this.sign; };
+ recipTmp.neg();
+ // First establish approximate result
+ mantissa = start-(mantissa>>>2);
+ int expRmd = exponent==0 ? 2 : (exponent-1)%3;
+ exponent = 0x40000000-(exponent-0x40000000-expRmd)/3;
+ normalize();
+ if (expRmd>0) {
+ { recipTmp2.sign=(byte)0; recipTmp2.exponent=0x3fffffff; recipTmp2.mantissa=0x6597fa94f5b8f20bL; }; // cbrt(1/2)
+ mul(recipTmp2);
+ if (expRmd>1)
+ mul(recipTmp2);
+ }
+ // Now perform Newton-Raphson iteration
+ // Xn+1 = (4*Xn - A*Xn**4)/3
+ for (int i=0; i<4; i++) {
+ { recipTmp2.mantissa = this.mantissa; recipTmp2.exponent = this.exponent; recipTmp2.sign = this.sign; };
+ sqr();
+ sqr();
+ mul(recipTmp);
+ recipTmp2.scalbn(2);
+ add(recipTmp2);
+ mul(THIRD);
+ }
+ recip();
+ if (!(this.exponent < 0 && this.mantissa != 0))
+ sign = s;
+ }
+ /**
+ * Calculates the n'th root of this Real.
+ * Replaces the contents of this Real with the result.
+ * For odd integer n, the n'th root of a negative value is the
+ * negative of the n'th root of that value's magnitude.
+ *
+ *
+ * Equivalent double code: |
+ * this = Math.{@link Math#pow(double,double)
+ * pow}(this,1/a);
+ * |
| Approximate error bound: |
+ * 2 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 110
+ * |
+ *
+ * @param n the Real argument.
+ */
+ public void nroot(Real n) {
+ if ((n.exponent < 0 && n.mantissa != 0)) {
+ makeNan();
+ return;
+ }
+ if (n.compare(THREE)==0) {
+ cbrt(); // Most probable application of nroot...
+ return;
+ } else if (n.compare(TWO)==0) {
+ sqrt(); // Also possible, should be optimized like this
+ return;
+ }
+ boolean negative = false;
+ if ((this.sign!=0) && n.isIntegral() && n.isOdd()) {
+ negative = true;
+ abs();
+ }
+ { tmp2.mantissa = n.mantissa; tmp2.exponent = n.exponent; tmp2.sign = n.sign; }; // Copy to temporary location in case of x.nroot(x)
+ tmp2.recip();
+ pow(tmp2);
+ if (negative)
+ neg();
+ }
+ /**
+ * Calculates sqrt(this*this+a*a).
+ * Replaces the contents of this Real with the result.
+ *
+ *
+ * Equivalent double code: |
+ * this = Math.{@link Math#hypot(double,double)
+ * hypot}(this,a);
+ * |
| Approximate error bound: |
+ * 1 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 24
+ * |
+ *
+ * @param a the Real argument.
+ */
+ public void hypot(Real a) {
+ { tmp1.mantissa = a.mantissa; tmp1.exponent = a.exponent; tmp1.sign = a.sign; }; // Copy to temporary location in case of x.hypot(x)
+ tmp1.sqr();
+ sqr();
+ add(tmp1);
+ sqrt();
+ }
+ private void exp2Internal(long extra) {
+ if ((this.exponent < 0 && this.mantissa != 0))
+ return;
+ if ((this.exponent < 0 && this.mantissa == 0)) {
+ if ((this.sign!=0))
+ makeZero(0);
+ return;
+ }
+ if ((this.exponent == 0 && this.mantissa == 0)) {
+ { this.mantissa = ONE.mantissa; this.exponent = ONE.exponent; this.sign = ONE.sign; };
+ return;
+ }
+ // Extract integer part
+ { expTmp.mantissa = this.mantissa; expTmp.exponent = this.exponent; expTmp.sign = this.sign; };
+ expTmp.add(HALF);
+ expTmp.floor();
+ int exp = expTmp.toInteger();
+ if (exp > 0x40000000) {
+ makeInfinity(sign);
+ return;
+ }
+ if (exp < -0x40000000) {
+ makeZero(sign);
+ return;
+ }
+ // Subtract integer part (this is where we need the extra accuracy)
+ expTmp.neg();
+ add128(extra,expTmp,0);
+ /*
+ * Adapted from:
+ * Cephes Math Library Release 2.7: May, 1998
+ * Copyright 1984, 1991, 1998 by Stephen L. Moshier
+ *
+ * exp2l.c
+ *
+ * long double exp2l(long double x);
+ */
+ // Now -0.5e raised to the power of this Real.
+ * Replaces the contents of this Real with the result.
+ *
+ *
+ * Equivalent double code: |
+ * this = Math.{@link Math#exp(double) exp}(this);
+ * |
| Approximate error bound: |
+ * 1 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 31
+ * |
+ */
+ public void exp() {
+ { expTmp.sign=(byte)0; expTmp.exponent=0x40000000; expTmp.mantissa=0x5c551d94ae0bf85dL; }; // log2(e)
+ long extra = mul128(0,expTmp,0xdf43ff68348e9f44L);
+ exp2Internal(extra);
+ }
+ /**
+ * Calculates 2 raised to the power of this Real.
+ * Replaces the contents of this Real with the result.
+ *
+ *
+ * Equivalent double code: |
+ * this = Math.{@link Math#exp(double) exp}(this *
+ * Math.{@link Math#log(double) log}(2));
+ * |
| Approximate error bound: |
+ * 1 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 27
+ * |
+ */
+ public void exp2() {
+ exp2Internal(0);
+ }
+ /**
+ * Calculates 10 raised to the power of this Real.
+ * Replaces the contents of this Real with the result.
+ *
+ *
+ * Equivalent double code: |
+ * this = Math.{@link Math#exp(double) exp}(this *
+ * Math.{@link Math#log(double) log}(10));
+ * |
| Approximate error bound: |
+ * 1 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 31
+ * |
+ */
+ public void exp10() {
+ { expTmp.sign=(byte)0; expTmp.exponent=0x40000001; expTmp.mantissa=0x6a4d3c25e68dc57fL; }; // log2(10)
+ long extra = mul128(0,expTmp,0x2495fb7fa6d7eda6L);
+ exp2Internal(extra);
+ }
+ private int lnInternal()
+ {
+ if ((this.exponent < 0 && this.mantissa != 0))
+ return 0;
+ if ((this.sign!=0)) {
+ makeNan();
+ return 0;
+ }
+ if ((this.exponent == 0 && this.mantissa == 0)) {
+ makeInfinity(1);
+ return 0;
+ }
+ if ((this.exponent < 0 && this.mantissa == 0))
+ return 0;
+ /*
+ * Adapted from:
+ * Cephes Math Library Release 2.7: May, 1998
+ * Copyright 1984, 1990, 1998 by Stephen L. Moshier
+ *
+ * logl.c
+ *
+ * long double logl(long double x);
+ */
+ // normalize to range [0.5, 1)
+ int e = exponent-0x3fffffff;
+ exponent = 0x3fffffff;
+ // rational appriximation
+ // log(1+x) = x - x²/2 + x³ P(x)/Q(x)
+ if (this.compare(SQRT1_2) < 0) {
+ e--;
+ exponent++;
+ }
+ sub(ONE);
+ { expTmp2.mantissa = this.mantissa; expTmp2.exponent = this.exponent; expTmp2.sign = this.sign; };
+ // P(x)
+ { this.sign=(byte)0; this.exponent=0x3ffffff1; this.mantissa=0x5ef0258ace5728ddL; };//4.5270000862445199635215E-5
+ mul(expTmp2);
+ { expTmp3.sign=(byte)0; expTmp3.exponent=0x3ffffffe; expTmp3.mantissa=0x7fa06283f86a0ce8L; };//0.4985410282319337597221
+ add(expTmp3);
+ mul(expTmp2);
+ { expTmp3.sign=(byte)0; expTmp3.exponent=0x40000002; expTmp3.mantissa=0x69427d1bd3e94ca1L; };//6.5787325942061044846969
+ add(expTmp3);
+ mul(expTmp2);
+ { expTmp3.sign=(byte)0; expTmp3.exponent=0x40000004; expTmp3.mantissa=0x77a5ce2e32e7256eL; };//29.911919328553073277375
+ add(expTmp3);
+ mul(expTmp2);
+ { expTmp3.sign=(byte)0; expTmp3.exponent=0x40000005; expTmp3.mantissa=0x79e63ae1b0cd4222L; };//60.949667980987787057556
+ add(expTmp3);
+ mul(expTmp2);
+ { expTmp3.sign=(byte)0; expTmp3.exponent=0x40000005; expTmp3.mantissa=0x7239d65d1e6840d6L; };//57.112963590585538103336
+ add(expTmp3);
+ mul(expTmp2);
+ { expTmp3.sign=(byte)0; expTmp3.exponent=0x40000004; expTmp3.mantissa=0x502880b6660c265fL; };//20.039553499201281259648
+ add(expTmp3);
+ // Q(x)
+ { expTmp.mantissa = expTmp2.mantissa; expTmp.exponent = expTmp2.exponent; expTmp.sign = expTmp2.sign; };
+ { expTmp3.sign=(byte)0; expTmp3.exponent=0x40000003; expTmp3.mantissa=0x7880d67a40f8dc5cL; };//15.062909083469192043167
+ expTmp.add(expTmp3);
+ expTmp.mul(expTmp2);
+ { expTmp3.sign=(byte)0; expTmp3.exponent=0x40000006; expTmp3.mantissa=0x530c2d4884d25e18L; };//83.047565967967209469434
+ expTmp.add(expTmp3);
+ expTmp.mul(expTmp2);
+ { expTmp3.sign=(byte)0; expTmp3.exponent=0x40000007; expTmp3.mantissa=0x6ee19643f3ed5776L; };//221.76239823732856465394
+ expTmp.add(expTmp3);
+ expTmp.mul(expTmp2);
+ { expTmp3.sign=(byte)0; expTmp3.exponent=0x40000008; expTmp3.mantissa=0x4d465177242295efL; };//309.09872225312059774938
+ expTmp.add(expTmp3);
+ expTmp.mul(expTmp2);
+ { expTmp3.sign=(byte)0; expTmp3.exponent=0x40000007; expTmp3.mantissa=0x6c36c4f923819890L; };//216.42788614495947685003
+ expTmp.add(expTmp3);
+ expTmp.mul(expTmp2);
+ { expTmp3.sign=(byte)0; expTmp3.exponent=0x40000005; expTmp3.mantissa=0x783cc111991239a3L; };//60.118660497603843919306
+ expTmp.add(expTmp3);
+ div(expTmp);
+ { expTmp3.mantissa = expTmp2.mantissa; expTmp3.exponent = expTmp2.exponent; expTmp3.sign = expTmp2.sign; };
+ expTmp3.sqr();
+ mul(expTmp3);
+ mul(expTmp2);
+ expTmp3.scalbn(-1);
+ sub(expTmp3);
+ add(expTmp2);
+ return e;
+ }
+ /**
+ * Calculates the natural logarithm (base-e) of this
+ * Real.
+ * Replaces the contents of this Real with the result.
+ *
+ *
+ * Equivalent double code: |
+ * this = Math.{@link Math#log(double) log}(this);
+ * |
| Approximate error bound: |
+ * 2 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 51
+ * |
+ */
+ public void ln() {
+ int exp = lnInternal();
+ expTmp.assign(exp);
+ expTmp.mul(LN2);
+ add(expTmp);
+ }
+ /**
+ * Calculates the base-2 logarithm of this Real.
+ * Replaces the contents of this Real with the result.
+ *
+ *
+ * Equivalent double code: |
+ * this = Math.{@link Math#log(double) log}(this)/Math.{@link
+ * Math#log(double) log}(2);
+ * |
| Approximate error bound: |
+ * 1 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 51
+ * |
+ */
+ public void log2() {
+ int exp = lnInternal();
+ mul(LOG2E);
+ add(exp);
+ }
+ /**
+ * Calculates the base-10 logarithm of this Real.
+ * Replaces the contents of this Real with the result.
+ *
+ *
+ * Equivalent double code: |
+ * this = Math.{@link Math#log10(double) log10}(this);
+ * |
| Approximate error bound: |
+ * 2 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 53
+ * |
+ */
+ public void log10() {
+ int exp = lnInternal();
+ expTmp.assign(exp);
+ expTmp.mul(LN2);
+ add(expTmp);
+ mul(LOG10E);
+ }
+ /**
+ * Calculates the closest power of 10 that is less than or equal to this
+ * Real.
+ * Replaces the contents of this Real with the result.
+ * The base-10 exponent of the result is returned.
+ *
+ *
+ * Equivalent double code: |
+ * int exp = (int)(Math.{@link Math#floor(double)
+ * floor}(Math.{@link Math#log10(double) log10}(this)));
+ *
this = Math.{@link Math#pow(double,double) pow}(10, exp);
+ * return exp;
+ * |
| Error bound: |
+ * ½ ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 3.6
+ * |
+ *
+ * @return the base-10 exponent
+ */
+ public int lowPow10() {
+ if (!(this.exponent >= 0 && this.mantissa != 0))
+ return 0;
+ { tmp2.mantissa = this.mantissa; tmp2.exponent = this.exponent; tmp2.sign = this.sign; };
+ // Approximate log10 using exponent only
+ int e = exponent - 0x40000000;
+ if (e<0) // it's important to achieve floor(exponent*ln2/ln10)
+ e = -(int)(((-e)*0x4d104d43L+((1L<<32)-1)) >> 32);
+ else
+ e = (int)(e*0x4d104d43L >> 32);
+ // Now, e < log10(this) < e+1
+ { this.mantissa = TEN.mantissa; this.exponent = TEN.exponent; this.sign = TEN.sign; };
+ pow(e);
+ if ((this.exponent == 0 && this.mantissa == 0)) { // A *really* small number, then
+ { tmp3.mantissa = TEN.mantissa; tmp3.exponent = TEN.exponent; tmp3.sign = TEN.sign; };
+ tmp3.pow(e+1);
+ } else {
+ { tmp3.mantissa = this.mantissa; tmp3.exponent = this.exponent; tmp3.sign = this.sign; };
+ tmp3.mul10();
+ }
+ if (tmp3.compare(tmp2) <= 0) {
+ // First estimate of log10 was too low
+ e++;
+ { this.mantissa = tmp3.mantissa; this.exponent = tmp3.exponent; this.sign = tmp3.sign; };
+ }
+ return e;
+ }
+ /**
+ * Calculates the value of this Real raised to the power of
+ * a.
+ * Replaces the contents of this Real with the result.
+ *
+ * Special cases:
+ *
+ * - if a is 0.0 or -0.0 then result is 1.0
+ *
- if a is NaN then result is NaN
+ *
- if this is NaN and a is not zero then result is NaN
+ *
- if a is 1.0 then result is this
+ *
- if |this| > 1.0 and a is +Infinity then result is +Infinity
+ *
- if |this| < 1.0 and a is -Infinity then result is +Infinity
+ *
- if |this| > 1.0 and a is -Infinity then result is +0
+ *
- if |this| < 1.0 and a is +Infinity then result is +0
+ *
- if |this| = 1.0 and a is ±Infinity then result is NaN
+ *
- if this = +0 and a > 0 then result is +0
+ *
- if this = +0 and a < 0 then result is +Inf
+ *
- if this = -0 and a > 0, and odd integer then result is -0
+ *
- if this = -0 and a < 0, and odd integer then result is -Inf
+ *
- if this = -0 and a > 0, not odd integer then result is +0
+ *
- if this = -0 and a < 0, not odd integer then result is +Inf
+ *
- if this = +Inf and a > 0 then result is +Inf
+ *
- if this = +Inf and a < 0 then result is +0
+ *
- if this = -Inf and a not integer then result is NaN
+ *
- if this = -Inf and a > 0, and odd integer then result is -Inf
+ *
- if this = -Inf and a > 0, not odd integer then result is +Inf
+ *
- if this = -Inf and a < 0, and odd integer then result is -0
+ *
- if this = -Inf and a < 0, not odd integer then result is +0
+ *
- if this < 0 and a not integer then result is NaN
+ *
- if this < 0 and a odd integer then result is -(|this|a)
+ *
- if this < 0 and a not odd integer then result is |this|a
+ *
- else result is exp(ln(this)*a)
+ *
+ *
+ *
+ * Equivalent double code: |
+ * this = Math.{@link Math#pow(double,double) pow}(this, a);
+ * |
| Approximate error bound: |
+ * 2 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 110
+ * |
+ *
+ * @param a the Real argument.
+ */
+ public void pow(Real a) {
+ if ((a.exponent == 0 && a.mantissa == 0)) {
+ { this.mantissa = ONE.mantissa; this.exponent = ONE.exponent; this.sign = ONE.sign; };
+ return;
+ }
+ if ((this.exponent < 0 && this.mantissa != 0) || (a.exponent < 0 && a.mantissa != 0)) {
+ makeNan();
+ return;
+ }
+ if (a.compare(ONE)==0)
+ return;
+ if ((a.exponent < 0 && a.mantissa == 0)) {
+ { tmp1.mantissa = this.mantissa; tmp1.exponent = this.exponent; tmp1.sign = this.sign; };
+ tmp1.abs();
+ int test = tmp1.compare(ONE);
+ if (test>0) {
+ if ((a.sign==0))
+ makeInfinity(0);
+ else
+ makeZero();
+ } else if (test<0) {
+ if ((a.sign!=0))
+ makeInfinity(0);
+ else
+ makeZero();
+ } else {
+ makeNan();
+ }
+ return;
+ }
+ if ((this.exponent == 0 && this.mantissa == 0)) {
+ if ((this.sign==0)) {
+ if ((a.sign==0))
+ makeZero();
+ else
+ makeInfinity(0);
+ } else {
+ if (a.isIntegral() && a.isOdd()) {
+ if ((a.sign==0))
+ makeZero(1);
+ else
+ makeInfinity(1);
+ } else {
+ if ((a.sign==0))
+ makeZero();
+ else
+ makeInfinity(0);
+ }
+ }
+ return;
+ }
+ if ((this.exponent < 0 && this.mantissa == 0)) {
+ if ((this.sign==0)) {
+ if ((a.sign==0))
+ makeInfinity(0);
+ else
+ makeZero();
+ } else {
+ if (a.isIntegral()) {
+ if (a.isOdd()) {
+ if ((a.sign==0))
+ makeInfinity(1);
+ else
+ makeZero(1);
+ } else {
+ if ((a.sign==0))
+ makeInfinity(0);
+ else
+ makeZero();
+ }
+ } else {
+ makeNan();
+ }
+ }
+ return;
+ }
+ if (a.isIntegral() && a.exponent <= 0x4000001e) {
+ pow(a.toInteger());
+ return;
+ }
+ byte s=0;
+ if ((this.sign!=0)) {
+ if (a.isIntegral()) {
+ if (a.isOdd())
+ s = 1;
+ } else {
+ makeNan();
+ return;
+ }
+ sign = 0;
+ }
+ { tmp1.mantissa = a.mantissa; tmp1.exponent = a.exponent; tmp1.sign = a.sign; };
+ if (tmp1.exponent <= 0x4000001e) {
+ // For increased accuracy, exponentiate with integer part of
+ // exponent by successive squaring
+ // (I really don't know why this works)
+ { tmp2.mantissa = tmp1.mantissa; tmp2.exponent = tmp1.exponent; tmp2.sign = tmp1.sign; };
+ tmp2.floor();
+ { tmp3.mantissa = this.mantissa; tmp3.exponent = this.exponent; tmp3.sign = this.sign; };
+ tmp3.pow(tmp2.toInteger());
+ tmp1.sub(tmp2);
+ } else {
+ { tmp3.mantissa = ONE.mantissa; tmp3.exponent = ONE.exponent; tmp3.sign = ONE.sign; };
+ }
+ // Do log2 and maintain accuracy
+ int e = lnInternal();
+ { tmp2.sign=(byte)0; tmp2.exponent=0x40000000; tmp2.mantissa=0x5c551d94ae0bf85dL; }; // log2(e)
+ long extra = mul128(0,tmp2,0xdf43ff68348e9f44L);
+ tmp2.assign(e);
+ extra = add128(extra,tmp2,0);
+ // Do exp2 of this multiplied by (fractional part of) exponent
+ extra = tmp1.mul128(0,this,extra);
+ tmp1.exp2Internal(extra);
+ { this.mantissa = tmp1.mantissa; this.exponent = tmp1.exponent; this.sign = tmp1.sign; };
+ mul(tmp3);
+ if (!(this.exponent < 0 && this.mantissa != 0))
+ sign = s;
+ }
+ /**
+ * Calculates the value of this Real raised to the power of
+ * the integer a.
+ * Replaces the contents of this Real with the result.
+ *
+ *
+ * Equivalent double code: |
+ * this = Math.{@link Math#pow(double,double) pow}(this, a);
+ * |
| Error bound: |
+ * ½ ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 84
+ * |
+ *
+ * @param a the integer argument.
+ */
+ public void pow(int a) {
+ // Calculate power of integer by successive squaring
+ boolean recp=false;
+ if (a < 0) {
+ a = -a; // Also works for 0x80000000
+ recp = true;
+ }
+ long extra = 0, expTmpExtra = 0;
+ { expTmp.mantissa = this.mantissa; expTmp.exponent = this.exponent; expTmp.sign = this.sign; };
+ { this.mantissa = ONE.mantissa; this.exponent = ONE.exponent; this.sign = ONE.sign; };
+ for (; a!=0; a>>>=1) {
+ if ((a & 1) != 0)
+ extra = mul128(extra,expTmp,expTmpExtra);
+ expTmpExtra = expTmp.mul128(expTmpExtra,expTmp,expTmpExtra);
+ }
+ if (recp)
+ extra = recip128(extra);
+ roundFrom128(extra);
+ }
+ private void sinInternal() {
+ /*
+ * Adapted from:
+ * Cephes Math Library Release 2.7: May, 1998
+ * Copyright 1985, 1990, 1998 by Stephen L. Moshier
+ *
+ * sinl.c
+ *
+ * long double sinl(long double x);
+ */
+ // XReal.
+ * Replaces the contents of this Real with the result.
+ * The input value is treated as an angle measured in radians.
+ *
+ *
+ * Equivalent double code: |
+ * this = Math.{@link Math#sin(double) sin}(this);
+ * |
| Approximate error bound: |
+ * 1 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 28
+ * |
+ */
+ public void sin() {
+ if (!(this.exponent >= 0 && this.mantissa != 0)) {
+ if (!(this.exponent == 0 && this.mantissa == 0))
+ makeNan();
+ return;
+ }
+ // Since sin(-x) = -sin(x) we can make sure that x > 0
+ boolean negative = false;
+ if ((this.sign!=0)) {
+ abs();
+ negative = true;
+ }
+ // Then reduce the argument to the range of 0 < x < pi*2
+ if (this.compare(PI2) > 0)
+ modInternal(PI2,0x62633145c06e0e69L);
+ // Since sin(pi*2 - x) = -sin(x) we can reduce the range 0 < x < pi
+ if (this.compare(PI) > 0) {
+ sub(PI2);
+ neg();
+ negative = !negative;
+ }
+ // Since sin(x) = sin(pi - x) we can reduce the range to 0 < x < pi/2
+ if (this.compare(PI_2) > 0) {
+ sub(PI);
+ neg();
+ }
+ // Since sin(x) = cos(pi/2 - x) we can reduce the range to 0 < x < pi/4
+ if (this.compare(PI_4) > 0) {
+ sub(PI_2);
+ neg();
+ cosInternal();
+ } else {
+ sinInternal();
+ }
+ if (negative)
+ neg();
+ if ((this.exponent == 0 && this.mantissa == 0))
+ abs(); // Remove confusing "-"
+ }
+ /**
+ * Calculates the trigonometric cosine of this Real.
+ * Replaces the contents of this Real with the result.
+ * The input value is treated as an angle measured in radians.
+ *
+ *
+ * Equivalent double code: |
+ * this = Math.{@link Math#cos(double) cos}(this);
+ * |
| Approximate error bound: |
+ * 1 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 37
+ * |
+ */
+ public void cos() {
+ if ((this.exponent == 0 && this.mantissa == 0)) {
+ { this.mantissa = ONE.mantissa; this.exponent = ONE.exponent; this.sign = ONE.sign; };
+ return;
+ }
+ if ((this.sign!=0))
+ abs();
+ if (this.compare(PI_4) < 0) {
+ cosInternal();
+ } else {
+ add(PI_2);
+ sin();
+ }
+ }
+ /**
+ * Calculates the trigonometric tangent of this Real.
+ * Replaces the contents of this Real with the result.
+ * The input value is treated as an angle measured in radians.
+ *
+ *
+ * Equivalent double code: |
+ * this = Math.{@link Math#tan(double) tan}(this);
+ * |
| Approximate error bound: |
+ * 2 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 70
+ * |
+ */
+ public void tan() {
+ { tmp4.mantissa = this.mantissa; tmp4.exponent = this.exponent; tmp4.sign = this.sign; };
+ tmp4.cos();
+ sin();
+ div(tmp4);
+ }
+ /**
+ * Calculates the trigonometric arc sine of this Real,
+ * in the range -π/2 to π/2.
+ * Replaces the contents of this Real with the result.
+ *
+ *
+ * Equivalent double code: |
+ * this = Math.{@link Math#asin(double) asin}(this);
+ * |
| Approximate error bound: |
+ * 3 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 68
+ * |
+ */
+ public void asin() {
+ { tmp1.mantissa = this.mantissa; tmp1.exponent = this.exponent; tmp1.sign = this.sign; };
+ sqr();
+ neg();
+ add(ONE);
+ rsqrt();
+ mul(tmp1);
+ atan();
+ }
+ /**
+ * Calculates the trigonometric arc cosine of this Real,
+ * in the range 0.0 to π.
+ * Replaces the contents of this Real with the result.
+ *
+ *
+ * Equivalent double code: |
+ * this = Math.{@link Math#acos(double) acos}(this);
+ * |
| Approximate error bound: |
+ * 2 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 67
+ * |
+ */
+ public void acos() {
+ boolean negative = (this.sign!=0);
+ abs();
+ { tmp1.mantissa = this.mantissa; tmp1.exponent = this.exponent; tmp1.sign = this.sign; };
+ sqr();
+ neg();
+ add(ONE);
+ sqrt();
+ div(tmp1);
+ atan();
+ if (negative) {
+ neg();
+ add(PI);
+ }
+ }
+ /**
+ * Calculates the trigonometric arc tangent of this Real,
+ * in the range -π/2 to π/2.
+ * Replaces the contents of this Real with the result.
+ *
+ *
+ * Equivalent double code: |
+ * this = Math.{@link Math#atan(double) atan}(this);
+ * |
| Approximate error bound: |
+ * 2 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 37
+ * |
+ */
+ public void atan() {
+ /*
+ * Adapted from:
+ * Cephes Math Library Release 2.7: May, 1998
+ * Copyright 1984, 1990, 1998 by Stephen L. Moshier
+ *
+ * atanl.c
+ *
+ * long double atanl(long double x);
+ */
+ if ((this.exponent == 0 && this.mantissa == 0) || (this.exponent < 0 && this.mantissa != 0))
+ return;
+ if ((this.exponent < 0 && this.mantissa == 0)) {
+ byte s = sign;
+ { this.mantissa = PI_2.mantissa; this.exponent = PI_2.exponent; this.sign = PI_2.sign; };
+ sign = s;
+ return;
+ }
+ byte s = sign;
+ sign = 0;
+ // range reduction
+ boolean addPI_2 = false;
+ boolean addPI_4 = false;
+ { tmp1.mantissa = SQRT2.mantissa; tmp1.exponent = SQRT2.exponent; tmp1.sign = SQRT2.sign; };
+ tmp1.add(ONE);
+ if (this.compare(tmp1) > 0) {
+ addPI_2 = true;
+ recip();
+ neg();
+ } else {
+ tmp1.sub(TWO);
+ if (this.compare(tmp1) > 0) {
+ addPI_4 = true;
+ { tmp1.mantissa = this.mantissa; tmp1.exponent = this.exponent; tmp1.sign = this.sign; };
+ tmp1.add(ONE);
+ sub(ONE);
+ div(tmp1);
+ }
+ }
+ // Now |X|Real divided by x, in the range -π
+ * to π. The signs of both arguments are used to determine the
+ * quadrant of the result. Replaces the contents of this
+ * Real with the result.
+ *
+ *
+ * Equivalent double code: |
+ * this = Math.{@link Math#atan2(double,double)
+ * atan2}(this,x);
+ * |
| Approximate error bound: |
+ * 2 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 48
+ * |
+ *
+ * @param x the Real argument.
+ */
+ public void atan2(Real x) {
+ if ((this.exponent < 0 && this.mantissa != 0) || (x.exponent < 0 && x.mantissa != 0) || ((this.exponent < 0 && this.mantissa == 0) && (x.exponent < 0 && x.mantissa == 0))) {
+ makeNan();
+ return;
+ }
+ if ((this.exponent == 0 && this.mantissa == 0) && (x.exponent == 0 && x.mantissa == 0))
+ return;
+ byte s = sign;
+ byte s2 = x.sign;
+ sign = 0;
+ x.sign = 0;
+ div(x);
+ atan();
+ if (s2 != 0) {
+ neg();
+ add(PI);
+ }
+ sign = s;
+ }
+ /**
+ * Calculates the hyperbolic sine of this Real.
+ * Replaces the contents of this Real with the result.
+ *
+ *
+ * Equivalent double code: |
+ * this = Math.{@link Math#sinh(double) sinh}(this);
+ * |
| Approximate error bound: |
+ * 2 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 67
+ * |
+ */
+ public void sinh() {
+ { tmp1.mantissa = this.mantissa; tmp1.exponent = this.exponent; tmp1.sign = this.sign; };
+ tmp1.neg();
+ tmp1.exp();
+ exp();
+ sub(tmp1);
+ scalbn(-1);
+ }
+ /**
+ * Calculates the hyperbolic cosine of this Real.
+ * Replaces the contents of this Real with the result.
+ *
+ *
+ * Equivalent double code: |
+ * this = Math.{@link Math#cosh(double) cosh}(this);
+ * |
| Approximate error bound: |
+ * 2 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 66
+ * |
+ */
+ public void cosh() {
+ { tmp1.mantissa = this.mantissa; tmp1.exponent = this.exponent; tmp1.sign = this.sign; };
+ tmp1.neg();
+ tmp1.exp();
+ exp();
+ add(tmp1);
+ scalbn(-1);
+ }
+ /**
+ * Calculates the hyperbolic tangent of this Real.
+ * Replaces the contents of this Real with the result.
+ *
+ *
+ * Equivalent double code: |
+ * this = Math.{@link Math#tanh(double) tanh}(this);
+ * |
| Approximate error bound: |
+ * 2 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 70
+ * |
+ */
+ public void tanh() {
+ { tmp1.mantissa = this.mantissa; tmp1.exponent = this.exponent; tmp1.sign = this.sign; };
+ tmp1.neg();
+ tmp1.exp();
+ exp();
+ { tmp2.mantissa = this.mantissa; tmp2.exponent = this.exponent; tmp2.sign = this.sign; };
+ tmp2.add(tmp1);
+ sub(tmp1);
+ div(tmp2);
+ }
+ /**
+ * Calculates the hyperbolic arc sine of this Real.
+ * Replaces the contents of this Real with the result.
+ *
+ *
+ * Equivalent double code: |
+ * none
+ * |
| Approximate error bound: |
+ * 2 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 77
+ * |
+ */
+ public void asinh() {
+ if (!(this.exponent >= 0 && this.mantissa != 0))
+ return;
+ // Use symmetry to prevent underflow error for very large negative
+ // values
+ byte s = sign;
+ sign = 0;
+ { tmp1.mantissa = this.mantissa; tmp1.exponent = this.exponent; tmp1.sign = this.sign; };
+ tmp1.sqr();
+ tmp1.add(ONE);
+ tmp1.sqrt();
+ add(tmp1);
+ ln();
+ if (!(this.exponent < 0 && this.mantissa != 0))
+ sign = s;
+ }
+ /**
+ * Calculates the hyperbolic arc cosine of this Real.
+ * Replaces the contents of this Real with the result.
+ *
+ *
+ * Equivalent double code: |
+ * none
+ * |
| Approximate error bound: |
+ * 2 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 75
+ * |
+ */
+ public void acosh() {
+ { tmp1.mantissa = this.mantissa; tmp1.exponent = this.exponent; tmp1.sign = this.sign; };
+ tmp1.sqr();
+ tmp1.sub(ONE);
+ tmp1.sqrt();
+ add(tmp1);
+ ln();
+ }
+ /**
+ * Calculates the hyperbolic arc tangent of this Real.
+ * Replaces the contents of this Real with the result.
+ *
+ *
+ * Equivalent double code: |
+ * none
+ * |
| Approximate error bound: |
+ * 2 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 57
+ * |
+ */
+ public void atanh() {
+ { tmp1.mantissa = this.mantissa; tmp1.exponent = this.exponent; tmp1.sign = this.sign; };
+ tmp1.neg();
+ tmp1.add(ONE);
+ add(ONE);
+ div(tmp1);
+ ln();
+ scalbn(-1);
+ }
+ /**
+ * Calculates the factorial of this Real.
+ * Replaces the contents of this Real with the result.
+ * The definition is generalized to all real numbers (not only integers),
+ * by using the fact that (n!)={@link #gamma() gamma}(n+1).
+ *
+ *
+ * Equivalent double code: |
+ * none
+ * |
| Approximate error bound: |
+ * 15 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 8-190
+ * |
+ */
+ public void fact() {
+ if (!(this.exponent >= 0))
+ return;
+ if (!this.isIntegral() || this.compare(ZERO)<0 || this.compare(200)>0)
+ {
+ // x<0, x>200 or not integer: fact(x) = gamma(x+1)
+ add(ONE);
+ gamma();
+ return;
+ }
+ { tmp1.mantissa = this.mantissa; tmp1.exponent = this.exponent; tmp1.sign = this.sign; };
+ { this.mantissa = ONE.mantissa; this.exponent = ONE.exponent; this.sign = ONE.sign; };
+ while (tmp1.compare(ONE) > 0) {
+ mul(tmp1);
+ tmp1.sub(ONE);
+ }
+ }
+ /**
+ * Calculates the gamma function for this Real.
+ * Replaces the contents of this Real with the result.
+ *
+ *
+ * Equivalent double code: |
+ * none
+ * |
| Approximate error bound: |
+ * 100+ ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 190
+ * |
+ */
+ public void gamma() {
+ if (!(this.exponent >= 0))
+ return;
+ // x<0: gamma(-x) = -pi/(x*gamma(x)*sin(pi*x))
+ boolean negative = (this.sign!=0);
+ abs();
+ { tmp1.mantissa = this.mantissa; tmp1.exponent = this.exponent; tmp1.sign = this.sign; };
+ // xn: gamma(x) = exp((x-1/2)*ln(x) - x + ln(2*pi)/2 + 1/12x - 1/360x³
+ // + 1/1260x**5 - 1/1680x**7+1/1188x**9)
+ { tmp3.mantissa = this.mantissa; tmp3.exponent = this.exponent; tmp3.sign = this.sign; }; // x
+ { tmp4.mantissa = this.mantissa; tmp4.exponent = this.exponent; tmp4.sign = this.sign; };
+ tmp4.sqr(); // x²
+ // (x-1/2)*ln(x)-x
+ ln(); { tmp5.mantissa = tmp3.mantissa; tmp5.exponent = tmp3.exponent; tmp5.sign = tmp3.sign; }; tmp5.sub(HALF); mul(tmp5); sub(tmp3);
+ // + ln(2*pi)/2
+ { tmp5.sign=(byte)0; tmp5.exponent=0x3fffffff; tmp5.mantissa=0x759fc72192fad29aL; }; add(tmp5);
+ // + 1/12x
+ tmp5.assign( 12); tmp5.mul(tmp3); tmp5.recip(); add(tmp5); tmp3.mul(tmp4);
+ // - 1/360x³
+ tmp5.assign( 360); tmp5.mul(tmp3); tmp5.recip(); sub(tmp5); tmp3.mul(tmp4);
+ // + 1/1260x**5
+ tmp5.assign(1260); tmp5.mul(tmp3); tmp5.recip(); add(tmp5); tmp3.mul(tmp4);
+ // - 1/1680x**7
+ tmp5.assign(1680); tmp5.mul(tmp3); tmp5.recip(); sub(tmp5); tmp3.mul(tmp4);
+ // + 1/1188x**9
+ tmp5.assign(1188); tmp5.mul(tmp3); tmp5.recip(); add(tmp5);
+ exp();
+ if (divide)
+ div(tmp2);
+ if (negative) {
+ { tmp5.mantissa = tmp1.mantissa; tmp5.exponent = tmp1.exponent; tmp5.sign = tmp1.sign; }; // sin() uses tmp1
+ // -pi/(x*gamma(x)*sin(pi*x))
+ mul(tmp5);
+ tmp5.scalbn(-1); tmp5.frac(); tmp5.mul(PI2); // Fixes integer inaccuracy
+ tmp5.sin(); mul(tmp5); recip(); mul(PI); neg();
+ }
+ }
+ private void erfc1Internal() {
+ // 3 5 7 9
+ // 2 / x x x x // erfc(x) = 1 - ------ | x - --- + ---- - ---- + ---- - ... |
+ // sqrt(pi)\ 3 2!*5 3!*7 4!*9 /
+ //
+ long extra=0,tmp1Extra,tmp2Extra,tmp3Extra,tmp4Extra;
+ { tmp1.mantissa = this.mantissa; tmp1.exponent = this.exponent; tmp1.sign = this.sign; }; tmp1Extra = 0;
+ { tmp2.mantissa = this.mantissa; tmp2.exponent = this.exponent; tmp2.sign = this.sign; };
+ tmp2Extra = tmp2.mul128(0,tmp2,0);
+ tmp2.neg();
+ { tmp3.mantissa = ONE.mantissa; tmp3.exponent = ONE.exponent; tmp3.sign = ONE.sign; }; tmp3Extra = 0;
+ int i=1;
+ do {
+ tmp1Extra = tmp1.mul128(tmp1Extra,tmp2,tmp2Extra);
+ tmp4.assign(i);
+ tmp3Extra = tmp3.mul128(tmp3Extra,tmp4,0);
+ tmp4.assign(2*i+1);
+ tmp4Extra = tmp4.mul128(0,tmp3,tmp3Extra);
+ tmp4Extra = tmp4.recip128(tmp4Extra);
+ tmp4Extra = tmp4.mul128(tmp4Extra,tmp1,tmp1Extra);
+ extra = add128(extra,tmp4,tmp4Extra);
+ i++;
+ } while (exponent - tmp4.exponent < 128);
+ { tmp1.sign=(byte)1; tmp1.exponent=0x40000000; tmp1.mantissa=0x48375d410a6db446L; }; // -2/sqrt(pi)
+ extra = mul128(extra,tmp1,0xb8ea453fb5ff61a2L);
+ extra = add128(extra,ONE,0);
+ roundFrom128(extra);
+ }
+ private void erfc2Internal() {
+ // -x² -1
+ // e x / 1 3 3*5 3*5*7 // erfc(x) = -------- | 1 - --- + ------ - ------ + ------ - ... |
+ // sqrt(pi) \ 2x² 2 3 4 /
+ // (2x²) (2x²) (2x²)
+ // Calculate iteration stop criteria
+ { tmp1.mantissa = this.mantissa; tmp1.exponent = this.exponent; tmp1.sign = this.sign; };
+ tmp1.sqr();
+ { tmp2.sign=(byte)0; tmp2.exponent=0x40000000; tmp2.mantissa=0x5c3811b4bfd0c8abL; }; // 1/0.694
+ tmp2.mul(tmp1);
+ tmp2.sub(HALF);
+ int digits = tmp2.toInteger(); // number of accurate digits = x*x/0.694-0.5
+ if (digits > 64)
+ digits = 64;
+ tmp1.scalbn(1);
+ int dxq = tmp1.toInteger()+1;
+ { tmp1.mantissa = this.mantissa; tmp1.exponent = this.exponent; tmp1.sign = this.sign; };
+ recip();
+ { tmp2.mantissa = this.mantissa; tmp2.exponent = this.exponent; tmp2.sign = this.sign; };
+ { tmp3.mantissa = this.mantissa; tmp3.exponent = this.exponent; tmp3.sign = this.sign; };
+ tmp3.sqr();
+ tmp3.neg();
+ tmp3.scalbn(-1);
+ { this.mantissa = ONE.mantissa; this.exponent = ONE.exponent; this.sign = ONE.sign; };
+ { tmp4.mantissa = ONE.mantissa; tmp4.exponent = ONE.exponent; tmp4.sign = ONE.sign; };
+ int i=1;
+ do {
+ tmp4.mul(2*i-1);
+ tmp4.mul(tmp3);
+ add(tmp4);
+ i++;
+ } while (tmp4.exponent-0x40000000>-(digits+2) && 2*i-1Real.
+ * Replaces the contents of this Real with the result.
+ *
+ * The complementary error function is defined as the integral from
+ * x to infinity of 2/√π ·e-t² dt. It is
+ * related to the error function, erf, by the formula
+ * erfc(x)=1-erf(x).
+ *
+ *
+ * Equivalent double code: |
+ * none
+ * |
| Approximate error bound: |
+ * 219 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 80-4900
+ * |
+ */
+ public void erfc() {
+ if ((this.exponent < 0 && this.mantissa != 0))
+ return;
+ if ((this.exponent == 0 && this.mantissa == 0)) {
+ { this.mantissa = ONE.mantissa; this.exponent = ONE.exponent; this.sign = ONE.sign; };
+ return;
+ }
+ if ((this.exponent < 0 && this.mantissa == 0) || toInteger()>27281) {
+ if ((this.sign!=0)) {
+ { this.mantissa = TWO.mantissa; this.exponent = TWO.exponent; this.sign = TWO.sign; };
+ } else
+ makeZero(0);
+ return;
+ }
+ byte s = sign;
+ sign = 0;
+ { tmp1.sign=(byte)0; tmp1.exponent=0x40000002; tmp1.mantissa=0x570a3d70a3d70a3dL; }; // 5.44
+ if (this.lessThan(tmp1))
+ erfc1Internal();
+ else
+ erfc2Internal();
+ if (s != 0) {
+ neg();
+ add(TWO);
+ }
+ }
+ /**
+ * Calculates the inverse complementary error function for this
+ * Real.
+ * Replaces the contents of this Real with the result.
+ *
+ *
+ * Equivalent double code: |
+ * none
+ * |
| Approximate error bound: |
+ * 219 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 240-5100
+ * |
+ */
+ public void inverfc() {
+ if ((this.exponent < 0 && this.mantissa != 0) || (this.sign!=0) || this.greaterThan(TWO)) {
+ makeNan();
+ return;
+ }
+ if ((this.exponent == 0 && this.mantissa == 0)) {
+ makeInfinity(0);
+ return;
+ }
+ if (this.equalTo(TWO)) {
+ makeInfinity(1);
+ return;
+ }
+ int sign = ONE.compare(this);
+ if (sign==0) {
+ makeZero();
+ return;
+ }
+ if (sign<0) {
+ neg();
+ add(TWO);
+ }
+ // Using invphi to calculate inverfc, like this
+ // inverfc(x) = -invphi(x/2)/(sqrt(2))
+ scalbn(-1);
+ // Inverse Phi Algorithm (phi(Z)=P, so invphi(P)=Z)
+ // ------------------------------------------------
+ // Part 1: Numerical Approximation Method for Inverse Phi
+ // This accepts input of P and outputs approximate Z as Y
+ // Source:Odeh & Evans. 1974. AS 70. Applied Statistics.
+ // R = sqrt(Ln(1/(Q²)))
+ { tmp1.mantissa = this.mantissa; tmp1.exponent = this.exponent; tmp1.sign = this.sign; };
+ tmp1.ln();
+ tmp1.mul(-2);
+ tmp1.sqrt();
+ // Y = -(R+((((P4*R+P3)*R+P2)*R+P1)*R+P0)/((((Q4*R+Q3)*R*Q2)*R+Q1)*R+Q0))
+ { tmp2.sign=(byte)1; tmp2.exponent=0x3ffffff1; tmp2.mantissa=0x5f22bb0fb4698674L; }; // P4=-0.0000453642210148
+ tmp2.mul(tmp1);
+ { tmp3.sign=(byte)1; tmp3.exponent=0x3ffffffa; tmp3.mantissa=0x53a731ce1ea0be15L; }; // P3=-0.0204231210245
+ tmp2.add(tmp3);
+ tmp2.mul(tmp1);
+ { tmp3.sign=(byte)1; tmp3.exponent=0x3ffffffe; tmp3.mantissa=0x579d2d719fc517f3L; }; // P2=-0.342242088547
+ tmp2.add(tmp3);
+ tmp2.mul(tmp1);
+ tmp2.add(-1); // P1=-1
+ tmp2.mul(tmp1);
+ { tmp3.sign=(byte)1; tmp3.exponent=0x3ffffffe; tmp3.mantissa=0x527dd3193bc8dd4cL; }; // P0=-0.322232431088
+ tmp2.add(tmp3);
+ { tmp3.sign=(byte)0; tmp3.exponent=0x3ffffff7; tmp3.mantissa=0x7e5b0f681d161e7dL; }; // Q4=0.0038560700634
+ tmp3.mul(tmp1);
+ { tmp4.sign=(byte)0; tmp4.exponent=0x3ffffffc; tmp4.mantissa=0x6a05ccf9917da0a8L; }; // Q3=0.103537752850
+ tmp3.add(tmp4);
+ tmp3.mul(tmp1);
+ { tmp4.sign=(byte)0; tmp4.exponent=0x3fffffff; tmp4.mantissa=0x43fb32c0d3c14ec4L; }; // Q2=0.531103462366
+ tmp3.add(tmp4);
+ tmp3.mul(tmp1);
+ { tmp4.sign=(byte)0; tmp4.exponent=0x3fffffff; tmp4.mantissa=0x4b56a41226f4ba95L; }; // Q1=0.588581570495
+ tmp3.add(tmp4);
+ tmp3.mul(tmp1);
+ { tmp4.sign=(byte)0; tmp4.exponent=0x3ffffffc; tmp4.mantissa=0x65bb9a7733dd5062L; }; // Q0=0.0993484626060
+ tmp3.add(tmp4);
+ tmp2.div(tmp3);
+ tmp1.add(tmp2);
+ tmp1.neg();
+ { sqrtTmp.mantissa = tmp1.mantissa; sqrtTmp.exponent = tmp1.exponent; sqrtTmp.sign = tmp1.sign; }; // sqrtTmp and tmp5 not used by erfc() and exp()
+ // Part 2: Refine to accuracy of erfc Function
+ // This accepts inputs Y and P (from above) and outputs Z
+ // (Using Halley's third order method for finding roots of equations)
+ // Q = erfc(-Y/sqrt(2))/2-P
+ { tmp5.mantissa = sqrtTmp.mantissa; tmp5.exponent = sqrtTmp.exponent; tmp5.sign = sqrtTmp.sign; };
+ tmp5.mul(SQRT1_2);
+ tmp5.neg();
+ tmp5.erfc();
+ tmp5.scalbn(-1);
+ tmp5.sub(this);
+ // R = Q*sqrt(2*pi)*e^(Y²/2)
+ { tmp3.mantissa = sqrtTmp.mantissa; tmp3.exponent = sqrtTmp.exponent; tmp3.sign = sqrtTmp.sign; };
+ tmp3.sqr();
+ tmp3.scalbn(-1);
+ tmp3.exp();
+ tmp5.mul(tmp3);
+ { tmp3.sign=(byte)0; tmp3.exponent=0x40000001; tmp3.mantissa=0x50364c7fd89c1659L; }; // sqrt(2*pi)
+ tmp5.mul(tmp3);
+ // Z = Y-R/(1+R*Y/2)
+ { this.mantissa = sqrtTmp.mantissa; this.exponent = sqrtTmp.exponent; this.sign = sqrtTmp.sign; };
+ mul(tmp5);
+ scalbn(-1);
+ add(ONE);
+ rdiv(tmp5);
+ neg();
+ add(sqrtTmp);
+ // calculate inverfc(x) = -invphi(x/2)/(sqrt(2))
+ mul(SQRT1_2);
+ if (sign>0)
+ neg();
+ }
+ //*************************************************************************
+ // Calendar conversions taken from
+ // http://www.fourmilab.ch/documents/calendar/
+ private static int floorDiv(int a, int b) {
+ if (a>=0)
+ return a/b;
+ return -((-a+b-1)/b);
+ }
+ private static int floorMod(int a, int b) {
+ if (a>=0)
+ return a%b;
+ return a+((-a+b-1)/b)*b;
+ }
+ private static boolean leap_gregorian(int year) {
+ return ((year % 4) == 0) &&
+ (!(((year % 100) == 0) && ((year % 400) != 0)));
+ }
+ // GREGORIAN_TO_JD -- Determine Julian day number from Gregorian
+ // calendar date -- Except that we use 1/1-0 as day 0
+ private static int gregorian_to_jd(int year, int month, int day) {
+ return ((366 - 1) +
+ (365 * (year - 1)) +
+ (floorDiv(year - 1, 4)) +
+ (-floorDiv(year - 1, 100)) +
+ (floorDiv(year - 1, 400)) +
+ ((((367 * month) - 362) / 12) +
+ ((month <= 2) ? 0 : (leap_gregorian(year) ? -1 : -2)) + day));
+ }
+ // JD_TO_GREGORIAN -- Calculate Gregorian calendar date from Julian
+ // day -- Except that we use 1/1-0 as day 0
+ private static int jd_to_gregorian(int jd) {
+ int wjd, depoch, quadricent, dqc, cent, dcent, quad, dquad,
+ yindex, year, yearday, leapadj, month, day;
+ wjd = jd;
+ depoch = wjd - 366;
+ quadricent = floorDiv(depoch, 146097);
+ dqc = floorMod(depoch, 146097);
+ cent = floorDiv(dqc, 36524);
+ dcent = floorMod(dqc, 36524);
+ quad = floorDiv(dcent, 1461);
+ dquad = floorMod(dcent, 1461);
+ yindex = floorDiv(dquad, 365);
+ year = (quadricent * 400) + (cent * 100) + (quad * 4) + yindex;
+ if (!((cent == 4) || (yindex == 4)))
+ year++;
+ yearday = wjd - gregorian_to_jd(year, 1, 1);
+ leapadj = ((wjd < gregorian_to_jd(year, 3, 1)) ? 0
+ : (leap_gregorian(year) ? 1 : 2));
+ month = floorDiv(((yearday + leapadj) * 12) + 373, 367);
+ day = (wjd - gregorian_to_jd(year, month, 1)) + 1;
+ return (year*100+month)*100+day;
+ }
+ /**
+ * Converts this Real from "hours" to "days, hours,
+ * minutes and seconds".
+ * Replaces the contents of this Real with the result.
+ *
+ * The format converted to is encoded into the digits of the
+ * number (in decimal form):
+ * "DDDDhh.mmss". Here "DDDD," is number
+ * of days, "hh" is hours (0-23), "mm" is
+ * minutes (0-59) and "ss" is seconds
+ * (0-59). Additional digits represent fractions of a second.
+ *
+ *
If the number of hours of the input is greater or equal to
+ * 8784 (number of hours in year 0), the format
+ * converted to is instead "YYYYMMDDhh.mmss". Here
+ * "YYYY" is the number of years since the imaginary
+ * year 0 in the Gregorian calendar, extrapolated back
+ * from year 1582. "MM" is the month (1-12) and
+ * "DD" is the day of the month (1-31). See a thorough
+ * discussion of date calculations here.
+ *
+ *
+ * Equivalent double code: |
+ * none
+ * |
| Approximate error bound: |
+ * ?
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 19
+ * |
+ */
+ public void toDHMS() {
+ if (!(this.exponent >= 0 && this.mantissa != 0))
+ return;
+ boolean negative = (this.sign!=0);
+ abs();
+ int D,m;
+ long h;
+ h = toLong();
+ frac();
+ tmp1.assign(60);
+ mul(tmp1);
+ m = toInteger();
+ frac();
+ mul(tmp1);
+ // MAGIC ROUNDING: Check if we are 2**-16 sec short of a whole minute
+ // i.e. "seconds" > 59.999985
+ { tmp2.mantissa = ONE.mantissa; tmp2.exponent = ONE.exponent; tmp2.sign = ONE.sign; };
+ tmp2.scalbn(-16);
+ add(tmp2);
+ if (this.compare(tmp1) >= 0) {
+ // Yes. So set zero secs instead and carry over to mins and hours
+ { this.mantissa = ZERO.mantissa; this.exponent = ZERO.exponent; this.sign = ZERO.sign; };
+ m++;
+ if (m >= 60) {
+ m -= 60;
+ h++;
+ }
+ // Phew! That was close. From now on it is integer arithmetic...
+ } else {
+ // Nope. So try to undo the damage...
+ sub(tmp2);
+ }
+ D = (int)(h/24);
+ h %= 24;
+ if (D >= 366)
+ D = jd_to_gregorian(D);
+ add(m*100);
+ div(10000);
+ tmp1.assign(D*100L+h);
+ add(tmp1);
+ if (negative)
+ neg();
+ }
+ /**
+ * Converts this Real from "days, hours, minutes and
+ * seconds" to "hours".
+ * Replaces the contents of this Real with the result.
+ *
+ * The format converted from is encoded into the digits of the
+ * number (in decimal form):
+ * "DDDDhh.mmss". Here "DDDD" is number of
+ * days, "hh" is hours (0-23), "mm" is
+ * minutes (0-59) and "ss" is seconds
+ * (0-59). Additional digits represent fractions of a second.
+ *
+ *
If the number of days in the input is greater than or equal to
+ * 10000, the format converted from is instead
+ * "YYYYMMDDhh.mmss". Here "YYYY" is the
+ * number of years since the imaginary year 0 in the
+ * Gregorian calendar, extrapolated back from year
+ * 1582. "MM" is the month (1-12) and
+ * "DD" is the day of the month (1-31). If month or day
+ * is 0 it is treated as 1. See a thorough discussion of date
+ * calculations here.
+ *
+ *
+ * Equivalent double code: |
+ * none
+ * |
| Approximate error bound: |
+ * ?
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 19
+ * |
+ */
+ public void fromDHMS() {
+ if (!(this.exponent >= 0 && this.mantissa != 0))
+ return;
+ boolean negative = (this.sign!=0);
+ abs();
+ int Y,M,D,m;
+ long h;
+ h = toLong();
+ frac();
+ tmp1.assign(100);
+ mul(tmp1);
+ m = toInteger();
+ frac();
+ mul(tmp1);
+ // MAGIC ROUNDING: Check if we are 2**-10 second short of 100 seconds
+ // i.e. "seconds" > 99.999
+ { tmp2.mantissa = ONE.mantissa; tmp2.exponent = ONE.exponent; tmp2.sign = ONE.sign; };
+ tmp2.scalbn(-10);
+ add(tmp2);
+ if (this.compare(tmp1) >= 0) {
+ // Yes. So set zero secs instead and carry over to mins and hours
+ { this.mantissa = ZERO.mantissa; this.exponent = ZERO.exponent; this.sign = ZERO.sign; };
+ m++;
+ if (m >= 100) {
+ m -= 100;
+ h++;
+ }
+ // Phew! That was close. From now on it is integer arithmetic...
+ } else {
+ // Nope. So try to undo the damage...
+ sub(tmp2);
+ }
+ D = (int)(h/100);
+ h %= 100;
+ if (D>=10000) {
+ M = D/100;
+ D %= 100;
+ if (D==0) D=1;
+ Y = M/100;
+ M %= 100;
+ if (M==0) M=1;
+ D = gregorian_to_jd(Y,M,D);
+ }
+ add(m*60);
+ div(3600);
+ tmp1.assign(D*24L+h);
+ add(tmp1);
+ if (negative)
+ neg();
+ }
+ /**
+ * Assigns this Real the current time. The time is
+ * encoded into the digits of the number (in decimal form), using the
+ * format "hh.mmss", where "hh" is hours,
+ * "mm" is minutes and "code>ss" is seconds.
+ *
+ *
+ * Equivalent double code: |
+ * none
+ * |
| Error bound: |
+ * ½ ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 8.9
+ * |
+ */
+ public void time() {
+ long now = System.currentTimeMillis();
+ int h,m,s;
+ now /= 1000;
+ s = (int)(now % 60);
+ now /= 60;
+ m = (int)(now % 60);
+ now /= 60;
+ h = (int)(now % 24);
+ assign((h*100+m)*100+s);
+ div(10000);
+ }
+ /**
+ * Assigns this Real the current date. The date is
+ * encoded into the digits of the number (in decimal form), using
+ * the format "YYYYMMDD00", where "YYYY"
+ * is the year, "MM" is the month (1-12) and
+ * "DD" is the day of the month (1-31). The
+ * "00" in this format is a sort of padding to make it
+ * compatible with the format used by {@link #toDHMS()} and {@link
+ * #fromDHMS()}.
+ *
+ *
+ * Equivalent double code: |
+ * none
+ * |
| Error bound: |
+ * 0 ULPs
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 30
+ * |
+ */
+ public void date() {
+ long now = System.currentTimeMillis();
+ now /= 86400000; // days
+ now *= 24; // hours
+ assign(now);
+ add(719528*24); // 1970-01-01 era
+ toDHMS();
+ }
+ //*************************************************************************
+ /**
+ * The seed of the first 64-bit CRC generator of the random
+ * routine. Set this value to control the generated sequence of random
+ * numbers. Should never be set to 0. See {@link #random()}.
+ * Initialized to mantissa of pi.
+ */
+ public static long randSeedA = 0x6487ed5110b4611aL;
+ /**
+ * The seed of the second 64-bit CRC generator of the random
+ * routine. Set this value to control the generated sequence of random
+ * numbers. Should never be set to 0. See {@link #random()}.
+ * Initialized to mantissa of e.
+ */
+ public static long randSeedB = 0x56fc2a2c515da54dL;
+ // 64 Bit CRC Generators
+ //
+ // The generators used here are not cryptographically secure, but
+ // two weak generators are combined into one strong generator by
+ // skipping bits from one generator whenever the other generator
+ // produces a 0-bit.
+ private static void advanceBit() {
+ randSeedA = (randSeedA<<1)^(randSeedA<0?0x1b:0);
+ randSeedB = (randSeedB<<1)^(randSeedB<0?0xb000000000000001L:0);
+ }
+ // Get next bits from the pseudo-random sequence
+ private static long nextBits(int bits) {
+ long answer = 0;
+ while (bits-- > 0) {
+ while (randSeedA >= 0)
+ advanceBit();
+ answer = (answer<<1) + (randSeedB < 0 ? 1 : 0);
+ advanceBit();
+ }
+ return answer;
+ }
+ /**
+ * Accumulate more randomness into the random number generator, to
+ * decrease the predictability of the output from {@link
+ * #random()}. The input should contain data with some form of
+ * inherent randomness e.g. System.currentTimeMillis().
+ *
+ * @param seed some extra randomness for the random number generator.
+ */
+ public static void accumulateRandomness(long seed) {
+ randSeedA ^= seed & 0x5555555555555555L;
+ randSeedB ^= seed & 0xaaaaaaaaaaaaaaaaL;
+ nextBits(63);
+ }
+ /**
+ * Calculates a pseudorandom number in the range [0, 1).
+ * Replaces the contents of this Real with the result.
+ *
+ * The algorithm used is believed to be cryptographically secure,
+ * combining two relatively weak 64-bit CRC generators into a strong
+ * generator by skipping bits from one generator whenever the other
+ * generator produces a 0-bit. The algorithm passes the ent test.
+ *
+ *
+ * Equivalent double code: |
+ * this = Math.{@link Math#random() random}();
+ * |
| Approximate error bound: |
+ * -
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 81
+ * |
+ */
+ public void random() {
+ sign = 0;
+ exponent = 0x3fffffff;
+ while (nextBits(1) == 0)
+ exponent--;
+ mantissa = 0x4000000000000000L+nextBits(62);
+ }
+ //*************************************************************************
+ private int digit(char a, int base, boolean twosComplement) {
+ int digit = -1;
+ if (a>='0' && a<='9')
+ digit = a-'0';
+ else if (a>='A' && a<='F')
+ digit = a-'A'+10;
+ if (digit >= base)
+ return -1;
+ if (twosComplement)
+ digit ^= base-1;
+ return digit;
+ }
+ private void shiftUp(int base) {
+ if (base==2)
+ scalbn(1);
+ else if (base==8)
+ scalbn(3);
+ else if (base==16)
+ scalbn(4);
+ else
+ mul10();
+ }
+ private void atof(String a, int base) {
+ makeZero();
+ int length = a.length();
+ int index = 0;
+ byte tmpSign = 0;
+ boolean compl = false;
+ while (index=0) {
+ shiftUp(base);
+ add(d);
+ index++;
+ }
+ int exp=0;
+ if (index=0) {
+ shiftUp(base);
+ add(d);
+ exp--;
+ index++;
+ }
+ }
+ if (compl)
+ add(ONE);
+ while (index='0' &&
+ a.charAt(index)<='9')
+ {
+ // This takes care of overflows and makes inf or 0
+ if (exp2 < 400000000)
+ exp2 = exp2*10 + a.charAt(index) - '0';
+ index++;
+ }
+ if (expNeg)
+ exp2 = -exp2;
+ exp += exp2;
+ }
+ if (base==2)
+ scalbn(exp);
+ else if (base==8)
+ scalbn(exp*3);
+ else if (base==16)
+ scalbn(exp*4);
+ else {
+ if (exp > 300000000 || exp < -300000000) {
+ // Kludge to be able to enter very large and very small
+ // numbers without causing over/underflows
+ { tmp1.mantissa = TEN.mantissa; tmp1.exponent = TEN.exponent; tmp1.sign = TEN.sign; };
+ if (exp<0) {
+ tmp1.pow(-exp/2);
+ div(tmp1);
+ } else {
+ tmp1.pow(exp/2);
+ mul(tmp1);
+ }
+ exp -= exp/2;
+ }
+ { tmp1.mantissa = TEN.mantissa; tmp1.exponent = TEN.exponent; tmp1.sign = TEN.sign; };
+ if (exp<0) {
+ tmp1.pow(-exp);
+ div(tmp1);
+ } else if (exp>0) {
+ tmp1.pow(exp);
+ mul(tmp1);
+ }
+ }
+ sign = tmpSign;
+ }
+ //*************************************************************************
+ private void normalizeDigits(byte[] digits, int nDigits, int base) {
+ byte carry = 0;
+ boolean isZero = true;
+ for (int i=nDigits-1; i>=0; i--) {
+ if (digits[i] != 0)
+ isZero = false;
+ digits[i] += carry;
+ carry = 0;
+ if (digits[i] >= base) {
+ digits[i] -= base;
+ carry = 1;
+ }
+ }
+ if (isZero) {
+ exponent = 0;
+ return;
+ }
+ if (carry != 0) {
+ if (digits[nDigits-1] >= base/2)
+ digits[nDigits-2] ++; // Rounding, may be inaccurate
+ System.arraycopy(digits, 0, digits, 1, nDigits-1);
+ digits[0] = carry;
+ exponent++;
+ if (digits[nDigits-1] >= base) {
+ // Oh, no, not again!
+ normalizeDigits(digits, nDigits, base);
+ }
+ }
+ while (digits[0] == 0) {
+ System.arraycopy(digits, 1, digits, 0, nDigits-1);
+ digits[nDigits-1] = 0;
+ exponent--;
+ }
+ }
+ private int getDigits(byte[] digits, int base) {
+ if (base == 10)
+ {
+ { tmp1.mantissa = this.mantissa; tmp1.exponent = this.exponent; tmp1.sign = this.sign; };
+ tmp1.abs();
+ { tmp2.mantissa = tmp1.mantissa; tmp2.exponent = tmp1.exponent; tmp2.sign = tmp1.sign; };
+ int exp = exponent = tmp1.lowPow10();
+ exp -= 18;
+ boolean exp_neg = exp <= 0;
+ exp = Math.abs(exp);
+ if (exp > 300000000) {
+ // Kludge to be able to print very large and very small numbers
+ // without causing over/underflows
+ { tmp1.mantissa = TEN.mantissa; tmp1.exponent = TEN.exponent; tmp1.sign = TEN.sign; };
+ tmp1.pow(exp/2); // So, divide twice by not-so-extreme numbers
+ if (exp_neg)
+ tmp2.mul(tmp1);
+ else
+ tmp2.div(tmp1);
+ { tmp1.mantissa = TEN.mantissa; tmp1.exponent = TEN.exponent; tmp1.sign = TEN.sign; };
+ tmp1.pow(exp-(exp/2));
+ } else {
+ { tmp1.mantissa = TEN.mantissa; tmp1.exponent = TEN.exponent; tmp1.sign = TEN.sign; };
+ tmp1.pow(exp);
+ }
+ if (exp_neg)
+ tmp2.mul(tmp1);
+ else
+ tmp2.div(tmp1);
+ long a;
+ if (tmp2.exponent > 0x4000003e) {
+ tmp2.exponent--;
+ tmp2.round();
+ a = tmp2.toLong();
+ if (a >= 5000000000000000000L) { // Rounding up gave 20 digits
+ exponent++;
+ a /= 5;
+ digits[18] = (byte)(a%10);
+ a /= 10;
+ } else {
+ digits[18] = (byte)((a%5)*2);
+ a /= 5;
+ }
+ } else {
+ tmp2.round();
+ a = tmp2.toLong();
+ digits[18] = (byte)(a%10);
+ a /= 10;
+ }
+ for (int i=17; i>=0; i--) {
+ digits[i] = (byte)(a%10);
+ a /= 10;
+ }
+ digits[19] = 0;
+ return 19;
+ }
+ int accurateBits = 64;
+ int bitsPerDigit = base == 2 ? 1 : base == 8 ? 3 : 4;
+ if ((this.exponent == 0 && this.mantissa == 0)) {
+ sign = 0; // Two's complement cannot display -0
+ } else {
+ if ((this.sign!=0)) {
+ mantissa = -mantissa;
+ if (((mantissa >> 62)&3) == 3) {
+ mantissa <<= 1;
+ exponent--;
+ accurateBits--; // ?
+ }
+ }
+ exponent -= 0x40000000-1;
+ int shift = bitsPerDigit-1 -
+ floorMod(exponent, bitsPerDigit);
+ exponent = floorDiv(exponent, bitsPerDigit);
+ if (shift == bitsPerDigit-1) {
+ // More accurate to shift up instead
+ mantissa <<= 1;
+ exponent--;
+ accurateBits--;
+ }
+ else if (shift>0) {
+ mantissa = (mantissa+(1L<<(shift-1)))>>>shift;
+ if ((this.sign!=0)) {
+ // Need to fill in some 1's at the top
+ // (">>", not ">>>")
+ mantissa |= 0x8000000000000000L>>(shift-1);
+ }
+ }
+ }
+ int accurateDigits = (accurateBits+bitsPerDigit-1)/bitsPerDigit;
+ for (int i=0; i>>(64-bitsPerDigit));
+ mantissa <<= bitsPerDigit;
+ }
+ digits[accurateDigits] = 0;
+ return accurateDigits;
+ }
+ private boolean carryWhenRounded(byte[] digits, int nDigits, int base) {
+ if (digits[nDigits] < base/2)
+ return false; // no rounding up, no carry
+ for (int i=nDigits-1; i>=0; i--)
+ if (digits[i] < base-1)
+ return false; // carry would not propagate
+ exponent++;
+ digits[0] = 1;
+ for (int i=1; i= base/2) {
+ digits[nDigits-1]++;
+ normalizeDigits(digits, nDigits, base);
+ }
+ }
+ /**
+ * The number format used to convert Real values to
+ * String using {@link Real#toString(Real.NumberFormat)
+ * Real.toString()}. The default number format uses base-10, maximum
+ * precision, removal of trailing zeros and '.' as radix point.
+ *
+ * Note that the fields of NumberFormat are not
+ * protected in any way, the user is responsible for setting the
+ * correct values to get a correct result.
+ */
+ public static class NumberFormat
+ {
+ /**
+ * The number base of the conversion. The default value is 10,
+ * valid options are 2, 8, 10 and 16. See {@link Real#and(Real)
+ * Real.and()} for an explanation of the interpretation of a
+ * Real in base 2, 8 and 16.
+ *
+ *
Negative numbers output in base-2, base-8 and base-16 are
+ * shown in two's complement form. This form guarantees that a
+ * negative number starts with at least one digit that is the
+ * maximum digit for that base, i.e. '1', '7', and 'F',
+ * respectively. A positive number is guaranteed to start with at
+ * least one '0'. Both positive and negative numbers are extended
+ * to the left using this digit, until {@link #maxwidth} is
+ * reached.
+ */
+ public int base = 10;
+ /**
+ * Maximum width of the converted string. The default value is 30.
+ * If the conversion of a Real with a given {@link
+ * #precision} would produce a string wider than
+ * maxwidth, precision is reduced until
+ * the number fits within the given width. If
+ * maxwidth is too small to hold the number with its
+ * sign, exponent and a precision of 1 digit, the
+ * string may become wider than maxwidth.
+ *
+ *
If align is set to anything but
+ * ALIGN_NONE and the converted string is shorter
+ * than maxwidth, the resulting string is padded with
+ * spaces to the specified width according to the alignment.
+ */
+ public int maxwidth = 30;
+ /**
+ * The precision, or number of digits after the radix point in the
+ * converted string when using the FIX, SCI or
+ * ENG format (see {@link #fse}). The default value is 16,
+ * valid values are 0-16 for base-10 and base-16 conversion, 0-21
+ * for base-8 conversion, and 0-63 for base-2 conversion.
+ *
+ *
The precision may be reduced to make the number
+ * fit within {@link #maxwidth}. The precision is
+ * also reduced if it is set higher than the actual numbers of
+ * significant digits in a Real. When
+ * fse is set to FSE_NONE, i.e. "normal"
+ * output, the precision is always at maximum, but trailing zeros
+ * are removed.
+ */
+ public int precision = 16;
+ /**
+ * The special output formats FIX, SCI or ENG
+ * are enabled with this field. The default value is
+ * FSE_NONE. Valid options are listed below.
+ *
+ *
Numbers are output in one of two main forms, according to
+ * this setting. The normal form has an optional sign, one or more
+ * digits before the radix point, and zero or more digits after the
+ * radix point, for example like this:
+ * 3.14159
+ * The exponent form is like the normal form, followed by an
+ * exponent marker 'e', an optional sign and one or more exponent
+ * digits, for example like this:
+ * -3.4753e-13
+ *
+ *
+ * - {@link #FSE_NONE}
+ *
- Normal output. Numbers are output with maximum precision,
+ * trailing zeros are removed. The format is changed to
+ * exponent form if the number is larger than the number of
+ * significant digits allows, or if the resulting string would
+ * exceed
maxwidth without the exponent form.
+ *
+ * - {@link #FSE_FIX}
+ *
- Like normal output, but the numbers are output with a
+ * fixed number of digits after the radix point, according to
+ * {@link #precision}. Trailing zeros are not removed.
+ *
+ *
- {@link #FSE_SCI}
+ *
- The numbers are always output in the exponent form, with
+ * one digit before the radix point, and a fixed number of
+ * digits after the radix point, according to
+ *
precision. Trailing zeros are not removed.
+ *
+ * - {@link #FSE_ENG}
+ *
- Like the SCI format, but the output shows one to
+ * three digits before the radix point, so that the exponent is
+ * always divisible by 3.
+ *
+ */
+ public int fse = FSE_NONE;
+ /**
+ * The character used as the radix point. The default value is
+ * '.'. Theoretcally any character that does not
+ * otherwise occur in the output can be used, such as
+ * ','.
+ *
+ * Note that setting this to anything but '.' and
+ * ',' is not supported by any conversion method from
+ * String back to Real.
+ */
+ public char point = '.';
+ /**
+ * Set to true to remove the radix point if this is
+ * the last character in the converted string. This is the
+ * default.
+ */
+ public boolean removePoint = true;
+ /**
+ * The character used as the thousands separator. The default
+ * value is the character code 0, which disables
+ * thousands-separation. Theoretcally any character that does not
+ * otherwise occur in the output can be used, such as
+ * ',' or ' '.
+ *
+ *
When thousand!=0, this character is inserted
+ * between every 3rd digit to the left of the radix point in
+ * base-10 conversion. In base-16 conversion, the separator is
+ * inserted between every 4th digit, and in base-2 conversion the
+ * separator is inserted between every 8th digit. In base-8
+ * conversion, no separator is ever inserted.
+ *
+ *
Note that tousands separators are not supported by any
+ * conversion method from String back to
+ * Real, so use of a thousands separator is meant
+ * only for the presentation of numbers.
+ */
+ public char thousand = 0;
+ /**
+ * The alignment of the output string within a field of {@link
+ * #maxwidth} characters. The default value is
+ * ALIGN_NONE. Valid options are defined as follows:
+ *
+ *
+ * - {@link #ALIGN_NONE}
+ *
- The resulting string is not padded with spaces.
+ *
+ *
- {@link #ALIGN_LEFT}
+ *
- The resulting string is padded with spaces on the right side
+ * until a width of
maxwidth is reached, making the
+ * number left-aligned within the field.
+ *
+ * - {@link #ALIGN_RIGHT}
+ *
- The resulting string is padded with spaces on the left side
+ * until a width of
maxwidth is reached, making the
+ * number right-aligned within the field.
+ *
+ * - {@link #ALIGN_CENTER}
+ *
- The resulting string is padded with spaces on both sides
+ * until a width of
maxwidth is reached, making the
+ * number center-aligned within the field.
+ *
+ */
+ public int align = ALIGN_NONE;
+ /** Normal output {@linkplain #fse format} */
+ public static final int FSE_NONE = 0;
+ /** FIX output {@linkplain #fse format} */
+ public static final int FSE_FIX = 1;
+ /** SCI output {@linkplain #fse format} */
+ public static final int FSE_SCI = 2;
+ /** ENG output {@linkplain #fse format} */
+ public static final int FSE_ENG = 3;
+ /** No {@linkplain #align alignment} */
+ public static final int ALIGN_NONE = 0;
+ /** Left {@linkplain #align alignment} */
+ public static final int ALIGN_LEFT = 1;
+ /** Right {@linkplain #align alignment} */
+ public static final int ALIGN_RIGHT = 2;
+ /** Center {@linkplain #align alignment} */
+ public static final int ALIGN_CENTER = 3;
+ }
+ private String align(StringBuffer s, NumberFormat format) {
+ if (format.align == NumberFormat.ALIGN_LEFT) {
+ while (s.length()"0123456789ABCDEF".
+ * See {@link #assign(String,int)}.
+ */
+ public static final String hexChar = "0123456789ABCDEF";
+ private String ftoa(NumberFormat format) {
+ ftoaBuf.setLength(0);
+ if ((this.exponent < 0 && this.mantissa != 0)) {
+ ftoaBuf.append("nan");
+ return align(ftoaBuf,format);
+ }
+ if ((this.exponent < 0 && this.mantissa == 0)) {
+ ftoaBuf.append((this.sign!=0) ? "-inf":"inf");
+ return align(ftoaBuf,format);
+ }
+ int digitsPerThousand;
+ switch (format.base) {
+ case 2:
+ digitsPerThousand = 8;
+ break;
+ case 8:
+ digitsPerThousand = 1000; // Disable thousands separator
+ break;
+ case 16:
+ digitsPerThousand = 4;
+ break;
+ case 10:
+ default:
+ digitsPerThousand = 3;
+ break;
+ }
+ if (format.thousand == 0)
+ digitsPerThousand = 1000; // Disable thousands separator
+ { tmp4.mantissa = this.mantissa; tmp4.exponent = this.exponent; tmp4.sign = this.sign; };
+ int accurateDigits = tmp4.getDigits(ftoaDigits, format.base);
+ if (format.base == 10 && (exponent > 0x4000003e || !isIntegral()))
+ accurateDigits = 16; // Only display 16 digits for non-integers
+ int precision;
+ int pointPos = 0;
+ do
+ {
+ int width = format.maxwidth-1; // subtract 1 for decimal point
+ int prefix = 0;
+ if (format.base != 10)
+ prefix = 1; // want room for at least one "0" or "f/7/1"
+ else if ((tmp4.sign!=0))
+ width--; // subtract 1 for sign
+ boolean useExp = false;
+ switch (format.fse) {
+ case NumberFormat.FSE_SCI:
+ precision = format.precision+1;
+ useExp = true;
+ break;
+ case NumberFormat.FSE_ENG:
+ pointPos = floorMod(tmp4.exponent,3);
+ precision = format.precision+1+pointPos;
+ useExp = true;
+ break;
+ case NumberFormat.FSE_FIX:
+ case NumberFormat.FSE_NONE:
+ default:
+ precision = 1000;
+ if (format.fse == NumberFormat.FSE_FIX)
+ precision = format.precision+1;
+ if (tmp4.exponent+1 >
+ width-(tmp4.exponent+prefix)/digitsPerThousand-prefix+
+ (format.removePoint ? 1:0) ||
+ tmp4.exponent+1 > accurateDigits ||
+ -tmp4.exponent >= width ||
+ -tmp4.exponent >= precision)
+ {
+ useExp = true;
+ } else {
+ pointPos = tmp4.exponent;
+ precision += tmp4.exponent;
+ if (tmp4.exponent > 0)
+ width -= (tmp4.exponent+prefix)/digitsPerThousand;
+ if (format.removePoint && tmp4.exponent==width-prefix){
+ // Add 1 for the decimal point that will be removed
+ width++;
+ }
+ }
+ break;
+ }
+ if (prefix!=0 && pointPos>=0)
+ width -= prefix;
+ ftoaExp.setLength(0);
+ if (useExp) {
+ ftoaExp.append('e');
+ ftoaExp.append(tmp4.exponent-pointPos);
+ width -= ftoaExp.length();
+ }
+ if (precision > accurateDigits)
+ precision = accurateDigits;
+ if (precision > width)
+ precision = width;
+ if (precision > width+pointPos) // In case of negative pointPos
+ precision = width+pointPos;
+ if (precision <= 0)
+ precision = 1;
+ }
+ while (tmp4.carryWhenRounded(ftoaDigits,precision,format.base));
+ tmp4.round(ftoaDigits,precision,format.base);
+ // Start generating the string. First the sign
+ if ((tmp4.sign!=0) && format.base == 10)
+ ftoaBuf.append('-');
+ // Save pointPos for hex/oct/bin prefixing with thousands-sep
+ int pointPos2 = pointPos < 0 ? 0 : pointPos;
+ // Add leading zeros (or f/7/1)
+ char prefixChar = (format.base==10 || (tmp4.sign==0)) ? '0' :
+ hexChar.charAt(format.base-1);
+ if (pointPos < 0) {
+ ftoaBuf.append(prefixChar);
+ ftoaBuf.append(format.point);
+ while (pointPos < -1) {
+ ftoaBuf.append(prefixChar);
+ pointPos++;
+ }
+ }
+ // Add fractional part
+ for (int i=0; i0 && pointPos%digitsPerThousand==0)
+ ftoaBuf.append(format.thousand);
+ if (pointPos == 0)
+ ftoaBuf.append(format.point);
+ pointPos--;
+ }
+ if (format.fse == NumberFormat.FSE_NONE) {
+ // Remove trailing zeros
+ while (ftoaBuf.charAt(ftoaBuf.length()-1) == '0')
+ ftoaBuf.setLength(ftoaBuf.length()-1);
+ }
+ if (format.removePoint) {
+ // Remove trailing point
+ if (ftoaBuf.charAt(ftoaBuf.length()-1) == format.point)
+ ftoaBuf.setLength(ftoaBuf.length()-1);
+ }
+ // Add exponent
+ ftoaBuf.append(ftoaExp.toString());
+ // In case hex/oct/bin number, prefix with 0's or f/7/1's
+ if (format.base!=10) {
+ while (ftoaBuf.length()0 && pointPos2%digitsPerThousand==0)
+ ftoaBuf.insert(0,format.thousand);
+ if (ftoaBuf.length()Real to a String using
+ * the default NumberFormat.
+ *
+ * See {@link Real.NumberFormat NumberFormat} for a description
+ * of the default way that numbers are formatted.
+ *
+ *
+ * Equivalent double code: |
+ * this.toString()
+ * |
|
+ * Execution time relative to add:
+ * |
+ * 130
+ * |
+ *
+ * @return a String representation of this Real.
+ */
+ public String toString() {
+ tmpFormat.base = 10;
+ return ftoa(tmpFormat);
+ }
+ /**
+ * Converts this Real to a String using
+ * the default NumberFormat with base set
+ * according to the argument.
+ *
+ * See {@link Real.NumberFormat NumberFormat} for a description
+ * of the default way that numbers are formatted.
+ *
+ *
+ * Equivalent double code: |
+ *
+ * this.toString() // Works only for base-10
+ * |
|
+ * Execution time relative to add:
+ * | base-2 |
+ * 120
+ * |
| base-8 |
+ * 110
+ * |
| base-10 |
+ * 130
+ * |
| base-16 |
+ * 120
+ * |
+ *
+ * @param base the base for the conversion. Valid base values are
+ * 2, 8, 10 and 16.
+ * @return a String representation of this Real.
+ */
+ public String toString(int base) {
+ tmpFormat.base = base;
+ return ftoa(tmpFormat);
+ }
+ /**
+ * Converts this Real to a String using
+ * the given NumberFormat.
+ *
+ * See {@link Real.NumberFormat NumberFormat} for a description of the
+ * various ways that numbers may be formatted.
+ *
+ *
+ * Equivalent double code: |
+ *
+ * String.format("%...g",this); // Works only for base-10
+ * |
|
+ * Execution time relative to add:
+ * | base-2 |
+ * 120
+ * |
| base-8 |
+ * 110
+ * |
| base-10 |
+ * 130
+ * |
| base-16 |
+ * 120
+ * |
+ *
+ * @param format the number format to use in the conversion.
+ * @return a String representation of this Real.
+ */
+ public String toString(NumberFormat format) {
+ return ftoa(format);
+ }
+}