[ODE] What is the ODE approach to denormals?

[GARY VANSICKLE]

[snip]
> -    if (v[i] > 0) {
> +    //  Denormals are a problem, because we divide by v[i], and then
> +    //  multiply that by 0. Alas, infinity times 0 is infinity (!)
> +    //  We also use v2[i], which is v[i] squared. Here's how the epsilons
> +    //  are chosen:
> +    //  float epsilon = 1.175494e-038 (smallest non-denormal number)
> +    //  double epsilon = 2.225074e-308 (smallest non-denormal number)
> +    //  For single precision, choose an epsilon such that v[i] squared is
> +    //  not a denormal; this is for performance.
> +    //  For double precision, choose an epsilon such that v[i] is not a
> +    //  denormal; this is for correctness.
> +#if defined( dSINGLE )
> +    if (v[i] > dReal(1e-19))
> +#else
> +    if (v[i] > dReal(1e-307))
> +#endif
> +    {
> 
> 
> 
> Is this the right approach?

In concept, probably yes.  In particulars, probably no.  I'd have to
educatedly guess that the smallest non-denormal numbers are going to vary
depending on the underlying hardware, not simply whether we're using float
or double.  This looks like a job for <limits>.  In particular, check this
bit of one popular <limits> out:

enum float_denorm_style {
   denorm_indeterminate = -1,
   denorm_absent = 0,
   denorm_present = 1
   };

Ouch: we might not even know if we have denormals!