Topic — Precision: how the PE rounds, and references that match

Needs: Lesson 2 — Values, control flow and functions, Lesson 3 — How a PE runs a loop: vector instructions and cycles.

Lesson 3 showed that the compiler fuses a multiply and an add into one instruction, and that a reference must round the same way. This topic shows where the PE rounds in each floating-point form a kernel commonly uses: conversions between hp, ep and sp, fused and split multiply-adds, dot products, and the math fast and math ehalf pragmas. It then builds a host reference that matches every one of them bit for bit.

What you will be able to do

• Say how many significant bits hp, ep and sp keep, and what a conversion between them does.
• Predict where a multiply-add, a dot product and an accumulation round, from the instructions in the listing.
• Use math fast and math ehalf, knowing exactly which bits they drop and what they gain.
• Write a host reference that rounds where the hardware rounds, and so matches it exactly.

The program

One PE. precision.w receives four input arrays from the image (miset): ah and bh, hp values in [0.5, 1), and as and bs, sp values in [0.5, 2). It computes five groups of results from them:

          ┌──────────────────────── one PE ───────────────────────────┐
 ah, bh ─>│ formats      xh = as as hp        xe = as as ep           │
 as, bs ─>│ multiply-add fh = bh + kh·ah      gh = the same, split    │
          │              fs = bs + ks·as      gs = the same, split    │
          │ dot products dh[0..1] (hp)        ds[0..3] (sp)           │
          │ pragmas      mf, mfa (math fast)  ea, ef (math ehalf)     │
          │ ep           em = xe·xe                                   │
          └───────────────────────────────────────────────────────────┘

kh and ks are the constant 0.3 in hp and in sp. check.py compares every result with an exact model in ref.py, and the test passes only if all of them match.

The three formats

Type	Sign, exponent, mantissa bits	Significant bits	Largest value	Smallest normal	Words
hp	1, 5, 10	11	65504	2⁻¹⁴	1
ep	1, 8, 10	11	about 3.4 × 10³⁸	2⁻¹²⁶	2
sp	1, 8, 23	24	about 3.4 × 10³⁸	2⁻¹²⁶	2

hp is IEEE half precision. ep has the precision of hp and the range of sp. It is stored exactly like an sp, with the low 13 mantissa bits zero, so it saves no memory. It exists because the multiplier does four ep products per cycle and one sp product per cycle.

∀i xh[i] <- as[i];
∀i xe[i] <- as[i];

Both conversions round to nearest, with ties to even. The listing shows that the ep conversion is an add of zero whose result is rounded to ep:

	fs2h [T24] = [T19];  //@../precision.w:41
	fadde [T25] = [T03], D6;  //@../precision.w:42

The zero it adds is +0.0, so the conversion turns -0.0 into +0.0. The conversion to hp keeps the sign of -0.0.

A multiply-add: one statement or two

∀i fh[i] <- bh[i] + kh * ah[i];
t[1] <- tsc[0];
∀i ph[i] <- kh * ah[i];
∀i gh[i] <- bh[i] + ph[i];

t[1] <- tsc[0]; reads the timestamp counter, as in lesson 3; the kernel reads it between the statements it times. The first statement is one fused multiply-add. It adds the exact product to bh[i] and rounds once. The two statements round twice: the product when it is stored in ph, and then the sum:

	fmach [T26] = [T04], [T1a], R3;  //@../precision.w:47
	fmulh [T27] = R3, [T1b];  //@../precision.w:49
	faddh [T28] = [T05], [T1c];  //@../precision.w:50

On the tutorial row, 6 of the 62 hp results and 10 of the 62 sp results differ between the two forms. Neither form is wrong; the reference must follow the form the kernel uses.

Dot products

dh[0] <- ∑ k in [0, N) ah[k] * bh[k];
t[5] <- tsc[0];
dh[1] <- 0.0:hp;
∀i dh[1] <+- ah[i] * bh[i];

