MATLAB Central Newsreader - bug in rem

Re: bug in rem

roberson@ibd.nrc-cnrc.gc.ca (Walter Roberson) — Tue, 31 Jul 2007 19:55:35 -0400

In article <46ae685a$1@news.eftel.com.au>, ross <ross@nospam.com> wrote:

>rem and mod should not, by definition, return values that exceed the magnitude
>of the modulus or radix, and the software should be written to reflect that.
>This is an "absolute" consideration. IMO the software is broken if it doesn't
>conform to the expected mathematical behaviour, when that is possible.

The documentation for rem and mod do not indicate that the functions
are the abstract algebraic functions of mathematics. Instead,
the documentation gives explicit formulae for how the result is
produced, and those explicit formulae lead to the results you observed.

This is no different than understanding that the '+' operator
does not implement the abstract algebraic addition of mathematics.
Operators in a computer language are defined by their documented
behaviour, not by the behaviour of mathematical constructs
with similar names. Sometimes the difference is startling,
but I don't think it rises to the level of,
'This is an "absolute" consideration.'.

If you need a language with complete abstract mathematical fidelity
or one in which all computations are arranged to maximize
retention of precision, then Matlab is not the right language
for your purposes.
--
Is there any thing whereof it may be said, See, this is new? It hath
been already of old time, which was before us. -- Ecclesiastes

Re: bug in rem

ross — Tue, 31 Jul 2007 06:38:18 -0400

"us " <us@neurol.unizh.ch> wrote:
>ross:
><SNIP FP headache evergreen...
>
>> x = 1.5511623754564467e105;
>> r = rem(x,10)
>> r = 2.54629497041811e+089
>
>...should not come as a surprise, because your <x> is represented with the
>uncertainty of
>
> eps(x)
>
>that happens to be ---
>
> 2.546294970418108e+089
>
>us

Look, the round-off error isn't the issue. You will note that I said "whatever
the cause". That phrase was meant to indicate that I know that the value
returned is due to roundoff error, but the real problem is different. In fact
the real problem is that a rem routine should never return that value. This is
not an evergreen newbie problem, it's a comment on the inadequacy of the
programming.

rem and mod should not, by definition, return values that exceed the magnitude
of the modulus or radix, and the software should be written to reflect that.
This is an "absolute" consideration. IMO the software is broken if it doesn't
conform to the expected mathematical behaviour, when that is possible. In other
words, if I want to express pi to 100 sig figs in a Matlab double, it can't be
done. However it is entirely possible to ensure that a remainder doesn't
exceed the divisor, and it absolutely should be done.

Secondly, I believe that there is a more "relative" consideration, which is
that the numerical routines should protect the user from the effects of
round-off error where possible. None of us disagree that one should understand
the effects of round off error in order to program in an informed way.
However, it is absurd to have to put all the error trapping code into every
application routine. It should be in the basic numerical functions. There are
many applications where certain inputs might generate high disparaties in the
values at some point in the routine. The routine should not crash at that
point. It's true that if one does a careful error analysis one might expect
those situations. But this is not the 1950s anymore guys - the user can expect
to have routines that dealt with out-of-range cases gracefully, not by
returning absurd values. If you don't include sensible error trapping in the
routines, you force the user to be watching the details all the time instead of
concentrating on the main objectives. It's the old mainframe mentality where
the programming issues were allowed to take precedence over what was important
to the user, which was why programming was such a black art then. It's very
important to get away from this oldfashioned point of view. (And for God's
sake, before you send a bunch of irrelevant flames -- I am not arguing that an
understanding of round-off error is oldfashioned).

Ross

bug in rem

ross — Tue, 31 Jul 2007 04:23:58 -0400

Can you confirm the following (R2007a):

x = 1.5511623754564467e105;
r = rem(x,10)

r =

2.54629497041811e+089

If you can, there is a bug in rem( ), since the result should never be > 10,
whatever the cause.

ross

Re: bug in rem

David Wilson — Tue, 31 Jul 2007 01:08:27 -0400

