Polynomial Arithmetic

[smithnormform]

Currently Cryptol has built-in polynomial arithmetic functions like pmult, pmod, and pdiv. However, pmult is restricted to polynomials with Bit coefficients. It would nice if Cryptol supported optimized primitive functions like pmult for polynomials with Z p or Integer coefficients.

[robdockins]

Related discussion here: #763

[brianhuffman]

A generic polynomial multiplication function for arbitrary rings can be defined in Cryptol like this:

polymul : {u, v, a} (fin u, fin v, Ring a) => [1 + u]a -> [1 + v]a -> [1 + u + v]a
polymul x y = map head (take`{v} sums) # last sums
  where
    rows : [1 + v][1 + u]a
    rows = [ [ xi * yi | xi <- x ] | yi <- y ]
    next : [1 + u]a -> [1 + u]a -> [1 + u]a
    next acc row = zipWith (+) (tail acc) (take`{u} row) # drop`{u} row
    sums : [1 + v][1 + u]a
    sums = scanl next (head rows) (tail rows)

Here's a correctness property to go with it:

polyeval : {n, a} (fin n, Ring a) => [n]a -> a -> a
polyeval p x = foldl (\y c -> x * y + c) zero p

polyeval_polymul :
  {u, v, a} (fin u, fin v, Ring a, Eq a) => [1 + u]a -> [1 + v]a -> a -> Bit
polyeval_polymul p q x =
  polyeval (polymul p q) x == polyeval p x * polyeval q x

@smithnormform: Would something like this be sufficient for your needs? If so, then no changes to Cryptol itself are necessary. If you need something with better asymptotic performance, then we could probably implement a more optimized version using a Karatsuba algorithm. (The code for this would be nicer if we implemented #701, but would probably still be possible with Cryptol as-is.) I think it's likely that we can get something that works well enough without adding any new primitives.

[smithnormform]

This is faster, but unfortunately not fast enough for the application. I think a Karatsuba algorithm is worth a shot! If the ring is Z or Z p then Fast Fourier Transform methods are often used for speedups, although these methods might cause issues with the SMT solvers. Whatever we come up with, it may be beneficial to place all of our results in a Polynomial package.

[weaversa]

I think built-in support for polynomials is necessary to support some PQ crypto such as NTRU and BFV. Both of those algorithms are elegant to spec in Cryptol excepting the polynomial arithmetic. Also, the solutions provided here are prohibitively slow on real parameters.

I like the suggestion of beefing-up the current polynomial arithmetic commands (pmod, pmul, pdiv, and adding padd) to work with Z types (and perhaps Bit for Z 2, which would need #1339). I'm also curious about @smithnormform's import polynomial suggestion.

(@brianhuffman this is an opportunity to take back <| and |>)

[robdockins]

Some time ago, I worked up a basic Karatsuba multiplier for bitvectors. It would probably not be too hard to rework it for polynomial rings.

https://github.com/GaloisInc/cryptol/blob/master/examples/Karatsuba.cry

It would be interesting to know if the asymptotic behavior or the constant factors are a bigger contributor at realistic parameter sizes.

[brianhuffman]

Here's a Karatsuba polynomial multiplier that I put together:

kmul : {u, v, a} (fin u, fin v, Ring a) => [1 + u]a -> [1 + v]a -> [1 + u + v]a
kmul x y = drop (karatsuba`{1 + max u v} (zero # x) (zero # y))

/** Multiply two big-endian polynomials, using a Karatsuba algorithm. */
karatsuba : {n, a} (fin n, Ring a) => [n]a -> [n]a -> [2 * n]a
karatsuba x y =
  if `n <= 50 then tail (polymul ([zero] # x) ([zero] # y)) else
  zipWith (+) (z1 # s # z1) (p1 # p2)
  where
    type n1 = n / 2
    type n2 = n - n / 2
    x1, y1 : [n1]a
    x2, y2 : [n2]a
    x1 # x2 = x
    y1 # y2 = y
    x', y' : [n2]a
    x' = padd x1 x2
    y' = padd y1 y2
    p1 : [2 * n1]a
    p1 = karatsuba x1 y1
    p2 : [2 * n2]a
    p2 = karatsuba x2 y2
    p' : [2 * n2]a
    p' = karatsuba x' y'
    s : [2 * n2]a
    s = psub (psub p' p1) p2
    z1 : [n1]a
    z1 = zero
    z2 : [n2]a
    z2 = zero

I chose a cutoff of 50 to switch to the naive polymul algorithm, based on a bit of experimentation. Performance doesn't seem to be too sensitive to the exact cutoff, but 10 is definitely too low, and 100 is too high. On my machine, it can multiply two polynomials of type [2000][64] in about 12.5 seconds, compared to over a minute for polymul. For type [4000][64], it's about 45 seconds, so it seems to be achieving the asymptotic behavior you'd expect for Karatsuba.

FFT-based algorithms could certainly be faster with better asymptotic behavior for very large polynomials, but unfortunately such algorithms can't be used polynomial multiplication over an arbitrary Ring, as FFT (or NTT, rather) algorithms only work on rings with suitable roots of unity.