A picture is worth a thousand LLM prompts:

Horizontal value is the input value x for cosf(x), ranging from 0 to 2π. The orange line shows the result of evaluating cosf(x), which looks exactly as expected. The blue line shows the time it took to evaluate cosf(x) for each x, in nanoseconds.

I discovered this while investigating why my main loop would occasionally stall for whole milliseconds (that’s hundreds of thousands of cycles for my RP2350!), causing USB HID reports to be dropped.

After a while, I figured out that the culprit was this code (slightly simplified):

level += level_step;
const float rampup = 0.5f - 0.5f * cosf(level);
phase += phase_step;
tx[i] = sinf(phase) * volume * rampup;

This generates a pure sine tone with a raised-cosine filter applied over it to avoid clicks from harmonics. Specifically, this happened when level would cross over π/2, and would go away entirely if I stubbed the raised-cosine filter with identity instead of cosf(), for instance.

With some help from Claude, I located the issue in the ARM assembly code of the implementation of cosf() in the Pico SDK. First, below are the relevant bits of the main block of cosf():

wrapper_func cosf
[…]
 bls 10b                     @ argument ≥1? needs reduction; also Inf/NaN […]
5:
@ here we have a quadrant count in r12 and a signed offset r0 from r12*π/2
[…]
 bx r14

The interesting part is that it handles quadrants.

In short, trigonometric functions are periodic, and it is most convenient to evaluate them with a small argument. So, instead of evaluating cosf(2π + 0.1f), it will actually evaluate (simplifying again) cosf(0.1f). Further, as you can see below, the range [π; 2π] is just the mirror image of the range [0;π] (line symmetry around x = π). And [π/2;π] is itself a mirror image of [0;π/2] (point symmetry around (π/2, 0)). That means cosf() does a small transformation depending on whether you are in [0;π/2], [π/2;π], [π; 3π/2] or [3π/2; 2π]; these are the quadrants.

Thus, most implementations will reduce the argument x to the first quadrant, evaluate cosf(x), and transform the result depending on what quadrant x was in.

In the code above, bls 10b jumps to the code labelled 10 backwards, if the argument is not in the first quadrant. Or, more exactly, if it is larger than 1.0f, due to optimization details. The corresponding code is shown (again, simplified) below.

The “range reduction” is finding the quadrant. The important thing to notice is that it finishes by jumping to the code label 5 forwards, which actually just brings us back to the main block.

@ light-duty range reduction
[…]
10:
[…]
 blo 40b                     @ very small result? use heavy-duty
[…]
 b 5f

In summary:

float cosf(float x) {
    if (x > 1.0f) {
        // find quadrant
    }
    // actually evaluate cosf(x) for first quadrant
    // adjust for quadrant
    return res;
}

From this, you clearly understand the step in the timing graph (reproduced below).

Now, the other thing to notice in that “light-duty range reduction” is the “very small result? use heavy-duty” comment, associated with a jump to code label 40 backwards. That code is responsible for the huge spikes.

There is no point in looking at it, but we should consider why this special case exists. As hinted by the comment, cos(π/2) = 0, so any argument that is very close to π/2 (such as M_PI_2) will evaluate to something that is very close to 0.

Small errors look large compared to small values. And one such small error is in the representation of π itself: C’s M_PI cannot be exactly equal to π. If the code just did cosf(x - M_PI_2), it would be introducing a relatively significant error into the calculation. So, the authors of the Pico SDK decided to automatically switch to another approach that evaluates cos(x – π/2) without just subtratcing M_PI_2 from x. Even if that means taking surprisingly longer to run the computation.

Now, then, it makes sense at π/2 and 3π/2 because cos(π/2) and cos(3π/2) are both zero. But what about the two other spikes? cos(π) = -1 and cos(2π) = 1, not 0, so there is no relative error explosion.

Simple: because that same code is also used for sinf(), which does have zeros at π and 2π. In other words: cosf() is unnecessarily slow near π and 2π and sinf() is unnecessarily slow near π/2 and 3π/2.

I would also argue that such a behavior is a bad trade-off on a microcontroller. There, you typically care much more about predictable timings than exact arithmetic. In fact, glibc’s cos() and sin() just don’t do that and are just as precise.

In any case, that’s something to be aware of when writing software for the RP2040 or the RP2350.