Skip to content

trig recurrence accuracy (and speed) vs. precomputed table #55

@stevengj

Description

@stevengj

Note that your trigonometric recurrence has $O(\sqrt{n})$ root-mean-square error growth, which is much worse than the $O(\sqrt{\log n})$ mean error growth of the Cooley–Tukey algorithm (see e.g. the references here).

In particular, if I pull out your trig recurrence code and compare it to a function that simply calls cispi to compute $e^{-2\pi k /n}$ for $k = 0,\ldots,n-1$, the code (for BigFloat) looks like:

function trigtable(n)
    n *= 2
    θ=big"-2"/n
    wtemp = sinpi/2)
    wpr, wpi = -2wtemp^2, sinpi(θ)
    wr, wi = one(BigFloat), zero(BigFloat)
    w = [complex(wr, wi)]
    for m=1:2:n-2
        wr = (wtemp=wr)*wpr-wi*wpi+wr
        wi = wi*wpr+wtemp*wpi+wi
        push!(w, complex(wr, wi))
    end
    return w
end

# a more accurate alternative:
trigtable2(n) = cispi.(.- (0:n-1) ./ BigInt(n)) # simple direct-call algorithm

For the default 256-bit BigFloat precision ($\varepsilon \approx 10^{-77}$), one can see the $O(\sqrt{n})$ error growth in the relative $L_2$ error $\Vert x - y \Vert_2 / \Vert y \Vert_2$ between these two methods:

Image

The easiest accurate alternative is to use a precomputed trig table such as that returned by trigtable2(n), which could be stored in a "plan" object for efficient repeated FFTs.

Maybe it's not such a big concern if the precision is large and $n$ is small, but it might be nice to have the option for a more accurate result. Also, precomputing the trig table might be faster if you re-use it for a lot of FFTs.

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions