# FastTransforms.jl

`FastTransforms.jl`

allows the user to conveniently work with orthogonal polynomials with degrees well into the millions.

This package provides a Julia wrapper for the C library of the same name. Additionally, all three types of nonuniform fast Fourier transforms are available, as well as the Padua transform.

## Installation

Installation, which uses BinaryBuilder for all of Julia's supported platforms (in particular Sandybridge Intel processors and beyond), may be as straightforward as:

pkg> add FastTransforms
julia> using FastTransforms, LinearAlgebra

## Fast orthogonal polynomial transforms

The orthogonal polynomial transforms are listed in `FastTransforms.Transforms`

or `FastTransforms.kind2string.(instances(FastTransforms.Transforms))`

. Univariate transforms may be planned with the standard normalization or with orthonormalization. For multivariate transforms, the standard normalization may be too severe for floating-point computations, so it is omitted. Here are two examples:

### The Chebyshev--Legendre transform

julia> c = rand(8192);
julia> leg2cheb(c);
julia> cheb2leg(c);
julia> norm(cheb2leg(leg2cheb(c; normcheb=true); normcheb=true)-c)/norm(c)
1.1866591414786334e-14

The implementation separates pre-computation into an `FTPlan`

. This type is constructed with either `plan_leg2cheb`

or `plan_cheb2leg`

. Let's see how much faster it is if we pre-compute.

julia> p1 = plan_leg2cheb(c);
julia> p2 = plan_cheb2leg(c);
julia> @time leg2cheb(c);
0.433938 seconds (9 allocations: 64.641 KiB)
julia> @time p1*c;
0.005713 seconds (8 allocations: 64.594 KiB)
julia> @time cheb2leg(c);
0.423865 seconds (9 allocations: 64.641 KiB)
julia> @time p2*c;
0.005829 seconds (8 allocations: 64.594 KiB)

Furthermore, for orthogonal polynomial connection problems that are degree-preserving, we should expect to be able to apply the transforms in-place:

julia> lmul!(p1, c);
julia> lmul!(p2, c);
julia> ldiv!(p1, c);
julia> ldiv!(p2, c);

### The spherical harmonic transform

Let `F`

be an array of spherical harmonic expansion coefficients with columns arranged by increasing order in absolute value, alternating between negative and positive orders. Then `sph2fourier`

converts the representation into a bivariate Fourier series, and `fourier2sph`

converts it back. Once in a bivariate Fourier series on the sphere, `plan_sph_synthesis`

converts the coefficients to function samples on an equiangular grid that does not sample the poles, and `plan_sph_analysis`

converts them back.

julia> F = sphrandn(Float64, 1024, 2047); # convenience method
julia> P = plan_sph2fourier(F);
julia> PS = plan_sph_synthesis(F);
julia> PA = plan_sph_analysis(F);
julia> @time G = PS*(P*F);
0.090767 seconds (12 allocations: 31.985 MiB, 1.46% gc time)
julia> @time H = P\(PA*G);
0.092726 seconds (12 allocations: 31.985 MiB, 1.63% gc time)
julia> norm(F-H)/norm(F)
2.1541073345177038e-15

Due to the structure of the spherical harmonic connection problem, these transforms may also be performed in-place with `lmul!`

and `ldiv!`

.

See also FastSphericalHarmonics.jl for a simpler interface to the spherical harmonic transforms defined in this package.

## Nonuniform fast Fourier transforms

The NUFFTs are implemented thanks to Alex Townsend:

`nufft1`

assumes uniform samples and noninteger frequencies;`nufft2`

assumes nonuniform samples and integer frequencies;`nufft3 ( = nufft)`

assumes nonuniform samples and noninteger frequencies;`inufft1`

inverts an`nufft1`

; and,`inufft2`

inverts an`nufft2`

.

Here is an example:

julia> n = 10^4;
julia> c = complex(rand(n));
julia> ω = collect(0:n-1) + rand(n);
julia> nufft1(c, ω, eps());
julia> p1 = plan_nufft1(ω, eps());
julia> @time p1*c;
0.002383 seconds (6 allocations: 156.484 KiB)
julia> x = (collect(0:n-1) + 3rand(n))/n;
julia> nufft2(c, x, eps());
julia> p2 = plan_nufft2(x, eps());
julia> @time p2*c;
0.001478 seconds (6 allocations: 156.484 KiB)
julia> nufft3(c, x, ω, eps());
julia> p3 = plan_nufft3(x, ω, eps());
julia> @time p3*c;
0.058999 seconds (6 allocations: 156.484 KiB)

## The Padua Transform

The Padua transform and its inverse are implemented thanks to Michael Clarke. These are optimized methods designed for computing the bivariate Chebyshev coefficients by interpolating a bivariate function at the Padua points on `[-1,1]^2`

.

julia> n = 200;
julia> N = div((n+1)*(n+2), 2);
julia> v = rand(N); # The length of v is the number of Padua points
julia> @time norm(ipaduatransform(paduatransform(v)) - v)/norm(v)
0.007373 seconds (543 allocations: 1.733 MiB)
3.925164683252905e-16

# References

[1] D. Ruiz—Antolín and A. Townsend, A nonuniform fast Fourier transform based on low rank approximation, *SIAM J. Sci. Comput.*, **40**:A529–A547, 2018.

[2] T. S. Gutleb, S. Olver and R. M. Slevinsky, Polynomial and rational measure modifications of orthogonal polynomials via infinite-dimensional banded matrix factorizations, arXiv:2302.08448, 2023.

[3] S. Olver, R. M. Slevinsky, and A. Townsend, Fast algorithms using orthogonal polynomials, *Acta Numerica*, **29**:573—699, 2020.

[4] R. M. Slevinsky, Fast and backward stable transforms between spherical harmonic expansions and bivariate Fourier series, *Appl. Comput. Harmon. Anal.*, **47**:585—606, 2019.

[5] R. M. Slevinsky, Conquering the pre-computation in two-dimensional harmonic polynomial transforms, arXiv:1711.07866, 2017.