Re[5]: [Help-glpk] proposal for gnu lp/mip low-level format

[Michael Hennebry]

On Fri, 30 Jul 2004, Andrew Makhorin wrote:

> >FORTRAN can't do symbol tables?
>
> I meant that one needs no symbol tables to read files in the proposed
> formats.
>
> >[regarding rationals, e.g. 2/3 != .666666667]
>
> As well as 3*(2/3) != (3*2)/3 in finite precision arithimetic.
>
> >That would be determined by the rounding mode,
>
> I think that changing the rounding mode is a bad idea (though I agree

Yup.  Absent special circumstances, round to nearest is best.

> that some specific applications, for example, for interval arithmetic,
> needs to control the rounding mode). All experts in numerical analysis
> suggest that only the "correct" rounding (i.e. rounding to the nearest
> representable floating-point number) must be used in all numerical
> computations, because it provides highest accuracy unlike other rounding
> modes.
>
> > not the optimization option.
>
> It must be so, but it is not. For example, the following fragment
>
>    eps = 1.0;
>    while (1.0 + eps > 1.0) eps *= 0.5;
>    eps += eps;
>
> allows computing the machine precision. Many "clever" compilers trying
> to "optimize" this fragment think that the condition 1.0 + eps > 1.0
> is always true as in exact arithmetic, ...

Ouch.  I'd known that some compilers are willing to optimize in ways
that might change the result, but I hadn't known about that one.

It doesn't affect the computation of 2/3 though.

> ...                                    other compilers, again due to
> optimization options, hold temporary results in floating-point registers
> which might be wider than the declared precision of eps that involves a
> wrong result, etc. I met this many times in my practice. The only right

gcc on Intel cpus for example.
This might affect the computation of 2/3.
Double rounding might cause the result to be off by more than half an ulp.
Converting from a decimal fraction presents the same problem
and could be worse if the reader is not well-written.

> way is *not to change* the precision and the order of computations at
> all, however, this rule is ignored by most compiler developers.

Most, but not all.  gcc has flags that distinguish between optimizations
that could change the result and those that won't.

> As to using rational numbers as input data, I think this would be
> actual, say, in a computer algebra system. In case of lp/mip it is
> always possible to reduce rational numbers to integer ones by
> appropriate scaling constraint and objective coefficients and bounds
> of variables. For example, let there be a constraint:
>
>    (1/3) * x1 + (1/2) * x2 + (2/7) * x3 >= 5/6
>
> Multiplying its both sides by lcm(3, 2, 7, 6) = 42 we have:
>
>    14 * x1 + 21 * x2 + 12 * x3 >= 35
>
> i.e. using integer numbers (which are always read exactly even being
> floating-point numbers) is sufficient.

In the vast majority of cases, I would expect the situation to be even easier.
There would be a common denominator with its origin
in the original problem and no need for lcm.

I see two situations in which the all integer technique would not work.
In one, lcm returns a number too large to be reliably represented as a double.
In the other, the user wants to do his own scaling, if any.
Allowing an octal or hexadecimal radix would alleviate that somewhat.

-- 
Mike   address@hidden
"Nothing says it like words if you know how to use them."
                    --  the Professional Organization of English Majors