From d6e71a4a1f351fbb9f06a633573cd698180015f8 Mon Sep 17 00:00:00 2001 From: Jonathan Schilling Date: Fri, 1 May 2026 00:11:29 +0200 Subject: [PATCH] review using Claude Code and update for consistency and correctness --- .gitignore | 4 +- CLAUDE.md | 54 ++++ README.md | 618 ++++++++++++++++------------------------- src/app_flux_surface.c | 25 +- src/test_1d_redft01.c | 80 +++--- src/test_1d_redft10.c | 75 +++-- src/test_1d_redft11.c | 81 +++--- src/test_1d_rodft01.c | 87 +++--- src/test_1d_rodft10.c | 79 +++--- src/test_1d_rodft11.c | 77 ++--- src/util.h | 7 + 11 files changed, 564 insertions(+), 623 deletions(-) create mode 100644 CLAUDE.md diff --git a/.gitignore b/.gitignore index fa5fd25..4f4323e 100644 --- a/.gitignore +++ b/.gitignore @@ -39,4 +39,6 @@ test_util test_util.dat # Jupyter Notebook checkpoints -*.ipynb_checkpoints \ No newline at end of file +*.ipynb_checkpoints + +arXiv-0708.4399v2 \ No newline at end of file diff --git a/CLAUDE.md b/CLAUDE.md new file mode 100644 index 0000000..e9f965a --- /dev/null +++ b/CLAUDE.md @@ -0,0 +1,54 @@ +# CLAUDE.md + +This file provides guidance to Claude Code (claude.ai/code) when working with code in this repository. + +## Project purpose + +A teaching collection of FFTW 3.3.9 examples in C. Each example computes a 1D or 2D DFT both via FFTW and via a hand-rolled naive reference, then compares the two — the goal is to make explicit the symmetries, scaling factors, and array layouts that FFTW assumes. The README is the long-form tutorial; the `src/` programs are the runnable companions to it. + +## Build / run / test + +Top-level Makefile builds with `gcc` and links `-lfftw3 -lm`. Each target is a single-file executable that `#include`s `src/util.h`. + +``` +make # build all test_* and app_* binaries (compiled into repo root, not src/) +make test # build all targets, then run each one via the %_test pattern rule +make clean # remove binaries and generated *.dat files +make # build a single target, e.g. `make test_1d_redft00` or `make app_flux_surface` +./ # run a single test/app (binaries land in the repo root) +``` + +There is no separate test framework — each `test_*.c` has its own `main` that returns non-zero on mismatch, and `make test` simply runs them in sequence. + +The `eqn/` subdirectory has its own Makefile that renders `.tex` snippets to `.png` via `./pnglatex` for the README; ignore it unless touching the rendered equations. + +## Architecture + +**Header-only utility layer (`src/util.h`).** All shared helpers — random/zero array fills, element-wise comparison with tolerance, file dumpers in NumPy-loadable formats, a generalized-DFT reference (`dft_1d_cplx`), and the LCFS file readers — live as `static`-style definitions directly in the header. Each `.c` file is compiled into its own executable, so there is no link-time conflict. Adding a new helper means editing `util.h` directly; there is no `util.c`. + +**Per-example pattern.** Every `test_*.c` follows the same skeleton: +1. Allocate input/output arrays with `fftw_alloc_real` / `fftw_alloc_complex`. +2. Create the FFTW plan **before** filling input (FFTW_ESTIMATE; some plans destroy input on execute). +3. Fill input with `fill_random_*`. +4. Compute a reference output by a direct, naive transform formula. +5. `fftw_execute`, then `compare_*` with `eps = 1e-12`. +6. Free plan and arrays. +For r2r tests (REDFT/RODFT), the reference is built by constructing the "logically equivalent" full DFT input (mirrored/anti-mirrored according to the transform's implicit symmetry) and comparing FFTW's output to a slice of that full DFT. + +**Application layer (`src/app_*.c`).** Three real-world stellarator-physics examples that consume the same primitives: `app_magn_axis` (1D c2r), `app_magn_axis_stellsym` (1D REDFT01, exploits stellarator symmetry), `app_flux_surface` (2D c2r, reads `lcfs.dat`). They write `axis_R*.dat` / `lcfs_R.dat` / `lcfs_Z.dat` to the repo root for plotting; `clean` removes these. + +**Python companions.** `src/plot_lcfs_realspace.py` (Mayavi) reads the `app_flux_surface` output. `src/redft_symmetries.py` and `src/redft_symmetries.ipynb` are exploratory NumPy notebooks on even/odd-harmonic DCT decomposition. + +## Non-obvious conventions + +- **2D arrays are row-major flat `double*` / `fftw_complex*`**, indexed as `arr[i*cols + j]`. The dumpers and comparators all assume this. +- **VMEC vs FFTW sign convention** (relevant in `app_flux_surface.c` and the corresponding README section): VMEC uses angle argument `(m θ − n ζ)` while FFTW's 2D DFT uses `(m θ + n ζ)`. Coefficients from VMEC must be placed at index `-n` in the FFTW input — see the `idx_in` computation. `m=0` row is special-cased because VMEC keeps the positive-`n` sign there. +- **Hermitian-symmetry scaling for c2r/r2c.** For real-input DFTs the complex array has `n/2+1` entries; all but DC (and Nyquist when `n` even) need a factor-of-2 weighting in the manual reference. Both array-layout figure (`img/array_structures.png`) and `test_1d_c2r.c` document this. +- **R\*DFT00 caveat.** Per the FFTW manual, `REDFT00`/`RODFT00` are numerically less stable; `app_magn_axis_stellsym` deliberately uses `REDFT01` instead. Don't "fix" a stellsym example to use 00. +- **Output `.dat` files are gitignored** but consumed by Python plotters and by some tests' own assertions — don't assume `make clean` is safe to run mid-debugging session if a Python plot is open. + +## When extending + +- New transform example → add `src/test_.c`, add the bare name to `TARGETS` in the Makefile. The pattern rule `$(TARGETS): % : src/%.c src/util.h` handles compilation; no other wiring needed. +- New utility → add it to `src/util.h` as a definition (not just a declaration). Keep it dependency-free beyond `` and libc. +- The repo also carries Eclipse CDT (`.cproject`, `.project`) and PyDev (`.pydevproject`) project files; leave them alone unless asked. diff --git a/README.md b/README.md index b063b66..258f0ee 100644 --- a/README.md +++ b/README.md @@ -69,7 +69,9 @@ The size of the DFT is specified via the variable `n` and the direction (forward Complex-valued arrays `in`, `ref_out` and `fftw_out` are allocated to hold the input (*X_k*) and outputs (*Y_k*) of the DFT. A plan for the corresponding DFT is created using `fftw_plan_dft_1d`. Only after this, the input data is written into the `in` array. -The reference output is computed now (before calling `fftw_execute`), since in some later examples (most notably multi-dimensional `c2r` transforms), +The reference output is computed now (before calling `fftw_execute`), since in some later examples +(`c2r` transforms in particular: 1D `c2r` destroys its input by default but supports `FFTW_PRESERVE_INPUT`, +while multi-dimensional `c2r` cannot generally be requested to preserve its input at all) FFTW overwrites the input data and for good practise, we keep this in mind already now. Note that this is an out-of-place transform, since `in` and `fftw_out` are allocated to be separate arrays. Next the FFTW transform can be executed via `fftw_execute`. @@ -243,20 +245,49 @@ the sign of a `r2c` DFT is always `FFTW_FORWARD`. ### 1D real-to-real Certain symmetries can be assumed for a given real input array of which a DFT is to be computed -that lead to the output array being purely real as well. This is another gain of a factor of 2 in speed -and memory usage over `r2c`/`c2r` transforms. +that lead to the output array being purely real as well. +Code that knows its data has these symmetries can use a real-to-real plan and avoid +the complex-array bookkeeping of an `r2c`/`c2r` transform; it does not, however, generally make +the transform itself faster, since FFTW implements the `r2r` transforms internally as a real DFT +of the *logical* size — see the +[FFTW manual](http://fftw.org/fftw3_doc/Real-even_002fodd-DFTs-_0028cosine_002fsine-transforms_0029.html#DOCF4), +which notes in particular that `REDFT00`/`RODFT00` are slower than the corresponding `r2c` because +of numerical-stability mitigations. Depending on even (odd) parity of the input array, the transform outputs have even (odd) parity. They are therefore called Discrete Cosine Transform (DCT) and Discrete Sine Transfrom (DST), respectively. -The logical size of the corresponding DFT is denoted as *N*. -The actual array size given to FFTW is denoted by `n`. -For the DCT/DST types implemented in FFTW, *N* is always even. -Note that this does not pose any restrictions on the input array sizes `n`. -One can think of the logical DFT input array as one that FFTW 'sees' internally and computes a regular DFT of. -The resulting output array is purely real and features the named symmetry properties, -since the logical input array was 'constructed' from the given input array to have the desired symmetry properties. - -In below example plots used to illustrate these symmetry properties/assumptions, +**Logical-equivalent DFT view.** Each `r2r` transform is mathematically equivalent to a regular (unshifted) +complex DFT of an enlarged real-valued sequence *X̃* of size *Ñ*, where *Ñ* depends on the transform type: + +| type | logical DFT size *Ñ* | one-sentence description of *X̃* | +|:-------:|:--------------------:|:---------------------------------| +| REDFT00 | *Ñ* = 2(`n`−1) | input data with even-around-0 and even-around-(n−1) extension | +| REDFT10 | *Ñ* = 4`n` | input data placed at odd positions, even-around-0 and even-around-2n | +| REDFT01 | *Ñ* = 4`n` | input data with doubled boundary at *X̃[0]=½X₀*, *X̃[2n]=−½X₀*, even-around-0 / odd-around-n | +| REDFT11 | *Ñ* = 8`n` | input data placed at odd positions, even-around-0 / odd-around-2n / even-around-4n / odd-around-6n | +| RODFT00 | *Ñ* = 2(`n`+1) | input data with explicit zeros at *j=0* and *j=n+1*, odd-around-0 and odd-around-(n+1) | +| RODFT10 | *Ñ* = 4`n` | input data placed at odd positions, odd-around-0 and odd-around-2n | +| RODFT01 | *Ñ* = 4`n` | input data with explicit zero at *X̃[n]=0* and *X̃[3n]=0*, odd-around-0 / even-around-n | +| RODFT11 | *Ñ* = 8`n` | input data placed at odd positions, odd-around-0 / even-around-2n / odd-around-4n / even-around-6n | + +The actual array size given to FFTW is denoted by `n`. For the unshifted types (REDFT00, RODFT00) the logical size *Ñ* +is the full size of an unshifted DFT whose input has the named symmetry. For the half-sample-shifted types +(REDFT10/01, RODFT10/01) the *Ñ* = 4`n` picture absorbs the half-sample shift into a *zero-interleaved* extension +(real samples sit at odd positions of *X̃*, all even positions are zero). For the doubly-shifted type-IV transforms +(REDFT11, RODFT11) the same idea is applied twice and *Ñ* = 8`n`. This is the construction used +by Shao & Johnson (arXiv [0708.4399](https://arxiv.org/abs/0708.4399), 2008, equations 5–7 for DCT-IV +and equations 22–24 for DCT-III); the type-II / type-III transforms follow the same template at half the logical size, +and DST variants flip the corresponding symmetries (odd ↔ even). + +The output values *Y_k* of the FFTW `r2r` transform are then recovered as a particular subsequence of the +size-*Ñ* DFT of *X̃* (the first *n* entries for the type-00 transforms, the *(2k+1)*-th entries for the +half-sample-shifted types, with an overall prefactor of 2 due to FFTW's normalization convention — see the per-type sections below). + +Each of the example programs `src/test_1d_*.c` builds the corresponding *X̃* explicitly, computes the size-*Ñ* +DFT via the helper `dft_1d_cplx_unshifted` (defined in `src/util.h`), extracts the relevant subset of outputs, +and checks that the result agrees with FFTW's `r2r` output to within `1e-12`. + +In the example plots used to illustrate these symmetry properties, random Fourier coefficients have been sampled and transformed back to real-space using the named inverse transforms. This allows to 'evaluate' the input data also in the samples given in the (small) input arrays. In these plots, red dashed vertical lines indicate even symmetry (*f(x)=f(-x)*) about the indicated position @@ -274,20 +305,7 @@ The shifting becomes necessary when formulating the symmetry properties over sam vs. symmetry axis that are possibly located at half-integer locations. The parity and symmetry properties of the output array are those of the input array for the inverse transform. -For all transforms, a periodicity of *N* is assumed for the *logical* input array as *X_j = X_{N+j}* where *X* is the input data array. - -Here is a quick overview table to indicate the assumed symmetries in the input array for the following types of `r2r` DFTs: - -| type | actual `r2r` input | logically-equivalent DFT input | -| :------:| :---------: | :-------------------: | -| REDFT00 | `a b c d e` | `a b c d e d c b` | -| REDFT10 | `a b c d ` | `a b c d d c b a` | -| REDFT01 | `a b c d ` | `a b c d 0 -d -c -b` | -| REDFT11 | `a b c d ` | `a b c d -d -c -b -a` | -| RODFT00 | `a b c ` | `0 a b c 0 -c -b -a` | -| RODFT10 | `a b c d ` | `a b c d -d -c -b -a` | -| RODFT01 | `a b c d ` | `a b c d c b a 0` | -| RODFT11 | `a b c d ` | `a b c d d c b a` | +For all transforms, a periodicity of *Ñ* is assumed for the *logical* input array as *X̃_j = X̃_{Ñ+j}*. For the DCTs, please also consider https://en.wikipedia.org/wiki/Discrete_cosine_transform#Formal_definition and in particular https://upload.wikimedia.org/wikipedia/commons/a/ae/DCT-symmetries.svg . @@ -382,7 +400,9 @@ The full example can be found in [`src/test_1d_redft00.c`](src/test_1d_redft00.c #### REDFT10 (DCT-II) In case of the real-valued even-parity DFT with shifted input data (REDFT10), -also called the DCT-II, the corresponding logical DFT size is given by *N* = 2`n`, corresponding to `n` = *N*/2. +also called the DCT-II, the appropriate logical-equivalent (unshifted) DFT size is *Ñ* = 4`n`, +following the same family of constructions as Shao & Johnson use for DCT-III but with the +half-sample shift on the input (rather than the output) side. This function is commonly known as "the" DCT. The formal definition of the REDFT10 is given below: @@ -394,161 +414,94 @@ The input array is assumed to have even symmetry around *j=-0.5* and even symmet ![REDFT10](img/redft10.png) -In above figure, the lowercase letters *a* to *e* refer to the input data *abcd* for the size-4 REDFT10, -which is logically equivalent to a size-8 DFT with real-valued input data *abcddcba*. - -In order to demonstrate the use of this method, -the logically equivalent DFT input is filled appropriately and its output is checked against `REDFT10`. -In the following code, `in` is the input array (size `n`) given to `REDFT10` -and `in_logical` is the (complex-valued) input array (size *N*) handed to a -[generic 1D DFT](https://en.wikipedia.org/wiki/Discrete_Fourier_transform#Generalized_DFT_(shifted_and_non-linear_phase)). -Similarly, `out` is the output array (size `n`) from `REDFT10` -and `out_logical` is the output array (size *N*) from a generic 1D DFT. - -Here is how the symmetric input is generated: +We construct the size-4`n` logical input *X̃* by placing each *X_j* at odd index *2j+1* and at its +mirror *4n−2j−1* (so the four-fold symmetry pattern is even-around-0 and even-around-2`n`, +with all even-index entries zero): ```C -// the first half of the array is identical -for (int i = 0; i < n; ++i) { - in_logical[i] = in[i]; -} - -// second half is filled according to even symmetry around n-0.5 -for (int i = 0; i < n; ++i) { - in_logical[n + i] = in[n - 1 - i]; +int Nlog = 4 * n; +fill_zero_1d_cplx(Nlog, in_logical); +for (int j = 0; j < n; ++j) { + in_logical[2*j + 1] = in[j]; + in_logical[Nlog - 2*j - 1] = in[j]; } ``` -The checks are a little bit more involved. -The logically equivalent DFT output should be purely real-valued: - -```C -for (int i = 0; i < N; ++i) { - if (fabs(cimag(out_logical[i])) > eps) { - printf("error: imag of [%d] is %g\n", i, cimag(out_logical[i])); - status = 1; - } else { - printf("imag of [%d] is %g\n", i, cimag(out_logical[i])); - } -} -``` - -The first `n` values should be identical between `REDFT10` and the generalized DFT: +After running `dft_1d_cplx_unshifted(Nlog, in_logical, out_logical)`, every DFT output is real +(by even-around-0) and the FFTW REDFT10 output is recovered as *Y_k = Re(X̃_k)* for *k = 0, …, n−1*. +The FFTW prefactor of 2 is built into the construction: each *X_j* is written into *X̃* twice +(at *2j+1* and at *4n−2j−1*), and the two copies' cosine contributions add to give the FFTW formula. ```C -for (int i = 0; i < n; ++i) { - delta = creal(out_logical[i]) - out[i]; +double delta; +for (int k = 0; k < n; ++k) { + delta = creal(out_logical[k]) - out[k]; if (fabs(delta) > eps) { - printf("error: delta of [%d] is %g\n", i, delta); + printf("error: delta of [%d] is %g\n", k, delta); status = 1; } else { - printf("match of [%d] (delta=%g)\n", i, delta); + printf("match of [%d] (delta=%g)\n", k, delta); } } ``` -Odd symmetry around `n` implies that the value at `n` is zero: - -```C -if (fabs(creal(out_logical[n])) > eps) { - printf("error: delta of [%d] is %g\n", n, creal(out_logical[n])); - status = 1; -} else { - printf("match of [%d] (delta=%g)\n", n, creal(out_logical[n])); -} -``` - -The remaining values should have odd symmetry around `n`: - -```C -for (int i = 1; i < n; ++i) { - delta = creal(out_logical[n + i]) - (-out[n - i]); - if (fabs(delta) > eps) { - printf("error: delta of [%d] is %g\n", n + i, delta); - status = 1; - } else { - printf("match of [%d] (delta=%g)\n", n + i, delta); - } -} -``` +The four-fold symmetry of *X̃* also forces *X̃_n = 0* (the DCT-II "Nyquist" boundary, which the +test verifies separately). The full example can be found in [`src/test_1d_redft10.c`](src/test_1d_redft10.c). #### REDFT01 (DCT-III) In case of the real-valued even-parity DFT with shifted output data (REDFT01), -also called the DCT-III, the corresponding logical DFT size is given by *N* = 2`n`, corresponding to `n` = *N*/2. +also called the DCT-III, the appropriate logical-equivalent (unshifted) DFT size is *Ñ* = 4`n`, following +the construction of [Shao & Johnson 2008 (arXiv 0708.4399), eqs. 22–24](https://arxiv.org/abs/0708.4399). This function is commonly known as "the" inverse DCT (IDCT). The formal definition of the REDFT01 is given below: ![REDFT01 formula](eqn/redft01.png) +This is *FFTW's* form: *Y_k = X_0 + 2 Σ_{j≥1} X_j cos[π j (k+½)/n]*. The Shao-Johnson DCT-III +(paper eq. 22) is *C_k^III = Σ_{j=0}^{n-1} X_j cos[π j (k+½)/n]*, summing uniformly from *j=0*. +The two forms are related by *Y_k = 2 · C_k^III − X_0*. + The inverse of this transform is REDFT10. The input array is assumed to have even symmetry around *j=0* and odd symmetry around *j=n*. ![REDFT01](img/redft01.png) -In order to demonstrate the use of this method, -the logically equivalent DFT input is filled appropriately and its output is checked against `REDFT01`. -In the following code, `in` is the input array (size `n`) given to `REDFT01` -and `in_logical` is the (complex-valued) input array (size *N*) handed to a -[generic 1D DFT](https://en.wikipedia.org/wiki/Discrete_Fourier_transform#Generalized_DFT_(shifted_and_non-linear_phase)). -Similarly, `out` is the output array (size `n`) from `REDFT01` -and `out_logical` is the output array (size *N*) from a generic 1D DFT. - -Here is how the symmetric input is generated: - -```C -// the first half of the array is identical -for (int i = 0; i < n; ++i) { - in_logical[i] = in[i]; -} - -// second half is filled according to odd symmetry around n -for (int i = 1; i < n; ++i) { - in_logical[n + i] = -in[n - i]; -} -``` - -The checks are a little bit more involved. -The logically equivalent DFT output should be purely real-valued: +We construct the size-4`n` logical input *X̃* per Shao-Johnson eqs. 23–24 (paper line 705-715): +the four symmetry axes are even at *j=0* and *j=2n*, odd at *j=n* and *j=3n*. The data fills *X̃* in full +(no zero interleaving), with two boundary terms *X̃_0 = X_0/2* and *X̃_{2n} = −X_0/2*: ```C -for (int i = 0; i < N; ++i) { - if (fabs(cimag(out_logical[i])) > eps) { - printf("error: imag of [%d] is %g\n", i, cimag(out_logical[i])); - status = 1; - } else { - printf("imag of [%d] is %g\n", i, cimag(out_logical[i])); - } +int Nlog = 4 * n; +fill_zero_1d_cplx(Nlog, in_logical); +in_logical[0] = 0.5 * in[0]; +in_logical[2*n] = -0.5 * in[0]; +// in_logical[n] = 0; in_logical[3*n] = 0; (already zero) +for (int j = 1; j < n; ++j) { + double v = 0.25 * in[j]; + in_logical[j] = v; // x̃_j = X_j/4 + in_logical[Nlog - j] = v; // x̃_{4n-j} = X_j/4 (even-around-0) + in_logical[2*n - j] = -v; // x̃_{2n-j} = -X_j/4 + in_logical[2*n + j] = -v; // x̃_{2n+j} = -X_j/4 (even-around-2n with sign flip) } ``` -The first `n` values should be identical between `REDFT01` and the generalized DFT: +After running `dft_1d_cplx_unshifted(Nlog, in_logical, out_logical)`, the (2k+1)-th DFT outputs +recover the unnormalized Shao-Johnson DCT-III: *X̃_{2k+1} = C_k^III*. FFTW's REDFT01 is then +*Y_k = 2 · X̃_{2k+1} − X_0*: ```C -for (int i = 0; i < n; ++i) { - delta = creal(out_logical[i]) - out[i]; - if (fabs(delta) > eps) { - printf("error: delta of [%d] is %g\n", i, delta); - status = 1; - } else { - printf("match of [%d] (delta=%g)\n", i, delta); - } -} -``` - -The remaining values should have even symmetry around `n-0.5`: - -```C -for (int i = 0; i < n; ++i) { - delta = creal(out_logical[n + i]) - out[n - 1 - i]; +double delta; +for (int k = 0; k < n; ++k) { + delta = (2.0 * creal(out_logical[2*k + 1]) - in[0]) - out[k]; if (fabs(delta) > eps) { - printf("error: delta of [%d] is %g\n", n + i, delta); + printf("error: delta of [%d] is %g\n", k, delta); status = 1; } else { - printf("match of [%d] (delta=%g)\n", n + i, delta); + printf("match of [%d] (delta=%g)\n", k, delta); } } ``` @@ -558,77 +511,79 @@ The full example can be found in [`src/test_1d_redft01.c`](src/test_1d_redft01.c #### REDFT11 (DCT-IV) In case of the real-valued even-parity DFT with both shifted input and output data (REDFT11), -also called the DCT-IV, the corresponding logical DFT size is given by *N* = 2`n`, corresponding to `n` = *N*/2. +also called the DCT-IV, the appropriate logical-equivalent (unshifted) DFT size is *Ñ* = 8`n`, following +the construction of [Shao & Johnson 2008 (arXiv 0708.4399), eqs. 5–7](https://arxiv.org/abs/0708.4399). +The half-sample shifts in the input and output of REDFT11 are absorbed into a *zero-interleaved* + +quadruple-symmetric extension of the original `n` samples. The formal definition of the REDFT11 is given below: ![REDFT11 formula](eqn/redft11.png) +This is *FFTW's* form: the prefactor is **2**. The Shao-Johnson DCT-IV (paper eq. 4) is +*C_k^IV = Σ X_j cos[π(j+½)(k+½)/n]* — the same kernel without the prefactor; +the relationship is therefore *Y_k = 2 · C_k^IV*. The paper notes (line 417-419): +*"The transform can be made orthogonal by multiplying with √(2/N), but for our purposes the unnormalized form is more convenient"* — +so both prefactor 1 (paper) and prefactor 2 (FFTW) are simply different unnormalized conventions for the same kernel. + The inverse of this transform is REDFT11 itself. The input array is assumed to have even symmetry around *j=-0.5* and odd symmetry around *j=n-0.5*. ![REDFT11](img/redft11.png) -In order to demonstrate the use of this method, -the logically equivalent DFT input is filled appropriately and its output is checked against `REDFT11`. -In the following code, `in` is the input array (size `n`) given to `REDFT11` -and `in_logical` is the (complex-valued) input array (size *N*) handed to a -[generic 1D DFT](https://en.wikipedia.org/wiki/Discrete_Fourier_transform#Generalized_DFT_(shifted_and_non-linear_phase)). -Similarly, `out` is the output array (size `n`) from `REDFT11` -and `out_logical` is the output array (size *N*) from a generic 1D DFT. - -Here is how the symmetric input is generated: +To verify, we explicitly construct the size-8`n` logical input *X̃* and compare its (unshifted) DFT +against FFTW's REDFT11 output. Per Shao-Johnson eqs. 6–7, *X̃* has zero at every even index, and at the +odd indices it carries quadruple copies of *X_j* with the four-fold symmetry pattern +(even-around-0, odd-around-2`n`, even-around-4`n`, odd-around-6`n`): ```C -// the first half of the array is identical -for (int i = 0; i < n; ++i) { - in_logical[i] = in[i]; -} - -// second half is filled according to odd symmetry around n-0.5 -for (int i = 0; i < n; ++i) { - in_logical[n + i] = -in[n - 1 - i]; +int Nlog = 8 * n; +fill_zero_1d_cplx(Nlog, in_logical); +for (int j = 0; j < n; ++j) { + double v = 0.25 * in[j]; + in_logical[2*j + 1] = v; // +X_j/4 at j=2j+1 (even-around-0) + in_logical[Nlog - 2*j - 1] = v; // +X_j/4 at j=8n-2j-1 (mirror) + in_logical[4*n - 2*j - 1] = -v; // -X_j/4 at j=4n-2j-1 (odd-around-2n) + in_logical[4*n + 2*j + 1] = -v; // -X_j/4 at j=4n+2j+1 (mirror) } ``` -The checks are a little bit more involved. -The logically equivalent DFT output should be purely real-valued: +The factor of *¼* in front of every entry is part of the SJ construction; it is absorbed into the +identity *cos(πℓ/N) = ¼(ω⁴ℓ − ω⁴ᴺ⁻⁴ℓ − ω⁴ᴺ⁺⁴ℓ + ω⁸ᴺ⁻⁴ℓ)* used to derive eq. 5 of the paper. + +The unshifted DFT of *X̃* is computed via `dft_1d_cplx_unshifted(Nlog, in_logical, out_logical)` +(see [`src/util.h`](src/util.h)); FFTW's REDFT11 is run on the original `n` samples in parallel. + +The first verification is that the entire DFT of *X̃* is purely real-valued +(a consequence of the even-around-0 symmetry of *X̃*) and that *only the odd-index outputs are non-zero* +(a consequence of the zero-interleaved structure of *X̃*): ```C -for (int i = 0; i < N; ++i) { - if (fabs(cimag(out_logical[i])) > eps) { - printf("error: imag of [%d] is %g\n", i, cimag(out_logical[i])); +for (int k = 0; k < n; ++k) { + if (fabs(cimag(out_logical[2*k + 1])) > eps) { + printf("error: imag of Y_logical[%d] is %g\n", 2*k + 1, cimag(out_logical[2*k + 1])); status = 1; - } else { - printf("imag of [%d] is %g\n", i, cimag(out_logical[i])); } } -``` - -The first `n` values should be identical between `REDFT11` and the generalized DFT: - -```C -for (int i = 0; i < n; ++i) { - delta = creal(out_logical[i]) - out[i]; - if (fabs(delta) > eps) { - printf("error: delta of [%d] is %g\n", i, delta); +for (int k = 0; k < 4*n; ++k) { + if (fabs(creal(out_logical[2*k])) > eps || fabs(cimag(out_logical[2*k])) > eps) { + printf("error: even-index Y_logical[%d] is non-zero\n", 2*k); status = 1; - } else { - printf("match of [%d] (delta=%g)\n", i, delta); } } ``` -The remaining values should have odd symmetry around `n-0.5`: +Finally, *Y_k* (FFTW REDFT11) is recovered as twice the *(2k+1)*-th DFT output of *X̃*: ```C -for (int i = 0; i < n; ++i) { - delta = creal(out_logical[n + i]) - (-out[n - 1 - i]); +double delta; +for (int k = 0; k < n; ++k) { + delta = 2.0 * creal(out_logical[2*k + 1]) - out[k]; if (fabs(delta) > eps) { - printf("error: delta of [%d] is %g\n", n + i, delta); + printf("error: delta of [%d] is %g\n", k, delta); status = 1; } else { - printf("match of [%d] (delta=%g)\n", n + i, delta); + printf("match of [%d] (delta=%g)\n", k, delta); } } ``` @@ -639,10 +594,11 @@ The full example can be found in [`src/test_1d_redft11.c`](src/test_1d_redft11.c In case of the real-valued odd-parity DFT with no shifts in either input or output array (RODFT00), also called the DST-I, the corresponding logical DFT size is given by *N* = 2(`n`+1), corresponding to `n` = *N*/2-1. -Note that the periodicity of *N* of the logical input array in combination with odd symmetry *X_j = -X_{n-j}* -leads to *X_0 = -X_0* which is equivalent to *X_0 = 0*. -This first always-zero element of the input array is not explicitly included in the input to FFTW -and the input array thus has a size of one less and the indices of the symmetry axis shift by 1. +The input array is assumed to have odd symmetry around *j=-1* and odd symmetry also around *j=n* (see below). +Odd symmetry around *j=-1* directly implies *X_{-1} = -X_{-1}*, hence *X_{-1} = 0*; the same argument at *j=n* gives *X_n = 0*. +FFTW therefore omits these always-zero boundary elements and only stores *X_0, ..., X_{n-1}*; in the size-*N* logical-equivalent array +written down explicitly below, these zero entries are placed at indices *0* and *n+1* respectively (i.e. the symmetry-axis indices +of the logical array shift by 1 relative to FFTW's convention). The formal definition of the RODFT00 is given below: @@ -752,7 +708,8 @@ The full example can be found in [`src/test_1d_rodft00.c`](src/test_1d_rodft00.c #### RODFT10 (DST-II) In case of the real-valued odd-parity DFT with shifted input data (RODFT10), -also called the DST-II, the corresponding logical DFT size is given by *N* = 2`n`, corresponding to `n` = *N*/2. +also called the DST-II, the appropriate logical-equivalent (unshifted) DFT size is *Ñ* = 4`n`, +the sine analog of REDFT10's construction. This function is commonly known as "the" DST. The formal definition of the RODFT10 is given below: @@ -764,92 +721,47 @@ The input array is assumed to have odd symmetry around *j=-0.5* and odd symmetry ![RODFT10](img/rodft10.png) -In order to demonstrate the use of this method, -the logically equivalent DFT input is filled appropriately and its output is checked against `RODFT10`. -In the following code, `in` is the input array (size `n`) given to `RODFT10` -and `in_logical` is the (complex-valued) input array (size *N*) handed to a -[generic 1D DFT](https://en.wikipedia.org/wiki/Discrete_Fourier_transform#Generalized_DFT_(shifted_and_non-linear_phase)). -Similarly, `out` is the output array (size `n`) from `RODFT10` -and `out_logical` is the output array (size *N*) from a generic 1D DFT. - -Here is how the symmetric input is generated: +We construct the size-4`n` logical input *X̃* by placing each *X_j* at odd index *2j+1* and *−X_j* +at the mirror position *4n−2j−1* (so the symmetries are odd-around-0 and odd-around-2`n`): ```C -// the first half of the array is identical -for (int i = 0; i < n; ++i) { - in_logical[i] = in[i]; -} - -// second half is filled according to odd symmetry around (n-0.5) -for (int i = 0; i < n; ++i) { - in_logical[n + i] = -in[n - 1 - i]; -} -``` - -The checks are a little bit more involved. -The logically equivalent DFT output should be purely imaginary-valued: - -```C -for (int i = 0; i < N; ++i) { - if (fabs(creal(out_logical[i])) > eps) { - printf("error: real of [%d] is %g\n", i, creal(out_logical[i])); - status = 1; - } else { - printf("real of [%d] is %g\n", i, creal(out_logical[i])); - } -} -``` - -The output of `RODFT10` is shifted by one index to the left with respect to the logically-equivalent DFT, -since the first input is constrained to be zero (resulting in odd parity about `-1` in the input to `RODFT10`). -Odd symmetry about `0` implies that the imaginary part of the first output should be zero: - -```C -if (fabs(cimag(out_logical[0])) > eps) { - printf("error: delta of [%d] is %g\n", 0, -cimag(out_logical[0])); - status = 1; -} else { - printf("match of [%d] (delta=%g)\n", 0, -cimag(out_logical[0])); +int Nlog = 4 * n; +fill_zero_1d_cplx(Nlog, in_logical); +for (int j = 0; j < n; ++j) { + in_logical[2*j + 1] = in[j]; + in_logical[Nlog - 2*j - 1] = -in[j]; } ``` -The next `n` values should have the output of `RODFT10` in their negative imaginary parts -with one index offset to account for the first zero in the input: +After running `dft_1d_cplx_unshifted(Nlog, in_logical, out_logical)`, every DFT output is purely imaginary +(odd-around-0 of *X̃*). One can show analytically that +*X̃_k = −2i Σ_j X_j sin[π(j+½) k / n]*; comparing with the FFTW formula +*Y_k = 2 Σ_j X_j sin[π(j+½)(k+1)/n]* shows that the FFTW RODFT10 output is recovered as +*Y_k = −Im(X̃_{k+1})* for *k = 0, …, n−1* (note the *(k+1)* offset reflecting RODFT10's right-boundary +anti-periodicity that produces the *(k+1)* in the FFTW formula). ```C -for (int i = 0; i < n; ++i) { - delta = -cimag(out_logical[i + 1]) - out[i]; +double delta; +for (int k = 0; k < n; ++k) { + delta = -cimag(out_logical[k + 1]) - out[k]; if (fabs(delta) > eps) { - printf("error: delta of [%d] is %g\n", i + 1, delta); + printf("error: delta of [%d] is %g\n", k, delta); status = 1; } else { - printf("match of [%d] (delta=%g)\n", i + 1, delta); + printf("match of [%d] (delta=%g)\n", k, delta); } } ``` -The remaining values have even symmetry around `n-1` -(note again that the real output from `RODFT10` needs to be compared against -the negative imaginary parts in the output of the logically equivalent DFT shifted by one full index): - -```C -for (int i = 0; i < n - 1; ++i) { - delta = -cimag(out_logical[n + 1 + i]) - out[n - 2 - i]; - if (fabs(delta) > eps) { - printf("error: delta of [%d] is %g\n", n + 1 + i, delta); - status = 1; - } else { - printf("match of [%d] (delta=%g)\n", n + 1 + i, delta); - } -} -``` +The DC bin *X̃_0 = 0* by the odd-around-0 symmetry of *X̃* (verified separately in the test). The full example can be found in [`src/test_1d_rodft10.c`](src/test_1d_rodft10.c). #### RODFT01 (DST-III) In case of the real-valued odd-parity DFT with shifted output data (RODFT01), -also called the DST-III, the corresponding logical DFT size is given by *N* = 2`n`, corresponding to `n` = *N*/2. +also called the DST-III, the appropriate logical-equivalent (unshifted) DFT size is *Ñ* = 4`n`, +the same as for REDFT01. This function is commonly known as "the" inverse DST (IDST). The formal definition of the RODFT01 is given below: @@ -861,74 +773,41 @@ The input array is assumed to have odd symmetry around *j=-1* and even symmetry ![RODFT01](img/rodft01.png) -In order to demonstrate the use of this method, -the logically equivalent DFT input is filled appropriately and its output is checked against `RODFT01`. -In the following code, `in` is the input array (size `n`) given to `RODFT01` -and `in_logical` is the (complex-valued) input array (size *N*) handed to a -[generic 1D DFT](https://en.wikipedia.org/wiki/Discrete_Fourier_transform#Generalized_DFT_(shifted_and_non-linear_phase)). -Similarly, `out` is the output array (size `n`) from `RODFT01` -and `out_logical` is the output array (size *N*) from a generic 1D DFT. - -The input of `RODFT01` is shifted by one index to the left with respect to the logically-equivalent DFT, -since the first input is constrained to be zero (resulting in odd parity about `-1` in the input to `RODFT01`). -Here is how the symmetric input is generated: - -```C -// a zero in the first entry is needed to satisfy -// odd symmetry about -1 in the shifted input array -in_logical[0] = 0.0; - -// the first half of the array is identical up to shift -for (int i = 0; i < n; ++i) { - in_logical[i + 1] = in[i]; -} - -// second half is filled according to even symmetry around (n-1) -for (int i = 0; i < n - 1; ++i) { - in_logical[n + 1 + i] = in[n - 2 - i]; -} -``` +The Shao-Johnson paper §V.C +([arXiv 0708.4399, eq. 32](https://arxiv.org/abs/0708.4399)) +shows that a DST-III is reducible to a DCT-III on the input-reversed array, with a sign-flip on every other output: -The checks are a little bit more involved. -The logically equivalent DFT output should be purely imaginary-valued: +> *S_k^III(X) = (-1)^k · Σ_{m=0}^{N−1} X_{N−m} cos[π m (k+½)/N] = (-1)^k · C_k^III(X_reversed)* -```C -for (int i = 0; i < N; ++i) { - if (fabs(creal(out_logical[i])) > eps) { - printf("error: real of [%d] is %g\n", i, creal(out_logical[i])); - status = 1; - } else { - printf("real of [%d] is %g\n", i, creal(out_logical[i])); - } -} -``` +We therefore reuse the size-4`n` REDFT01 construction (eqs. 22–24 of the paper) on the *reversed* input +*X_reversed[j] = X[n−1−j]*, and recover the FFTW RODFT01 output as -The first `n` values should have the output of `RODFT01` in their negative imaginary parts: +*Y_k_FFTW(RODFT01) = 2 · S_k^III(X) − (-1)^k · X_{n−1} = (-1)^k · (2 · X̃_{2k+1} − X_{n−1})*. ```C -for (int i = 0; i < n; ++i) { - delta = -cimag(out_logical[i]) - out[i]; - if (fabs(delta) > eps) { - printf("error: delta of [%d] is %g\n", i, delta); - status = 1; - } else { - printf("match of [%d] (delta=%g)\n", i, delta); - } +int Nlog = 4 * n; +fill_zero_1d_cplx(Nlog, in_logical); +in_logical[0] = 0.5 * in[n - 1]; // X_reversed_0 = X_{n-1} +in_logical[2*n] = -0.5 * in[n - 1]; +for (int j = 1; j < n; ++j) { + double v = 0.25 * in[n - 1 - j]; // X_reversed_j = X_{n-1-j} + in_logical[j] = v; + in_logical[Nlog - j] = v; + in_logical[2*n - j] = -v; + in_logical[2*n + j] = -v; } -``` -The remaining values have odd symmetry around `n-0.5` -(note again that the real output from `RODFT01` needs to be compared against -the negative imaginary parts in the output of the logically equivalent DFT shifted by one full index): +dft_1d_cplx_unshifted(Nlog, in_logical, out_logical); -```C -for (int i = 0; i < n; ++i) { - delta = -cimag(out_logical[n + i]) - (-out[n - 1 - i]); +double delta; +for (int k = 0; k < n; ++k) { + double sign = (k % 2 == 0) ? 1.0 : -1.0; + delta = sign * (2.0 * creal(out_logical[2*k + 1]) - in[n - 1]) - out[k]; if (fabs(delta) > eps) { - printf("error: delta of [%d] is %g\n", n + i, delta); + printf("error: delta of [%d] is %g\n", k, delta); status = 1; } else { - printf("match of [%d] (delta=%g)\n", n + i, delta); + printf("match of [%d] (delta=%g)\n", k, delta); } } ``` @@ -938,79 +817,60 @@ The full example can be found in [`src/test_1d_rodft01.c`](src/test_1d_rodft01.c #### RODFT11 (DST-IV) In case of the real-valued odd-parity DFT with both shifted input and output data (RODFT11), -also called the DST-IV, the corresponding logical DFT size is given by *N* = 2`n`, corresponding to `n` = *N*/2. +also called the DST-IV, the appropriate logical-equivalent (unshifted) DFT size is *Ñ* = 8`n`. The formal definition of the RODFT11 is given below: ![RODFT11 formula](eqn/rodft11.png) +This is *FFTW's* form (prefactor 2). The Shao-Johnson DST-IV +([arXiv 0708.4399, eq. 19](https://arxiv.org/abs/0708.4399)) +is *S_k^IV = Σ X_j sin[π(j+½)(k+½)/n]* with prefactor 1; *Y_k = 2 · S_k^IV*. + The inverse of this transform is RODFT11 itself. The input array is assumed to have odd symmetry around *j=-0.5* and even symmetry around *j=n-0.5*. ![RODFT11](img/rodft11.png) -In order to demonstrate the use of this method, -the logically equivalent DFT input is filled appropriately and its output is checked against `RODFT11`. -In the following code, `in` is the input array (size `n`) given to `RODFT11` -and `in_logical` is the (complex-valued) input array (size *N*) handed to a -[generic 1D DFT](https://en.wikipedia.org/wiki/Discrete_Fourier_transform#Generalized_DFT_(shifted_and_non-linear_phase)). -Similarly, `out` is the output array (size `n`) from `RODFT11` -and `out_logical` is the output array (size *N*) from a generic 1D DFT. - -Here is how the symmetric input is generated: +The Shao-Johnson paper (§VII, eqs. 20–21) shows that the DST-IV is *exactly equivalent* to a DCT-IV +of a sign-flipped, output-reversed input — i.e. -```C -// the first half of the array is identical up to shift -for (int i = 0; i < n; ++i) { - in_logical[i] = in[i]; -} +> *S_{N-1-k}^IV(X) = Σ (-1)^n X_n cos[π(n+½)(k+½)/N] = C_k^IV(X')* with *X'_n = (-1)^n X_n*. -// second half is filled according to even symmetry around (n-0.5) -for (int i = 0; i < n; ++i) { - in_logical[n + i] = in[n - 1 - i]; -} -``` +(NB: the prefactor 2 printed in the published version of eq. 21 is a typo; verifying directly at +*N=2*, e.g. *X_0=1, X_1=0*, recovers *S_1 = cos(π/8) ≈ 0.924*, not *2 cos(π/8)*.) -The checks are a little bit more involved. -The logically equivalent DFT output should be purely imaginary-valued: +We therefore reuse the size-8`n` REDFT11 construction on *X'_n = (-1)^n X_n*, then read out +*Y_k_RODFT11 = 2·Y_{n−1−k}_REDFT11(X') = 2·X̃_{2(n−1−k)+1}*: ```C -for (int i = 0; i < N; ++i) { - if (fabs(creal(out_logical[i])) > eps) { - printf("error: real of [%d] is %g\n", i, creal(out_logical[i])); - status = 1; - } else { - printf("real of [%d] is %g\n", i, creal(out_logical[i])); - } +int Nlog = 8 * n; +fill_zero_1d_cplx(Nlog, in_logical); +for (int j = 0; j < n; ++j) { + double sign = (j % 2 == 0) ? 1.0 : -1.0; + double v = sign * 0.25 * in[j]; + in_logical[2*j + 1] = v; + in_logical[Nlog - 2*j - 1] = v; + in_logical[4*n - 2*j - 1] = -v; + in_logical[4*n + 2*j + 1] = -v; } -``` - -The first `n` values should have the output of `RODFT11` in their negative imaginary parts: -```C -for (int i = 0; i < n; ++i) { - delta = -cimag(out_logical[i]) - out[i]; - if (fabs(delta) > eps) { - printf("error: delta of [%d] is %g\n", i, delta); - status = 1; - } else { - printf("match of [%d] (delta=%g)\n", i, delta); - } -} +dft_1d_cplx_unshifted(Nlog, in_logical, out_logical); ``` -The remaining values have even symmetry around `n-0.5` -(note again that the real output from `RODFT11` needs to be compared against -the negative imaginary parts in the output of the logically equivalent DFT): +The DFT output is purely real (same four-fold symmetry as REDFT11) and only odd-index entries are non-zero. +The FFTW RODFT11 output is recovered by reverse-ordering: ```C -for (int i = 0; i < n; ++i) { - delta = -cimag(out_logical[n + i]) - out[n - 1 - i]; +double delta; +for (int k = 0; k < n; ++k) { + int idx = 2*(n - 1 - k) + 1; + delta = 2.0 * creal(out_logical[idx]) - out[k]; if (fabs(delta) > eps) { - printf("error: delta of [%d] is %g\n", n + i, delta); + printf("error: delta of [%d] is %g\n", k, delta); status = 1; } else { - printf("match of [%d] (delta=%g)\n", n + i, delta); + printf("match of [%d] (delta=%g)\n", k, delta); } } ``` @@ -1617,28 +1477,36 @@ if (n <= 0) { This code is abbreviated in the following as `/* compute idx_in */`. For the first coefficients with `m=0`, -the copying is performed as follows: +the copying is performed as follows. +Note that for `m=0` the VMEC basis collapses to *cos(nζ)* / *sin(nζ)* (with positive *nζ*, in contrast to *m>0* which uses *(mθ − nζ)*), +which is exactly the same situation as in [`src/app_magn_axis.c`](src/app_magn_axis.c) (the 1D magnetic-axis example), +so we apply the same Euler decomposition: +*cos(nζ) = ½ e^{inζ} + ½ e^{-inζ}*, *sin(nζ) = −(i/2) e^{inζ} + (i/2) e^{-inζ}*. +For *n=0* the DC entry is just *R^cos_{0,0}* (and 0 for `in_Z` since *sin(0)=0*). +For *n≥1*, because the toroidal dimension is the *full* (unhalved) one in this 2D `c2r` +(only the poloidal dimension is halved), we write the Hermitian-conjugate pair *(+n, −n)* explicitly: ```C int idx_vmec = 0; int m = 0; for (int n = 0; n <= ntor; ++n) { - /* compute idx_in */ - - in_R[idx_in] = rmnc[idx_vmec]; - in_Z[idx_in] = I * zmns[idx_vmec]; + if (n == 0) { + in_R[m] = rmnc[idx_vmec]; + in_Z[m] = 0.0; // sin(0) = 0 + } else { + int idx_pos = n * nyq_pol + m; // +n frequency + int idx_neg = (n_zeta - n) * nyq_pol + m; // -n frequency + in_R[idx_pos] = 0.5 * rmnc[idx_vmec]; + in_R[idx_neg] = 0.5 * rmnc[idx_vmec]; // conj of real = real + in_Z[idx_pos] = -0.5 * I * zmns[idx_vmec]; + in_Z[idx_neg] = 0.5 * I * zmns[idx_vmec]; // conj of -iZ/2 = +iZ/2 + } idx_vmec++; } ``` -Note that VMEC weights the coefficients for `m=0` with positive `n ζ` -in contrast to the following coefficients with `m>0`. -The cosine-parity quantities like *R* are not influence by this, -but the sine-parity quantities like *Z* need to get the sign of their value flipped -due to the odd symmetry `sin(-n zeta) = -sin(n zeta)`. - The second part of the copying adresses the Fourier coefficients with `m>0`: @@ -1657,9 +1525,13 @@ for (m = 1; m < mpol; ++m) { ``` Here, the coefficients are scaled by a factor of `0.5`. -FFTW implies that coefficientss for `m<0` were present in the logically equivalent DFT input, -which is not the case for VMEC output. Thus, they are (in this case errornously) -scaled internally by a factor of `2`, which needs to be counteracted by a scaling factor of `0.5` in the input. +The reason is the standard Euler decomposition of the basis function: +*cos(mθ - nζ) = ½ e^{i(mθ - nζ)} + ½ e^{-i(mθ - nζ)}*. +For *m>0*, only the first half — at FFTW input index *(m, -n)* — is stored explicitly; +the conjugate counterpart at *(-m, n)* is implicit because the poloidal dimension is the halved +(`n_theta/2+1`) one in the 2D `c2r`, and FFTW reconstructs the upper half via Hermitian conjugation. +Both halves contribute equally to the `c2r` output, so each stored coefficient must carry the prefactor +*0.5* in order for the sum of the two halves to recover the full *R^cos_{m,n} cos(mθ - nζ)*. Assuming [stellarator symmetry](https://doi.org/10.1016/S0167-2789(97)00216-9), half of the Fourier coefficients can be omitted and the transform reduces to the two-dimensional IDCT and IDST: diff --git a/src/app_flux_surface.c b/src/app_flux_surface.c index 6a0867d..34919e6 100644 --- a/src/app_flux_surface.c +++ b/src/app_flux_surface.c @@ -78,20 +78,23 @@ void app_flux_surface(int n_theta, int n_zeta) { int m = 0; for (int n = 0; n <= ntor; ++n) { - // compute target index for input array - if (n <= 0) { - idx_in = -n * nyq_pol + m; + // m=0 row: VMEC sums only over positive n with basis cos(nζ) / sin(nζ). + // Use the same Euler decomposition as in app_magn_axis.c (1D c2r): + // cos(nζ) = ½ e^{inζ} + ½ e^{-inζ} sin(nζ) = -i/2 e^{inζ} + i/2 e^{-inζ} + // The toroidal dimension n_zeta is the *full* (unhalved) one in this 2D c2r, + // so we explicitly set up Hermitian symmetry by writing both the +n and -n bins. + if (n == 0) { + in_R[m] = rmnc[idx_vmec]; + in_Z[m] = 0.0; // sin(0) = 0 } else { - idx_in = (n_zeta - n) * nyq_pol + m; + int idx_pos = n * nyq_pol + m; // +n frequency + int idx_neg = (n_zeta - n) * nyq_pol + m; // -n frequency + in_R[idx_pos] = 0.5 * rmnc[idx_vmec]; + in_R[idx_neg] = 0.5 * rmnc[idx_vmec]; // conj of real = real + in_Z[idx_pos] = -0.5 * I * zmns[idx_vmec]; + in_Z[idx_neg] = 0.5 * I * zmns[idx_vmec]; // conj of -iZ/2 = +iZ/2 } - // VMEC: for m=0, only positive ntor --> summed up within VMEC definition with negative-n --> 2*0.5 = 1 - // also for sin-parity quantities: m=0 would imply (m theta -n zeta) == (-n zeta), but VMEC uses (n zeta) for m=0 - // --> does not matter for cos-parity since cos(-n zeta) = cos(n zeta), but sin(-n zeta) = -sin(n zeta) - // and hence the additional -1 cancels out by compensating the -1 from the complex multiply in the DFT definition - in_R[idx_in] = rmnc[idx_vmec]; - in_Z[idx_in] = I * zmns[idx_vmec]; - idx_vmec++; } diff --git a/src/test_1d_redft01.c b/src/test_1d_redft01.c index a5a4759..58b2c03 100644 --- a/src/test_1d_redft01.c +++ b/src/test_1d_redft01.c @@ -9,75 +9,81 @@ int test_1d_redft01(int n) { - int N = 2 * n; + // logical-equivalent DFT size, per Shao & Johnson 2008 (arXiv 0708.4399, eqs. 22-24) + int Nlog = 4 * n; double *in = fftw_alloc_real(n); - fftw_complex *in_logical = fftw_alloc_complex(N); + fftw_complex *in_logical = fftw_alloc_complex(Nlog); double *out = fftw_alloc_real(n); - fftw_complex *out_logical = fftw_alloc_complex(N); + fftw_complex *out_logical = fftw_alloc_complex(Nlog); fftw_plan p = fftw_plan_r2r_1d(n, in, out, FFTW_REDFT01, FFTW_ESTIMATE); fill_random_1d_real(n, in); - // the first half of the array is identical - for (int i = 0; i < n; ++i) { - in_logical[i] = in[i]; - } - - // second half is filled according to odd symmetry around n - for (int i = 1; i < n; ++i) { - in_logical[n + i] = -in[n - i]; + // Build the size-4n logical input x̃ per Shao & Johnson 2008 (paper line 705-715): + // x̃_0 = X_0/2, x̃_n = 0, x̃_{2n} = -X_0/2, x̃_{3n} = 0 + // x̃_j = x̃_{4n - j} = +(1/4) X_j for 0 < j < n + // x̃_{2n - j} = x̃_{2n + j} = -(1/4) X_j for 0 < j < n + // Symmetry axes: even at j=0 and j=2n; odd at j=n and j=3n. + fill_zero_1d_cplx(Nlog, in_logical); + in_logical[0] = 0.5 * in[0]; + in_logical[2*n] = -0.5 * in[0]; + // in_logical[n] = 0; (already zero) + // in_logical[3*n] = 0; (already zero) + for (int j = 1; j < n; ++j) { + double v = 0.25 * in[j]; + in_logical[j] = v; + in_logical[Nlog - j] = v; + in_logical[2*n - j] = -v; + in_logical[2*n + j] = -v; } dump_1d_real("test_1d_redft01_in.dat", n, in); - dump_1d_cplx("test_1d_redft01_in_logical.dat", N, in_logical); + dump_1d_cplx("test_1d_redft01_in_logical.dat", Nlog, in_logical); - // manual implementation of logically-equivalent DFT - double a = 0.0; - double b = 0.5; - dft_1d_cplx(N, in_logical, out_logical, a, b); + // unshifted DFT of size 4n + dft_1d_cplx_unshifted(Nlog, in_logical, out_logical); - // use FFTW to execute REDFT01 + // FFTW REDFT01 of the n original samples fftw_execute(p); dump_1d_real("test_1d_redft01_out.dat", n, out); - dump_1d_cplx("test_1d_redft01_out_logical.dat", N, out_logical); + dump_1d_cplx("test_1d_redft01_out_logical.dat", Nlog, out_logical); - // check output of logically equivalent DFT double eps = 1e-12; int status = 0; - // 1. logically equivalent output should be purely real-valued - for (int i = 0; i < N; ++i) { - if (fabs(cimag(out_logical[i])) > eps) { - printf("error: imag of [%d] is %g\n", i, cimag(out_logical[i])); + // 1. (2k+1)-th DFT outputs are real (the four-fold symmetry of x̃ implies this) + for (int k = 0; k < n; ++k) { + if (fabs(cimag(out_logical[2*k + 1])) > eps) { + printf("error: imag of Y_logical[%d] is %g\n", 2*k + 1, cimag(out_logical[2*k + 1])); status = 1; } else { - printf("imag of [%d] is %g\n", i, cimag(out_logical[i])); + printf("imag of Y_logical[%d] is %g\n", 2*k + 1, cimag(out_logical[2*k + 1])); } } - // 2. first n values should have identical real values - double delta; - for (int i = 0; i < n; ++i) { - delta = creal(out_logical[i]) - out[i]; - if (fabs(delta) > eps) { - printf("error: delta of [%d] is %g\n", i, delta); + // 2. even-index DFT outputs vanish + for (int k = 0; k < 2*n; ++k) { + if (fabs(creal(out_logical[2*k])) > eps || fabs(cimag(out_logical[2*k])) > eps) { + printf("error: Y_logical[%d] (even-index) is (%g, %g) -- expected zero\n", + 2*k, creal(out_logical[2*k]), cimag(out_logical[2*k])); status = 1; - } else { - printf("match of [%d] (delta=%g)\n", i, delta); } } - // 3. even symmetry of output values around n-0.5 - for (int i = 0; i < n; ++i) { - delta = creal(out_logical[n + i]) - out[n - 1 - i]; + // 3. FFTW REDFT01 output Y_k equals 2 * X̃_{2k+1} - X_0. + // (FFTW form: Y_k = X_0 + 2 sum_{j>=1} X_j cos[...] = 2 C_k^III - X_0, + // while the Shao-Johnson DCT-III sums uniformly from j=0: C_k^III = X̃_{2k+1}.) + double delta; + for (int k = 0; k < n; ++k) { + delta = (2.0 * creal(out_logical[2*k + 1]) - in[0]) - out[k]; if (fabs(delta) > eps) { - printf("error: delta of [%d] is %g\n", n + i, delta); + printf("error: delta of [%d] is %g\n", k, delta); status = 1; } else { - printf("match of [%d] (delta=%g)\n", n + i, delta); + printf("match of [%d] (delta=%g)\n", k, delta); } } diff --git a/src/test_1d_redft10.c b/src/test_1d_redft10.c index 207087a..8a6c328 100644 --- a/src/test_1d_redft10.c +++ b/src/test_1d_redft10.c @@ -9,84 +9,75 @@ int test_1d_redft10(int n) { - int N = 2 * n; + // logical-equivalent DFT size for the half-sample-input-shifted DCT-II: + // Ñ = 4n, with data zero-interleaved (real samples at odd positions only) + // and even symmetry around 0 and 2n. Same family as Shao & Johnson eqs. 22-24 + // for DCT-III but with the shift on the input side rather than the output side. + int Nlog = 4 * n; double *in = fftw_alloc_real(n); - fftw_complex *in_logical = fftw_alloc_complex(N); + fftw_complex *in_logical = fftw_alloc_complex(Nlog); double *out = fftw_alloc_real(n); - fftw_complex *out_logical = fftw_alloc_complex(N); + fftw_complex *out_logical = fftw_alloc_complex(Nlog); fftw_plan p = fftw_plan_r2r_1d(n, in, out, FFTW_REDFT10, FFTW_ESTIMATE); fill_random_1d_real(n, in); - // the first half of the array is identical - for (int i = 0; i < n; ++i) { - in_logical[i] = in[i]; - } - - // second half is filled according to even symmetry around n-0.5 - for (int i = 0; i < n; ++i) { - in_logical[n + i] = in[n - 1 - i]; + // Build x̃ of size 4n with X_j at positions 2j+1 and 4n-2j-1 (mirror), zero elsewhere. + // Symmetry: even around j=0 and even around j=2n. + fill_zero_1d_cplx(Nlog, in_logical); + for (int j = 0; j < n; ++j) { + in_logical[2*j + 1] = in[j]; + in_logical[Nlog - 2*j - 1] = in[j]; } dump_1d_real("test_1d_redft10_in.dat", n, in); - dump_1d_cplx("test_1d_redft10_in_logical.dat", N, in_logical); + dump_1d_cplx("test_1d_redft10_in_logical.dat", Nlog, in_logical); - // manual implementation of logically-equivalent DFT - double a = 0.5; - double b = 0.0; - dft_1d_cplx(N, in_logical, out_logical, a, b); + // unshifted DFT of size 4n + dft_1d_cplx_unshifted(Nlog, in_logical, out_logical); - // use FFTW to execute REDFT10 + // FFTW REDFT10 of the original n samples fftw_execute(p); dump_1d_real("test_1d_redft10_out.dat", n, out); - dump_1d_cplx("test_1d_redft10_out_logical.dat", N, out_logical); + dump_1d_cplx("test_1d_redft10_out_logical.dat", Nlog, out_logical); - // check output of logically equivalent DFT double eps = 1e-12; int status = 0; - // 1. logically equivalent output should be purely real-valued - for (int i = 0; i < N; ++i) { - if (fabs(cimag(out_logical[i])) > eps) { - printf("error: imag of [%d] is %g\n", i, cimag(out_logical[i])); + // 1. all DFT outputs are real (consequence of even symmetry of x̃ around 0) + for (int k = 0; k < Nlog; ++k) { + if (fabs(cimag(out_logical[k])) > eps) { + printf("error: imag of Y_logical[%d] is %g\n", k, cimag(out_logical[k])); status = 1; - } else { - printf("imag of [%d] is %g\n", i, cimag(out_logical[i])); } } - // 2. first n values should have identical real values + // 2. FFTW REDFT10 output Y_k equals Re(X̃_k) for k = 0, ..., n-1. + // The FFTW prefactor of 2 is built into the construction (each X_j appears + // twice in x̃: once at 2j+1 and once at the mirror 4n-2j-1; their cosine + // contributions combine to give the FFTW formula). double delta; - for (int i = 0; i < n; ++i) { - delta = creal(out_logical[i]) - out[i]; + for (int k = 0; k < n; ++k) { + delta = creal(out_logical[k]) - out[k]; if (fabs(delta) > eps) { - printf("error: delta of [%d] is %g\n", i, delta); + printf("error: delta of [%d] is %g\n", k, delta); status = 1; } else { - printf("match of [%d] (delta=%g)\n", i, delta); + printf("match of [%d] (delta=%g)\n", k, delta); } } - // 3. odd symmetry of output values around n - // (implies that value at n is zero) + // 3. X̃_n is zero (DCT-II "Nyquist" boundary, follows from even-around-2n) if (fabs(creal(out_logical[n])) > eps) { - printf("error: delta of [%d] is %g\n", n, creal(out_logical[n])); + printf("error: Re(Y_logical[%d]) = %g (expected 0 by DCT-II symmetry)\n", + n, creal(out_logical[n])); status = 1; } else { printf("match of [%d] (delta=%g)\n", n, creal(out_logical[n])); } - for (int i = 1; i < n; ++i) { - delta = creal(out_logical[n + i]) - (-out[n - i]); - if (fabs(delta) > eps) { - printf("error: delta of [%d] is %g\n", n + i, delta); - status = 1; - } else { - printf("match of [%d] (delta=%g)\n", n + i, delta); - } - } if (status == 0) { printf("=> all ok\n"); diff --git a/src/test_1d_redft11.c b/src/test_1d_redft11.c index 75453b8..adb45ee 100644 --- a/src/test_1d_redft11.c +++ b/src/test_1d_redft11.c @@ -9,75 +9,82 @@ int test_1d_redft11(int n) { - int N = 2 * n; + // logical-equivalent DFT size, per Shao & Johnson 2008 (arXiv 0708.4399, eq. 5) + int Nlog = 8 * n; double *in = fftw_alloc_real(n); - fftw_complex *in_logical = fftw_alloc_complex(N); + fftw_complex *in_logical = fftw_alloc_complex(Nlog); double *out = fftw_alloc_real(n); - fftw_complex *out_logical = fftw_alloc_complex(N); + fftw_complex *out_logical = fftw_alloc_complex(Nlog); fftw_plan p = fftw_plan_r2r_1d(n, in, out, FFTW_REDFT11, FFTW_ESTIMATE); fill_random_1d_real(n, in); - // the first half of the array is identical - for (int i = 0; i < n; ++i) { - in_logical[i] = in[i]; - } - - // second half is filled according to odd symmetry around n-0.5 - for (int i = 0; i < n; ++i) { - in_logical[n + i] = -in[n - 1 - i]; + // Build the size-8n logical input x̃ per Shao & Johnson 2008 + // (arXiv 0708.4399, eqs. 6-7): + // x̃[2j+1] = x̃[8n - 2j - 1] = +(1/4) X_j for 0 <= j < n + // x̃[4n - 2j - 1] = x̃[4n + 2j + 1] = -(1/4) X_j for 0 <= j < n + // x̃[2j] = 0 for 0 <= j < 4n + // Symmetry axes: even at j=0, odd at j=2n, even at j=4n, odd at j=6n + // (paper line 720-721: "even around 0 and N̄/2, odd around N̄/4 and 3N̄/4"). + fill_zero_1d_cplx(Nlog, in_logical); + for (int j = 0; j < n; ++j) { + double v = 0.25 * in[j]; + in_logical[2*j + 1] = v; + in_logical[Nlog - 2*j - 1] = v; + in_logical[4*n - 2*j - 1] = -v; + in_logical[4*n + 2*j + 1] = -v; } dump_1d_real("test_1d_redft11_in.dat", n, in); - dump_1d_cplx("test_1d_redft11_in_logical.dat", N, in_logical); + dump_1d_cplx("test_1d_redft11_in_logical.dat", Nlog, in_logical); - // manual implementation of logically-equivalent DFT - double a = 0.5; - double b = 0.5; - dft_1d_cplx(N, in_logical, out_logical, a, b); + // logically-equivalent (unshifted) DFT of size 8n + dft_1d_cplx_unshifted(Nlog, in_logical, out_logical); - // use FFTW to execute REDFT11 + // FFTW REDFT11 of the n original samples fftw_execute(p); dump_1d_real("test_1d_redft11_out.dat", n, out); - dump_1d_cplx("test_1d_redft11_out_logical.dat", N, out_logical); + dump_1d_cplx("test_1d_redft11_out_logical.dat", Nlog, out_logical); - // check output of logically equivalent DFT double eps = 1e-12; int status = 0; - // 1. logically equivalent output should be purely real-valued - for (int i = 0; i < N; ++i) { - if (fabs(cimag(out_logical[i])) > eps) { - printf("error: imag of [%d] is %g\n", i, cimag(out_logical[i])); + // 1. the (2k+1)-th DFT outputs are purely real + // (the four-fold symmetry of x̃ implies that the entire DFT is real-valued; + // we report the imaginary parts to make this visible, and the + // even-index outputs are zero by the same symmetry — checked below.) + for (int k = 0; k < n; ++k) { + if (fabs(cimag(out_logical[2*k + 1])) > eps) { + printf("error: imag of Y_logical[%d] is %g\n", 2*k + 1, cimag(out_logical[2*k + 1])); status = 1; } else { - printf("imag of [%d] is %g\n", i, cimag(out_logical[i])); + printf("imag of Y_logical[%d] is %g\n", 2*k + 1, cimag(out_logical[2*k + 1])); } } - // 2. first n values should have identical real values - double delta; - for (int i = 0; i < n; ++i) { - delta = creal(out_logical[i]) - out[i]; - if (fabs(delta) > eps) { - printf("error: delta of [%d] is %g\n", i, delta); + // 2. even-index DFT outputs vanish (a consequence of the zero-interleaved x̃) + for (int k = 0; k < 4*n; ++k) { + if (fabs(creal(out_logical[2*k])) > eps || fabs(cimag(out_logical[2*k])) > eps) { + printf("error: Y_logical[%d] (even-index) is (%g, %g) -- expected zero\n", + 2*k, creal(out_logical[2*k]), cimag(out_logical[2*k])); status = 1; - } else { - printf("match of [%d] (delta=%g)\n", i, delta); } } - // 3. odd symmetry of output values around n-0.5 - for (int i = 0; i < n; ++i) { - delta = creal(out_logical[n + i]) - (-out[n - 1 - i]); + // 3. FFTW REDFT11 output Y_k equals 2 * X̃_{2k+1}. + // The unnormalized Shao-Johnson DCT-IV is C_k^IV = X̃_{2k+1}; + // FFTW prepends a factor of 2 (uniform convention across all REDFT/RODFT types). + double delta; + for (int k = 0; k < n; ++k) { + delta = 2.0 * creal(out_logical[2*k + 1]) - out[k]; if (fabs(delta) > eps) { - printf("error: delta of [%d] is %g\n", n + i, delta); + printf("error: delta of [%d] is %g\n", k, delta); status = 1; } else { - printf("match of [%d] (delta=%g)\n", n + i, delta); + printf("match of [%d] (delta=%g)\n", k, delta); } } diff --git a/src/test_1d_rodft01.c b/src/test_1d_rodft01.c index ef77042..723cb5f 100644 --- a/src/test_1d_rodft01.c +++ b/src/test_1d_rodft01.c @@ -9,78 +9,85 @@ int test_1d_rodft01(int n) { - int N = 2 * n; + // logical-equivalent DFT size, same Ñ = 4n picture as REDFT01 (DCT-III). + // RODFT01 is reduced to a DCT-III via Shao & Johnson eq. 32 (paper line 906-908): + // S_k^III(X) = (-1)^k * sum_{m=0}^{N-1} X_{N-m} cos[π m (k+1/2)/N] + // = (-1)^k * C_k^III(X_reversed) + // where X_reversed[m] = X[n-1-m] (using FFTW's 0-based indexing). + int Nlog = 4 * n; double *in = fftw_alloc_real(n); - fftw_complex *in_logical = fftw_alloc_complex(N); + fftw_complex *in_logical = fftw_alloc_complex(Nlog); double *out = fftw_alloc_real(n); - fftw_complex *out_logical = fftw_alloc_complex(N); + fftw_complex *out_logical = fftw_alloc_complex(Nlog); fftw_plan p = fftw_plan_r2r_1d(n, in, out, FFTW_RODFT01, FFTW_ESTIMATE); fill_random_1d_real(n, in); - // a zero in the first entry is needed to satisfy - // odd symmetry about -1 in the shifted input array - in_logical[0] = 0.0; - - // the first half of the array is identical up to shift - for (int i = 0; i < n; ++i) { - in_logical[i + 1] = in[i]; - } - - // second half is filled according to even symmetry around (n-1) - for (int i = 0; i < n - 1; ++i) { - in_logical[n + 1 + i] = in[n - 2 - i]; + // Build the size-4n logical input x̃ for the REVERSED input array, using the + // Shao & Johnson DCT-III construction (eqs. 22-24). + fill_zero_1d_cplx(Nlog, in_logical); + in_logical[0] = 0.5 * in[n - 1]; // x̃_0 = X_reversed_0 / 2 = X_{n-1} / 2 + in_logical[2*n] = -0.5 * in[n - 1]; + for (int j = 1; j < n; ++j) { + double v = 0.25 * in[n - 1 - j]; // X_reversed_j = X_{n-1-j} + in_logical[j] = v; + in_logical[Nlog - j] = v; + in_logical[2*n - j] = -v; + in_logical[2*n + j] = -v; } dump_1d_real("test_1d_rodft01_in.dat", n, in); - dump_1d_cplx("test_1d_rodft01_in_logical.dat", N, in_logical); + dump_1d_cplx("test_1d_rodft01_in_logical.dat", Nlog, in_logical); - // manual implementation of logically-equivalent DFT - double a = 0.0; - double b = 0.5; - dft_1d_cplx(N, in_logical, out_logical, a, b); + // unshifted DFT of size 4n + dft_1d_cplx_unshifted(Nlog, in_logical, out_logical); + // FFTW RODFT01 on the original n samples fftw_execute(p); dump_1d_real("test_1d_rodft01_out.dat", n, out); - dump_1d_cplx("test_1d_rodft01_out_logical.dat", N, out_logical); + dump_1d_cplx("test_1d_rodft01_out_logical.dat", Nlog, out_logical); - // check output of logically equivalent DFT double eps = 1e-12; int status = 0; - // 1. logically equivalent output should be purely imaginary-valued - for (int i = 0; i < N; ++i) { - if (fabs(creal(out_logical[i])) > eps) { - printf("error: real of [%d] is %g\n", i, creal(out_logical[i])); + // 1. (2k+1)-th DFT outputs are real + for (int k = 0; k < n; ++k) { + if (fabs(cimag(out_logical[2*k + 1])) > eps) { + printf("error: imag of Y_logical[%d] is %g\n", 2*k + 1, cimag(out_logical[2*k + 1])); status = 1; } else { - printf("real of [%d] is %g\n", i, creal(out_logical[i])); + printf("imag of Y_logical[%d] is %g\n", 2*k + 1, cimag(out_logical[2*k + 1])); } } - // 2. next n values should have the output of RODFT01 in the negative imaginary part - double delta; - for (int i = 0; i < n; ++i) { - delta = -cimag(out_logical[i]) - out[i]; - if (fabs(delta) > eps) { - printf("error: delta of [%d] is %g\n", i, delta); + // 2. even-index DFT outputs vanish + for (int k = 0; k < 2*n; ++k) { + if (fabs(creal(out_logical[2*k])) > eps || fabs(cimag(out_logical[2*k])) > eps) { + printf("error: Y_logical[%d] (even-index) is (%g, %g) -- expected zero\n", + 2*k, creal(out_logical[2*k]), cimag(out_logical[2*k])); status = 1; - } else { - printf("match of [%d] (delta=%g)\n", i, delta); } } - // 3. odd symmetry around (n-0.5) of neg. imag. part vs. output of RODFT01 - for (int i = 0; i < n; ++i) { - delta = -cimag(out_logical[n + i]) - (-out[n - 1 - i]); + // 3. FFTW RODFT01 output reconstructed from the DCT-III of the reversed input. + // Per paper eq. 32 (using X_paper_n = X_FFTW_{n-1} so the paper's sum_{n=1}^{N} + // becomes the FFTW j=0..n-1 sum): + // S_k^III(X) = (-1)^k * C_k^III(X_reversed) = (-1)^k * X̃_{2k+1} + // FFTW's RODFT01 has the same prefactor relationship to S_k^III as REDFT01 has + // to C_k^III: Y_k_FFTW = 2 S_k^III - (-1)^k X_{n-1} + // = (-1)^k * (2 * X̃_{2k+1} - X_{n-1}). + double delta; + for (int k = 0; k < n; ++k) { + double sign = (k % 2 == 0) ? 1.0 : -1.0; + delta = sign * (2.0 * creal(out_logical[2*k + 1]) - in[n - 1]) - out[k]; if (fabs(delta) > eps) { - printf("error: delta of [%d] is %g\n", n + i, delta); + printf("error: delta of [%d] is %g\n", k, delta); status = 1; } else { - printf("match of [%d] (delta=%g)\n", n + i, delta); + printf("match of [%d] (delta=%g)\n", k, delta); } } diff --git a/src/test_1d_rodft10.c b/src/test_1d_rodft10.c index 1a31f26..c196209 100644 --- a/src/test_1d_rodft10.c +++ b/src/test_1d_rodft10.c @@ -9,83 +9,70 @@ int test_1d_rodft10(int n) { - int N = 2 * n; + // logical-equivalent DFT size for the half-sample-input-shifted DST-II: + // Ñ = 4n, with data zero-interleaved (real samples at odd positions only) + // and odd-around-0 / odd-around-2n symmetry. Sine analog of REDFT10. + int Nlog = 4 * n; double *in = fftw_alloc_real(n); - fftw_complex *in_logical = fftw_alloc_complex(N); + fftw_complex *in_logical = fftw_alloc_complex(Nlog); double *out = fftw_alloc_real(n); - fftw_complex *out_logical = fftw_alloc_complex(N); + fftw_complex *out_logical = fftw_alloc_complex(Nlog); fftw_plan p = fftw_plan_r2r_1d(n, in, out, FFTW_RODFT10, FFTW_ESTIMATE); fill_random_1d_real(n, in); - // the first half of the array is identical - for (int i = 0; i < n; ++i) { - in_logical[i] = in[i]; - } - - // second half is filled according to odd symmetry around (n-0.5) - for (int i = 0; i < n; ++i) { - in_logical[n + i] = -in[n - 1 - i]; + // Build x̃ of size 4n with X_j at position 2j+1 and -X_j at the mirror 4n-2j-1. + // Symmetries: odd-around-0 and odd-around-2n; all even positions zero. + fill_zero_1d_cplx(Nlog, in_logical); + for (int j = 0; j < n; ++j) { + in_logical[2*j + 1] = in[j]; + in_logical[Nlog - 2*j - 1] = -in[j]; } dump_1d_real("test_1d_rodft10_in.dat", n, in); - dump_1d_cplx("test_1d_rodft10_in_logical.dat", N, in_logical); + dump_1d_cplx("test_1d_rodft10_in_logical.dat", Nlog, in_logical); - // manual implementation of logically-equivalent DFT - double a = 0.5; - double b = 0.0; - dft_1d_cplx(N, in_logical, out_logical, a, b); + // unshifted DFT of size 4n + dft_1d_cplx_unshifted(Nlog, in_logical, out_logical); + // FFTW RODFT10 of the original n samples fftw_execute(p); dump_1d_real("test_1d_rodft10_out.dat", n, out); - dump_1d_cplx("test_1d_rodft10_out_logical.dat", N, out_logical); + dump_1d_cplx("test_1d_rodft10_out_logical.dat", Nlog, out_logical); - // check output of logically equivalent DFT double eps = 1e-12; int status = 0; - // 1. logically equivalent output should be purely imaginary-valued - for (int i = 0; i < N; ++i) { - if (fabs(creal(out_logical[i])) > eps) { - printf("error: real of [%d] is %g\n", i, creal(out_logical[i])); + // 1. all DFT outputs are purely imaginary (consequence of odd symmetry of x̃ around 0) + for (int k = 0; k < Nlog; ++k) { + if (fabs(creal(out_logical[k])) > eps) { + printf("error: real of Y_logical[%d] is %g\n", k, creal(out_logical[k])); status = 1; - } else { - printf("real of [%d] is %g\n", i, creal(out_logical[i])); } } - // 2. first output value is zero (?) - if (fabs(cimag(out_logical[0])) > eps) { - printf("error: delta of [%d] is %g\n", 0, -cimag(out_logical[0])); - status = 1; - } else { - printf("match of [%d] (delta=%g)\n", 0, -cimag(out_logical[0])); - } - - // 3. next n values should have the output of RODFT10 in the negative imaginary part + // 2. FFTW RODFT10 output Y_k equals -Im(X̃_{k+1}) for k = 0, ..., n-1. + // Derivation: X̃_k = -2i Σ_j X_j sin(π(j+1/2) k / n), so X̃_{k+1} carries the FFTW + // sin[π(j+1/2)(k+1)/n] factor; Y_k_FFTW (prefactor 2) = -Im(X̃_{k+1}). double delta; - for (int i = 0; i < n; ++i) { - delta = -cimag(out_logical[i + 1]) - out[i]; + for (int k = 0; k < n; ++k) { + delta = -cimag(out_logical[k + 1]) - out[k]; if (fabs(delta) > eps) { - printf("error: delta of [%d] is %g\n", i + 1, delta); + printf("error: delta of [%d] is %g\n", k, delta); status = 1; } else { - printf("match of [%d] (delta=%g)\n", i + 1, delta); + printf("match of [%d] (delta=%g)\n", k, delta); } } - // 4. even symmetry around (n-1) of neg. imag. part vs. output of RODFT10 - for (int i = 0; i < n - 1; ++i) { - delta = -cimag(out_logical[n + 1 + i]) - out[n - 2 - i]; - if (fabs(delta) > eps) { - printf("error: delta of [%d] is %g\n", n + 1 + i, delta); - status = 1; - } else { - printf("match of [%d] (delta=%g)\n", n + 1 + i, delta); - } + // 3. X̃_0 = 0 (the DC bin vanishes by the odd-around-0 symmetry of x̃) + if (fabs(cimag(out_logical[0])) > eps) { + printf("error: Im(Y_logical[0]) = %g (expected 0 by odd-around-0)\n", + cimag(out_logical[0])); + status = 1; } if (status == 0) { diff --git a/src/test_1d_rodft11.c b/src/test_1d_rodft11.c index c9408d7..5f19dbc 100644 --- a/src/test_1d_rodft11.c +++ b/src/test_1d_rodft11.c @@ -9,74 +9,79 @@ int test_1d_rodft11(int n) { - int N = 2 * n; + // logical-equivalent DFT size, same Ñ = 8n picture as REDFT11 + int Nlog = 8 * n; double *in = fftw_alloc_real(n); - fftw_complex *in_logical = fftw_alloc_complex(N); + fftw_complex *in_logical = fftw_alloc_complex(Nlog); double *out = fftw_alloc_real(n); - fftw_complex *out_logical = fftw_alloc_complex(N); + fftw_complex *out_logical = fftw_alloc_complex(Nlog); fftw_plan p = fftw_plan_r2r_1d(n, in, out, FFTW_RODFT11, FFTW_ESTIMATE); fill_random_1d_real(n, in); - // the first half of the array is identical up to shift - for (int i = 0; i < n; ++i) { - in_logical[i] = in[i]; - } - - // second half is filled according to even symmetry around (n-0.5) - for (int i = 0; i < n; ++i) { - in_logical[n + i] = in[n - 1 - i]; + // Shao & Johnson 2008 (arXiv 0708.4399), §VII: a DST-IV is exactly equivalent to + // a DCT-IV in which the outputs are reversed and every other input is multiplied by -1 + // (see eq. 20 of the paper, which derives S_{N-1-k}^IV = sum_n (-1)^n X_n cos[...]). + // We therefore reuse the DCT-IV construction (eqs. 5-7) on the sign-flipped input + // X'_j = (-1)^j X_j + // and recover Y_k_RODFT11 = 2 * X̃_{2(n-1-k)+1} = 2 * X̃_{2n-1-2k}. + fill_zero_1d_cplx(Nlog, in_logical); + for (int j = 0; j < n; ++j) { + double sign = (j % 2 == 0) ? 1.0 : -1.0; + double v = sign * 0.25 * in[j]; + in_logical[2*j + 1] = v; + in_logical[Nlog - 2*j - 1] = v; + in_logical[4*n - 2*j - 1] = -v; + in_logical[4*n + 2*j + 1] = -v; } dump_1d_real("test_1d_rodft11_in.dat", n, in); - dump_1d_cplx("test_1d_rodft11_in_logical.dat", N, in_logical); + dump_1d_cplx("test_1d_rodft11_in_logical.dat", Nlog, in_logical); - // manual implementation of logically-equivalent DFT - double a = 0.5; - double b = 0.5; - dft_1d_cplx(N, in_logical, out_logical, a, b); + // unshifted size-8n DFT + dft_1d_cplx_unshifted(Nlog, in_logical, out_logical); + // FFTW RODFT11 on the original n samples fftw_execute(p); dump_1d_real("test_1d_rodft11_out.dat", n, out); - dump_1d_cplx("test_1d_rodft11_out_logical.dat", N, out_logical); + dump_1d_cplx("test_1d_rodft11_out_logical.dat", Nlog, out_logical); - // check output of logically equivalent DFT double eps = 1e-12; int status = 0; - // 1. logically equivalent output should be purely imaginary-valued - for (int i = 0; i < N; ++i) { - if (fabs(creal(out_logical[i])) > eps) { - printf("error: real of [%d] is %g\n", i, creal(out_logical[i])); + // 1. (2k+1)-th DFT outputs are purely real (same four-fold symmetry as DCT-IV) + for (int k = 0; k < n; ++k) { + if (fabs(cimag(out_logical[2*k + 1])) > eps) { + printf("error: imag of Y_logical[%d] is %g\n", 2*k + 1, cimag(out_logical[2*k + 1])); status = 1; } else { - printf("real of [%d] is %g\n", i, creal(out_logical[i])); + printf("imag of Y_logical[%d] is %g\n", 2*k + 1, cimag(out_logical[2*k + 1])); } } - // 2. next n values should have the output of RODFT11 in the negative imaginary part - double delta; - for (int i = 0; i < n; ++i) { - delta = -cimag(out_logical[i]) - out[i]; - if (fabs(delta) > eps) { - printf("error: delta of [%d] is %g\n", i, delta); + // 2. even-index DFT outputs vanish + for (int k = 0; k < 4*n; ++k) { + if (fabs(creal(out_logical[2*k])) > eps || fabs(cimag(out_logical[2*k])) > eps) { + printf("error: Y_logical[%d] (even-index) is (%g, %g) -- expected zero\n", + 2*k, creal(out_logical[2*k]), cimag(out_logical[2*k])); status = 1; - } else { - printf("match of [%d] (delta=%g)\n", i, delta); } } - // 3. even symmetry around (n-0.5) of neg. imag. part vs. output of RODFT11 - for (int i = 0; i < n; ++i) { - delta = -cimag(out_logical[n + i]) - out[n - 1 - i]; + // 3. FFTW RODFT11 output Y_k equals 2 * X̃_{2(n-1-k)+1} + // (DST-IV recovered from DCT-IV via output reversal; FFTW's prefactor of 2 vs. SJ's prefactor 1). + double delta; + for (int k = 0; k < n; ++k) { + int idx = 2*(n - 1 - k) + 1; + delta = 2.0 * creal(out_logical[idx]) - out[k]; if (fabs(delta) > eps) { - printf("error: delta of [%d] is %g\n", n + i, delta); + printf("error: delta of [%d] is %g\n", k, delta); status = 1; } else { - printf("match of [%d] (delta=%g)\n", n + i, delta); + printf("match of [%d] (delta=%g)\n", k, delta); } } diff --git a/src/util.h b/src/util.h index ee0c221..263765e 100644 --- a/src/util.h +++ b/src/util.h @@ -291,6 +291,13 @@ void dft_1d_cplx(int n, fftw_complex *in, fftw_complex *out, double a, double b) } } +// plain unshifted DFT (a = b = 0 in dft_1d_cplx) +// used by the REDFT/RODFT examples to compute the logical-equivalent DFT +// of a zero-interleaved + 4-fold symmetric extension of size N (4n or 8n) +void dft_1d_cplx_unshifted(int N, fftw_complex *in, fftw_complex *out) { + dft_1d_cplx(N, in, out, 0.0, 0.0); +} + int read_lcfs_header(char *filename, int *mpol, int *ntor, int *mnmax, int *nfp) { char *line = NULL; char *mpol_line = "# mpol\n";