001/**
002 * Copyright (c) 2011, The University of Southampton and the individual contributors.
003 * All rights reserved.
004 *
005 * Redistribution and use in source and binary forms, with or without modification,
006 * are permitted provided that the following conditions are met:
007 *
008 *   *  Redistributions of source code must retain the above copyright notice,
009 *      this list of conditions and the following disclaimer.
010 *
011 *   *  Redistributions in binary form must reproduce the above copyright notice,
012 *      this list of conditions and the following disclaimer in the documentation
013 *      and/or other materials provided with the distribution.
014 *
015 *   *  Neither the name of the University of Southampton nor the names of its
016 *      contributors may be used to endorse or promote products derived from this
017 *      software without specific prior written permission.
018 *
019 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
020 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
021 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
022 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
023 * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
024 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
025 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
026 * ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
027 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
028 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
029 */
030package org.openimaj.math.util;
031
032
033/**
034 * A collection of maths functions not available anywhere else ive seen
035 *
036 * @author Sina Samangooei (ss@ecs.soton.ac.uk)
037 *
038 */
039public class MathUtils {
040
041        /**
042         * Given log(a) and log(b) calculate log(a + b) boils down to log(
043         * exp(log_a) + exp(log_b) ) but this might overflow, so we turn this into
044         * log([exp(log_a - log_c) + exp(log_b - log_c)]exp(log_c)) and we set log_c
045         * == max(log_a,log_b) and so it becomes: LARGE + log(1 + exp(SMALL -
046         * LARGE)) == LARGE + log(1 + SMALL) ~= large the whole idea being to avoid
047         * an overflow (exp(LARGE) == VERY LARGE == overflow)
048         *
049         * @param log_a
050         * @param log_b
051         * @return log(a+b)
052         */
053        public static double logSum(final double log_a, final double log_b) {
054                double v;
055
056                if (log_a < log_b) {
057                        v = log_b + Math.log(1 + Math.exp(log_a - log_b));
058                } else {
059                        v = log_a + Math.log(1 + Math.exp(log_b - log_a));
060                }
061                return (v);
062        }
063
064        /**
065         * Returns the next power of 2 superior to n.
066         *
067         * @param n
068         *            The value to find the next power of 2 above
069         * @return The next power of 2
070         */
071        public static int nextPowerOf2(final int n) {
072                return (int) Math.pow(2, 32 - Integer.numberOfLeadingZeros(n - 1));
073        }
074
075        /**
076         * Implementation of the C <code>frexp</code> function to break
077         * floating-point number into normalized fraction and power of 2.
078         *
079         * @see "http://stackoverflow.com/questions/1552738/is-there-a-java-equivalent-of-frexp"
080         *
081         * @param value
082         *            the value
083         * @return the exponent and mantissa of the input value
084         */
085        public static ExponentAndMantissa frexp(double value) {
086                final ExponentAndMantissa ret = new ExponentAndMantissa();
087
088                ret.exponent = 0;
089                ret.mantissa = 0;
090
091                if (value == 0.0 || value == -0.0) {
092                        return ret;
093                }
094
095                if (Double.isNaN(value)) {
096                        ret.mantissa = Double.NaN;
097                        ret.exponent = -1;
098                        return ret;
099                }
100
101                if (Double.isInfinite(value)) {
102                        ret.mantissa = value;
103                        ret.exponent = -1;
104                        return ret;
105                }
106
107                ret.mantissa = value;
108                ret.exponent = 0;
109                int sign = 1;
110
111                if (ret.mantissa < 0f) {
112                        sign--;
113                        ret.mantissa = -(ret.mantissa);
114                }
115                while (ret.mantissa < 0.5f) {
116                        ret.mantissa *= 2.0f;
117                        ret.exponent -= 1;
118                }
119                while (ret.mantissa >= 1.0f) {
120                        ret.mantissa *= 0.5f;
121                        ret.exponent++;
122                }
123                ret.mantissa *= sign;
124                return ret;
125        }
126
127        /**
128         * Class to hold an exponent and mantissa
129         *
130         * @author Jonathon Hare (jsh2@ecs.soton.ac.uk)
131         */
132        public static class ExponentAndMantissa {
133                /**
134                 * The exponent
135                 */
136                public int exponent;
137
138                /**
139                 * The mantissa
140                 */
141                public double mantissa;
142        }
143}