What algorithm is used to calculate $\exp(i 2 \pi/1024)$ in FFT?

Question

For the fast Fourier transform algorithm, a prerequisite is calculating quantities like $\exp(i 2\pi/ 2^n)$. How is this done?

Is there any specific algorithm tailored for this kind of special angles?

Calcuting $\exp(i\theta)$ at the precision of a float isn't complicated with the Taylor series of $\exp(z)$. — reuns, Dec 05 '17 at 12:13

Somos · Answer 1 · 2017-12-05T13:38:52.637

1

Let $a_n:=\exp(i2\pi/2^n).$ Then $a_1=-1,a_2=i,$ and $a_{n+1}=\sqrt{a_n}.$ There are good algorithms to compute square roots of complex numbers, especially roots of unity, and they can be used here. For example, if $x+yi=(u+vi)^2, x^2+y^2=1,$ then $u=\sqrt{(1+x)/2},v=y/(2u).$

edited Dec 05 '17 at 13:38

answered Dec 05 '17 at 13:04

Somos

35,251
3
30
76

1

Don't you think it starts with the $6$-terms Taylor polynomial, then improves the result using an iterative algorithm for the square root ? – reuns Dec 05 '17 at 13:13

user5713492 · Answer 2 · 2017-12-08T23:02:08.537

FFT algorithms don't necessarily do any complex multiplications because they imply unpacking and repacking SIMD data. A table of $\cos(2\pi ik/2^N)$ for $0<k<2^{N-2}$ is sufficient to perform an FFT. My recollection is that it's best to maintain the table in bit-reversed order. Whether it's efficient to keep that much data around depends on the size of the problem.

EDIT: I changed my minmax polynomial program to use Hormer's method instead of Chebyshev polynomials. Not only faster but a little more accurate with worst-case error of $1.12$ ulps for $\cos\left(\frac{2\pi\cdot176}{1024}\right)$.

! poly2.f90 -- Test minmax formulas for N = 2**10
program poly2
   implicit none
   integer, parameter :: dp = kind(1.0d0)
   integer, parameter :: qp = selected_real_kind(33,4931)
   integer, parameter :: wp = dp
   integer, parameter :: logN = 10
   integer, parameter :: N = 2**logN
   real(qp), parameter :: pi = 4*atan(1.0_qp)
! Parameters for 8 term sin and cos approximations
! cos(pi*x/4) = sum([(params(2*j)*x**(2*j),j=0,7)])
! sin(pi*x/4) = sum([(params(2*j+1)*x**(2*j+1),j=0,7)])
! -1 <= x <= 1
   real(dp), parameter :: params(0:15) = [ &
           1.0000000000000000_wp, &
          0.78539816339744828_wp, &
         -0.30842513753404244_wp, &
         -0.80745512188280785E-01_wp, &
          0.15854344243815419E-01_wp, &
          0.24903945701927120E-02_wp, &
         -0.32599188692673781E-03_wp, &
         -0.36576204182126918E-04_wp, &
          0.35908604460283630E-05_wp, &
          0.31336168886997841E-06_wp, &
         -0.24611364034255195E-07_wp, &
         -0.17572473559409242E-08_wp, &
          0.11500512072041326E-09_wp, &
          0.69481110892351598E-11_wp, &
         -0.38581902594191821E-12_wp, &
         -0.20214433162524027E-13_wp]
! Bit reversal of indices
   integer args(0:2**(logN-2)-1)
   integer j, L
! Argument and squares
   real(dp) Z(0:2**(logN-3)),W(0:2**(logN-3)),T(1:2**(logN-3)-1)
   real(dp) C(0:2**(logN-2)-1)
   real(dp) ulps(2**(logN-2)-1)

! Compute bit reversed array
   args(0) = 0
   do L = 1, logN-2
      args(0:2**(L-1)-1) = 2*args(0:2**(L-1)-1)
      args(2**(L-1):2**L-1) = 1+args(0:2**(L-1)-1)
   end do
! Compute arguments to trig functions
   Z = [(j,j=0,2**(logN-3))]*8.0_dp/N
   W = Z**2
   T = W(2**(logN-3)-1:1:-1)
   C = 0
! Add up series via Horner's method
   do L = ubound(params,1)-1, 0, -2
      C(0:2**(logN-3)) = C(0:2**(logN-3))*W+params(L)
      C(2**(logN-3)+1:2**(logN-2)-1) = C(2**(logN-3)+1:2**(logN-2)-1)*T+params(L+1)
   end do
! Need last factor for sines
      C(2**(logN-3)+1:2**(logN-2)-1) = C(2**(logN-3)+1:2**(logN-2)-1)*Z(2**(logN-3)-1:1:-1)
! Bit reverse array of cosines
   C = C(args)
! Get error in ulps
   do j = 1, 2**(logN-2)-1
      ulps(j) = abs(C(j)-cos(2*pi*args(j)/N))/spacing(C(j))
   end do
