Source Code Cross Referenced for MathLib.java in  » Science » javolution-5.2 » javolution » lang » Java Source Code / Java DocumentationJava Source Code and Java Documentation

Java Source Code / Java Documentation
1. 6.0 JDK Core
2. 6.0 JDK Modules
3. 6.0 JDK Modules com.sun
4. 6.0 JDK Modules com.sun.java
5. 6.0 JDK Modules sun
6. 6.0 JDK Platform
7. Ajax
8. Apache Harmony Java SE
9. Aspect oriented
10. Authentication Authorization
11. Blogger System
12. Build
13. Byte Code
14. Cache
15. Chart
16. Chat
17. Code Analyzer
18. Collaboration
19. Content Management System
20. Database Client
21. Database DBMS
22. Database JDBC Connection Pool
23. Database ORM
24. Development
25. EJB Server geronimo
26. EJB Server GlassFish
27. EJB Server JBoss 4.2.1
28. EJB Server resin 3.1.5
29. ERP CRM Financial
30. ESB
31. Forum
32. GIS
33. Graphic Library
34. Groupware
35. HTML Parser
36. IDE
37. IDE Eclipse
38. IDE Netbeans
39. Installer
40. Internationalization Localization
41. Inversion of Control
42. Issue Tracking
43. J2EE
44. JBoss
45. JMS
46. JMX
47. Library
48. Mail Clients
49. Net
50. Parser
51. PDF
52. Portal
53. Profiler
54. Project Management
55. Report
56. RSS RDF
57. Rule Engine
58. Science
59. Scripting
60. Search Engine
61. Security
62. Sevlet Container
63. Source Control
64. Swing Library
65. Template Engine
66. Test Coverage
67. Testing
68. UML
69. Web Crawler
70. Web Framework
71. Web Mail
72. Web Server
73. Web Services
74. Web Services apache cxf 2.0.1
75. Web Services AXIS2
76. Wiki Engine
77. Workflow Engines
78. XML
79. XML UI
Java
Java Tutorial
Java Open Source
Jar File Download
Java Articles
Java Products
Java by API
Photoshop Tutorials
Maya Tutorials
Flash Tutorials
3ds-Max Tutorials
Illustrator Tutorials
GIMP Tutorials
C# / C Sharp
C# / CSharp Tutorial
C# / CSharp Open Source
ASP.Net
ASP.NET Tutorial
JavaScript DHTML
JavaScript Tutorial
JavaScript Reference
HTML / CSS
HTML CSS Reference
C / ANSI-C
C Tutorial
C++
C++ Tutorial
Ruby
PHP
Python
Python Tutorial
Python Open Source
SQL Server / T-SQL
SQL Server / T-SQL Tutorial
Oracle PL / SQL
Oracle PL/SQL Tutorial
PostgreSQL
SQL / MySQL
MySQL Tutorial
VB.Net
VB.Net Tutorial
Flash / Flex / ActionScript
VBA / Excel / Access / Word
XML
XML Tutorial
Microsoft Office PowerPoint 2007 Tutorial
Microsoft Office Excel 2007 Tutorial
Microsoft Office Word 2007 Tutorial
Java Source Code / Java Documentation » Science » javolution 5.2 » javolution.lang 
Source Cross Referenced  Class Diagram Java Document (Java Doc) 


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 &lt;= 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 &lt;= 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:        }
www.java2java.com | Contact Us
Copyright 2009 - 12 Demo Source and Support. All rights reserved.
All other trademarks are property of their respective owners.