Skip to content

Commit 5c6877b

Browse files
jartjart
committed
Introduce support for trapping math
The feenableexcept() and fedisableexcept() APIs are now provided which let you detect when NaNs appear the moment it happens from anywhere in your program. Tests have also been added for the mission critical math functions expf() and erff(), whose perfect operation has been assured. See examples/trapping.c to see how to use this powerful functionality.
1 parent 403bc25 commit 5c6877b

File tree

13 files changed

+576
-2
lines changed

13 files changed

+576
-2
lines changed

‎examples/trapping.c

Lines changed: 142 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,142 @@
1+
#include <fenv.h>
2+
#include <math.h>
3+
#include <signal.h>
4+
#include <stdio.h>
5+
#include <string.h>
6+
#include <ucontext.h>
7+
#include <unistd.h>
8+
#include "libc/calls/struct/aarch64.internal.h"
9+
10+
/*
11+
Do you put lots of assert(!isnan(x)) in your code??
12+
Your microprocessor has a feature to automate this.
13+
14+
Uncaught SIGFPE (FPE_FLTINV)
15+
__math_invalidf at libc/tinymath/math_errf.c:88
16+
logf at libc/tinymath/logf.c:100
17+
main at examples/trapping.c:29
18+
cosmo at libc/runtime/cosmo.S:105
19+
_start at libc/crt/crt.S:116
20+
21+
This file shows how to use floating point exception
22+
trapping with Cosmopolitan Libc.
23+
*/
24+
25+
#define TRAPS (FE_INVALID | FE_DIVBYZERO | FE_OVERFLOW | FE_UNDERFLOW)
26+
27+
void spring_trap(int sig, siginfo_t *si, void *arg) {
28+
29+
// print signal safely
30+
const char *msg;
31+
int sic = si->si_code;
32+
if (sic == FPE_INTDIV)
33+
msg = "FPE_INTDIV: "; // integer divide by zero
34+
else if (sic == FPE_INTOVF)
35+
msg = "FPE_INTOVF: "; // integer overflow
36+
else if (sic == FPE_FLTDIV)
37+
msg = "FPE_FLTDIV: "; // floating point divide by zero
38+
else if (sic == FPE_FLTOVF)
39+
msg = "FPE_FLTOVF: "; // floating point overflow
40+
else if (sic == FPE_FLTUND)
41+
msg = "FPE_FLTUND: "; // floating point underflow
42+
else if (sic == FPE_FLTRES)
43+
msg = "FPE_FLTRES: "; // floating point inexact
44+
else if (sic == FPE_FLTINV)
45+
msg = "FPE_FLTINV: "; // invalid floating point operation
46+
else if (sic == FPE_FLTSUB)
47+
msg = "FPE_FLTSUB: "; // subscript out of range
48+
else
49+
msg = "SIGFPE: ";
50+
write(1, msg, strlen(msg));
51+
52+
// recover from trap so that execution may resume
53+
// without this the same signal will just keep getting raised
54+
ucontext_t *ctx = arg;
55+
#ifdef __x86_64__
56+
if (ctx->uc_mcontext.fpregs) {
57+
ctx->uc_mcontext.fpregs->mxcsr |= TRAPS << 7; // disable traps
58+
ctx->uc_mcontext.fpregs->mxcsr &= ~TRAPS; // clear cages
59+
return;
60+
}
61+
#elif defined(__aarch64__)
62+
struct _aarch64_ctx *ac;
63+
for (ac = (struct _aarch64_ctx *)ctx->uc_mcontext.__reserved; ac->magic;
64+
ac = (struct _aarch64_ctx *)((char *)ac + ac->size)) {
65+
if (ac->magic == FPSIMD_MAGIC) {
66+
struct fpsimd_context *sm = (struct fpsimd_context *)ac;
67+
sm->fpcr &= ~(TRAPS << 8); // disable traps
68+
sm->fpsr &= ~TRAPS; // clear cages
69+
return;
70+
}
71+
}
72+
#endif
73+
74+
// exit if we can't recover execution
75+
msg = "cannot recover from signal\n";
76+
write(1, msg, strlen(msg));
77+
_exit(1);
78+
}
79+
80+
void setup_trap(void) {
81+
struct sigaction sa;
82+
sigemptyset(&sa.sa_mask);
83+
sa.sa_flags = SA_SIGINFO;
84+
sa.sa_sigaction = spring_trap;
85+
sigaction(SIGFPE, &sa, 0);
86+
}
87+
88+
void activate_trap(void) {
89+
feclearexcept(TRAPS);
90+
if (feenableexcept(TRAPS)) {
91+
static bool once;
92+
if (!once) {
93+
fprintf(stderr, "warning: trapping math isn't supported on this cpu\n");
94+
once = true;
95+
}
96+
}
97+
}
98+
99+
float ident(float x) {
100+
return x;
101+
}
102+
float (*veil)(float) = ident;
103+
104+
int main(int argc, char *argv[]) {
105+
float x;
106+
setup_trap();
107+
108+
// test illegal math
109+
activate_trap();
110+
x = 0 / veil(0);
111+
printf("0/0 = %g\n", x);
112+
113+
// test divide by zero
114+
activate_trap();
115+
x = 1 / veil(0);
116+
printf("1/0 = %g\n", x);
117+
118+
// test divide by zero again
119+
activate_trap();
120+
x = -1 / veil(0);
121+
printf("-1/0 = %g\n", x);
122+
123+
// test domain error
124+
activate_trap();
125+
x = logf(veil(-1));
126+
printf("log(-1) = %g\n", x);
127+
128+
// test imaginary number
129+
activate_trap();
130+
x = sqrtf(veil(-1));
131+
printf("sqrt(-1) = %g\n", x);
132+
133+
// test overflow
134+
activate_trap();
135+
x = expf(veil(88.8));
136+
printf("expf(88.8) = %g\n", x);
137+
138+
// test underflow
139+
activate_trap();
140+
x = expf(veil(-104));
141+
printf("expf(-104) = %g\n", x);
142+
}

