Re: [lmi] Is bourn_cast demonstrably correct? (+PATCH)

[Greg Chicares]

On 2017-03-20 20:46, Vadim Zeitlin wrote:
> On Mon, 20 Mar 2017 20:04:03 +0000 Greg Chicares <address@hidden> wrote:
[...]
> GC> BTW, let me share my short-term plans for this facility. First, I asked
> GC> myself whether it should just forbid floating point; the answer is no,
> GC> because we really do have code that converts double to integral types
> GC> using value_cast<>(), and the whole point of bourn_cast<>() is to
> GC> replace boost::numeric_cast<>() in value_cast<>().
> 
>  Is the code doing this really supposed to lose precision? I.e. should it
> convert M_PI to 3? Or is it only meant to convert 3.0 to 3 and throw for
> the numbers with non zero fractional part?

A long time ago, I wrote this:

/// Function template numeric_value_cast() wraps boost::numeric_cast
/// to make it DWISOTT according to the boost-1.31.0 documentation:
///   "An exception is thrown when a runtime value-preservation
///   check fails."
/// The problem is that
///   boost::numeric_cast<int>(2.71828);
/// returns the integer 2 without throwing, but 2.71828 and 2 are
/// different values. It seems unreasonable to call truncation a
/// value-preserving relation.

template<typename To, typename From>
To numeric_value_cast(From const& from)
{
    To result = boost::numeric_cast<To>(from);
    if(result == from)
        {
        return result;
        }
    else
[...throw]

Now I believe the flaw was in the documentation, not the code.
I can't be certain what the documentation should have said, but,
as Cacciola's paper:
  http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2005/n1879.htm
explains, there are two problems in mapping between floating and
integral types:

  (1) loss of range: e.g., DBL_MAX --> int, where static_cast produces
undefined behavior [4.9/1] and definitely gives a wrong answer, because
there really is no right answer;

  (2) loss of precision: e.g., M_PI --> int, or ULLONG_MAX --> float,
where static_cast gives a well-defined result [4.9/2] that is as good
an approximation as it can be, plus or minus one ulp;

and Henney's real intention was probably to throw on (1) and pass (2)
through without complaint: IOW, to write a static_cast wrapper that
simply blocks obvious pitfalls like casting -1U into FFFFFFFF.

lmi uses the numeric_value_cast() above in 'calendar_date.cpp' in
birthdate_limit::operator()(), which attacks this question...

/// Determine a birthdate limit, iteratively.
///
/// To be age A on date D, one must have been born on a date B in
/// [Bmin, Bmax]. Problem: to find Bmin or Bmax, given A and D.

...with a polynomial solver (presumably a closed-form solution exists,
but the arithmetic would be daunting). The solver is tuned in this case
to produce a result rounded to zero decimals, so numeric_value_cast()
converts it to integer correctly, without loss of precision (because
both int and double can represent any JDN within lmi's lifecycle).

numeric_value_cast() is also used indirectly, through value_cast(),
in situations like:
  ledger_text_formats.cpp:    int age = value_cast<int>(invar().Age);
in the context of:
  class LedgerInvariant { ... double          Age; ... }
where type double is treated as almost a union of all arithmetic
types (so that we can apply algorithms to a vector<double>), but
many of the members have integer semantics (like "int age" above).

IOW, we've used numeric_value_cast() broadly, but we have never seen
it throw an "inexact" exception, so, arguably, the runtime value-
preservation check provides no actual benefit, and might be eliminated;
but OTOH it's cheap insurance, so why throw it away?

I think the best answer is to keep bourn_cast as a "guarded" wrapper
for static_cast, which throws only when the behavior of static_cast
is surprising (-1U --> FFFFFFFF) or undefined (DBL_MAX --> int); and
to provide a wrapper-wrapper like numeric_value_cast() above if value
preservation is important.

Hmmm...in that wrapper-wrapper:
    To result = boost::numeric_cast<To>(from);
    if(result == from)
is "result == from" the best we can do, or would
    double x;
    if(0.0 == std::modf(from, x))
be better? No, I guess not: either way, we're going to compare two
floating-point values for exact equality (and this is one case where
we really do want to do that); but with modf(), we potentially call
fesetround() twice, and we have to choose among {modf, modff, modfl}.
Testing "result == from" is simpler and better.