ross <ross@nospam.com> wrote in message <46ae685a$1@news.eftel.com.au>...
> "us " <us@neurol.unizh.ch> wrote:
> >ross:
> ><SNIP FP headache evergreen...
> >
> >> x = 1.5511623754564467e105;
> >> r = rem(x,10)
> >> r = 2.54629497041811e+089
> >
> >...should not come as a surprise, because your <x> is represented with the
> >uncertainty of
> >
> > eps(x)
> >
> >that happens to be ---
> >
> > 2.546294970418108e+089
> >
> >us
>
>
> Look, the round-off error isn't the issue. You will note that I said "whatever
> the cause". That phrase was meant to indicate that I know that the value
> returned is due to roundoff error, but the real problem is different. In fact
> the real problem is that a rem routine should never return that value. This is
> not an evergreen newbie problem, it's a comment on the inadequacy of the
> programming.
>
> rem and mod should not, by definition, return values that exceed the magnitude
> of the modulus or radix, and the software should be written to reflect that.
> This is an "absolute" consideration. IMO the software is broken if it doesn't
> conform to the expected mathematical behaviour, when that is possible. In other
> words, if I want to express pi to 100 sig figs in a Matlab double, it can't be
> done. However it is entirely possible to ensure that a remainder doesn't
> exceed the divisor, and it absolutely should be done.
>
> Secondly, I believe that there is a more "relative" consideration, which is
> that the numerical routines should protect the user from the effects of
> round-off error where possible. None of us disagree that one should understand
> the effects of round off error in order to program in an informed way.
> However, it is absurd to have to put all the error trapping code into every
> application routine. It should be in the basic numerical functions. There are
> many applications where certain inputs might generate high disparaties in the
> values at some point in the routine. The routine should not crash at that
> point. It's true that if one does a careful error analysis one might expect
> those situations. But this is not the 1950s anymore guys - the user can expect
> to have routines that dealt with out-of-range cases gracefully, not by
> returning absurd values. If you don't include sensible error trapping in the
> routines, you force the user to be watching the details all the time instead of
> concentrating on the main objectives. It's the old mainframe mentality where
> the programming issues were allowed to take precedence over what was important
> to the user, which was why programming was such a black art then. It's very
> important to get away from this oldfashioned point of view. (And for God's
> sake, before you send a bunch of irrelevant flames -- I am not arguing that an
> understanding of round-off error is oldfashioned).
>
> Ross
>

While your argument is well constructed regarding the fact that rem(a,b) should never return a value outside [0,b-1],
the problem remains that what, in your opinion, should rem then return in your somewhat extreme example?

We should not really return NaN or Inf simply because we can't compute it (& it is already used), and we cannot really return anything in the range [0,b-1] since that seems to indicate that that is infact the answer. We could throw up a warning I suppose, but that's not done in other similar overflow cases.

Re: bug in rem

ellieandrogerxyzzy@mindspring.com.invalid (Roger Stafford) — Tue, 31 Jul 2007 00:04:38 -0400

In article <46ae685a$1@news.eftel.com.au>, ross <ross@nospam.com> wrote:

