All Downloads are FREE. Search and download functionalities are using the official Maven repository.

javolution.lang.MathLib Maven / Gradle / Ivy

Go to download

Only the Java Core part of Javolution library, with slight modifications for use in MSFTBX.

There is a newer version: 6.11.8
Show newest version
/*
 * Javolution - Java(TM) Solution for Real-Time and Embedded Systems
 * Copyright (C) 2012 - Javolution (http://javolution.org/)
 * All rights reserved.
 * 
 * Permission to use, copy, modify, and distribute this software is
 * freely granted, provided that this notice is preserved.
 */
package javolution.lang;


/**
 * 

An utility class providing a {@link Realtime} implementation of * the math library.

* * @author Jean-Marie Dautelle * @version 4.2, January 6, 2007 */ @Realtime public final class MathLib { /** * Default constructor. */ private MathLib() {} /** * Returns the number of bits in the minimal two's-complement representation * of the specified int, excluding a sign bit. * For positive int, this is equivalent to the number of bits * in the ordinary binary representation. For negative int, * it is equivalent to the number of bits of the positive value * -(i + 1). * * @param i the int value for which the bit length is returned. * @return the bit length of i. */ public static int bitLength(int i) { if (i < 0) i = -++i; return (i < 1 << 16) ? (i < 1 << 8) ? BIT_LENGTH[i] : BIT_LENGTH[i >>> 8] + 8 : (i < 1 << 24) ? BIT_LENGTH[i >>> 16] + 16 : BIT_LENGTH[i >>> 24] + 24; } private static final byte[] BIT_LENGTH = { 0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8 }; /** * Returns the number of bits in the minimal two's-complement representation * of the specified long, excluding a sign bit. * For positive long, this is equivalent to the number of bits * in the ordinary binary representation. For negative long, * it is equivalent to the number of bits of the positive value * -(l + 1). * * @param l the long value for which the bit length is returned. * @return the bit length of l. */ public static int bitLength(long l) { int i = (int) (l >> 32); if (i > 0) return (i < 1 << 16) ? (i < 1 << 8) ? BIT_LENGTH[i] + 32 : BIT_LENGTH[i >>> 8] + 40 : (i < 1 << 24) ? BIT_LENGTH[i >>> 16] + 48 : BIT_LENGTH[i >>> 24] + 56; if (i < 0) return bitLength(-++l); i = (int) l; return (i < 0) ? 32 : (i < 1 << 16) ? (i < 1 << 8) ? BIT_LENGTH[i] : BIT_LENGTH[i >>> 8] + 8 : (i < 1 << 24) ? BIT_LENGTH[i >>> 16] + 16 : BIT_LENGTH[i >>> 24] + 24; } /** * Returns the number of one-bits in the two's complement binary * representation of the specified long value. * This function is sometimes referred to as the population count. * * @param longValue the long value. * @return the number of one-bits in the two's complement binary * representation of the specified longValue. */ public static int bitCount(long longValue) { longValue = longValue - ((longValue >>> 1) & 0x5555555555555555L); longValue = (longValue & 0x3333333333333333L) + ((longValue >>> 2) & 0x3333333333333333L); longValue = (longValue + (longValue >>> 4)) & 0x0f0f0f0f0f0f0f0fL; longValue = longValue + (longValue >>> 8); longValue = longValue + (longValue >>> 16); longValue = longValue + (longValue >>> 32); return (int) longValue & 0x7f; } /** * Returns the number of zero bits preceding the highest-order * ("leftmost") one-bit in the two's complement binary representation * of the specified long value. Returns 64 if the specifed * value is zero. * * @param longValue the long value. * @return the number of leading zero bits. */ public static int numberOfLeadingZeros(long longValue) { // From Hacker's Delight if (longValue == 0) return 64; int n = 1; int x = (int)(longValue >>> 32); if (x == 0) { n += 32; x = (int)longValue; } if (x >>> 16 == 0) { n += 16; x <<= 16; } if (x >>> 24 == 0) { n += 8; x <<= 8; } if (x >>> 28 == 0) { n += 4; x <<= 4; } if (x >>> 30 == 0) { n += 2; x <<= 2; } n -= x >>> 31; return n; } /** * Returns the number of zero bits following the lowest-order ("rightmost") * one-bit in the two's complement binary representation of the specified * long value. Returns 64 if the specifed value is zero. * * @param longValue the long value. * @return the number of trailing zero bits. */ public static int numberOfTrailingZeros(long longValue) { // From Hacker's Delight int x, y; if (longValue == 0) return 64; int n = 63; y = (int)longValue; if (y != 0) { n = n -32; x = y; } else x = (int)(longValue>>>32); y = x <<16; if (y != 0) { n = n -16; x = y; } y = x << 8; if (y != 0) { n = n - 8; x = y; } y = x << 4; if (y != 0) { n = n - 4; x = y; } y = x << 2; if (y != 0) { n = n - 2; x = y; } return n - ((x << 1) >>> 31); } /** * Returns the number of digits of the decimal representation of the * specified int value, excluding the sign character if any. * * @param i the int value for which the digit length is returned. * @return String.valueOf(i).length() for zero or positive values; * String.valueOf(i).length() - 1 for negative values. */ public static int digitLength(int i) { if (i >= 0) return (i >= 100000) ? (i >= 10000000) ? (i >= 1000000000) ? 10 : (i >= 100000000) ? 9 : 8 : (i >= 1000000) ? 7 : 6 : (i >= 100) ? (i >= 10000) ? 5 : (i >= 1000) ? 4 : 3 : (i >= 10) ? 2 : 1; if (i == Integer.MIN_VALUE) return 10; // "2147483648".length() return digitLength(-i); // No overflow possible. } /** * Returns the number of digits of the decimal representation of the * the specified long, excluding the sign character if any. * * @param l the long value for which the digit length is returned. * @return String.valueOf(l).length() for zero or positive values; * String.valueOf(l).length() - 1 for negative values. */ public static int digitLength(long l) { if (l >= 0) return (l <= Integer.MAX_VALUE) ? digitLength((int) l) : // At least 10 digits or more. (l >= 100000000000000L) ? (l >= 10000000000000000L) ? (l >= 1000000000000000000L) ? 19 : (l >= 100000000000000000L) ? 18 : 17 : (l >= 1000000000000000L) ? 16 : 15 : (l >= 100000000000L) ? (l >= 10000000000000L) ? 14 : (l >= 1000000000000L) ? 13 : 12 : (l >= 10000000000L) ? 11 : 10; if (l == Long.MIN_VALUE) return 19; // "9223372036854775808".length() return digitLength(-l); } /** * Returns the closest double representation of the * specified long number multiplied by a power of two. * * @param m the long multiplier. * @param n the power of two exponent. * @return m * 2n. */ public static double toDoublePow2(long m, int n) { if (m == 0) return 0.0; if (m == Long.MIN_VALUE) return toDoublePow2(Long.MIN_VALUE >> 1, n + 1); if (m < 0) return -toDoublePow2(-m, n); int bitLength = MathLib.bitLength(m); int shift = bitLength - 53; long exp = 1023L + 52 + n + shift; // Use long to avoid overflow. if (exp >= 0x7FF) return Double.POSITIVE_INFINITY; if (exp <= 0) { // Degenerated number (subnormal, assume 0 for bit 52) if (exp <= -54) return 0.0; return toDoublePow2(m, n + 54) / 18014398509481984L; // 2^54 Exact. } // Normal number. long bits = (shift > 0) ? (m >> shift) + ((m >> (shift - 1)) & 1) : // Rounding. m << -shift; if (((bits >> 52) != 1) && (++exp >= 0x7FF)) return Double.POSITIVE_INFINITY; bits &= 0x000fffffffffffffL; // Clears MSB (bit 52) bits |= exp << 52; return Double.longBitsToDouble(bits); } /**/ /** * Returns the closest double representation of the * specified long number multiplied by a power of ten. * * @param m the long multiplier. * @param n the power of ten exponent. * @return multiplier * 10n. **/ public static double toDoublePow10(long m, int n) { if (m == 0) return 0.0; if (m == Long.MIN_VALUE) return toDoublePow10(Long.MIN_VALUE / 10, n + 1); if (m < 0) return -toDoublePow10(-m, n); if (n >= 0) { // Positive power. if (n > 308) return Double.POSITIVE_INFINITY; // Works with 4 x 32 bits registers (x3:x2:x1:x0) long x0 = 0; // 32 bits. long x1 = 0; // 32 bits. long x2 = m & MASK_32; // 32 bits. long x3 = m >>> 32; // 32 bits. int pow2 = 0; while (n != 0) { int i = (n >= POW5_INT.length) ? POW5_INT.length - 1 : n; int coef = POW5_INT[i]; // 31 bits max. if (((int) x0) != 0) x0 *= coef; // 63 bits max. if (((int) x1) != 0) x1 *= coef; // 63 bits max. x2 *= coef; // 63 bits max. x3 *= coef; // 63 bits max. x1 += x0 >>> 32; x0 &= MASK_32; x2 += x1 >>> 32; x1 &= MASK_32; x3 += x2 >>> 32; x2 &= MASK_32; // Adjusts powers. pow2 += i; n -= i; // Normalizes (x3 should be 32 bits max). long carry = x3 >>> 32; if (carry != 0) { // Shift. x0 = x1; x1 = x2; x2 = x3 & MASK_32; x3 = carry; pow2 += 32; } } // Merges registers to a 63 bits mantissa. int shift = 31 - MathLib.bitLength(x3); // -1..30 pow2 -= shift; long mantissa = (shift < 0) ? (x3 << 31) | (x2 >>> 1) : // x3 is 32 bits. (((x3 << 32) | x2) << shift) | (x1 >>> (32 - shift)); return toDoublePow2(mantissa, pow2); } else { // n < 0 if (n < -324 - 20) return 0.0; // Works with x1:x0 126 bits register. long x1 = m; // 63 bits. long x0 = 0; // 63 bits. int pow2 = 0; while (true) { // Normalizes x1:x0 int shift = 63 - MathLib.bitLength(x1); x1 <<= shift; x1 |= x0 >>> (63 - shift); x0 = (x0 << shift) & MASK_63; pow2 -= shift; // Checks if division has to be performed. if (n == 0) break; // Done. // Retrieves power of 5 divisor. int i = (-n >= POW5_INT.length) ? POW5_INT.length - 1 : -n; int divisor = POW5_INT[i]; // Performs the division (126 bits by 31 bits). long wh = (x1 >>> 32); long qh = wh / divisor; long r = wh - qh * divisor; long wl = (r << 32) | (x1 & MASK_32); long ql = wl / divisor; r = wl - ql * divisor; x1 = (qh << 32) | ql; wh = (r << 31) | (x0 >>> 32); qh = wh / divisor; r = wh - qh * divisor; wl = (r << 32) | (x0 & MASK_32); ql = wl / divisor; x0 = (qh << 32) | ql; // Adjusts powers. n += i; pow2 -= i; } return toDoublePow2(x1, pow2); } } private static final long MASK_63 = 0x7FFFFFFFFFFFFFFFL; private static final long MASK_32 = 0xFFFFFFFFL; private static final int[] POW5_INT = { 1, 5, 25, 125, 625, 3125, 15625, 78125, 390625, 1953125, 9765625, 48828125, 244140625, 1220703125 }; /**/ /** * Returns the closest long representation of the * specified double number multiplied by a power of two. * * @param d the double multiplier. * @param n the power of two exponent. * @return d * 2n * @throws ArithmeticException if the conversion cannot be performed * (NaN, Infinity or overflow). **/ public static long toLongPow2(double d, int n) { long bits = Double.doubleToLongBits(d); boolean isNegative = (bits >> 63) != 0; int exp = ((int) (bits >> 52)) & 0x7FF; long m = bits & 0x000fffffffffffffL; if (exp == 0x7FF) throw new ArithmeticException( "Cannot convert to long (Infinity or NaN)"); if (exp == 0) { if (m == 0) return 0L; return toLongPow2(d * 18014398509481984L, n - 54); // 2^54 Exact. } m |= 0x0010000000000000L; // Sets MSB (bit 52) long shift = exp - 1023L - 52 + n; // Use long to avoid overflow. if (shift <= -64) return 0L; if (shift >= 11) throw new ArithmeticException("Cannot convert to long (overflow)"); m = (shift >= 0) ? m << shift : (m >> -shift) + ((m >> -(shift + 1)) & 1); // Rounding. return isNegative ? -m : m; } /**/ /** * Returns the closest long representation of the * specified double number multiplied by a power of ten. * * @param d the double multiplier. * @param n the power of two exponent. * @return d * 10n. */ public static long toLongPow10(double d, int n) { long bits = Double.doubleToLongBits(d); boolean isNegative = (bits >> 63) != 0; int exp = ((int) (bits >> 52)) & 0x7FF; long m = bits & 0x000fffffffffffffL; if (exp == 0x7FF) throw new ArithmeticException( "Cannot convert to long (Infinity or NaN)"); if (exp == 0) { if (m == 0) return 0L; return toLongPow10(d * 1E16, n - 16); } m |= 0x0010000000000000L; // Sets MSB (bit 52) int pow2 = exp - 1023 - 52; // Retrieves 63 bits m with n == 0. if (n >= 0) { // Works with 4 x 32 bits registers (x3:x2:x1:x0) long x0 = 0; // 32 bits. long x1 = 0; // 32 bits. long x2 = m & MASK_32; // 32 bits. long x3 = m >>> 32; // 32 bits. while (n != 0) { int i = (n >= POW5_INT.length) ? POW5_INT.length - 1 : n; int coef = POW5_INT[i]; // 31 bits max. if (((int) x0) != 0) x0 *= coef; // 63 bits max. if (((int) x1) != 0) x1 *= coef; // 63 bits max. x2 *= coef; // 63 bits max. x3 *= coef; // 63 bits max. x1 += x0 >>> 32; x0 &= MASK_32; x2 += x1 >>> 32; x1 &= MASK_32; x3 += x2 >>> 32; x2 &= MASK_32; // Adjusts powers. pow2 += i; n -= i; // Normalizes (x3 should be 32 bits max). long carry = x3 >>> 32; if (carry != 0) { // Shift. x0 = x1; x1 = x2; x2 = x3 & MASK_32; x3 = carry; pow2 += 32; } } // Merges registers to a 63 bits mantissa. int shift = 31 - MathLib.bitLength(x3); // -1..30 pow2 -= shift; m = (shift < 0) ? (x3 << 31) | (x2 >>> 1) : // x3 is 32 bits. (((x3 << 32) | x2) << shift) | (x1 >>> (32 - shift)); } else { // n < 0 // Works with x1:x0 126 bits register. long x1 = m; // 63 bits. long x0 = 0; // 63 bits. while (true) { // Normalizes x1:x0 int shift = 63 - MathLib.bitLength(x1); x1 <<= shift; x1 |= x0 >>> (63 - shift); x0 = (x0 << shift) & MASK_63; pow2 -= shift; // Checks if division has to be performed. if (n == 0) break; // Done. // Retrieves power of 5 divisor. int i = (-n >= POW5_INT.length) ? POW5_INT.length - 1 : -n; int divisor = POW5_INT[i]; // Performs the division (126 bits by 31 bits). long wh = (x1 >>> 32); long qh = wh / divisor; long r = wh - qh * divisor; long wl = (r << 32) | (x1 & MASK_32); long ql = wl / divisor; r = wl - ql * divisor; x1 = (qh << 32) | ql; wh = (r << 31) | (x0 >>> 32); qh = wh / divisor; r = wh - qh * divisor; wl = (r << 32) | (x0 & MASK_32); ql = wl / divisor; x0 = (qh << 32) | ql; // Adjusts powers. n += i; pow2 -= i; } m = x1; } if (pow2 > 0) throw new ArithmeticException("Overflow"); if (pow2 < -63) return 0; m = (m >> -pow2) + ((m >> -(pow2 + 1)) & 1); // Rounding. return isNegative ? -m : m; } /**/ /** * Returns the largest power of 2 that is less than or equal to the * the specified positive value. * * @param d the double number. * @return floor(Log2(abs(d))) * @throws ArithmeticException if d <= 0 or d * is NaN or Infinity. **/ public static int floorLog2(double d) { if (d <= 0) throw new ArithmeticException("Negative number or zero"); long bits = Double.doubleToLongBits(d); int exp = ((int) (bits >> 52)) & 0x7FF; if (exp == 0x7FF) throw new ArithmeticException("Infinity or NaN"); if (exp == 0) return floorLog2(d * 18014398509481984L) - 54; // 2^54 Exact. return exp - 1023; } /**/ /** * Returns the largest power of 10 that is less than or equal to the * the specified positive value. * * @param d the double number. * @return floor(Log10(abs(d))) * @throws ArithmeticException if d <= 0 or d * is NaN or Infinity. **/ public static int floorLog10(double d) { int guess = (int) (LOG2_DIV_LOG10 * MathLib.floorLog2(d)); double pow10 = MathLib.toDoublePow10(1, guess); if ((pow10 <= d) && (pow10 * 10 > d)) return guess; if (pow10 > d) return guess - 1; return guess + 1; } private static final double LOG2_DIV_LOG10 = 0.3010299956639811952137388947; /** * The natural logarithm. **/ public static final double E = 2.71828182845904523536028747135266; /** * The ratio of the circumference of a circle to its diameter. **/ public static final double PI = 3.1415926535897932384626433832795; /** * Half the ratio of the circumference of a circle to its diameter. **/ public static final double HALF_PI = 1.5707963267948966192313216916398; /** * Twice the ratio of the circumference of a circle to its diameter. **/ public static final double TWO_PI = 6.283185307179586476925286766559; /** * Four time the ratio of the circumference of a circle to its diameter. **/ public static final double FOUR_PI = 12.566370614359172953850573533118; /** * Holds {@link #PI} * {@link #PI}. **/ public static final double PI_SQUARE = 9.8696044010893586188344909998762; /** * The natural logarithm of two. **/ public static final double LOG2 = 0.69314718055994530941723212145818; /** * The natural logarithm of ten. **/ public static final double LOG10 = 2.3025850929940456840179914546844; /** * The square root of two. **/ public static final double SQRT2 = 1.4142135623730950488016887242097; /** * Not-A-Number. **/ public static final double NaN = 0.0 / 0.0; /** * Infinity. **/ public static final double Infinity = 1.0 / 0.0; /**/ /** * Converts an angle in degrees to radians. * * @param degrees the angle in degrees. * @return the specified angle in radians. **/ public static double toRadians(double degrees) { return degrees * (PI / 180.0); } /**/ /** * Converts an angle in radians to degrees. * * @param radians the angle in radians. * @return the specified angle in degrees. **/ public static double toDegrees(double radians) { return radians * (180.0 / PI); } /**/ /** * Returns the positive square root of the specified value. * * @param x the value. * @return java.lang.Math.sqrt(x) **/ public static double sqrt(double x) { return Math.sqrt(x); // CLDC 1.1 } /**/ /** * Returns the remainder of the division of the specified two arguments. * * @param x the dividend. * @param y the divisor. * @return x - round(x / y) * y **/ public static double rem(double x, double y) { double tmp = x / y; if (MathLib.abs(tmp) <= Long.MAX_VALUE) return x - MathLib.round(tmp) * y; else return NaN; } /**/ /** * Returns the smallest (closest to negative infinity) * double value that is not less than the argument and is * equal to a mathematical integer. * * @param x the value. * @return java.lang.Math.ceil(x) **/ public static double ceil(double x) { return Math.ceil(x); // CLDC 1.1 } /**/ /** * Returns the largest (closest to positive infinity) * double value that is not greater than the argument and * is equal to a mathematical integer. * * @param x the value. * @return java.lang.Math.ceil(x) **/ public static double floor(double x) { return Math.floor(x); // CLDC 1.1 } /**/ /** * Returns the trigonometric sine of the specified angle in radians. * * @param radians the angle in radians. * @return java.lang.Math.sin(radians) **/ public static double sin(double radians) { return Math.sin(radians); // CLDC 1.1 } /**/ /** * Returns the trigonometric cosine of the specified angle in radians. * * @param radians the angle in radians. * @return java.lang.Math.cos(radians) **/ public static double cos(double radians) { return Math.cos(radians); // CLDC 1.1 } /**/ /** * Returns the trigonometric tangent of the specified angle in radians. * * @param radians the angle in radians. * @return java.lang.Math.tan(radians) **/ public static double tan(double radians) { return Math.tan(radians); // CLDC 1.1 } /**/ /** * Returns the arc sine of the specified value, * in the range of -pi/2 through pi/2. * * @param x the value whose arc sine is to be returned. * @return the arc sine in radians for the specified value. **/ public static double asin(double x) { if (x < -1.0 || x > 1.0) return MathLib.NaN; if (x == -1.0) return -HALF_PI; if (x == 1.0) return HALF_PI; return MathLib.atan(x / MathLib.sqrt(1.0 - x * x)); } /**/ /** * Returns the arc cosine of the specified value, * in the range of 0.0 through pi. * * @param x the value whose arc cosine is to be returned. * @return the arc cosine in radians for the specified value. **/ public static double acos(double x) { return HALF_PI - MathLib.asin(x); } /**/ /** * Returns the arc tangent of the specified value, * in the range of -pi/2 through pi/2. * * @param x the value whose arc tangent is to be returned. * @return the arc tangent in radians for the specified value. * @see * Inverse Tangent -- from MathWorld **/ public static double atan(double x) { return MathLib._atan(x); } /**/ /** * Returns the angle theta such that * (x == cos(theta)) && (y == sin(theta)). * * @param y the y value. * @param x the x value. * @return the angle theta in radians. * @see Wikipedia: Atan2 **/ public static double atan2(double y, double x) { // From Wikipedia. if (x > 0) return MathLib.atan(y / x); if ((y >= 0) && (x < 0)) return MathLib.atan(y / x) + PI; if ((y < 0) && (x < 0)) return MathLib.atan(y / x) - PI; if ((y > 0) && (x == 0)) return PI / 2; if ((y < 0) && (x == 0)) return -PI / 2; return Double.NaN; // ((y == 0) && (x == 0)) } /**/ /** * Returns the hyperbolic sine of x. * * @param x the value for which the hyperbolic sine is calculated. * @return (exp(x) - exp(-x)) / 2 **/ public static double sinh(double x) { return (MathLib.exp(x) - MathLib.exp(-x)) * 0.5; } /**/ /** * Returns the hyperbolic cosine of x. * * @param x the value for which the hyperbolic cosine is calculated. * @return (exp(x) + exp(-x)) / 2 **/ public static double cosh(double x) { return (MathLib.exp(x) + MathLib.exp(-x)) * 0.5; } /**/ /** * Returns the hyperbolic tangent of x. * * @param x the value for which the hyperbolic tangent is calculated. * @return (exp(2 * x) - 1) / (exp(2 * x) + 1) **/ public static double tanh(double x) { return (MathLib.exp(2 * x) - 1) / (MathLib.exp(2 * x) + 1); } /**/ /** * Returns {@link #E e} raised to the specified power. * * @param x the exponent. * @return ex * @see * Exponential Function -- from MathWorld **/ public static double exp(double x) { return MathLib._ieee754_exp(x); } /**/ /** * Returns the natural logarithm (base {@link #E e}) of the specified * value. * * @param x the value greater than 0.0. * @return the value y such as ey == x **/ public static double log(double x) { return MathLib._ieee754_log(x); } /**/ /** * Returns the decimal logarithm of the specified value. * * @param x the value greater than 0.0. * @return the value y such as 10y == x **/ public static double log10(double x) { return log(x) * INV_LOG10; } private static double INV_LOG10 = 0.43429448190325182765112891891661; /** * Returns the value of the first argument raised to the power of the * second argument. * * @param x the base. * @param y the exponent. * @return xy **/ public static double pow(double x, double y) { // Use close approximation (+/- LSB) if ((x < 0) && (y == (int) y)) return (((int) y) & 1) == 0 ? pow(-x, y) : -pow(-x, y); return MathLib.exp(y * MathLib.log(x)); } /** * Returns the closest int to the specified argument. * * @param f the float value to be rounded to a int * @return the nearest int value. **/ public static int round(float f) { return (int) floor(f + 0.5f); } /**/ /** * Returns the closest long to the specified argument. * * @param d the double value to be rounded to a * long * @return the nearest long value. **/ public static long round(double d) { return (long) floor(d + 0.5d); } /** * Returns the absolute value of the specified int argument. * * @param i the int value. * @return i or -i */ public static int abs(int i) { return (i < 0) ? -i : i; } /** * Returns the absolute value of the specified long argument. * * @param l the long value. * @return l or -l */ public static long abs(long l) { return (l < 0) ? -l : l; } /** * Returns the absolute value of the specified float argument. * * @param f the float value. * @return f or -f **/ public static float abs(float f) { return (f < 0) ? -f : f; } /** * Returns the absolute value of the specified double argument. * * @param d the double value. * @return d or -d **/ public static double abs(double d) { return (d < 0) ? -d : d; } /**/ /** * Returns the greater of two int values. * * @param x the first value. * @param y the second value. * @return the larger of x and y. */ public static int max(int x, int y) { return (x >= y) ? x : y; } /** * Returns the greater of two long values. * * @param x the first value. * @param y the second value. * @return the larger of x and y. */ public static long max(long x, long y) { return (x >= y) ? x : y; } /** * Returns the greater of two float values. * * @param x the first value. * @param y the second value. * @return the larger of x and y. **/ public static float max(float x, float y) { return (x >= y) ? x : y; } /** * Returns the greater of two double values. * * @param x the first value. * @param y the second value. * @return the larger of x and y. **/ public static double max(double x, double y) { return (x >= y) ? x : y; } /**/ /** * Returns the smaller of two int values. * * @param x the first value. * @param y the second value. * @return the smaller of x and y. */ public static int min(int x, int y) { return (x < y) ? x : y; } /** * Returns the smaller of two long values. * * @param x the first value. * @param y the second value. * @return the smaller of x and y. */ public static long min(long x, long y) { return (x < y) ? x : y; } /** * Returns the smaller of two float values. * * @param x the first value. * @param y the second value. * @return the smaller of x and y. **/ public static float min(float x, float y) { return (x < y) ? x : y; } /**/ /** * Returns the smaller of two double values. * * @param x the first value. * @param y the second value. * @return the smaller of x and y. **/ public static double min(double x, double y) { return (x < y) ? x : y; } //////////////////////////////////////////////////////////////////////////// /* @(#)s_atan.c 1.3 95/01/18 */ /* * ==================================================== * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * * Developed at SunSoft, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this * software is freely granted, provided that this notice * is preserved. * ==================================================== * */ /* atan(x) * Method * 1. Reduce x to positive by atan(x) = -atan(-x). * 2. According to the integer k=4t+0.25 chopped, t=x, the argument * is further reduced to one of the following intervals and the * arctangent of t is evaluated by the corresponding formula: * * [0,7/16] atan(x) = t-t^3*(a1+t^2*(a2+...(a10+t^2*a11)...) * [7/16,11/16] atan(x) = atan(1/2) + atan( (t-0.5)/(1+t/2) ) * [11/16.19/16] atan(x) = atan( 1 ) + atan( (t-1)/(1+t) ) * [19/16,39/16] atan(x) = atan(3/2) + atan( (t-1.5)/(1+1.5t) ) * [39/16,INF] atan(x) = atan(INF) + atan( -1/t ) * * Constants: * The hexadecimal values are the intended ones for the following * constants. The decimal values may be used, provided that the * compiler will convert from decimal to binary accurately enough * to produce the hexadecimal values shown. */ static final double atanhi[] = { 4.63647609000806093515e-01, // atan(0.5)hi 0x3FDDAC67, 0x0561BB4F 7.85398163397448278999e-01, // atan(1.0)hi 0x3FE921FB, 0x54442D18 9.82793723247329054082e-01, // atan(1.5)hi 0x3FEF730B, 0xD281F69B 1.57079632679489655800e+00, // atan(inf)hi 0x3FF921FB, 0x54442D18 }; static final double atanlo[] = { 2.26987774529616870924e-17, // atan(0.5)lo 0x3C7A2B7F, 0x222F65E2 3.06161699786838301793e-17, // atan(1.0)lo 0x3C81A626, 0x33145C07 1.39033110312309984516e-17, // atan(1.5)lo 0x3C700788, 0x7AF0CBBD 6.12323399573676603587e-17, // atan(inf)lo 0x3C91A626, 0x33145C07 }; static final double aT[] = { 3.33333333333329318027e-01, // 0x3FD55555, 0x5555550D -1.99999999998764832476e-01, // 0xBFC99999, 0x9998EBC4 1.42857142725034663711e-01, // 0x3FC24924, 0x920083FF -1.11111104054623557880e-01, // 0xBFBC71C6, 0xFE231671 9.09088713343650656196e-02, // 0x3FB745CD, 0xC54C206E -7.69187620504482999495e-02, // 0xBFB3B0F2, 0xAF749A6D 6.66107313738753120669e-02, // 0x3FB10D66, 0xA0D03D51 -5.83357013379057348645e-02, // 0xBFADDE2D, 0x52DEFD9A 4.97687799461593236017e-02, // 0x3FA97B4B, 0x24760DEB -3.65315727442169155270e-02, // 0xBFA2B444, 0x2C6A6C2F 1.62858201153657823623e-02, // 0x3F90AD3A, 0xE322DA11 }; static final double one = 1.0, huge = 1.0e300; static double _atan(double x) { double w, s1, s2, z; int ix, hx, id; long xBits = Double.doubleToLongBits(x); int __HIx = (int) (xBits >> 32); int __LOx = (int) xBits; hx = __HIx; ix = hx & 0x7fffffff; if (ix >= 0x44100000) { // if |x| >= 2^66 if (ix > 0x7ff00000 || (ix == 0x7ff00000 && (__LOx != 0))) return x + x; // NaN if (hx > 0) return atanhi[3] + atanlo[3]; else return -atanhi[3] - atanlo[3]; } if (ix < 0x3fdc0000) { // |x| < 0.4375 if (ix < 0x3e200000) // |x| < 2^-29 if (huge + x > one) return x; id = -1; } else { x = MathLib.abs(x); if (ix < 0x3ff30000) // |x| < 1.1875 if (ix < 0x3fe60000) { // 7/16 <=|x|<11/16 id = 0; x = (2.0 * x - one) / (2.0 + x); } else { // 11/16<=|x|< 19/16 id = 1; x = (x - one) / (x + one); } else if (ix < 0x40038000) { // |x| < 2.4375 id = 2; x = (x - 1.5) / (one + 1.5 * x); } else { // 2.4375 <= |x| < 2^66 id = 3; x = -1.0 / x; } } // end of argument reduction z = x * x; w = z * z; // break sum from i=0 to 10 aT[i]z**(i+1) into odd and even poly s1 = z * (aT[0] + w * (aT[2] + w * (aT[4] + w * (aT[6] + w * (aT[8] + w * aT[10]))))); s2 = w * (aT[1] + w * (aT[3] + w * (aT[5] + w * (aT[7] + w * aT[9])))); if (id < 0) return x - x * (s1 + s2); else { z = atanhi[id] - ((x * (s1 + s2) - atanlo[id]) - x); return (hx < 0) ? -z : z; } } /**/ //////////////////////////////////////////////////////////////////////////// /* @(#)e_log.c 1.3 95/01/18 */ /* * ==================================================== * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * * Developed at SunSoft, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this * software is freely granted, provided that this notice * is preserved. * ==================================================== */ /* __ieee754_log(x) * Return the logrithm of x * * Method : * 1. Argument Reduction: find k and f such that * x = 2^k * (1+f), * where sqrt(2)/2 < 1+f < sqrt(2) . * * 2. Approximation of log(1+f). * Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s) * = 2s + 2/3 s**3 + 2/5 s**5 + ....., * = 2s + s*R * We use a special Reme algorithm on [0,0.1716] to generate * a polynomial of degree 14 to approximate R The maximum error * of this polynomial approximation is bounded by 2**-58.45. In * other words, * 2 4 6 8 10 12 14 * R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s +Lg6*s +Lg7*s * (the values of Lg1 to Lg7 are listed in the program) * and * | 2 14 | -58.45 * | Lg1*s +...+Lg7*s - R(z) | <= 2 * | | * Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2. * In order to guarantee error in log below 1ulp, we compute log * by * log(1+f) = f - s*(f - R) (if f is not too large) * log(1+f) = f - (hfsq - s*(hfsq+R)). (better accuracy) * * 3. Finally, log(x) = k*ln2 + log(1+f). * = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo))) * Here ln2 is split into two floating point number: * ln2_hi + ln2_lo, * where n*ln2_hi is always exact for |n| < 2000. * * Special cases: * log(x) is NaN with signal if x < 0 (including -INF) ; * log(+INF) is +INF; log(0) is -INF with signal; * log(NaN) is that NaN with no signal. * * Accuracy: * according to an error analysis, the error is always less than * 1 ulp (unit in the last place). * * Constants: * The hexadecimal values are the intended ones for the following * constants. The decimal values may be used, provided that the * compiler will convert from decimal to binary accurately enough * to produce the hexadecimal values shown. */ static final double ln2_hi = 6.93147180369123816490e-01, // 3fe62e42 fee00000 ln2_lo = 1.90821492927058770002e-10, // 3dea39ef 35793c76 two54 = 1.80143985094819840000e+16, // 43500000 00000000 Lg1 = 6.666666666666735130e-01, // 3FE55555 55555593 Lg2 = 3.999999999940941908e-01, // 3FD99999 9997FA04 Lg3 = 2.857142874366239149e-01, // 3FD24924 94229359 Lg4 = 2.222219843214978396e-01, // 3FCC71C5 1D8E78AF Lg5 = 1.818357216161805012e-01, // 3FC74664 96CB03DE Lg6 = 1.531383769920937332e-01, // 3FC39A09 D078C69F Lg7 = 1.479819860511658591e-01; // 3FC2F112 DF3E5244 static final double zero = 0.0; static double _ieee754_log(double x) { double hfsq, f, s, z, R, w, t1, t2, dk; int k, hx, i, j; int lx; // unsigned long xBits = Double.doubleToLongBits(x); hx = (int) (xBits >> 32); lx = (int) xBits; k = 0; if (hx < 0x00100000) { // x < 2**-1022 if (((hx & 0x7fffffff) | lx) == 0) return -two54 / zero; // log(+-0)=-inf if (hx < 0) return (x - x) / zero; // log(-#) = NaN k -= 54; x *= two54; // subnormal number, scale up x xBits = Double.doubleToLongBits(x); hx = (int) (xBits >> 32); // high word of x } if (hx >= 0x7ff00000) return x + x; k += (hx >> 20) - 1023; hx &= 0x000fffff; i = (hx + 0x95f64) & 0x100000; xBits = Double.doubleToLongBits(x); int HIx = hx | (i ^ 0x3ff00000); // normalize x or x/2 xBits = ((HIx & 0xFFFFFFFFL) << 32) | (xBits & 0xFFFFFFFFL); x = Double.longBitsToDouble(xBits); k += (i >> 20); f = x - 1.0; if ((0x000fffff & (2 + hx)) < 3) { // |f| < 2**-20 if (f == zero) if (k == 0) return zero; else { dk = (double) k; return dk * ln2_hi + dk * ln2_lo; } R = f * f * (0.5 - 0.33333333333333333 * f); if (k == 0) return f - R; else { dk = (double) k; return dk * ln2_hi - ((R - dk * ln2_lo) - f); } } s = f / (2.0 + f); dk = (double) k; z = s * s; i = hx - 0x6147a; w = z * z; j = 0x6b851 - hx; t1 = w * (Lg2 + w * (Lg4 + w * Lg6)); t2 = z * (Lg1 + w * (Lg3 + w * (Lg5 + w * Lg7))); i |= j; R = t2 + t1; if (i > 0) { hfsq = 0.5 * f * f; if (k == 0) return f - (hfsq - s * (hfsq + R)); else return dk * ln2_hi - ((hfsq - (s * (hfsq + R) + dk * ln2_lo)) - f); } else if (k == 0) return f - s * (f - R); else return dk * ln2_hi - ((s * (f - R) - dk * ln2_lo) - f); } /**/ //////////////////////////////////////////////////////////////////////////// /* @(#)e_exp.c 1.6 04/04/22 */ /* * ==================================================== * Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved. * * Permission to use, copy, modify, and distribute this * software is freely granted, provided that this notice * is preserved. * ==================================================== */ /* __ieee754_exp(x) * Returns the exponential of x. * * Method * 1. Argument reduction: * Reduce x to an r so that |r| <= 0.5*ln2 ~ 0.34658. * Given x, find r and integer k such that * * x = k*ln2 + r, |r| <= 0.5*ln2. * * Here r will be represented as r = hi-lo for better * accuracy. * * 2. Approximation of exp(r) by a special rational function on * the interval [0,0.34658]: * Write * R(r**2) = r*(exp(r)+1)/(exp(r)-1) = 2 + r*r/6 - r**4/360 + ... * We use a special Remes algorithm on [0,0.34658] to generate * a polynomial of degree 5 to approximate R. The maximum error * of this polynomial approximation is bounded by 2**-59. In * other words, * R(z) ~ 2.0 + P1*z + P2*z**2 + P3*z**3 + P4*z**4 + P5*z**5 * (where z=r*r, and the values of P1 to P5 are listed below) * and * | 5 | -59 * | 2.0+P1*z+...+P5*z - R(z) | <= 2 * | | * The computation of exp(r) thus becomes * 2*r * exp(r) = 1 + ------- * R - r * r*R1(r) * = 1 + r + ----------- (for better accuracy) * 2 - R1(r) * where * 2 4 10 * R1(r) = r - (P1*r + P2*r + ... + P5*r ). * * 3. Scale back to obtain exp(x): * From step 1, we have * exp(x) = 2^k * exp(r) * * Special cases: * exp(INF) is INF, exp(NaN) is NaN; * exp(-INF) is 0, and * for finite argument, only exp(0)=1 is exact. * * Accuracy: * according to an error analysis, the error is always less than * 1 ulp (unit in the last place). * * Misc. info. * For IEEE double * if x > 7.09782712893383973096e+02 then exp(x) overflow * if x < -7.45133219101941108420e+02 then exp(x) underflow * * Constants: * The hexadecimal values are the intended ones for the following * constants. The decimal values may be used, provided that the * compiler will convert from decimal to binary accurately enough * to produce the hexadecimal values shown. */ static final double halF[] = { 0.5, -0.5, }, twom1000 = 9.33263618503218878990e-302, // 2**-1000=0x01700000,0 o_threshold = 7.09782712893383973096e+02, // 0x40862E42, 0xFEFA39EF u_threshold = -7.45133219101941108420e+02, // 0xc0874910, 0xD52D3051 ln2HI[] = { 6.93147180369123816490e-01, // 0x3fe62e42, 0xfee00000 -6.93147180369123816490e-01, },// 0xbfe62e42, 0xfee00000 ln2LO[] = { 1.90821492927058770002e-10, // 0x3dea39ef, 0x35793c76 -1.90821492927058770002e-10, },// 0xbdea39ef, 0x35793c76 invln2 = 1.44269504088896338700e+00, // 0x3ff71547, 0x652b82fe P1 = 1.66666666666666019037e-01, // 0x3FC55555, 0x5555553E P2 = -2.77777777770155933842e-03, // 0xBF66C16C, 0x16BEBD93 P3 = 6.61375632143793436117e-05, // 0x3F11566A, 0xAF25DE2C P4 = -1.65339022054652515390e-06, // 0xBEBBBD41, 0xC5D26BF1 P5 = 4.13813679705723846039e-08; // 0x3E663769, 0x72BEA4D0 static double _ieee754_exp(double x) // default IEEE double exp { double y, hi = 0, lo = 0, c, t; int k = 0, xsb; int hx; // Unsigned. long xBits = Double.doubleToLongBits(x); int __HIx = (int) (xBits >> 32); int __LOx = (int) xBits; hx = __HIx; // high word of x xsb = (hx >> 31) & 1; // sign bit of x hx &= 0x7fffffff; // high word of |x| // filter out non-finite argument if (hx >= 0x40862E42) { // if |x|>=709.78... if (hx >= 0x7ff00000) if (((hx & 0xfffff) | __LOx) != 0) return x + x; // NaN else return (xsb == 0) ? x : 0.0; if (x > o_threshold) return huge * huge; // overflow if (x < u_threshold) return twom1000 * twom1000; // underflow } // argument reduction if (hx > 0x3fd62e42) { // if |x| > 0.5 ln2 if (hx < 0x3FF0A2B2) { // and |x| < 1.5 ln2 hi = x - ln2HI[xsb]; lo = ln2LO[xsb]; k = 1 - xsb - xsb; } else { k = (int) (invln2 * x + halF[xsb]); t = k; hi = x - t * ln2HI[0]; // t*ln2HI is exact here lo = t * ln2LO[0]; } x = hi - lo; } else if (hx < 0x3e300000) { // when |x|<2**-28 if (huge + x > one) return one + x;// trigger inexact } else k = 0; // x is now in primary range t = x * x; c = x - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t * P5)))); if (k == 0) return one - ((x * c) / (c - 2.0) - x); else y = one - ((lo - (x * c) / (2.0 - c)) - hi); long yBits = Double.doubleToLongBits(y); int __HIy = (int) (yBits >> 32); if (k >= -1021) { __HIy += (k << 20); // add k to y's exponent yBits = ((__HIy & 0xFFFFFFFFL) << 32) | (yBits & 0xFFFFFFFFL); y = Double.longBitsToDouble(yBits); return y; } else { __HIy += ((k + 1000) << 20);// add k to y's exponent yBits = ((__HIy & 0xFFFFFFFFL) << 32) | (yBits & 0xFFFFFFFFL); y = Double.longBitsToDouble(yBits); return y * twom1000; } } }