! Print out index of max relative error, k value, and error
   write(*,*) maxloc(ulps),args(maxloc(ulps,1)),maxval(ulps)
end program poly2

One of the big problems with computing trig functions via minmax polynomials is the reduction of the argument to the domain of validity of the polynomial, but in producing a table like this the argument reduction has been performed in advance so the polynomials can be evaluated in parallel using the SIMD units and pipelining.

Another method is to use angle sum formulas considering that the angles are in arithmetic progression, but there is the potential for roundoff errors in this method.

EDIT: Here is the arithmetic progression method. As anticipated, the errors really blew up, growing to $34.56$ ulps at $\cos\left(\frac{2\pi\cdot128}{1024}\right)$.

! arithmetic.f90 -- Test arithmetic progression formulas for N = 2**10
program arithmetic
   implicit none
   integer, parameter :: dp = kind(1.0d0)
   integer, parameter :: qp = selected_real_kind(33,4931)
   integer, parameter :: logN = 10
   integer, parameter :: N = 2**logN
   integer, parameter :: kmax = 2**(logN-2)-1
   real(dp), parameter :: pi = 4*atan(1.0_dp)
   real(qp), parameter :: qpi = 4*atan(1.0_qp)
! sin(2*pi/N), cos(2*pi/N)
   real(dp) sin_theta, cos_theta
   real(dp) x, xsquared
   real(dp) C(1:kmax)
   integer i
   real(dp) ulps(1:kmax)
! Compute trig functions of smalles angles via Taylor series
   x = 2*pi*1/N
   xsquared = x**2
   sin_theta = 1
   cos_theta = 1
   do i = 4, 0, -2
      cos_theta = 1-cos_theta*xsquared/((i+1)*(i+2))
      sin_theta = 1-sin_theta*xsquared/((i+2)*(i+3))
   end do
   sin_theta = sin_theta*x
! Get the rest of the cosines as
! cos((i+1)*x) = cos(i*x)*cos(x)-sin(i*x)*sin(x)
! sin((i+1)*x) = sin(i*x)*cos(x)+cos(i*x)*sin(x)
   C(1) = cos_theta
   C(kmax) = sin_theta
   do i = 2,2**(logN-3)
      C(kmax+1-i) = C(kmax+2-i)*cos_theta+C(i-1)*sin_theta
      C(i) = C(i-1)*cos_theta-C(kmax+2-i)*sin_theta
   end do
! Get error in ulps
   do i = 1, kmax
      ulps(i) = abs(C(i)-cos(2*qpi*i/N))/spacing(C(i))
   end do
! Print out index of max relative error, and error
   write(*,*) maxloc(ulps),maxval(ulps)
end program arithmetic

If one doesn't have another algorithm coded, the half-angle formulas can produce the necessary cosines at the cost of a division or square root per table entry.

EDIT: Here is my program using the half-angle formulas. I would have thought it to be quite accurate but in fact it was off by $1.44$ ulps for $\cos\left(\frac{2\pi\cdot187}{1024}\right)$.

! sqrts.f90 -- Test half-angle formulas for N = 2**10
program sqrts
   implicit none
   integer, parameter :: dp = kind(1.0d0)
   integer, parameter :: qp = selected_real_kind(33,4931)
   integer, parameter :: logN = 10
   integer, parameter :: N = 2**logN
   integer, parameter :: kmax = 2**(logN-2)-1
   integer L, k
   real(qp), parameter :: pi = 4*atan(1.0_qp)
   real(dp) ulps(kmax)
! Array of cosines in bit-reversed order
   real(dp) C(kmax)
! Get cos(pi/4)
   C(1) = sqrt(0.5_dp)
   do L = 1, logN-3
! cos(theta/2) = sqrt((1+cos(theta))/2)
      C(2**L:2**(L+1)-2:2) = sqrt((1+C(2**(L-1):2**L-1))/2)
! sin(theta/2) = sin(theta)/(2*cos(theta))
      C(2**L+1:2**(L+1)-1:2) = C(2**(L-1):2**L-1)/(2*C(2**(L+1)-2:2**L:-2))
   end do
! Get errors in ulps
   do k = 1, kmax
      ulps(k) = abs(C(k)-cos(2*pi*bitrev(k,logN-2)/N))/spacing(C(k))
   end do
! Print out index of max relative error, k value, and error
   write(*,*) maxloc(ulps),bitrev(maxloc(ulps,1),logN-2),maxval(ulps)
   contains
! Computes bit reversal of nbits-long number num.
      function bitrev(num,nbits)
         integer bitrev
         integer num
         integer nbits
         integer i
         bitrev = sum([(2**(nbits-1-i)*ibits(num,i,1),i=0,nbits-1)])
      end function bitrev
end program sqrts

What algorithm is used to calculate $\exp(i 2 \pi/1024)$ in FFT?

2 Answers2