/*
* Javolution - Java(TM) Solution for Real-Time and Embedded Systems
* Copyright (C) 2005 - 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 javolution37.javolution.lang;
/**
* This utility class ensures cross-platform portability of the math
* library. Functions not supported by the platform are emulated.
* Developers may replace the current implementation with native
* implementations for faster execution.
*
* @author Jean-Marie Dautelle
* @version 1.1, November 6, 2004
*/
public final class MathLib {
/**
* Default constructor.
*/
private MathLib() {}
/**
* Returns a pseudo random int value in the range
* [min, max].
* @param min the minimum value inclusive.
* @param max the maximum value inclusive.
* @return a pseudo random number in the range [min, max].
*/
public synchronized static int randomInt(int min, int max) {
// This is a linear congruential pseudorandom number generator, as
// defined by D. H. Lehmer and described by Donald E. Knuth in The
// Art of Computer Programming, Volume 2: Seminumerical
// Algorithms, section 3.2.1.
MathLib.Seed = (MathLib.Seed * 0x5DEECE66DL + 0xBL) & ((1L << 48) - 1);
int i = (int)(MathLib.Seed >>> (48 - 32)); // 48 - bits
return (i < 0) ? min - i % (max - min + 1) : min + i % (max - min + 1);
}
private static long Seed = System.currentTimeMillis() + 8682522807148012L;
/**
* 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.
*/
public static double atan2(double y, double x) {
final double epsilon = 1E-128;
if (MathLib.abs(x) > epsilon) {
double temp = MathLib.atan(MathLib.abs(y) / MathLib.abs(x));
if( x < 0.0 ) temp = PI - temp;
if( y < 0.0 ) temp = TWO_PI - temp;
return temp;
} else if( y > epsilon) {
return HALF_PI;
} else if( y < -epsilon) {
return 3 * HALF_PI;
} else {
return 0.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) {
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 a random number between zero and one.
*
* @return a double greater than or equal
* to 0.0 and less than 1.0.
*/
public static double random() {
return randomInt(0, Integer.MAX_VALUE) * INV_MAX_VALUE;
}
private static final double INV_MAX_VALUE = 1.0 / Integer.MAX_VALUE;
/**/
/**
* 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;
}
/**/
/**
* Returns the best double approximation of the specified
* long mantissa multiplied by a power of 10
* (scientific notation).
*
* @param mantissa the mantissa as long.
* @param n the power of 10 exponent.
* @return mantissaÃ10n
**/
public static double toDouble(long mantissa, int n) {
return (n < 0) ? ((n >= -323) ? mantissa * POW_10_NEG[-n] : mantissa * 0.0) :
((n <= 308) ? mantissa * POW_10_POS[n] : mantissa * Double.POSITIVE_INFINITY);
}
// Note: Approximation for exponents > 21. This may introduce round-off
// errors (e.g. 1E23 represented as "9.999999999999999E22").
private static final double[] POW_10_POS = new double[] { 1E000, 1E001,
1E002, 1E003, 1E004, 1E005, 1E006, 1E007, 1E008, 1E009, 1E010,
1E011, 1E012, 1E013, 1E014, 1E015, 1E016, 1E017, 1E018, 1E019,
1E020, 1E021, 1E022, 1E023, 1E024, 1E025, 1E026, 1E027, 1E028,
1E029, 1E030, 1E031, 1E032, 1E033, 1E034, 1E035, 1E036, 1E037,
1E038, 1E039, 1E040, 1E041, 1E042, 1E043, 1E044, 1E045, 1E046,
1E047, 1E048, 1E049, 1E050, 1E051, 1E052, 1E053, 1E054, 1E055,
1E056, 1E057, 1E058, 1E059, 1E060, 1E061, 1E062, 1E063, 1E064,
1E065, 1E066, 1E067, 1E068, 1E069, 1E070, 1E071, 1E072, 1E073,
1E074, 1E075, 1E076, 1E077, 1E078, 1E079, 1E080, 1E081, 1E082,
1E083, 1E084, 1E085, 1E086, 1E087, 1E088, 1E089, 1E090, 1E091,
1E092, 1E093, 1E094, 1E095, 1E096, 1E097, 1E098, 1E099,
1E100, 1E101, 1E102, 1E103, 1E104, 1E105, 1E106, 1E107, 1E108,
1E109, 1E110, 1E111, 1E112, 1E113, 1E114, 1E115, 1E116, 1E117,
1E118, 1E119, 1E120, 1E121, 1E122, 1E123, 1E124, 1E125, 1E126,
1E127, 1E128, 1E129, 1E130, 1E131, 1E132, 1E133, 1E134, 1E135,
1E136, 1E137, 1E138, 1E139, 1E140, 1E141, 1E142, 1E143, 1E144,
1E145, 1E146, 1E147, 1E148, 1E149, 1E150, 1E151, 1E152, 1E153,
1E154, 1E155, 1E156, 1E157, 1E158, 1E159, 1E160, 1E161, 1E162,
1E163, 1E164, 1E165, 1E166, 1E167, 1E168, 1E169, 1E170, 1E171,
1E172, 1E173, 1E174, 1E175, 1E176, 1E177, 1E178, 1E179, 1E180,
1E181, 1E182, 1E183, 1E184, 1E185, 1E186, 1E187, 1E188, 1E189,
1E190, 1E191, 1E192, 1E193, 1E194, 1E195, 1E196, 1E197, 1E198,
1E199,
1E200, 1E201, 1E202, 1E203, 1E204, 1E205, 1E206, 1E207, 1E208,
1E209, 1E210, 1E211, 1E212, 1E213, 1E214, 1E215, 1E216, 1E217,
1E218, 1E219, 1E220, 1E221, 1E222, 1E223, 1E224, 1E225, 1E226,
1E227, 1E228, 1E229, 1E230, 1E231, 1E232, 1E233, 1E234, 1E235,
1E236, 1E237, 1E238, 1E239, 1E240, 1E241, 1E242, 1E243, 1E244,
1E245, 1E246, 1E247, 1E248, 1E249, 1E250, 1E251, 1E252, 1E253,
1E254, 1E255, 1E256, 1E257, 1E258, 1E259, 1E260, 1E261, 1E262,
1E263, 1E264, 1E265, 1E266, 1E267, 1E268, 1E269, 1E270, 1E271,
1E272, 1E273, 1E274, 1E275, 1E276, 1E277, 1E278, 1E279, 1E280,
1E281, 1E282, 1E283, 1E284, 1E285, 1E286, 1E287, 1E288, 1E289,
1E290, 1E291, 1E292, 1E293, 1E294, 1E295, 1E296, 1E297, 1E298,
1E299,
1E300, 1E301, 1E302, 1E303, 1E304, 1E305, 1E306, 1E307, 1E308 };
private static final double[] POW_10_NEG = new double[] { 1E-000, 1E-001,
1E-002, 1E-003, 1E-004, 1E-005, 1E-006, 1E-007, 1E-008, 1E-009, 1E-010,
1E-011, 1E-012, 1E-013, 1E-014, 1E-015, 1E-016, 1E-017, 1E-018, 1E-019,
1E-020, 1E-021, 1E-022, 1E-023, 1E-024, 1E-025, 1E-026, 1E-027, 1E-028,
1E-029, 1E-030, 1E-031, 1E-032, 1E-033, 1E-034, 1E-035, 1E-036, 1E-037,
1E-038, 1E-039, 1E-040, 1E-041, 1E-042, 1E-043, 1E-044, 1E-045, 1E-046,
1E-047, 1E-048, 1E-049, 1E-050, 1E-051, 1E-052, 1E-053, 1E-054, 1E-055,
1E-056, 1E-057, 1E-058, 1E-059, 1E-060, 1E-061, 1E-062, 1E-063, 1E-064,
1E-065, 1E-066, 1E-067, 1E-068, 1E-069, 1E-070, 1E-071, 1E-072, 1E-073,
1E-074, 1E-075, 1E-076, 1E-077, 1E-078, 1E-079, 1E-080, 1E-081, 1E-082,
1E-083, 1E-084, 1E-085, 1E-086, 1E-087, 1E-088, 1E-089, 1E-090, 1E-091,
1E-092, 1E-093, 1E-094, 1E-095, 1E-096, 1E-097, 1E-098, 1E-099,
1E-100, 1E-101, 1E-102, 1E-103, 1E-104, 1E-105, 1E-106, 1E-107, 1E-108,
1E-109, 1E-110, 1E-111, 1E-112, 1E-113, 1E-114, 1E-115, 1E-116, 1E-117,
1E-118, 1E-119, 1E-120, 1E-121, 1E-122, 1E-123, 1E-124, 1E-125, 1E-126,
1E-127, 1E-128, 1E-129, 1E-130, 1E-131, 1E-132, 1E-133, 1E-134, 1E-135,
1E-136, 1E-137, 1E-138, 1E-139, 1E-140, 1E-141, 1E-142, 1E-143, 1E-144,
1E-145, 1E-146, 1E-147, 1E-148, 1E-149, 1E-150, 1E-151, 1E-152, 1E-153,
1E-154, 1E-155, 1E-156, 1E-157, 1E-158, 1E-159, 1E-160, 1E-161, 1E-162,
1E-163, 1E-164, 1E-165, 1E-166, 1E-167, 1E-168, 1E-169, 1E-170, 1E-171,
1E-172, 1E-173, 1E-174, 1E-175, 1E-176, 1E-177, 1E-178, 1E-179, 1E-180,
1E-181, 1E-182, 1E-183, 1E-184, 1E-185, 1E-186, 1E-187, 1E-188, 1E-189,
1E-190, 1E-191, 1E-192, 1E-193, 1E-194, 1E-195, 1E-196, 1E-197, 1E-198,
1E-199,
1E-200, 1E-201, 1E-202, 1E-203, 1E-204, 1E-205, 1E-206, 1E-207, 1E-208,
1E-209, 1E-210, 1E-211, 1E-212, 1E-213, 1E-214, 1E-215, 1E-216, 1E-217,
1E-218, 1E-219, 1E-220, 1E-221, 1E-222, 1E-223, 1E-224, 1E-225, 1E-226,
1E-227, 1E-228, 1E-229, 1E-230, 1E-231, 1E-232, 1E-233, 1E-234, 1E-235,
1E-236, 1E-237, 1E-238, 1E-239, 1E-240, 1E-241, 1E-242, 1E-243, 1E-244,
1E-245, 1E-246, 1E-247, 1E-248, 1E-249, 1E-250, 1E-251, 1E-252, 1E-253,
1E-254, 1E-255, 1E-256, 1E-257, 1E-258, 1E-259, 1E-260, 1E-261, 1E-262,
1E-263, 1E-264, 1E-265, 1E-266, 1E-267, 1E-268, 1E-269, 1E-270, 1E-271,
1E-272, 1E-273, 1E-274, 1E-275, 1E-276, 1E-277, 1E-278, 1E-279, 1E-280,
1E-281, 1E-282, 1E-283, 1E-284, 1E-285, 1E-286, 1E-287, 1E-288, 1E-289,
1E-290, 1E-291, 1E-292, 1E-293, 1E-294, 1E-295, 1E-296, 1E-297, 1E-298,
1E-299,
1E-300, 1E-301, 1E-302, 1E-303, 1E-304, 1E-305, 1E-306, 1E-307, 1E-308,
1E-309, 1E-310, 1E-311, 1E-312, 1E-313, 1E-314, 1E-315, 1E-316, 1E-317,
1E-318, 1E-319, 1E-320, 1E-321, 1E-322, 1E-323 };
/**/
////////////////////////////////////////////////////////////////////////////
/* @(#)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; // raise inexact
}
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; // exp(+-inf)={inf,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;
}
}
/**/
}
Ads help maintain this website.