[wenzeslaus]

Integer modulo by zero (r.mapcalc "x = 5 % 0") hangs r.mapcalc. This adds a zero-divisor check and removes the SIGFPE handler and its per-cell flag. With floating point traps off by default, the handler only ever ran for this integer case, where it caused the hang, while resetting the shared flag for every cell slowed down computation with multiple threads.

Fun fact: I used grass-legacy repo to investigate this, so there is some code history junkie stuff to ponder in the details block.

Edit: To clarify the history (thanks @nilason): the signal handler is older than the missing check. It dates back to 1990 and was meant to catch floating point errors as they happen during the computation and turn the affected cell into NULL, instead of relying on each function to check its inputs. Integer % had its own zero check until the 2002 rewrite dropped it. Today, floating point errors don't raise the signal at all: they produce infinity or NaN, the functions check their inputs where they want NULL, and NaN is NULL in GRASS. So the handler wasn't needed for floating point, and this integer case was the only thing it caught.

Details: cause, behavior, benchmarks, and tests (AI generated)

Cause: The integer % did not check for a zero divisor. On x86, integer division by zero traps, the SIGFPE handler returns, and the CPU executes the same division again, so r.mapcalc spins forever. Integer / already had the check. Integer % checked for a zero divisor until the 2002 rewrite (de8f7a4fd9 in grass-legacy) dropped the check.

Flags: Floating point exceptions are masked by default, so the handler, which dates from 1990 (in its sigaction() form from 2002), never runs for floating point operations, and on Windows it is not installed at all. Since the OpenMP support, all threads reset the one shared volatile flag to 0 before every cell and read it after.

Division by zero, domain errors, and overflow without the flag: Floating point exceptions are masked by default (IEEE 754), so an invalid operation, division by zero, or overflow does not raise SIGFPE; it returns NaN or infinity and the computation continues. What a cell gets is the same before and after this PR, except integer %:

• Checks before computing (GRASS code): zero divisor in integer and floating point / and now in integer %, sqrt of a negative, log of a non-positive value or base, and pow/exp of a negative base with a non-integer exponent.
• NaN results (C library): acos/asin outside [-1, 1], floating point % by zero (fmod), sin/cos/tan of infinity, and 0/0 such as log(1, 1). GRASS treats any NaN as NULL (Rast_is_d_null_value() is x != x): writing a row records it in the null file, and reading gives NULL back.
• Infinity is a value, not NULL: overflow such as exp(1000) or pow(10.0, 400), and log(x, 1) for x other than 1, give inf or -inf. GRASS writes and reads infinity like any other value, which the new test_overflow_is_infinity checks on all CI platforms.

The log(x, 1) case does not match the documentation, before or after this PR: the manual says that invalid arguments should give NULL, and the log() source comment says that a base of 1.0 gives NULL, but the code gives inf or -inf (NULL only for log(1, 1)). This PR keeps that behavior and the comment's claim, and has no test for it; it is left for a separate fix.

API: This removes pre_exec(), post_exec(), floating_point_exception, and floating_point_exception_occurred from the C API. grass/calc.h and grass/defs/calc.h are installed with both build systems and libgrass_calc exported the symbols, so this is a change in the public API, although nothing marks the library as internal or external. In this repository, only r.mapcalc and lib/calc include grass/calc.h, grass-addons has no users (checked in a clone from July 2026), and there are no Python bindings for the library.

History: From grass-legacy and this repository's history:

• 1990, GRASS 4 (src/Gmapcalc/mapcalc/fpe.c): the SIGFPE handler sets the two flags and re-installs itself with signal().
• 1999, the oldest GRASS 5 revision in grass-legacy (src/raster/r.mapcalc/mapcalc/math.c): floating point operations clear the flag, compute, and set NULL if the flag was set. Integer modulus checks *b == 0 before %=, so it never relied on the flag. A floating point error printed a note and gave exit code 1, with the map still written (exit(allok), "fixed 11/99").
• 2002, "New version of r.mapcalc" (de8f7a4fd9): volatile flags, sigaction() in pre_exec(), and the same per-cell pattern in the x*.c functions; integer % lost the zero check.
• 2016 (87274dfb09): the functions move to lib/calc.
• 2025, OpenMP support (d5c51924ea): threads compute rows in parallel and share the flag.

The per-cell pattern can only work if the error is reported after the operation completes, so that execution reaches the flag check after the handler returns. An x86 integer division by zero is a fault that re-runs the instruction instead, so no check placed after the division is ever reached.

Uncertainties: I did not verify which platforms GRASS ran on between 1990 and 2002, or whether they trapped floating point errors by default. My understanding is that some architectures reported arithmetic errors as traps after the instruction (e.g., VAX), and that x87 reports unmasked exceptions at the next floating point instruction; whether that made the flags work in practice for GRASS is not known. Exit codes before November 1999 were not checked.

Behavior: Old is current main, likely intended is what the SIGFPE handler and flag checks were written for.

case	new (this PR)	old	likely intended
integer % by zero (5 % 0)	NULL	hangs on x86 Linux, see Platforms	NULL (manual; integer / checks for zero)
float % by zero (5.0 % 0.0)	NULL (NaN)	NULL (NaN)	NULL
domain errors (acos(2), sqrt(-1), log(0), pow(-8.0, 0.5))	NULL (explicit check or NaN)	NULL (explicit check or NaN)	NULL
overflow (exp(1000), pow(10.0, 400))	infinity	infinity	NULL
"Floating point error(s) occurred" warning and exit code 1	removed	never (only after a trap, and the % trap hangs)	map written with NULL cells, then warning and exit 1
per-cell flag with multiple threads	none	one shared flag reset by all threads for every cell	single-threaded code, one flag

