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 GSL 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 static 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 static void fft_test_error (int ret)
48 {
49 if (ret != 0) {
50 fft_error_handler ("math error", "", 0, ret);
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 // Pop arguments.
59 A68_REF desc;
60 POP_REF (p, &desc);
61 CHECK_REF (p, desc, M_ROW_REAL);
62 A68_ARRAY *arr; A68_TUPLE *tup;
63 GET_DESCRIPTOR (arr, tup, &desc);
64 *len = ROW_SIZE (tup);
65 if ((*len) <= 0) {
66 return NO_REAL;
67 }
68 REAL_T *v = (REAL_T *) get_heap_space (2 * (size_t) (*len) * sizeof (REAL_T));
69 fft_test_error (v == NO_REAL ? GSL_ENOMEM : GSL_SUCCESS);
70 BYTE_T *base = DEREF (BYTE_T, &ARRAY (arr));
71 int index = VECTOR_OFFSET (arr, tup);
72 int inc = SPAN (tup) * ELEM_SIZE (arr);
73 for (int k = 0; k < (*len); k++, index += inc) {
74 A68_REAL *x = (A68_REAL *) (base + index);
75 CHECK_INIT (p, INITIALISED (x), M_REAL);
76 v[2 * k] = VALUE (x);
77 v[2 * k + 1] = 0.0;
78 }
79 return v;
80 }
81
82 //! @brief Push REAL_T [] on the stack as [] REAL.
83
84 void push_array_real (NODE_T * p, REAL_T * v, int len)
85 {
86 A68_REF desc, row; A68_ARRAY arr; A68_TUPLE tup;
87 NEW_ROW_1D (desc, row, arr, tup, M_ROW_REAL, M_REAL, len);
88 BYTE_T *base = DEREF (BYTE_T, &ARRAY (&arr));
89 int index = VECTOR_OFFSET (&arr, &tup);
90 int inc = SPAN (&tup) * ELEM_SIZE (&arr);
91 for (int k = 0; k < len; k++, index += inc) {
92 A68_REAL *x = (A68_REAL *) (base + index);
93 STATUS (x) = INIT_MASK;
94 VALUE (x) = v[2 * k];
95 CHECK_REAL (p, VALUE (x));
96 }
97 PUSH_REF (p, desc);
98 }
99
100 //! @brief Pop [] COMPLEX on the stack as REAL_T [].
101
102 REAL_T *pop_array_complex (NODE_T * p, int *len)
103 {
104 // Pop arguments.
105 A68_REF desc;
106 POP_REF (p, &desc);
107 CHECK_REF (p, desc, M_ROW_COMPLEX);
108 A68_ARRAY *arr; A68_TUPLE *tup;
109 GET_DESCRIPTOR (arr, tup, &desc);
110 *len = ROW_SIZE (tup);
111 if ((*len) <= 0) {
112 return NO_REAL;
113 }
114 REAL_T *v = (REAL_T *) get_heap_space (2 * (size_t) (*len) * sizeof (REAL_T));
115 fft_test_error (v == NO_REAL ? GSL_ENOMEM : GSL_SUCCESS);
116 BYTE_T *base = DEREF (BYTE_T, &ARRAY (arr));
117 int index = VECTOR_OFFSET (arr, tup);
118 int inc = SPAN (tup) * ELEM_SIZE (arr);
119 for (int k = 0; k < (*len); k++, index += inc) {
120 A68_REAL *re = (A68_REAL *) (base + index);
121 A68_REAL *im = (A68_REAL *) (base + index + SIZE (M_REAL));
122 CHECK_INIT (p, INITIALISED (re), M_COMPLEX);
123 CHECK_INIT (p, INITIALISED (im), M_COMPLEX);
124 v[2 * k] = VALUE (re);
125 v[2 * k + 1] = VALUE (im);
126 }
127 return v;
128 }
129
130 //! @brief Push REAL_T [] on the stack as [] COMPLEX.
131
132 void push_array_complex (NODE_T * p, REAL_T * v, int len)
133 {
134 A68_REF desc, row; A68_ARRAY arr; A68_TUPLE tup;
135 NEW_ROW_1D (desc, row, arr, tup, M_ROW_COMPLEX, M_COMPLEX, len);
136 BYTE_T *base = DEREF (BYTE_T, &ARRAY (&arr));
137 int index = VECTOR_OFFSET (&arr, &tup);
138 int inc = SPAN (&tup) * ELEM_SIZE (&arr);
139 for (int k = 0; k < len; k++, index += inc) {
140 A68_REAL *re = (A68_REAL *) (base + index);
141 A68_REAL *im = (A68_REAL *) (base + index + SIZE (M_REAL));
142 STATUS (re) = INIT_MASK;
143 VALUE (re) = v[2 * k];
144 STATUS (im) = INIT_MASK;
145 VALUE (im) = v[2 * k + 1];
146 CHECK_COMPLEX (p, VALUE (re), VALUE (im));
147 }
148 PUSH_REF (p, desc);
149 }
150
151 //! @brief Push prime factorisation on the stack as [] INT.
152
153 void genie_prime_factors (NODE_T * p)
154 {
155 gsl_error_handler_t *save_handler = gsl_set_error_handler (fft_error_handler);
156 A68_INT n;
157 POP_OBJECT (p, &n, A68_INT);
158 CHECK_INIT (p, INITIALISED (&n), M_INT);
159 gsl_fft_complex_wavetable *wt = gsl_fft_complex_wavetable_alloc ((size_t) (VALUE (&n)));
160 int len = (int) (NF (wt));
161 A68_REF desc, row; A68_ARRAY arr; A68_TUPLE tup;
162 NEW_ROW_1D (desc, row, arr, tup, M_ROW_INT, M_INT, len);
163 BYTE_T *base = DEREF (BYTE_T, &ARRAY (&arr));
164 int index = VECTOR_OFFSET (&arr, &tup);
165 int inc = SPAN (&tup) * ELEM_SIZE (&arr);
166 for (int k = 0; k < len; k++, index += inc) {
167 A68_INT *x = (A68_INT *) (base + index);
168 STATUS (x) = INIT_MASK;
169 VALUE (x) = (int) ((FACTOR (wt))[k]);
170 }
171 gsl_fft_complex_wavetable_free (wt);
172 PUSH_REF (p, desc);
173 (void) gsl_set_error_handler (save_handler);
174 }
175
176 //! @brief PROC ([] COMPLEX) [] COMPLEX fft complex forward
177
178 void genie_fft_complex_forward (NODE_T * p)
179 {
180 gsl_error_handler_t *save_handler = gsl_set_error_handler (fft_error_handler);
181 int len;
182 REAL_T *data = pop_array_complex (p, &len);
183 fft_test_error (len == 0 ? GSL_EDOM : GSL_SUCCESS);
184 gsl_fft_complex_wavetable *wt = gsl_fft_complex_wavetable_alloc ((size_t) len);
185 gsl_fft_complex_workspace *ws = gsl_fft_complex_workspace_alloc ((size_t) len);
186 int ret = gsl_fft_complex_forward (data, 1, (size_t) len, wt, ws);
187 fft_test_error (ret);
188 push_array_complex (p, data, len);
189 gsl_fft_complex_wavetable_free (wt);
190 gsl_fft_complex_workspace_free (ws);
191 a68_free (data);
192 (void) gsl_set_error_handler (save_handler);
193 }
194
195 //! @brief PROC ([] COMPLEX) [] COMPLEX fft complex backward
196
197 void genie_fft_complex_backward (NODE_T * p)
198 {
199 gsl_error_handler_t *save_handler = gsl_set_error_handler (fft_error_handler);
200 int len;
201 REAL_T *data = pop_array_complex (p, &len);
202 fft_test_error (len == 0 ? GSL_EDOM : GSL_SUCCESS);
203 gsl_fft_complex_wavetable *wt = gsl_fft_complex_wavetable_alloc ((size_t) len);
204 gsl_fft_complex_workspace *ws = gsl_fft_complex_workspace_alloc ((size_t) len);
205 int ret = gsl_fft_complex_backward (data, 1, (size_t) len, wt, ws);
206 fft_test_error (ret);
207 push_array_complex (p, data, len);
208 gsl_fft_complex_wavetable_free (wt);
209 gsl_fft_complex_workspace_free (ws);
210 a68_free (data);
211 (void) gsl_set_error_handler (save_handler);
212 }
213
214 //! @brief PROC ([] COMPLEX) [] COMPLEX fft complex inverse
215
216 void genie_fft_complex_inverse (NODE_T * p)
217 {
218 gsl_error_handler_t *save_handler = gsl_set_error_handler (fft_error_handler);
219 int len;
220 REAL_T *data = pop_array_complex (p, &len);
221 fft_test_error (len == 0 ? GSL_EDOM : GSL_SUCCESS);
222 gsl_fft_complex_wavetable *wt = gsl_fft_complex_wavetable_alloc ((size_t) len);
223 gsl_fft_complex_workspace *ws = gsl_fft_complex_workspace_alloc ((size_t) len);
224 int ret = gsl_fft_complex_inverse (data, 1, (size_t) len, wt, ws);
225 fft_test_error (ret);
226 push_array_complex (p, data, len);
227 gsl_fft_complex_wavetable_free (wt);
228 gsl_fft_complex_workspace_free (ws);
229 a68_free (data);
230 (void) gsl_set_error_handler (save_handler);
231 }
232
233 //! @brief PROC ([] REAL) [] COMPLEX fft forward
234
235 void genie_fft_forward (NODE_T * p)
236 {
237 gsl_error_handler_t *save_handler = gsl_set_error_handler (fft_error_handler);
238 int len;
239 REAL_T *data = pop_array_real (p, &len);
240 fft_test_error (len == 0 ? GSL_EDOM : GSL_SUCCESS);
241 gsl_fft_complex_wavetable *wt = gsl_fft_complex_wavetable_alloc ((size_t) len);
242 gsl_fft_complex_workspace *ws = gsl_fft_complex_workspace_alloc ((size_t) len);
243 int ret = gsl_fft_complex_forward (data, 1, (size_t) len, wt, ws);
244 fft_test_error (ret);
245 push_array_complex (p, data, len);
246 gsl_fft_complex_wavetable_free (wt);
247 gsl_fft_complex_workspace_free (ws);
248 a68_free (data);
249 (void) gsl_set_error_handler (save_handler);
250 }
251
252 //! @brief PROC ([] COMPLEX) [] REAL fft backward
253
254 void genie_fft_backward (NODE_T * p)
255 {
256 gsl_error_handler_t *save_handler = gsl_set_error_handler (fft_error_handler);
257 int len;
258 REAL_T *data = pop_array_complex (p, &len);
259 fft_test_error (len == 0 ? GSL_EDOM : GSL_SUCCESS);
260 gsl_fft_complex_wavetable *wt = gsl_fft_complex_wavetable_alloc ((size_t) len);
261 gsl_fft_complex_workspace *ws = gsl_fft_complex_workspace_alloc ((size_t) len);
262 int ret = gsl_fft_complex_backward (data, 1, (size_t) len, wt, ws);
263 fft_test_error (ret);
264 push_array_real (p, data, len);
265 gsl_fft_complex_wavetable_free (wt);
266 gsl_fft_complex_workspace_free (ws);
267 a68_free (data);
268 (void) gsl_set_error_handler (save_handler);
269 }
270
271 //! @brief PROC ([] COMPLEX) [] REAL fft inverse
272
273 void genie_fft_inverse (NODE_T * p)
274 {
275 gsl_error_handler_t *save_handler = gsl_set_error_handler (fft_error_handler);
276 int len;
277 REAL_T *data = pop_array_complex (p, &len);
278 fft_test_error (len == 0 ? GSL_EDOM : GSL_SUCCESS);
279 gsl_fft_complex_wavetable *wt = gsl_fft_complex_wavetable_alloc ((size_t) len);
280 gsl_fft_complex_workspace *ws = gsl_fft_complex_workspace_alloc ((size_t) len);
281 int ret = gsl_fft_complex_inverse (data, 1, (size_t) len, wt, ws);
282 fft_test_error (ret);
283 push_array_real (p, data, len);
284 gsl_fft_complex_wavetable_free (wt);
285 gsl_fft_complex_workspace_free (ws);
286 a68_free (data);
287 (void) gsl_set_error_handler (save_handler);
288 }
289
290 #endif
