[GRASS-dev] Overflow warning in r.mapcalc calculation

[Glynn Clements]

Rainer M. Krug wrote:

> >> Sounds like a sensible approach without adding to many new
> >> functions. But I would actually split the two, i.e. have two more
> >> arguments, where one specifies the type, 
> >> and the other one the number of decimals to round to, i.e.
> >> 
> >> round(x, 0, "I") would be the default, rounding to whole number and
> >> return an integer
> >
> > That interface isn't possible. r.mapcalc doesn't have a string type,
> > and there's no way that a function's return type can depend upon the
> > value of a parameter.
> 
> 
> In this case I would suggest different functions, e.g.
> 
> round(x) - as it is now
> 
> roundI(x, dec) - round to integer, dec>0 is ignored, dec<0 as below
> 
> roundD(x, dec) - round to double, where dec is the number of decimals it
> should be rounded to
> 
> roundF(x, dec) - round to float, where dec is the number of decimals it
> should be rounded to
> 
> if dec is positive, it is the number of decimals, if negative the
> ... (no idea what it is called, but if it is -1, 111 is rounded to 110,
> -2 111 is rounded to 100)
> 
> default from dec would be 0.
> 
> These different functions would also have the advantage of probably
> being faster then one single function.

1. I wouldn't bother with separate functions, but instead use the type
of the second argument to determine the return type, i.e.

	round(x, 1)    -> CELL
	round(x, 1.0f) -> FCELL
	round(x, 1.0)  -> DCELL

The value of the second argument can be ignored, or it can be used for
other purposes.

2. Rounding to a given number of decimals is unnecessarily limiting.
The algorithm for generalised rounding is essentially:

	roundTo(x, k) = round(x / k) * k.

Rounding to N decimal places is just a case of using k=1/10^N. If you
allow k to be specified directly, then you can round to arbitrary
steps (e.g. k=5 would round to the nearest multiple of 5, etc).

3. However: there's a slight problem with doing it that way: 0.1 isn't
exactly representable in binary, so e.g. x/0.1 isn't equal to x*10; it
would be more accurate to use:

	roundTo(x, k) = round(x * k) / k

where k is the reciprocal of the step, so k=10^N to round to N decimal
places (or k=2 to round to 1/2).

The downside is that the interface is less useful if you want to round
to something other than a fixed number of decimal places. E.g. if you
wanted to round to the nearest multiple of 45 degrees, you'd need to
use k=1.0/45, which isn't exactly representable.

-- 
Glynn Clements <glynn at gclements.plus.com>