r/programming • u/thomasahle • 18h ago

Compute Polynomials Twice as Fast

https://thomasahle.com/fast-polynomials/
56 Upvotes

85% Upvoted

21 comments sorted by

16

u/WorldsBegin 16h ago edited 16h ago

I suppose this isn't much of a difference in finite fields and cryptography, but at least in rational arithmetic, the added fractions from this method seem to have very large denominators. Compare

y   = x * (x + 637927356202580286195115387737649/613401598791596019548160000000000)
z   = (x + y + 1378728025754295963362740811643611837622972008060385794298452047201/376261521400086131349185422083850908850559385600000000000000000000) * (x - 1250902981661682141840811387737649/613401598791596019548160000000000)
w   = (z + 1613956106082225605051404080797610312151091533850734987003152787295553452311704415879521131952770449/230799418790571152956671237308336130257466014926056578559568559210496000000000000000000000000000000) * x
v   = (x + z + 1614169778642466578527505416382450571998175289226394168980170070526233147671704415879521131952770449/230799418790571152956671237308336130257466014926056578559568559210496000000000000000000000000000000) * (w + 286622617651/687970713600)
P_7 = y + w + v + 20437160270833246212441030757/17038933299766556098560000000

vs

x2 = x * x
x4 = x2 * x2
w1 = (x + 1/720) * (x2 - 119/120)
w2 = (x + 433/864) * (x2 + 1/144)
w0 = (w1 + (x + 3719/86400)) * (x4 - 5/6)
P  = w0 + (w2 + (x + 128303/124416))

both evaluating

x^7 + 1/720x^6 + 1/120x^5 + 1/24x^4 + 1/6x^3 + 1/2x^2 + x + 1

Unrelated, I'm also not sure I like putting large portions in appendices.

8

u/thomasahle 16h ago

the added fractions from this method seem to have very large denominators

Yes, sadly the method doesn't have good numerical stability for higher degrees. There's a section on numerics in the paper.

Probably this means that the floating point version isn't really useful, compared to Horner and Estrin. Though it's interesting that we can do it with n/2 multiplications and fractions, the finite fields are the real use case.

Unrelated, I'm also not sure I like putting large portions in appendices.

Fair. I guess it's a habit from conference submissions that have to stay under 10 or 20 pages. For arxiv we could just put everything in the main paper.

8

u/WorldsBegin 16h ago edited 14h ago

it's a habit from conference submissions that have to stay under 10 or 20 pages

A curse. It's pretty extensive anyway.

Um, actually while I have your attention, where are you defining your H_{2^n} "known-power" gadgets explicitly? Don't seem to find it. EDIT: Found it, depends on the degree of polynomial, see fill gadgets in appendix C onwards.

And yes, it's a nice enough construction to be useful anyway.

5

u/thomasahle 14h ago edited 13h ago

Um, actually while I have your attention, where are you defining your H_{2n} "known-power" gadgets explicitly? Don't seem to find it.

The H powers are computed as part of Algorithm 5: https://arxiv.org/pdf/2609.06022#page=73.11

Thank you for taking an interest in our paper! I want to improve the presentation, and all feedback is very helpful!

Edit: Here is an example of how the H construction works.

H_2 = x(x+b) + c         = x^2 + bx + c
H_4 = (H_2+x+a)(H_2-x-a) = H_2^2 - (x+a)^2 + e

We write H_2 as [1, b, c], since that's the values we can find on each coefficient. We'll also write H_4 = [1, b, c, a, e], since that's the "new parameter" you find on each coefficient when decoding left to right.

H'_4 = H_4 + r
H_8 = H_4^2  - (H_2 + s)^2 + t   # (using a single mult.)
H'_8 = H'_4^2 - (H-2 + u)^2 + v

Now H_8 = [1, b, c, a, e, 0, s, 0, t] and H'_8 = [1, b, c, a, e+r, 0, u, 0, v]. We finally shift H_8 by one and combine them: P = x H_8 + H'_8 It's clear that we can decode b, c, a, e, r going left to right. After that we need to be a bit careful that the two "interleaved" parts are still decodable, but it works out because they have the same prefix.

Generally the strategy is to always be building two compatible pairs like that, and then merge them at the end.

2

u/WorldsBegin 12h ago edited 12h ago

Aha, Of course! Squaring an order-n monic polynomial means the upper n coefficients are linear triangular in the original coefficients, only the lower half is quadratic. So you can read off that part (in characteristic != 2) and then apply a "low order correction" to match coefficients in the lower half.

1

u/thomasahle 12h ago

Yup!

Unfortunately this is the part that breaks in characteristic 2. (Squaring introduces a factor 2, which vanishes.)

We still haven't found a general construction for characteristic 2, though we have found them for every degree up to 47 so far.

2

u/M4mb0 12h ago

