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 &mdash; <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 &lt; 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 &lt; 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 &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 - 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 *
"line.728"> class="sourceLineNo">729 * public static long seed() {"line.729"> class="sourceLineNo">730 * return seed(System.nanoTime());"line.730"> class="sourceLineNo">731 * }"line.731">
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 *
"line.754"> class="sourceLineNo">755 * public static long seed(final long base) {"line.755"> class="sourceLineNo">756 * final long objectHashSeed = ((long)(new Object().hashCode()) &lt;&lt; 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 &lt;&lt; 17;"line.759"> class="sourceLineNo">760 * seed ^= seed &gt;&gt;&gt; 31;"line.760"> class="sourceLineNo">761 * seed ^= seed &lt;&lt; 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}