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.
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:
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.
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.