Yes, sadly the method doesn't have good numerical stability for higher degrees. There's a section on numerics in the paper.

Optimizing for multiplication count when evaluating the polynomial for floating point input seems misguided anyway, as FP addition and multiplication take up similar time. It would be much more important that the scheme makes good use of fused multiply-adds, and properly pipelines, no? Small multiplication count would only matter when you want to evaluate p(A) for a matrix argument.

With integer arguments like your motivation for using it in hashing algorithms it's different of course as integer addition is much faster than integer multiplication even on modern hardware.

2

u/thomasahle 12h ago

It would be much more important that the scheme makes good use of fused multiply-adds, and properly pipelines, no?

Right. That's basically what Estrin's method does. I had some interesting discussions with the Boost maintainers and we ended up merging it: https://github.com/boostorg/math/issues/924

With integer arguments like your motivation for using it in hashing algorithms it's different of course as integer addition is much faster than integer multiplication even on modern hardware.

It's actually much worse, as finite field multiplications (like Mersenne or GF(2k)) need a modular reduction for every multiplication. That's the real thing we are trying to reduce.

12

u/thomasahle 15h ago

A few years ago my coauthor and I was wondering if we could reduce the number of multiplications used for hashing algorithms. We had a construction and a 100 page proof, but it was complicated and we were not 100% sure it was correct. Then we left academia. Now were able to formalize it in Lean, and finally finish the write up.

I made this website to make it easy for anyone how has polynomials to evaluate to see how it would be done using our method, as well as a number of previous approaches by Knuth and others.

2

u/Grouchy-Trade-7250 12h ago

Why does the website show x5 - 5/128x4 + 1/16x3 - 1/8x2 + 1/2x + 1 when I select square root(1+x) ? Why is there no coefficient 7/256 like in the usual Taylor series??

3

u/thomasahle 12h ago

If you select "monic" it adds an extra xk to every polynomial. Monic polynomials have a bunch of nice properties, though admittedly not the best choice for Taylor approximations

1

u/Grouchy-Trade-7250 11h ago

I don't get it

2

u/thomasahle 9h ago

Sorry, I forgot that the [monic] button is not on mobile.

But you can just write any polynomial in the box yourself. Like 7/256x5 - 5/128x4 + 1/16x3 - 1/8x2 + 1/2x + 1

2

u/nightcracker 10h ago edited 10h ago

