x87 Floating Point Unit

On x86, floating-point operations require a bit of compiler and OS support. This is primarily due to the CPU using extra registers for floating-point operations. These registers have to be saved and restored for context switches. We focus on the x87 Floating Point Unit, as it has a relatively small register file (108 byte). This FPU is not entirely accurate and some operations (srqt, sin, cos) might be approximated.

1. Compiler Flags

The following flags tell the compiler that we want to use the x87 FPU.

C++: tools/build.mk

# Replace this variable with the following line
CXXFLAGS_NOFPU = -mno-mmx -mno-sse -mfpmath=387 -mpc32 -m80387 -mno-fp-ret-in-387

2. Thread Context

Additionally, we have to ensure that we do not lose any FPU registers during a context switch. Similar to the non-scratch registers, we have to manually save and restore them.

machine/context.h

struct Context {
    intptr_t ebx;
    // ...
    uint8_t fpu[108]; 
} __attribute__((packed));

‍Reserving enough space for the x87 FPU.

machine/context.asm

align 8
context_switch:
    ; Save non-scratch register of current context
    ; ...
    fsave [eax + 20]
    ; Save stack pointer of current context
    mov [eax + 16], esp
 
    ; Load stack pointer of next context
    mov esp, [ecx + 16]
    ; Restore register of next context
    ; ...
    frstor [ecx + 20]
 
    ret
 
align 8
context_launch:
    ; Load stack pointer of next context
    mov esp, [eax + 16]
    ; Restore register of next context
    ; ...
    finit
 
    ret

‍Just add the finit, fsave, and frstor instructions to the context switch functions.

3. Bonus: Sqrt, Sin, Cos

Now, all basic floating-point operations should work as expected (add, sub, mul, div). However, some special operations, like sqare-root or sine, are still missing. Luckily, the x87 FPU has special instructions for most of them. Still, we have to implement them manually.

utils/math.h

const float PI = 3.1415926;
inline float sqrt(float x) {
    asm("fld %0; fsqrt; fstp %0" : "=m"(x) : "m"(x));
    return x;
}
inline float sin(float x) {
    asm("fld %0; fsin; fstp %0" : "=m"(x) : "m"(x));
    return x;
}
inline float cos(float x) {
    asm("fld %0; fcos; fstp %0" : "=m"(x) : "m"(x));
    return x;
}

‍Of course, you can also use the Fast Inverse Sqrt from Quake III.