Infinitely small numbers

[Stewbond]

This evening I wrote a function to calculate a derivative:

typedef double (*func)double;

double fprime( func f , double x )
{
    const double dx = 0.0000000001;
    return ( f(x + dx) - f(x) ) / dx;
}

I don't like line 5 above. dx should be an infinitely small number.

The value chosen above may work for a sin function, but if I'm calculating molar mass, it's way too big. Even worse, if I'm doing orbital mechanics, it's so small that it will be less than the LSB of the mantissa of a floating point number which will have the disastrous effects of not incremental x at all.

Is there a function available in std:: that will return the value of the LSB of the mantissa? This would let me calculate dx properly for any value of x.

I've been checking in <limits> but I haven't found anything and I don't want to dig into IEEE754 to do it bitwise myself.

Edit: I've continued researching and have found that what I am looking for is the ULP (Unit in last place) of the floating point number.

[Athar]

There is epsilon() (in <limits>), however you might find that it's actually too small to be useful in this case.

[Stewbond]

Ah perfect!

I was starting to look at boost::math::float_advance() in boost/math/special_functions.

I'm glad there's something in std. I will certainly multiply epsilon() by a factor as minimum resolution operations will be prone to larger errors than slightly larger operations.

Edit: I checked out epsilon() and it only gives me the difference between 1.0 and the next largest number. This would be different than the difference between 1000000 and the next largest number so I can't use this function.
http://en.cppreference.com/w/cpp/types/numeric_limits/epsilon

I think that boost::math::float_advance() is the best way to do it. I took a quick peak at what it does. If I move it forward by 3 ULPs, this is what I get:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21

#include <boost/math/special_functions/next.hpp> #include <iomanip> union ftoi { float f; unsigned int i; }; int main() { ftoi x0, x1; x0.f = 1.0f; std::cout << std::setw(8) << std::setfill('0') << std::hex << x0.i << std::endl; x1.f = boost::math::float_advance( x0.f, 3 ); std::cout << std::setw(8) << std::setfill('0') << std::hex << x1.i << std::endl; }

3f800000 3f800003 Press any key to continue . . .

Which is very cool. My function now looks like this:

double fprime( func f , double x)
{
    const double xn = boost::math::float_advance( x, 1 );
    return ( f(xn) - f(x) ) / (xn - x);
}

I know that 1 should make this prone to errors, but I just evaluated the cos(x) = 0 (which calculates to pi/2) and got it accurate to 17 places (which was the most I could display in the console). I'm happy with that.

[Cubbi]

Just a comment, the example on that cppreference page shows exactly how to scale the epsilon to get the desired ulps. And the seealso points at the <cmath> functions to get the next/previous representible float/double given any starting point.

[Stewbond]

Oh, I see. You're right!

template<class T>
typename std::enable_if<!std::numeric_limits<T>::is_integer, T>::type
    small_delta(T x, int ulp)
{
    return abs(x) * std::numeric_limits<T>::epsilon() * ulp;
}

double fprime( func f , double x)
{
    const double dx = small_delta( x, 5 );
    return ( f(x + dx) - f(x) ) / (dx);
}

I saw those cmath functions: nextafter, nexttoward. But they are C++11 and not in my implementation (VS2010).