> Look, the round-off error isn't the issue. You will note that I said
"whatever
> the cause". That phrase was meant to indicate that I know that the value
> returned is due to roundoff error, but the real problem is different.
In fact
> the real problem is that a rem routine should never return that value.
This is
> not an evergreen newbie problem, it's a comment on the inadequacy of the
> programming.
>
> rem and mod should not, by definition, return values that exceed the
magnitude
> of the modulus or radix, and the software should be written to reflect that.
> This is an "absolute" consideration. IMO the software is broken if it
doesn't
> conform to the expected mathematical behaviour, when that is possible.
In other
> words, if I want to express pi to 100 sig figs in a Matlab double, it
can't be
> done. However it is entirely possible to ensure that a remainder doesn't
> exceed the divisor, and it absolutely should be done.
>
> Secondly, I believe that there is a more "relative" consideration, which is
> that the numerical routines should protect the user from the effects of
> round-off error where possible. None of us disagree that one should
understand
> the effects of round off error in order to program in an informed way.
> However, it is absurd to have to put all the error trapping code into every
> application routine. It should be in the basic numerical functions.
There are
> many applications where certain inputs might generate high disparaties in the
> values at some point in the routine. The routine should not crash at that
> point. It's true that if one does a careful error analysis one might expect
> those situations. But this is not the 1950s anymore guys - the user can
expect
> to have routines that dealt with out-of-range cases gracefully, not by
> returning absurd values. If you don't include sensible error trapping in the
> routines, you force the user to be watching the details all the time
instead of
> concentrating on the main objectives. It's the old mainframe mentality where
> the programming issues were allowed to take precedence over what was
important
> to the user, which was why programming was such a black art then. It's very
> important to get away from this oldfashioned point of view. (And for God's
> sake, before you send a bunch of irrelevant flames -- I am not arguing
that an
> understanding of round-off error is oldfashioned).
>
> Ross
--------------------
There is an old computer adage that is pertinent here, Ross. "Garbage
In, Garbage Out". And, given the 53-bit accuracy limitation of IEEE
binary double precision floating points numbers, "garbage" is precisely
what you have put into the 'rem' function input when you write the absurd
expression rem(1.5511623754564467e105,10)! Even with the very first step
of finding the nearest representable number to that decimal string, an
integer in this case, the correct answer for the remainder with ten as
divisor will have been irretrievably lost. IEEE 754 binary numbers almost
certainly cannot represent 1.5511623754564467e105 exactly even if it is an
integer with 89 trailing zeros. Matlab has to select the nearest integer
it can realize and the difference in this case will be astronomically
large, which will of course change the remainder completely. Would you
rather the 'rem' function issue the message, "You've to be kidding,
fella!"? In my opinion it would almost be justified. The lesson to be
learned here is that there is no substitute for using the little gray
cells when it comes to writing good programs. Matlab, and indeed all such
numerical computation programs, simply won't do all your thinking for you
- perhaps by the next millennium, but certainly not yet. You as a
programmer have to be constantly aware of such limitations when dealing
with a finite set of numbers and attempting to approximate the
mathematical ideal of the infinite continuum of real numbers. If you do
otherwise, you put yourself in the same class as people who ask
indignantly why 100*.07==7 (or numerous similar equalities) come out false
in matlab.

Roger Stafford

Re: bug in rem

us — Mon, 30 Jul 2007 20:44:26 -0400

ross:
<SNIP FP headache evergreen...

> x = 1.5511623754564467e105;
> r = rem(x,10)
> r = 2.54629497041811e+089

...should not come as a surprise, because your <x> is represented with the uncertainty of

eps(x)

that happens to be ---

2.546294970418108e+089

us

Re: bug in rem

roberson@ibd.nrc-cnrc.gc.ca (Walter Roberson) — Mon, 30 Jul 2007 20:35:18 -0400

In article <46ae48de$1@news.eftel.com.au>, ross <ross@nospam.com> wrote:
>Can you confirm the following (R2007a):

>x = 1.5511623754564467e105;
>r = rem(x,10)

>r =
>
> 2.54629497041811e+089

>If you can, there is a bug in rem( ), since the result should never be > 10,
>whatever the cause.

>> help rem
REM Remainder after division.
    REM(x,y) is x - n.*y where n = fix(x./y) if y ~= 0. If y is not an
    integer and the quotient x./y is within roundoff error of an integer,
    then n is that integer. The inputs x and y must be real arrays of the
    same size, or real scalars.

If the inputs are floating point, then you are dealing with the
limits of floating point representation. There simply aren't enough
bits available in a Matlab float to calculate the answer properly.

1.5511623754564467e105 does not have enough bits to represent
right down to the units position. 1.5511623754564467e105 is
really 1.5511623754564467e105 +/- one part in 2^(-53)
which is 1.5511623754564467e105 +/- 1.722136184164e+89 or so.
The answer you got back was within margin of error for the
question.
--
  "law -- it's a commodity"
                         -- Andrew Ryan (The Globe and Mail, 2005/11/26)