//
// NumericUtilities.h
//
// This file contains the NumericUtilities template class. 
// This code is part of the Bioinformatics Template Library (BTL).
//
// Copyright (C) 1997 Birkbeck College, Malet Street, London WC1E 7HX, U.K.
// (classlib@mail.cryst.bbk.ac.uk)
// 
// This library is free software; you can
// redistribute it and/or modify it under the terms of
// the GNU Library General Public License as published by the Free
// Software Foundation; either version 2 of the License, or (at your
// option) any later version.  This library 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 Library General Public License for more details.
// You should have received a copy of the GNU Library General Public
// License along with this library; if not, write to the Free Software
// Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
///////////////////////////////////////////////////////////////////////////

#if !defined(NUMERICUTILITIES_H)
#define NUMERICUTILITIES_H

#include "NumericLimits.h"

/**#: [Description="This template class contains members which preform certain
     numerical functions. The design of this class is not object-oriented and
     the member functions could easily be provided as stand alone functions. The
     reasons for containing the functions within a class is to put them into a
     namespace and to make documentation easier. Standand C++  namespaces are
     not still not implemented in many C++ compilers but otherwise this method
     may well be better than putting the functions into a class. Most C++ design
     and automatic documentation tools do not deal with non-class parts of a
     library. We use Together/C++ for our documentation and this tool only
     produces documentation for classes and class members."]
    [Usage="e.g. unsigned int error =
     NumericUtilities&ltfloat&gt::AreNotEqual(a,b);"]
    [Restrictions="This template class only works for float, double, and long
       double types."]
    [Authors = "W.R.Pitt"]
    [Files = "<A HREF=./btl/NumericUtilities.h>NumericUtilities.h</A>"]
    [Friends="None"]
    [Dependencies="<A HREF=#NumericLimits>NumericLimits</A>"]*/

template <class T>
class NumericUtilities
{
public:

    /**#: [Description="Returns i where x = n * 10^i and 10>n>0"]*/
    
    static int 
    PowerOf10(const T& x);

    /**#: [Description="Finds the difference between two real numbers - the
           outputed result is the number of digits after the decimal point of
	   the difference divided by the smallest detectable difference. The
	   output is zero if the difference was smaller than the smallest
	   detectable difference. The maximum output depends on the precision of
	   the input numbers e.g for two float numbers the maximum might be 6
	   (FLT_DIG) and for double's it might be 15 (DBL_DIG).<P>
	   N.B. The result can be out by 1 due to rounding errors in the
	   calculation which can occur for example when a float is converted to
	   a double. However, the result will never be greater than the
	   precision given in the constants mentioned above."] */

    static unsigned int 
    AreNotEqual(const T& a, const T& b);
};

//.............................................................................


template <class T>
int NumericUtilities<T>::PowerOf10(const T& x)
{
    if (x == 0.0) return 1;

    int power10 = (int) log10(fabs(x));
    int lastPreciseValue;


    // Check answer and adjust it if necessary
    if (power10 != 0) lastPreciseValue =(int) floor(fabs(x * pow(10,-power10)));
    else lastPreciseValue = (int) floor(fabs(x));

    if (lastPreciseValue == 0) power10--;
    if (lastPreciseValue > 9) power10++;

    return power10;
}

//.............................................................................

template <class T>
unsigned int NumericUtilities<T>::AreNotEqual(const T& a, const T& b)
{
    const int DIGITS10 = NumericLimits<T>::digits10();

    // Deal with special cases

    // There is an infinity difference between any non-zero number and zero
    if (a == 0.0 && b == 0.0)
    	return 0;
    else if (a == 0.0 || b == 0.0)
    	return DIGITS10;

    // Find the absolute difference between the two input values (DELTA)
    T diff = fabs(a-b);

    // Differences of zero (possibly resulting from underflow)
    if (diff == 0.0) return 0;

    // Find largest absolute value between the 2 inputs and DELTA.
    // The last decimal place used in the subtraction is dependent on this
    // value.
    T maxVal = diff;
    if (fabs(a) > maxVal) maxVal = a;
    if (fabs(b) > maxVal) maxVal = b;

    // Find the value of the exponent (i) of this number (maxVal) in
    // the formula: maxVal = n * 10^i where 10>n>0
    int power10 = PowerOf10(maxVal);

    // The position of the last significant digit DELTA is dependent on
    // the precision used and the previous calculated value.
    int lastPreciseDigit = power10 - DIGITS10 + 1;

    // The calculated error is DELTA divided by the smallest detectable DELTA.
    // The error value that is output is the number of digits after the
    // decimal point of this calculated error.
    unsigned int error;
    if (lastPreciseDigit != 0)
    {
    	int errorPower = PowerOf10(diff * pow(10,-lastPreciseDigit)) + 1;
    	if (errorPower < 0)
    	    error = 0;
    	else if (errorPower>DIGITS10)
          // Can get rounding errors that result error being out by 1, for
          // example when float values are converted to double for fabs function
          error = DIGITS10;
      else
    	    error = errorPower;
    }
    else // Smallest detectable DELTA = 1.
    	error = PowerOf10(diff) + 1;
    return error;
}
#endif