From 6134ff7476933984b1877753527e15e582058f0f Mon Sep 17 00:00:00 2001 From: serso Date: Sat, 2 Apr 2016 22:56:39 +0200 Subject: [PATCH] Automatic code reformat --- app/src/main/java/midpcalc/Real.java | 4657 ++++++++++++++++---------- 1 file changed, 2795 insertions(+), 1862 deletions(-) diff --git a/app/src/main/java/midpcalc/Real.java b/app/src/main/java/midpcalc/Real.java index d6a6e0d0..d8dfe4d2 100755 --- a/app/src/main/java/midpcalc/Real.java +++ b/app/src/main/java/midpcalc/Real.java @@ -1,109 +1,299 @@ - - - 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 -{ +@SuppressWarnings({"unused", "StatementWithEmptyBody"}) +public final class Real { + /** + * 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); + /** + * This string holds the only valid characters to use in hexadecimal + * numbers. Equals "0123456789ABCDEF". + * See {@link #assign(String, int)}. + */ + public static final String hexChar = "0123456789ABCDEF"; + 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}; + /** + * 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; + /** + * 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; + // 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(); + private static byte[] ftoaDigits = new byte[65]; + private static StringBuffer ftoaBuf = new StringBuffer(40); + private static StringBuffer ftoaExp = new StringBuffer(15); + private static NumberFormat tmpFormat = new NumberFormat(); /** * 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 - * + *      -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. @@ -114,11 +304,11 @@ public final class Real * 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 @@ -130,178 +320,23 @@ public final class Real * 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)}. @@ -309,8 +344,13 @@ public final class Real * @param a the Real to assign. */ public Real(Real a) { - { this.mantissa = a.mantissa; this.exponent = a.exponent; this.sign = a.sign; }; + { + 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)}. @@ -320,6 +360,7 @@ public final class Real public Real(int a) { assign(a); } + /** * Creates a new Real, assigning the value of a long * integer. See {@link #assign(long)}. @@ -329,6 +370,7 @@ public final class Real public Real(long a) { assign(a); } + /** * Creates a new Real, assigning the value encoded in a * String using base-10. See {@link #assign(String)}. @@ -336,59 +378,203 @@ public final class Real * @param a the String to assign. */ public Real(String a) { - assign(a,10); + assign(a, 10); } + /** * Creates a new Real, assigning the value encoded in a * String using the specified number base. See {@link - * #assign(String,int)}. + * #assign(String, int)}. * - * @param a the String to assign. + * @param a the String to assign. * @param base the number base of a. Valid base values are 2, - * 8, 10 and 16. + * 8, 10 and 16. */ public Real(String a, int base) { - assign(a,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)}. + * 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; }; + { + 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)}. + * #toBytes(byte[], int) toBytes}. See {@link #assign(byte[], int)}. * - * @param data byte array to decode into this Real. + * @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. + * data[offset]...data[offset+11] will be + * read. */ - public Real(byte [] data, int offset) { - assign(data,offset); + public Real(byte[] data, int offset) { + assign(data, offset); } + + 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; + } + + //************************************************************************* + // 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; + } + + // 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); + } + /** * Assigns this Real the value of another Real. - * + *

*

