single-fft.c
1 //! @file single-fft.c
2 //! @author J. Marcel van der Veer
3 //!
4 //! @section Copyright
5 //!
6 //! This file is part of Algol68G - an Algol 68 compiler-interpreter.
7 //! Copyright 2001-2023 J. Marcel van der Veer [algol68g@xs4all.nl].
8 //!
9 //! @section License
10 //!
11 //! This program is free software; you can redistribute it and/or modify it
12 //! under the terms of the GNU General Public License as published by the
13 //! Free Software Foundation; either version 3 of the License, or
14 //! (at your option) any later version.
15 //!
16 //! This program is distributed in the hope that it will be useful, but
17 //! WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
18 //! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
19 //! more details. You should have received a copy of the GNU General Public
20 //! License along with this program. If not, see [http://www.gnu.org/licenses/].
21
22 //! @section Synopsis
23 //!
24 //! REAL, COMPLEX fast fourier transform.
25
26 #include "a68g.h"
27 #include "a68g-genie.h"
28 #include "a68g-prelude.h"
29
30 #if defined (HAVE_GSL)
31
32 //! @brief Map GSL error handler onto a68g error handler.
33
34 void fft_error_handler (const char *reason, const char *file, int line, int gsl_errno)
35 {
36 if (line != 0) {
37 ASSERT (snprintf (A68 (edit_line), SNPRINTF_SIZE, "%s in line %d of file %s", reason, line, file) >= 0);
38 } else {
39 ASSERT (snprintf (A68 (edit_line), SNPRINTF_SIZE, "%s", reason) >= 0);
40 }
41 diagnostic (A68_RUNTIME_ERROR, A68 (f_entry), ERROR_FFT, A68 (edit_line), gsl_strerror (gsl_errno));
42 exit_genie (A68 (f_entry), A68_RUNTIME_ERROR);
43 }
44
45 //! @brief Detect math errors.
46
47 void fft_test_error (int rc)
48 {
49 if (rc != 0) {
50 fft_error_handler ("math error", "", 0, rc);
51 }
52 }
53
54 //! @brief Pop [] REAL on the stack as complex REAL_T [].
55
56 REAL_T *pop_array_real (NODE_T * p, int *len)
57 {
58 A68_REF desc;
59 A68_ARRAY *arr;
60 A68_TUPLE *tup;
61 int inc, iindex, k;
62 BYTE_T *base;
63 REAL_T *v;
64 A68 (f_entry) = p;
65 // Pop arguments.
66 POP_REF (p, &desc);
67 CHECK_REF (p, desc, M_ROW_REAL);
68 GET_DESCRIPTOR (arr, tup, &desc);
69 *len = ROW_SIZE (tup);
70 if ((*len) <= 0) {
71 return NO_REAL;
72 }
73 v = (REAL_T *) get_heap_space (2 * (size_t) (*len) * sizeof (REAL_T));
74 fft_test_error (v == NO_REAL ? GSL_ENOMEM : GSL_SUCCESS);
75 base = DEREF (BYTE_T, &ARRAY (arr));
76 iindex = VECTOR_OFFSET (arr, tup);
77 inc = SPAN (tup) * ELEM_SIZE (arr);
78 for (k = 0; k < (*len); k++, iindex += inc) {
79 A68_REAL *x = (A68_REAL *) (base + iindex);
80 CHECK_INIT (p, INITIALISED (x), M_REAL);
81 v[2 * k] = VALUE (x);
82 v[2 * k + 1] = 0.0;
83 }
84 return v;
85 }
86
87 //! @brief Push REAL_T [] on the stack as [] REAL.
88
89 void push_array_real (NODE_T * p, REAL_T * v, int len)
90 {
91 A68_REF desc, row;
92 A68_ARRAY arr;
93 A68_TUPLE tup;
94 int inc, iindex, k;
95 BYTE_T *base;
96 A68 (f_entry) = p;
97 NEW_ROW_1D (desc, row, arr, tup, M_ROW_REAL, M_REAL, len);
98 base = DEREF (BYTE_T, &ARRAY (&arr));
99 iindex = VECTOR_OFFSET (&arr, &tup);
100 inc = SPAN (&tup) * ELEM_SIZE (&arr);
101 for (k = 0; k < len; k++, iindex += inc) {
102 A68_REAL *x = (A68_REAL *) (base + iindex);
103 STATUS (x) = INIT_MASK;
104 VALUE (x) = v[2 * k];
105 CHECK_REAL (p, VALUE (x));
106 }
107 PUSH_REF (p, desc);
108 }
109
110 //! @brief Pop [] COMPLEX on the stack as REAL_T [].
111
112 REAL_T *pop_array_complex (NODE_T * p, int *len)
113 {
114 A68_REF desc;
115 A68_ARRAY *arr;
116 A68_TUPLE *tup;
117 int inc, iindex, k;
118 BYTE_T *base;
119 REAL_T *v;
120 A68 (f_entry) = p;
121 // Pop arguments.
122 POP_REF (p, &desc);
123 CHECK_REF (p, desc, M_ROW_COMPLEX);
124 GET_DESCRIPTOR (arr, tup, &desc);
125 *len = ROW_SIZE (tup);
126 if ((*len) <= 0) {
127 return NO_REAL;
128 }
129 v = (REAL_T *) get_heap_space (2 * (size_t) (*len) * sizeof (REAL_T));
130 fft_test_error (v == NO_REAL ? GSL_ENOMEM : GSL_SUCCESS);
131 base = DEREF (BYTE_T, &ARRAY (arr));
132 iindex = VECTOR_OFFSET (arr, tup);
133 inc = SPAN (tup) * ELEM_SIZE (arr);
134 for (k = 0; k < (*len); k++, iindex += inc) {
135 A68_REAL *re = (A68_REAL *) (base + iindex);
136 A68_REAL *im = (A68_REAL *) (base + iindex + SIZE (M_REAL));
137 CHECK_INIT (p, INITIALISED (re), M_COMPLEX);
138 CHECK_INIT (p, INITIALISED (im), M_COMPLEX);
139 v[2 * k] = VALUE (re);
140 v[2 * k + 1] = VALUE (im);
141 }
142 return v;
143 }
144
145 //! @brief Push REAL_T [] on the stack as [] COMPLEX.
146
147 void push_array_complex (NODE_T * p, REAL_T * v, int len)
148 {
149 A68_REF desc, row;
150 A68_ARRAY arr;
151 A68_TUPLE tup;
152 int inc, iindex, k;
153 BYTE_T *base;
154 A68 (f_entry) = p;
155 NEW_ROW_1D (desc, row, arr, tup, M_ROW_COMPLEX, M_COMPLEX, len);
156 base = DEREF (BYTE_T, &ARRAY (&arr));
157 iindex = VECTOR_OFFSET (&arr, &tup);
158 inc = SPAN (&tup) * ELEM_SIZE (&arr);
159 for (k = 0; k < len; k++, iindex += inc) {
160 A68_REAL *re = (A68_REAL *) (base + iindex);
161 A68_REAL *im = (A68_REAL *) (base + iindex + SIZE (M_REAL));
162 STATUS (re) = INIT_MASK;
163 VALUE (re) = v[2 * k];
164 STATUS (im) = INIT_MASK;
165 VALUE (im) = v[2 * k + 1];
166 CHECK_COMPLEX (p, VALUE (re), VALUE (im));
167 }
168 PUSH_REF (p, desc);
169 }
170
171 //! @brief Push prime factorisation on the stack as [] INT.
172
173 void genie_prime_factors (NODE_T * p)
174 {
175 gsl_error_handler_t *save_handler = gsl_set_error_handler (fft_error_handler);
176 A68_INT n;
177 A68_REF desc, row;
178 A68_ARRAY arr;
179 A68_TUPLE tup;
180 int len, inc, iindex, k;
181 BYTE_T *base;
182 gsl_fft_complex_wavetable *wt;
183 A68 (f_entry) = p;
184 POP_OBJECT (p, &n, A68_INT);
185 CHECK_INIT (p, INITIALISED (&n), M_INT);
186 wt = gsl_fft_complex_wavetable_alloc ((size_t) (VALUE (&n)));
187 len = (int) (NF (wt));
188 NEW_ROW_1D (desc, row, arr, tup, M_ROW_INT, M_INT, len);
189 base = DEREF (BYTE_T, &ARRAY (&arr));
190 iindex = VECTOR_OFFSET (&arr, &tup);
191 inc = SPAN (&tup) * ELEM_SIZE (&arr);
192 for (k = 0; k < len; k++, iindex += inc) {
193 A68_INT *x = (A68_INT *) (base + iindex);
194 STATUS (x) = INIT_MASK;
195 VALUE (x) = (int) ((FACTOR (wt))[k]);
196 }
197 gsl_fft_complex_wavetable_free (wt);
198 PUSH_REF (p, desc);
199 (void) gsl_set_error_handler (save_handler);
200 }
201
202 //! @brief PROC ([] COMPLEX) [] COMPLEX fft complex forward
203
204 void genie_fft_complex_forward (NODE_T * p)
205 {
206 gsl_error_handler_t *save_handler = gsl_set_error_handler (fft_error_handler);
207 int len, rc;
208 REAL_T *data;
209 gsl_fft_complex_wavetable *wt;
210 gsl_fft_complex_workspace *ws;
211 A68 (f_entry) = p;
212 data = pop_array_complex (p, &len);
213 fft_test_error (len == 0 ? GSL_EDOM : GSL_SUCCESS);
214 wt = gsl_fft_complex_wavetable_alloc ((size_t) len);
215 ws = gsl_fft_complex_workspace_alloc ((size_t) len);
216 rc = gsl_fft_complex_forward (data, 1, (size_t) len, wt, ws);
217 fft_test_error (rc);
218 push_array_complex (p, data, len);
219 gsl_fft_complex_wavetable_free (wt);
220 gsl_fft_complex_workspace_free (ws);
221 if (data != NO_REAL) {
222 a68_free (data);
223 }
224 (void) gsl_set_error_handler (save_handler);
225 }
226
227 //! @brief PROC ([] COMPLEX) [] COMPLEX fft complex backward
228
229 void genie_fft_complex_backward (NODE_T * p)
230 {
231 gsl_error_handler_t *save_handler = gsl_set_error_handler (fft_error_handler);
232 int len, rc;
233 REAL_T *data;
234 gsl_fft_complex_wavetable *wt;
235 gsl_fft_complex_workspace *ws;
236 A68 (f_entry) = p;
237 data = pop_array_complex (p, &len);
238 fft_test_error (len == 0 ? GSL_EDOM : GSL_SUCCESS);
239 wt = gsl_fft_complex_wavetable_alloc ((size_t) len);
240 ws = gsl_fft_complex_workspace_alloc ((size_t) len);
241 rc = gsl_fft_complex_backward (data, 1, (size_t) len, wt, ws);
242 fft_test_error (rc);
243 push_array_complex (p, data, len);
244 gsl_fft_complex_wavetable_free (wt);
245 gsl_fft_complex_workspace_free (ws);
246 if (data != NO_REAL) {
247 a68_free (data);
248 }
249 (void) gsl_set_error_handler (save_handler);
250 }
251
252 //! @brief PROC ([] COMPLEX) [] COMPLEX fft complex inverse
253
254 void genie_fft_complex_inverse (NODE_T * p)
255 {
256 gsl_error_handler_t *save_handler = gsl_set_error_handler (fft_error_handler);
257 int len, rc;
258 REAL_T *data;
259 gsl_fft_complex_wavetable *wt;
260 gsl_fft_complex_workspace *ws;
261 A68 (f_entry) = p;
262 data = pop_array_complex (p, &len);
263 fft_test_error (len == 0 ? GSL_EDOM : GSL_SUCCESS);
264 wt = gsl_fft_complex_wavetable_alloc ((size_t) len);
265 ws = gsl_fft_complex_workspace_alloc ((size_t) len);
266 rc = gsl_fft_complex_inverse (data, 1, (size_t) len, wt, ws);
267 fft_test_error (rc);
268 push_array_complex (p, data, len);
269 gsl_fft_complex_wavetable_free (wt);
270 gsl_fft_complex_workspace_free (ws);
271 if (data != NO_REAL) {
272 a68_free (data);
273 }
274 (void) gsl_set_error_handler (save_handler);
275 }
276
277 //! @brief PROC ([] REAL) [] COMPLEX fft forward
278
279 void genie_fft_forward (NODE_T * p)
280 {
281 gsl_error_handler_t *save_handler = gsl_set_error_handler (fft_error_handler);
282 int len, rc;
283 REAL_T *data;
284 gsl_fft_complex_wavetable *wt;
285 gsl_fft_complex_workspace *ws;
286 A68 (f_entry) = p;
287 data = pop_array_real (p, &len);
288 fft_test_error (len == 0 ? GSL_EDOM : GSL_SUCCESS);
289 wt = gsl_fft_complex_wavetable_alloc ((size_t) len);
290 ws = gsl_fft_complex_workspace_alloc ((size_t) len);
291 rc = gsl_fft_complex_forward (data, 1, (size_t) len, wt, ws);
292 fft_test_error (rc);
293 push_array_complex (p, data, len);
294 gsl_fft_complex_wavetable_free (wt);
295 gsl_fft_complex_workspace_free (ws);
296 if (data != NO_REAL) {
297 a68_free (data);
298 }
299 (void) gsl_set_error_handler (save_handler);
300 }
301
302 //! @brief PROC ([] COMPLEX) [] REAL fft backward
303
304 void genie_fft_backward (NODE_T * p)
305 {
306 gsl_error_handler_t *save_handler = gsl_set_error_handler (fft_error_handler);
307 int len, rc;
308 REAL_T *data;
309 gsl_fft_complex_wavetable *wt;
310 gsl_fft_complex_workspace *ws;
311 A68 (f_entry) = p;
312 data = pop_array_complex (p, &len);
313 fft_test_error (len == 0 ? GSL_EDOM : GSL_SUCCESS);
314 wt = gsl_fft_complex_wavetable_alloc ((size_t) len);
315 ws = gsl_fft_complex_workspace_alloc ((size_t) len);
316 rc = gsl_fft_complex_backward (data, 1, (size_t) len, wt, ws);
317 fft_test_error (rc);
318 push_array_real (p, data, len);
319 gsl_fft_complex_wavetable_free (wt);
320 gsl_fft_complex_workspace_free (ws);
321 if (data != NO_REAL) {
322 a68_free (data);
323 }
324 (void) gsl_set_error_handler (save_handler);
325 }
326
327 //! @brief PROC ([] COMPLEX) [] REAL fft inverse
328
329 void genie_fft_inverse (NODE_T * p)
330 {
331 gsl_error_handler_t *save_handler = gsl_set_error_handler (fft_error_handler);
332 int len, rc;
333 REAL_T *data;
334 gsl_fft_complex_wavetable *wt;
335 gsl_fft_complex_workspace *ws;
336 A68 (f_entry) = p;
337 data = pop_array_complex (p, &len);
338 fft_test_error (len == 0 ? GSL_EDOM : GSL_SUCCESS);
339 wt = gsl_fft_complex_wavetable_alloc ((size_t) len);
340 ws = gsl_fft_complex_workspace_alloc ((size_t) len);
341 rc = gsl_fft_complex_inverse (data, 1, (size_t) len, wt, ws);
342 fft_test_error (rc);
343 push_array_real (p, data, len);
344 gsl_fft_complex_wavetable_free (wt);
345 gsl_fft_complex_workspace_free (ws);
346 if (data != NO_REAL) {
347 a68_free (data);
348 }
349 (void) gsl_set_error_handler (save_handler);
350 }
351
352 #endif
