001 package com.softnetConsult.utils.math;
002
003 import java.math.BigDecimal;
004 import java.math.BigInteger;
005
006
007 /**
008 * This class is a collection of static mathematical utility methods and functions.
009 *
010 * <p style="font-size:smaller;">This product includes software developed by the
011 * <strong>SoftNet-Consult Java Utility Library</strong> project and its contributors.<br />
012 * (<a href="http://java-tools.sourceforge.net" target="_blank">http://java-tools.sourceforge.net</a>)<br />
013 * Copyright (c) 2007-2008 SoftNet-Consult.<br />
014 * Copyright (c) 2007-2008 G. Paperin.<br />
015 * All rights reserved.
016 * </p>
017 * <p style="font-size:smaller;">File: MathTools.java<br />
018 * Library API version: {@value com.softnetConsult.utils.APIProperties#apiVersion}<br />
019 * Java compliance version: {@value com.softnetConsult.utils.APIProperties#javaComplianceVersion}
020 * </p>
021 * <p style="font-size:smaller;">Redistribution and use in source and binary forms, with or
022 * without modification, are permitted provided that the following terms and conditions are met:
023 * </p>
024 * <p style="font-size:smaller;">1. Redistributions of source code must retain the above
025 * acknowledgement of the SoftNet-Consult Java Utility Library project, the above copyright
026 * notice, this list of conditions and the following disclaimer.<br />
027 * 2. Redistributions in binary form must reproduce the above acknowledgement of the
028 * SoftNet-Consult Java Utility Library project, the above copyright notice, this list of
029 * conditions and the following disclaimer in the documentation and/or other materials
030 * provided with the distribution.<br />
031 * 3. All advertising materials mentioning features or use of this software or any derived
032 * software must display the following acknowledgement:<br />
033 * <em>This product includes software developed by the SoftNet-Consult Java Utility Library
034 * project and its contributors.</em>
035 * </p>
036 * <p style="font-size:smaller;">THIS SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY
037 * OF ANY KIND, EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
038 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
039 * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
040 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR
041 * IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
042 * </p>
043 * @author Greg Paperin (<a href="http://www.paperin.org" target="_blank">http://www.paperin.org</a>)
044 * @version {@value com.softnetConsult.utils.APIProperties#apiVersion}
045 *
046 */
047 public final class MathTools {
048
049 /**
050 * Prevents instances of this class from being created
051 * as this class contains only static utility methods.
052 */
053 private MathTools() {}
054
055 /**
056 * Powers of 10 as integers for fast fetching.
057 */
058 private static final int[] pow10 = new int[] {1, // 10^0
059 10, // 10^1
060 100, // 10^2
061 1000, // 10^3
062 10000, // 10^4
063 100000, // 10^5
064 1000000, // 10^6
065 10000000, // 10^7
066 100000000, // 10^8
067 1000000000}; // 10^9
068
069 /**
070 * Powers of 10 as longs for fast fetchng.
071 */
072 private static final long[] pow10L = new long[] {1L, // 10^0
073 10L, // 10^1
074 100L, // 10^2
075 1000L, // 10^3
076 10000L, // 10^4
077 100000L, // 10^5
078 1000000L, // 10^6
079 10000000L, // 10^7
080 100000000L, // 10^8
081 1000000000L, // 10^9
082 10000000000L, // 10^10
083 100000000000L, // 10^11
084 1000000000000L, // 10^12
085 10000000000000L, // 10^13
086 100000000000000L, // 10^14
087 1000000000000000L, // 10^15
088 10000000000000000L, // 10^16
089 100000000000000000L, // 10^17
090 1000000000000000000L}; // 10^18
091
092
093 /**
094 * Lookup tables of factorials that can be represented as an int.
095 */
096 private static final int[] factorialsInt = new int[] {1, 1, 2, 6, 24, 120, 720, 5040, 40320,
097 362880, 3628800, 39916800, 479001600};
098
099 /**
100 * Lookup tables of factorials that can be represented as a long.
101 */
102 private static final long[] factorialsLong = new long[] {1L, 1L, 2L, 6L, 24L, 120L, 720L, 5040L, 40320L,
103 362880L, 3628800L, 39916800L, 479001600L,
104 6227020800L, 87178291200L, 1307674368000L,
105 20922789888000L, 355687428096000L, 6402373705728000L,
106 121645100408832000L, 2432902008176640000L};
107
108 /**
109 * Used for fast calculation of the gamma function.
110 */
111 private static final double[] gammaCoefficients = new double[] {0.99999999999980993, 676.5203681218851,
112 -1259.1392167224028, 771.32342877765313, -176.61502916214059, 12.507343278686905,
113 -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7};
114
115 /**
116 * Constant used for the calculation of the gamma function.
117 */
118 private static final double sqrt2Pi = Math.sqrt(2. * Math.PI);
119
120 /**
121 * Constant used for the canculation of the natural logarithm of the gamma function.
122 */
123 private static final double ln2Pi = Math.log(2. * Math.PI);
124
125
126 /**
127 * Computes an approximation to the gamma function.
128 * In most cases, this approximation is correct for at least 15 significant digits.
129 *
130 * @param x The real parameter to the gamma function.
131 * @return The approximate value of the gamma function at {@code x}.
132 */
133 public static double gamma(final double x) {
134
135 if (x < 0.5)
136 return Math.PI / (Math.sin(Math.PI * x) * gamma(1. - x));
137
138 final double c1 = x - 1;
139 double s = gammaCoefficients[0];
140 for (int i = 1; i <= 8; i++)
141 s += gammaCoefficients[i] / (c1 + (double) i);
142 final double c2 = x + 6.5;
143 return sqrt2Pi * Math.pow(c2, c1 + 0.5) * Math.exp(-c2) * s;
144 }
145
146
147 /**
148 * Computes an approximation to the natiral logarithm of the gamma function.
149 * In most cases, this approximation is correct for at least 8 significant digits.
150 * This approximation is designed for positive arguments only; the return value for
151 * other agruments is unspecified and is most likely to be {@code Double.NaN}.
152 *
153 * @param x A positive real number.
154 * @return The approximate value of {@code ln(gamma(x))}.
155 */
156 public static double lnGamma(final double x) {
157 double x6 = x * x * x; x6 *= x6;
158 final double lnX = Math.log(x);
159 double _2lnG = ln2Pi - lnX
160 + x * (2. * lnX + Math.log(x * Math.sinh(1./x) + 1./(810. * x6)) - 2.);
161 return 0.5 * _2lnG;
162 }
163
164
165 /**
166 * Computes the factorial of a specified number.
167 *
168 * @param x A number between 0 and 12 inclusive (the factorial of larger numbers cannot be represented as int).
169 * @return {@code x!}.
170 */
171 public static int factorial(final int x) {
172
173 if (0 > x)
174 throw new IllegalArgumentException("Cannot compute an int-type factorial of a negative number");
175 if (12 < x)
176 throw new IllegalArgumentException("Cannot compute an int-type factorial of a number larger"
177 + " than 12 due to overflow.");
178
179 return factorialsInt[x];
180 }
181
182 /**
183 * Computes the factorial of a specified number.
184 *
185 * @param x A number between 0 and 20 inclusive (the factorial of larger numbers cannot be represented as long).
186 * @return {@code x!}.
187 */
188 public static long factorial(final long x) {
189
190 if (0 > x)
191 throw new IllegalArgumentException("Cannot compute a long-type factorial of a negative number");
192 if (20 < x)
193 throw new IllegalArgumentException("Cannot compute a long-type factorial of a number larger"
194 + " than 20 due to overflow.");
195
196 return factorialsLong[(int) x];
197 }
198
199 /**
200 * Computes an approximation to the factorial of a specified number
201 * using the fact that {@code x! = gamma(x + 1)}.
202 *
203 * @param x A number.
204 * @return {@code x!} calculated as {@code MathTools.gamma(x + 1.0)}.
205 * @see #gamma(double)
206 */
207 public static double factorial(final double x) {
208 return gamma(x + 1.0);
209 }
210
211 /**
212 * Computes an approximation to the natural logarithm of the factorial of a specified number.
213 * For values larger than {@code 1.0} this method uses a direct approximation to the natural
214 * logarithm of the gamma function (i.e. {@code ln(x!) = MathTools.lnGamma(x + 1.0))};
215 * for {@code 1.0} and smaller values this method explicitly
216 * computes the factorial using {@code MathTools.factorial(x)} and then takes the natural
217 * logarithm using {@code java.lang.Math.log()}.
218 *
219 * @param x A number larger than {@code -1}.
220 * @return THe approximate value of {@code ln(x!)}.
221 * @see #lnGamma(double)
222 * @see #factorial(double)
223 */
224 public static double lnFactorial(final double x) {
225 if (x <= 1.)
226 return Math.log(factorial(x));
227 return lnGamma(x + 1.0);
228 }
229
230
231 /**
232 * Probability mass function of the Poisson distribution.
233 *
234 * @param lambda LAMBDA-Poisson parameter.
235 * @param k K-Poisson parameter.
236 * @return Approximate value of the probability mass function of the Poisson distribution at given parameters.
237 */
238 public static double poisson(final double lambda, final int k) {
239
240 if (0 > lambda)
241 throw new IllegalArgumentException("Poisson parameter may not be negative.");
242
243 if (0 > k)
244 throw new IllegalArgumentException("Poisson parameter may not be negative.");
245
246 final double _k = (double) k;
247 final double nom = Math.exp(- lambda) * Math.pow(lambda, _k);
248 final double denom = MathTools.factorial(_k);
249 final double p = nom / denom;
250 return p;
251 }
252
253
254 /**
255 * Linearly transforms the specified value {@code value} from the specified value range
256 * interval {@code [valueRangeMin, valueRangeMax]} into the unit interval {@code [0, 1]}.
257 *
258 * @param value The value to transform.
259 * @param valueRangeMin Inclusive minimum of the value's original range.
260 * @param valueRangeMax Inclusive maximum of the value's original range.
261 * @return Image of the linear transform of the specified value into the unit inerval.
262 */
263 public static double transToUnitInterval(final double value,
264 final double valueRangeMin, final double valueRangeMax) {
265 final double valueRangeSize = valueRangeMax - valueRangeMin;
266 final double factor = (0. == valueRangeSize ? 0. : 1. / valueRangeSize);
267 final double transVal = (value - valueRangeMin) * factor;
268 return transVal;
269 }
270
271
272 /**
273 * Returns the logarithm of the specified value to the specified base.
274 *
275 * @param base The base for the logarithm.
276 * @param x The value to take the logarithm of.
277 * @return The logarithm of {@code x} to base {@code base}.
278 */
279 public static double log(final double base, final double x) {
280 return Math.log(x) / Math.log(base);
281 }
282
283
284 /**
285 * Compares the specified parameters.
286 * Returns a value smaller than 0 if the first parameter is smaller than the second;
287 * returns a value larger than 0 if the second parameter is smaller than the first;
288 * returns 0 if the parameters represent an equal value.<br />
289 * This method attempts to use the objects' own {@code compareTo} methods if both
290 * parameters are of same type and if that type is a standard Java type that inherits
291 * directly from {@code Number} and implements {@code Comparable}. If this fails,
292 * then the values returned by {@code doubleValue()} method are compared.
293 *
294 * @param <TX> A {@code Number} type.
295 * @param <TY> A {@code Number} type.
296 * @param x A number
297 * @param y A number.
298 * @return A value smaller than 0 if the first parameter is smaller than the second,
299 * a value larger than 0 if the second parameter is smaller than the first,
300 * 0 if the parameters represent an equal value.
301 */
302 public static <TX extends Number, TY extends Number> int compare(final TX x, final TY y) {
303
304 if (null == x && null == y)
305 return 0;
306
307 if (null == x || null == y)
308 throw new NullPointerException("Cannot comapre null to not null");
309
310 if (x instanceof Byte && y instanceof Byte)
311 return ((Byte) x).compareTo((Byte) y);
312
313 if (x instanceof Short && y instanceof Short)
314 return ((Short) x).compareTo((Short) y);
315
316 if (x instanceof Integer && y instanceof Integer)
317 return ((Integer) x).compareTo((Integer) y);
318
319 if (x instanceof Long && y instanceof Long)
320 return ((Long) x).compareTo((Long) y);
321
322 if (x instanceof Float && y instanceof Float)
323 return ((Float) x).compareTo((Float) y);
324
325 if (x instanceof Double && y instanceof Double)
326 return ((Double) x).compareTo((Double) y);
327
328 if (x instanceof BigInteger && y instanceof BigInteger)
329 return ((BigInteger) x).compareTo((BigInteger) y);
330
331 if (x instanceof BigDecimal && y instanceof BigDecimal)
332 return ((BigDecimal) x).compareTo((BigDecimal) y);
333
334 return Double.compare(x.doubleValue(), y.doubleValue());
335 }
336
337
338 /**
339 * Raises 10 to a positive integer power.
340 *
341 * @param exp Power to which to raise 10.
342 * @return If {@code exp} is negative - {@code 0}, otherwise {@code 10^exp}.
343 * @throws IllegalArgumentException If {@code exp > 9} as the result would be too
344 * large to be represented in a Java {@code int}.
345 */
346 public static int pow10(int exp) {
347 if (0 > exp)
348 return 0;
349 if (pow10.length <= exp)
350 throw new IllegalArgumentException("Cannot compute 10^" + exp + ": exponent too large");
351 return pow10[exp];
352 }
353
354
355 /**
356 * Raises 10 to a positive integer power.
357 *
358 * @param exp Power to which to raise 10.
359 * @return If {@code exp} is negative - {@code 0}, otherwise {@code 10^exp}.
360 * @throws IllegalArgumentException If {@code exp > 18} as the result would be too
361 * large to be represented in a Java {@code long}.
362 */
363 public static long pow10L(int exp) {
364 if (0 > exp)
365 return 0L;
366 if (pow10L.length <= exp)
367 throw new IllegalArgumentException("Cannot compute 10^" + exp + ": exponent too large");
368 return pow10L[exp];
369 }
370
371
372 /**
373 * Rounds the specified integer to the nearest specified power of 10.<br />
374 * For example:<br />
375 * {@code round(12345, 0)} = <em>12345</em> <br />
376 * {@code round(12345, 1)} = <em>12340</em> <br />
377 * {@code round(12345, 2)} = <em>12300</em> <br />
378 * {@code round(12345, 3)} = <em>12000</em> <br />
379 * {@code round(12345, 4)} = <em>10000</em> <br />
380 * {@code round(12345, 5)} = <em>0</em> <br />
381 * {@code round(54321, 5)} = <em>100000</em> <br />
382 *
383 * @param val A number.
384 * @param pow10 The power of 10 to which to round,
385 * only non-negative positive values are permitted.
386 * @return The closest integer to {@code val} that is a multiple of {@code 10^pow10}.
387 * Midway rounding is towards the even number.
388 */
389 public static int round(final int val, final int pow10) {
390
391 if (0 >= pow10)
392 return val;
393
394 final int scale = pow10(pow10);
395 final int quot = val / scale;
396 final int mod = val % scale;
397
398 // Round:
399 final int half = scale / 2;
400
401 if (mod > half)
402 return (quot + 1) * scale;
403
404 if (mod < half)
405 return quot * scale;
406
407 if (0 == quot % 2)
408 return quot * scale;
409 else
410 return (quot + 1) * scale;
411 }
412
413
414 /**
415 * Rounds the specified long integer to the nearest specified power of 10.<br />
416 * For example:<br />
417 * {@code round(12345, 0)} = <em>12345</em> <br />
418 * {@code round(12345, 1)} = <em>12340</em> <br />
419 * {@code round(12345, 2)} = <em>12300</em> <br />
420 * {@code round(12345, 3)} = <em>12000</em> <br />
421 * {@code round(12345, 4)} = <em>10000</em> <br />
422 * {@code round(12345, 5)} = <em>0</em> <br />
423 * {@code round(54321, 5)} = <em>100000</em> <br />
424 *
425 * @param val A number.
426 * @param pow10 The power of 10 to which to round,
427 * only non-negative positive values are permitted.
428 * @return The closest long integer to {@code val} that is a multiple of {@code 10^pow10}.
429 * Midway rounding is towards the even number.
430 */
431 public static long round(final long val, final int pow10) {
432
433 if (0 >= pow10)
434 return val;
435
436 final long scale = pow10L(pow10);
437 final long quot = val / scale;
438 final long mod = val % scale;
439
440 // Round:
441 final long half = scale / 2L;
442
443 if (mod > half)
444 return (quot + 1L) * scale;
445
446 if (mod < half)
447 return quot * scale;
448
449 if (0 == quot % 2L)
450 return quot * scale;
451 else
452 return (quot + 1L) * scale;
453 }
454
455
456 /**
457 * Rounds the specified number to the nearest specified power of 10.
458 * For example:<br />
459 * {@code round(12345.12345, 0)} = <em>12345.0</em> <br />
460 * {@code round(12345.12345, 1)} = <em>12340.0</em> <br />
461 * {@code round(12345.12345, 2)} = <em>12300.0</em> <br />
462 * {@code round(12345.12345, 3)} = <em>12000.0</em> <br />
463 * {@code round(12345.12345, 4)} = <em>10000.0</em> <br />
464 * {@code round(12345.12345, 5)} = <em>0.0</em> <br />
465 * {@code round(54321.0, 5)} = <em>100000.0</em> <br />
466 * {@code round(12345.12345, -1)} = <em>12345.1</em> <br />
467 * {@code round(12345.12345, -2)} = <em>12345.12</em> <br />
468 * {@code round(12345.12345, -3)} = <em>12345.123</em> <br />
469 * {@code round(12345.12345, -4)} = <em>12345.1234</em> <br />
470 * {@code round(12345.12345, -5)} = <em>12345.12345</em> <br />
471 * {@code round(12345.12345, -6)} = <em>12345.12345</em> <br />
472 *
473 * @param val A number.
474 * @param pow10 The power of 10 to which to round.
475 * @return The closest number to {@code val} that is a multiple of {@code 10^pow10}.
476 * Midway rounding is towards the even number.
477 */
478 public static double round(final double val, final int pow10) {
479 final double scale = Math.pow(10.0, pow10);
480 return Math.rint(val / scale) * scale;
481 }
482
483
484 /**
485 * Rounds the specified number to the nearest specified power of 10.
486 * For example:<br />
487 * {@code round(12345.12345, 0)} = <em>12345.0</em> <br />
488 * {@code round(12345.12345, 1)} = <em>12340.0</em> <br />
489 * {@code round(12345.12345, 2)} = <em>12300.0</em> <br />
490 * {@code round(12345.12345, 3)} = <em>12000.0</em> <br />
491 * {@code round(12345.12345, 4)} = <em>10000.0</em> <br />
492 * {@code round(12345.12345, 5)} = <em>0.0</em> <br />
493 * {@code round(54321.0, 5)} = <em>100000.0</em> <br />
494 * {@code round(12345.12345, -1)} = <em>12345.1</em> <br />
495 * {@code round(12345.12345, -2)} = <em>12345.12</em> <br />
496 * {@code round(12345.12345, -3)} = <em>12345.123</em> <br />
497 * {@code round(12345.12345, -4)} = <em>12345.1234</em> <br />
498 * {@code round(12345.12345, -5)} = <em>12345.12345</em> <br />
499 * {@code round(12345.12345, -6)} = <em>12345.12345</em> <br />
500 *
501 * @param val A number.
502 * @param pow10 The power of 10 to which to round.
503 * @return The closest number to {@code val} that is a multiple of {@code 10^pow10}.
504 * Midway rounding is towards the even number.
505 */
506 public static float round(final float val, final int pow10) {
507 final double scale = Math.pow(10.0, pow10);
508 return (float) (Math.rint(val / scale) * scale);
509 }
510
511 } // public class MathTools