∑ k in [0, N) e is the sum of e for k from 0 to N − 1; it declares k itself. a <+- e adds e to a, so ∀i dh[1] <+- ah[i] * bh[i]; adds each product to dh[1].

The PE has a dot-product unit with four accumulators. An accumulator keeps 24 significant bits, as sp does, with a 10-bit exponent, so a sum in it cannot overflow or underflow. ∑ compiles to it:

	fh2s [T2c] = [T10];  //@../precision.w:59
	fh2s [T2d] = [T11];  //@../precision.w:59
	fdps DPR0.clear = [T08], [T12];  //@../precision.w:59
	fstdpas De = DPR0;  //@../precision.w:59
	fs2h Rd = De;  //@../precision.w:59

For hp operands, Carbide first copies both arrays to sp (fh2s), then runs the sp dot product, fdps, which takes one element at a time. fdps rounds each product to the accumulator's 24 bits, then adds it and rounds again. A product of two hp values is exact in 24 bits, so here only the additions round. fstdpas reads the accumulator as an sp, and fs2h rounds that once to hp.

<+- into an array element is one fused multiply-add per element, fmach, each rounded to the element's type. It runs element by element: about 9 cycles per element, against under 2 for ∑ (see the cost table below). On the tutorial row:

Result	Value	Rounded
exact sum of the 62 products	35.44886326789856
dh[0], ∑ into hp	35.4375	to 24 bits at each step, then once to hp
dh[1], <+- into an hp element	35.46875	to hp at each of the 62 steps

dh[0] is the hp value nearest the exact sum. The chain ends one hp step (2⁻⁵, near 35) above it.

The same pair in sp, ds[0] and ds[1], can differ: the dot product rounds each product before it adds it, and the fused multiply-add does not. On the tutorial row they happen to agree in every bit, 91.30236. On the bank row (N = 64) they are 98.222702 and 98.222717, two units in the last place apart.

Into a scalar variable, <+- can use the dot-product unit as well. It does so for sp operands into an sp scalar and for hp operands into an hp scalar. With hp operands into an sp scalar it became a chain of fused multiply-adds, one per element. Write ∑ when you want the dot-product unit.

ds[2] <- ∑ k in [0, N) ah[k] * bh[k];
ds[3] <- ∑ k in [0, N) (sp)ah[k] * (sp)bh[k];

A ∑ of hp operands is an hp value, even when its destination is an sp. ds[2] is 35.4375, the hp result widened. Casting the operands makes the sum an sp value: ds[3] is 35.448868.

math fast and math ehalf

[[math.fast(on)]] ∀i mf[i] <- as[i] * bs[i];
t[10] <- tsc[0];
∀i ps[i] <- as[i] * bs[i];
t[11] <- tsc[0];
[[math.fast(on)]] ∀i mfa[i] <- bs[i] + ks * as[i];
[[math.ehalf(on)]] ∀i ea[i] <- as[i] + bs[i];
[[math.ehalf(on)]] ∀i ef[i] <- bs[i] + ks * as[i];
∀i em[i] <- xe[i] * xe[i];

A pragma in [[ ]] applies to the statement after it. ps is the plain product, for comparison. math fast and math ehalf change sp arithmetic to the ep instructions:

Statement	Instruction	The factors	The result
a * b, plain	fmuls	used in full	rounded to sp
a * b, math fast	fmule	cut to ep	rounded to ep
c + k * a, math fast	fmaces	cut to ep	exact product added, rounded to sp
a + b, math ehalf	fadde	(an add: used in full)	rounded to ep
c + k * a, math ehalf	fmace	cut to ep	exact product added, rounded to ep

"Cut to ep" means that the low 13 mantissa bits of each factor are dropped. That rounds the factor toward zero, not to nearest: with math fast, 1.00146484375 × 1.0 gives 1.0009765625, where rounding the factor to nearest even would give 1.001953125. The result then rounds to nearest. So 3.14159265 × 2.0 gives 6.28125: about three significant digits survive.

