0001: /*
0002: * Javolution - Java(TM) Solution for Real-Time and Embedded Systems
0003: * Copyright (C) 2006 - Javolution (http://javolution.org/)
0004: * All rights reserved.
0005: *
0006: * Permission to use, copy, modify, and distribute this software is
0007: * freely granted, provided that this notice is preserved.
0008: */
0009: package javolution.lang;
0010:
0011: import java.util.Random;
0012:
0013: /**
0014: * <p> This utility class ensures cross-platform portability of the math
0015: * library. Functions not supported by the platform are emulated.
0016: * Developers may replace the current implementation with native
0017: * implementations for faster execution (unlike
0018: * <code>java.lang.Math</code> this class can be customized).<p>
0019: *
0020: * @author <a href="mailto:jean-marie@dautelle.com">Jean-Marie Dautelle</a>
0021: * @version 4.2, January 6, 2007
0022: */
0023: public final class MathLib {
0024:
0025: /**
0026: * Default constructor.
0027: */
0028: private MathLib() {
0029: }
0030:
0031: /**
0032: * Returns a pseudo random <code>int</code> value in the range
0033: * <code>[min, max]</code> (fast and thread-safe without synchronization).
0034: *
0035: * @param min the minimum value inclusive.
0036: * @param max the maximum value exclusive.
0037: * @return a pseudo random number in the range <code>[min, max]</code>.
0038: */
0039: public static int random(int min, int max) {
0040: int next = RANDOM.nextInt();
0041: if ((next >= min) && (next <= max))
0042: return next;
0043: next += Integer.MIN_VALUE;
0044: if ((next >= min) && (next <= max))
0045: return next;
0046: // We should not have interval overflow as the interval has to be less
0047: // or equal to Integer.MAX_VALUE (otherwise we would have exited before).
0048: final int interval = 1 + max - min; // Positive.
0049: if (interval <= 0)
0050: throw new Error("Interval [" + min + ".." + max + "] error"); // In case.
0051: return MathLib.abs(next % interval) + min;
0052: }
0053:
0054: private static final Random RANDOM = new Random();
0055:
0056: /**
0057: * Returns a pseudo random <code>long</code> value in the range
0058: * <code>[min, max]</code> (fast and thread-safe without synchronization).
0059: *
0060: * @param min the minimum value inclusive.
0061: * @param max the maximum value inclusive.
0062: * @return a pseudo random number in the range <code>[min, max]</code>.
0063: */
0064: public static long random(long min, long max) {
0065: long next = RANDOM.nextLong();
0066: if ((next >= min) && (next <= max))
0067: return next;
0068: next += Long.MIN_VALUE;
0069: if ((next >= min) && (next <= max))
0070: return next;
0071: // We should not have interval overflow as the interval has to be less
0072: // or equal to Long.MAX_VALUE (otherwise we would have exited before).
0073: final long interval = 1L + max - min; // Positive.
0074: if (interval <= 0)
0075: throw new Error("Interval error"); // In case.
0076: return MathLib.abs(next % interval) + min;
0077: }
0078:
0079: /**
0080: * Returns a pseudo random <code>float</code> value in the range
0081: * <code>[min, max]</code> (fast and thread-safe without synchronization).
0082: *
0083: * @param min the minimum value inclusive.
0084: * @param max the maximum value inclusive.
0085: * @return a pseudo random number in the range <code>[min, max]</code>.
0086: * @JVM-1.1+@
0087: public static float random(float min, float max) {
0088: return (float) random((double)min, (double) max);
0089: }
0090: /**/
0091:
0092: /**
0093: * Returns a pseudo random <code>double</code> value in the range
0094: * <code>[min, max]</code> (fast and thread-safe without synchronization).
0095: *
0096: * @param min the minimum value inclusive.
0097: * @param max the maximum value inclusive.
0098: * @return a pseudo random number in the range <code>[min, max]</code>.
0099: * @JVM-1.1+@
0100: public static double random(double min, double max) {
0101: double next = RANDOM.nextDouble(); // 0.0 .. 1.0
0102: return min + (next * max) - (next * min);
0103: // Due to potential numeric errors, both min and max are possible.
0104: }
0105: /**/
0106:
0107: /**
0108: * Returns the number of bits in the minimal two's-complement representation
0109: * of the specified <code>int</code>, excluding a sign bit.
0110: * For positive <code>int</code>, this is equivalent to the number of bits
0111: * in the ordinary binary representation. For negative <code>int</code>,
0112: * it is equivalent to the number of bits of the positive value
0113: * <code>-(i + 1)</code>.
0114: *
0115: * @param i the <code>int</code> value for which the bit length is returned.
0116: * @return the bit length of <code>i</code>.
0117: */
0118: public static int bitLength(int i) {
0119: if (i < 0)
0120: i = -++i;
0121: return (i < 1 << 16) ? (i < 1 << 8) ? BIT_LENGTH[i]
0122: : BIT_LENGTH[i >>> 8] + 8
0123: : (i < 1 << 24) ? BIT_LENGTH[i >>> 16] + 16
0124: : BIT_LENGTH[i >>> 24] + 24;
0125: }
0126:
0127: private static final byte[] BIT_LENGTH = { 0, 1, 2, 2, 3, 3, 3, 3,
0128: 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
0129: 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
0130: 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7,
0131: 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
0132: 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
0133: 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
0134: 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
0135: 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
0136: 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
0137: 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
0138: 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
0139: 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
0140: 8, 8, 8, 8, 8, 8, 8, 8 };
0141:
0142: /**
0143: * Returns the number of bits in the minimal two's-complement representation
0144: * of the specified <code>long</code>, excluding a sign bit.
0145: * For positive <code>long</code>, this is equivalent to the number of bits
0146: * in the ordinary binary representation. For negative <code>long</code>,
0147: * it is equivalent to the number of bits of the positive value
0148: * <code>-(l + 1)</code>.
0149: *
0150: * @param l the <code>long</code> value for which the bit length is returned.
0151: * @return the bit length of <code>l</code>.
0152: */
0153: public static int bitLength(long l) {
0154: int i = (int) (l >> 32);
0155: if (i > 0)
0156: return (i < 1 << 16) ? (i < 1 << 8) ? BIT_LENGTH[i] + 32
0157: : BIT_LENGTH[i >>> 8] + 40
0158: : (i < 1 << 24) ? BIT_LENGTH[i >>> 16] + 48
0159: : BIT_LENGTH[i >>> 24] + 56;
0160: if (i < 0)
0161: return bitLength(-++l);
0162: i = (int) l;
0163: return (i < 0) ? 32
0164: : (i < 1 << 16) ? (i < 1 << 8) ? BIT_LENGTH[i]
0165: : BIT_LENGTH[i >>> 8] + 8
0166: : (i < 1 << 24) ? BIT_LENGTH[i >>> 16] + 16
0167: : BIT_LENGTH[i >>> 24] + 24;
0168: }
0169:
0170: /**
0171: * Returns the number of digits in the minimal ten's-complement
0172: * representation of the specified <code>int</code>, excluding a sign bit.
0173: * For positive <code>int</code>, this is equivalent to the number of digit
0174: * in the ordinary decimal representation. For negative <code>int</code>,
0175: * it is equivalent to the number of digit of the positive value
0176: * <code>-(i + 1)</code>.
0177: *
0178: * @param i the <code>int</code> value for which the digit length is returned.
0179: * @return the digit length of <code>i</code>.
0180: */
0181: public static int digitLength(int i) {
0182: if (i < 0)
0183: i = -++i;
0184: return (i >= 100000) ? (i >= 10000000) ? (i >= 1000000000) ? 10
0185: : (i >= 100000000) ? 9 : 8 : (i >= 1000000) ? 7 : 6
0186: : (i >= 100) ? (i >= 10000) ? 5 : (i >= 1000) ? 4 : 3
0187: : (i >= 10) ? 2 : (i >= 1) ? 1 : 0;
0188: }
0189:
0190: /**
0191: * Returns the number of digits in the minimal ten's-complement
0192: * representation of the specified <code>long</code>, excluding a sign bit.
0193: * For positive <code>long</code>, this is equivalent to the number of digit
0194: * in the ordinary decimal representation. For negative <code>long</code>,
0195: * it is equivalent to the number of digit of the positive value
0196: * <code>-(value + 1)</code>.
0197: *
0198: * @param l the <code>long</code> value for which the digit length is returned.
0199: * @return the digit length of <code>l</code>.
0200: */
0201: public static int digitLength(long l) {
0202: if (l < 0)
0203: l = -++l;
0204: if (l <= Integer.MAX_VALUE) {
0205: int i = (int) l;
0206: return (i >= 100000) ? (i >= 10000000) ? (i >= 1000000000) ? 10
0207: : (i >= 100000000) ? 9 : 8
0208: : (i >= 1000000) ? 7 : 6
0209: : (i >= 100) ? (i >= 10000) ? 5 : (i >= 1000) ? 4
0210: : 3 : (i >= 10) ? 2 : (i >= 1) ? 1 : 0;
0211: }
0212: // At least 10 digits or more.
0213: return (l >= 100000000000000L) ? (l >= 10000000000000000L) ? (l >= 1000000000000000000L) ? 19
0214: : (l >= 100000000000000000L) ? 18 : 17
0215: : (l >= 1000000000000000L) ? 16 : 15
0216: : (l >= 100000000000L) ? (l >= 10000000000000L) ? 14
0217: : (l >= 1000000000000L) ? 13 : 12
0218: : (l >= 10000000000L) ? 11 : 10;
0219:
0220: }
0221:
0222: /**
0223: * Returns the closest <code>double</code> representation of the
0224: * specified <code>long</code> number multiplied by a power of two.
0225: *
0226: * @param m the <code>long</code> multiplier.
0227: * @param n the power of two exponent.
0228: * @return <code>m * 2<sup>n</sup></code>.
0229: @JVM-1.1+@
0230: public static double toDoublePow2(long m, int n) {
0231: if (m == 0)
0232: return 0.0;
0233: if (m == Long.MIN_VALUE)
0234: return toDoublePow2(Long.MIN_VALUE >> 1, n + 1);
0235: if (m < 0)
0236: return -toDoublePow2(-m, n);
0237: int bitLength = MathLib.bitLength(m);
0238: int shift = bitLength - 53;
0239: long exp = 1023L + 52 + n + shift; // Use long to avoid overflow.
0240: if (exp >= 0x7FF)
0241: return Double.POSITIVE_INFINITY;
0242: if (exp <= 0) { // Degenerated number (subnormal, assume 0 for bit 52)
0243: if (exp <= -54)
0244: return 0.0;
0245: return toDoublePow2(m, n + 54) / 18014398509481984L; // 2^54 Exact.
0246: }
0247: // Normal number.
0248: long bits = (shift > 0) ? (m >> shift)
0249: + ((m >> (shift - 1)) & 1) : // Rounding.
0250: m << -shift;
0251: if (((bits >> 52) != 1) && (++exp >= 0x7FF)) // Rare case where rounding push to next power of 2.
0252: return Double.POSITIVE_INFINITY;
0253: bits &= 0x000fffffffffffffL; // Clears MSB (bit 52)
0254: bits |= exp << 52;
0255: return Double.longBitsToDouble(bits);
0256: }
0257: /**/
0258:
0259: /**
0260: * Returns the closest <code>double</code> representation of the
0261: * specified <code>long</code> number multiplied by a power of ten.
0262: *
0263: * @param m the <code>long</code> multiplier.
0264: * @param n the power of ten exponent.
0265: * @return <code>multiplier * 10<sup>n</sup></code>.
0266: * @JVM-1.1+@
0267: public static double toDoublePow10(long m, int n) {
0268: if (m == 0)
0269: return 0.0;
0270: if (m == Long.MIN_VALUE)
0271: return toDoublePow10(Long.MIN_VALUE / 10, n + 1);
0272: if (m < 0)
0273: return -toDoublePow10(-m, n);
0274: if (n >= 0) { // Positive power.
0275: if (n > 308)
0276: return Double.POSITIVE_INFINITY;
0277: // Works with 4 x 32 bits registers (x3:x2:x1:x0)
0278: long x0 = 0; // 32 bits.
0279: long x1 = 0; // 32 bits.
0280: long x2 = m & MASK_32; // 32 bits.
0281: long x3 = m >>> 32; // 32 bits.
0282: int pow2 = 0;
0283: while (n != 0) {
0284: int i = (n >= POW5_INT.length) ? POW5_INT.length - 1 : n;
0285: int coef = POW5_INT[i]; // 31 bits max.
0286:
0287: if (((int)x0) != 0) x0 *= coef; // 63 bits max.
0288: if (((int)x1) != 0) x1 *= coef; // 63 bits max.
0289: x2 *= coef; // 63 bits max.
0290: x3 *= coef; // 63 bits max.
0291:
0292: x1 += x0 >>> 32;
0293: x0 &= MASK_32;
0294:
0295: x2 += x1 >>> 32;
0296: x1 &= MASK_32;
0297:
0298: x3 += x2 >>> 32;
0299: x2 &= MASK_32;
0300:
0301: // Adjusts powers.
0302: pow2 += i;
0303: n -= i;
0304:
0305: // Normalizes (x3 should be 32 bits max).
0306: long carry = x3 >>> 32;
0307: if (carry != 0) { // Shift.
0308: x0 = x1;
0309: x1 = x2;
0310: x2 = x3 & MASK_32;
0311: x3 = carry;
0312: pow2 += 32;
0313: }
0314: }
0315:
0316: // Merges registers to a 63 bits mantissa.
0317: int shift = 31 - MathLib.bitLength(x3); // -1..30
0318: pow2 -= shift;
0319: long mantissa = (shift < 0) ?
0320: (x3 << 31) | (x2 >>> 1) : // x3 is 32 bits.
0321: (((x3 << 32) | x2) << shift) | (x1 >>> (32 - shift));
0322: return toDoublePow2(mantissa, pow2);
0323:
0324: } else { // n < 0
0325: if (n < -324 - 20) // m less than 20 digits.
0326: return 0.0;
0327:
0328: // Works with x1:x0 126 bits register.
0329: long x1 = m; // 63 bits.
0330: long x0 = 0; // 63 bits.
0331: int pow2 = 0;
0332: while (true) {
0333:
0334: // Normalizes x1:x0
0335: int shift = 63 - MathLib.bitLength(x1);
0336: x1 <<= shift;
0337: x1 |= x0 >>> (63 - shift);
0338: x0 = (x0 << shift) & MASK_63;
0339: pow2 -= shift;
0340:
0341: // Checks if division has to be performed.
0342: if (n == 0)
0343: break; // Done.
0344:
0345: // Retrieves power of 5 divisor.
0346: int i = (-n >= POW5_INT.length) ? POW5_INT.length - 1 : -n;
0347: int divisor = POW5_INT[i];
0348:
0349: // Performs the division (126 bits by 31 bits).
0350: long wh = (x1 >>> 32);
0351: long qh = wh / divisor;
0352: long r = wh - qh * divisor;
0353: long wl = (r << 32) | (x1 & MASK_32);
0354: long ql = wl / divisor;
0355: r = wl - ql * divisor;
0356: x1 = (qh << 32) | ql;
0357:
0358: wh = (r << 31) | (x0 >>> 32);
0359: qh = wh / divisor;
0360: r = wh - qh * divisor;
0361: wl = (r << 32) | (x0 & MASK_32);
0362: ql = wl / divisor;
0363: x0 = (qh << 32) | ql;
0364:
0365: // Adjusts powers.
0366: n += i;
0367: pow2 -= i;
0368: }
0369: return toDoublePow2(x1, pow2);
0370: }
0371: }
0372:
0373: private static final long MASK_63 = 0x7FFFFFFFFFFFFFFFL;
0374:
0375: private static final long MASK_32 = 0xFFFFFFFFL;
0376:
0377: private static final int[] POW5_INT = { 1, 5, 25, 125, 625, 3125, 15625,
0378: 78125, 390625, 1953125, 9765625, 48828125, 244140625, 1220703125 };
0379: /**/
0380:
0381: /**
0382: * Returns the closest <code>long</code> representation of the
0383: * specified <code>double</code> number multiplied by a power of two.
0384: *
0385: * @param d the <code>double</code> multiplier.
0386: * @param n the power of two exponent.
0387: * @return <code>d * 2<sup>n</sup></code>
0388: * @throws ArithmeticException if the conversion cannot be performed
0389: * (NaN, Infinity or overflow).
0390: * @JVM-1.1+@
0391: public static long toLongPow2(double d, int n) {
0392: long bits = Double.doubleToLongBits(d);
0393: boolean isNegative = (bits >> 63) != 0;
0394: int exp = ((int)(bits >> 52)) & 0x7FF;
0395: long m = bits & 0x000fffffffffffffL;
0396: if (exp == 0x7FF) throw new ArithmeticException(
0397: "Cannot convert to long (Infinity or NaN)");
0398: if (exp == 0) {
0399: if (m == 0) return 0L;
0400: return toLongPow2(d * 18014398509481984L, n - 54); // 2^54 Exact.
0401: }
0402: m |= 0x0010000000000000L; // Sets MSB (bit 52)
0403: long shift = exp - 1023L - 52 + n; // Use long to avoid overflow.
0404: if (shift <= -64) return 0L;
0405: if (shift >= 11) throw new ArithmeticException(
0406: "Cannot convert to long (overflow)");
0407: m = (shift >= 0) ? m << shift :
0408: (m >> -shift) + ((m >> -(shift + 1)) & 1) ; // Rounding.
0409: return isNegative ? -m : m;
0410: }
0411: /**/
0412:
0413: /**
0414: * Returns the closest <code>long</code> representation of the
0415: * specified <code>double</code> number multiplied by a power of ten.
0416: *
0417: * @param d the <code>double</code> multiplier.
0418: * @param n the power of two exponent.
0419: * @return <code>d * 10<sup>n</sup></code>.
0420: @JVM-1.1+@
0421: public static long toLongPow10(double d, int n) {
0422: long bits = Double.doubleToLongBits(d);
0423: boolean isNegative = (bits >> 63) != 0;
0424: int exp = ((int)(bits >> 52)) & 0x7FF;
0425: long m = bits & 0x000fffffffffffffL;
0426: if (exp == 0x7FF) throw new ArithmeticException(
0427: "Cannot convert to long (Infinity or NaN)");
0428: if (exp == 0) {
0429: if (m == 0) return 0L;
0430: return toLongPow10(d * 1E16, n - 16);
0431: }
0432: m |= 0x0010000000000000L; // Sets MSB (bit 52)
0433: int pow2 = exp - 1023 - 52;
0434: // Retrieves 63 bits m with n == 0.
0435: if (n >= 0) {
0436: // Works with 4 x 32 bits registers (x3:x2:x1:x0)
0437: long x0 = 0; // 32 bits.
0438: long x1 = 0; // 32 bits.
0439: long x2 = m & MASK_32; // 32 bits.
0440: long x3 = m >>> 32; // 32 bits.
0441: while (n != 0) {
0442: int i = (n >= POW5_INT.length) ? POW5_INT.length - 1 : n;
0443: int coef = POW5_INT[i]; // 31 bits max.
0444:
0445: if (((int)x0) != 0) x0 *= coef; // 63 bits max.
0446: if (((int)x1) != 0) x1 *= coef; // 63 bits max.
0447: x2 *= coef; // 63 bits max.
0448: x3 *= coef; // 63 bits max.
0449:
0450: x1 += x0 >>> 32;
0451: x0 &= MASK_32;
0452:
0453: x2 += x1 >>> 32;
0454: x1 &= MASK_32;
0455:
0456: x3 += x2 >>> 32;
0457: x2 &= MASK_32;
0458:
0459: // Adjusts powers.
0460: pow2 += i;
0461: n -= i;
0462:
0463: // Normalizes (x3 should be 32 bits max).
0464: long carry = x3 >>> 32;
0465: if (carry != 0) { // Shift.
0466: x0 = x1;
0467: x1 = x2;
0468: x2 = x3 & MASK_32;
0469: x3 = carry;
0470: pow2 += 32;
0471: }
0472: }
0473:
0474: // Merges registers to a 63 bits mantissa.
0475: int shift = 31 - MathLib.bitLength(x3); // -1..30
0476: pow2 -= shift;
0477: m = (shift < 0) ?
0478: (x3 << 31) | (x2 >>> 1) : // x3 is 32 bits.
0479: (((x3 << 32) | x2) << shift) | (x1 >>> (32 - shift));
0480:
0481: } else { // n < 0
0482:
0483: // Works with x1:x0 126 bits register.
0484: long x1 = m; // 63 bits.
0485: long x0 = 0; // 63 bits.
0486: while (true) {
0487:
0488: // Normalizes x1:x0
0489: int shift = 63 - MathLib.bitLength(x1);
0490: x1 <<= shift;
0491: x1 |= x0 >>> (63 - shift);
0492: x0 = (x0 << shift) & MASK_63;
0493: pow2 -= shift;
0494:
0495: // Checks if division has to be performed.
0496: if (n == 0)
0497: break; // Done.
0498:
0499: // Retrieves power of 5 divisor.
0500: int i = (-n >= POW5_INT.length) ? POW5_INT.length - 1 : -n;
0501: int divisor = POW5_INT[i];
0502:
0503: // Performs the division (126 bits by 31 bits).
0504: long wh = (x1 >>> 32);
0505: long qh = wh / divisor;
0506: long r = wh - qh * divisor;
0507: long wl = (r << 32) | (x1 & MASK_32);
0508: long ql = wl / divisor;
0509: r = wl - ql * divisor;
0510: x1 = (qh << 32) | ql;
0511:
0512: wh = (r << 31) | (x0 >>> 32);
0513: qh = wh / divisor;
0514: r = wh - qh * divisor;
0515: wl = (r << 32) | (x0 & MASK_32);
0516: ql = wl / divisor;
0517: x0 = (qh << 32) | ql;
0518:
0519: // Adjusts powers.
0520: n += i;
0521: pow2 -= i;
0522: }
0523: m = x1;
0524: }
0525: if (pow2 > 0) throw new ArithmeticException("Overflow");
0526: if (pow2 < -63) return 0;
0527: m = (m >> -pow2) + ((m >> -(pow2 + 1)) & 1) ; // Rounding.
0528: return isNegative ? -m : m;
0529: }
0530: /**/
0531:
0532: /**
0533: * Returns the largest power of 2 that is less than or equal to the
0534: * the specified positive value.
0535: *
0536: * @param d the <code>double</code> number.
0537: * @return <code>floor(Log2(abs(d)))</code>
0538: * @throws ArithmeticException if <code>d <= 0<code> or <code>d</code>
0539: * is <code>NaN</code> or <code>Infinity</code>.
0540: * @JVM-1.1+@
0541: public static int floorLog2(double d) {
0542: if (d <= 0) throw new ArithmeticException("Negative number or zero");
0543: long bits = Double.doubleToLongBits(d);
0544: int exp = ((int)(bits >> 52)) & 0x7FF;
0545: if (exp == 0x7FF) throw new ArithmeticException("Infinity or NaN");
0546: if (exp == 0) // Degenerated.
0547: return floorLog2(d * 18014398509481984L) - 54; // 2^54 Exact.
0548: return exp - 1023;
0549: }
0550: /**/
0551:
0552: /**
0553: * Returns the largest power of 10 that is less than or equal to the
0554: * the specified positive value.
0555: *
0556: * @param d the <code>double</code> number.
0557: * @return <code>floor(Log10(abs(d)))</code>
0558: * @throws ArithmeticException if <code>d <= 0<code> or <code>d</code>
0559: * is <code>NaN</code> or <code>Infinity</code>.
0560: * @JVM-1.1+@
0561: public static int floorLog10(double d) {
0562: int guess = (int) (LOG2_DIV_LOG10 * MathLib.floorLog2(d));
0563: double pow10 = MathLib.toDoublePow10(1, guess);
0564: if ((pow10 <= d) && (pow10 * 10 > d)) // Good guess.
0565: return guess;
0566: if (pow10 > d) // Too large.
0567: return guess - 1;
0568: return guess + 1;
0569: }
0570: private static final double LOG2_DIV_LOG10 = 0.3010299956639811952137388947;
0571: /**/
0572:
0573: /**
0574: * The natural logarithm.
0575: * @JVM-1.1+@
0576: public static final double E = 2.71828182845904523536028747135266;
0577:
0578: /**
0579: * The ratio of the circumference of a circle to its diameter.
0580: * @JVM-1.1+@
0581: public static final double PI = 3.1415926535897932384626433832795;
0582:
0583: /**
0584: * Half the ratio of the circumference of a circle to its diameter.
0585: * @JVM-1.1+@
0586: public static final double HALF_PI = 1.5707963267948966192313216916398;
0587:
0588: /**
0589: * Twice the ratio of the circumference of a circle to its diameter.
0590: * @JVM-1.1+@
0591: public static final double TWO_PI = 6.283185307179586476925286766559;
0592:
0593: /**
0594: * Four time the ratio of the circumference of a circle to its diameter.
0595: * @JVM-1.1+@
0596: public static final double FOUR_PI = 12.566370614359172953850573533118;
0597:
0598: /**
0599: * Holds {@link #PI} * {@link #PI}.
0600: * @JVM-1.1+@
0601: public static final double PI_SQUARE = 9.8696044010893586188344909998762;
0602:
0603: /**
0604: * The natural logarithm of two.
0605: * @JVM-1.1+@
0606: public static final double LOG2 = 0.69314718055994530941723212145818;
0607:
0608: /**
0609: * The natural logarithm of ten.
0610: * @JVM-1.1+@
0611: public static final double LOG10 = 2.3025850929940456840179914546844;
0612:
0613: /**
0614: * The square root of two.
0615: * @JVM-1.1+@
0616: public static final double SQRT2 = 1.4142135623730950488016887242097;
0617:
0618: /**
0619: * Not-A-Number.
0620: * @JVM-1.1+@
0621: public static final double NaN = 0.0 / 0.0;
0622:
0623: /**
0624: * Infinity.
0625: * @JVM-1.1+@
0626: public static final double Infinity = 1.0 / 0.0;
0627:
0628: /**/
0629:
0630: /**
0631: * Converts an angle in degrees to radians.
0632: *
0633: * @param degrees the angle in degrees.
0634: * @return the specified angle in radians.
0635: * @JVM-1.1+@
0636: public static double toRadians(double degrees) {
0637: return degrees * (PI / 180.0);
0638: }
0639: /**/
0640:
0641: /**
0642: * Converts an angle in radians to degrees.
0643: *
0644: * @param radians the angle in radians.
0645: * @return the specified angle in degrees.
0646: * @JVM-1.1+@
0647: public static double toDegrees(double radians) {
0648: return radians * (180.0 / PI);
0649: }
0650: /**/
0651:
0652: /**
0653: * Returns the positive square root of the specified value.
0654: *
0655: * @param x the value.
0656: * @return <code>java.lang.Math.sqrt(x)</code>
0657: * @JVM-1.1+@
0658: public static double sqrt(double x) {
0659: return Math.sqrt(x); // CLDC 1.1
0660: }
0661: /**/
0662:
0663: /**
0664: * Returns the remainder of the division of the specified two arguments.
0665: *
0666: * @param x the dividend.
0667: * @param y the divisor.
0668: * @return <code>x - round(x / y) * y</code>
0669: * @JVM-1.1+@
0670: public static double rem(double x, double y) {
0671: double tmp = x / y;
0672: if (MathLib.abs(tmp) <= Long.MAX_VALUE) {
0673: return x - MathLib.round(tmp) * y;
0674: } else {
0675: return NaN;
0676: }
0677: }
0678: /**/
0679:
0680: /**
0681: * Returns the smallest (closest to negative infinity)
0682: * <code>double</code> value that is not less than the argument and is
0683: * equal to a mathematical integer.
0684: *
0685: * @param x the value.
0686: * @return <code>java.lang.Math.ceil(x)</code>
0687: * @JVM-1.1+@
0688: public static double ceil(double x) {
0689: return Math.ceil(x); // CLDC 1.1
0690: }
0691: /**/
0692:
0693: /**
0694: * Returns the largest (closest to positive infinity)
0695: * <code>double</code> value that is not greater than the argument and
0696: * is equal to a mathematical integer.
0697: *
0698: * @param x the value.
0699: * @return <code>java.lang.Math.ceil(x)</code>
0700: * @JVM-1.1+@
0701: public static double floor(double x) {
0702: return Math.floor(x); // CLDC 1.1
0703: }
0704: /**/
0705:
0706: /**
0707: * Returns the trigonometric sine of the specified angle in radians.
0708: *
0709: * @param radians the angle in radians.
0710: * @return <code>java.lang.Math.sin(radians)</code>
0711: * @JVM-1.1+@
0712: public static double sin(double radians) {
0713: return Math.sin(radians); // CLDC 1.1
0714: }
0715: /**/
0716:
0717: /**
0718: * Returns the trigonometric cosine of the specified angle in radians.
0719: *
0720: * @param radians the angle in radians.
0721: * @return <code>java.lang.Math.cos(radians)</code>
0722: * @JVM-1.1+@
0723: public static double cos(double radians) {
0724: return Math.cos(radians); // CLDC 1.1
0725: }
0726: /**/
0727:
0728: /**
0729: * Returns the trigonometric tangent of the specified angle in radians.
0730: *
0731: * @param radians the angle in radians.
0732: * @return <code>java.lang.Math.tan(radians)</code>
0733: * @JVM-1.1+@
0734: public static double tan(double radians) {
0735: return Math.tan(radians); // CLDC 1.1
0736: }
0737: /**/
0738:
0739: /**
0740: * Returns the arc sine of the specified value,
0741: * in the range of -<i>pi</i>/2 through <i>pi</i>/2.
0742: *
0743: * @param x the value whose arc sine is to be returned.
0744: * @return the arc sine in radians for the specified value.
0745: * @JVM-1.1+@
0746: public static double asin(double x) {
0747: if (x < -1.0 || x > 1.0) return MathLib.NaN;
0748: if (x == -1.0) return - HALF_PI;
0749: if (x == 1.0) return HALF_PI;
0750: return MathLib.atan(x / MathLib.sqrt(1.0 - x * x));
0751: }
0752: /**/
0753:
0754: /**
0755: * Returns the arc cosine of the specified value,
0756: * in the range of 0.0 through <i>pi</i>.
0757: *
0758: * @param x the value whose arc cosine is to be returned.
0759: * @return the arc cosine in radians for the specified value.
0760: * @JVM-1.1+@
0761: public static double acos(double x) {
0762: return HALF_PI - MathLib.asin(x);
0763: }
0764: /**/
0765:
0766: /**
0767: * Returns the arc tangent of the specified value,
0768: * in the range of -<i>pi</i>/2 through <i>pi</i>/2.
0769: *
0770: * @param x the value whose arc tangent is to be returned.
0771: * @return the arc tangent in radians for the specified value.
0772: * @see <a href="http://mathworld.wolfram.com/InverseTangent.html">
0773: * Inverse Tangent -- from MathWorld</a>
0774: * @JVM-1.1+@
0775: public static double atan(double x) {
0776: return MathLib._atan(x);
0777: }
0778: /**/
0779:
0780: /**
0781: * Returns the angle theta such that
0782: * <code>(x == cos(theta)) && (y == sin(theta))</code>.
0783: *
0784: * @param y the y value.
0785: * @param x the x value.
0786: * @return the angle theta in radians.
0787: * @JVM-1.1+@
0788: public static double atan2(double y, double x) {
0789: final double epsilon = 1E-128;
0790: if (MathLib.abs(x) > epsilon) {
0791: double temp = MathLib.atan(MathLib.abs(y) / MathLib.abs(x));
0792: if( x < 0.0 ) temp = PI - temp;
0793: if( y < 0.0 ) temp = TWO_PI - temp;
0794: return temp;
0795: } else if( y > epsilon) {
0796: return HALF_PI;
0797: } else if( y < -epsilon) {
0798: return 3 * HALF_PI;
0799: } else {
0800: return 0.0;
0801: }
0802: }
0803: /**/
0804:
0805: /**
0806: * Returns the hyperbolic sine of x.
0807: *
0808: * @param x the value for which the hyperbolic sine is calculated.
0809: * @return <code>(exp(x) - exp(-x)) / 2</code>
0810: * @JVM-1.1+@
0811: public static double sinh(double x) {
0812: return (MathLib.exp(x) - MathLib.exp(-x)) * 0.5;
0813: }
0814: /**/
0815:
0816: /**
0817: * Returns the hyperbolic cosine of x.
0818: *
0819: * @param x the value for which the hyperbolic cosine is calculated.
0820: * @return <code>(exp(x) + exp(-x)) / 2</code>
0821: * @JVM-1.1+@
0822: public static double cosh(double x) {
0823: return (MathLib.exp(x) + MathLib.exp(-x)) * 0.5;
0824: }
0825: /**/
0826:
0827: /**
0828: * Returns the hyperbolic tangent of x.
0829: *
0830: * @param x the value for which the hyperbolic tangent is calculated.
0831: * @return <code>(exp(2 * x) - 1) / (exp(2 * x) + 1)</code>
0832: * @JVM-1.1+@
0833: public static double tanh(double x) {
0834: return (MathLib.exp(2 * x) - 1) / (MathLib.exp(2 * x) + 1);
0835: }
0836: /**/
0837:
0838: /**
0839: * Returns <i>{@link #E e}</i> raised to the specified power.
0840: *
0841: * @param x the exponent.
0842: * @return <code><i>e</i><sup>x</sup></code>
0843: * @see <a href="http://mathworld.wolfram.com/ExponentialFunction.html">
0844: * Exponential Function -- from MathWorld</a>
0845: * @JVM-1.1+@
0846: public static double exp(double x) {
0847: return MathLib._ieee754_exp(x);
0848: }
0849: /**/
0850:
0851: /**
0852: * Returns the natural logarithm (base <i>{@link #E e}</i>) of the specified
0853: * value.
0854: *
0855: * @param x the value greater than <code>0.0</code>.
0856: * @return the value y such as <code><i>e</i><sup>y</sup> == x</code>
0857: * @JVM-1.1+@
0858: public static double log(double x) {
0859: return MathLib._ieee754_log(x);
0860: }
0861: /**/
0862:
0863: /**
0864: * Returns the decimal logarithm of the specified value.
0865: *
0866: * @param x the value greater than <code>0.0</code>.
0867: * @return the value y such as <code>10<sup>y</sup> == x</code>
0868: * @JVM-1.1+@
0869: public static double log10(double x) {
0870: return log(x) * INV_LOG10;
0871: }
0872: private static double INV_LOG10 = 0.43429448190325182765112891891661;
0873: /**/
0874:
0875: /**
0876: * Returns the value of the first argument raised to the power of the
0877: * second argument.
0878: *
0879: * @param x the base.
0880: * @param y the exponent.
0881: * @return <code>x<sup>y</sup></code>
0882: * @JVM-1.1+@
0883: public static double pow(double x, double y) {
0884: /**/
0885: /* @JVM-1.4+@ // Use java.lang.Math value.
0886: if (true) return Math.pow(x, y);
0887: /**/
0888: /* @JVM-1.1+@ // Else (J2ME) use close approximation (+/- LSB)
0889: if ((x < 0) && (y == (int)y)) return
0890: (((int)y) & 1) == 0 ? pow(-x, y) : -pow(-x, y);
0891: return MathLib.exp(y * MathLib.log(x));
0892: }
0893: /**/
0894:
0895: /**
0896: * Returns the closest <code>int</code> to the specified argument.
0897: *
0898: * @param f the <code>float</code> value to be rounded to a <code>int</code>
0899: * @return the nearest <code>int</code> value.
0900: * @JVM-1.1+@
0901: public static int round(float f) {
0902: return (int) floor(f + 0.5f);
0903: }
0904: /**/
0905:
0906: /**
0907: * Returns the closest <code>long</code> to the specified argument.
0908: *
0909: * @param d the <code>double</code> value to be rounded to a
0910: * <code>long</code>
0911: * @return the nearest <code>long</code> value.
0912: * @JVM-1.1+@
0913: public static long round(double d) {
0914: return (long) floor(d + 0.5d);
0915: }
0916: /**/
0917:
0918: /**
0919: * Returns a random number between zero and one.
0920: *
0921: * @return a <code>double</code> greater than or equal
0922: * to <code>0.0</code> and less than <code>1.0</code>.
0923: * @JVM-1.1+@
0924: public static double random() {
0925: return random(0, Integer.MAX_VALUE) * NORMALIZATION_FACTOR;
0926: }
0927: private static final double NORMALIZATION_FACTOR = -1.0 / Integer.MIN_VALUE;
0928: /**/
0929:
0930: /**
0931: * Returns the absolute value of the specified <code>int</code> argument.
0932: *
0933: * @param i the <code>int</code> value.
0934: * @return <code>i</code> or <code>-i</code>
0935: */
0936: public static int abs(int i) {
0937: return (i < 0) ? -i : i;
0938: }
0939:
0940: /**
0941: * Returns the absolute value of the specified <code>long</code> argument.
0942: *
0943: * @param l the <code>long</code> value.
0944: * @return <code>l</code> or <code>-l</code>
0945: */
0946: public static long abs(long l) {
0947: return (l < 0) ? -l : l;
0948: }
0949:
0950: /**
0951: * Returns the absolute value of the specified <code>float</code> argument.
0952: *
0953: * @param f the <code>float</code> value.
0954: * @return <code>f</code> or <code>-f</code>
0955: * @JVM-1.1+@
0956: public static float abs(float f) {
0957: return (f < 0) ? -f : f;
0958: }
0959: /**/
0960:
0961: /**
0962: * Returns the absolute value of the specified <code>double</code> argument.
0963: *
0964: * @param d the <code>double</code> value.
0965: * @return <code>d</code> or <code>-d</code>
0966: * @JVM-1.1+@
0967: public static double abs(double d) {
0968: return (d < 0) ? -d : d;
0969: }
0970: /**/
0971:
0972: /**
0973: * Returns the greater of two <code>int</code> values.
0974: *
0975: * @param x the first value.
0976: * @param y the second value.
0977: * @return the larger of <code>x</code> and <code>y</code>.
0978: */
0979: public static int max(int x, int y) {
0980: return (x >= y) ? x : y;
0981: }
0982:
0983: /**
0984: * Returns the greater of two <code>long</code> values.
0985: *
0986: * @param x the first value.
0987: * @param y the second value.
0988: * @return the larger of <code>x</code> and <code>y</code>.
0989: */
0990: public static long max(long x, long y) {
0991: return (x >= y) ? x : y;
0992: }
0993:
0994: /**
0995: * Returns the greater of two <code>float</code> values.
0996: *
0997: * @param x the first value.
0998: * @param y the second value.
0999: * @return the larger of <code>x</code> and <code>y</code>.
1000: * @JVM-1.1+@
1001: public static float max(float x, float y) {
1002: return (x >= y) ? x : y;
1003: }
1004: /**/
1005:
1006: /**
1007: * Returns the greater of two <code>double</code> values.
1008: *
1009: * @param x the first value.
1010: * @param y the second value.
1011: * @return the larger of <code>x</code> and <code>y</code>.
1012: * @JVM-1.1+@
1013: public static double max(double x, double y) {
1014: return (x >= y) ? x : y;
1015: }
1016: /**/
1017:
1018: /**
1019: * Returns the smaller of two <code>int</code> values.
1020: *
1021: * @param x the first value.
1022: * @param y the second value.
1023: * @return the smaller of <code>x</code> and <code>y</code>.
1024: */
1025: public static int min(int x, int y) {
1026: return (x < y) ? x : y;
1027: }
1028:
1029: /**
1030: * Returns the smaller of two <code>long</code> values.
1031: *
1032: * @param x the first value.
1033: * @param y the second value.
1034: * @return the smaller of <code>x</code> and <code>y</code>.
1035: */
1036: public static long min(long x, long y) {
1037: return (x < y) ? x : y;
1038: }
1039:
1040: /**
1041: * Returns the smaller of two <code>float</code> values.
1042: *
1043: * @param x the first value.
1044: * @param y the second value.
1045: * @return the smaller of <code>x</code> and <code>y</code>.
1046: * @JVM-1.1+@
1047: public static float min(float x, float y) {
1048: return (x < y) ? x : y;
1049: }
1050: /**/
1051:
1052: /**
1053: * Returns the smaller of two <code>double</code> values.
1054: *
1055: * @param x the first value.
1056: * @param y the second value.
1057: * @return the smaller of <code>x</code> and <code>y</code>.
1058: * @JVM-1.1+@
1059: public static double min(double x, double y) {
1060: return (x < y) ? x : y;
1061: }
1062: /**/
1063:
1064: ////////////////////////////////////////////////////////////////////////////
1065: /* @(#)s_atan.c 1.3 95/01/18 */
1066: /*
1067: * ====================================================
1068: * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
1069: *
1070: * Developed at SunSoft, a Sun Microsystems, Inc. business.
1071: * Permission to use, copy, modify, and distribute this
1072: * software is freely granted, provided that this notice
1073: * is preserved.
1074: * ====================================================
1075: *
1076: */
1077:
1078: /* atan(x)
1079: * Method
1080: * 1. Reduce x to positive by atan(x) = -atan(-x).
1081: * 2. According to the integer k=4t+0.25 chopped, t=x, the argument
1082: * is further reduced to one of the following intervals and the
1083: * arctangent of t is evaluated by the corresponding formula:
1084: *
1085: * [0,7/16] atan(x) = t-t^3*(a1+t^2*(a2+...(a10+t^2*a11)...)
1086: * [7/16,11/16] atan(x) = atan(1/2) + atan( (t-0.5)/(1+t/2) )
1087: * [11/16.19/16] atan(x) = atan( 1 ) + atan( (t-1)/(1+t) )
1088: * [19/16,39/16] atan(x) = atan(3/2) + atan( (t-1.5)/(1+1.5t) )
1089: * [39/16,INF] atan(x) = atan(INF) + atan( -1/t )
1090: *
1091: * Constants:
1092: * The hexadecimal values are the intended ones for the following
1093: * constants. The decimal values may be used, provided that the
1094: * compiler will convert from decimal to binary accurately enough
1095: * to produce the hexadecimal values shown.
1096: @JVM-1.1+@
1097: static final double atanhi[] = {
1098: 4.63647609000806093515e-01, // atan(0.5)hi 0x3FDDAC67, 0x0561BB4F
1099: 7.85398163397448278999e-01, // atan(1.0)hi 0x3FE921FB, 0x54442D18
1100: 9.82793723247329054082e-01, // atan(1.5)hi 0x3FEF730B, 0xD281F69B
1101: 1.57079632679489655800e+00, // atan(inf)hi 0x3FF921FB, 0x54442D18
1102: };
1103: static final double atanlo[] = {
1104: 2.26987774529616870924e-17, // atan(0.5)lo 0x3C7A2B7F, 0x222F65E2
1105: 3.06161699786838301793e-17, // atan(1.0)lo 0x3C81A626, 0x33145C07
1106: 1.39033110312309984516e-17, // atan(1.5)lo 0x3C700788, 0x7AF0CBBD
1107: 6.12323399573676603587e-17, // atan(inf)lo 0x3C91A626, 0x33145C07
1108: };
1109: static final double aT[] = {
1110: 3.33333333333329318027e-01, // 0x3FD55555, 0x5555550D
1111: -1.99999999998764832476e-01, // 0xBFC99999, 0x9998EBC4
1112: 1.42857142725034663711e-01, // 0x3FC24924, 0x920083FF
1113: -1.11111104054623557880e-01, // 0xBFBC71C6, 0xFE231671
1114: 9.09088713343650656196e-02, // 0x3FB745CD, 0xC54C206E
1115: -7.69187620504482999495e-02, // 0xBFB3B0F2, 0xAF749A6D
1116: 6.66107313738753120669e-02, // 0x3FB10D66, 0xA0D03D51
1117: -5.83357013379057348645e-02, // 0xBFADDE2D, 0x52DEFD9A
1118: 4.97687799461593236017e-02, // 0x3FA97B4B, 0x24760DEB
1119: -3.65315727442169155270e-02, // 0xBFA2B444, 0x2C6A6C2F
1120: 1.62858201153657823623e-02, // 0x3F90AD3A, 0xE322DA11
1121: };
1122: static final double
1123: one = 1.0,
1124: huge = 1.0e300;
1125: static double _atan(double x)
1126: {
1127: double w,s1,s2,z;
1128: int ix,hx,id;
1129: long xBits = Double.doubleToLongBits(x);
1130: int __HIx = (int) (xBits >> 32);
1131: int __LOx = (int) xBits;
1132:
1133: hx = __HIx;
1134: ix = hx&0x7fffffff;
1135: if(ix>=0x44100000) { // if |x| >= 2^66
1136: if(ix>0x7ff00000||
1137: (ix==0x7ff00000&&(__LOx!=0)))
1138: return x+x; // NaN
1139: if(hx>0) return atanhi[3]+atanlo[3];
1140: else return -atanhi[3]-atanlo[3];
1141: } if (ix < 0x3fdc0000) { // |x| < 0.4375
1142: if (ix < 0x3e200000) { // |x| < 2^-29
1143: if(huge+x>one) return x; // raise inexact
1144: }
1145: id = -1;
1146: } else {
1147: x = MathLib.abs(x);
1148: if (ix < 0x3ff30000) { // |x| < 1.1875
1149: if (ix < 0x3fe60000) { // 7/16 <=|x|<11/16
1150: id = 0; x = (2.0*x-one)/(2.0+x);
1151: } else { // 11/16<=|x|< 19/16
1152: id = 1; x = (x-one)/(x+one);
1153: }
1154: } else {
1155: if (ix < 0x40038000) { // |x| < 2.4375
1156: id = 2; x = (x-1.5)/(one+1.5*x);
1157: } else { // 2.4375 <= |x| < 2^66
1158: id = 3; x = -1.0/x;
1159: }
1160: }}
1161: // end of argument reduction
1162: z = x*x;
1163: w = z*z;
1164: // break sum from i=0 to 10 aT[i]z**(i+1) into odd and even poly
1165: s1 = z*(aT[0]+w*(aT[2]+w*(aT[4]+w*(aT[6]+w*(aT[8]+w*aT[10])))));
1166: s2 = w*(aT[1]+w*(aT[3]+w*(aT[5]+w*(aT[7]+w*aT[9]))));
1167: if (id<0) return x - x*(s1+s2);
1168: else {
1169: z = atanhi[id] - ((x*(s1+s2) - atanlo[id]) - x);
1170: return (hx<0)? -z:z;
1171: }
1172: }
1173: /**/
1174:
1175: ////////////////////////////////////////////////////////////////////////////
1176: /* @(#)e_log.c 1.3 95/01/18 */
1177: /*
1178: * ====================================================
1179: * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
1180: *
1181: * Developed at SunSoft, a Sun Microsystems, Inc. business.
1182: * Permission to use, copy, modify, and distribute this
1183: * software is freely granted, provided that this notice
1184: * is preserved.
1185: * ====================================================
1186: */
1187:
1188: /* __ieee754_log(x)
1189: * Return the logrithm of x
1190: *
1191: * Method :
1192: * 1. Argument Reduction: find k and f such that
1193: * x = 2^k * (1+f),
1194: * where sqrt(2)/2 < 1+f < sqrt(2) .
1195: *
1196: * 2. Approximation of log(1+f).
1197: * Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
1198: * = 2s + 2/3 s**3 + 2/5 s**5 + .....,
1199: * = 2s + s*R
1200: * We use a special Reme algorithm on [0,0.1716] to generate
1201: * a polynomial of degree 14 to approximate R The maximum error
1202: * of this polynomial approximation is bounded by 2**-58.45. In
1203: * other words,
1204: * 2 4 6 8 10 12 14
1205: * R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s +Lg6*s +Lg7*s
1206: * (the values of Lg1 to Lg7 are listed in the program)
1207: * and
1208: * | 2 14 | -58.45
1209: * | Lg1*s +...+Lg7*s - R(z) | <= 2
1210: * | |
1211: * Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2.
1212: * In order to guarantee error in log below 1ulp, we compute log
1213: * by
1214: * log(1+f) = f - s*(f - R) (if f is not too large)
1215: * log(1+f) = f - (hfsq - s*(hfsq+R)). (better accuracy)
1216: *
1217: * 3. Finally, log(x) = k*ln2 + log(1+f).
1218: * = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo)))
1219: * Here ln2 is split into two floating point number:
1220: * ln2_hi + ln2_lo,
1221: * where n*ln2_hi is always exact for |n| < 2000.
1222: *
1223: * Special cases:
1224: * log(x) is NaN with signal if x < 0 (including -INF) ;
1225: * log(+INF) is +INF; log(0) is -INF with signal;
1226: * log(NaN) is that NaN with no signal.
1227: *
1228: * Accuracy:
1229: * according to an error analysis, the error is always less than
1230: * 1 ulp (unit in the last place).
1231: *
1232: * Constants:
1233: * The hexadecimal values are the intended ones for the following
1234: * constants. The decimal values may be used, provided that the
1235: * compiler will convert from decimal to binary accurately enough
1236: * to produce the hexadecimal values shown.
1237: @JVM-1.1+@
1238: static final double
1239: ln2_hi = 6.93147180369123816490e-01, // 3fe62e42 fee00000
1240: ln2_lo = 1.90821492927058770002e-10, // 3dea39ef 35793c76
1241: two54 = 1.80143985094819840000e+16, // 43500000 00000000
1242: Lg1 = 6.666666666666735130e-01, // 3FE55555 55555593
1243: Lg2 = 3.999999999940941908e-01, // 3FD99999 9997FA04
1244: Lg3 = 2.857142874366239149e-01, // 3FD24924 94229359
1245: Lg4 = 2.222219843214978396e-01, // 3FCC71C5 1D8E78AF
1246: Lg5 = 1.818357216161805012e-01, // 3FC74664 96CB03DE
1247: Lg6 = 1.531383769920937332e-01, // 3FC39A09 D078C69F
1248: Lg7 = 1.479819860511658591e-01; // 3FC2F112 DF3E5244
1249: static final double zero = 0.0;
1250: static double _ieee754_log(double x)
1251: {
1252: double hfsq,f,s,z,R,w,t1,t2,dk;
1253: int k,hx,i,j;
1254: int lx; // unsigned
1255:
1256: long xBits = Double.doubleToLongBits(x);
1257: hx = (int) (xBits >> 32);
1258: lx = (int) xBits;
1259:
1260: k=0;
1261: if (hx < 0x00100000) { // x < 2**-1022
1262: if (((hx&0x7fffffff)|lx)==0)
1263: return -two54/zero; // log(+-0)=-inf
1264: if (hx<0) return (x-x)/zero; // log(-#) = NaN
1265: k -= 54; x *= two54; // subnormal number, scale up x
1266: xBits = Double.doubleToLongBits(x);
1267: hx = (int) (xBits >> 32); // high word of x
1268: }
1269: if (hx >= 0x7ff00000) return x+x;
1270: k += (hx>>20)-1023;
1271: hx &= 0x000fffff;
1272: i = (hx+0x95f64)&0x100000;
1273: xBits = Double.doubleToLongBits(x);
1274: int HIx = hx|(i^0x3ff00000); // normalize x or x/2
1275: xBits = ((HIx & 0xFFFFFFFFL) << 32) | (xBits & 0xFFFFFFFFL);
1276: x = Double.longBitsToDouble(xBits);
1277: k += (i>>20);
1278: f = x-1.0;
1279: if((0x000fffff&(2+hx))<3) { // |f| < 2**-20
1280: if(f==zero) if(k==0) return zero; else {dk=(double)k;
1281: return dk*ln2_hi+dk*ln2_lo;}
1282: R = f*f*(0.5-0.33333333333333333*f);
1283: if(k==0) return f-R; else {dk=(double)k;
1284: return dk*ln2_hi-((R-dk*ln2_lo)-f);}
1285: }
1286: s = f/(2.0+f);
1287: dk = (double)k;
1288: z = s*s;
1289: i = hx-0x6147a;
1290: w = z*z;
1291: j = 0x6b851-hx;
1292: t1= w*(Lg2+w*(Lg4+w*Lg6));
1293: t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7)));
1294: i |= j;
1295: R = t2+t1;
1296: if(i>0) {
1297: hfsq=0.5*f*f;
1298: if(k==0) return f-(hfsq-s*(hfsq+R)); else
1299: return dk*ln2_hi-((hfsq-(s*(hfsq+R)+dk*ln2_lo))-f);
1300: } else {
1301: if(k==0) return f-s*(f-R); else
1302: return dk*ln2_hi-((s*(f-R)-dk*ln2_lo)-f);
1303: }
1304: }
1305: /**/
1306:
1307: ////////////////////////////////////////////////////////////////////////////
1308: /* @(#)e_exp.c 1.6 04/04/22 */
1309: /*
1310: * ====================================================
1311: * Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved.
1312: *
1313: * Permission to use, copy, modify, and distribute this
1314: * software is freely granted, provided that this notice
1315: * is preserved.
1316: * ====================================================
1317: */
1318:
1319: /* __ieee754_exp(x)
1320: * Returns the exponential of x.
1321: *
1322: * Method
1323: * 1. Argument reduction:
1324: * Reduce x to an r so that |r| <= 0.5*ln2 ~ 0.34658.
1325: * Given x, find r and integer k such that
1326: *
1327: * x = k*ln2 + r, |r| <= 0.5*ln2.
1328: *
1329: * Here r will be represented as r = hi-lo for better
1330: * accuracy.
1331: *
1332: * 2. Approximation of exp(r) by a special rational function on
1333: * the interval [0,0.34658]:
1334: * Write
1335: * R(r**2) = r*(exp(r)+1)/(exp(r)-1) = 2 + r*r/6 - r**4/360 + ...
1336: * We use a special Remes algorithm on [0,0.34658] to generate
1337: * a polynomial of degree 5 to approximate R. The maximum error
1338: * of this polynomial approximation is bounded by 2**-59. In
1339: * other words,
1340: * R(z) ~ 2.0 + P1*z + P2*z**2 + P3*z**3 + P4*z**4 + P5*z**5
1341: * (where z=r*r, and the values of P1 to P5 are listed below)
1342: * and
1343: * | 5 | -59
1344: * | 2.0+P1*z+...+P5*z - R(z) | <= 2
1345: * | |
1346: * The computation of exp(r) thus becomes
1347: * 2*r
1348: * exp(r) = 1 + -------
1349: * R - r
1350: * r*R1(r)
1351: * = 1 + r + ----------- (for better accuracy)
1352: * 2 - R1(r)
1353: * where
1354: * 2 4 10
1355: * R1(r) = r - (P1*r + P2*r + ... + P5*r ).
1356: *
1357: * 3. Scale back to obtain exp(x):
1358: * From step 1, we have
1359: * exp(x) = 2^k * exp(r)
1360: *
1361: * Special cases:
1362: * exp(INF) is INF, exp(NaN) is NaN;
1363: * exp(-INF) is 0, and
1364: * for finite argument, only exp(0)=1 is exact.
1365: *
1366: * Accuracy:
1367: * according to an error analysis, the error is always less than
1368: * 1 ulp (unit in the last place).
1369: *
1370: * Misc. info.
1371: * For IEEE double
1372: * if x > 7.09782712893383973096e+02 then exp(x) overflow
1373: * if x < -7.45133219101941108420e+02 then exp(x) underflow
1374: *
1375: * Constants:
1376: * The hexadecimal values are the intended ones for the following
1377: * constants. The decimal values may be used, provided that the
1378: * compiler will convert from decimal to binary accurately enough
1379: * to produce the hexadecimal values shown.
1380: @JVM-1.1+@
1381: static final double
1382: halF[] = {0.5,-0.5,},
1383: twom1000= 9.33263618503218878990e-302, // 2**-1000=0x01700000,0
1384: o_threshold= 7.09782712893383973096e+02, // 0x40862E42, 0xFEFA39EF
1385: u_threshold= -7.45133219101941108420e+02, // 0xc0874910, 0xD52D3051
1386: ln2HI[] ={ 6.93147180369123816490e-01, // 0x3fe62e42, 0xfee00000
1387: -6.93147180369123816490e-01,},// 0xbfe62e42, 0xfee00000
1388: ln2LO[] ={ 1.90821492927058770002e-10, // 0x3dea39ef, 0x35793c76
1389: -1.90821492927058770002e-10,},// 0xbdea39ef, 0x35793c76
1390: invln2 = 1.44269504088896338700e+00, // 0x3ff71547, 0x652b82fe
1391: P1 = 1.66666666666666019037e-01, // 0x3FC55555, 0x5555553E
1392: P2 = -2.77777777770155933842e-03, // 0xBF66C16C, 0x16BEBD93
1393: P3 = 6.61375632143793436117e-05, // 0x3F11566A, 0xAF25DE2C
1394: P4 = -1.65339022054652515390e-06, // 0xBEBBBD41, 0xC5D26BF1
1395: P5 = 4.13813679705723846039e-08; // 0x3E663769, 0x72BEA4D0
1396: static double _ieee754_exp(double x) // default IEEE double exp
1397: {
1398: double y,hi = 0,lo = 0,c,t;
1399: int k = 0,xsb;
1400: int hx; // Unsigned.
1401: long xBits = Double.doubleToLongBits(x);
1402: int __HIx = (int) (xBits >> 32);
1403: int __LOx = (int) xBits;
1404:
1405: hx = __HIx; // high word of x
1406: xsb = (hx>>31)&1; // sign bit of x
1407: hx &= 0x7fffffff; // high word of |x|
1408:
1409: // filter out non-finite argument
1410: if(hx >= 0x40862E42) { // if |x|>=709.78...
1411: if(hx>=0x7ff00000) {
1412: if(((hx&0xfffff)|__LOx)!=0)
1413: return x+x; // NaN
1414: else return (xsb==0)? x:0.0; // exp(+-inf)={inf,0}
1415: }
1416: if(x > o_threshold) return huge*huge; // overflow
1417: if(x < u_threshold) return twom1000*twom1000; // underflow
1418: }
1419:
1420: // argument reduction
1421: if(hx > 0x3fd62e42) { // if |x| > 0.5 ln2
1422: if(hx < 0x3FF0A2B2) { // and |x| < 1.5 ln2
1423: hi = x-ln2HI[xsb]; lo=ln2LO[xsb]; k = 1-xsb-xsb;
1424: } else {
1425: k = (int)(invln2*x+halF[xsb]);
1426: t = k;
1427: hi = x - t*ln2HI[0]; // t*ln2HI is exact here
1428: lo = t*ln2LO[0];
1429: }
1430: x = hi - lo;
1431: }
1432: else if(hx < 0x3e300000) { // when |x|<2**-28
1433: if(huge+x>one) return one+x;// trigger inexact
1434: }
1435: else k = 0;
1436:
1437: // x is now in primary range
1438: t = x*x;
1439: c = x - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
1440: if(k==0) return one-((x*c)/(c-2.0)-x);
1441: else y = one-((lo-(x*c)/(2.0-c))-hi);
1442: long yBits = Double.doubleToLongBits(y);
1443: int __HIy = (int) (yBits >> 32);
1444: if(k >= -1021) {
1445: __HIy += (k<<20); // add k to y's exponent
1446: yBits = ((__HIy & 0xFFFFFFFFL) << 32) | (yBits & 0xFFFFFFFFL);
1447: y = Double.longBitsToDouble(yBits);
1448: return y;
1449: } else {
1450: __HIy += ((k+1000)<<20);// add k to y's exponent
1451: yBits = ((__HIy & 0xFFFFFFFFL) << 32) | (yBits & 0xFFFFFFFFL);
1452: y = Double.longBitsToDouble(yBits);
1453: return y*twom1000;
1454: }
1455: }
1456: /**/
1457:
1458: }
|