/*
* DrMath.java
*
* Copyright (c) 2002-2015 Alexei Drummond, Andrew Rambaut and Marc Suchard
*
* This file is part of BEAST.
* See the NOTICE file distributed with this work for additional
* information regarding copyright ownership and licensing.
*
* BEAST is free software; you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as
* published by the Free Software Foundation; either version 2
* of the License, or (at your option) any later version.
*
* BEAST is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with BEAST; if not, write to the
* Free Software Foundation, Inc., 51 Franklin St, Fifth Floor,
* Boston, MA 02110-1301 USA
*/
package dr.math.functionEval;
import java.io.PrintStream;
/**
* This class implements additional mathematical functions
* and determines the parameters of the floating point representation.
*
* @author Didier H. Besset
*/
public final class DrMath
{
/**
* Typical meaningful precision for numerical calculations.
*/
static private double defaultNumericalPrecision = 0;
/**
* Typical meaningful small number for numerical calculations.
*/
static private double smallNumber = 0;
/**
* Radix used by floating-point numbers.
*/
static private int radix = 0;
/**
* Largest positive value which, when added to 1.0, yields 0.
*/
static private double machinePrecision = 0;
/**
* Largest positive value which, when subtracted to 1.0, yields 0.
*/
static private double negativeMachinePrecision = 0;
/**
* Smallest number different from zero.
*/
static private double smallestNumber = 0;
/**
* Largest possible number
*/
static private double largestNumber = 0;
/**
* Largest argument for the exponential
*/
static private double largestExponentialArgument = 0;
/**
* Values used to compute human readable scales.
*/
private static final double scales[] = {1.25, 2, 2.5, 4, 5, 7.5, 8, 10};
private static final double semiIntegerScales[] = {2, 2.5, 4, 5, 7.5, 8, 10};
private static final double integerScales[] = {2, 4, 5, 8, 10};
private static void computeLargestNumber()
{
double floatingRadix = getRadix();
double fullMantissaNumber = 1.0d -
floatingRadix * getNegativeMachinePrecision();
while ( !Double.isInfinite( fullMantissaNumber) )
{
largestNumber = fullMantissaNumber;
fullMantissaNumber *= floatingRadix;
}
}
private static void computeMachinePrecision()
{
double floatingRadix = getRadix();
double inverseRadix = 1.0d / floatingRadix;
machinePrecision = 1.0d;
double tmp = 1.0d + machinePrecision;
while ( tmp - 1.0d != 0.0d)
{
machinePrecision *= inverseRadix;
tmp = 1.0d + machinePrecision;
}
}
private static void computeNegativeMachinePrecision()
{
double floatingRadix = getRadix();
double inverseRadix = 1.0d / floatingRadix;
negativeMachinePrecision = 1.0d;
double tmp = 1.0d - negativeMachinePrecision;
while ( tmp - 1.0d != 0.0d)
{
negativeMachinePrecision *= inverseRadix;
tmp = 1.0d - negativeMachinePrecision;
}
}
private static void computeRadix()
{
double a = 1.0d;
double tmp1, tmp2;
do { a += a;
tmp1 = a + 1.0d;
tmp2 = tmp1 - a;
} while ( tmp2 - 1.0d != 0.0d);
double b = 1.0d;
while ( radix == 0)
{
b += b;
tmp1 = a + b;
radix = (int) ( tmp1 - a);
}
}
private static void computeSmallestNumber()
{
double floatingRadix = getRadix();
double inverseRadix = 1.0d / floatingRadix;
double fullMantissaNumber = 1.0d - floatingRadix * getNegativeMachinePrecision();
while ( fullMantissaNumber != 0.0d )
{
smallestNumber = fullMantissaNumber;
fullMantissaNumber *= inverseRadix;
}
}
public static double defaultNumericalPrecision()
{
if ( defaultNumericalPrecision == 0 )
defaultNumericalPrecision = Math.sqrt( getMachinePrecision());
return defaultNumericalPrecision;
}
/**
* @return boolean true if the difference between a and b is
* less than the default numerical precision
* @param a double
* @param b double
*/
public static boolean equal( double a, double b)
{
return equal( a, b, defaultNumericalPrecision());
}
/**
* @return boolean true if the relative difference between a and b
* is less than precision
* @param a double
* @param b double
* @param precision double
*/
public static boolean equal( double a, double b, double precision)
{
double norm = Math.max( Math.abs(a), Math.abs( b));
return norm < precision || Math.abs( a - b) < precision * norm;
}
public static double getLargestExponentialArgument()
{
if ( largestExponentialArgument == 0 )
largestExponentialArgument = Math.log(getLargestNumber());
return largestExponentialArgument;
}
/**
* (c) Copyrights Didier BESSET, 1999, all rights reserved.
*/
public static double getLargestNumber()
{
if ( largestNumber == 0 )
computeLargestNumber();
return largestNumber;
}
public static double getMachinePrecision()
{
if ( machinePrecision == 0 )
computeMachinePrecision();
return machinePrecision;
}
public static double getNegativeMachinePrecision()
{
if ( negativeMachinePrecision == 0 )
computeNegativeMachinePrecision();
return negativeMachinePrecision;
}
public static int getRadix()
{
if ( radix == 0 )
computeRadix();
return radix;
}
public static double getSmallestNumber()
{
if ( smallestNumber == 0 )
computeSmallestNumber();
return smallestNumber;
}
public static void printParameters( PrintStream printStream)
{
printStream.println( "Floating-point machine parameters");
printStream.println( "---------------------------------");
printStream.println( " ");
printStream.println( "radix = "+ getRadix());
printStream.println( "Machine precision = "
+ getMachinePrecision());
printStream.println( "Negative machine precision = "
+ getNegativeMachinePrecision());
printStream.println( "Smallest number = "+ getSmallestNumber());
printStream.println( "Largest number = "+ getLargestNumber());
return;
}
public static void reset()
{
defaultNumericalPrecision = 0;
smallNumber = 0;
radix = 0;
machinePrecision = 0;
negativeMachinePrecision = 0;
smallestNumber = 0;
largestNumber = 0;
}
/**
* This method returns the specified value rounded to
* the nearest integer multiple of the specified scale.
*
* @param value number to be rounded
* @param scale defining the rounding scale
* @return rounded value
*/
public static double roundTo( double value, double scale)
{
return Math.round( value / scale) * scale;
}
/**
* Round the specified value upward to the next scale value.
* @param the value to be rounded.
* @param a fag specified whether integer scale are used, otherwise double scale is used.
* @return a number rounded upward to the next scale value.
*/
public static double roundToScale( double value, boolean integerValued)
{
double[] scaleValues;
int orderOfMagnitude = (int) Math.floor( Math.log( value) / Math.log( 10.0));
if ( integerValued )
{
orderOfMagnitude = Math.max( 1, orderOfMagnitude);
if ( orderOfMagnitude == 1)
scaleValues = integerScales;
else if ( orderOfMagnitude == 2)
scaleValues = semiIntegerScales;
else
scaleValues = scales;
}
else
scaleValues = scales;
double exponent = Math.pow( 10.0, orderOfMagnitude);
double rValue = value / exponent;
for ( int n = 0; n < scaleValues.length; n++)
{
if ( rValue <= scaleValues[n])
return scaleValues[n] * exponent;
}
return exponent; // Should never reach here
}
/**
* (c) Copyrights Didier BESSET, 1999, all rights reserved.
* @return double
*/
public static double smallNumber()
{
if ( smallNumber == 0 )
smallNumber = Math.sqrt( getSmallestNumber());
return smallNumber;
}
}