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 */
020package org.jenetics.util;
021
022import static java.lang.String.format;
023import static java.util.Objects.requireNonNull;
024
025import 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 — <em>$Date: 2014-04-03 $</em>
033 */
034public 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 * "line.181">
class="sourceLineNo">182 * double a = 0.0;"line.182">
class="sourceLineNo">183 * for (int i = 0; i < 10; ++i) {"line.183">
class="sourceLineNo">184 * a = Math.nextAfter(a, Double.POSITIVE_INFINITY);"line.184">
class="sourceLineNo">185 * }"line.185">
class="sourceLineNo">186 *"line.186">
class="sourceLineNo">187 * for (int i = 0; i < 19; ++i) {"line.187">
class="sourceLineNo">188 * a = Math.nextAfter(a, Double.NEGATIVE_INFINITY);"line.188">
class="sourceLineNo">189 * System.out.println("line.189">
class="sourceLineNo">190 * a + "\t" + ulpPosition(a) + "\t" + ulpDistance(0.0, a)"line.190">
class="sourceLineNo">191 * );"line.191">
class="sourceLineNo">192 * }"line.192">
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[l] + 1;
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[i] = 0;
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 + 1 - i;
401 l = 1 + (sub[ip - 1]*k - 1)/n;
402 ids = sub[ip - 1] - ((l - 1)*n)/k;
403 sub[ip - 1] = 0;
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 + 1 - ll;
412
413 if (sub[l - 1] != 0) {
414 ir = l;
415 m0 = 1 + ((sub[l - 1] - 1)*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 sub[ i- 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 - 1 : 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 — <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 — <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 - min) + 1;
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 - 1) < 0);
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 *
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—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 * "line.754">
class="sourceLineNo">755 * public static long seed(final long base) {"line.755">
class="sourceLineNo">756 * final long objectHashSeed = ((long)(new Object().hashCode()) << 32) |"line.756">
class="sourceLineNo">757 * new Object().hashCode();"line.757">
class="sourceLineNo">758 * long seed = base ^ objectHashSeed;"line.758">
class="sourceLineNo">759 * seed ^= seed << 17;"line.759">
class="sourceLineNo">760 * seed ^= seed >>> 31;"line.760">
class="sourceLineNo">761 * seed ^= seed << 8;"line.761">
class="sourceLineNo">762 * return seed;"line.762">
class="sourceLineNo">763 * }"line.763">
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}