From cba12afc83b6ba9f80cc85b7e39d07170b5cc00a Mon Sep 17 00:00:00 2001 From: serso Date: Sat, 2 Apr 2016 22:52:54 +0200 Subject: [PATCH] Import Real from midpcalc --- app/src/main/java/midpcalc/Real.java | 6021 ++++++++++++++++++++++++++ 1 file changed, 6021 insertions(+) create mode 100755 app/src/main/java/midpcalc/Real.java 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 + *

+ */ +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: + * + *

+ * + *

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: + *

+ * + *

+ * 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); + } +}