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