‎libc/calls/ucontext.h

Lines changed: 41 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -43,14 +43,55 @@ struct FpuStackEntry {
4343
};
4444

4545
struct thatispacked FpuState {
46+
47+
/* 8087 FPU Control Word
48+
IM: Invalid Operation ───────────────┐
49+
DM: Denormal Operand ───────────────┐│
50+
ZM: Zero Divide ───────────────────┐││
51+
OM: Overflow ─────────────────────┐│││
52+
UM: Underflow ───────────────────┐││││
53+
PM: Precision ──────────────────┐│││││
54+
PC: Precision Control ───────┐ ││││││
55+
{float,∅,double,long double}│ ││││││
56+
RC: Rounding Control ──────┐ │ ││││││
57+
{even, →-∞, →+∞, →0} │┌┤ ││││││
58+
┌┤││ ││││││
59+
d││││rr││││││
60+
0b0000001001111111 */
4661
uint16_t cwd;
62+
63+
/* 8087 FPU Status Word */
4764
uint16_t swd;
65+
4866
uint16_t ftw;
4967
uint16_t fop;
5068
uint64_t rip;
5169
uint64_t rdp;
70+
71+
/* SSE CONTROL AND STATUS REGISTER
72+
IE: Invalid Operation Flag ──────────────┐
73+
DE: Denormal Flag ──────────────────────┐│
74+
ZE: Divide-by-Zero Flag ───────────────┐││
75+
OE: Overflow Flag ────────────────────┐│││
76+
UE: Underflow Flag ──────────────────┐││││
77+
PE: Precision Flag ─────────────────┐│││││
78+
DAZ: Denormals Are Zeros ──────────┐││││││
79+
IM: Invalid Operation Mask ───────┐│││││││
80+
DM: Denormal Operation Mask ─────┐││││││││
81+
ZM: Divide-by-Zero Mask ────────┐│││││││││
82+
OM: Overflow Mask ─────────────┐││││││││││
83+
UM: Underflow Mask ───────────┐│││││││││││
84+
PM: Precision Mask ──────────┐││││││││││││
85+
RC: Rounding Control ───────┐│││││││││││││
86+
{even, →-∞, →+∞, →0} ││││││││││││││
87+
FTZ: Flush To Zero ───────┐ ││││││││││││││
88+
│┌┤│││││││││││││
89+
┌──────────────┐││││││││││││││││
90+
│ reserved │││││││││││││││││
91+
0b00000000000000000001111110000000 */
5292
uint32_t mxcsr;
5393
uint32_t mxcr_mask;
94+
5495
struct FpuStackEntry st[8];
5596
struct XmmRegister xmm[16];
5697
uint32_t __padding[24];

‎libc/cosmo.h

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -5,6 +5,7 @@ COSMOPOLITAN_C_START_
55
errno_t cosmo_once(_Atomic(uint32_t) *, void (*)(void));
66
int systemvpe(const char *, char *const[], char *const[]) libcesque;
77
char *GetProgramExecutableName(void);
8+
void unleaf(void);
89

910
COSMOPOLITAN_C_END_
1011
#endif /* COSMOPOLITAN_LIBC_COSMO_H_ */

‎libc/intrin/fedisableexcept.c

Lines changed: 92 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,92 @@
1+
/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
2+
│ vi: set et ft=c ts=2 sts=2 sw=2 fenc=utf-8 :vi │
3+
╞══════════════════════════════════════════════════════════════════════════════╡
4+
│ Copyright 2024 Justine Alexandra Roberts Tunney │
5+
│ │
6+
│ Permission to use, copy, modify, and/or distribute this software for │
7+
│ any purpose with or without fee is hereby granted, provided that the │
8+
│ above copyright notice and this permission notice appear in all copies. │
9+
│ │
10+
│ THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL │
11+
│ WARRANTIES WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED │
12+
│ WARRANTIES OF MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE │
13+
│ AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL │
14+
│ DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR │
15+
│ PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER │
16+
│ TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR │
17+
│ PERFORMANCE OF THIS SOFTWARE. │
18+
╚─────────────────────────────────────────────────────────────────────────────*/
19+
#include "libc/runtime/fenv.h"
20+
21+
/**
22+
* Disables floating point exception trapping, e.g.
23+
*
24+
* feenableexcept(FE_INVALID | FE_DIVBYZERO |
25+
* FE_OVERFLOW | FE_UNDERFLOW);
26+
*
27+
* When trapping is enabled, something should handle SIGFPE. Calling
28+
* ShowCrashReports() at startup will install a generic handler with
29+
* backtraces and the symbol of the `si->si_code` which UNIX defines
30+
*
31+
* - `FPE_INTOVF`: integer overflow
32+
* - `FPE_INTDIV`: integer divide by zero
33+
* - `FPE_FLTDIV`: floating point divide by zero
34+
* - `FPE_FLTOVF`: floating point overflow
35+
* - `FPE_FLTUND`: floating point underflow
36+
* - `FPE_FLTRES`: floating point inexact
37+
* - `FPE_FLTINV`: invalid floating point operation
38+
* - `FPE_FLTSUB`: subscript out of range
39+
*
40+
* It's important to not use the `-ffast-math` or `-Ofast` flags when
41+
* compiling code that needs to be debugged. Using `-fsignaling-nans`
42+
* will also help, since GCC doesn't enable that by default.
43+
*
44+
* @param excepts may bitwise-or the following:
45+
* - `FE_INVALID`
46+
* - `FE_DIVBYZERO`
47+
* - `FE_OVERFLOW`
48+
* - `FE_UNDERFLOW`
49+
* - `FE_INEXACT`
50+
* - `FE_ALL_EXCEPT` (all of the above)
51+
* @see fetestexcept() if you don't want to deal with signals
52+
* @see feenableexcept() to turn it on in the first place
53+
*/
54+
int fedisableexcept(int excepts) {
55+
56+
// limit to what we know
57+
excepts &= FE_ALL_EXCEPT;
58+
59+
#ifdef __x86_64__
60+
61+
#ifndef NOX87
62+
// configure 8087 fpu control word
63+
// setting the bits enables suppression
64+
unsigned short x87cw;
65+
asm("fstcw\t%0" : "=m"(x87cw));
66+
x87cw |= excepts;
67+
asm("fldcw\t%0" : /* no inputs */ : "m"(x87cw));
68+
#endif
69+
70+
// configure modern sse control word
71+
// setting the bits enables suppression
72+
unsigned mxcsr;
73+
asm("stmxcsr\t%0" : "=m"(mxcsr));
74+
mxcsr |= excepts << 7;
75+
asm("ldmxcsr\t%0" : /* no inputs */ : "m"(mxcsr));
76+
77+
return 0;
78+
79+
#elif defined(__aarch64__)
80+
81+
unsigned fpcr;
82+
unsigned fpcr2;
83+
fpcr = __builtin_aarch64_get_fpcr();
84+
fpcr2 = fpcr & ~(excepts << 8);
85+
if (fpcr != fpcr2)
86+
__builtin_aarch64_set_fpcr(fpcr2);
87+
return (fpcr >> 8) & FE_ALL_EXCEPT;
88+
89+
#else
90+
return -1;
91+
#endif
92+
}

‎libc/intrin/feenableexcept.c

Lines changed: 98 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,98 @@
1+
/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
2+
│ vi: set et ft=c ts=2 sts=2 sw=2 fenc=utf-8 :vi │
3+
╞══════════════════════════════════════════════════════════════════════════════╡
4+
│ Copyright 2024 Justine Alexandra Roberts Tunney │
5+
│ │
6+
│ Permission to use, copy, modify, and/or distribute this software for │
7+
│ any purpose with or without fee is hereby granted, provided that the │
8+
│ above copyright notice and this permission notice appear in all copies. │
9+
│ │
10+
│ THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL │
11+
│ WARRANTIES WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED │
12+
│ WARRANTIES OF MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE │
13+
│ AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL │
14+
│ DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR │
15+
│ PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER │
16+
│ TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR │
17+
│ PERFORMANCE OF THIS SOFTWARE. │
18+
╚─────────────────────────────────────────────────────────────────────────────*/
19+
#include "libc/runtime/fenv.h"
20+
21+
/**
22+
* Enables floating point exception trapping, e.g.
23+
*
24+
* feenableexcept(FE_INVALID | FE_DIVBYZERO |
25+
* FE_OVERFLOW | FE_UNDERFLOW);
26+
*
27+
* When trapping is enabled, something should handle SIGFPE. Calling
28+
* ShowCrashReports() at startup will install a generic handler with
29+
* backtraces and the symbol of the `si->si_code` which UNIX defines
30+
*
31+
* - `FPE_INTOVF`: integer overflow
32+
* - `FPE_INTDIV`: integer divide by zero
33+
* - `FPE_FLTDIV`: floating point divide by zero
34+
* - `FPE_FLTOVF`: floating point overflow
35+
* - `FPE_FLTUND`: floating point underflow
36+
* - `FPE_FLTRES`: floating point inexact
37+
* - `FPE_FLTINV`: invalid floating point operation
38+
* - `FPE_FLTSUB`: subscript out of range
39+
*
40+
* It's important to not use the `-ffast-math` or `-Ofast` flags when
41+
* compiling code that needs to be debugged. Using `-fsignaling-nans`
42+
* will also help, since GCC doesn't enable that by default.
43+
*
44+
* @param excepts may bitwise-or the following:
45+
* - `FE_INVALID`
46+
* - `FE_DIVBYZERO`
47+
* - `FE_OVERFLOW`
48+
* - `FE_UNDERFLOW`
49+
* - `FE_INEXACT`
50+
* - `FE_ALL_EXCEPT` (all of the above)
51+
* @see fetestexcept() if you don't want to deal with signals
52+
* @see fedisableexcept() to turn it back off again
53+
*/
54+
int feenableexcept(int excepts) {
55+
56+
// limit to what we know
57+
excepts &= FE_ALL_EXCEPT;
58+
59+
#ifdef __x86_64__
60+
61+
#ifndef NOX87
62+
// configure 8087 fpu control word
63+
// celaring the bits disables suppression
64+
unsigned short x87cw;
65+
asm("fstcw\t%0" : "=m"(x87cw));
66+
x87cw &= ~excepts;
67+
asm("fldcw\t%0" : /* no inputs */ : "m"(x87cw));
68+
#endif
69+
70+
// configure modern sse control word
71+
// clearing the bits disables suppression
72+
unsigned mxcsr;
73+
asm("stmxcsr\t%0" : "=m"(mxcsr));
74+
mxcsr &= ~(excepts << 7);
75+
asm("ldmxcsr\t%0" : /* no inputs */ : "m"(mxcsr));
76+
77+
return 0;
78+
79+
#elif defined(__aarch64__)
80+
81+
unsigned fpcr;
82+
unsigned fpcr2;
83+
unsigned updated_fpcr;
84+
fpcr = __builtin_aarch64_get_fpcr();
85+
fpcr2 = fpcr | (excepts << 8);
86+
if (fpcr != fpcr2) {
87+
__builtin_aarch64_set_fpcr(fpcr2);
88+
// floating point exception trapping is optional in aarch64
89+
updated_fpcr = __builtin_aarch64_get_fpsr();
90+
if (fpcr2 & ~updated_fpcr)
91+
return -1;
92+
}
93+
return (fpcr >> 8) & FE_ALL_EXCEPT;
94+
95+
#else
96+
return -1;
97+
#endif
98+
}

0 commit comments

Comments
 (0)