math.java
001 /*
002  * Java Genetic Algorithm Library (jenetics-2.0.2).
003  * Copyright (c) 2007-2014 Franz Wilhelmstötter
004  *
005  * Licensed under the Apache License, Version 2.0 (the "License");
006  * you may not use this file except in compliance with the License.
007  * You may obtain a copy of the License at
008  *
009  *      http://www.apache.org/licenses/LICENSE-2.0
010  *
011  * Unless required by applicable law or agreed to in writing, software
012  * distributed under the License is distributed on an "AS IS" BASIS,
013  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014  * See the License for the specific language governing permissions and
015  * limitations under the License.
016  *
017  * Author:
018  *    Franz Wilhelmstötter (franz.wilhelmstoetter@gmx.at)
019  */
020 package org.jenetics.util;
021 
022 import static java.lang.String.format;
023 import static java.util.Objects.requireNonNull;
024 
025 import java.util.Random;
026 
027 /**
028  * This object contains mathematical helper functions.
029  *
030  @author <a href="mailto:franz.wilhelmstoetter@gmx.at">Franz Wilhelmstötter</a>
031  @since 1.0
032  @version 1.4 &mdash; <em>$Date: 2014-04-03 $</em>
033  */
034 public final class math extends StaticObject {
035     private math() {}
036 
037     /**
038      * Add to long values and throws an ArithmeticException in the case of an
039      * overflow.
040      *
041      @param a the first summand.
042      @param b the second summand.
043      @return the sum of the given values.
044      @throws ArithmeticException if the summation would lead to an overflow.
045      */
046     public static long plus(final long a, final long b) {
047         final long z = a + b;
048         if (((a^z(b^z)) 0) {
049             throw new ArithmeticException(format("Overflow: %d + %d", a, b));
050         }
051 
052         return z;
053     }
054 
055     /**
056      * Subtracts to long values and throws an ArithmeticException in the case of
057      * an overflow.
058      *
059      @param a the minuend.
060      @param b the subtrahend.
061      @return the difference of the given values.
062      @throws ArithmeticException if the subtraction would lead to an overflow.
063      */
064     public static long minus(final long a, final long b) {
065         final long z = a - b;
066         if (((a^b(a^z)) 0) {
067             throw new ArithmeticException(format("Overflow: %d - %d", a, b));
068         }
069 
070         return z;
071     }
072 
073     /**
074      * Normalize the given double array, so that it sum to one. The
075      * normalization is performed in place and the same {@code values} are
076      * returned.
077      *
078      @param values the values to normalize.
079      @return the {@code values} array.
080      @throws NullPointerException if the given double array is {@code null}.
081      */
082     public static double[] normalize(final double[] values) {
083         final double sum = 1.0/statistics.sum(values);
084         for (int i = values.length; --i >= 0;) {
085             values[i= values[i]*sum;
086         }
087 
088         return values;
089     }
090 
091     /**
092      <i>Clamping</i> a value between a pair of boundary values.
093      <i>Note: using clamp with floating point numbers may give unexpected
094      * results if one of the values is {@code NaN}.</i>
095      *
096      @param v the value to <i>clamp</i>
097      @param lo the lower bound.
098      @param hi the upper bound.
099      @return The clamped value:
100      *        <ul>
101      *            <li>{@code lo if v < lo}</li>
102      *            <li>{@code hi if hi < v}</li>
103      *            <li>{@code otherwise, v}</li>
104      *        </ul>
105      */
106     public static double clamp(final double v, final double lo, final double hi) {
107         return v < lo ? lo : (v > hi ? hi : v);
108     }
109 
110     /**
111      * Component wise multiplication of the given double array.
112      *
113      @param values the double values to multiply.
114      @param multiplier the multiplier.
115      @throws NullPointerException if the given double array is {@code null}.
116      */
117     public static void times(final double[] values, final double multiplier) {
118         for (int i = values.length; --i >= 0;) {
119             values[i*= multiplier;
120         }
121     }
122 
123     /**
124      * Component wise division of the given double array.
125      *
126      @param values the double values to divide.
127      @param divisor the divisor.
128      @throws NullPointerException if the given double array is {@code null}.
129      */
130     public static void divide(final double[] values, final double divisor) {
131         for (int i = values.length; --i >= 0;) {
132             values[i/= divisor;
133         }
134     }
135 
136     /**
137      * Binary exponentiation algorithm.
138      *
139      @param b the base number.
140      @param e the exponent.
141      @return {@code b^e}.
142      */
143     public static long pow(final long b, final long e) {
144         long base = b;
145         long exp = e;
146         long result = 1;
147 
148         while (exp != 0) {
149             if ((exp & 1!= 0) {
150                 result *= base;
151             }
152             exp >>>= 1;
153             base *= base;
154         }
155 
156         return result;
157     }
158 
159     static boolean isMultiplicationSave(final int a, final int b) {
160         final long m = (long)a*(long)b;
161         return ((int)m== m;
162     }
163 
164     /**
165      * Return the <a href="http://en.wikipedia.org/wiki/Unit_in_the_last_place">ULP</a>
166      * distance of the given two double values.
167      *
168      @param a first double.
169      @param b second double.
170      @return the ULP distance.
171      @throws ArithmeticException if the distance doesn't fit in a long value.
172      */
173     public static long ulpDistance(final double a, final double b) {
174         return minus(ulpPosition(a), ulpPosition(b));
175     }
176 
177     /**
178      * Calculating the <a href="http://en.wikipedia.org/wiki/Unit_in_the_last_place">ULP</a>
179      * position of a double number.
180      *
181      * [code]
182      * double a = 0.0;
183      * for (int i = 0; i &lt; 10; ++i) {
184      *     a = Math.nextAfter(a, Double.POSITIVE_INFINITY);
185      * }
186      *
187      * for (int i = 0; i &lt; 19; ++i) {
188      *     a = Math.nextAfter(a, Double.NEGATIVE_INFINITY);
189      *     System.out.println(
190      *          a + "\t" + ulpPosition(a) + "\t" + ulpDistance(0.0, a)
191      *     );
192      * }
193      * [/code]
194      *
195      * The code fragment above will create the following output:
196      <pre>
197      *   4.4E-323    9  9
198      *   4.0E-323    8  8
199      *   3.5E-323    7  7
200      *   3.0E-323    6  6
201      *   2.5E-323    5  5
202      *   2.0E-323    4  4
203      *   1.5E-323    3  3
204      *   1.0E-323    2  2
205      *   4.9E-324    1  1
206      *   0.0         0  0
207      *  -4.9E-324   -1  1
208      *  -1.0E-323   -2  2
209      *  -1.5E-323   -3  3
210      *  -2.0E-323   -4  4
211      *  -2.5E-323   -5  5
212      *  -3.0E-323   -6  6
213      *  -3.5E-323   -7  7
214      *  -4.0E-323   -8  8
215      *  -4.4E-323   -9  9
216      </pre>
217      *
218      @param a the double number.
219      @return the ULP position.
220      */
221     public static long ulpPosition(final double a) {
222         long t = Double.doubleToLongBits(a);
223         if (t < 0) {
224             t = Long.MIN_VALUE - t;
225         }
226         return t;
227     }
228 
229     /**
230      * Selects a random subset of size {@code k} from a set of size {@code n}.
231      *
232      @see #subset(int, int[])
233      *
234      @param n the size of the set.
235      @param k the size of the subset.
236      @throws IllegalArgumentException if {@code n < k}, {@code k == 0} or if
237      *          {@code n*k} will cause an integer overflow.
238      @return the subset array.
239      */
240     public static int[] subset(final int n, final int k) {
241         return subset(n, k, RandomRegistry.getRandom());
242     }
243 
244     /**
245      * Selects a random subset of size {@code k} from a set of size {@code n}.
246      *
247      @see #subset(int, int[], Random)
248      *
249      @param n the size of the set.
250      @param k the size of the subset.
251      @param random the random number generator used.
252      @throws NullPointerException if {@code random} is {@code null}.
253      @throws IllegalArgumentException if {@code n < k}, {@code k == 0} or if
254      *         {@code n*k} will cause an integer overflow.
255      @return the subset array.
256      */
257     public static int[] subset(final int n, final int k, final Random random) {
258         requireNonNull(random, "Random");
259         if (k <= 0) {
260             throw new IllegalArgumentException(format(
261                     "Subset size smaller or equal zero: %s", k
262                 ));
263         }
264         if (n < k) {
265             throw new IllegalArgumentException(format(
266                     "n smaller than k: %s < %s.", n, k
267                 ));
268         }
269 
270         final int[] sub = new int[k];
271         subset(n, sub,random);
272         return sub;
273     }
274 
275     /**
276      <p>
277      * Selects a random subset of size {@code sub.length} from a set of size
278      * {@code n}.
279      </p>
280      *
281      <p>
282      <em>Authors:</em>
283      *      FORTRAN77 original version by Albert Nijenhuis, Herbert Wilf. This
284      *      version based on the  C++ version by John Burkardt.
285      </p>
286      *
287      <p><em><a href="https://people.scs.fsu.edu/~burkardt/c_src/subset/subset.html">
288      *  Reference:</a></em>
289      *      Albert Nijenhuis, Herbert Wilf,
290      *      Combinatorial Algorithms for Computers and Calculators,
291      *      Second Edition,
292      *      Academic Press, 1978,
293      *      ISBN: 0-12-519260-6,
294      *      LC: QA164.N54.
295      </p>
296      *
297      @param n the size of the set.
298      @param sub the sub set array.
299      @throws NullPointerException if {@code sub} is {@code null}.
300      @throws IllegalArgumentException if {@code n < sub.length},
301      *         {@code sub.length == 0} or {@code n*sub.length} will cause an
302      *         integer overflow.
303      */
304     public static void subset(final int n, final int sub[]) {
305         subset(n, sub, RandomRegistry.getRandom());
306     }
307 
308     /**
309      <p>
310      * Selects a random subset of size {@code sub.length} from a set of size
311      * {@code n}.
312      </p>
313      *
314      <p>
315      <em>Authors:</em>
316      *      FORTRAN77 original version by Albert Nijenhuis, Herbert Wilf. This
317      *      version based on the  C++ version by John Burkardt.
318      </p>
319      *
320      <p><em><a href="https://people.scs.fsu.edu/~burkardt/c_src/subset/subset.html">
321      *  Reference:</a></em>
322      *      Albert Nijenhuis, Herbert Wilf,
323      *      Combinatorial Algorithms for Computers and Calculators,
324      *      Second Edition,
325      *      Academic Press, 1978,
326      *      ISBN: 0-12-519260-6,
327      *      LC: QA164.N54.
328      </p>
329      *
330      @param n the size of the set.
331      @param sub the sub set array.
332      @param random the random number generator used.
333      @return the sub-set array for the given parameter
334      @throws NullPointerException if {@code sub} or {@code random} is
335      *         {@code null}.
336      @throws IllegalArgumentException if {@code n < sub.length},
337      *         {@code sub.length == 0} or {@code n*sub.length} will cause an
338      *         integer overflow.
339      */
340     public static int[] subset(final int n, final int sub[]final Random random) {
341         requireNonNull(random, "Random");
342         requireNonNull(sub, "Sub set array");
343 
344         final int k = sub.length;
345         if (k <= 0) {
346             throw new IllegalArgumentException(format(
347                 "Subset size smaller or equal zero: %s", k
348             ));
349         }
350         if (n < k) {
351             throw new IllegalArgumentException(format(
352                 "n smaller than k: %s < %s.", n, k
353             ));
354         }
355         if (!math.isMultiplicationSave(n, k)) {
356             throw new IllegalArgumentException(format(
357                 "n*sub.length > Integer.MAX_VALUE (%s*%s = %s > %s)",
358                 n, sub.length, (long)n*(long)k, Integer.MAX_VALUE
359             ));
360         }
361 
362         if (sub.length == n) {
363             for (int i = 0; i < sub.length; ++i) {
364                 sub[i= i;
365             }
366             return sub;
367         }
368 
369         for (int i = 0; i < k; ++i) {
370             sub[i(i*n)/k;
371         }
372 
373         int l = 0;
374         int ix = 0;
375         for (int i = 0; i < k; ++i) {
376             do {
377                 ix = nextInt(random, 1, n);
378                 l = (ix*k - 1)/n;
379             while (sub[l>= ix);
380 
381             sub[l= sub[l1;
382         }
383 
384         int m = 0;
385         int ip = 0;
386         int is = k;
387         for (int i = 0; i < k; ++i) {
388             m = sub[i];
389             sub[i0;
390 
391             if (m != (i*n)/k) {
392                 ip = ip + 1;
393                 sub[ip - 1= m;
394             }
395         }
396 
397         int ihi = ip;
398         int ids = 0;
399         for (int i = 1; i <= ihi; ++i) {
400             ip = ihi + - i;
401             l = (sub[ip - 1]*k - 1)/n;
402             ids = sub[ip - 1((l - 1)*n)/k;
403             sub[ip - 10;
404             sub[is - 1= l;
405             is = is - ids;
406         }
407 
408         int ir = 0;
409         int m0 = 0;
410         for (int ll = 1; ll <= k; ++ll) {
411             l = k + - ll;
412 
413             if (sub[l - 1!= 0) {
414                 ir = l;
415                 m0 = ((sub[l - 11)*n)/k;
416                 m = (sub[l-1]*n)/k - m0 + 1;
417             }
418 
419             ix = nextInt(random, m0, m0 + m - 1);
420 
421             int i = l + 1;
422             while (i <= ir && ix >= sub[i - 1]) {
423                 ix = ix + 1;
424                 subi- 2= sub[i - 1];
425                 i = i + 1;
426             }
427 
428             sub[i - 2= ix;
429             --m;
430         }
431 
432         return sub;
433     }
434 
435     private static int nextInt(final Random random, final int a, final int b) {
436         return a == b ? a - : random.nextInt(b - a+ a;
437     }
438 
439     /**
440      * Some helper method concerning statistics.
441      *
442      @author <a href="mailto:franz.wilhelmstoetter@gmx.at">Franz Wilhelmstötter</a>
443      @since 1.3
444      @version 1.3 &mdash; <em>$Date: 2014-04-03 $</em>
445      */
446     public static final class statistics extends StaticObject {
447         private statistics() {}
448 
449         /**
450          * Return the minimum value of the given double array.
451          *
452          @param values the double array.
453          @return the minimum value or {@link Double#NaN} if the given array is
454          *         empty.
455          @throws NullPointerException if the given array is {@code null}.
456          */
457         public static double min(final double[] values) {
458             double min = Double.NaN;
459             if (values.length > 0) {
460                 min = values[0];
461 
462                 for (int i = values.length; --i >= 1;) {
463                     if (values[i< min) {
464                         min = values[i];
465                     }
466                 }
467             }
468 
469             return min;
470         }
471 
472         /**
473          * Return the maximum value of the given double array.
474          *
475          @param values the double array.
476          @return the maximum value or {@link Double#NaN} if the given array is
477          *         empty.
478          @throws NullPointerException if the given array is {@code null}.
479          */
480         public static double max(final double[] values) {
481             double max = Double.NaN;
482             if (values.length > 0) {
483                 max = values[0];
484 
485                 for (int i = values.length; --i >= 1;) {
486                     if (values[i> max) {
487                         max = values[i];
488                     }
489                 }
490             }
491 
492             return max;
493         }
494 
495         /**
496          * Implementation of the <a href="http://en.wikipedia.org/wiki/Kahan_summation_algorithm">
497          * Kahan summation algorithm</a>.
498          *
499          @param values the values to sum up.
500          @return the sum of the given {@code values}.
501          @throws NullPointerException if the given array is {@code null}.
502          */
503         public static double sum(final double[] values) {
504             double sum = 0.0;
505             double c = 0.0;
506             double y = 0.0;
507             double t = 0.0;
508 
509             for (int i = values.length; --i >= 0;) {
510                 y = values[i- c;
511                 t = sum + y;
512                 c = t - sum - y;
513                 sum = t;
514             }
515 
516             return sum;
517         }
518 
519         /**
520          * Add the values of the given array.
521          *
522          @param values the values to add.
523          @return the values sum.
524          @throws NullPointerException if the values are null;
525          */
526         public static long sum(final long[] values) {
527             long sum = 0;
528             for (int i = values.length; --i >= 0;) {
529                 sum += values[i];
530             }
531             return sum;
532         }
533 
534     }
535 
536     /**
537      * Some helper method concerning random numbers and random seed generation.
538      *
539      @author <a href="mailto:franz.wilhelmstoetter@gmx.at">Franz Wilhelmstötter</a>
540      @since 1.1
541      @version 1.2 &mdash; <em>$Date: 2014-04-03 $</em>
542      */
543     public static final class random extends StaticObject {
544         private random() {}
545 
546         /**
547          * Returns a pseudo-random, uniformly distributed int value between min
548          * and max (min and max included).
549          *
550          @param random the random engine to use for calculating the random
551          *        int value
552          @param min lower bound for generated integer
553          @param max upper bound for generated integer
554          @return a random integer greater than or equal to {@code min} and
555          *         less than or equal to {@code max}
556          @throws IllegalArgumentException if {@code min >= max}
557          */
558         public static int nextInt(
559             final Random random,
560             final int min, final int max
561         ) {
562             if (min >= max) {
563                 throw new IllegalArgumentException(format(
564                     "Min >= max: %d >= %d", min, max
565                 ));
566             }
567 
568             final int diff = max - min + 1;
569             int result = 0;
570 
571             if (diff <= 0) {
572                 do {
573                     result = random.nextInt();
574                 while (result < min || result > max);
575             else {
576                 result = random.nextInt(diff+ min;
577             }
578 
579             return result;
580         }
581 
582         /**
583          * Returns a pseudo-random, uniformly distributed int value between min
584          * and max (min and max included).
585          *
586          @param random the random engine to use for calculating the random
587          *        long value
588          @param min lower bound for generated long integer
589          @param max upper bound for generated long integer
590          @return a random long integer greater than or equal to {@code min}
591          *         and less than or equal to {@code max}
592          @throws IllegalArgumentException if {@code min >= max}
593          */
594         public static long nextLong(
595             final Random random,
596             final long min, final long max
597         ) {
598             if (min >= max) {
599                 throw new IllegalArgumentException(format(
600                     "min >= max: %d >= %d.", min, max
601                 ));
602             }
603 
604             final long diff = (max - min1;
605             long result = 0;
606 
607             if (diff <= 0) {
608                 do {
609                     result = random.nextLong();
610                 while (result < min || result > max);
611             else if (diff < Integer.MAX_VALUE) {
612                 result = random.nextInt((int)diff+ min;
613             else {
614                 result = nextLong(random, diff+ min;
615             }
616 
617             return result;
618         }
619 
620         /**
621          * Returns a pseudo-random, uniformly distributed int value between 0
622          * (inclusive) and the specified value (exclusive), drawn from the given
623          * random number generator's sequence.
624          *
625          @param random the random engine used for creating the random number.
626          @param n the bound on the random number to be returned. Must be
627          *        positive.
628          @return the next pseudo-random, uniformly distributed int value
629          *         between 0 (inclusive) and n (exclusive) from the given random
630          *         number generator's sequence
631          @throws IllegalArgumentException if n is smaller than 1.
632          */
633         public static long nextLong(final Random random, final long n) {
634             if (n <= 0) {
635                 throw new IllegalArgumentException(format(
636                     "n is smaller than one: %d", n
637                 ));
638             }
639 
640             long bits;
641             long result;
642             do {
643                 bits = random.nextLong() 0x7fffffffffffffffL;
644                 result = bits%n;
645             while (bits - result + (n - 10);
646 
647             return result;
648         }
649 
650         /**
651          * Returns a pseudo-random, uniformly distributed double value between
652          * min (inclusively) and max (exclusively).
653          *
654          @param random the random engine used for creating the random number.
655          @param min lower bound for generated float value
656          @param max upper bound for generated float value
657          @return a random float greater than or equal to {@code min} and less
658          *         than to {@code max}
659          */
660         public static float nextFloat(
661             final Random random,
662             final float min, final float max
663         ) {
664             return random.nextFloat()*(max - min+ min;
665         }
666 
667         /**
668          * Returns a pseudo-random, uniformly distributed double value between
669          * min (inclusively) and max (exclusively).
670          *
671          @param random the random engine used for creating the random number.
672          @param min lower bound for generated double value
673          @param max upper bound for generated double value
674          @return a random double greater than or equal to {@code min} and less
675          *         than to {@code max}
676          */
677         public static double nextDouble(
678             final Random random,
679             final double min, final double max
680         ) {
681             return random.nextDouble()*(max - min+ min;
682         }
683 
684         /**
685          * Create a new <em>seed</em> byte array of the given length.
686          *
687          @see #seed(byte[])
688          @see #seed()
689          *
690          @param length the length of the returned byte array.
691          @return a new <em>seed</em> byte array of the given length
692          @throws NegativeArraySizeException if the given length is smaller
693          *         than zero.
694          */
695         public static byte[] seedBytes(final int length) {
696             return seed(new byte[length]);
697         }
698 
699         /**
700          * Fills the given byte array with random bytes, created by successive
701          * calls of the {@link #seed()} method.
702          *
703          @see #seed()
704          *
705          @param seed the byte array seed to fill with random bytes.
706          @return the given byte array, for method chaining.
707          @throws NullPointerException if the {@code seed} array is
708          *         {@code null}.
709          */
710         public static byte[] seed(final byte[] seed) {
711             for (int i = 0, len = seed.length; i < len;) {
712                 int n = Math.min(len - i, Long.SIZE/Byte.SIZE);
713 
714                 for (long x = seed(); n-- > 0; x >>= Byte.SIZE) {
715                     seed[i++(byte)x;
716                 }
717             }
718 
719             return seed;
720         }
721 
722         /**
723          * Calculating a 64 bit seed value which can be used for initializing
724          * PRNGs. This method uses a combination of {@code System.nanoTime()}
725          * and {@code new Object().hashCode()} calls to create a reasonable safe
726          * seed value:
727          <p>
728          * [code]
729          * public static long seed() {
730          *     return seed(System.nanoTime());
731          * }
732          * [/code]
733          <p>
734          * This method passes all of the statistical tests of the
735          * <a href="http://www.phy.duke.edu/~rgb/General/dieharder.php">
736          * dieharder</a> test suite&mdash;executed on a linux machine with
737          * JDK version 1.7. <em>Since there is no prove that this will the case
738          * for every Java version and OS, it is recommended to only use this
739          * method for seeding other PRNGs.</em>
740          *
741          @see #seed(long)
742          *
743          @return the random seed value.
744          */
745         public static long seed() {
746             return seed(System.nanoTime());
747         }
748 
749         /**
750          * Uses the given {@code base} value to create a reasonable safe seed
751          * value. This is done by combining it with values of
752          * {@code new Object().hashCode()}:
753          <p>
754          * [code]
755          * public static long seed(final long base) {
756          *     final long objectHashSeed = ((long)(new Object().hashCode()) &lt;&lt; 32) |
757          *                                         new Object().hashCode();
758          *     long seed = base ^ objectHashSeed;
759          *     seed ^= seed &lt;&lt; 17;
760          *     seed ^= seed &gt;&gt;&gt; 31;
761          *     seed ^= seed &lt;&lt; 8;
762          *     return seed;
763          * }
764          * [/code]
765          *
766          @param base the base value of the seed to create
767          @return the created seed value.
768          */
769         public static long seed(final long base) {
770             return mix(base, objectHashSeed());
771         }
772 
773         private static long mix(final long a, final long b) {
774             long c = a^b;
775             c ^= c << 17;
776             c ^= c >>> 31;
777             c ^= c << 8;
778             return c;
779         }
780 
781         private static long objectHashSeed() {
782             return ((long)(new Object().hashCode()) << 32|
783                             new Object().hashCode();
784         }
785 
786     }
787 
788 }