An ep variable's value already has its low 13 bits zero, so em, a product of two ep values, compiles to the same fmule without a pragma and loses nothing in its factors.

What math fast buys is throughput. On the tutorial row (N = 62) and the bank row (N = 64):

Statement	N = 62: cycles	per element	N = 64: cycles	per element
hp b + k*a, one statement (fmach)	20	0.32	22	0.34
hp b + k*a, two statements	31	0.50	32	0.50
sp b + k*a, one statement (fmacs)	64	1.03	130	2.03
sp b + k*a, two statements	95	1.53	98	1.53
sp a*b (fmuls)	64	1.03	130	2.03
sp a*b, math fast (fmule)	38	0.61	39	0.61
hp ∑	112	1.81	177	2.77
hp <+- into an element	564	9.10	582	9.09
sp ∑	74	1.19	140	2.19
sp <+- into an element	562	9.06	580	9.06

On the tutorial row, math fast makes the sp multiply 1.7 times faster, and one fused statement is faster than two. fmule does four multiplies per cycle, but it does not reach four elements per cycle here: memory delivers at most four 32-bit words per cycle, which is two elements of each of two sp arrays (lesson 3).

The bank row shows a memory-bank conflict, the one lesson 3 describes. At N = 64, as and bs are 128 words apart, a multiple of 8, so as[i] and bs[i] always lie in the same memory bank. fmacs, fmuls and the ∑s, whose operands Carbide also places 128 words apart, take one element of each array at a time, so each waits a cycle per element. That doubles their cost, and makes the split sp statement look faster than the fused one. fmule and the add of the split statement take the same time on both rows, although the add's operands bs and ps are 640 words apart, also a multiple of 8. Like the sp addition in lesson 3, they read four elements of one array at a time, one word from every memory bank, so their two operands take turns at the memory banks wherever the arrays lie. The address rule tells you where a conflict can happen, not which instruction waits. Two sp arrays declared one after the other are a multiple of 8 words apart whenever N is a multiple of 4.

cycles.py prints this table on every run. The numbers depend on the compiler build and on where the arrays lie, so the test never checks them.

A reference that matches

