0

I have been studying this topic in an Algorithms textbook.

The clever usage of the complex roots of unity seems to be mathematically working. However, I do not understand how one could actually represent this in a computer.

I can think of two things:

  • Use the real/imaginary decomposition to represent the complex numbers. But this means using floats, which means I open up my algorithm to numerical error and I would loose out precision even if I want to multiply two polynomials with integer coefficients.
  • Represent exp(i 2pi/n) as n. So, I'd eventually get a tuple in omega, and if I have to keep it in this form, I'd essentially be doing polynomial multiplication in omega again, taking us back to square one.

I'd really like to see an implementation of this algorithm in a familiar programming language.

  • "Numerical Recipes in C" has an FFT implementation in all its glory. NumPy/SciPy would have Python versions. There are Java open source libraries. Find the language you know and download the source. – duffymo Aug 10 '17 at 15:42
  • Perhaps https://math.stackexchange.com/questions/764727/concrete-fft-polynomial-multiplication-example ? – Severin Pappadeux Aug 10 '17 at 23:18

1 Answers1

3

Indeed as you identify, the roots of unity are typically not nice numbers that can be stored well in a computer. Since the numerical error is small, if you know the output should be integers, rounding usually produces the right result.

If you don't want to (or cannot) rely on that, an exact option is the Number Theoretic Transform. It substitutes the roots of unity in the complex plane with roots of unity in a finite field ℤ/pℤ where p is a suitable prime. p has to be large enough for all the necessary roots to exist, and the efficiency is affected by properties of p. If you choose a Fermat prime then the roots of unity have convenient forms and there is a trick to do reduction modulo p more efficiently than usual. That is all exact integer arithmetic and the values stay small, so there is no problem implementing it in a computer.

That technique is used in the Schönhage–Strassen algorithm so you can look up the specifics there.

harold
  • 61,398
  • 6
  • 86
  • 164
  • $p$ does not need to be a prime, and indeed Schönhage does not discuss primality. Only the form to look like Fermat primes is important, as that readily supplies the equivalent of the root of unity. – Lutz Lehmann Aug 13 '17 at 09:44
  • Yes some composites are fine, all we really need is for the roots to exist. But I think such details are going needlessly further down the rabbit hole here. – harold Aug 13 '17 at 13:03
  • Yes NTT is possible way for integers see [Modular arithmetics and NTT (finite field DFT) optimizations](https://stackoverflow.com/q/18577076/2521214) for C++ example for strings multiplication. – Spektre Aug 18 '17 at 08:09