exp, log, sqrt and sin from adds and multiplies

Why not just call expf?

C has expf, logf, sqrtf and sinf. But each C library computes them its own way, and two machines can disagree in the last bit. The grader compares f32 results bit for bit, so the course needs one definition everyone can match. The spec gives one: each math op is a specific program of core ops. If your compiler runs exactly those steps, your bits match the oracle's on any machine.

There's a second reason, and it's the one the syllabus cares about. A GPU, or a small accelerator, often has no exp instruction at all. Knowing how to build one from multiplies is how compilers for such hardware work.

How an f32 is stored

Every f32 is 32 bits: a sign bit, an 8-bit exponent field E, and a 23-bit mantissa M. For normal numbers the value is (-1)^sign · 1.M · 2^(E - 127). The 127 is a bias, so exponents from -126 to 127 fit in a field from 1 to 254. When the exponent field is 0, the number is subnormal: its value is (-1)^sign · 0.M · 2^-126. These are the tiny numbers below 2^-126, with fewer bits of precision.

The new op bitcast x dtype=i32 gives you those 32 bits as an integer, and bitcast n dtype=f32 goes the other way. Nothing is converted: the same bits are read as the other type. That's what lets the math ops pull a number apart and put one together.

exp2: split, then a polynomial

2^x for any x comes from two easy cases. Round x to the nearest integer n and let f = x - n, so f is in [-0.5, 0.5]. Then 2^x = 2^f · 2^n.

• 2^n needs no arithmetic: it's the f32 with exponent field n + 127 and mantissa 0. Build the integer (n + 127) · 2^23 and bitcast it.
• 2^f on a small interval is a polynomial. 2^f = e^(f ln 2) = Σ (f ln 2)^k / k!. With |f| ≤ 0.5, eight terms (k = 0 to 7) are more precise than f32 can store.

The spec evaluates the polynomial with Horner's rule: c0 + f(c1 + f(c2 + ...)) with c_k = (ln 2)^k / k! (1.0, 0.6931472, 0.2402265, ...), starting from the highest coefficient. That's one multiply and one add per term, and it's also more accurate than adding up powers.

Two details make the edges right. First, NaN is caught before anything else, because casting NaN to i32 is an error. Second, x is clamped to [-150, 128]. 2^128 is already infinity in f32, and 2^-150 already rounds to 0, so clamping changes no answer and keeps n small. The bit trick (n + 127) · 2^23 only works for n from -126 to 127, and after the clamp n runs from -150 to 128. Below -126 the result is subnormal, and at 128 a single 2^128 would already be infinity even when the final result (say exp2 127.9) is finite. So the spec multiplies by 2^h and then 2^k, with h + k = n and both in [-75, 64], safely normal.

log2: read the exponent, fix up the rest

log2 runs the same idea backward. Bitcast x to an integer b. Its exponent field is b div 2^23, so e = b div 2^23 - 127. Its mantissa bits are b mod 2^23. Adding the bits of 1.0 (1065353216) gives the f32 with that mantissa and exponent 0, which is m, in [1, 2), with x = m · 2^e. So log2 x = e + log2 m.

A series for log2 m converges fastest when m is close to 1, so if m > 1.414 the spec halves it and adds 1 to e. Now m is in [0.707, 1.414]. With s = (m - 1) / (m + 1), |s| ≤ 0.172 and

log2 m = (2 / ln 2) · (s + s³/3 + s⁵/5 + s⁷/7 + s⁹/9 + ...)

Five terms give full f32 precision. A subnormal x has exponent field 0, so it's multiplied by 2^23 first and 23 is subtracted at the end. The special cases are checked at the end: zero gives -inf, a negative gives NaN, and inf and NaN pass through.

sqrt: a guess, then Newton

sqrt x = 2^(log2(x) / 2) is a fine first guess, but two rounded steps leave it a few units in the last place (ulps) off. One ulp is the gap between an f32 and the next one. One step of Newton's method fixes that. If y is close to √x, then (y + x/y) / 2 is much closer, because the error roughly squares each step. Perfect squares come out exact: sqrt 16.0 is exactly 4.0.

sin: subtract turns, fold, then Taylor

sin repeats every 2π, so first subtract the nearest whole number of turns: k = floor(x / 2π + 0.5), then r = x - k · 2π. 2π isn't an f32, so the spec subtracts it in two parts: 6.28125, which is exact, and a small correction. That keeps r accurate for x up to about 10,000.

Now r is in about [-π, π]. Since sin(π - r) = sin r, folding brings it to [-π/2, π/2], where the Taylor series r - r³/3! + r⁵/5! - ... up to r^13 is accurate. For a huge x, f32 can't even say which turn x is on, so the result is meaningless. The spec still clamps it so it stays in [-1, 1].

Special values fall out

Nobody special-cases exp2(-inf) = 0: the clamp sends it to -150, and 2^-150 rounds to 0. Most of the table in the spec works like that. The cases that need an explicit where are the ones where a step would otherwise fail or lie: NaN before a cast, zero and negatives before a log, and infinity, whose bits look like a huge exponent.

One more rule makes NaN predictable. There are millions of NaN bit patterns, and different hardware makes different ones (x86 produces a negative NaN for 0 · inf). bitcast of any NaN gives 0x7fc00000, so bitcast never exposes which NaN you had. In C, check isnan first.