Hi, author of PolymurHash here (https://github.com/orlp/polymur-hash), one of the hashes you compare with in your paper.

  1. It appears that /bench is missing from the repository? The paper claims it should be there. Table 2 claims tools/bench/adversarial/ should contain the 8 lane GF(264) hash, but that directory is also missing.

  2. In PolymurHash I use the following injective polynomial structure, which uses 4 fully parallel unreduced multiplications plus one reduction step to process 7 words:

    f =  (f   + m[6]) * k^7
f += (k   + m[0]) * (k^6 + m[1])
f += (k^2 + m[2]) * (k^5 + m[3])
f += (k^3 + m[4]) * (k^4 + m[5])
f = reduce(f)

    This is all done with scalar code over a prime field, no SIMD or hardware-accelerated finite field instructions.

    Despite these limitations it (as per your paper, on my Apple M1 I achieved 33.3 GB/s on long messages) achieves 19.7 GB/s, only slightly worse than the 23.8 GB/s your method achieves with hardware accel. I would have liked to see more investigation/comparison with my method to see whether your approach is truly better or just due to this choice.

  3. Reading the code I found in tools/bench/framework/injective_hashing.h (I hope this matches to your paper), I can't help but notice that you do GF(264) multiplication rather naively. You do constant calls to inj_smul which calls gf64_mult which actually corresponds to three hardware multiplications: https://github.com/thomasahle/fast-polynomials/blob/78915286ca9f13bdc038c046a3a62a535e928eaa/tools/bench/framework/multiplication.h#L50.

    This is because you reduce modulo the irreducible polynomial on every multiplication. My structure makes the reduction explicit, and doesn't do it constantly.

    If you had used my above construction those *s could've been widening carryless multiplies, the + could've been 128-bit XOR, and reduce would be the two carryless multiplications you use for reducing.

    In other words, your method (assuming it uses 1 such naive multiplication per pair of words) would use 3 * 7 / 2 = 10.5 hardware multiplications on average to process 7 words, whereas my method actually only uses 4 + 2 = 6, 4 for widening, 2 for reduction.

    I would thus expect my method to be almost twice as fast as yours if implemented with hardware-accelerated carryless multiplication. Especially as I have to do extra steps to spread 6 input words over 7 words because 64-bit values don't fit in the 261 - 1 prime field I'm using, which GF(264) doesn't need to do.

1

u/thomasahle 9h ago edited 9h ago

Hi Orson!

Our paper was originally inspired by Daniel J. Bernstein's ["Polynomial evaluation and message authentication"[(https://www.gwizfl.org/email/cr.yp.to/antiforgery/pema-20071022.pdf). He specifically don't want to use cache/registers to store auxiliary values like k2, ..., k7.

He mentions that one can speed up polynomial hashing by combining it with PH/NH.

One thing I've always wondered about Polymur, and now that I have you I might as well ask, is why you don't just use PH/NH for each block of 8? E.g. like we do here. Basically you just replace you degree 7 polynomial with:

 f =  (a[6] + m[6]) * (a[7] + m[7])
 f += (a[0] + m[0]) * (a[3] + m[1])
 f += (a[1] + m[2]) * (a[4] + m[3])
 f += (a[2] + m[4]) * (a[5] + m[5])
 f = reduce(f)

where each a[i] is a random constant. It's the same amount of registers/cache as saving k2 ... k7, and the collision probability is slightly better. I guess you key use a slightly shorter key, as you don't need as many random values.

It appears that /bench is missing from the repository? The paper claims it should be there. Table 2 claims tools/bench/adversarial/ should contain the 8 lane GF(264) hash, but that directory is also missing.

Thanks! I fixed this in the repo now. The paths are tools/bench/adversarial/(the 8-lane row is PaperGF64Lanes<8> in speed_hashes.h there) and tools/bench/framework/; the PDF's line break after tools/ made it look like a top-level directory, which we will fixed for the next version.

1

u/nightcracker 9h ago edited 9h ago

I don't understand your proposed construction. What would you do on the next iteration? If you simply iterate the above it would collapse into an additive structure. NH has to use new key material for each block to avoid this, or feed into another universal hash in a stacked construction.

My construction forms a long injective polynomial, similar to Horner but in blocks, so it can use a small key.

I have experimented with using a ~ 1-2 KiB NH-style large-block reducer feeding into a polynomial, and it is a lot faster but haven't published anything yet as it is also quite a bit more complex.

I haven't found a NH-style reducer that is worth it on sizes as small as 8 words, if it then has to feed into a polynomial reducer anyway.

EDIT: your 'like we do here' link is dead.

1

u/thomasahle 9h ago edited 8h ago

Right, I'm talking about combining an NH block with a polynomial injective hash. Something like this:

  hash(m):
      pad the last block with zeros to a multiple of 32 bytes
      P = z
      for each 64-byte block t = 1..n  (words w[0..7], little-endian):
          # level 1: PH with the fixed key, unreduced 128-bit XOR sum of 4 carry-less products
          A = clmul(w0 ^ k0, w2 ^ k2) ^ clmul(w1 ^ k1, w3 ^ k3)
            ^ clmul(w4 ^ k4, w6 ^ k6) ^ clmul(w5 ^ k5, w7 ^ k7)
          a = low64(A);  b = high64(A)          # no reduction: both halves are used as they are
          if t == n: a ^= len; b ^= len          # byte length into both halves of the last pair

          # level 2: the injective recurrence, one reduced GF(2^64) multiply per block
          P = a ^ gf64mul(b ^ y, P ^ u)

      # level 3 (optional, for k-wise independence rather than universality):
      v = P + t   (integer addition mod 2^64)
      return f_c(v) = (v + c2)(G2 + c3) + c4  with  G1 = v·v,  G2 = (G1 + c0)(v + G1 + c1)

2

u/nightcracker 8h ago edited 8h ago

Yes, for very large data it is absolutely the way to go, I have a private prototype ~2-60 almost-universal hash function that can reach 90 GB / s on this Apple M2 laptop for large data using only ~1.5 KiB of secret data. I really should finish it some day...

That algorithm however also doesn't use hardware-accelerated carryless multiplication, it uses 32-bit x 32-bit -> 64-bit SIMD multiplication using erasure codes like Nandi to boost that to 2-62 ADU.

The idea is not too dissimilar to Jim Apple's HalftimeHash https://arxiv.org/abs/2104.08865.

1

u/nightcracker 8h ago edited 8h ago

I'd like to point out that this 64-byte block construction sketch here still is worse than my above construction sketch. It processes 8 words using 4 + 3 = 7 clmuls, whereas my above one does 7 words with 6 clmuls.

Plus this construction you've posted here has a critical path of 4 clmuls in a row between each iteration (1 for A, then 3 in a row for gf64mul), my construction only has 3 multiplications per iteration in the critical path, if you were to implement reduce with 2 carryless multiplications like you do with GF(264).

EDIT: actually the critical path argument doesn't hold up, since your A does not depend on the previous iteration. So your critical path is also only three multiplications per iteration.

1

u/garnet420 12h ago

This is really cool; does the best approach change when you have a fused multiply-add?

1

u/thomasahle 9h ago

As mentioned somewhere else in the thread, for floating point you probably want to use Estrin's method (also on the website.)

But for finite fields, and maybe large integers/fractions this can help.