Java > Open Source Codes > javolution > lang > MathLib


1 /*
2  * Javolution - Java(TM) Solution for Real-Time and Embedded Systems
3  * Copyright (C) 2006 - Javolution (http://javolution.org/)
4  * All rights reserved.
5  *
6  * Permission to use, copy, modify, and distribute this software is
7  * freely granted, provided that this notice is preserved.
8  */
9 package javolution.lang;
10
11 import java.util.Random ;
12
13 /**
14  * <p> This utility class ensures cross-platform portability of the math
15  * library. Functions not supported by the platform are emulated.
16  * Developers may replace the current implementation with native
17  * implementations for faster execution (unlike
18  * <code>java.lang.Math</code> this class can be customized).<p>
19  *
20  * @author <a HREF="mailto:jean-marie@dautelle.com">Jean-Marie Dautelle</a>
21  * @version 4.2, January 6, 2007
22  */
23 public final class MathLib {
24
25     /**
26      * Default constructor.
27      */
28     private MathLib() {
29     }
30
31     /**
32      * Returns a pseudo random <code>int</code> value in the range
33      * <code>[min, max]</code> (fast and thread-safe without synchronization).
34      *
35      * @param min the minimum value inclusive.
36      * @param max the maximum value exclusive.
37      * @return a pseudo random number in the range <code>[min, max]</code>.
38      */
39     public static int random(int min, int max) {
40         int next = RANDOM.nextInt();
41         if ((next >= min) && (next <= max))
42             return next;
43         next += Integer.MIN_VALUE;
44         if ((next >= min) && (next <= max))
45             return next;
46         // We should not have interval overflow as the interval has to be less
47 // or equal to Integer.MAX_VALUE (otherwise we would have exited before).
48 final int interval = 1 + max - min; // Positive.
49 if (interval <= 0)
50             throw new Error ("Interval error"); // In case.
51 return MathLib.abs(next % interval) + min;
52     }
53
54     private static final Random RANDOM = new Random ();
55
56     /**
57      * Returns a pseudo random <code>long</code> value in the range
58      * <code>[min, max]</code> (fast and thread-safe without synchronization).
59      *
60      * @param min the minimum value inclusive.
61      * @param max the maximum value inclusive.
62      * @return a pseudo random number in the range <code>[min, max]</code>.
63      */
64     public static long random(long min, long max) {
65         long next = RANDOM.nextLong();
66         if ((next >= min) && (next <= max))
67             return next;
68         next += Long.MIN_VALUE;
69         if ((next >= min) && (next <= max))
70             return next;
71         // We should not have interval overflow as the interval has to be less
72 // or equal to Long.MAX_VALUE (otherwise we would have exited before).
73 final long interval = 1L + max - min; // Positive.
74 if (interval <= 0)
75             throw new Error ("Interval error"); // In case.
76 return MathLib.abs(next % interval) + min;
77     }
78
79     /**
80      * Returns a pseudo random <code>float</code> value in the range
81      * <code>[min, max]</code> (fast and thread-safe without synchronization).
82      *
83      * @param min the minimum value inclusive.
84      * @param max the maximum value inclusive.
85      * @return a pseudo random number in the range <code>[min, max]</code>.
86      * @JVM-1.1+@
87      public static float random(float min, float max) {
88      return (float) random((double)min, (double) max);
89      }
90      /**/
91
92     /**
93      * Returns a pseudo random <code>double</code> value in the range
94      * <code>[min, max]</code> (fast and thread-safe without synchronization).
95      *
96      * @param min the minimum value inclusive.
97      * @param max the maximum value inclusive.
98      * @return a pseudo random number in the range <code>[min, max]</code>.
99      * @JVM-1.1+@
100      public static double random(double min, double max) {
101      double next = RANDOM.nextDouble(); // 0.0 .. 1.0
102      return min + (next * max) - (next * min);
103      // Due to potential numeric errors, both min and max are possible.
104      }
105      /**/
106
107     /**
108      * Returns the number of bits in the minimal two's-complement representation
109      * of the specified <code>int</code>, excluding a sign bit.
110      * For positive <code>int</code>, this is equivalent to the number of bits
111      * in the ordinary binary representation. For negative <code>int</code>,
112      * it is equivalent to the number of bits of the positive value
113      * <code>-(i + 1)</code>.
114      *
115      * @param i the <code>int</code> value for which the bit length is returned.
116      * @return the bit length of <code>i</code>.
117      */
118     public static int bitLength(int i) {
119         if (i < 0)
120             i = -++i;
121         return (i < 1 << 16) ? (i < 1 << 8) ? BIT_LENGTH[i]
122                 : BIT_LENGTH[i >>> 8] + 8
123                 : (i < 1 << 24) ? BIT_LENGTH[i >>> 16] + 16
124                         : BIT_LENGTH[i >>> 24] + 24;
125     }
126
127     private static final byte[] BIT_LENGTH = { 0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4,
128             4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6,
129             6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
130             6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
131             7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
132             7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
133             7, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
134             8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
135             8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
136             8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
137             8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
138             8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
139             8, 8, 8 };
140
141     /**
142      * Returns the number of bits in the minimal two's-complement representation
143      * of the specified <code>long</code>, excluding a sign bit.
144      * For positive <code>long</code>, this is equivalent to the number of bits
145      * in the ordinary binary representation. For negative <code>long</code>,
146      * it is equivalent to the number of bits of the positive value
147      * <code>-(l + 1)</code>.
148      *
149      * @param l the <code>long</code> value for which the bit length is returned.
150      * @return the bit length of <code>l</code>.
151      */
152     public static int bitLength(long l) {
153         int i = (int) (l >> 32);
154         if (i > 0)
155             return (i < 1 << 16) ? (i < 1 << 8) ? BIT_LENGTH[i] + 32
156                     : BIT_LENGTH[i >>> 8] + 40
157                     : (i < 1 << 24) ? BIT_LENGTH[i >>> 16] + 48
158                             : BIT_LENGTH[i >>> 24] + 56;
159         if (i < 0)
160             return bitLength(-++l);
161         i = (int) l;
162         return (i < 0) ? 32 : (i < 1 << 16) ? (i < 1 << 8) ? BIT_LENGTH[i]
163                 : BIT_LENGTH[i >>> 8] + 8
164                 : (i < 1 << 24) ? BIT_LENGTH[i >>> 16] + 16
165                         : BIT_LENGTH[i >>> 24] + 24;
166     }
167
168     /**
169      * Returns the number of digits in the minimal ten's-complement
170      * representation of the specified <code>int</code>, excluding a sign bit.
171      * For positive <code>int</code>, this is equivalent to the number of digit
172      * in the ordinary decimal representation. For negative <code>int</code>,
173      * it is equivalent to the number of digit of the positive value
174      * <code>-(i + 1)</code>.
175      *
176      * @param i the <code>int</code> value for which the digit length is returned.
177      * @return the digit length of <code>i</code>.
178      */
179     public static int digitLength(int i) {
180         if (i < 0)
181             i = -++i;
182         return (i >= 100000) ? (i >= 10000000) ? (i >= 1000000000) ? 10
183                 : (i >= 100000000) ? 9 : 8 : (i >= 1000000) ? 7 : 6
184                 : (i >= 100) ? (i >= 10000) ? 5 : (i >= 1000) ? 4 : 3
185                         : (i >= 10) ? 2 : (i >= 1) ? 1 : 0;
186     }
187
188     /**
189      * Returns the number of digits in the minimal ten's-complement
190      * representation of the specified <code>long</code>, excluding a sign bit.
191      * For positive <code>long</code>, this is equivalent to the number of digit
192      * in the ordinary decimal representation. For negative <code>long</code>,
193      * it is equivalent to the number of digit of the positive value
194      * <code>-(value + 1)</code>.
195      *
196      * @param l the <code>long</code> value for which the digit length is returned.
197      * @return the digit length of <code>l</code>.
198      */
199     public static int digitLength(long l) {
200         if (l < 0)
201             l = -++l;
202         if (l <= Integer.MAX_VALUE) {
203             int i = (int) l;
204             return (i >= 100000) ? (i >= 10000000) ? (i >= 1000000000) ? 10
205                     : (i >= 100000000) ? 9 : 8 : (i >= 1000000) ? 7 : 6
206                     : (i >= 100) ? (i >= 10000) ? 5 : (i >= 1000) ? 4 : 3
207                             : (i >= 10) ? 2 : (i >= 1) ? 1 : 0;
208         }
209         // At least 10 digits or more.
210 return (l >= 100000000000000L) ? (l >= 10000000000000000L) ? (l >= 1000000000000000000L) ? 19
211                 : (l >= 100000000000000000L) ? 18 : 17
212                 : (l >= 1000000000000000L) ? 16 : 15
213                 : (l >= 100000000000L) ? (l >= 10000000000000L) ? 14
214                         : (l >= 1000000000000L) ? 13 : 12
215                         : (l >= 10000000000L) ? 11 : 10;
216
217     }
218
219     /**
220      * Returns the closest <code>double</code> representation of the
221      * specified <code>long</code> number multiplied by a power of two.
222      *
223      * @param m the <code>long</code> multiplier.
224      * @param n the power of two exponent.
225      * @return <code>m * 2<sup>n</sup></code>.
226      @JVM-1.1+@
227      public static double toDoublePow2(long m, int n) {
228      if (m == 0)
229      return 0.0;
230      if (m == Long.MIN_VALUE)
231      return toDoublePow2(Long.MIN_VALUE >> 1, n + 1);
232      if (m < 0)
233      return -toDoublePow2(-m, n);
234      int bitLength = MathLib.bitLength(m);
235      int shift = bitLength - 53;
236      long exp = 1023L + 52 + n + shift; // Use long to avoid overflow.
237      if (exp >= 0x7FF)
238      return Double.POSITIVE_INFINITY;
239      if (exp <= 0) { // Degenerated number (subnormal, assume 0 for bit 52)
240      if (exp <= -54)
241      return 0.0;
242      return toDoublePow2(m, n + 54) / 18014398509481984L; // 2^54 Exact.
243      }
244      // Normal number.
245      long bits = (shift > 0) ? (m >> shift)
246      + ((m >> (shift - 1)) & 1) : // Rounding.
247      m << -shift;
248      if (((bits >> 52) != 1) && (++exp >= 0x7FF)) // Rare case where rounding push to next power of 2.
249      return Double.POSITIVE_INFINITY;
250      bits &= 0x000fffffffffffffL; // Clears MSB (bit 52)
251      bits |= exp << 52;
252      return Double.longBitsToDouble(bits);
253      }
254      /**/
255
256     /**
257      * Returns the closest <code>double</code> representation of the
258      * specified <code>long</code> number multiplied by a power of ten.
259      *
260      * @param m the <code>long</code> multiplier.
261      * @param n the power of ten exponent.
262      * @return <code>multiplier * 10<sup>n</sup></code>.
263      * @JVM-1.1+@
264      public static double toDoublePow10(long m, int n) {
265      if (m == 0)
266      return 0.0;
267      if (m == Long.MIN_VALUE)
268      return toDoublePow10(Long.MIN_VALUE / 10, n + 1);
269      if (m < 0)
270      return -toDoublePow10(-m, n);
271      if (n >= 0) { // Positive power.
272      if (n > 308)
273      return Double.POSITIVE_INFINITY;
274      // Works with 4 x 32 bits registers (x3:x2:x1:x0)
275      long x0 = 0; // 32 bits.
276      long x1 = 0; // 32 bits.
277      long x2 = m & MASK_32; // 32 bits.
278      long x3 = m >>> 32; // 32 bits.
279      int pow2 = 0;
280      while (n != 0) {
281      int i = (n >= POW5_INT.length) ? POW5_INT.length - 1 : n;
282      int coef = POW5_INT[i]; // 31 bits max.
283      
284      if (((int)x0) != 0) x0 *= coef; // 63 bits max.
285      if (((int)x1) != 0) x1 *= coef; // 63 bits max.
286      x2 *= coef; // 63 bits max.
287      x3 *= coef; // 63 bits max.
288      
289      x1 += x0 >>> 32;
290      x0 &= MASK_32;
291      
292      x2 += x1 >>> 32;
293      x1 &= MASK_32;
294      
295      x3 += x2 >>> 32;
296      x2 &= MASK_32;
297      
298      // Adjusts powers.
299      pow2 += i;
300      n -= i;
301      
302      // Normalizes (x3 should be 32 bits max).
303      long carry = x3 >>> 32;
304      if (carry != 0) { // Shift.
305      x0 = x1;
306      x1 = x2;
307      x2 = x3 & MASK_32;
308      x3 = carry;
309      pow2 += 32;
310      }
311      }
312      
313      // Merges registers to a 63 bits mantissa.
314      int shift = 31 - MathLib.bitLength(x3); // -1..30
315      pow2 -= shift;
316      long mantissa = (shift < 0) ?
317      (x3 << 31) | (x2 >>> 1) : // x3 is 32 bits.
318      (((x3 << 32) | x2) << shift) | (x1 >>> (32 - shift));
319      return toDoublePow2(mantissa, pow2);
320      
321      } else { // n < 0
322      if (n < -324 - 20) // m less than 20 digits.
323      return 0.0;
324      
325      // Works with x1:x0 126 bits register.
326      long x1 = m; // 63 bits.
327      long x0 = 0; // 63 bits.
328      int pow2 = 0;
329      while (true) {
330      
331      // Normalizes x1:x0
332      int shift = 63 - MathLib.bitLength(x1);
333      x1 <<= shift;
334      x1 |= x0 >>> (63 - shift);
335      x0 = (x0 << shift) & MASK_63;
336      pow2 -= shift;
337      
338      // Checks if division has to be performed.
339      if (n == 0)
340      break; // Done.
341      
342      // Retrieves power of 5 divisor.
343      int i = (-n >= POW5_INT.length) ? POW5_INT.length - 1 : -n;
344      int divisor = POW5_INT[i];
345      
346      // Performs the division (126 bits by 31 bits).
347      long wh = (x1 >>> 32);
348      long qh = wh / divisor;
349      long r = wh - qh * divisor;
350      long wl = (r << 32) | (x1 & MASK_32);
351      long ql = wl / divisor;
352      r = wl - ql * divisor;
353      x1 = (qh << 32) | ql;
354      
355      wh = (r << 31) | (x0 >>> 32);
356      qh = wh / divisor;
357      r = wh - qh * divisor;
358      wl = (r << 32) | (x0 & MASK_32);
359      ql = wl / divisor;
360      x0 = (qh << 32) | ql;
361      
362      // Adjusts powers.
363      n += i;
364      pow2 -= i;
365      }
366      return toDoublePow2(x1, pow2);
367      }
368      }
369      
370      private static final long MASK_63 = 0x7FFFFFFFFFFFFFFFL;
371      
372      private static final long MASK_32 = 0xFFFFFFFFL;
373      
374      private static final int[] POW5_INT = { 1, 5, 25, 125, 625, 3125, 15625,
375      78125, 390625, 1953125, 9765625, 48828125, 244140625, 1220703125 };
376      /**/
377
378     /**
379      * Returns the closest <code>long</code> representation of the
380      * specified <code>double</code> number multiplied by a power of two.
381      *
382      * @param d the <code>double</code> multiplier.
383      * @param n the power of two exponent.
384      * @return <code>d * 2<sup>n</sup></code>
385      * @throws ArithmeticException if the conversion cannot be performed
386      * (NaN, Infinity or overflow).
387      * @JVM-1.1+@
388      public static long toLongPow2(double d, int n) {
389      long bits = Double.doubleToRawLongBits(d);
390      boolean isNegative = (bits >> 63) != 0;
391      int exp = ((int)(bits >> 52)) & 0x7FF;
392      long m = bits & 0x000fffffffffffffL;
393      if (exp == 0x7FF) throw new ArithmeticException(
394      "Cannot convert to long (Infinity or NaN)");
395      if (exp == 0) {
396      if (m == 0) return 0L;
397      return toLongPow2(d * 18014398509481984L, n - 54); // 2^54 Exact.
398      }
399      m |= 0x0010000000000000L; // Sets MSB (bit 52)
400      long shift = exp - 1023L - 52 + n; // Use long to avoid overflow.
401      if (shift <= -64) return 0L;
402      if (shift >= 11) throw new ArithmeticException(
403      "Cannot convert to long (overflow)");
404      m = (shift >= 0) ? m << shift :
405      (m >> -shift) + ((m >> -(shift + 1)) & 1) ; // Rounding.
406      return isNegative ? -m : m;
407      }
408      /**/
409
410     /**
411      * Returns the closest <code>long</code> representation of the
412      * specified <code>double</code> number multiplied by a power of ten.
413      *
414      * @param d the <code>double</code> multiplier.
415      * @param n the power of two exponent.
416      * @return <code>d * 10<sup>n</sup></code>.
417      @JVM-1.1+@
418      public static long toLongPow10(double d, int n) {
419      long bits = Double.doubleToRawLongBits(d);
420      boolean isNegative = (bits >> 63) != 0;
421      int exp = ((int)(bits >> 52)) & 0x7FF;
422      long m = bits & 0x000fffffffffffffL;
423      if (exp == 0x7FF) throw new ArithmeticException(
424      "Cannot convert to long (Infinity or NaN)");
425      if (exp == 0) {
426      if (m == 0) return 0L;
427      return toLongPow10(d * 1E16, n - 16);
428      }
429      m |= 0x0010000000000000L; // Sets MSB (bit 52)
430      int pow2 = exp - 1023 - 52;
431      // Retrieves 63 bits m with n == 0.
432      if (n >= 0) {
433      // Works with 4 x 32 bits registers (x3:x2:x1:x0)
434      long x0 = 0; // 32 bits.
435      long x1 = 0; // 32 bits.
436      long x2 = m & MASK_32; // 32 bits.
437      long x3 = m >>> 32; // 32 bits.
438      while (n != 0) {
439      int i = (n >= POW5_INT.length) ? POW5_INT.length - 1 : n;
440      int coef = POW5_INT[i]; // 31 bits max.
441      
442      if (((int)x0) != 0) x0 *= coef; // 63 bits max.
443      if (((int)x1) != 0) x1 *= coef; // 63 bits max.
444      x2 *= coef; // 63 bits max.
445      x3 *= coef; // 63 bits max.
446      
447      x1 += x0 >>> 32;
448      x0 &= MASK_32;
449      
450      x2 += x1 >>> 32;
451      x1 &= MASK_32;
452      
453      x3 += x2 >>> 32;
454      x2 &= MASK_32;
455      
456      // Adjusts powers.
457      pow2 += i;
458      n -= i;
459      
460      // Normalizes (x3 should be 32 bits max).
461      long carry = x3 >>> 32;
462      if (carry != 0) { // Shift.
463      x0 = x1;
464      x1 = x2;
465      x2 = x3 & MASK_32;
466      x3 = carry;
467      pow2 += 32;
468      }
469      }
470      
471      // Merges registers to a 63 bits mantissa.
472      int shift = 31 - MathLib.bitLength(x3); // -1..30
473      pow2 -= shift;
474      m = (shift < 0) ?
475      (x3 << 31) | (x2 >>> 1) : // x3 is 32 bits.
476      (((x3 << 32) | x2) << shift) | (x1 >>> (32 - shift));
477      
478      } else { // n < 0
479      
480      // Works with x1:x0 126 bits register.
481      long x1 = m; // 63 bits.
482      long x0 = 0; // 63 bits.
483      while (true) {
484      
485      // Normalizes x1:x0
486      int shift = 63 - MathLib.bitLength(x1);
487      x1 <<= shift;
488      x1 |= x0 >>> (63 - shift);
489      x0 = (x0 << shift) & MASK_63;
490      pow2 -= shift;
491      
492      // Checks if division has to be performed.
493      if (n == 0)
494      break; // Done.
495      
496      // Retrieves power of 5 divisor.
497      int i = (-n >= POW5_INT.length) ? POW5_INT.length - 1 : -n;
498      int divisor = POW5_INT[i];
499      
500      // Performs the division (126 bits by 31 bits).
501      long wh = (x1 >>> 32);
502      long qh = wh / divisor;
503      long r = wh - qh * divisor;
504      long wl = (r << 32) | (x1 & MASK_32);
505      long ql = wl / divisor;
506      r = wl - ql * divisor;
507      x1 = (qh << 32) | ql;
508      
509      wh = (r << 31) | (x0 >>> 32);
510      qh = wh / divisor;
511      r = wh - qh * divisor;
512      wl = (r << 32) | (x0 & MASK_32);
513      ql = wl / divisor;
514      x0 = (qh << 32) | ql;
515      
516      // Adjusts powers.
517      n += i;
518      pow2 -= i;
519      }
520      m = x1;
521      }
522      if (pow2 > 0) throw new ArithmeticException("Overflow");
523      if (pow2 < -63) return 0;
524      m = (m >> -pow2) + ((m >> -(pow2 + 1)) & 1) ; // Rounding.
525      return isNegative ? -m : m;
526      }
527      /**/
528
529     /**
530      * Returns the largest power of 2 that is less than or equal to the
531      * the specified positive value.
532      *
533      * @param d the <code>double</code> number.
534      * @return <code>floor(Log2(abs(d)))</code>
535      * @throws ArithmeticException if <code>d &lt;= 0<code> or <code>d</code>
536      * is <code>NaN</code> or <code>Infinity</code>.
537      * @JVM-1.1+@
538      public static int floorLog2(double d) {
539      if (d <= 0) throw new ArithmeticException("Negative number or zero");
540      long bits = Double.doubleToRawLongBits(d);
541      int exp = ((int)(bits >> 52)) & 0x7FF;
542      if (exp == 0x7FF) throw new ArithmeticException("Infinity or NaN");
543      if (exp == 0) // Degenerated.
544      return floorLog2(d * 18014398509481984L) - 54; // 2^54 Exact.
545      return exp - 1023;
546      }
547      /**/
548
549     /**
550      * Returns the largest power of 10 that is less than or equal to the
551      * the specified positive value.
552      *
553      * @param d the <code>double</code> number.
554      * @return <code>floor(Log10(abs(d)))</code>
555      * @throws ArithmeticException if <code>d &lt;= 0<code> or <code>d</code>
556      * is <code>NaN</code> or <code>Infinity</code>.
557      * @JVM-1.1+@
558      public static int floorLog10(double d) {
559      int guess = (int) (LOG2_DIV_LOG10 * MathLib.floorLog2(d));
560      double pow10 = MathLib.toDoublePow10(1, guess);
561      if ((pow10 <= d) && (pow10 * 10 > d)) // Good guess.
562      return guess;
563      if (pow10 > d) // Too large.
564      return guess - 1;
565      return guess + 1;
566      }
567      private static final double LOG2_DIV_LOG10 = 0.3010299956639811952137388947;
568      /**/
569
570     /**
571      * The natural logarithm.
572      * @JVM-1.1+@
573      public static final double E = 2.71828182845904523536028747135266;
574      
575      /**
576      * The ratio of the circumference of a circle to its diameter.
577      * @JVM-1.1+@
578      public static final double PI = 3.1415926535897932384626433832795;
579
580      /**
581      * Half the ratio of the circumference of a circle to its diameter.
582      * @JVM-1.1+@
583      public static final double HALF_PI = 1.5707963267948966192313216916398;
584
585      /**
586      * Twice the ratio of the circumference of a circle to its diameter.
587      * @JVM-1.1+@
588      public static final double TWO_PI = 6.283185307179586476925286766559;
589      
590      /**
591      * Four time the ratio of the circumference of a circle to its diameter.
592      * @JVM-1.1+@
593      public static final double FOUR_PI = 12.566370614359172953850573533118;
594      
595      /**
596      * Holds {@link #PI} * {@link #PI}.
597      * @JVM-1.1+@
598      public static final double PI_SQUARE = 9.8696044010893586188344909998762;
599      
600      /**
601      * The natural logarithm of two.
602      * @JVM-1.1+@
603      public static final double LOG2 = 0.69314718055994530941723212145818;
604      
605      /**
606      * The natural logarithm of ten.
607      * @JVM-1.1+@
608      public static final double LOG10 = 2.3025850929940456840179914546844;
609      
610      /**
611      * The square root of two.
612      * @JVM-1.1+@
613      public static final double SQRT2 = 1.4142135623730950488016887242097;
614      
615      /**
616      * Not-A-Number.
617      * @JVM-1.1+@
618      public static final double NaN = 0.0 / 0.0;
619
620      /**
621      * Infinity.
622      * @JVM-1.1+@
623      public static final double Infinity = 1.0 / 0.0;
624
625      /**/
626
627     /**
628      * Converts an angle in degrees to radians.
629      *
630      * @param degrees the angle in degrees.
631      * @return the specified angle in radians.
632      * @JVM-1.1+@
633      public static double toRadians(double degrees) {
634      return degrees * (PI / 180.0);
635      }
636      /**/
637
638     /**
639      * Converts an angle in radians to degrees.
640      *
641      * @param radians the angle in radians.
642      * @return the specified angle in degrees.
643      * @JVM-1.1+@
644      public static double toDegrees(double radians) {
645      return radians * (180.0 / PI);
646      }
647      /**/
648
649     /**
650      * Returns the positive square root of the specified value.
651      *
652      * @param x the value.
653      * @return <code>java.lang.Math.sqrt(x)</code>
654      * @JVM-1.1+@
655      public static double sqrt(double x) {
656      return Math.sqrt(x); // CLDC 1.1
657      }
658      /**/
659
660     /**
661      * Returns the remainder of the division of the specified two arguments.
662      *
663      * @param x the dividend.
664      * @param y the divisor.
665      * @return <code>x - round(x / y) * y</code>
666      * @JVM-1.1+@
667      public static double rem(double x, double y) {
668      double tmp = x / y;
669      if (MathLib.abs(tmp) <= Long.MAX_VALUE) {
670      return x - MathLib.round(tmp) * y;
671      } else {
672      return NaN;
673      }
674      }
675      /**/
676
677     /**
678      * Returns the smallest (closest to negative infinity)
679      * <code>double</code> value that is not less than the argument and is
680      * equal to a mathematical integer.
681      *
682      * @param x the value.
683      * @return <code>java.lang.Math.ceil(x)</code>
684      * @JVM-1.1+@
685      public static double ceil(double x) {
686      return Math.ceil(x); // CLDC 1.1
687      }
688      /**/
689
690     /**
691      * Returns the largest (closest to positive infinity)
692      * <code>double</code> value that is not greater than the argument and
693      * is equal to a mathematical integer.
694      *
695      * @param x the value.
696      * @return <code>java.lang.Math.ceil(x)</code>
697      * @JVM-1.1+@
698      public static double floor(double x) {
699      return Math.floor(x); // CLDC 1.1
700      }
701      /**/
702
703     /**
704      * Returns the trigonometric sine of the specified angle in radians.
705      *
706      * @param radians the angle in radians.
707      * @return <code>java.lang.Math.sin(radians)</code>
708      * @JVM-1.1+@
709      public static double sin(double radians) {
710      return Math.sin(radians); // CLDC 1.1
711      }
712      /**/
713
714     /**
715      * Returns the trigonometric cosine of the specified angle in radians.
716      *
717      * @param radians the angle in radians.
718      * @return <code>java.lang.Math.cos(radians)</code>
719      * @JVM-1.1+@
720      public static double cos(double radians) {
721      return Math.cos(radians); // CLDC 1.1
722      }
723      /**/
724
725     /**
726      * Returns the trigonometric tangent of the specified angle in radians.
727      *
728      * @param radians the angle in radians.
729      * @return <code>java.lang.Math.tan(radians)</code>
730      * @JVM-1.1+@
731      public static double tan(double radians) {
732      return Math.tan(radians); // CLDC 1.1
733      }
734      /**/
735
736     /**
737      * Returns the arc sine of the specified value,
738      * in the range of -<i>pi</i>/2 through <i>pi</i>/2.
739      *
740      * @param x the value whose arc sine is to be returned.
741      * @return the arc sine in radians for the specified value.
742      * @JVM-1.1+@
743      public static double asin(double x) {
744      if (x < -1.0 || x > 1.0) return MathLib.NaN;
745      if (x == -1.0) return - HALF_PI;
746      if (x == 1.0) return HALF_PI;
747      return MathLib.atan(x / MathLib.sqrt(1.0 - x * x));
748      }
749      /**/
750
751     /**
752      * Returns the arc cosine of the specified value,
753      * in the range of 0.0 through <i>pi</i>.
754      *
755      * @param x the value whose arc cosine is to be returned.
756      * @return the arc cosine in radians for the specified value.
757      * @JVM-1.1+@
758      public static double acos(double x) {
759      return HALF_PI - MathLib.asin(x);
760      }
761      /**/
762
763     /**
764      * Returns the arc tangent of the specified value,
765      * in the range of -<i>pi</i>/2 through <i>pi</i>/2.
766      *
767      * @param x the value whose arc tangent is to be returned.
768      * @return the arc tangent in radians for the specified value.
769      * @see <a HREF="http://mathworld.wolfram.com/InverseTangent.html">
770      * Inverse Tangent -- from MathWorld</a>
771      * @JVM-1.1+@
772      public static double atan(double x) {
773      return MathLib._atan(x);
774      }
775      /**/
776
777     /**
778      * Returns the angle theta such that
779      * <code>(x == cos(theta)) && (y == sin(theta))</code>.
780      *
781      * @param y the y value.
782      * @param x the x value.
783      * @return the angle theta in radians.
784      * @JVM-1.1+@
785      public static double atan2(double y, double x) {
786      final double epsilon = 1E-128;
787      if (MathLib.abs(x) > epsilon) {
788      double temp = MathLib.atan(MathLib.abs(y) / MathLib.abs(x));
789      if( x < 0.0 ) temp = PI - temp;
790      if( y < 0.0 ) temp = TWO_PI - temp;
791      return temp;
792      } else if( y > epsilon) {
793      return HALF_PI;
794      } else if( y < -epsilon) {
795      return 3 * HALF_PI;
796      } else {
797      return 0.0;
798      }
799      }
800      /**/
801
802     /**
803      * Returns the hyperbolic sine of x.
804      *
805      * @param x the value for which the hyperbolic sine is calculated.
806      * @return <code>(exp(x) - exp(-x)) / 2</code>
807      * @JVM-1.1+@
808      public static double sinh(double x) {
809      return (MathLib.exp(x) - MathLib.exp(-x)) * 0.5;
810      }
811      /**/
812
813     /**
814      * Returns the hyperbolic cosine of x.
815      *
816      * @param x the value for which the hyperbolic cosine is calculated.
817      * @return <code>(exp(x) + exp(-x)) / 2</code>
818      * @JVM-1.1+@
819      public static double cosh(double x) {
820      return (MathLib.exp(x) + MathLib.exp(-x)) * 0.5;
821      }
822      /**/
823
824     /**
825      * Returns the hyperbolic tangent of x.
826      *
827      * @param x the value for which the hyperbolic tangent is calculated.
828      * @return <code>(exp(2 * x) - 1) / (exp(2 * x) + 1)</code>
829      * @JVM-1.1+@
830      public static double tanh(double x) {
831      return (MathLib.exp(2 * x) - 1) / (MathLib.exp(2 * x) + 1);
832      }
833      /**/
834
835     /**
836      * Returns <i>{@link #E e}</i> raised to the specified power.
837      *
838      * @param x the exponent.
839      * @return <code><i>e</i><sup>x</sup></code>
840      * @see <a HREF="http://mathworld.wolfram.com/ExponentialFunction.html">
841      * Exponential Function -- from MathWorld</a>
842      * @JVM-1.1+@
843      public static double exp(double x) {
844      return MathLib._ieee754_exp(x);
845      }
846      /**/
847
848     /**
849      * Returns the natural logarithm (base <i>{@link #E e}</i>) of the specified
850      * value.
851      *
852      * @param x the value greater than <code>0.0</code>.
853      * @return the value y such as <code><i>e</i><sup>y</sup> == x</code>
854      * @JVM-1.1+@
855      public static double log(double x) {
856      return MathLib._ieee754_log(x);
857      }
858      /**/
859
860     /**
861      * Returns the decimal logarithm of the specified value.
862      *
863      * @param x the value greater than <code>0.0</code>.
864      * @return the value y such as <code>10<sup>y</sup> == x</code>
865      * @JVM-1.1+@
866      public static double log10(double x) {
867      return log(x) * INV_LOG10;
868      }
869      private static double INV_LOG10 = 0.43429448190325182765112891891661;
870      /**/
871
872     /**
873      * Returns the value of the first argument raised to the power of the
874      * second argument.
875      *
876      * @param x the base.
877      * @param y the exponent.
878      * @return <code>x<sup>y</sup></code>
879      * @JVM-1.1+@
880      public static double pow(double x, double y) {
881      /**/
882     /* @JVM-1.4+@ // Use java.lang.Math value.
883      if (true) return Math.pow(x, y);
884      /**/
885     /* @JVM-1.1+@ // Else (J2ME) use close approximation (+/- LSB)
886      if ((x < 0) && (y == (int)y)) return
887      (((int)y) & 1) == 0 ? pow(-x, y) : -pow(-x, y);
888      return MathLib.exp(y * MathLib.log(x));
889      }
890      /**/
891
892     /**
893      * Returns the closest <code>int</code> to the specified argument.
894      *
895      * @param f the <code>float</code> value to be rounded to a <code>int</code>
896      * @return the nearest <code>int</code> value.
897      * @JVM-1.1+@
898      public static int round(float f) {
899      return (int) floor(f + 0.5f);
900      }
901      /**/
902
903     /**
904      * Returns the closest <code>long</code> to the specified argument.
905      *
906      * @param d the <code>double</code> value to be rounded to a
907      * <code>long</code>
908      * @return the nearest <code>long</code> value.
909      * @JVM-1.1+@
910      public static long round(double d) {
911      return (long) floor(d + 0.5d);
912      }
913      /**/
914
915     /**
916      * Returns a random number between zero and one.
917      *
918      * @return a <code>double</code> greater than or equal
919      * to <code>0.0</code> and less than <code>1.0</code>.
920      * @JVM-1.1+@
921      public static double random() {
922      return random(0, Integer.MAX_VALUE) * NORMALIZATION_FACTOR;
923      }
924      private static final double NORMALIZATION_FACTOR = -1.0 / Integer.MIN_VALUE;
925      /**/
926
927     /**
928      * Returns the absolute value of the specified <code>int</code> argument.
929      *
930      * @param i the <code>int</code> value.
931      * @return <code>i</code> or <code>-i</code>
932      */
933     public static int abs(int i) {
934         return (i < 0) ? -i : i;
935     }
936
937     /**
938      * Returns the absolute value of the specified <code>long</code> argument.
939      *
940      * @param l the <code>long</code> value.
941      * @return <code>l</code> or <code>-l</code>
942      */
943     public static long abs(long l) {
944         return (l < 0) ? -l : l;
945     }
946
947     /**
948      * Returns the absolute value of the specified <code>float</code> argument.
949      *
950      * @param f the <code>float</code> value.
951      * @return <code>f</code> or <code>-f</code>
952      * @JVM-1.1+@
953      public static float abs(float f) {
954      return (f < 0) ? -f : f;
955      }
956      /**/
957
958     /**
959      * Returns the absolute value of the specified <code>double</code> argument.
960      *
961      * @param d the <code>double</code> value.
962      * @return <code>d</code> or <code>-d</code>
963      * @JVM-1.1+@
964      public static double abs(double d) {
965      return (d < 0) ? -d : d;
966      }
967      /**/
968
969     /**
970      * Returns the greater of two <code>int</code> values.
971      *
972      * @param x the first value.
973      * @param y the second value.
974      * @return the larger of <code>x</code> and <code>y</code>.
975      */
976     public static int max(int x, int y) {
977         return (x >= y) ? x : y;
978     }
979
980     /**
981      * Returns the greater of two <code>long</code> values.
982      *
983      * @param x the first value.
984      * @param y the second value.
985      * @return the larger of <code>x</code> and <code>y</code>.
986      */
987     public static long max(long x, long y) {
988         return (x >= y) ? x : y;
989     }
990
991     /**
992      * Returns the greater of two <code>float</code> values.
993      *
994      * @param x the first value.
995      * @param y the second value.
996      * @return the larger of <code>x</code> and <code>y</code>.
997      * @JVM-1.1+@
998      public static float max(float x, float y) {
999      return (x >= y) ? x : y;
1000     }
1001     /**/
1002
1003    /**
1004     * Returns the greater of two <code>double</code> values.
1005     *
1006     * @param x the first value.
1007     * @param y the second value.
1008     * @return the larger of <code>x</code> and <code>y</code>.
1009     * @JVM-1.1+@
1010     public static double max(double x, double y) {
1011     return (x >= y) ? x : y;
1012     }
1013     /**/
1014
1015    /**
1016     * Returns the smaller of two <code>int</code> values.
1017     *
1018     * @param x the first value.
1019     * @param y the second value.
1020     * @return the smaller of <code>x</code> and <code>y</code>.
1021     */
1022    public static int min(int x, int y) {
1023        return (x < y) ? x : y;
1024    }
1025
1026    /**
1027     * Returns the smaller of two <code>long</code> values.
1028     *
1029     * @param x the first value.
1030     * @param y the second value.
1031     * @return the smaller of <code>x</code> and <code>y</code>.
1032     */
1033    public static long min(long x, long y) {
1034        return (x < y) ? x : y;
1035    }
1036
1037    /**
1038     * Returns the smaller of two <code>float</code> values.
1039     *
1040     * @param x the first value.
1041     * @param y the second value.
1042     * @return the smaller of <code>x</code> and <code>y</code>.
1043     * @JVM-1.1+@
1044     public static float min(float x, float y) {
1045     return (x < y) ? x : y;
1046     }
1047     /**/
1048
1049    /**
1050     * Returns the smaller of two <code>double</code> values.
1051     *
1052     * @param x the first value.
1053     * @param y the second value.
1054     * @return the smaller of <code>x</code> and <code>y</code>.
1055     * @JVM-1.1+@
1056     public static double min(double x, double y) {
1057     return (x < y) ? x : y;
1058     }
1059     /**/
1060
1061    ////////////////////////////////////////////////////////////////////////////
1062 /* @(#)s_atan.c 1.3 95/01/18 */
1063    /*
1064     * ====================================================
1065     * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
1066     *
1067     * Developed at SunSoft, a Sun Microsystems, Inc. business.
1068     * Permission to use, copy, modify, and distribute this
1069     * software is freely granted, provided that this notice
1070     * is preserved.
1071     * ====================================================
1072     *
1073     */
1074
1075    /* atan(x)
1076     * Method
1077     * 1. Reduce x to positive by atan(x) = -atan(-x).
1078     * 2. According to the integer k=4t+0.25 chopped, t=x, the argument
1079     * is further reduced to one of the following intervals and the
1080     * arctangent of t is evaluated by the corresponding formula:
1081     *
1082     * [0,7/16] atan(x) = t-t^3*(a1+t^2*(a2+...(a10+t^2*a11)...)
1083     * [7/16,11/16] atan(x) = atan(1/2) + atan( (t-0.5)/(1+t/2) )
1084     * [11/16.19/16] atan(x) = atan( 1 ) + atan( (t-1)/(1+t) )
1085     * [19/16,39/16] atan(x) = atan(3/2) + atan( (t-1.5)/(1+1.5t) )
1086     * [39/16,INF] atan(x) = atan(INF) + atan( -1/t )
1087     *
1088     * Constants:
1089     * The hexadecimal values are the intended ones for the following
1090     * constants. The decimal values may be used, provided that the
1091     * compiler will convert from decimal to binary accurately enough
1092     * to produce the hexadecimal values shown.
1093     @JVM-1.1+@
1094     static final double atanhi[] = {
1095     4.63647609000806093515e-01, // atan(0.5)hi 0x3FDDAC67, 0x0561BB4F
1096     7.85398163397448278999e-01, // atan(1.0)hi 0x3FE921FB, 0x54442D18
1097     9.82793723247329054082e-01, // atan(1.5)hi 0x3FEF730B, 0xD281F69B
1098     1.57079632679489655800e+00, // atan(inf)hi 0x3FF921FB, 0x54442D18
1099     };
1100     static final double atanlo[] = {
1101     2.26987774529616870924e-17, // atan(0.5)lo 0x3C7A2B7F, 0x222F65E2
1102     3.06161699786838301793e-17, // atan(1.0)lo 0x3C81A626, 0x33145C07
1103     1.39033110312309984516e-17, // atan(1.5)lo 0x3C700788, 0x7AF0CBBD
1104     6.12323399573676603587e-17, // atan(inf)lo 0x3C91A626, 0x33145C07
1105     };
1106     static final double aT[] = {
1107     3.33333333333329318027e-01, // 0x3FD55555, 0x5555550D
1108     -1.99999999998764832476e-01, // 0xBFC99999, 0x9998EBC4
1109     1.42857142725034663711e-01, // 0x3FC24924, 0x920083FF
1110     -1.11111104054623557880e-01, // 0xBFBC71C6, 0xFE231671
1111     9.09088713343650656196e-02, // 0x3FB745CD, 0xC54C206E
1112     -7.69187620504482999495e-02, // 0xBFB3B0F2, 0xAF749A6D
1113     6.66107313738753120669e-02, // 0x3FB10D66, 0xA0D03D51
1114     -5.83357013379057348645e-02, // 0xBFADDE2D, 0x52DEFD9A
1115     4.97687799461593236017e-02, // 0x3FA97B4B, 0x24760DEB
1116     -3.65315727442169155270e-02, // 0xBFA2B444, 0x2C6A6C2F
1117     1.62858201153657823623e-02, // 0x3F90AD3A, 0xE322DA11
1118     };
1119     static final double
1120     one = 1.0,
1121     huge = 1.0e300;
1122     static double _atan(double x)
1123     {
1124     double w,s1,s2,z;
1125     int ix,hx,id;
1126     long xBits = Double.doubleToLongBits(x);
1127     int __HIx = (int) (xBits >> 32);
1128     int __LOx = (int) xBits;
1129
1130     hx = __HIx;
1131     ix = hx&0x7fffffff;
1132     if(ix>=0x44100000) { // if |x| >= 2^66
1133     if(ix>0x7ff00000||
1134     (ix==0x7ff00000&&(__LOx!=0)))
1135     return x+x; // NaN
1136     if(hx>0) return atanhi[3]+atanlo[3];
1137     else return -atanhi[3]-atanlo[3];
1138     } if (ix < 0x3fdc0000) { // |x| < 0.4375
1139     if (ix < 0x3e200000) { // |x| < 2^-29
1140     if(huge+x>one) return x; // raise inexact
1141     }
1142     id = -1;
1143     } else {
1144     x = MathLib.abs(x);
1145     if (ix < 0x3ff30000) { // |x| < 1.1875
1146     if (ix < 0x3fe60000) { // 7/16 <=|x|<11/16
1147     id = 0; x = (2.0*x-one)/(2.0+x);
1148     } else { // 11/16<=|x|< 19/16
1149     id = 1; x = (x-one)/(x+one);
1150     }
1151     } else {
1152     if (ix < 0x40038000) { // |x| < 2.4375
1153     id = 2; x = (x-1.5)/(one+1.5*x);
1154     } else { // 2.4375 <= |x| < 2^66
1155     id = 3; x = -1.0/x;
1156     }
1157     }}
1158     // end of argument reduction
1159     z = x*x;
1160     w = z*z;
1161     // break sum from i=0 to 10 aT[i]z**(i+1) into odd and even poly
1162     s1 = z*(aT[0]+w*(aT[2]+w*(aT[4]+w*(aT[6]+w*(aT[8]+w*aT[10])))));
1163     s2 = w*(aT[1]+w*(aT[3]+w*(aT[5]+w*(aT[7]+w*aT[9]))));
1164     if (id<0) return x - x*(s1+s2);
1165     else {
1166     z = atanhi[id] - ((x*(s1+s2) - atanlo[id]) - x);
1167     return (hx<0)? -z:z;
1168     }
1169     }
1170     /**/
1171
1172    ////////////////////////////////////////////////////////////////////////////
1173 /* @(#)e_log.c 1.3 95/01/18 */
1174    /*
1175     * ====================================================
1176     * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
1177     *
1178     * Developed at SunSoft, a Sun Microsystems, Inc. business.
1179     * Permission to use, copy, modify, and distribute this
1180     * software is freely granted, provided that this notice
1181     * is preserved.
1182     * ====================================================
1183     */
1184
1185    /* __ieee754_log(x)
1186     * Return the logrithm of x
1187     *
1188     * Method :
1189     * 1. Argument Reduction: find k and f such that
1190     * x = 2^k * (1+f),
1191     * where sqrt(2)/2 < 1+f < sqrt(2) .
1192     *
1193     * 2. Approximation of log(1+f).
1194     * Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
1195     * = 2s + 2/3 s**3 + 2/5 s**5 + .....,
1196     * = 2s + s*R
1197     * We use a special Reme algorithm on [0,0.1716] to generate
1198     * a polynomial of degree 14 to approximate R The maximum error
1199     * of this polynomial approximation is bounded by 2**-58.45. In
1200     * other words,
1201     * 2 4 6 8 10 12 14
1202     * R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s +Lg6*s +Lg7*s
1203     * (the values of Lg1 to Lg7 are listed in the program)
1204     * and
1205     * | 2 14 | -58.45
1206     * | Lg1*s +...+Lg7*s - R(z) | <= 2
1207     * | |
1208     * Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2.
1209     * In order to guarantee error in log below 1ulp, we compute log
1210     * by
1211     * log(1+f) = f - s*(f - R) (if f is not too large)
1212     * log(1+f) = f - (hfsq - s*(hfsq+R)). (better accuracy)
1213     *
1214     * 3. Finally, log(x) = k*ln2 + log(1+f).
1215     * = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo)))
1216     * Here ln2 is split into two floating point number:
1217     * ln2_hi + ln2_lo,
1218     * where n*ln2_hi is always exact for |n| < 2000.
1219     *
1220     * Special cases:
1221     * log(x) is NaN with signal if x < 0 (including -INF) ;
1222     * log(+INF) is +INF; log(0) is -INF with signal;
1223     * log(NaN) is that NaN with no signal.
1224     *
1225     * Accuracy:
1226     * according to an error analysis, the error is always less than
1227     * 1 ulp (unit in the last place).
1228     *
1229     * Constants:
1230     * The hexadecimal values are the intended ones for the following
1231     * constants. The decimal values may be used, provided that the
1232     * compiler will convert from decimal to binary accurately enough
1233     * to produce the hexadecimal values shown.
1234     @JVM-1.1+@
1235     static final double
1236     ln2_hi = 6.93147180369123816490e-01, // 3fe62e42 fee00000
1237     ln2_lo = 1.90821492927058770002e-10, // 3dea39ef 35793c76
1238     two54 = 1.80143985094819840000e+16, // 43500000 00000000
1239     Lg1 = 6.666666666666735130e-01, // 3FE55555 55555593
1240     Lg2 = 3.999999999940941908e-01, // 3FD99999 9997FA04
1241     Lg3 = 2.857142874366239149e-01, // 3FD24924 94229359
1242     Lg4 = 2.222219843214978396e-01, // 3FCC71C5 1D8E78AF
1243     Lg5 = 1.818357216161805012e-01, // 3FC74664 96CB03DE
1244     Lg6 = 1.531383769920937332e-01, // 3FC39A09 D078C69F
1245     Lg7 = 1.479819860511658591e-01; // 3FC2F112 DF3E5244
1246     static final double zero = 0.0;
1247     static double _ieee754_log(double x)
1248     {
1249     double hfsq,f,s,z,R,w,t1,t2,dk;
1250     int k,hx,i,j;
1251     int lx; // unsigned
1252
1253     long xBits = Double.doubleToLongBits(x);
1254     hx = (int) (xBits >> 32);
1255     lx = (int) xBits;
1256
1257     k=0;
1258     if (hx < 0x00100000) { // x < 2**-1022
1259     if (((hx&0x7fffffff)|lx)==0)
1260     return -two54/zero; // log(+-0)=-inf
1261     if (hx<0) return (x-x)/zero; // log(-#) = NaN
1262     k -= 54; x *= two54; // subnormal number, scale up x
1263     xBits = Double.doubleToLongBits(x);
1264     hx = (int) (xBits >> 32); // high word of x
1265     }
1266     if (hx >= 0x7ff00000) return x+x;
1267     k += (hx>>20)-1023;
1268     hx &= 0x000fffff;
1269     i = (hx+0x95f64)&0x100000;
1270     xBits = Double.doubleToLongBits(x);
1271     int HIx = hx|(i^0x3ff00000); // normalize x or x/2
1272     xBits = ((HIx & 0xFFFFFFFFL) << 32) | (xBits & 0xFFFFFFFFL);
1273     x = Double.longBitsToDouble(xBits);
1274     k += (i>>20);
1275     f = x-1.0;
1276     if((0x000fffff&(2+hx))<3) { // |f| < 2**-20
1277     if(f==zero) if(k==0) return zero; else {dk=(double)k;
1278     return dk*ln2_hi+dk*ln2_lo;}
1279     R = f*f*(0.5-0.33333333333333333*f);
1280     if(k==0) return f-R; else {dk=(double)k;
1281     return dk*ln2_hi-((R-dk*ln2_lo)-f);}
1282     }
1283     s = f/(2.0+f);
1284     dk = (double)k;
1285     z = s*s;
1286     i = hx-0x6147a;
1287     w = z*z;
1288     j = 0x6b851-hx;
1289     t1= w*(Lg2+w*(Lg4+w*Lg6));
1290     t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7)));
1291     i |= j;
1292     R = t2+t1;
1293     if(i>0) {
1294     hfsq=0.5*f*f;
1295     if(k==0) return f-(hfsq-s*(hfsq+R)); else
1296     return dk*ln2_hi-((hfsq-(s*(hfsq+R)+dk*ln2_lo))-f);
1297     } else {
1298     if(k==0) return f-s*(f-R); else
1299     return dk*ln2_hi-((s*(f-R)-dk*ln2_lo)-f);
1300     }
1301     }
1302     /**/
1303
1304    ////////////////////////////////////////////////////////////////////////////
1305 /* @(#)e_exp.c 1.6 04/04/22 */
1306    /*
1307     * ====================================================
1308     * Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved.
1309     *
1310     * Permission to use, copy, modify, and distribute this
1311     * software is freely granted, provided that this notice
1312     * is preserved.
1313     * ====================================================
1314     */
1315
1316    /* __ieee754_exp(x)
1317     * Returns the exponential of x.
1318     *
1319     * Method
1320     * 1. Argument reduction:
1321     * Reduce x to an r so that |r| <= 0.5*ln2 ~ 0.34658.
1322     * Given x, find r and integer k such that
1323     *
1324     * x = k*ln2 + r, |r| <= 0.5*ln2.
1325     *
1326     * Here r will be represented as r = hi-lo for better
1327     * accuracy.
1328     *
1329     * 2. Approximation of exp(r) by a special rational function on
1330     * the interval [0,0.34658]:
1331     * Write
1332     * R(r**2) = r*(exp(r)+1)/(exp(r)-1) = 2 + r*r/6 - r**4/360 + ...
1333     * We use a special Remes algorithm on [0,0.34658] to generate
1334     * a polynomial of degree 5 to approximate R. The maximum error
1335     * of this polynomial approximation is bounded by 2**-59. In
1336     * other words,
1337     * R(z) ~ 2.0 + P1*z + P2*z**2 + P3*z**3 + P4*z**4 + P5*z**5
1338     * (where z=r*r, and the values of P1 to P5 are listed below)
1339     * and
1340     * | 5 | -59
1341     * | 2.0+P1*z+...+P5*z - R(z) | <= 2
1342     * | |
1343     * The computation of exp(r) thus becomes
1344     * 2*r
1345     * exp(r) = 1 + -------
1346     * R - r
1347     * r*R1(r)
1348     * = 1 + r + ----------- (for better accuracy)
1349     * 2 - R1(r)
1350     * where
1351     * 2 4 10
1352     * R1(r) = r - (P1*r + P2*r + ... + P5*r ).
1353     *
1354     * 3. Scale back to obtain exp(x):
1355     * From step 1, we have
1356     * exp(x) = 2^k * exp(r)
1357     *
1358     * Special cases:
1359     * exp(INF) is INF, exp(NaN) is NaN;
1360     * exp(-INF) is 0, and
1361     * for finite argument, only exp(0)=1 is exact.
1362     *
1363     * Accuracy:
1364     * according to an error analysis, the error is always less than
1365     * 1 ulp (unit in the last place).
1366     *
1367     * Misc. info.
1368     * For IEEE double
1369     * if x > 7.09782712893383973096e+02 then exp(x) overflow
1370     * if x < -7.45133219101941108420e+02 then exp(x) underflow
1371     *
1372     * Constants:
1373     * The hexadecimal values are the intended ones for the following
1374     * constants. The decimal values may be used, provided that the
1375     * compiler will convert from decimal to binary accurately enough
1376     * to produce the hexadecimal values shown.
1377     @JVM-1.1+@
1378     static final double
1379     halF[] = {0.5,-0.5,},
1380     twom1000= 9.33263618503218878990e-302, // 2**-1000=0x01700000,0
1381     o_threshold= 7.09782712893383973096e+02, // 0x40862E42, 0xFEFA39EF
1382     u_threshold= -7.45133219101941108420e+02, // 0xc0874910, 0xD52D3051
1383     ln2HI[] ={ 6.93147180369123816490e-01, // 0x3fe62e42, 0xfee00000
1384     -6.93147180369123816490e-01,},// 0xbfe62e42, 0xfee00000
1385     ln2LO[] ={ 1.90821492927058770002e-10, // 0x3dea39ef, 0x35793c76
1386     -1.90821492927058770002e-10,},// 0xbdea39ef, 0x35793c76
1387     invln2 = 1.44269504088896338700e+00, // 0x3ff71547, 0x652b82fe
1388     P1 = 1.66666666666666019037e-01, // 0x3FC55555, 0x5555553E
1389     P2 = -2.77777777770155933842e-03, // 0xBF66C16C, 0x16BEBD93
1390     P3 = 6.61375632143793436117e-05, // 0x3F11566A, 0xAF25DE2C
1391     P4 = -1.65339022054652515390e-06, // 0xBEBBBD41, 0xC5D26BF1
1392     P5 = 4.13813679705723846039e-08; // 0x3E663769, 0x72BEA4D0
1393     static double _ieee754_exp(double x) // default IEEE double exp
1394     {
1395     double y,hi = 0,lo = 0,c,t;
1396     int k = 0,xsb;
1397     int hx; // Unsigned.
1398     long xBits = Double.doubleToLongBits(x);
1399     int __HIx = (int) (xBits >> 32);
1400     int __LOx = (int) xBits;
1401
1402     hx = __HIx; // high word of x
1403     xsb = (hx>>31)&1; // sign bit of x
1404     hx &= 0x7fffffff; // high word of |x|
1405
1406     // filter out non-finite argument
1407     if(hx >= 0x40862E42) { // if |x|>=709.78...
1408     if(hx>=0x7ff00000) {
1409     if(((hx&0xfffff)|__LOx)!=0)
1410     return x+x; // NaN
1411     else return (xsb==0)? x:0.0; // exp(+-inf)={inf,0}
1412     }
1413     if(x > o_threshold) return huge*huge; // overflow
1414     if(x < u_threshold) return twom1000*twom1000; // underflow
1415     }
1416
1417     // argument reduction
1418     if(hx > 0x3fd62e42) { // if |x| > 0.5 ln2
1419     if(hx < 0x3FF0A2B2) { // and |x| < 1.5 ln2
1420     hi = x-ln2HI[xsb]; lo=ln2LO[xsb]; k = 1-xsb-xsb;
1421     } else {
1422     k = (int)(invln2*x+halF[xsb]);
1423     t = k;
1424     hi = x - t*ln2HI[0]; // t*ln2HI is exact here
1425     lo = t*ln2LO[0];
1426     }
1427     x = hi - lo;
1428     }
1429     else if(hx < 0x3e300000) { // when |x|<2**-28
1430     if(huge+x>one) return one+x;// trigger inexact
1431     }
1432     else k = 0;
1433
1434     // x is now in primary range
1435     t = x*x;
1436     c = x - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
1437     if(k==0) return one-((x*c)/(c-2.0)-x);
1438     else y = one-((lo-(x*c)/(2.0-c))-hi);
1439     long yBits = Double.doubleToLongBits(y);
1440     int __HIy = (int) (yBits >> 32);
1441     if(k >= -1021) {
1442     __HIy += (k<<20); // add k to y's exponent
1443     yBits = ((__HIy & 0xFFFFFFFFL) << 32) | (yBits & 0xFFFFFFFFL);
1444     y = Double.longBitsToDouble(yBits);
1445     return y;
1446     } else {
1447     __HIy += ((k+1000)<<20);// add k to y's exponent
1448     yBits = ((__HIy & 0xFFFFFFFFL) << 32) | (yBits & 0xFFFFFFFFL);
1449     y = Double.longBitsToDouble(yBits);
1450     return y*twom1000;
1451     }
1452     }
1453     /**/
1454
1455}

A to Z: JavaDoc & Examples

---

Daily Java News & Articles

---

Open Source Projects

---

Open Source Codes

---

Free Computer Books

---

Remove Frame

---

Popular Tags