Hi to everybody,
I am looking for a way to determine the floating-point epsilon value (i.e., the minimum number that should be added to 1.0 in order to get something more than 1.0, like the constants FLT_EPSILON and DBL_EPSILON in float.h). I would like to use it for a generic function, something like a generic epsilon[T] function:
import math
type
Func[T] = proc(x : T) : T
proc der[T : SomeReal](f : Func, pt: T) : T =
let eps = sqrt(epsilon[T]())
result = (f(pt + eps) - f(pt)) / eps
(My real case does not involve the calculation of derivatives, but something more complex; however, this example should give the idea). I am looking for something like Ada's Machine_Epsilon attribute: does Nim provide something similar?
Varriount: Thanks for the reply, this is exactly what I did at the end.
Sadly, I think this imply that the answer to my question "does Nim provide something similar?" is… no. I was looking for a solution in pure Nim that can be used with generics, as in my example above.
In Julia, I would implement it in the following way:
myEps(::Type{Float32}) = FLT_EPSILON # Whatever myEps(::Type{Float64}) = DBL_EPSILON # Whatever
This has the advantage that Julia will complain whenever I call e.g. myEps(Int32).
My point is to have a generic function epsilon[T], like the one in the example I gave above, which gets specialized for T equal to float32 and float64 but causes a compilation error when used with other types (e.g., int16).
So far, I've been able to put together this implementation:
var FLT_EPSILON {. importc: "FLT_EPSILON" header: "float.h" .} : cfloat
var DBL_EPSILON {. importc: "DBL_EPSILON" header: "float.h" .} : cdouble
proc epsilon[T: SomeReal]() : T =
when T is cfloat:
result = FLT_EPSILON
when T is cdouble:
result = DBL_EPSILON
echo("FLT_EPSILON = " & $FLT_EPSILON)
echo("DBL_EPSILON = " & $DBL_EPSILON)
echo("epsilon with 64-bit FP: " & $float(epsilon[cfloat]()))
echo("epsilon with 64-bit FP: " & $float(epsilon[cdouble]()))
Sadly, it does not work:
FLT_EPSILON = 1.192092895507812e-07 DBL_EPSILON = 2.220446049250313e-16 epsilon with 32-bit FP: 2.220446049250313e-16 epsilon with 64-bit FP: 2.220446049250313e-16
Moreover, there are a few things that annoy me:
If I am able to put together some working implementation for epsilon[T], would Araq be interested in a PR for lib/pure/math.nim? (Also, would this file be the right place?) I would like to add other FP properties from float.h as well (e.g., *_MAX, *_MANT_DIG, and so on).
But is there any correct solution that can be exhaustive
Come on, be creative, when T is float32 ... elif T is float64 else: {.error: "unsupported floating point type".} does the trick.
Araq, thanks for your answers. I keep forgetting that elif can be used with while too. (And I wasn't aware of the existence of the error pragma.)
Regarding the problem in the code above, I found that the same applies when I use a template:
var FLT_EPSILON {. importc: "FLT_EPSILON" header: "float.h" .} : cfloat
var DBL_EPSILON {. importc: "DBL_EPSILON" header: "float.h" .} : cdouble
template eps(t : typedesc) : expr =
when t is float32:
FLT_EPSILON
elif t is float64:
DBL_EPSILON
else:
{. error: "unsupported floating-point type" .}
echo("FLT_EPSILON = " & $FLT_EPSILON)
echo("DBL_EPSILON = " & $DBL_EPSILON)
echo("eps with 32-bit FP: " & $eps(float32))
echo("eps with 64-bit FP: " & $eps(float64))
The output is still wrong:
FLT_EPSILON = 1.192092895507812e-07 DBL_EPSILON = 2.220446049250313e-16 eps with 32-bit FP: 1.192092895507812e-07 eps with 64-bit FP: 1.192092895507812e-07
However, at the end I got a right way to make the Nim compiler do what I want:
template epsilon(T : typedesc[float32]): float32 = FLT_EPSILON
template epsilon(T : typedesc[float64]): float64 = DBL_EPSILON
This solution finally works! And the two lines above are ok even if I call epsilon(cfloat) or epsilon(cdouble). I really like this solution, as it is shorter and much similar to what I would have done using Julia. (Another big plus is that, in the generated C code, no epsilon function needs to be called, because FLT_EPSILON and DBL_EPSILON are substituted directly in the call site.)
Many, many thanks for this. I'll create a PR soon.