Hi, I've been thinking a bit more on the structure of polynomial evaluation, which at a high level is rather similar for ghash and gcm.
** Intro **
The function to be computed is
R_j = K (R_{j-1} + M_j)
where K is the secret key and M_j are the message blocks. Operations take place in some finite field. With n message blocks, M_0, ..., M_{n-1}, and initial R_{-1} = 0, we get
R_{n-1} = m_{n-1} K + m_{n-2} K^2 + ... m_0 K^n
I.e., a degree n polynomial with coefficients M_j, constant term = 0, evaluated at the point K.
To be concrete, consider a 64-bit architecture, and a finite field where elements are represented as two words (for poly1305, we need two words plus a few extra bits, but ignore that now, for simplicity. And for ghash, let's also ignore the complications from bit-reversal). Let B represent the bignum base. I'm going to make this a bit handwavy, but we'll have B = 2^64 or B = x^64 depending on type of field.
The finite field is represented by a mod P operation, where the structure of P is nice, with a leading one bit followed by more than 64 zeros, and then a few more non-zero bits at the end. This implies that we can define a multiply operation
Y_2 B^2 + Y_1 B + Y_0 = (X_1 B + X_0) (K_1 B + K_0) (mod P)
by four independent independent multiplication instructions (widening, 64x64 --> 128) involving X and some precomputed values depending on K, and accumulation only involving shifts/adds that have low latency. The result is not fully reduced mod P; it consists of three words.
We can then reduce this to two words by multiplying Y_2 by a suitable single-word constant, with final two-word result
R = C Y_2 + Y_1 B + Y_0
Applying this to the original recurrency, R_j = K (R_{j-1} + M_j), the X input corresponds to R_{j-1} + M_j, so the critical dependency path from R_{j-1} to R_j includes *two* multiply latencies. E.g, if multiply latency is 5 cycles, it's not possible to get this evaluation scheme to run faster than 10 cycles per block. (In practice, the accumulation will also contribute one or a few cycles to the critical path, but *much* less than those two multiplies).
So question is, how can we do better?
** Postponing reduction **
The approach taken in the new x86_64 poly1305 code I just pushed is to skip the final reduction, and let the state be one word larger (and this is particularly cheap for poly1305, because we don't quite increase state size by a word, but from two words + 3 bits to two words plus ~60 bits). The multiply operation becomes
Y_2 B^2 + Y_1 B + Y_0 = (X_2 B^2 + X_1 B + X_0) (K_1 B + K_0) (mod P)
This can be arranged with 6 independent multiply instructions + cheap accumulation. (I haven't worked out the details for the ghash case, but I do expect that it's rather practival there too).
Then the dependency chain from one block to the next is reduced to one multiply latency, 5 cycles in our example. In case all other needed instructions can be scheduled (manually, or by the processor's out-of-order machinery) to run in 5 cycles in parallel with the multiplies, we would get a running time of 5 cycles per block.
** Interleaving **
The other approach, used in the recent powerpc gcm code, is to interleave multiple blocks. For simplicity, only consider 2-way interleaving here. The key thing is that if we expand te recurrency once, we get
R_j = K (M_j + K (R_{j-2} + M_{j-1})) = K M_j + K^2 (R_{j-2} + M_{j-1})
We get two field multiplications, but one of them, K M_j, is completely independent of previous blocks (R_{j-2}), and can be computed in parallel. It may add a cycle or so to accumulation latency, but we can do essentially twice as much work without making the critical path longer. We get 8 independent multiply instructions, and one dependent for the final folding. Can be extended to more than two blocks if needed (depending on number of available registers).
Another variant could be to separate even and odd parts of the polynomial being evaluated, and evaluate both parts at K^2. We can then compute the two recurrencies
E_j = K^2 (E_{j-1} + M_{2j}) O_j = K^2 (O_{j-1} + M_{2j+1})
in parallel. It's unclear to me what the pros and cons are compared to previous variant. One may get some advantage from both multiplies using the same factor K^2. On the other hand, each recurrency has to be accumulated and folded separately, which costs instructions and registers. Maybe more useful for hardware implementation? This variant is currently not used in Nettle.
** Doing both **
It's possible to combine those two tricks. Processing of two blocks would then be an operation of the form
Z_2 B^2 + Z_1 B + Z_0 = (X_2 B^2 + X_1 B + X_0) K^2 + (Y_1 B + Y_0) K
Here, the Xs (three words) represent R_{j-2} + M_{j-1}, the Ys repreent M_{j-2}, and the Zs represents R_j, as three words (without final folding). We would need 10 independent multiples, one more than with plain interleaving, but critical path includes only one multiply latency.
I think this is a promising alternative, if one would otherwise need to interleaving a large number of blocks to get full utilization of the multipliers.
** How to choose **
When implementing one of those schemes, different processor resources may be the bottleneck. I'd expect it to be one of
o Multiply latency, i.e, latency of the dependency chain from one block to the next (including also a few additions, but multiply latency willb e the main part). If this is the bottleneck, it means all other instructions can be scheduled in parallel, and the processor will sit idle for some cycles, waiting for a multiply to complete. Typical latency for multiply is 5 times longer than for an addition (but ratio difers quite a bit between processors, of course)
o Multiply throughput, i.e., the maximum number of (independent) multiply instructions that can be run per cycle. Typical number is 0.5 -- 2. If this is the bottleneck, the processor will spend some cycles idle, waiting for a multiplier to be ready to accept a new input.
o A superscalar processor can issue several instructinos in the same cycle, but there's a fix small limit. Typical number is 2 -- 6. So, e.g., if the processor can issue maximum 4 instructions per cycle, the evaluation loop consists of 40 instructions, and the loop actually runs in close to 10 cycles per iteration, then instruction issue is the bottleneck.
The tricks discussed in this note are useful for finding an evaluation scheme where multiply latency isn't a bottleneck. But once a loop hits the limit on multiply throughput or instructions per cycle, other tricks are needed to optimize further. In particular, the postponed reduction has a cost in multiply throughput, since it needs some additional multiply instruction.
I think one should aim to hit the limit on multiply throughput; that one is hard to negotiate (it's possible to reduce the number of multiply instructions somewhat, by the Karatsuba trick, but due to the additional overhead, likely to be useful only on processors with particularly low multiply throughput).
Regards, /Niels
nisse@lysator.liu.se (Niels Möller) writes:
Y_2 B^2 + Y_1 B + Y_0 = (X_2 B^2 + X_1 B + X_0) (K_1 B + K_0) (mod P)
This can be arranged with 6 independent multiply instructions + cheap accumulation. (I haven't worked out the details for the ghash case, but I do expect that it's rather practival there too).
I've found a rather straight forward way to express that.
Recall that for ghash, due to the bit-reversal, the multiply operation of interest is
M H x^{-128} mod P
where structure of P means that x^{-64} = x^{64] + P_1, and P_1 is a single word. Split M and H into halves,
M = M_1 x^{64} + M_0 H = H_1 x^{64} + H_0
The previous notes defines the precomputation of
D_1 x^{64} + D_0 = H_0 x^{64} + H_1 + H_0 P_1
Alternatively, D can be defined as D = x^{-64} H. And the accumulation part can then be written as
(M_1 x^64 + M_0) H x^{-128} = (M_1 H + M_0 D) x^{-64}
As before, accumulate this in two 128-bit registers R and F, as
(M_1 x^64 + M_0) H x^{-128} = R + F x^{-64}
with
R = M_1 H_1 + M_0 D_1 F = M_1 H_0 + M_0 D_0
If we add one more unreduced word to M,
M = M_1 x^{64} + M_0 + M_{-1} x^{-64}
all we need is to precompute one more constant E = H x^{-128} = D x^{-64}, in the same way
E_1 x^{64} + E_0 = D x^{-64} = D_0 x^{64} + D_1 + D0 P_1
and we get one more term each for R and F,
R = M_1 H_1 + M_0 D_1 + M_{-1} E_1 F = M_1 H_0 + M_0 D_0 + M_{-1} E_2
At the end of the iteration, just add the high half of F into R, but keep F_0 as an input (the place of the M_{-1}) for the next iteration.
Regards, /Niels
On Sat, Jan 29, 2022 at 4:29 PM Niels Möller nisse@lysator.liu.se wrote:
Hi, I've been thinking a bit more on the structure of polynomial evaluation, which at a high level is rather similar for ghash and gcm.
** Intro **
The function to be computed is
R_j = K (R_{j-1} + M_j)
where K is the secret key and M_j are the message blocks. Operations take place in some finite field. With n message blocks, M_0, ..., M_{n-1}, and initial R_{-1} = 0, we get
R_{n-1} = m_{n-1} K + m_{n-2} K^2 + ... m_0 K^n
I.e., a degree n polynomial with coefficients M_j, constant term = 0, evaluated at the point K.
To be concrete, consider a 64-bit architecture, and a finite field where elements are represented as two words (for poly1305, we need two words plus a few extra bits, but ignore that now, for simplicity. And for ghash, let's also ignore the complications from bit-reversal). Let B represent the bignum base. I'm going to make this a bit handwavy, but we'll have B = 2^64 or B = x^64 depending on type of field.
The finite field is represented by a mod P operation, where the structure of P is nice, with a leading one bit followed by more than 64 zeros, and then a few more non-zero bits at the end. This implies that we can define a multiply operation
Y_2 B^2 + Y_1 B + Y_0 = (X_1 B + X_0) (K_1 B + K_0) (mod P)
by four independent independent multiplication instructions (widening, 64x64 --> 128) involving X and some precomputed values depending on K, and accumulation only involving shifts/adds that have low latency. The result is not fully reduced mod P; it consists of three words.
We can then reduce this to two words by multiplying Y_2 by a suitable single-word constant, with final two-word result
R = C Y_2 + Y_1 B + Y_0
Applying this to the original recurrency, R_j = K (R_{j-1} + M_j), the X input corresponds to R_{j-1} + M_j, so the critical dependency path from R_{j-1} to R_j includes *two* multiply latencies. E.g, if multiply latency is 5 cycles, it's not possible to get this evaluation scheme to run faster than 10 cycles per block. (In practice, the accumulation will also contribute one or a few cycles to the critical path, but *much* less than those two multiplies).
So question is, how can we do better?
** Postponing reduction **
The approach taken in the new x86_64 poly1305 code I just pushed is to skip the final reduction, and let the state be one word larger (and this is particularly cheap for poly1305, because we don't quite increase state size by a word, but from two words + 3 bits to two words plus ~60 bits). The multiply operation becomes
Y_2 B^2 + Y_1 B + Y_0 = (X_2 B^2 + X_1 B + X_0) (K_1 B + K_0) (mod P)
This can be arranged with 6 independent multiply instructions + cheap accumulation. (I haven't worked out the details for the ghash case, but I do expect that it's rather practival there too).
Then the dependency chain from one block to the next is reduced to one multiply latency, 5 cycles in our example. In case all other needed instructions can be scheduled (manually, or by the processor's out-of-order machinery) to run in 5 cycles in parallel with the multiplies, we would get a running time of 5 cycles per block.
** Interleaving **
The other approach, used in the recent powerpc gcm code, is to interleave multiple blocks. For simplicity, only consider 2-way interleaving here. The key thing is that if we expand te recurrency once, we get
R_j = K (M_j + K (R_{j-2} + M_{j-1})) = K M_j + K^2 (R_{j-2} + M_{j-1})
We get two field multiplications, but one of them, K M_j, is completely independent of previous blocks (R_{j-2}), and can be computed in parallel. It may add a cycle or so to accumulation latency, but we can do essentially twice as much work without making the critical path longer. We get 8 independent multiply instructions, and one dependent for the final folding.
But in radix 2^64, doesn't K^2 need to be represented in 3 words as described in a previous note so considering that we get 10 independent multiply instructions rather than 8?
Can be extended to more than two blocks if needed (depending on number of available registers).
Another variant could be to separate even and odd parts of the polynomial being evaluated, and evaluate both parts at K^2. We can then compute the two recurrencies
E_j = K^2 (E_{j-1} + M_{2j}) O_j = K^2 (O_{j-1} + M_{2j+1})
in parallel. It's unclear to me what the pros and cons are compared to previous variant. One may get some advantage from both multiplies using the same factor K^2. On the other hand, each recurrency has to be accumulated and folded separately, which costs instructions and registers. Maybe more useful for hardware implementation? This variant is currently not used in Nettle.
** Doing both **
It's possible to combine those two tricks. Processing of two blocks would then be an operation of the form
Z_2 B^2 + Z_1 B + Z_0 = (X_2 B^2 + X_1 B + X_0) K^2 + (Y_1 B + Y_0) K
Here, the Xs (three words) represent R_{j-2} + M_{j-1}, the Ys repreent M_{j-2}, and the Zs represents R_j, as three words (without final folding). We would need 10 independent multiples, one more than with plain interleaving, but critical path includes only one multiply latency.
Same matter here, if we consider 3 words for K^3 this sums up the number of independent multiples to 13.
I think this is a promising alternative, if one would otherwise need to interleaving a large number of blocks to get full utilization of the multipliers.
** How to choose **
When implementing one of those schemes, different processor resources may be the bottleneck. I'd expect it to be one of
o Multiply latency, i.e, latency of the dependency chain from one block to the next (including also a few additions, but multiply latency willb e the main part). If this is the bottleneck, it means all other instructions can be scheduled in parallel, and the processor will sit idle for some cycles, waiting for a multiply to complete. Typical latency for multiply is 5 times longer than for an addition (but ratio difers quite a bit between processors, of course)
o Multiply throughput, i.e., the maximum number of (independent) multiply instructions that can be run per cycle. Typical number is 0.5 -- 2. If this is the bottleneck, the processor will spend some cycles idle, waiting for a multiplier to be ready to accept a new input.
o A superscalar processor can issue several instructinos in the same cycle, but there's a fix small limit. Typical number is 2 -- 6. So, e.g., if the processor can issue maximum 4 instructions per cycle, the evaluation loop consists of 40 instructions, and the loop actually runs in close to 10 cycles per iteration, then instruction issue is the bottleneck.
The tricks discussed in this note are useful for finding an evaluation scheme where multiply latency isn't a bottleneck. But once a loop hits the limit on multiply throughput or instructions per cycle, other tricks are needed to optimize further. In particular, the postponed reduction has a cost in multiply throughput, since it needs some additional multiply instruction.
I think one should aim to hit the limit on multiply throughput; that one is hard to negotiate (it's possible to reduce the number of multiply instructions somewhat, by the Karatsuba trick, but due to the additional overhead, likely to be useful only on processors with particularly low multiply throughput).
Regards, /Niels
-- Niels Möller. PGP key CB4962D070D77D7FCB8BA36271D8F1FF368C6677. Internet email is subject to wholesale government surveillance.
nettle-bugs mailing list nettle-bugs@lists.lysator.liu.se http://lists.lysator.liu.se/mailman/listinfo/nettle-bugs
Maamoun TK maamoun.tk@googlemail.com writes:
On Sat, Jan 29, 2022 at 4:29 PM Niels Möller nisse@lysator.liu.se wrote:
** Interleaving **
The other approach, used in the recent powerpc gcm code, is to interleave multiple blocks. For simplicity, only consider 2-way interleaving here. The key thing is that if we expand te recurrency once, we get
R_j = K (M_j + K (R_{j-2} + M_{j-1})) = K M_j + K^2 (R_{j-2} + M_{j-1})
We get two field multiplications, but one of them, K M_j, is completely independent of previous blocks (R_{j-2}), and can be computed in parallel. It may add a cycle or so to accumulation latency, but we can do essentially twice as much work without making the critical path longer. We get 8 independent multiply instructions, and one dependent for the final folding.
But in radix 2^64, doesn't K^2 need to be represented in 3 words as described in a previous note so considering that we get 10 independent multiply instructions rather than 8?
For poly1305, yes. For gcm, I don't think so.
** Doing both **
It's possible to combine those two tricks. Processing of two blocks would then be an operation of the form
Z_2 B^2 + Z_1 B + Z_0 = (X_2 B^2 + X_1 B + X_0) K^2 + (Y_1 B + Y_0) K
Here, the Xs (three words) represent R_{j-2} + M_{j-1}, the Ys repreent M_{j-2}, and the Zs represents R_j, as three words (without final folding). We would need 10 independent multiples, one more than with plain interleaving, but critical path includes only one multiply latency.
Same matter here, if we consider 3 words for K^3 this sums up the number of independent multiples to 13.
For poly1305, yes. This seems to be an unavoidable cost of doing interleaving with poly1305: Structure of K is choosen to make multiplication particularly convenient, but that doesn't carry over to K^2. And we'll also have an X_3, of just a few bits. Unclear if it helps that most significant words are small; clearly multplying top 3 bits of X by top 3 bits of K^2 can use a short multiply (imul, on x86_64), but it seems to me that all other products need the more expensive 64x64->128 multiply instruction.
Also, if we have to use an extra word for poly1305, we could perhaps also explore going to radix 2^44, if hardware multipliers can support that.
Regards, /Niels
nettle-bugs@lists.lysator.liu.se
-
Maamoun TK -
nisse＠lysator.liu.se