[GRASS-dev] r.mapcal rand() strangeness

[Glynn Clements]

Maciej Sieczka wrote:

> >>>> $ r.mapcalc 'map=rand(-2147483648,2147483647)'
> 
> >>>     res[i] = (lo == hi) ? lo : lo + x % (hi - lo);
> >>>
> >>> The expression (hi - lo) is overflowing the range of a signed integer.
> 
> > Are all GRASS modules limited to Int32, ie. -2147483648 to 214748364? I 
> > can see r.mapcalc silently forces 214748364 for anything bigger...
> 
> This is not entirely correct. More findings:
> 
> 
> 
> The following silently forces 2147483647 on anything bigger:
> 
> $ r.mapcalc 'map=9999999999999999999'
> $ r.info -rt map
> min=2147483647
> max=2147483647
> datatype=CELL

That is down to the atoi() function used to parse integers.

I could change it to use strtol(), which sets errno to ERANGE on
overflow.

> This one too, though I'm trying to enforce double:
> 
> $ r.mapcalc 'map=double(9999999999999999999)'
> $ r.info -rt map
> min=2147483647.000000
> max=2147483647.000000
> datatype=DCELL

The over-large integer constant is truncated by atoi(), then the
truncated constant is converted to a double.

> Trying the same by adding a decimal place separator, results in 
> something bizzare to me:
> 
> $ r.mapcalc 'map=99999999999999999999999.0'
> $ r.info -rt map
> min=99999999999999991611392.000000
> max=99999999999999991611392.000000
> datatype=DCELL

The presence of a dot causes the value to be parsed as a
floating-point value using atof().

A double has a precision of ~16 decimal digits, which matches what you
see above.

> $ r.mapcalc 'map=9999999999999999999.0'
> $ r.info -rt map
> min=10000000000000000000.000000
> max=10000000000000000000.000000
> datatype=DCELL

Ditto.

> At the edges of Float32 range other strange things happen:
> 
> Shouldn't the below be rather 
> min,max=340282000000000000000000000000000000000.0?:
> 
> $ r.mapcalc 'map=3.40282e38'
> $ r.info -rt map
> min=340282000000000014192072600942972764160.000000
> max=340282000000000014192072600942972764160.000000
> datatype=DCELL
> 
> The the one below should be -340282000000000000000000000000000000000.0 I 
> guess:
> 
> $ r.mapcalc 'map=-3.40282e38'
> $ r.info -rt map
> min=-340282000000000014192072600942972764160.000000
> max=-340282000000000014192072600942972764160.000000
> datatype=DCEL

Again, you're running into the limits of the precision of the double
type.

Also, bear in mind that the conversion from a binary FP value to a
string of decimal digits is being performed by r.info; r.mapcalc has
no control over that part of the process.

Essentially, the %f format behaves badly for large numbers. If you use
%e or %g, it will use exponential notation, with a specific relative
precision rather than a specific absolute precision.

> At bigger Float64 values r.mapcalc crashes:
> 
> $r.mapcalc 'map=3.40282e110'
> 
> *** stack smashing detected ***: r.mapcalc terminated
> Aborted (core dumped)
> 
> $ r.mapcalc 'map=1.79769e+308'
>   100%
> *** stack smashing detected ***: r.mapcalc terminated
> Aborted (core dumped

Right. It's sprintf()ing the value into a 64-byte buffer using %f.

I'll change it to use %g:

--- raster/r.mapcalc/expression.c	(revision 30318)
+++ raster/r.mapcalc/expression.c	(working copy)
@@ -312,7 +312,7 @@
 	if (e->res_type == CELL_TYPE)
 		sprintf(buff, "%d", e->data.con.ival);
 	else
-		sprintf(buff, "%f", e->data.con.fval);
+		sprintf(buff, "%.8g", e->data.con.fval);
 
 	return strdup(buff);
 }

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