Back to index

wims  3.65+svn20090927
MoreMath.java
Go to the documentation of this file.
00001 /*
00002  * Copyright (C) 2006-2008 Mihai Preda.
00003  *
00004  * Licensed under the Apache License, Version 2.0 (the "License");
00005  * you may not use this file except in compliance with the License.
00006  * You may obtain a copy of the License at
00007  *
00008  *      http://www.apache.org/licenses/LICENSE-2.0
00009  *
00010  * Unless required by applicable law or agreed to in writing, software
00011  * distributed under the License is distributed on an "AS IS" BASIS,
00012  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
00013  * See the License for the specific language governing permissions and
00014  * limitations under the License.
00015  *
00016  */
00017 
00018 package org.javia.arity;
00019 
00020 class MoreMath extends BaseMath {
00021     private static final double LOG2E  = 1.4426950408889634074;
00022 
00023     public static final double asinh(double x) {
00024         return (x < 0) ? -asinh(-x) : log(x + x + 1/(Math.sqrt(x*x + 1) + x));
00025     }
00026     
00027     public static final double acosh(double x) {
00028         return log(x + x - 1/(Math.sqrt(x*x - 1) + x));
00029     }
00030 
00031     public static final double atanh(double x) {
00032         return (x < 0) ? -atanh(-x) : 0.5 * log(1. + (x + x)/(1 - x));
00033     }
00034 
00035     public static final double trunc(double x) {
00036         return x >= 0 ? Math.floor(x) : Math.ceil(x);
00037     }
00038 
00039     public static final double gcd(double x, double y) {
00040         //double remainder = y;
00041         if (Double.isNaN(x) || Double.isNaN(y) ||
00042             Double.isInfinite(x) || Double.isInfinite(y)) {
00043             return Double.NaN;
00044         }
00045         x = Math.abs(x);
00046         y = Math.abs(y);
00047         double save;
00048         while (y > 1e-12) {
00049             save = y;
00050             y = x % y;
00051             x = save;
00052             //Log.log(y);
00053         } 
00054         return x > 1e-10 ? x : 0;
00055     }
00056  
00057     public static final double lgamma(double x) {
00058         double tmp = x + 5.2421875; //== 607/128. + .5;
00059         return 0.9189385332046727418 //LN_SQRT2PI, ln(sqrt(2*pi))
00060             + log(
00061                   0.99999999999999709182 +
00062                   57.156235665862923517     / ++x +
00063                   -59.597960355475491248    / ++x +
00064                   14.136097974741747174     / ++x +
00065                   -0.49191381609762019978   / ++x +
00066                   .33994649984811888699e-4  / ++x +
00067                   .46523628927048575665e-4  / ++x +
00068                   -.98374475304879564677e-4 / ++x +
00069                   .15808870322491248884e-3  / ++x +
00070                   -.21026444172410488319e-3 / ++x +
00071                   .21743961811521264320e-3  / ++x +
00072                   -.16431810653676389022e-3 / ++x +
00073                   .84418223983852743293e-4  / ++x +
00074                   -.26190838401581408670e-4 / ++x +
00075                   .36899182659531622704e-5  / ++x
00076                   )
00077             + (tmp-4.7421875)*log(tmp) - tmp
00078             ;
00079     }
00080 
00081     static final double FACT[] = {
00082         1.0,
00083         40320.0,
00084         2.0922789888E13,
00085         6.204484017332394E23,
00086         2.631308369336935E35,
00087         8.159152832478977E47,
00088         1.2413915592536073E61,
00089         7.109985878048635E74,
00090         1.2688693218588417E89,
00091         6.1234458376886085E103,
00092         7.156945704626381E118,
00093         1.8548264225739844E134,
00094         9.916779348709496E149,
00095         1.0299016745145628E166,
00096         1.974506857221074E182,
00097         6.689502913449127E198,
00098         3.856204823625804E215,
00099         3.659042881952549E232,
00100         5.5502938327393044E249,
00101         1.3113358856834524E267,
00102         4.7147236359920616E284,
00103         2.5260757449731984E302,
00104     };
00105 
00106     public static final double factorial(double x) {
00107         if (x < 0) { // x <= -1 ?
00108             return Double.NaN;
00109         }
00110         if (x <= 170) {
00111             if (Math.floor(x) == x) {
00112                 int n = (int)x;
00113                 double extra = x;
00114                 switch (n & 7) {
00115                 case 7: extra *= --x;
00116                 case 6: extra *= --x;
00117                 case 5: extra *= --x;
00118                 case 4: extra *= --x;
00119                 case 3: extra *= --x;
00120                 case 2: extra *= --x;
00121                 case 1: return FACT[n >> 3] * extra;
00122                 case 0: return FACT[n >> 3];
00123                 }
00124             }
00125         }
00126         return exp(lgamma(x));
00127     }
00128 
00129     public static final double comb(double n, double k) {
00130         if (n < 0 || k < 0) { return Double.NaN; }
00131         if (n < k) { return 0; }
00132         if (Math.floor(n) == n && Math.floor(k) == k) {
00133             k = Math.min(k, n-k);
00134             if (n <= 170 && 12 < k && k <= 170) {
00135                 return factorial(n)/factorial(k)/factorial(n-k);
00136             } else {
00137                 double r = 1, diff = n-k;
00138                 for (double i = k; i > .5 && r < Double.POSITIVE_INFINITY; --i) {
00139                     r *= (diff+i)/i;
00140                 }
00141                 return r;
00142             }
00143         } else {
00144             return exp(lgamma(n) - lgamma(k) - lgamma(n-k));
00145         }
00146     }
00147 
00148     public static final double perm(double n, double k) {
00149         if (n < 0 || k < 0) { return Double.NaN; }
00150         if (n < k) { return 0; }
00151         if (Math.floor(n) == n && Math.floor(k) == k) {
00152             if (n <= 170 && 10 < k && k <= 170) {
00153                 return factorial(n)/factorial(n-k);
00154             } else {
00155                 double r = 1, limit = n-k+.5;
00156                 for (double i = n; i > limit && r < Double.POSITIVE_INFINITY; --i) {
00157                     r *= i;
00158                 }
00159                 return r;
00160             }
00161         } else {
00162             return exp(lgamma(n) - lgamma(n-k));
00163         }
00164     }
00165 
00166     public static final double log2(double x) {
00167         return log(x) * LOG2E;
00168     }
00169 
00170     private static final boolean isPiMultiple(double x) {
00171         return x % Math.PI == 0;
00172     }
00173 
00174     public static final int intLog10(double x) {
00175         //an alternative implem is using a for loop.
00176         return (int)Math.floor(log10(x));
00177         //return (int)log10(x);
00178     }
00179 
00180     public static final double intExp10(int exp) {
00181         return Double.parseDouble("1E" + exp);
00182     }
00183 }