List of usage examples for java.lang Math expm1
public static double expm1(double x)
From source file:org.renjin.primitives.random.ChiSquare.java
public static double pnchisq_raw(double x, double f, double theta, double errmax, double reltol, int itrmax, boolean lower_tail) { double lam, x2, f2, term, bound, f_x_2n, f_2n; double l_lam = -1., l_x = -1.; /* initialized for -Wall */ int n;//from w ww .ja v a 2 s. com boolean lamSml, tSml, is_r, is_b, is_it; double ans, u, v, t, lt, lu = -1; final double _dbl_min_exp = Math.log(2) * Double.MIN_EXPONENT; /*= -708.3964 for IEEE double precision */ if (x <= 0.) { if (x == 0. && f == 0.) { return lower_tail ? Math.exp(-0.5 * theta) : -Math.expm1(-0.5 * theta); } /* x < 0 or {x==0, f > 0} */ return lower_tail ? 0. : 1.; } if (!DoubleVector.isFinite(x)) { return lower_tail ? 1. : 0.; } if (theta < 80) { /* use 110 for Inf, as ppois(110, 80/2, lower.tail=FALSE) is 2e-20 */ double sum = 0, sum2 = 0, lambda = 0.5 * theta, pr = Math.exp(-lambda); double ans_inner; int i; /* we need to renormalize here: the result could be very close to 1 */ for (i = 0; i < 110; pr *= lambda / ++i) { sum2 += pr; sum += pr * Distributions.pchisq(x, f + 2 * i, lower_tail, false); if (sum2 >= 1 - 1e-15) { break; } } ans_inner = sum / sum2; return ans_inner; } lam = .5 * theta; lamSml = (-lam < _dbl_min_exp); if (lamSml) { /* MATHLIB_ERROR( "non centrality parameter (= %g) too large for current algorithm", theta) */ u = 0; lu = -lam;/* == ln(u) */ l_lam = Math.log(lam); } else { u = Math.exp(-lam); } /* evaluate the first term */ v = u; x2 = .5 * x; f2 = .5 * f; f_x_2n = f - x; if (f2 * SignRank.DBL_EPSILON > 0.125 && /* very large f and x ~= f: probably needs */ Math .abs(t = x2 - f2) < /* another algorithm anyway */ Math.sqrt(SignRank.DBL_EPSILON) * f2) { /* evade cancellation error */ /* t = exp((1 - t)*(2 - t/(f2 + 1))) / sqrt(2*M_PI*(f2 + 1));*/ lt = (1 - t) * (2 - t / (f2 + 1)) - 0.5 * Math.log(2 * Math.PI * (f2 + 1)); } else { /* Usual case 2: careful not to overflow .. : */ lt = f2 * Math.log(x2) - x2 - org.apache.commons.math.special.Gamma.logGamma(f2 + 1); } tSml = (lt < _dbl_min_exp); if (tSml) { if (x > f + theta + 5 * Math.sqrt(2 * (f + 2 * theta))) { /* x > E[X] + 5* sigma(X) */ return lower_tail ? 1. : 0.; /* FIXME: We could be more accurate than 0. */ } /* else */ l_x = Math.log(x); ans = term = t = 0.; } else { t = Math.exp(lt); ans = term = v * t; } for (n = 1, f_2n = f + 2., f_x_2n += 2.;; n++, f_2n += 2, f_x_2n += 2) { /* f_2n === f + 2*n * f_x_2n === f - x + 2*n > 0 <==> (f+2n) > x */ if (f_x_2n > 0) { /* find the error bound and check for convergence */ bound = t * x / f_x_2n; is_r = is_it = false; /* convergence only if BOTH absolute and relative error < 'bnd' */ if (((is_b = (bound <= errmax)) && (is_r = (term <= reltol * ans))) || (is_it = (n > itrmax))) { break; /* out completely */ } } /* evaluate the next term of the */ /* expansion and then the partial sum */ if (lamSml) { lu += l_lam - Math.log(n); /* u = u* lam / n */ if (lu >= _dbl_min_exp) { /* no underflow anymore ==> change regime */ v = u = Math.exp(lu); /* the first non-0 'u' */ lamSml = false; } } else { u *= lam / n; v += u; } if (tSml) { lt += l_x - Math.log(f_2n);/* t <- t * (x / f2n) */ if (lt >= _dbl_min_exp) { /* no underflow anymore ==> change regime */ t = Math.exp(lt); /* the first non-0 't' */ tSml = false; } } else { t *= x / f_2n; } if (!lamSml && !tSml) { term = v * t; ans += term; } } /* for(n ...) */ if (is_it) { // How to alert this message without an exception? //(_("pnchisq(x=%g, ..): not converged in %d iter."),x, itrmax); } return lower_tail ? ans : 1 - ans; }
From source file:uk.ac.diamond.scisoft.analysis.dataset.Maths.java
/** * expm1 - evaluate the exponential function - 1 on each element of the dataset * @param a//from w ww .j ava 2s . co m * @return dataset */ @SuppressWarnings("cast") public static AbstractDataset expm1(final AbstractDataset a) { final int isize; final IndexIterator it = a.getIterator(); AbstractDataset ds; final int dt = a.getDtype(); switch (dt) { case AbstractDataset.INT8: ds = AbstractDataset.zeros(a, AbstractDataset.FLOAT32); final byte[] i8data = ((ByteDataset) a).data; final float[] oi8data = ((FloatDataset) ds).getData(); for (int i = 0; it.hasNext();) { final byte ix = i8data[it.index]; float ox; ox = (float) (Math.expm1(ix)); oi8data[i++] = ox; } break; case AbstractDataset.INT16: ds = AbstractDataset.zeros(a, AbstractDataset.FLOAT32); final short[] i16data = ((ShortDataset) a).data; final float[] oi16data = ((FloatDataset) ds).getData(); for (int i = 0; it.hasNext();) { final short ix = i16data[it.index]; float ox; ox = (float) (Math.expm1(ix)); oi16data[i++] = ox; } break; case AbstractDataset.INT32: ds = AbstractDataset.zeros(a, AbstractDataset.FLOAT64); final int[] i32data = ((IntegerDataset) a).data; final double[] oi32data = ((DoubleDataset) ds).getData(); for (int i = 0; it.hasNext();) { final int ix = i32data[it.index]; double ox; ox = (double) (Math.expm1(ix)); oi32data[i++] = ox; } break; case AbstractDataset.INT64: ds = AbstractDataset.zeros(a, AbstractDataset.FLOAT64); final long[] i64data = ((LongDataset) a).data; final double[] oi64data = ((DoubleDataset) ds).getData(); for (int i = 0; it.hasNext();) { final long ix = i64data[it.index]; double ox; ox = (double) (Math.expm1(ix)); oi64data[i++] = ox; } break; case AbstractDataset.ARRAYINT8: ds = AbstractDataset.zeros(a, AbstractDataset.ARRAYFLOAT32); isize = a.getElementsPerItem(); final byte[] ai8data = ((CompoundByteDataset) a).data; final float[] oai8data = ((CompoundFloatDataset) ds).getData(); for (int i = 0; it.hasNext();) { for (int j = 0; j < isize; j++) { final byte ix = ai8data[it.index + j]; float ox; ox = (float) (Math.expm1(ix)); oai8data[i++] = ox; } } break; case AbstractDataset.ARRAYINT16: ds = AbstractDataset.zeros(a, AbstractDataset.ARRAYFLOAT32); isize = a.getElementsPerItem(); final short[] ai16data = ((CompoundShortDataset) a).data; final float[] oai16data = ((CompoundFloatDataset) ds).getData(); for (int i = 0; it.hasNext();) { for (int j = 0; j < isize; j++) { final short ix = ai16data[it.index + j]; float ox; ox = (float) (Math.expm1(ix)); oai16data[i++] = ox; } } break; case AbstractDataset.ARRAYINT32: ds = AbstractDataset.zeros(a, AbstractDataset.ARRAYFLOAT64); isize = a.getElementsPerItem(); final int[] ai32data = ((CompoundIntegerDataset) a).data; final double[] oai32data = ((CompoundDoubleDataset) ds).getData(); for (int i = 0; it.hasNext();) { for (int j = 0; j < isize; j++) { final int ix = ai32data[it.index + j]; double ox; ox = (double) (Math.expm1(ix)); oai32data[i++] = ox; } } break; case AbstractDataset.ARRAYINT64: ds = AbstractDataset.zeros(a, AbstractDataset.ARRAYFLOAT64); isize = a.getElementsPerItem(); final long[] ai64data = ((CompoundLongDataset) a).data; final double[] oai64data = ((CompoundDoubleDataset) ds).getData(); for (int i = 0; it.hasNext();) { for (int j = 0; j < isize; j++) { final long ix = ai64data[it.index + j]; double ox; ox = (double) (Math.expm1(ix)); oai64data[i++] = ox; } } break; case AbstractDataset.FLOAT32: ds = AbstractDataset.zeros(a, AbstractDataset.FLOAT32); final float[] f32data = ((FloatDataset) a).data; final float[] of32data = ((FloatDataset) ds).getData(); for (int i = 0; it.hasNext();) { final float ix = f32data[it.index]; float ox; ox = (float) (Math.expm1(ix)); of32data[i++] = ox; } break; case AbstractDataset.FLOAT64: ds = AbstractDataset.zeros(a, AbstractDataset.FLOAT64); final double[] f64data = ((DoubleDataset) a).data; final double[] of64data = ((DoubleDataset) ds).getData(); for (int i = 0; it.hasNext();) { final double ix = f64data[it.index]; double ox; ox = (double) (Math.expm1(ix)); of64data[i++] = ox; } break; case AbstractDataset.ARRAYFLOAT32: ds = AbstractDataset.zeros(a, AbstractDataset.ARRAYFLOAT32); isize = a.getElementsPerItem(); final float[] af32data = ((CompoundFloatDataset) a).data; final float[] oaf32data = ((CompoundFloatDataset) ds).getData(); for (int i = 0; it.hasNext();) { for (int j = 0; j < isize; j++) { final float ix = af32data[it.index + j]; float ox; ox = (float) (Math.expm1(ix)); oaf32data[i++] = ox; } } break; case AbstractDataset.ARRAYFLOAT64: ds = AbstractDataset.zeros(a, AbstractDataset.ARRAYFLOAT64); isize = a.getElementsPerItem(); final double[] af64data = ((CompoundDoubleDataset) a).data; final double[] oaf64data = ((CompoundDoubleDataset) ds).getData(); for (int i = 0; it.hasNext();) { for (int j = 0; j < isize; j++) { final double ix = af64data[it.index + j]; double ox; ox = (double) (Math.expm1(ix)); oaf64data[i++] = ox; } } break; case AbstractDataset.COMPLEX64: ds = AbstractDataset.zeros(a, AbstractDataset.COMPLEX64); final float[] c64data = ((ComplexFloatDataset) a).data; final float[] oc64data = ((ComplexFloatDataset) ds).getData(); for (int i = 0; it.hasNext();) { final float ix = c64data[it.index]; final float iy = c64data[it.index + 1]; float tf; float ox; float oy; tf = (float) (Math.expm1(ix)); ox = (float) (tf * Math.cos(iy)); oy = (float) (tf * Math.sin(iy)); oc64data[i++] = ox; oc64data[i++] = oy; } break; case AbstractDataset.COMPLEX128: ds = AbstractDataset.zeros(a, AbstractDataset.COMPLEX128); final double[] c128data = ((ComplexDoubleDataset) a).data; final double[] oc128data = ((ComplexDoubleDataset) ds).getData(); for (int i = 0; it.hasNext();) { final double ix = c128data[it.index]; final double iy = c128data[it.index + 1]; double tf; double ox; double oy; tf = (double) (Math.expm1(ix)); ox = (double) (tf * Math.cos(iy)); oy = (double) (tf * Math.sin(iy)); oc128data[i++] = ox; oc128data[i++] = oy; } break; default: throw new IllegalArgumentException( "expm1 supports integer, compound integer, real, compound real, complex datasets only"); } ds.setName(a.getName()); addFunctionName(ds, "expm1"); return ds; }