* Equivalent double code: - * this = a; + * this = a; *
Error bound: - * 0 ULPs + * 0 ULPs *
* Execution time relative to add:   * - * 0.3 + * 0.3 *
* * @param a the Real to assign. @@ -402,61 +588,68 @@ public final class Real 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; + * this = (double)a; *
Error bound: - * 0 ULPs + * 0 ULPs *
* Execution time relative to add:   * - * 0.6 + * 0.6 *
* * @param a the int to assign. */ public void assign(int a) { - if (a==0) { + if (a == 0) { makeZero(); return; } sign = 0; - if (a<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); + 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; + * this = (double)a; *
Error bound: - * 0 ULPs + * 0 ULPs *
* Execution time relative to add:   * - * 1.0 + * 1.0 *
* * @param a the long to assign. */ public void assign(long a) { sign = 0; - if (a<0) { + if (a < 0) { sign = 1; a = -a; // Also works for 0x8000000000000000 } @@ -464,119 +657,122 @@ public final class Real mantissa = a; normalize(); } + /** * Assigns this Real a value encoded in a String - * using base-10, as specified in {@link #assign(String,int)}. - * + * using base-10, as specified in {@link #assign(String, int)}. + *

*

* Equivalent double code: - * this = Double.{@link Double#valueOf(String) valueOf}(a); + * this = Double.{@link Double#valueOf(String) valueOf}(a); *
Approximate error bound: - * ½-1 ULPs + * �-1 ULPs *
* Execution time relative to add:   * - * 80 + * 80 *
* * @param a the String to assign. */ public void assign(String a) { - assign(a,10); + 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 + * this = Double.{@link Double#valueOf(String) valueOf}(a); + * // Works only for base-10 *
* Approximate error bound: * base-10 - * ½-1 ULPs + * �-1 ULPs *
2/8/16 - * ½ ULPs + * � ULPs *
* Execution time relative to add:   * base-2 - * 54 + * 54 *
base-8 - * 60 + * 60 *
base-10 - * 80 + * 80 *
base-16   - * 60 + * 60 *
* - * @param a the String to assign. + * @param a the String to assign. * @param base the number base of a. Valid base values are - * 2, 8, 10 and 16. + * 2, 8, 10 and 16. */ public void assign(String a, int base) { - if (a==null || a.length()==0) { + if (a == null || a.length() == 0) { assign(ZERO); return; } - atof(a,base); + 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); + * this = (1-2*s) * m * + * Math.{@link Math#pow(double, double) + * pow}(2.0,e-0x400000e3); *
Error bound: - * 0 ULPs + * 0 ULPs *
* Execution time relative to add:   * - * 0.3 + * 0.3 *
* * @param s {@link #sign} bit, 0 for positive sign, 1 for negative sign @@ -584,113 +780,116 @@ public final class Real * @param m {@link #mantissa} */ public void assign(int s, int e, long m) { - sign = (byte)s; + 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}. - * + * #toBytes(byte[], int) toBytes}. + *

*

* Error bound: - * 0 ULPs + * 0 ULPs *
* Execution time relative to add:   * - * 1.2 + * 1.2 *
* - * @param data byte array to decode into this Real. + * @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. + * 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))); + 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)}. - * + * #assign(byte[], int)}. + *

*

* Execution time relative to add:   * - * 1.2 + * 1.2 *
* - * @param data byte array to save this Real in. + * @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. + * 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); + 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); + * this = Float.{@link Float#intBitsToFloat(int) + * intBitsToFloat}(bits); *
Error bound: - * 0 ULPs + * 0 ULPs *
* Execution time relative to add:   * - * 0.6 + * 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; + 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) { + if (exponent == 0) { // Degenerate small float - exponent = 0x40000000-126; + exponent = 0x40000000 - 126; normalize(); return; } if (exponent <= 254) { // Normal IEEE 754 float - exponent += 0x40000000-127; - mantissa |= 1L<<62; + exponent += 0x40000000 - 127; + mantissa |= 1L << 62; return; } if (mantissa == 0) @@ -698,43 +897,44 @@ public final class Real 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); + * this = Double.{@link Double#longBitsToDouble(long) + * longBitsToDouble}(bits); *
Error bound: - * 0 ULPs + * 0 ULPs *
* Execution time relative to add:   * - * 0.6 + * 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; + 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) { + if (exponent == 0) { // Degenerate small float - exponent = 0x40000000-1022; + exponent = 0x40000000 - 1022; normalize(); return; } if (exponent <= 2046) { // Normal IEEE 754 float - exponent += 0x40000000-1023; - mantissa |= 1L<<62; + exponent += 0x40000000 - 1023; + mantissa |= 1L << 62; return; } if (mantissa == 0) @@ -742,19 +942,20 @@ public final class Real 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) + * Float.{@link Float#floatToIntBits(float) + * floatToIntBits}(this) *
* Execution time relative to add:   * - * 0.7 + * 0.7 *
* * @return the bits that represent the floating-point number. @@ -762,38 +963,39 @@ public final class Real public int toFloatBits() { if ((this.exponent < 0 && this.mantissa != 0)) return 0x7fffffff; // nan - int e = exponent-0x40000000+127; + int e = exponent - 0x40000000 + 127; long m = mantissa; // Round properly! - m += 1L<<38; - if (m<0) { + m += 1L << 38; + if (m < 0) { m >>>= 1; e++; if (exponent < 0) // Overflow - return (sign<<31)|0x7f800000; // inf + return (sign << 31) | 0x7f800000; // inf } if ((this.exponent < 0 && this.mantissa == 0) || e > 254) - return (sign<<31)|0x7f800000; // inf + return (sign << 31) | 0x7f800000; // inf if ((this.exponent == 0 && this.mantissa == 0) || e < -22) - return (sign<<31); // zero + return (sign << 31); // zero if (e <= 0) // Degenerate small float - return (sign<<31)|((int)(m>>>(40-e))&0x007fffff); + return (sign << 31) | ((int) (m >>> (40 - e)) & 0x007fffff); // Normal IEEE 754 float - return (sign<<31)|(e<<23)|((int)(m>>>39)&0x007fffff); + 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) + * Double.{@link Double#doubleToLongBits(double) + * doubleToLongBits}(this) *
* Execution time relative to add:   * - * 0.7 + * 0.7 *
* * @return the bits that represent the floating-point number. @@ -801,36 +1003,37 @@ public final class Real public long toDoubleBits() { if ((this.exponent < 0 && this.mantissa != 0)) return 0x7fffffffffffffffL; // nan - int e = exponent-0x40000000+1023; + int e = exponent - 0x40000000 + 1023; long m = mantissa; // Round properly! - m += 1L<<9; - if (m<0) { + m += 1L << 9; + if (m < 0) { m >>>= 1; e++; if (exponent < 0) - return ((long)sign<<63)|0x7ff0000000000000L; // inf + return ((long) sign << 63) | 0x7ff0000000000000L; // inf } if ((this.exponent < 0 && this.mantissa == 0) || e > 2046) - return ((long)sign<<63)|0x7ff0000000000000L; // inf + return ((long) sign << 63) | 0x7ff0000000000000L; // inf if ((this.exponent == 0 && this.mantissa == 0) || e < -51) - return ((long)sign<<63); // zero + return ((long) sign << 63); // zero if (e <= 0) // Degenerate small double - return ((long)sign<<63)|((m>>>(11-e))&0x000fffffffffffffL); + return ((long) sign << 63) | ((m >>> (11 - e)) & 0x000fffffffffffffL); // Normal IEEE 754 double - return ((long)sign<<63)|((long)e<<52)|((m>>>10)&0x000fffffffffffffL); + return ((long) sign << 63) | ((long) e << 52) | ((m >>> 10) & 0x000fffffffffffffL); } + /** * Makes this Real the value of positive zero. - * + *

*

* Equivalent double code: - * this = 0; + * this = 0; *
* Execution time relative to add:   * - * 0.2 + * 0.2 *
*/ public void makeZero() { @@ -838,60 +1041,63 @@ public final class Real mantissa = 0; exponent = 0; } + /** * Makes this Real the value of zero with the specified sign. - * + *

*

* Equivalent double code: - * this = 0.0 * (1-2*s); + * this = 0.0 * (1-2*s); *
* Execution time relative to add:   * - * 0.2 + * 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; + 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); + * this = Double.{@link Double#POSITIVE_INFINITY POSITIVE_INFINITY} + * * (1-2*s); *
* Execution time relative to add:   * - * 0.3 + * 0.3 *
* * @param s sign bit, 0 to make positive infinity, 1 to make negative - * infinity + * infinity */ public void makeInfinity(int s) { - sign = (byte)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}; + * this = Double.{@link Double#NaN NaN}; *
* Execution time relative to add:   * - * 0.3 + * 0.3 *
*/ public void makeNan() { @@ -899,187 +1105,196 @@ public final class Real mantissa = 0x4000000000000000L; exponent = 0x80000000; } + /** * Returns true if the value of this Real is * zero, false otherwise. - * + *

*

* Equivalent double code: - * (this == 0) + * (this == 0) *
* Execution time relative to add:   * - * 0.3 + * 0.3 *
* * @return true if the value represented by this object is - * zero, false otherwise. + * 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) + * Double.{@link Double#isInfinite(double) isInfinite}(this) *
* Execution time relative to add:   * - * 0.3 + * 0.3 *
* * @return true if the value represented by this object is - * infinite, false if it is finite or NaN. + * 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) + * Double.{@link Double#isNaN(double) isNaN}(this) *
* Execution time relative to add:   * - * 0.3 + * 0.3 *
* * @return true if the value represented by this object is - * NaN, false otherwise. + * 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)) + * (!Double.{@link Double#isNaN(double) isNaN}(this) && + * !Double.{@link Double#isInfinite(double) + * isInfinite}(this)) *
* Execution time relative to add:   * - * 0.3 + * 0.3 *
* * @return true if the value represented by this object is - * finite, false if it is infinite or NaN. + * 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)) + * (!Double.{@link Double#isNaN(double) isNaN}(this) && + * !Double.{@link Double#isInfinite(double) isInfinite}(this) && + * (this!=0)) *
* Execution time relative to add:   * - * 0.3 + * 0.3 *
* * @return true if the value represented by this object is - * finite and nonzero, false if it is infinite, NaN or - * zero. + * 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) + * (this < 0) *
* Execution time relative to add:   * - * 0.3 + * 0.3 *
* * @return true if the value represented by this object - * is negative, false if it is positive or NaN. + * is negative, false if it is positive or NaN. */ public boolean isNegative() { - return sign!=0; + 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); + * this = Math.{@link Math#abs(double) abs}(this); *
Error bound: - * 0 ULPs + * 0 ULPs *
* Execution time relative to add:   * - * 0.2 + * 0.2 *
*/ public void abs() { sign = 0; } + /** * Negates this Real. * Replaces the contents of this Real with the result. - * + *

*

* Equivalent double code: - * this = -this; + * this = -this; *
Error bound: - * 0 ULPs + * 0 ULPs *
* Execution time relative to add:   * - * 0.2 + * 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); + * this = Math.{@link Math#abs(double) + * abs}(this)*Math.{@link Math#signum(double) signum}(a); *
Error bound: - * 0 ULPs + * 0 ULPs *
* Execution time relative to add:   * - * 0.2 + * 0.2 *
* * @param a the Real to copy the sign from. @@ -1091,6 +1306,7 @@ public final class Real } sign = a.sign; } + /** * Readjusts the mantissa of this Real. The exponent is * adjusted accordingly. This is necessary when the mantissa has been @@ -1099,48 +1315,52 @@ public final class 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 + * � ULPs *
* Execution time relative to add:   * - * 0.7 + * 0.7 *
*/ public void normalize() { if ((this.exponent >= 0)) { - if (mantissa > 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; + 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 ++; + } else if (mantissa < 0) { + mantissa = (mantissa + 1) >>> 1; + exponent++; if (mantissa == 0) { // Ooops, it was 0xffffffffffffffffL mantissa = 0x4000000000000000L; - exponent ++; + exponent++; } if (exponent < 0) // Overflow makeInfinity(sign); - } - else // mantissa == 0 + } else // mantissa == 0 { exponent = 0; } } } + /** * Readjusts the mantissa of a Real with extended * precision. The exponent is adjusted accordingly. This is necessary when @@ -1149,21 +1369,21 @@ public final class Real * 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) + * 2-64 ULPs (i.e. of a normal precision Real) *
* Execution time relative to add:   * - * 0.7 + * 0.7 *
* * @param extra the extra 64 bits of mantissa of this extended precision - * Real. + * Real. * @return the extra 64 bits of mantissa of the resulting extended - * precision Real. + * precision Real. */ public long normalize128(long extra) { if (!(this.exponent >= 0)) @@ -1182,9 +1402,9 @@ public final class Real } } if (mantissa < 0) { - extra = (mantissa<<63)+(extra>>>1); + extra = (mantissa << 63) + (extra >>> 1); mantissa >>>= 1; - exponent ++; + exponent++; if (exponent < 0) { // Overflow makeInfinity(sign); return 0; @@ -1192,13 +1412,20 @@ public final class Real 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; + 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)); + mantissa = (mantissa << clz) + (extra >>> (64 - clz)); extra <<= clz; exponent -= clz; if (exponent < 0) { // Underflow @@ -1207,28 +1434,30 @@ public final class Real } 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 + * � ULPs *
* Execution time relative to add:   * - * 1.0 + * 1.0 *
* * @param extra the extra 64 bits of mantissa of this extended precision - * Real. + * Real. */ public void roundFrom128(long extra) { - mantissa += (extra>>63)&1; + mantissa += (extra >> 63) & 1; normalize(); } + /** * Returns true if this Java object is the same * object as a. Since a Real should be @@ -1240,79 +1469,82 @@ public final class Real * @return true if this object is the same as a. */ public boolean equals(Object a) { - return this==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; + 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) + * (this == a) *
* Execution time relative to add:   * - * 1.0 + * 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. + * 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; + return !invalidCompare(a) && compare(a) == 0; } + /** * Returns true if this Real is equal to * the integer a. - * + *

*

* Equivalent double code: - * (this == a) + * (this == a) *
* Execution time relative to add:   * - * 1.7 + * 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. + * 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. @@ -1321,27 +1553,26 @@ public final class Real * returned. * This distinguishes notEqualTo(a) from the expression * !equalTo(a). - * + *

*

* Equivalent double code: - * (this != a) + * (this != a) *
* Execution time relative to add:   * - * 1.0 + * 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. + * 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; + return !invalidCompare(a) && compare(a) != 0; } + /** * Returns true if this Real is not equal to * the integer a. @@ -1349,257 +1580,258 @@ public final class Real * returned. * This distinguishes notEqualTo(a) from the expression * !equalTo(a). - * + *

*

* Equivalent double code: - * (this != a) + * (this != a) *
* Execution time relative to add:   * - * 1.7 + * 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. + * 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) + * (this < a) *
* Execution time relative to add:   * - * 1.0 + * 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. + * 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; + return !invalidCompare(a) && 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) + * (this < a) *
* Execution time relative to add:   * - * 1.7 + * 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. + * 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) + * (this <= a) *
* Execution time relative to add:   * - * 1.0 + * 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. + * 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; + return !invalidCompare(a) && 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) + * (this <= a) *
* Execution time relative to add:   * - * 1.7 + * 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. + * 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) + * (this > a) *
* Execution time relative to add:   * - * 1.0 + * 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. + * 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; + return !invalidCompare(a) && 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) + * (this > a) *
* Execution time relative to add:   * - * 1.7 + * 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. + * 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) + * (this >= a) *
* Execution time relative to add:   * - * 1.0 + * 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. + * 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; + return !invalidCompare(a) && 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) + * (this >= a) *
* Execution time relative to add:   * - * 1.7 + * 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. + * 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)) + * (Math.{@link Math#abs(double) abs}(this) < + * Math.{@link Math#abs(double) abs}(a)) *
* Execution time relative to add:   * - * 0.5 + * 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. + * 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)) @@ -1607,25 +1839,26 @@ public final class Real 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); + * this *= Math.{@link Math#pow(double, double) pow}(2.0,n); *
Error bound: - * 0 ULPs + * 0 ULPs *
* Execution time relative to add:   * - * 0.3 + * 0.3 *
* * @param n the integer argument. @@ -1635,29 +1868,30 @@ public final class Real return; exponent += n; if (exponent < 0) { - if (n<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); + * this += Math.{@link Math#ulp(double) ulp}(this)*Math.{@link + * Math#signum(double) signum}(a-this); *
Error bound: - * 0 ULPs + * 0 ULPs *
* Execution time relative to add:   * - * 0.8 + * 0.8 *
* * @param a the Real argument. @@ -1673,49 +1907,58 @@ public final class Real 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); + { + 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); + { + 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 ++; + if ((this.sign == 0) ^ dir < 0) { + mantissa++; } else { if (mantissa == 0x4000000000000000L) { mantissa <<= 1; exponent--; } - mantissa --; + 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); + * this = Math.{@link Math#floor(double) floor}(this); *
Approximate error bound: - * 0 ULPs + * 0 ULPs *
* Execution time relative to add:   * - * 0.5 + * 0.5 *
*/ public void floor() { if (!(this.exponent >= 0 && this.mantissa != 0)) return; if (exponent < 0x40000000) { - if ((this.sign==0)) + if ((this.sign == 0)) makeZero(sign); else { exponent = ONE.exponent; @@ -1724,31 +1967,32 @@ public final class Real } return; } - int shift = 0x4000003e-exponent; + 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); + * this = Math.{@link Math#ceil(double) ceil}(this); *
Approximate error bound: - * 0 ULPs + * 0 ULPs *
* Execution time relative to add:   * - * 1.8 + * 1.8 *
*/ public void ceil() { @@ -1756,23 +2000,24 @@ public final class Real 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); + * this = Math.{@link Math#rint(double) rint}(this); *
Approximate error bound: - * 0 ULPs + * 0 ULPs *
* Execution time relative to add:   * - * 1.3 + * 1.3 *
*/ public void round() { @@ -1782,28 +2027,29 @@ public final class Real makeZero(sign); return; } - int shift = 0x4000003e-exponent; + 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); + * this = (double)((long)this); *
Approximate error bound: - * 0 ULPs + * 0 ULPs *
* Execution time relative to add:   * - * 1.2 + * 1.2 *
*/ public void trunc() { @@ -1813,61 +2059,63 @@ public final class Real makeZero(sign); return; } - int shift = 0x4000003e-exponent; + 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); + * this -= (double)((long)this); *
Approximate error bound: - * 0 ULPs + * 0 ULPs *
* Execution time relative to add:   * - * 1.2 + * 1.2 *
*/ public void frac() { if (!(this.exponent >= 0 && this.mantissa != 0) || exponent < 0x40000000) return; - int shift = 0x4000003e-exponent; + 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 + * (int)this *
Approximate error bound: - * 0 ULPs + * 0 ULPs *
* Execution time relative to add:   * - * 0.6 + * 0.6 *
* * @return an int representation of this Real. @@ -1876,40 +2124,41 @@ public final class Real 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; + return ((this.sign == 0)) ? 0x7fffffff : 0x80000001; // 0x80000001, so that you can take -x.toInteger() } if (exponent < 0x40000000) return 0; - int shift = 0x4000003e-exponent; + int shift = 0x4000003e - exponent; if (shift < 32) { - return ((this.sign==0)) ? 0x7fffffff : 0x80000001; + return ((this.sign == 0)) ? 0x7fffffff : 0x80000001; // 0x80000001, so that you can take -x.toInteger() } - return (this.sign==0) ? - (int)(mantissa>>>shift) : -(int)(mantissa>>>shift); + 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 + * (long)this *
Approximate error bound: - * 0 ULPs + * 0 ULPs *
* Execution time relative to add:   * - * 0.5 + * 0.5 *
* * @return a long representation of this Real. @@ -1918,35 +2167,36 @@ public final class Real 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; + return ((this.sign == 0)) ? 0x7fffffffffffffffL : 0x8000000000000001L; // 0x8000000000000001L, so that you can take -x.toLong() } if (exponent < 0x40000000) return 0; - int shift = 0x4000003e-exponent; + int shift = 0x4000003e - exponent; if (shift < 0) { - return ((this.sign==0))? 0x7fffffffffffffffL:0x8000000000000001L; + return ((this.sign == 0)) ? 0x7fffffffffffffffL : 0x8000000000000001L; // 0x8000000000000001L, so that you can take -x.toLong() } - return (this.sign==0) ? (mantissa>>>shift) : -(mantissa>>>shift); + 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) + * (this == (long)this) *
* Execution time relative to add:   * - * 0.6 + * 0.6 *
* * @return true if the value represented by this object - * represents a mathematical integer, false otherwise. + * represents a mathematical integer, false otherwise. */ public boolean isIntegral() { if ((this.exponent < 0 && this.mantissa != 0)) @@ -1955,87 +2205,81 @@ public final class Real 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) + * ((((long)this)&1) == 1) *
* Execution time relative to add:   * - * 0.6 + * 0.6 *
* * @return true if the mathematical integer represented by - * this Real is odd, false otherwise. + * this Real is odd, false otherwise. */ public boolean isOdd() { if (!(this.exponent >= 0 && this.mantissa != 0) || - exponent < 0x40000000 || exponent > 0x4000003e) + exponent < 0x40000000 || exponent > 0x4000003e) return false; - int shift = 0x4000003e-exponent; - return ((mantissa>>>shift)&1) != 0; + 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; + * tmp=this; this=a; a=tmp; *
Error bound: - * 0 ULPs + * 0 ULPs *
* Execution time relative to add:   * - * 0.5 + * 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; + 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; + * this += a; *
Error bound: - * ½ ULPs + * � ULPs *
* Execution time relative to add:   * - * «« 1.0 »» + * �� 1.0 �� *
* * @param a the Real to add to this. @@ -2053,18 +2297,20 @@ public final class Real 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)) - { this.mantissa = a.mantissa; this.exponent = a.exponent; this.sign = a.sign; }; - if ((this.exponent == 0 && this.mantissa == 0)) - sign=0; + sign = 0; return; } byte s; int e; long m; if (exponent > a.exponent || - (exponent == a.exponent && mantissa>=a.mantissa)) - { + (exponent == a.exponent && mantissa >= a.mantissa)) { s = a.sign; e = a.exponent; m = a.mantissa; @@ -2076,38 +2322,38 @@ public final class Real exponent = a.exponent; mantissa = a.mantissa; } - int shift = exponent-e; - if (shift>=64) + 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 + 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 ++; + mantissa = (mantissa + 1) >>> 1; + exponent++; if (exponent < 0) { // Overflow makeInfinity(sign); return; } } } else { - if (shift>0) { + if (shift > 0) { // Shift mantissa up to increase accuracy mantissa <<= 1; - exponent --; - shift --; + 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 + 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) { + 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... @@ -2118,8 +2364,8 @@ public final class Real // 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; + (exponent == 0x4000003e && mantissa + m > 0)) { + long mask = (1 << (0x4000003e - exponent)) - 1; if ((mantissa & mask) != 0 || (m & mask) != 0) mantissa = 0; } else @@ -2129,23 +2375,24 @@ public final class Real } // else... if (shift>=1 && mantissa>=0) it should be a-ok } if ((this.exponent == 0 && this.mantissa == 0)) - sign=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; + * this += a; *
Error bound: - * ½ ULPs + * � ULPs *
* Execution time relative to add:   * - * 1.8 + * 1.8 *
* * @param a the int to add to this. @@ -2154,41 +2401,42 @@ public final class Real 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; + * this += a; *
Approximate error bound: - * 2-62 ULPs (i.e. of a normal precision Real) + * 2-62 ULPs (i.e. of a normal precision Real) *
* Execution time relative to add:   * - * 2.0 + * 2.0 *
* - * @param extra the extra 64 bits of mantissa of this extended precision - * Real. - * @param a the Real to add to this. + * @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. + * value a. * @return the extra 64 bits of mantissa of the resulting extended - * precision Real. + * precision Real. */ public long add128(long extra, Real a, long aExtra) { if ((this.exponent < 0 && this.mantissa != 0) || (a.exponent < 0 && a.mantissa != 0)) { @@ -2204,11 +2452,15 @@ public final class Real } 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; }; + { + this.mantissa = a.mantissa; + this.exponent = a.exponent; + this.sign = a.sign; + } extra = aExtra; } if ((this.exponent == 0 && this.mantissa == 0)) - sign=0; + sign = 0; return extra; } byte s; @@ -2216,10 +2468,9 @@ public final class Real long m; long x; if (exponent > a.exponent || - (exponent == a.exponent && mantissa>a.mantissa) || - (exponent == a.exponent && mantissa==a.mantissa && - (extra>>>1)>=(aExtra>>>1))) - { + (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; @@ -2234,25 +2485,25 @@ public final class Real mantissa = a.mantissa; extra = aExtra; } - int shift = exponent-e; - if (shift>=127) + int shift = exponent - e; + if (shift >= 127) return extra; - if (shift>=64) { - x = m>>>(shift-64); + if (shift >= 64) { + x = m >>> (shift - 64); m = 0; - } else if (shift>0) { - x = (x>>>shift)+(m<<(64-shift)); + } 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 += (extra >> 63) & 1; mantissa += m; } else { extra -= x; - mantissa -= (extra>>63)&1; + 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... @@ -2262,29 +2513,30 @@ public final class Real extra <<= 1; extra = normalize128(extra); if ((this.exponent == 0 && this.mantissa == 0)) - sign=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 + * #add128(long, Real, long) add128}(extra,a,aExtra), since only * the sign bit of a need to be changed.) - * + *

*

* Equivalent double code: - * this -= a; + * this -= a; *
Error bound: - * ½ ULPs + * � ULPs *
* Execution time relative to add:   * - * 2.0 + * 2.0 *
* * @param a the Real to subtract from this. @@ -2292,24 +2544,25 @@ public final class Real public void sub(Real a) { tmp0.mantissa = a.mantissa; tmp0.exponent = a.exponent; - tmp0.sign = (byte)(a.sign^1); + 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; + * this -= a; *
Error bound: - * ½ ULPs + * � ULPs *
* Execution time relative to add:   * - * 2.4 + * 2.4 *
* * @param a the int to subtract from this. @@ -2318,20 +2571,21 @@ public final class Real 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; + * this *= a; *
Error bound: - * ½ ULPs + * � ULPs *
* Execution time relative to add:   * - * 1.3 + * 1.3 *
* * @param a the Real to multiply to this. @@ -2357,12 +2611,12 @@ public final class Real long a1 = mantissa >>> 31; long b0 = a.mantissa & 0x7fffffff; long b1 = a.mantissa >>> 31; - mantissa = a1*b1; + 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 round = (mantissa < 0) ? 0 : 0x40000000; + mantissa += ((a0 * b1 + a1 * b0 + ((a0 * b0) >>> 31) + round) >>> 31); int aExp = a.exponent; - exponent += aExp-0x40000000; + exponent += aExp - 0x40000000; if (exponent < 0) { if (exponent == -1 && aExp < 0x40000000 && mantissa < 0) { // Not underflow after all, it will be corrected in the @@ -2377,27 +2631,28 @@ public final class Real } // Simplified normalize() if (mantissa < 0) { - mantissa = (mantissa+1)>>>1; - exponent ++; + 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; + * this *= a; *
Error bound: - * ½ ULPs + * � ULPs *
* Execution time relative to add:   * - * 1.3 + * 1.3 *
* * @param a the int to multiply to this. @@ -2405,11 +2660,11 @@ public final class Real public void mul(int a) { if ((this.exponent < 0 && this.mantissa != 0)) return; - if (a<0) { + if (a < 0) { sign ^= 1; a = -a; } - if ((this.exponent == 0 && this.mantissa == 0) || a==0) { + if ((this.exponent == 0 && this.mantissa == 0) || a == 0) { if ((this.exponent < 0 && this.mantissa == 0)) makeNan(); else @@ -2419,9 +2674,14 @@ public final class Real 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; + 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 @@ -2430,48 +2690,49 @@ public final class Real long a0 = mantissa & 0x7fffffff; long a1 = mantissa >>> 31; long b0 = a & 0xffffffffL; - mantissa = a1*b0; + 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); + int round = (mantissa < 0) ? 0 : 0x40000000; + mantissa += ((a0 * b0 + round) >>> 31); // Simplified normalize() if (mantissa < 0) { - mantissa = (mantissa+1)>>>1; - exponent ++; + 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)}. - * + * #add128(long, Real, long)}. + *

*

* Equivalent double code: - * this *= a; + * this *= a; *
Approximate error bound: - * 2-60 ULPs + * 2-60 ULPs *
* Execution time relative to add:   * - * 3.1 + * 3.1 *
* - * @param extra the extra 64 bits of mantissa of this extended precision - * Real. - * @param a the Real to multiply to this. + * @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. + * value a. * @return the extra 64 bits of mantissa of the resulting extended - * precision Real. + * precision Real. */ public long mul128(long extra, Real a, long aExtra) { if ((this.exponent < 0 && this.mantissa != 0) || (a.exponent < 0 && a.mantissa != 0)) { @@ -2491,7 +2752,7 @@ public final class Real return 0; } int aExp = a.exponent; - exponent += aExp-0x40000000; + exponent += aExp - 0x40000000; if (exponent < 0) { if (aExp < 0x40000000) makeZero(sign); // Underflow @@ -2508,52 +2769,54 @@ public final class Real 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; + 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; + 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 += ((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); + 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; + mantissa += (mantissa + 2) >>> 2; exponent += 3; if (mantissa < 0) { - mantissa = (mantissa+1)>>>1; + 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; + * this = this*this; *
Error bound: - * ½ ULPs + * � ULPs *
* Execution time relative to add:   * - * 1.1 + * 1.1 *
*/ public void sqr() { @@ -2561,7 +2824,7 @@ public final class Real if (!(this.exponent >= 0 && this.mantissa != 0)) return; int e = exponent; - exponent += exponent-0x40000000; + exponent += exponent - 0x40000000; if (exponent < 0) { if (e < 0x40000000) makeZero(sign); // Underflow @@ -2569,77 +2832,40 @@ public final class Real makeInfinity(sign); // Overflow return; } - long a0 = mantissa&0x7fffffff; - long a1 = mantissa>>>31; - mantissa = a1*a1; + 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); + 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 ++; + 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) + * calling {@link #mul128(long, Real, long) * mul128}(extra,a,aExtra).) - * + *

*

* Equivalent double code: - * this /= a; + * this /= a; *
Error bound: - * ½ ULPs + * � ULPs *
* Execution time relative to add:   * - * 2.6 + * 2.6 *
* * @param a the Real to divide this with. @@ -2668,7 +2894,7 @@ public final class Real makeInfinity(sign); return; } - exponent += 0x40000000-a.exponent; + exponent += 0x40000000 - a.exponent; if (mantissa < a.mantissa) { mantissa <<= 1; exponent--; @@ -2682,23 +2908,24 @@ public final class Real } if (a.mantissa == 0x4000000000000000L) return; - mantissa = ldiv(mantissa,a.mantissa); + 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; + * this /= a; *
Error bound: - * ½ ULPs + * � ULPs *
* Execution time relative to add:   * - * 2.6 + * 2.6 *
* * @param a the int to divide this with. @@ -2706,77 +2933,94 @@ public final class Real public void div(int a) { if ((this.exponent < 0 && this.mantissa != 0)) return; - if (a<0) { + 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) + if (a == 0) makeNan(); return; } - if (a==0) { + if (a == 0) { makeInfinity(sign); return; } long denom = a & 0xffffffffL; - long remainder = mantissa%denom; + 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; + 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; + 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; + * this = a/this; *
Error bound: - * ½ ULPs + * � ULPs *
* Execution time relative to add:   * - * 3.1 + * 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.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; }; + { + 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; + * this = a/this; *
Error bound: - * ½ ULPs + * � ULPs *
* Execution time relative to add:   * - * 3.9 + * 3.9 *
* * @param a the int to be divided by this. @@ -2785,20 +3029,21 @@ public final class Real 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; + * this = 1/this; *
Error bound: - * ½ ULPs + * � ULPs *
* Execution time relative to add:   * - * 2.3 + * 2.3 *
*/ public void recip() { @@ -2812,41 +3057,42 @@ public final class Real makeInfinity(sign); return; } - exponent = 0x80000000-exponent; + exponent = 0x80000000 - exponent; if (mantissa == 0x4000000000000000L) { if (exponent < 0) makeInfinity(sign); // Overflow return; } exponent--; - mantissa = ldiv(0x8000000000000000L,mantissa); + 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)}. - * + * for the extra mantissa. See also {@link #add128(long, Real, long)}. + *

*

* Equivalent double code: - * this = 1/this; + * this = 1/this; *
Approximate error bound: - * 2-60 ULPs + * 2-60 ULPs *
* Execution time relative to add:   * - * 17 + * 17 *
* * @param extra the extra 64 bits of mantissa of this extended precision - * Real. + * Real. * @return the extra 64 bits of mantissa of the resulting extended - * precision Real. + * precision Real. */ public long recip128(long extra) { if ((this.exponent < 0 && this.mantissa != 0)) @@ -2863,27 +3109,35 @@ public final class Real sign = 0; // Special case, simple power of 2 if (mantissa == 0x4000000000000000L && extra == 0) { - exponent = 0x80000000-exponent; - if (exponent<0) // Overflow + exponent = 0x80000000 - exponent; + if (exponent < 0) // Overflow makeInfinity(s); return 0; } // Normalize exponent - int exp = 0x40000000-exponent; + int exp = 0x40000000 - exponent; exponent = 0x40000000; // Save -A - { recipTmp.mantissa = this.mantissa; recipTmp.exponent = this.exponent; recipTmp.sign = this.sign; }; + { + 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); + { + 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 @@ -2891,21 +3145,22 @@ public final class Real 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); + * this = Math.{@link Math#floor(double) floor}(this/a); *
Error bound: - * 0 ULPs + * 0 ULPs *
* Execution time relative to add:   * - * 22 + * 22 *
* * @param a the Real argument. @@ -2933,59 +3188,69 @@ public final class Real makeInfinity(sign); return; } - { tmp0.mantissa = a.mantissa; tmp0.exponent = a.exponent; tmp0.sign = a.sign; }; // tmp0 should be free + { + 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)) - { + 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 + { + 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 + 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)) - { + 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); + 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); + * this = this - + * a*Math.{@link Math#floor(double) floor}(this/a); *
Error bound: - * 0 ULPs + * 0 ULPs *
* Execution time relative to add:   * - * 27 + * 27 *
* * @param a the Real argument. @@ -3015,19 +3280,20 @@ public final class Real makeZero(a.sign); return; } - modInternal(a,0); + 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 @@ -3041,20 +3307,20 @@ public final class Real * 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; + * this &= a; *
Error bound: - * 0 ULPs + * 0 ULPs *
* Execution time relative to add:   * - * 1.5 + * 1.5 *
* * @param a the Real argument @@ -3069,12 +3335,16 @@ public final class Real 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)) + 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)) + (this.sign != 0) && (a.sign != 0)) ; // makeInfinity(1) else makeZero(); @@ -3095,17 +3365,17 @@ public final class Real exponent = a.exponent; mantissa = a.mantissa; } - int shift = exponent-e; - if (shift>=64) { + int shift = exponent - e; + if (shift >= 64) { if (s == 0) makeZero(sign); return; } if (s != 0) m = -m; - if ((this.sign!=0)) + if ((this.sign != 0)) mantissa = -mantissa; - mantissa &= m>>shift; + mantissa &= m >> shift; sign = 0; if (mantissa < 0) { mantissa = -mantissa; @@ -3113,26 +3383,27 @@ public final class Real } 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; + * this |= a; *
Error bound: - * 0 ULPs + * 0 ULPs *
* Execution time relative to add:   * - * 1.6 + * 1.6 *
* * @param a the Real argument @@ -3143,15 +3414,22 @@ public final class Real 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)) { + 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)) + 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 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; @@ -3159,9 +3437,8 @@ public final class Real byte s; int e; long m; - if (((this.sign!=0) && exponent <= a.exponent) || - ((a.sign==0) && exponent >= a.exponent)) - { + if (((this.sign != 0) && exponent <= a.exponent) || + ((a.sign == 0) && exponent >= a.exponent)) { s = a.sign; e = a.exponent; m = a.mantissa; @@ -3173,17 +3450,17 @@ public final class Real exponent = a.exponent; mantissa = a.mantissa; } - int shift = exponent-e; - if (shift>=64 || shift<=-64) + int shift = exponent - e; + if (shift >= 64 || shift <= -64) return; if (s != 0) m = -m; - if ((this.sign!=0)) + if ((this.sign != 0)) mantissa = -mantissa; - if (shift>=0) - mantissa |= m>>shift; + if (shift >= 0) + mantissa |= m >> shift; else - mantissa |= m<<(-shift); + mantissa |= m << (-shift); sign = 0; if (mantissa < 0) { mantissa = -mantissa; @@ -3191,33 +3468,34 @@ public final class Real } 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; + * this ^= a; *
Error bound: - * 0 ULPs + * 0 ULPs *
* Execution time relative to add:   * - * 1.5 + * 1.5 *
* * @param a the Real argument @@ -3228,8 +3506,11 @@ public final class Real 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)) { + 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)) { @@ -3251,14 +3532,14 @@ public final class Real exponent = a.exponent; mantissa = a.mantissa; } - int shift = exponent-e; - if (shift>=64) + int shift = exponent - e; + if (shift >= 64) return; if (s != 0) m = -m; - if ((this.sign!=0)) + if ((this.sign != 0)) mantissa = -mantissa; - mantissa ^= m>>shift; + mantissa ^= m >> shift; sign = 0; if (mantissa < 0) { mantissa = -mantissa; @@ -3266,26 +3547,27 @@ public final class Real } 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; + * this &= ~a; *
Error bound: - * 0 ULPs + * 0 ULPs *
* Execution time relative to add:   * - * 1.5 + * 1.5 *
* * @param a the Real argument @@ -3299,12 +3581,12 @@ public final class Real 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)) + if ((this.sign != 0)) + if ((a.sign != 0)) makeInfinity(0); else makeInfinity(1); - } else if ((a.sign!=0)) { + } else if ((a.sign != 0)) { if ((a.exponent < 0 && a.mantissa == 0)) makeInfinity(0); else @@ -3312,25 +3594,25 @@ public final class Real } return; } - int shift = exponent-a.exponent; - if (shift>=64 || (shift<=-64 && (this.sign==0))) + int shift = exponent - a.exponent; + if (shift >= 64 || (shift <= -64 && (this.sign == 0))) return; long m = a.mantissa; - if ((a.sign!=0)) + if ((a.sign != 0)) m = -m; - if ((this.sign!=0)) + if ((this.sign != 0)) mantissa = -mantissa; - if (shift<0) { - if ((this.sign!=0)) { - if (shift<=-64) + if (shift < 0) { + if ((this.sign != 0)) { + if (shift <= -64) mantissa = ~m; else - mantissa = (mantissa>>(-shift)) & ~m; + mantissa = (mantissa >> (-shift)) & ~m; exponent = a.exponent; } else - mantissa &= ~(m<<(-shift)); + mantissa &= ~(m << (-shift)); } else - mantissa &= ~(m>>shift); + mantissa &= ~(m >> shift); sign = 0; if (mantissa < 0) { mantissa = -mantissa; @@ -3338,24 +3620,26 @@ public final class Real } 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); + * this = Math.{@link Math#sqrt(double) sqrt}(this); *
Approximate error bound: - * 1 ULPs + * 1 ULPs *
* Execution time relative to add:   * - * 19 + * 19 *
*/ public void sqrt() { @@ -3371,79 +3655,108 @@ public final class Real if ((this.exponent < 0 && this.mantissa != 0)) return; if ((this.exponent == 0 && this.mantissa == 0)) { - sign=0; + sign = 0; return; } - if ((this.sign!=0)) { + 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; }; + { + recipTmp.mantissa = this.mantissa; + recipTmp.exponent = this.exponent; + recipTmp.sign = this.sign; + } // normalize to range [0.5, 1) - int e = exponent-0x3fffffff; + 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 + { + 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 + { + 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 + { + sqrtTmp.sign = (byte) 0; + sqrtTmp.exponent = 0x3ffffffe; + sqrtTmp.mantissa = 0x5045ee6baf28677aL; + } + //0.31356706742295303132394 add(sqrtTmp); // adjust for odd powers of 2 - if ((e&1) != 0) + if ((e & 1) != 0) mul(SQRT2); // calculate exponent - exponent += e>>1; + 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; }; + 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); + * this = 1/Math.{@link Math#sqrt(double) sqrt}(this); *
Approximate error bound: - * 1 ULPs + * 1 ULPs *
* Execution time relative to add:   * - * 21 + * 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); + * this = Math.{@link Math#cbrt(double) cbrt}(this); *
Approximate error bound: - * 2 ULPs + * 2 ULPs *
* Execution time relative to add:   * - * 32 + * 32 *
*/ public void cbrt() { @@ -3455,23 +3768,36 @@ public final class 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.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; + 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) + if (expRmd > 0) { + { + recipTmp2.sign = (byte) 0; + recipTmp2.exponent = 0x3fffffff; + recipTmp2.mantissa = 0x6597fa94f5b8f20bL; + } + // cbrt(1/2) mul(recipTmp2); - if (expRmd>1) + 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; }; + for (int i = 0; i < 4; i++) { + { + recipTmp2.mantissa = this.mantissa; + recipTmp2.exponent = this.exponent; + recipTmp2.sign = this.sign; + } sqr(); sqr(); mul(recipTmp); @@ -3483,23 +3809,24 @@ public final class Real 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); + * this = Math.{@link Math#pow(double, double) + * pow}(this,1/a); *
Approximate error bound: - * 2 ULPs + * 2 ULPs *
* Execution time relative to add:   * - * 110 + * 110 *
* * @param n the Real argument. @@ -3509,64 +3836,84 @@ public final class Real makeNan(); return; } - if (n.compare(THREE)==0) { + if (n.compare(THREE) == 0) { cbrt(); // Most probable application of nroot... return; - } else if (n.compare(TWO)==0) { + } 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()) { + 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.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); + * this = Math.{@link Math#hypot(double, double) + * hypot}(this,a); *
Approximate error bound: - * 1 ULPs + * 1 ULPs *
* Execution time relative to add:   * - * 24 + * 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.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)) + 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; }; + { + 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.mantissa = this.mantissa; + expTmp.exponent = this.exponent; + expTmp.sign = this.sign; + } expTmp.add(HALF); expTmp.floor(); int exp = expTmp.toInteger(); @@ -3580,7 +3927,7 @@ public final class Real } // Subtract integer part (this is where we need the extra accuracy) expTmp.neg(); - add128(extra,expTmp,0); + add128(extra, expTmp, 0); /* * Adapted from: * Cephes Math Library Release 2.7: May, 1998 @@ -3592,27 +3939,61 @@ public final class Real */ // 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); + * this = Math.{@link Math#exp(double) exp}(this); *
Approximate error bound: - * 1 ULPs + * 1 ULPs *
* Execution time relative to add:   * - * 31 + * 31 *
*/ public void exp() { - { expTmp.sign=(byte)0; expTmp.exponent=0x40000000; expTmp.mantissa=0x5c551d94ae0bf85dL; }; // log2(e) - long extra = mul128(0,expTmp,0xdf43ff68348e9f44L); + { + 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)); + * this = Math.{@link Math#exp(double) exp}(this * + * Math.{@link Math#log(double) log}(2)); *
Approximate error bound: - * 1 ULPs + * 1 ULPs *
* Execution time relative to add:   * - * 27 + * 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)); + * this = Math.{@link Math#exp(double) exp}(this * + * Math.{@link Math#log(double) log}(10)); *
Approximate error bound: - * 1 ULPs + * 1 ULPs *
* Execution time relative to add:   * - * 31 + * 31 *
*/ public void exp10() { - { expTmp.sign=(byte)0; expTmp.exponent=0x40000001; expTmp.mantissa=0x6a4d3c25e68dc57fL; }; // log2(10) - long extra = mul128(0,expTmp,0x2495fb7fa6d7eda6L); + { + expTmp.sign = (byte) 0; + expTmp.exponent = 0x40000001; + expTmp.mantissa = 0x6a4d3c25e68dc57fL; + } + // log2(10) + long extra = mul128(0, expTmp, 0x2495fb7fa6d7eda6L); exp2Internal(extra); } - private int lnInternal() - { + + private int lnInternal() { if ((this.exponent < 0 && this.mantissa != 0)) return 0; - if ((this.sign!=0)) { + if ((this.sign != 0)) { makeNan(); return 0; } @@ -3708,57 +4102,134 @@ public final class Real * long double logl(long double x); */ // normalize to range [0.5, 1) - int e = exponent-0x3fffffff; + int e = exponent - 0x3fffffff; exponent = 0x3fffffff; // rational appriximation - // log(1+x) = x - x²/2 + x³ P(x)/Q(x) + // 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; }; + { + 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 + { + 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 + { + 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 + { + 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 + { + 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 + { + 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 + { + 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 + { + 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.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 + { + 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 + { + 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 + { + 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 + { + 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 + { + 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.mantissa = expTmp2.mantissa; + expTmp3.exponent = expTmp2.exponent; + expTmp3.sign = expTmp2.sign; + } expTmp3.sqr(); mul(expTmp3); mul(expTmp2); @@ -3767,21 +4238,22 @@ public final class Real 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); + * this = Math.{@link Math#log(double) log}(this); *
Approximate error bound: - * 2 ULPs + * 2 ULPs *
* Execution time relative to add:   * - * 51 + * 51 *
*/ public void ln() { @@ -3790,21 +4262,22 @@ public final class Real 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); + * this = Math.{@link Math#log(double) log}(this)/Math.{@link + * Math#log(double) log}(2); *
Approximate error bound: - * 1 ULPs + * 1 ULPs *
* Execution time relative to add:   * - * 51 + * 51 *
*/ public void log2() { @@ -3812,20 +4285,21 @@ public final class Real 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); + * this = Math.{@link Math#log10(double) log10}(this); *
Approximate error bound: - * 2 ULPs + * 2 ULPs *
* Execution time relative to add:   * - * 53 + * 53 *
*/ public void log10() { @@ -3835,25 +4309,26 @@ public final class Real 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; + * 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 + * � ULPs *
* Execution time relative to add:   * - * 3.6 + * 3.6 *
* * @return the base-10 exponent @@ -3861,35 +4336,56 @@ public final class Real public int lowPow10() { if (!(this.exponent >= 0 && this.mantissa != 0)) return 0; - { tmp2.mantissa = this.mantissa; tmp2.exponent = this.exponent; tmp2.sign = this.sign; }; + { + 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); + 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); + e = (int) (e * 0x4d104d43L >> 32); // Now, e < log10(this) < e+1 - { this.mantissa = TEN.mantissa; this.exponent = TEN.exponent; this.sign = TEN.sign; }; + { + 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); + { + 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.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; }; + { + 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); + * this = Math.{@link Math#pow(double, double) pow}(this, a); *
Approximate error bound: - * 2 ULPs + * 2 ULPs *
* Execution time relative to add:   * - * 110 + * 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; }; + { + 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) + 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.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)) + if (test > 0) { + if ((a.sign == 0)) makeInfinity(0); else makeZero(); - } else if (test<0) { - if ((a.sign!=0)) + } else if (test < 0) { + if ((a.sign != 0)) makeInfinity(0); else makeZero(); @@ -3965,19 +4469,19 @@ public final class Real return; } if ((this.exponent == 0 && this.mantissa == 0)) { - if ((this.sign==0)) { - if ((a.sign==0)) + if ((this.sign == 0)) { + if ((a.sign == 0)) makeZero(); else makeInfinity(0); } else { if (a.isIntegral() && a.isOdd()) { - if ((a.sign==0)) + if ((a.sign == 0)) makeZero(1); else makeInfinity(1); } else { - if ((a.sign==0)) + if ((a.sign == 0)) makeZero(); else makeInfinity(0); @@ -3986,20 +4490,20 @@ public final class Real return; } if ((this.exponent < 0 && this.mantissa == 0)) { - if ((this.sign==0)) { - if ((a.sign==0)) + if ((this.sign == 0)) { + if ((a.sign == 0)) makeInfinity(0); else makeZero(); } else { if (a.isIntegral()) { if (a.isOdd()) { - if ((a.sign==0)) + if ((a.sign == 0)) makeInfinity(1); else makeZero(1); } else { - if ((a.sign==0)) + if ((a.sign == 0)) makeInfinity(0); else makeZero(); @@ -4014,8 +4518,8 @@ public final class Real pow(a.toInteger()); return; } - byte s=0; - if ((this.sign!=0)) { + byte s = 0; + if ((this.sign != 0)) { if (a.isIntegral()) { if (a.isOdd()) s = 1; @@ -4025,71 +4529,106 @@ public final class Real } sign = 0; } - { tmp1.mantissa = a.mantissa; tmp1.exponent = a.exponent; tmp1.sign = a.sign; }; + { + 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.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.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; }; + { + 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.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); + extra = add128(extra, tmp2, 0); // Do exp2 of this multiplied by (fractional part of) exponent - extra = tmp1.mul128(0,this,extra); + extra = tmp1.mul128(0, this, extra); tmp1.exp2Internal(extra); - { this.mantissa = tmp1.mantissa; this.exponent = tmp1.exponent; this.sign = tmp1.sign; }; + { + 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); + * this = Math.{@link Math#pow(double, double) pow}(this, a); *
Error bound: - * ½ ULPs + * � ULPs *
* Execution time relative to add:   * - * 84 + * 84 *
* * @param a the integer argument. */ public void pow(int a) { // Calculate power of integer by successive squaring - boolean recp=false; + 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) { + { + 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); + extra = mul128(extra, expTmp, expTmpExtra); + expTmpExtra = expTmp.mul128(expTmpExtra, expTmp, expTmpExtra); } if (recp) extra = recip128(extra); roundFrom128(extra); } + private void sinInternal() { /* * Adapted from: @@ -4102,33 +4641,77 @@ public final class Real */ // 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); + * this = Math.{@link Math#sin(double) sin}(this); *
Approximate error bound: - * 1 ULPs + * 1 ULPs *
* Execution time relative to add:   * - * 28 + * 28 *
*/ public void sin() { if (!(this.exponent >= 0 && this.mantissa != 0)) { - if (!(this.exponent == 0 && this.mantissa == 0)) + if (!(this.exponent == 0)) makeNan(); return; } // Since sin(-x) = -sin(x) we can make sure that x > 0 boolean negative = false; - if ((this.sign!=0)) { + 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); + 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); @@ -4225,29 +4852,34 @@ public final class Real 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); + * this = Math.{@link Math#cos(double) cos}(this); *
Approximate error bound: - * 1 ULPs + * 1 ULPs *
* Execution time relative to add:   * - * 37 + * 37 *
*/ public void cos() { if ((this.exponent == 0 && this.mantissa == 0)) { - { this.mantissa = ONE.mantissa; this.exponent = ONE.exponent; this.sign = ONE.sign; }; + { + this.mantissa = ONE.mantissa; + this.exponent = ONE.exponent; + this.sign = ONE.sign; + } return; } - if ((this.sign!=0)) + if ((this.sign != 0)) abs(); if (this.compare(PI_4) < 0) { cosInternal(); @@ -4256,48 +4888,58 @@ public final class Real 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); + * this = Math.{@link Math#tan(double) tan}(this); *
Approximate error bound: - * 2 ULPs + * 2 ULPs *
* Execution time relative to add:   * - * 70 + * 70 *
*/ public void tan() { - { tmp4.mantissa = this.mantissa; tmp4.exponent = this.exponent; tmp4.sign = this.sign; }; + { + 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); + * this = Math.{@link Math#asin(double) asin}(this); *
Approximate error bound: - * 3 ULPs + * 3 ULPs *
* Execution time relative to add:   * - * 68 + * 68 *
*/ public void asin() { - { tmp1.mantissa = this.mantissa; tmp1.exponent = this.exponent; tmp1.sign = this.sign; }; + { + tmp1.mantissa = this.mantissa; + tmp1.exponent = this.exponent; + tmp1.sign = this.sign; + } sqr(); neg(); add(ONE); @@ -4305,27 +4947,32 @@ public final class Real 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); + * this = Math.{@link Math#acos(double) acos}(this); *
Approximate error bound: - * 2 ULPs + * 2 ULPs *
* Execution time relative to add:   * - * 67 + * 67 *
*/ public void acos() { - boolean negative = (this.sign!=0); + boolean negative = (this.sign != 0); abs(); - { tmp1.mantissa = this.mantissa; tmp1.exponent = this.exponent; tmp1.sign = this.sign; }; + { + tmp1.mantissa = this.mantissa; + tmp1.exponent = this.exponent; + tmp1.sign = this.sign; + } sqr(); neg(); add(ONE); @@ -4337,21 +4984,22 @@ public final class Real 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); + * this = Math.{@link Math#atan(double) atan}(this); *
Approximate error bound: - * 2 ULPs + * 2 ULPs *
* Execution time relative to add:   * - * 37 + * 37 *
*/ public void atan() { @@ -4368,7 +5016,11 @@ public final class Real 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; }; + { + this.mantissa = PI_2.mantissa; + this.exponent = PI_2.exponent; + this.sign = PI_2.sign; + } sign = s; return; } @@ -4377,7 +5029,11 @@ public final class Real // range reduction boolean addPI_2 = false; boolean addPI_4 = false; - { tmp1.mantissa = SQRT2.mantissa; tmp1.exponent = SQRT2.exponent; tmp1.sign = SQRT2.sign; }; + { + tmp1.mantissa = SQRT2.mantissa; + tmp1.exponent = SQRT2.exponent; + tmp1.sign = SQRT2.sign; + } tmp1.add(ONE); if (this.compare(tmp1) > 0) { addPI_2 = true; @@ -4387,7 +5043,11 @@ public final class Real tmp1.sub(TWO); if (this.compare(tmp1) > 0) { addPI_4 = true; - { tmp1.mantissa = this.mantissa; tmp1.exponent = this.exponent; tmp1.sign = this.sign; }; + { + tmp1.mantissa = this.mantissa; + tmp1.exponent = this.exponent; + tmp1.sign = this.sign; + } tmp1.add(ONE); sub(ONE); div(tmp1); @@ -4395,39 +5055,101 @@ public final class Real } // 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); + * this = Math.{@link Math#atan2(double, double) + * atan2}(this,x); *
Approximate error bound: - * 2 ULPs + * 2 ULPs *
* Execution time relative to add:   * - * 48 + * 48 *
* * @param x the Real argument. @@ -4479,94 +5202,114 @@ public final class Real } 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); + * this = Math.{@link Math#sinh(double) sinh}(this); *
Approximate error bound: - * 2 ULPs + * 2 ULPs *
* Execution time relative to add:   * - * 67 + * 67 *
*/ public void sinh() { - { tmp1.mantissa = this.mantissa; tmp1.exponent = this.exponent; tmp1.sign = this.sign; }; + { + 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); + * this = Math.{@link Math#cosh(double) cosh}(this); *
Approximate error bound: - * 2 ULPs + * 2 ULPs *
* Execution time relative to add:   * - * 66 + * 66 *
*/ public void cosh() { - { tmp1.mantissa = this.mantissa; tmp1.exponent = this.exponent; tmp1.sign = this.sign; }; + { + 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); + * this = Math.{@link Math#tanh(double) tanh}(this); *
Approximate error bound: - * 2 ULPs + * 2 ULPs *
* Execution time relative to add:   * - * 70 + * 70 *
*/ public void tanh() { - { tmp1.mantissa = this.mantissa; tmp1.exponent = this.exponent; tmp1.sign = this.sign; }; + { + 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.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 + * none *
Approximate error bound: - * 2 ULPs + * 2 ULPs *
* Execution time relative to add:   * - * 77 + * 77 *
*/ public void asinh() { @@ -4576,7 +5319,11 @@ public final class Real // values byte s = sign; sign = 0; - { tmp1.mantissa = this.mantissa; tmp1.exponent = this.exponent; tmp1.sign = this.sign; }; + { + tmp1.mantissa = this.mantissa; + tmp1.exponent = this.exponent; + tmp1.sign = this.sign; + } tmp1.sqr(); tmp1.add(ONE); tmp1.sqrt(); @@ -4585,48 +5332,58 @@ public final class Real 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 + * none *
Approximate error bound: - * 2 ULPs + * 2 ULPs *
* Execution time relative to add:   * - * 75 + * 75 *
*/ public void acosh() { - { tmp1.mantissa = this.mantissa; tmp1.exponent = this.exponent; tmp1.sign = this.sign; }; + { + 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 + * none *
Approximate error bound: - * 2 ULPs + * 2 ULPs *
* Execution time relative to add:   * - * 57 + * 57 *
*/ public void atanh() { - { tmp1.mantissa = this.mantissa; tmp1.exponent = this.exponent; tmp1.sign = this.sign; }; + { + tmp1.mantissa = this.mantissa; + tmp1.exponent = this.exponent; + tmp1.sign = this.sign; + } tmp1.neg(); tmp1.add(ONE); add(ONE); @@ -4634,210 +5391,349 @@ public final class Real 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 + * none *
Approximate error bound: - * 15 ULPs + * 15 ULPs *
* Execution time relative to add:   * - * 8-190 + * 8-190 *
*/ public void fact() { if (!(this.exponent >= 0)) return; - if (!this.isIntegral() || this.compare(ZERO)<0 || this.compare(200)>0) - { + 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; }; + { + 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 + * none *
Approximate error bound: - * 100+ ULPs + * 100+ ULPs *
* Execution time relative to add:   * - * 190 + * 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); + boolean negative = (this.sign != 0); abs(); - { tmp1.mantissa = this.mantissa; tmp1.exponent = this.exponent; tmp1.sign = this.sign; }; + { + 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³ + // x>n: 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² + { + 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(); + { + 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); + { + 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); + 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); + 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); + 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); + 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 + { + 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(); + 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); + 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; + { + tmp3.mantissa = ONE.mantissa; + tmp3.exponent = ONE.exponent; + tmp3.sign = ONE.sign; + } + tmp3Extra = 0; + int i = 1; do { - tmp1Extra = tmp1.mul128(tmp1Extra,tmp2,tmp2Extra); + 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); + 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); + 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); + { + 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 + // -x� -1 // e x / 1 3 3*5 3*5*7 // erfc(x) = -------- | 1 - --- + ------ - ------ + ------ - ... | - // sqrt(pi) \ 2x² 2 3 4 / - // (2x²) (2x²) (2x²) + // 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.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.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; }; + 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; }; + { + 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; + { + 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(2 * i - 1); tmp4.mul(tmp3); add(tmp4); i++; - } while (tmp4.exponent-0x40000000>-(digits+2) && 2*i-1 -(digits + 2) && 2 * i - 1 < dxq); mul(tmp2); tmp1.sqr(); tmp1.neg(); tmp1.exp(); mul(tmp1); - { tmp1.sign=(byte)0; tmp1.exponent=0x3fffffff; tmp1.mantissa=0x48375d410a6db447L; }; // 1/sqrt(pi) + { + tmp1.sign = (byte) 0; + tmp1.exponent = 0x3fffffff; + tmp1.mantissa = 0x48375d410a6db447L; + } + // 1/sqrt(pi) mul(tmp1); } + /** * Calculates the complementary error function for this Real. * 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 + * overline;">π ï¿½e-t� dt. It is * related to the error function, erf, by the formula * erfc(x)=1-erf(x). - * + *

*

* Equivalent double code: - * none + * none *
Approximate error bound: - * 219 ULPs + * 219 ULPs *
* Execution time relative to add:   * - * 80-4900 + * 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; }; + { + 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; }; + 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 + { + tmp1.sign = (byte) 0; + tmp1.exponent = 0x40000002; + tmp1.mantissa = 0x570a3d70a3d70a3dL; + } + // 5.44 if (this.lessThan(tmp1)) erfc1Internal(); else @@ -4847,25 +5743,26 @@ public final class Real add(TWO); } } + /** * Calculates the inverse complementary error function for this * Real. * Replaces the contents of this Real with the result. - * + *

*

* Equivalent double code: - * none + * none *
Approximate error bound: - * 219 ULPs + * 219 ULPs *
* Execution time relative to add:   * - * 240-5100 + * 240-5100 *
*/ public void inverfc() { - if ((this.exponent < 0 && this.mantissa != 0) || (this.sign!=0) || this.greaterThan(TWO)) { + if ((this.exponent < 0 && this.mantissa != 0) || (this.sign != 0) || this.greaterThan(TWO)) { makeNan(); return; } @@ -4878,11 +5775,11 @@ public final class Real return; } int sign = ONE.compare(this); - if (sign==0) { + if (sign == 0) { makeZero(); return; } - if (sign<0) { + if (sign < 0) { neg(); add(TWO); } @@ -4894,61 +5791,132 @@ public final class Real // 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; }; + // 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.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 + { + 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 + { + 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 + { + 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.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 + { + 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 + { + 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 + { + 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 + { + 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() + { + 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.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; }; + // 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) + { + 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; }; + { + this.mantissa = sqrtTmp.mantissa; + this.exponent = sqrtTmp.exponent; + this.sign = sqrtTmp.sign; + } mul(tmp5); scalbn(-1); add(ONE); @@ -4957,73 +5925,22 @@ public final class Real add(sqrtTmp); // calculate inverfc(x) = -invphi(x/2)/(sqrt(2)) mul(SQRT1_2); - if (sign>0) + 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 @@ -5033,25 +5950,25 @@ public final class Real * "DD" is the day of the month (1-31). See a thorough * discussion of date calculations here. - * + *

*

* Equivalent double code: - * none + * none *
Approximate error bound: - * ? + * ? *
* Execution time relative to add:   * - * 19 + * 19 *
*/ public void toDHMS() { if (!(this.exponent >= 0 && this.mantissa != 0)) return; - boolean negative = (this.sign!=0); + boolean negative = (this.sign != 0); abs(); - int D,m; + int D, m; long h; h = toLong(); frac(); @@ -5062,12 +5979,20 @@ public final class Real 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.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; }; + { + this.mantissa = ZERO.mantissa; + this.exponent = ZERO.exponent; + this.sign = ZERO.sign; + } m++; if (m >= 60) { m -= 60; @@ -5078,29 +6003,30 @@ public final class Real // Nope. So try to undo the damage... sub(tmp2); } - D = (int)(h/24); + D = (int) (h / 24); h %= 24; if (D >= 366) D = jd_to_gregorian(D); - add(m*100); + add(m * 100); div(10000); - tmp1.assign(D*100L+h); + 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 @@ -5111,25 +6037,25 @@ public final class Real * is 0 it is treated as 1. See a thorough discussion of date * calculations here. - * + *

*

* Equivalent double code: - * none + * none *
Approximate error bound: - * ? + * ? *
* Execution time relative to add:   * - * 19 + * 19 *
*/ public void fromDHMS() { if (!(this.exponent >= 0 && this.mantissa != 0)) return; - boolean negative = (this.sign!=0); + boolean negative = (this.sign != 0); abs(); - int Y,M,D,m; + int Y, M, D, m; long h; h = toLong(); frac(); @@ -5140,12 +6066,20 @@ public final class Real 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.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; }; + { + this.mantissa = ZERO.mantissa; + this.exponent = ZERO.exponent; + this.sign = ZERO.sign; + } m++; if (m >= 100) { m -= 100; @@ -5156,54 +6090,56 @@ public final class Real // Nope. So try to undo the damage... sub(tmp2); } - D = (int)(h/100); + D = (int) (h / 100); h %= 100; - if (D>=10000) { - M = D/100; + if (D >= 10000) { + M = D / 100; D %= 100; - if (D==0) D=1; - Y = M/100; + if (D == 0) D = 1; + Y = M / 100; M %= 100; - if (M==0) M=1; - D = gregorian_to_jd(Y,M,D); + if (M == 0) M = 1; + D = gregorian_to_jd(Y, M, D); } - add(m*60); + add(m * 60); div(3600); - tmp1.assign(D*24L+h); + 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 + * none *
Error bound: - * ½ ULPs + * � ULPs *
* Execution time relative to add:   * - * 8.9 + * 8.9 *
*/ public void time() { long now = System.currentTimeMillis(); - int h,m,s; + int h, m, s; now /= 1000; - s = (int)(now % 60); + s = (int) (now % 60); now /= 60; - m = (int)(now % 60); + m = (int) (now % 60); now /= 60; - h = (int)(now % 24); - assign((h*100+m)*100+s); + 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 @@ -5213,17 +6149,17 @@ public final class Real * "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 + * none *
Error bound: - * 0 ULPs + * 0 ULPs *
* Execution time relative to add:   * - * 30 + * 30 *
*/ public void date() { @@ -5231,78 +6167,30 @@ public final class Real now /= 86400000; // days now *= 24; // hours assign(now); - add(719528*24); // 1970-01-01 era + 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}(); + * this = Math.{@link Math#random() random}(); *
Approximate error bound: - * - + * - *
* Execution time relative to add:   * - * 81 + * 81 *
*/ public void random() { @@ -5310,60 +6198,63 @@ public final class Real exponent = 0x3fffffff; while (nextBits(1) == 0) exponent--; - mantissa = 0x4000000000000000L+nextBits(62); + 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 (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; + digit ^= base - 1; return digit; } + private void shiftUp(int base) { - if (base==2) + if (base == 2) scalbn(1); - else if (base==8) + else if (base == 8) scalbn(3); - else if (base==16) + 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) { + while (index < length && (d = digit(a.charAt(index), base, compl)) >= 0) { shiftUp(base); add(d); index++; } - int exp=0; - if (index=0) { + while (index < length && (d = digit(a.charAt(index), base, compl)) >= 0) { shiftUp(base); add(d); exp--; @@ -5372,66 +6263,74 @@ public final class Real } if (compl) add(ONE); - while (index='0' && - a.charAt(index)<='9') - { + while (index < length && a.charAt(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'; + exp2 = exp2 * 10 + a.charAt(index) - '0'; index++; } if (expNeg) exp2 = -exp2; exp += exp2; } - if (base==2) + if (base == 2) scalbn(exp); - else if (base==8) - scalbn(exp*3); - else if (base==16) - scalbn(exp*4); + 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); + { + 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); + tmp1.pow(exp / 2); mul(tmp1); } - exp -= exp/2; + exp -= exp / 2; } - { tmp1.mantissa = TEN.mantissa; tmp1.exponent = TEN.exponent; tmp1.sign = TEN.sign; }; - if (exp<0) { + { + tmp1.mantissa = TEN.mantissa; + tmp1.exponent = TEN.exponent; + tmp1.sign = TEN.sign; + } + if (exp < 0) { tmp1.pow(-exp); div(tmp1); - } else if (exp>0) { + } 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--) { + for (int i = nDigits - 1; i >= 0; i--) { if (digits[i] != 0) isZero = false; digits[i] += carry; @@ -5446,28 +6345,36 @@ public final class Real 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); + 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) { + 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; + 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; }; + 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; }; + { + tmp2.mantissa = tmp1.mantissa; + tmp2.exponent = tmp1.exponent; + tmp2.sign = tmp1.sign; + } int exp = exponent = tmp1.lowPow10(); exp -= 18; boolean exp_neg = exp <= 0; @@ -5475,16 +6382,28 @@ public final class Real 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 + { + 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)); + { + 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.mantissa = TEN.mantissa; + tmp1.exponent = TEN.exponent; + tmp1.sign = TEN.sign; + } tmp1.pow(exp); } if (exp_neg) @@ -5499,20 +6418,20 @@ public final class Real if (a >= 5000000000000000000L) { // Rounding up gave 20 digits exponent++; a /= 5; - digits[18] = (byte)(a%10); + digits[18] = (byte) (a % 10); a /= 10; } else { - digits[18] = (byte)((a%5)*2); + digits[18] = (byte) ((a % 5) * 2); a /= 5; } } else { tmp2.round(); a = tmp2.toLong(); - digits[18] = (byte)(a%10); + digits[18] = (byte) (a % 10); a /= 10; } - for (int i=17; i>=0; i--) { - digits[i] = (byte)(a%10); + for (int i = 17; i >= 0; i--) { + digits[i] = (byte) (a % 10); a /= 10; } digits[19] = 0; @@ -5523,274 +6442,87 @@ public final class Real if ((this.exponent == 0 && this.mantissa == 0)) { sign = 0; // Two's complement cannot display -0 } else { - if ((this.sign!=0)) { + if ((this.sign != 0)) { mantissa = -mantissa; - if (((mantissa >> 62)&3) == 3) { + if (((mantissa >> 62) & 3) == 3) { mantissa <<= 1; exponent--; accurateBits--; // ? } } - exponent -= 0x40000000-1; - int shift = bitsPerDigit-1 - - floorMod(exponent, bitsPerDigit); + exponent -= 0x40000000 - 1; + int shift = bitsPerDigit - 1 - + floorMod(exponent, bitsPerDigit); exponent = floorDiv(exponent, bitsPerDigit); - if (shift == bitsPerDigit-1) { + 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)) { + } 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); + mantissa |= 0x8000000000000000L >> (shift - 1); } } } - int accurateDigits = (accurateBits+bitsPerDigit-1)/bitsPerDigit; - for (int i=0; i>>(64-bitsPerDigit)); + int accurateDigits = (accurateBits + bitsPerDigit - 1) / bitsPerDigit; + for (int i = 0; i < accurateDigits; i++) { + digits[i] = (byte) (mantissa >>> (64 - bitsPerDigit)); mantissa <<= bitsPerDigit; } digits[accurateDigits] = 0; return accurateDigits; } + private boolean carryWhenRounded(byte[] digits, int nDigits, int base) { - if (digits[nDigits] < base/2) + 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) + 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]++; + if (digits[nDigits] >= 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); + return align(ftoaBuf, format); } if ((this.exponent < 0 && this.mantissa == 0)) { - ftoaBuf.append((this.sign!=0) ? "-inf":"inf"); - return align(ftoaBuf,format); + ftoaBuf.append((this.sign != 0) ? "-inf" : "inf"); + return align(ftoaBuf, format); } int digitsPerThousand; switch (format.base) { @@ -5810,29 +6542,32 @@ public final class Real } if (format.thousand == 0) digitsPerThousand = 1000; // Disable thousands separator - { tmp4.mantissa = this.mantissa; tmp4.exponent = this.exponent; tmp4.sign = this.sign; }; + { + 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 + 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)) + else if ((tmp4.sign != 0)) width--; // subtract 1 for sign boolean useExp = false; switch (format.fse) { case NumberFormat.FSE_SCI: - precision = format.precision+1; + precision = format.precision + 1; useExp = true; break; case NumberFormat.FSE_ENG: - pointPos = floorMod(tmp4.exponent,3); - precision = format.precision+1+pointPos; + pointPos = floorMod(tmp4.exponent, 3); + precision = format.precision + 1 + pointPos; useExp = true; break; case NumberFormat.FSE_FIX: @@ -5840,54 +6575,53 @@ public final class Real 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) - { + 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){ + 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) + if (prefix != 0 && pointPos >= 0) width -= prefix; ftoaExp.setLength(0); if (useExp) { ftoaExp.append('e'); - ftoaExp.append(tmp4.exponent-pointPos); + 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 > 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); + 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) + 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); + char prefixChar = (format.base == 10 || (tmp4.sign == 0)) ? '0' : + hexChar.charAt(format.base - 1); if (pointPos < 0) { ftoaBuf.append(prefixChar); ftoaBuf.append(format.point); @@ -5897,9 +6631,9 @@ public final class Real } } // Add fractional part - for (int i=0; i0 && pointPos%digitsPerThousand==0) + if (pointPos > 0 && pointPos % digitsPerThousand == 0) ftoaBuf.append(format.thousand); if (pointPos == 0) ftoaBuf.append(format.point); @@ -5907,46 +6641,46 @@ public final class Real } if (format.fse == NumberFormat.FSE_NONE) { // Remove trailing zeros - while (ftoaBuf.charAt(ftoaBuf.length()-1) == '0') - ftoaBuf.setLength(ftoaBuf.length()-1); + 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); + 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() 0 && pointPos2 % digitsPerThousand == 0) + ftoaBuf.insert(0, format.thousand); + if (ftoaBuf.length() < format.maxwidth) + ftoaBuf.insert(0, prefixChar); } if (ftoaBuf.charAt(0) == format.thousand) ftoaBuf.deleteCharAt(0); } - return align(ftoaBuf,format); + return align(ftoaBuf, format); } - private static NumberFormat tmpFormat = new NumberFormat(); + /** * Converts this 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() + * this.toString() *
* Execution time relative to add:   * - * 130 + * 130 *
* * @return a String representation of this Real. @@ -5955,61 +6689,63 @@ public final class Real 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 + * this.toString() // Works only for base-10 *
* Execution time relative to add:   * base-2 - * 120 + * 120 *
base-8 - * 110 + * 110 *
base-10 - * 130 + * 130 *
base-16   - * 120 + * 120 *
* * @param base the base for the conversion. Valid base values are - * 2, 8, 10 and 16. + * 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 + * String.format("%...g",this); // Works only for base-10 *
* Execution time relative to add:   * base-2 - * 120 + * 120 *
base-8 - * 110 + * 110 *
base-10 - * 130 + * 130 *
base-16   - * 120 + * 120 *
* * @param format the number format to use in the conversion. @@ -6018,4 +6754,201 @@ public final class Real public String toString(NumberFormat format) { return ftoa(format); } + + /** + * 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 { + /** + * 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; + /** + * 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; + } }