-
Notifications
You must be signed in to change notification settings - Fork 0
/
api.tex
440 lines (399 loc) · 19.3 KB
/
api.tex
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
\section{Detailed API}
Boost.Math.FFT brings two kinds of interfaces: one-time and planned transforms.
The first one is inspired by the standard library's \texttt{std::transform}
and the second by \fftw's plan interface. Listings
\ref{ls:complex_dft_transform} and \ref{ls:complex_dft_plan} are two examples of
the use of the one-time and planned interfaces respectively.
The planned transform is used with
the declaration of a plan variable and subsequently calling its methods.
The internal state of the plan is determined by the size of the \dft\ transform,
and throughout its life it caches the necessary data and memory allocation
needed to perform \dft s on request. In principle it should be faster to call a
method of a plan, than the equivalent routine from the one-time interface.
The library's interface can be further classified according to the domain where
the \dft\ is applied: let it be complex numbers, real and halfcomplex, or a
purely algebraic set satisfying the ring axioms. This distinction is necessary
because all three behave differently and the user expects different interfaces.
Consider the complex case. There are only two types of \dft\ we can perform on
a $\Complex^n$ vector, they are the \verb|forward| and the \verb|backward|,
because there are only two roots of unity in that domain:
$w = \exp( - 2 \pi i/n)$
and
$w^{-1} = \exp( 2 \pi i/n)$.
While in general for an algebraic transform on a type that satisfies the ring
axioms the \dft\ engine has no way of knowing \emph{a priori} what is the root
of unity $w$ and that one must be supplied by the user either to the constructor
of the plan or as an extra argument to the one-time transform routine. In this
case, once $w$ is known there are again only two types of \dft s: \verb|forward|
when we use $w$ or \verb|backward| when using $w^{-1}$.
The real/halfcomplex interface works very differently from the complex and the
algebraic. Since the result of any \dft\ on a $\Real^n$ vector yields always a
vector in $\hComplex^n$ and viceversa; this interface contains the methods:
\verb|real_to_halfcomplex| and \verb|halfcomplex_to_real| to emphasize that
property. Because of the trivial relations between the results of the forward
and backward transforms applied to real and/or halfcomplex data (see fig.
\ref{fig:halfcomplex}) we do not provide a choice for forward or backward
transform. \verb|real_to_halfcomplex| will in effect compute a forward \dft\ on
a $\Real^n$ vector and \verb|halfcomplex_to_real| will compute a backward \dft\ on a
$\hComplex^n$ vector. We also provide the methods \verb|real_to_complex|
and \verb|complex_to_real| similar to the previous \verb|real_to_halfcomplex|
and \verb|halfcomplex_to_real| respectively, where the halfcomplex part is
represented in an array of $n$ complex numbers, hence keeping redundant values.
The library computes the requested operations using one of three backends:
\gsl, \fftw\ or our own Boost licensed engine that we call \bsl.
\subsection{Complex \dft}
Complex valued \dft\ are best handled by Boost.Math.FFT complex oriented
interface. All three backends are accessible through the plan-like classes:
\verb|bsl_dft|, \verb|gsl_dft| and \verb|fftw_dft|, representing the three
engines \bsl, \gsl\ and \fftw\ respectively. They share the same interface hence
it is enough to describe only the member functions of \verb|bsl_dft|.
\begin{lstlisting}[language=C++,caption=Complex plan declaration.]
template< class Complex = std::complex<double>,
class Allocator = std::allocator<Complex> >
class bsl_dft;
\end{lstlisting}
Notice that even though the complex \dft\ plans are template on the complex
type, not all backends handle all complex types. \gsl\ only works for
\verb|std::complex<double>| and \fftw\ work with a richer collection of
precision types
\verb|std::complex<float>|, \verb|std::complex<double>|,
\verb|std::complex<long double>| and
\verb|boost::multiprecision::complex128|\footnote{or equivalently
\texttt{boost::multiprecision::complex<boost::multiprecision::float128>}.}.
On the other hand, the implementation of \verb|bsl_dft| allows, in principle, to
compute \dft s
for any complex type as long as it has the essential complex
algebraic methods plus sine and cosine functions for its internal real
\verb|value_type|.
\verb|bsl_dft| can handle, for instance besides the previously mentioned complex types,
also the family of \verb|boost::multiprecision|'s \verb|cpp_complex|
and \verb|mpc_complex| based on \verb|cpp_bin_float| and
\verb|mpfr_float_backend| real types respectively; eg.
\verb|cpp_complex_50|,
\verb|cpp_complex_quad| and
\verb|mpc_complex_50|.
\subsubsection*{Member functions}
\begin{lstlisting}[language=C++,caption=Complex plan constructor.,label=ls:complex_plan]
template< class Complex, class Allocator >
constexpr bsl_dft<Complex,Allocator>::bsl_dft(
unsigned int n, const Allocator& alloc = Allocator{} );
\end{lstlisting}
The constructor of a complex \dft\ plan takes two arguments (listing
\ref{ls:complex_plan}): the first is an integer number $N$ that specifies the size of
the transform, the second is a
reference to an allocator (allocators are discussed in section
\ref{subsec:allocators}).
\begin{lstlisting}[language=C++,caption=Complex plan forward function.]
template<class InputIt,
class OutputIt>
void forward(InputIt first1, InputIt last1, OutputIt d_first);
\end{lstlisting}
\begin{lstlisting}[language=C++,caption=Complex plan backward function.]
template<class InputIt,
class OutputIt>
void backward(InputIt first1, InputIt last1, OutputIt d_first);
\end{lstlisting}
\begin{lstlisting}[language=C++,caption=Complex plan size.]
constexpr std::size_t size() const;
\end{lstlisting}
\begin{lstlisting}[language=C++,caption=Resize.]
void resize(std::size_t new_size);
\end{lstlisting}
\subsubsection*{One-time transform methods}
The same functionality is accessible through the one-time transform methods
below.
\begin{lstlisting}[language=C++,caption=Forward complex transform.]
template<class InputIt,
class OutputIt>
void bsl_transform::forward(InputIt first1, InputIt last1, OutputIt d_first);
\end{lstlisting}
\begin{lstlisting}[language=C++,caption=Backward complex transform.]
template<class InputIt,
class OutputIt>
void bsl_transform::backward(InputIt first1, InputIt last1, OutputIt d_first);
\end{lstlisting}
\subsubsection*{Examples}
See listings \ref{ls:complex_dft_transform} and \ref{ls:complex_dft_plan}
for examples of the use of the complex one-time transform and plan interfaces.
\subsection{Algebraic \dft}
Algebraic type \dft\ is understood in the sense of definition \ref{def:dft}.
It can be applied to algebraic types that satisfy the ring axioms (see appendix
\ref{app:ring}) as long as a root of unity is known for the required \dft\ size.
Boost.Math.FFT has achieved this level of abtraction with the use of C++
templates. There is not equivalent functionality in \gsl\ or \fftw\ libraries
for this kind of transform, hence we only describe the \bsl\ plan-like backend
called \verb|bsl_algebraic_dft|.
\begin{lstlisting}[language=C++,caption=Algebraic plan.,label=ls:alg_plan]
template< class T = std::complex<double>,
class Allocator = std::allocator<T> >
class bsl_algebraic_dft;
\end{lstlisting}
The constructor of an algebraic \dft\ plan takes three arguments (listing
\ref{ls:alg_plan}): the first is an integer number $N$ that specifies the size of
the transform, the second is an $N$-root of unity $w$ and the third is a
reference to an allocator (allocators are discussed in section
\ref{subsec:allocators}).
\subsubsection*{Member functions}
\begin{lstlisting}[language=C++,caption=Algebraic plan constructor.]
template< class T, class Allocator >
constexpr bsl_algebraic_dft<T,Allocator>::bsl_algebraic_dft(
unsigned int n, T w, const Allocator& alloc = Allocator{} );
\end{lstlisting}
\begin{lstlisting}[language=C++,caption=Algebraic plan forward function.]
template<class InputIt,
class OutputIt>
void forward(InputIt first1, InputIt last1, OutputIt d_first);
\end{lstlisting}
\begin{lstlisting}[language=C++,caption=Algebraic plan backward function.]
template<class InputIt,
class OutputIt>
void backward(InputIt first1, InputIt last1, OutputIt d_first);
\end{lstlisting}
\begin{lstlisting}[language=C++,caption=Size.]
constexpr std::size_t size() const;
\end{lstlisting}
% The fact is that the algebraic_dft does have a resize method, but resizing the
% DFT here must be followed with a change of the internal root of unity w. A
% correct resize method here requires two arguments: the new size and the new
% value of w.
% \begin{lstlisting}[language=C++]
% void resize(std::size_t new_size);
% \end{lstlisting}
\subsubsection*{One-time transform methods}
The same functionality is accessible through the one-time transform methods
below.
\begin{lstlisting}[language=C++,caption=Forward algebraic transform.]
template<class InputIt,
class OutputIt,
class T>
void bsl_algebraic_transform::forward(
InputIt first1, InputIt last1, OutputIt d_first, const T w);
\end{lstlisting}
\begin{lstlisting}[language=C++,caption=Backward algebraic transform.]
template<class InputIt,
class OutputIt,
class T>
void bsl_algebraic_transform::backward(
InputIt first1, InputIt last1, OutputIt d_first, const T w);
\end{lstlisting}
\subsubsection*{Examples}
Complex numbers are a particular case of an algebraic ring, which fortunately we
know for every $N$ only two possible $N$-roots of unity; ie.
$\exp(-i \frac{2\pi}{N})$
and
$\exp(i \frac{2\pi}{N})$, where $i=\sqrt{-1}$. Hence we could use the algebraic
interface to compute complex \dft s. Listing \ref{ls:complex_alg_dft} shows how
could we compute the \dft\ of a complex array using the algebraic interface.
\begin{lstlisting}[language=C++,caption=Complex DFT using the generic interface.,label=ls:complex_alg_dft]
using Real = double;
using Complex = std::complex<Real>;
std::vector<Complex> A{1.,-1.,3.,.5,9.},B,C;
const int N = A.size();
const Real w_phase = 2*boost::math::constants::pi<Real>()/N;
const Complex w{cos(w_phase), -sin(w_phase)}; // w is an N-root of unity in the complex
// compute the complex forward DFT using the algebraic interface
bsl_algebraic_transform::forward(A.cbegin(),A.cend(),std::back_inserter(B), w);
// compute the complex forward DFT using the complex interface
bsl_transform::forward(A.cbegin(),A.cend(),std::back_inserter(C));
// B and C are equal up to floating point errors
\end{lstlisting}
\begin{lstlisting}[language=C++,caption=Integer fast multiplication.]
using M_int = fft::my_modulo_lib::mint<Z337>;
// 85 is a primitive root of 337,
// ie. the smallest number k such that 85^k mod 337 = 1 is k=336
const M_int w{85};
const M_int inv_8{M_int{8}.inverse()};
// Multiplying 1234 times 5678 = 7006652
std::vector<M_int> A{4, 3, 2, 1, 0, 0, 0, 0};
std::vector<M_int> B{8, 7, 6, 5, 0, 0, 0, 0};
// forward FFT
boost::math::fft::bsl_algebraic_dft<M_int> P(A.size(),w);
P.forward(A.cbegin(),A.cend(),A.begin());
P.forward(B.cbegin(),B.cend(),B.begin());
// convolution in Fourier space
std::vector<M_int> AB;
std::transform(A.begin(), A.end(), B.begin(),
std::back_inserter(AB),
[](M_int x, M_int y) { return x * y; });
// backwards FFT
P.backward(AB.cbegin(),AB.cend(),AB.begin());
std::transform(AB.begin(), AB.end(), AB.begin(),
[&inv_8](M_int x) { return x * inv_8; });
// carry the remainders in base 10
std::vector<int> C;
M_int r{0};
for (auto x : AB)
{
auto y = x + r;
C.emplace_back(int(y) % 10);
r = M_int(int(y) / 10);
}
// C[] = 7006652
\end{lstlisting}
\subsection{Real \dft}
Real valued \dft\ are best handled by Boost.Math.FFT real oriented
interface. The plan-like classes specialized for real numbers are denoted
\verb|bsl_rdft|, \verb|gsl_rdft| and \verb|fftw_rdft|,
representing the three
engines \bsl, \gsl\ and \fftw\ respectively. They share the same interface hence
it is enough to describe only the member functions of \verb|bsl_rdft|.
The real methods take advantage of the symmetries involved in the real valued
\dft s described in Subsection \ref{subsec:real}. In principle one could compute
those real transforms using complex numbers and the complex interface (eg. by
setting to zero the imaginary part of the input and ignoring the halfcomplex
symmetry), but the real routines will perform
half of the operations and occupy half the space in memory than the complex
unconstrained counterpart.
\begin{lstlisting}[language=C++,caption=Real plan.]
template< class Real = double,
class Allocator = std::allocator<Real> >
class bsl_rdft;
\end{lstlisting}
Real plans, as well as the complex plans seen previously, are templates on the
type that represents the real numbers. Again \gsl\ backend will work only for
\verb|double| representations while \fftw\ supports \verb|float|, \verb|double|,
and \verb|long double|, and \verb|boost::multiprecision::float128|.
The \bsl\ functions are all templates on the type,
hence it will accept any type representing real numbers.
\subsubsection*{Member functions}
\begin{lstlisting}[language=C++,caption=Real plan constructor.,label=ls:real_plan]
template< class Real, class Allocator >
constexpr bsl_rdft<Real,Allocator>::bsl_rdft(
unsigned int n, const Allocator& alloc = Allocator{} );
\end{lstlisting}
The constructor of a real \dft\ plan takes two arguments (listing
\ref{ls:real_plan}): the first is an integer number $N$ that specifies the size of
the transform, the second is a
reference to an allocator (allocators are discussed in section
\ref{subsec:allocators}).
\begin{lstlisting}[language=C++,caption=Real to halfcomplex.]
template<class InputIt,
class OutputIt>
void real_to_halfcomplex(InputIt first1, InputIt last1, OutputIt d_first);
\end{lstlisting}
\begin{lstlisting}[language=C++,caption=Real to complex.]
template<class InputIt,
class OutputIt>
void real_to_complex(InputIt first1, InputIt last1, OutputIt d_first);
\end{lstlisting}
\begin{lstlisting}[language=C++,caption=Halfcomplex to real.]
template<class InputIt,
class OutputIt>
void halfcomplex_to_real(InputIt first1, InputIt last1, OutputIt d_first);
\end{lstlisting}
\begin{lstlisting}[language=C++,caption=Complex to real.]
template<class InputIt,
class OutputIt>
void complex_to_real(InputIt first1, InputIt last1, OutputIt d_first);
\end{lstlisting}
\begin{lstlisting}[language=C++,caption=Size.]
constexpr std::size_t size() const;
\end{lstlisting}
\begin{lstlisting}[language=C++,caption=Resize.]
void resize(std::size_t new_size);
\end{lstlisting}
\subsubsection*{One-time transform methods}
The same functionality is accessible through the one-time transform methods
below.
\begin{lstlisting}[language=C++,caption=One-time transform real to halfcomplex.]
template<class InputIt,
class OutputIt>
void bsl_real_transform::real_to_halfcomplex(
InputIt first1, InputIt last1, OutputIt d_first);
\end{lstlisting}
\begin{lstlisting}[language=C++,caption=One-time transform real to complex.]
template<class InputIt,
class OutputIt>
void bsl_real_transform::real_to_complex(
InputIt first1, InputIt last1, OutputIt d_first);
\end{lstlisting}
\begin{lstlisting}[language=C++,caption=One-time transform halfcomplex to real.]
template<class InputIt,
class OutputIt>
void bsl_real_transform::halfcomplex_to_real(
InputIt first1, InputIt last1, OutputIt d_first);
\end{lstlisting}
\begin{lstlisting}[language=C++,caption=One-time transform complex to real.]
template<class InputIt,
class OutputIt>
void bsl_real_transform::complex_to_real(
InputIt first1, InputIt last1, OutputIt d_first);
\end{lstlisting}
\subsubsection*{Examples}
\begin{lstlisting}[language=C++,caption=Polynomial multiplication.]
using Real = double;
std::vector<Real> A{1.,4.,-5.,1.,0.,0.,0.,0.};
std::vector<Real> B{-1.,1.,2.,3.,0.,0.,0.,0.};
const std::size_t N = A.size();
std::vector< std::complex<Real> > TA(N),TB(N);
boost::math::fft::bsl_rdft<Real> P(N);
P.real_to_complex(A.begin(),A.end(),TA.begin());
P.real_to_complex(B.begin(),B.end(),TB.begin());
for(unsigned int i=0;i<N;++i)
TA[i]*=TB[i];
std::vector<Real> C(N);
P.complex_to_real(TA.begin(),TA.end(),C.begin());
std::transform(C.begin(), C.end(), C.begin(),
[N](Real x) { return x / N; });
// C[] = [-1, -3, 11, 5, 3, -13, 3, 0,]
\end{lstlisting}
\subsection{Convolution}
Currently, there is no convolution transform implemented as a member of any of
the plan-like classes. Instead the circular convolution (see definition
\ref{def:convolution}) of two series of complex
numbers is handled by the one-time transform described in listing~\ref{ls:convolution}.
\begin{lstlisting}[language=C++,caption=Convolution.,label=ls:convolution]
template<typename InputIt1,
typename InputIt2,
typename OutputIt>
void bsl_transform::convolution(
InputIt1 first1, InputIt1 last1, InputIt2 first2, OutputIt d_first);
\end{lstlisting}
\subsection{The one-time transform generalities}
For each plan type and member function there is an equivalent one-time transform
method. These are static methods of a template class denoted \verb|transform|.
This design choice allows to select one backends programmatically.
For instance \verb|bsl_transform| is an alias for \verb|transform<bsl_dft<>>|,
\verb|bsl_real_transform| is and alias for \verb|transform<bsl_rdft<>>|, etc.
Table \ref{tab:transform} summarizes the naming of the one-time transform family.
\begin{table}[h]
\centering
\begin{tabular}{ll}
Alias & Instantiation \\
\hline
\verb|bsl_transform| & \verb|transform< bsl_dft<> >|\\
\verb|bsl_algebraic_transform| & \verb|transform< bsl_algebraic_dft<> >|\\
\verb|bsl_real_transform| & \verb|transform< bsl_rdft<> >|\\
\verb|gsl_transform| & \verb|transform< gsl_dft<> >|\\
\verb|gsl_real_transform| & \verb|transform< gsl_rdft<> >|\\
\verb|fftw_transform| & \verb|transform< fftw_dft<> >|\\
\verb|fftw_real_transform| & \verb|transform< fftw_rdft<> >|\\
\end{tabular}
\caption{Summary of the one-time transform classes.}
\label{tab:transform}
\end{table}
\subsection{The allocator template parameter}\label{subsec:allocators}
The plan-like interface has a template parameter named \verb|Allocator|.
The plan uses the \emph{allocator} supplied through the constructor, to allocate
all the internal heap memory necessary to compute \dft s.
Notice that very often, \dft s cannot be performed in place, and scratch space
is used. Furthermore, the work space is required to be a contiguous array in
memory and there is no way to deduce that property for any iterator before
C++20's \emph{contiguous iterator} concept; hence scratch space is used even in
cases when the in-place \dft\ could be possible in principle.
The example in listing \ref{ls:pmr_allocator}, shows how could the allocator
parameter could be used to harness memory allocation.
\begin{lstlisting}[language=C++,caption=Example use of custom allocators.,label=ls:pmr_allocator]
const std::size_t N = 100;
std::array<char,30000> buf;
boost::container::pmr::monotonic_buffer_resource
pool{buf.data(),buf.size(),boost::container::pmr::null_memory_resource()};
using Real = double;
using Complex = boost::multiprecision::complex<Real>;
using Allocator = boost::container::pmr::polymorphic_allocator<Complex>;
std::vector<Complex,Allocator> A(N,Complex(),&pool);
std::vector<Complex,Allocator> B(N,Complex(),&pool);
// ...
bsl_dft<Complex,Allocator> plan(N,&pool);
plan.forward(A.begin(),A.end(),B.begin());
\end{lstlisting}