ref.py computes with exact rational numbers (Python's fractions.Fraction) and rounds only where the hardware rounds. One function does all the rounding:

HP = (11, -14)     # 10 stored mantissa bits; subnormal below 2**-14
EP = (11, -126)    # hp's significand, sp's exponent range
SP = (24, -126)
ACC = (24, -510)   # dot-product accumulator: sp's significand, a 10-bit exponent

rne(x, fmt) rounds the exact value x to nearest, ties to even, in the format fmt. Each model is then one line that says where the hardware rounds:

def fma(c, a, b, fmt):
    """c + a*b with the exact product, rounded once."""
    return rne(exact(c) + exact(a) * exact(b), fmt)


def mul_add(c, a, b, fmt):
    """c + a*b in two statements: the product is rounded, then the sum."""
    return rne(exact(c) + rne(exact(a) * exact(b), fmt), fmt)

check.py prints, for each output, how many elements match its model, and how many a reference that rounds differently would match. On the tutorial row:

array  model                                          match   a reference that uses            would match
xh     as rounded to hp                             62/62     as with the low bits dropped        31/62
xe     as rounded to ep                             62/62     as with the low bits dropped        31/62
fh     bh + kh*ah, rounded once                     62/62     product rounded first               56/62
gh     product rounded, then the sum                62/62     rounded once                        56/62
fs     bs + ks*as, rounded once                     62/62     product rounded first               52/62
gs     product rounded, then the sum                62/62     rounded once                        52/62
dh     dot unit, then hp; hp FMA chain               2/2      exact sum rounded to hp              1/2
ds     dot unit; sp FMA chain; hp sum; sp sum        4/4      all four rounded once from exact     2/4
mf     factors cut to ep, product rounded to ep     62/62     factors rounded to ep               23/62
ps     as*bs rounded to sp                          62/62     math fast's result                   0/62
mfa    factors cut to ep, sum rounded to sp         62/62     factors rounded to ep                0/62
ea     as + bs rounded to ep                        62/62     terms cut to ep first               30/62
ef     factors cut to ep, sum rounded to ep         62/62     sum rounded to sp                    0/62
em     xe*xe rounded to ep                          62/62     rounded to sp                        2/62

Three rules give a reference that matches:

1. Round where the hardware rounds. Read the listing: fmach rounds once, fmulh and then faddh round twice, fdps rounds each product and each sum.
2. Compute exactly between the roundings. NumPy's float16 and float32 arithmetic round after every operation, so it cannot model a fused multiply-add. Exact rationals or integers can.
3. Model each input the way the instruction reads it. An ep multiply drops the low bits of its factors; a conversion rounds them.

When an implementation is allowed to differ from the correctly rounded result, as the library's division is, compare within a tolerance in units in the last place instead (see the division topic).

What the hardware does

• hp, ep and sp instructions round to nearest even by default. The rounding mode is set in the floating-point control word, and the compiler never sets that word.
• A fused multiply-add adds the infinitely precise product, then rounds once.
• An instruction that takes ep factors drops their low 13 mantissa bits.
• The four dot-product accumulators keep 24 significant bits and a 10-bit exponent. The sp dot product, fdps, takes one element at a time. Storing an accumulator as sp rounds only when its exponent is out of sp's range.
• One instruction can do eight hp multiplies, four ep multiplies, or one sp multiply per cycle.

See the architecture notes for the PE as a whole.

Run it

In your copy of the tutorial (see Lesson 1):

cd tutorial/ext-numerics/precision
make TEST=tutorial     # N=62
make TEST='*'          # also N=64, N=16 and N=201

A passing row prints a line ending in PASS tutorial [ N=62 SEED=1 ARCH=sdr ], and test-tutorial/make-log holds the check table and the cycle table, the first ending in Ok!. If an output differs from its model, its match column shows fewer than all elements and the log ends with FAIL: an output differs from its model.

After you edit a source file, run make clean before the next test.

Pitfalls

• math fast costs digits silently. An sp product under math fast keeps about three significant digits. Use it where the data has no more than that. On run-time xp and u16 operands math fast selects imul11, which multiplies only the low 11 bits of each operand: variables holding 2048 and 2 multiply to 0, and -3 and 5 to 10225. A product of two literals is computed by the compiler and is not affected. Keep integer multiplies out of the statements that math fast covers.
• A misspelled pragma does nothing. [[math.fats(on)]] compiles without a warning and the multiply stays fmuls. Check the listing whenever a pragma matters.
• A pragma on a declaration with an initializer does nothing. [[math.fast(on)]] sp x = a * b; gives fmuls. Declare x first and put the pragma on the assignment.
• ∑ of hp values is hp. Storing it in an sp does not restore the lost bits. Cast the operands to sp inside the ∑.
• <+- into an array element rounds every step. In hp the error grows with the length of the sum, and the loop runs element by element. Use ∑ for a long sum.
• A cost measured at one size can be a memory-bank conflict. Before you compare two forms of an sp loop, time them at a size where the arrays do not collide, such as N = 62 here.
• A NumPy reference in float16 models two statements, not one. kh * ah + bh computed on NumPy float16 arrays disagreed with the fused fh in 6 of 62 elements, and agreed with gh.

Exercises

1. Change kh to 0.25:hp in precision.w and to np.float16(0.25) in check.py, run make clean and the test, and look at the fh and gh rows of the table. Each now matches the other row's model too: a product with a power of two is exact, so the extra rounding changes nothing.
2. Put [[math.fast(on)]] on the fs statement. The test fails. Change fs's model in check.py to the one mfa uses, and the test passes again.

Next

Topic — Bits and wide integers works with integer bit patterns: shifts, rotations and 64-bit arithmetic in 16-bit pieces.

References

• Extended precision (ep)
• Dot product
• Floating-point evaluation
• Floating-point precision pragmas
• Casts