001// ============================================================================
002//   Copyright 2006-2012 Daniel W. Dyer
003//
004//   Licensed under the Apache License, Version 2.0 (the "License");
005//   you may not use this file except in compliance with the License.
006//   You may obtain a copy of the License at
007//
008//       http://www.apache.org/licenses/LICENSE-2.0
009//
010//   Unless required by applicable law or agreed to in writing, software
011//   distributed under the License is distributed on an "AS IS" BASIS,
012//   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
013//   See the License for the specific language governing permissions and
014//   limitations under the License.
015// ============================================================================
016package squidpony.squidmath;
017
018import squidpony.annotation.GwtIncompatible;
019
020import java.math.BigInteger;
021
022/**
023 * Mathematical operations not provided by {@link Math java.lang.Math}.
024 * <br>
025 * Originally part of the <a href="http://maths.uncommons.org/">Uncommon Maths software package</a> as Maths.
026 * @author Daniel Dyer
027 */
028public final class MathExtras
029{
030    // The biggest factorial that can be calculated using 64-bit signed longs.
031    private static final int MAX_LONG_FACTORIAL = 20;
032
033    // Cache BigInteger factorial values because they are expensive to generate.
034    private static final int CACHE_SIZE = 256;
035    static final long[] factorialsStart = new long[]{
036            0, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800, 479001600,
037            6227020800L, 87178291200L, 1307674368000L, 20922789888000L, 355687428096000L,
038            6402373705728000L, 121645100408832000L, 2432902008176640000L
039    };
040    private static final OrderedMap<Integer, BigInteger> BIG_FACTORIALS
041        = new OrderedMap<Integer, BigInteger>();
042
043
044    private MathExtras()
045    {
046        // Prevent instantiation.
047    }
048
049
050    /**
051     * Calculates the factorial of n where n is a number in the
052     * range 0 - 20.  Zero factorial is equal to 1.  For values of
053     * n greater than 20 you must use {@link #bigFactorial(int)}.
054     * @param n The factorial to calculate.
055     * @return The factorial of n.
056     * @see #bigFactorial(int)
057     */
058    public static long factorial(int n)
059    {
060        if (n < 0 || n > MAX_LONG_FACTORIAL)
061        {
062            throw new IllegalArgumentException("Argument must be in the range 0 - 20.");
063        }
064        /*
065        long factorial = 1;
066        for (int i = n; i > 1; i--)
067        {
068            factorial *= i;
069        }*/
070        return factorialsStart[n];
071    }
072
073
074    /**
075     * Calculates the factorial of n where n is a positive integer.
076     * Zero factorial is equal to 1.  For values of n up to 20, consider
077     * using {@link #factorial(int)} instead since it uses a faster
078     * implementation.
079     * @param n The factorial to calculate.
080     * @return The factorial of n.
081     * @see #factorial(int)
082     */
083    public static BigInteger bigFactorial(int n)
084    {
085        if (n < 0)
086        {
087            throw new IllegalArgumentException("Argument must greater than or equal to zero.");
088        }
089
090        BigInteger factorial = null;
091        if (n < CACHE_SIZE) // Check for a cached value.
092        {
093            factorial = BIG_FACTORIALS.get(n);
094        }
095
096        if (factorial == null)
097        {
098            factorial = BigInteger.ONE;
099            for (int i = n; i > 1; i--)
100            {
101                factorial = factorial.multiply(BigInteger.valueOf(i));
102            }
103            if (n < CACHE_SIZE) // Cache value.
104            {
105                if(!BIG_FACTORIALS.containsKey(n))
106                    BIG_FACTORIALS.put(n, factorial);
107            }
108        }
109
110        return factorial;
111    }
112
113
114    /**
115     * Calculate the first argument raised to the power of the second.
116     * This method only supports non-negative powers.
117     * @param value The number to be raised.
118     * @param power The exponent (must be positive).
119     * @return {@code value} raised to {@code power}.
120     */
121    public static long raiseToPower(int value, int power)
122    {
123        if (power < 0)
124        {
125            throw new IllegalArgumentException("This method does not support negative powers.");
126        }
127        long result = 1;
128        for (int i = 0; i < power; i++)
129        {
130            result *= value;
131        }
132        return result;
133    }
134
135
136    /**
137     * Calculate logarithms for arbitrary bases.
138     * @param base The base for the logarithm.
139     * @param arg The value to calculate the logarithm for.
140     * @return The log of {@code arg} in the specified {@code base}.
141     */
142    public static double log(double base, double arg)
143    {
144        // Use natural logarithms and change the base.
145        return Math.log(arg) / Math.log(base);
146    }
147
148
149    /**
150     * Checks that two values are approximately equal (plus or minus a specified tolerance).
151     * @param value1 The first value to compare.
152     * @param value2 The second value to compare.
153     * @param tolerance How much (in percentage terms, as a percentage of the first value)
154     * the values are allowed to differ and still be considered equal.  Expressed as a value
155     * between 0 and 1.
156     * @return true if the values are approximately equal, false otherwise.
157     */
158    public static boolean approxEquals(double value1,
159                                       double value2,
160                                       double tolerance)
161    {
162        if (tolerance < 0 || tolerance > 1)
163        {
164            throw new IllegalArgumentException("Tolerance must be between 0 and 1.");
165        }
166        return Math.abs(value1 - value2) <= value1 * tolerance;
167    }
168
169
170    /**
171     * If the specified value is not greater than or equal to the specified minimum and
172     * less than or equal to the specified maximum, adjust it so that it is.
173     * @param value The value to check.
174     * @param min The minimum permitted value.
175     * @param max The maximum permitted value.
176     * @return {@code value} if it is between the specified limits, {@code min} if the value
177     * is too low, or {@code max} if the value is too high.
178     */
179    public static int clamp(int value, int min, int max)
180    {
181        return Math.min(Math.max(value, min), max);
182    }
183
184
185    /**
186     * If the specified value is not greater than or equal to the specified minimum and
187     * less than or equal to the specified maximum, adjust it so that it is.
188     * @param value The value to check.
189     * @param min The minimum permitted value.
190     * @param max The maximum permitted value.
191     * @return {@code value} if it is between the specified limits, {@code min} if the value
192     * is too low, or {@code max} if the value is too high.
193     */
194    public static long clamp(long value, long min, long max)
195    {
196        return Math.min(Math.max(value, min), max);
197    }
198
199
200    /**
201     * If the specified value is not greater than or equal to the specified minimum and
202     * less than or equal to the specified maximum, adjust it so that it is.
203     * @param value The value to check.
204     * @param min The minimum permitted value.
205     * @param max The maximum permitted value.
206     * @return {@code value} if it is between the specified limits, {@code min} if the value
207     * is too low, or {@code max} if the value is too high.
208     */
209    public static double clamp(double value, double min, double max)
210    {
211        return Math.min(Math.max(value, min), max);
212    }
213
214    /**
215     * If the specified value is not greater than or equal to the specified minimum and
216     * less than or equal to the specified maximum, adjust it so that it is.
217     * @param value The value to check.
218     * @param min The minimum permitted value.
219     * @param max The maximum permitted value.
220     * @return {@code value} if it is between the specified limits, {@code min} if the value
221     * is too low, or {@code max} if the value is too high.
222     */
223    public static float clamp(float value, float min, float max)
224    {
225        return Math.min(Math.max(value, min), max);
226    }
227
228    /**
229     * Like the modulo operator {@code %}, but the result will always match the sign of {@code d} instead of {@code op}.
230     * @param op the dividend; negative values are permitted and wrap instead of producing negative results
231     * @param d the divisor; if this is negative then the result will be negative, otherwise it will be positive
232     * @return the remainder of the division of op by d, with a sign matching d
233     */
234    public static double remainder(final double op, final double d) {
235        return (op % d + d) % d;
236    }
237
238    /**
239     * Determines the greatest common divisor of a pair of natural numbers
240     * using the Euclidean algorithm.  This method only works with natural
241     * numbers.  If negative integers are passed in, the absolute values will
242     * be used.  The return value is always positive.
243     * @param a The first value.
244     * @param b The second value.
245     * @return The greatest common divisor.
246     */
247    public static long greatestCommonDivisor(long a, long b)
248    {
249        a = Math.abs(a);
250        b = Math.abs(b);
251        while (b != 0)
252        {
253            long temp = b;
254            b = a % b;
255            a = temp;
256        }
257        return a;
258    }
259
260    /**
261     * Given any odd int {@code a}, this finds another odd int {@code b} such that {@code a * b == 1}.
262     * <br>
263     * This is incompatible with GWT, but it should usually only find uses in exploratory code or in tests anyway...
264     * It is only incompatible because it tends to rely on multiplication overflow to work.
265     * @param a any odd int; note that even numbers do not have inverses modulo 2 to the 32
266     * @return the multiplicative inverse of {@code a} modulo 4294967296 (or, 2 to the 32)
267     */
268    @GwtIncompatible
269    public static int modularMultiplicativeInverse(int a)
270    {
271        int x = 2 ^ a * 3;     //  5 bits
272        x *= 2 - a * x;        // 10
273        x *= 2 - a * x;        // 20
274        x *= 2 - a * x;        // 40 -- 32 low bits
275        return x;
276    }
277
278    /**
279     * Given any odd long {@code a}, this finds another odd long {@code b} such that {@code a * b == 1L}.
280     * @param a any odd long; note that even numbers do not have inverses modulo 2 to the 64
281     * @return the multiplicative inverse of {@code a} modulo 18446744073709551616 (or, 2 to the 64)
282     */
283    public static long modularMultiplicativeInverse(long a)
284    {
285        long x = 2 ^ a * 3;    //  5 bits
286        x *= 2 - a * x;        // 10
287        x *= 2 - a * x;        // 20
288        x *= 2 - a * x;        // 40
289        x *= 2 - a * x;        // 80 -- 64 low bits
290        return x;
291    }
292}