Platforms: Only the hang on x86-64 Linux was reproduced. The other platforms are expected from the code and the CPU behavior, not tested: x86 macOS installs the same handler, so it should hang too; Windows has no handler (#ifndef _WIN32), so the trap should crash r.mapcalc; on ARM64 (Apple Silicon, Linux aarch64), integer division by zero does not trap and gives 0, so 5 % 0 should silently give 5. In C, it is undefined behavior anyway. With the zero check, the result is NULL everywhere, which the new tests check on all CI platforms.

Benchmark: f_div on 3000x20000 cells called directly from OpenMP threads, median of 7 runs, i7-10510U, -O2:

threads	before	after
1	0.133 s	0.096 s
4	0.229 s	0.025 s
8	0.200 s	0.025 s

sqrt behaves the same.

Whole r.mapcalc, 6000x6000 cells, x = sqrt(a / b) + log(a / b) + a / b / b + ... with 25 such terms, one run each on the same machine:

nprocs	before	after
1	80.6 s	78.7 s
2	45.7 s	44.3 s
3	35.6 s	32.8 s
4	34.7 s	29.1 s
5	39.0 s	31.3 s
6	35.4 s	28.3 s
7	35.4 s	26.2 s
8	39.2 s	29.2 s

Before, the time stops improving at 3 threads and stays at 35-39 s; after, it keeps improving to 26 s at 7 threads, and is at least 16% faster from 4 threads up. The machine has 4 physical cores, so 5-8 threads use hyper-threading.

Benchmark at -O0 and -O2: Whole r.mapcalc on a Xeon W-2295 @ 3.00GHz (18 cores, 36 threads, shared with other users), 3000x3000 cells, x = sqrt(a / b) + log(a / b) + a / b / b + sin(a) * cos(b) + exp(a, 0.5) + sqrt(b / a) + log(b / a) with a and b random double maps, i.e., 13 operations per cell which reset the flag before this PR. Wall time in seconds, median (min-max) of 5 runs, the builds interleaved within each run (figure above):

	nprocs 1	nprocs 2	nprocs 4	nprocs 8	nprocs 16
main, -O2	2.01 (2.00-2.03)	1.15 (1.10-1.18)	0.70 (0.69-0.71)	0.55 (0.54-0.56)	0.58 (0.57-0.61)
this PR, -O2	2.02 (1.99-2.03)	1.08 (1.07-1.09)	0.59 (0.56-0.62)	0.36 (0.35-0.37)	0.39 (0.38-0.39)
main, -O0	4.57 (4.53-4.65)	2.42 (2.39-2.46)	1.31 (1.29-1.33)	0.82 (0.81-0.83)	0.77 (0.74-0.79)
this PR, -O0	4.54 (4.50-4.63)	2.34 (2.30-2.42)	1.24 (1.20-1.25)	0.72 (0.68-0.75)	0.52 (0.51-0.53)

At 1 thread the builds are equal within the min-max ranges at both optimization levels. At -O2, main stops improving at 8 threads and is no faster at 16, while this PR reaches 0.36 s at 8 threads and stays there at 16; this PR is 16% faster at 4 threads, 35% at 8, and 33% at 16. At -O0, the gap appears later, 13% at 8 threads and 32% at 16, because the unoptimized arithmetic hides the contention on the shared flag until more threads write it. The result maps are identical cell by cell across builds and thread counts. Beyond 18 threads the machine shares cores, so 16 threads is the maximum used.

Alternatives measured: Making the flag thread-local scaled with threads, but was 1.6x slower on one thread than the current code. Detecting errors with <fenv.h> per cell was 3-40x slower and would turn overflow into NULL with a non-zero exit.

Tests: The new pytest file checks that division and modulo by zero give NULL for all three types, and so do domain errors of each function whose flag check was removed and infinite input to sin, cos, and tan. It also checks correct values of those functions, overflow to infinity, and a zero divisor in only some cells, all with 1 and 4 threads. With the old library, exactly the four % cases fail by timing out.

The hang was reproduced, the alternatives were benchmarked, and the fix and tests were prepared with AI assistance (Claude Code with Opus 5.5). The -O0 and -O2 benchmark was run with Claude Code with Fable 5.1.

Suggested squash commit message

r.mapcalc: Fix modulo by zero hang and speed up threads

Check for a zero divisor in the integer modulo, so that 5 % 0 gives
NULL as documented. It hung r.mapcalc on x86 Linux because the SIGFPE
handler returned and the division was executed again.

Remove the SIGFPE handler and the floating_point_exception flags.
Floating point exceptions do not trap by default, so the handler never
ran for floating point operations and the flags never became 1. Still,
every thread reset the shared flag to 0 for every cell, which made the
computation scale poorly with threads. Domain errors still give NaN,
which is NULL, and overflow still gives infinity.

Remove pre_exec(), post_exec(), and the two floating_point_exception
variables from the installed grass/calc.h and from libgrass_calc, so
technically a change to the API. Only r.mapcalc used them.

Test division and modulo by zero, domain errors and infinite input of
the functions which had the flag check, their regular results, and
overflow, all with 1 and 4 threads.

The fix, tests, and benchmarks were prepared with AI assistance (Claude
Code with Opus 5.5).

[wenzeslaus]

I added explanation and a more detailed benchmark showing the speedup. I'm resolving the conversation and merging.

A Claude Fable agent discovered that for the log function, the code, comments and user documentation do not match. It is unrelated to the goal here, so I'm leaving this for a separate PR.