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
    No fields configured for issues without a type.

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions