Towards an iterative algorithm for the DTF

We first change the intermediate polynomials. Let n = 2^r mathend000# be a power of 2 mathend000#, let $\omega$ $\in$ R mathend000# be a primitive n mathend000#-th root of unity and A(x) = $\Sigma_{{{0 \, \leq \, i \, < \, n}}}^{{}}$ a_i xⁱ mathend000#. We define

A[0](x)	=	a0 + a2x + a4x2 + ... + an-2xn/2-1
A[1](x)	=	a1 + a3x + a5x2 + ... + an-1xn/2-1

(55)

such that we have

A(x)  =  A^[0](x²) + x A^[1](x²). (56)

Hence evaluating A mathend000# at 1,$\omega$,$\omega^{{2}}_{{}}$,...,$\omega^{{{n-1}}}_{{}}$ mathend000# is reduced to evaluate A^[0] mathend000# and A^[1] mathend000# at 1,$\omega^{{2}}_{{}}$,$\omega^{{4}}_{{}}$,...,$\omega^{{{2n-2}}}_{{}}$ mathend000#. But this second sequence of powers of $\omega$ mathend000# consists of two identical sequences of length n/2 mathend000#. Indeed, for every integer k mathend000# we have

$$\omega^{{{k+ n/2}}}_{{}} \omega^{{{k}}}_{{}}$$ (57)

Since $\omega^{{{k}}}_{{}}$(1 - $\omega^{{{n/2}}}_{{}}$) $\neq$ 0 mathend000# holds, we obtain

$$\omega^{{{k+ n/2}}}_{{}} \omega^{{{k}}}_{{}}$$ (58)

Hence we only need to evaluate A^[0] mathend000# and A^[1] mathend000# at

$$\omega^{{2}}_{{}} \omega^{{4}}_{{}} \omega^{{{n-2}}}_{{}}$$ (59)

This provides us directly with the values of A mathend000# at

$$\omega \omega^{{2}}_{{}} \omega^{{{n/2-1}}}_{{}}$$ (60)

and the values of A mathend000# at

$$\omega^{{{n/2}}}_{{}} \omega^{{{n/2+1}}}_{{}} \omega^{{{n-1}}}_{{}}$$ (61)

by using the fact these numbers are respectively equal to

$$\omega \omega^{{2}}_{{}} \omega^{{{n/2-1}}}_{{}}$$ (62)

To summarize, for 0 $\leq$ i $\leq$ n/2 - 1 mathend000# we have

A($$\omega^{{i}}_{{}}$$)	=	A[0]($$\omega^{{{2i}}}_{{}}$$) + $$\omega^{{i}}_{{}}$$ A[1]($$\omega^{{{2i}}}_{{}}$$)
A($$\omega^{{{i+ n/2}}}_{{}}$$)	=	A[0]($$\omega^{{{2i}}}_{{}}$$) - $$\omega^{{i}}_{{}}$$ A[1]($$\omega^{{{2i}}}_{{}}$$)

(63)

Algorithm 3  

mathend000#

The validity of Algorithm 3 follows from the two following observations.

• Formulas (63).
• If $\omega$ mathend000# is a primitive 2^s mathend000#-root of unity, with s > 1 mathend000#, then $\omega^{{2}}_{{}}$ mathend000# is a primitive 2^s-1 mathend000#-root of unity.

Observe that the recursive calls of DFT(n, A,$\Omega$) mathend000# defines an ordering of the coefficients of A mathend000# shown on Figure 2 for n = 8 mathend000#. Let us call this ordering the DFT ordering of A mathend000#.

Observe that if the input array A mathend000# is sorted w.r.t. this DFT ordering, one can describe easily the recursive calls with a binary tree. If this tree is traversed bottom-up, from left to right, we are led to an iterative algorithm working in place! For instance, for n = 8 mathend000#, assuming that the input array A mathend000# is

We obtain the desired result in A mathend000# using the following straigth-line program

1. a mathend000# := a₀ mathend000#
2. b mathend000# := a₁ mathend000#
3. a₀ mathend000# := a + b mathend000#
4. a₁ mathend000# := a - b mathend000#
5. a mathend000# := a₂ mathend000#
6. b mathend000# := a₃ mathend000#
7. a₂ mathend000# := a + b mathend000#
8. a₃ mathend000# := a - b mathend000#
9. a mathend000# := a₄ mathend000#
10. b mathend000# := a₅ mathend000#
11. a₄ mathend000# := a + b mathend000#
12. a₅ mathend000# := a - b mathend000#
13. a mathend000# := a₆ mathend000#
14. b mathend000# := a₇ mathend000#
15. a₆ mathend000# := a + b mathend000#
16. a₇ mathend000# := a - b mathend000#
17. a mathend000# := a₀ mathend000#
18. b mathend000# := a₁ mathend000#
19. c mathend000# := a₂ mathend000#
20. d mathend000# := a₃ mathend000#
21. a₀ mathend000# := a + c mathend000#
22. a₂ mathend000# := a - c mathend000#
23. a₁ mathend000# := b + $\omega$d mathend000#
24. a₃ mathend000# := b - $\omega$d mathend000#
25. a mathend000# := a₄ mathend000#
26. b mathend000# := a₅ mathend000#
27. c mathend000# := a₆ mathend000#
28. d mathend000# := a₇ mathend000#
29. a₀ mathend000# := a + c mathend000#
30. a₂ mathend000# := a - c mathend000#
31. a₁ mathend000# := b + $\omega$d mathend000#
32. a₃ mathend000# := b - $\omega$d mathend000#
33. a mathend000# := a₀ mathend000#
34. b mathend000# := a₁ mathend000#
35. c mathend000# := a₂ mathend000#
36. d mathend000# := a₃ mathend000#
37. e mathend000# := a₄ mathend000#
38. f mathend000# := a₅ mathend000#
39. g mathend000# := a₆ mathend000#
40. h mathend000# := a₇ mathend000#
41. a₀ mathend000# := a + e mathend000#
42. a₁ mathend000# := b + $\omega$f mathend000#
43. a₂ mathend000# := c + $\omega^{{2}}_{{}}$f mathend000#
44. a₃ mathend000# := d + $\omega^{{3}}_{{}}$h mathend000#
45. a₄ mathend000# := a + e mathend000#
46. a₅ mathend000# := b - $\omega$f mathend000#
47. a₆ mathend000# := c - $\omega^{{2}}_{{}}$f mathend000#
48. a₇ mathend000# := d - $\omega^{{3}}_{{}}$h mathend000#

Figure 2: Recursive calls for n = 8 mathend000#.

Algorithm 4  

mathend000#

Remark 8   To get the DFT ordering for A mathend000#, that is [a₀, a₄, a₂, a₆, a₁, a₅, a₃, a₇] mathend000#, observe that this ordering changes

000, 001, 010, 011, 100, 101, 011, 111 (64)

to

000, 100, 010, 110, 001, 101, 110, 111 (65)

Remark 9   The ring $\mbox{${\mathbb Z}$}$ mathend000# does not have any interesting roots of unity. In order to take advantage of the FFT-based multiplication one could use a modular approach as follows.

• Let f mathend000# and g mathend000# be two univariate polynomials in $\mbox{${\mathbb Z}$}$[x] mathend000#. We can assume that all their coefficients are positive. Otherwise, we can restrict to this case by setting

  f  =  f⁺ - f^-    and    g  =  g⁺ - g^- (66)

  
  
  where f⁺ mathend000#, f^- mathend000#, g⁺ mathend000#, g^- mathend000# are polynomials with positive coefficients and writing

  f g  =  f⁺g⁺ + f^-g^- - f⁺g^- - f^-g⁺ (67)

  
  
  and computing the products f⁺g⁺ mathend000#, f^-g^- mathend000#, f⁺g^- mathend000# and f^-g⁺ mathend000# one after another.
• Let | | f | | mathend000# and | | g | | mathend000# be the max-norms of f mathend000# and g mathend000#. Let d mathend000# be the degree of f g mathend000#. For k = 0 ^... d mathend000# a bound for the coefficient c_k mathend000# of x^k mathend000# in f g mathend000# is given by

  $$\leq$$ (68)

  
  

• Let M mathend000# be the biggest value of these bounds. Let p₁,..., p_r mathend000# be primes such
  • 2^d mathend000# divides each p₁ -1,..., p_r - 1 mathend000#
  • the product of p₁,..., p_r mathend000# is bigger than 2 M mathend000#.
• Then we can compute f g mathend000# in each of the $\mbox{${\mathbb Z}$}$/p_i$\mbox{${\mathbb Z}$}$[x] mathend000# using the FFT-based multiplication.
• Finally we can recombine these modular product by applying the CRA to each coefficient c_k mathend000# (of x^k mathend000#) in f g mathend000#.