The structure of fftrn is almost identical to fftr1 and the basic equations used are the same. The key difference comes from the fact that in more than one dimension the element can no longer simply be found at the location . Rather the index must be reversed from to separately for each dimension. Likewise in the one-dimensional case only needed to include non-negative value of whereas in this case that is true only in the last dimension.
This indexing is particularly tricky since fftrn treats the data as a one-dimensional array. (There is no practical way to treat the data as a multi-dimensional array without knowing in advance how many dimensions are going to be used.) Fortunately there is a simple solution to this problem. The function fftrn defines an array called that keeps track of the separate indices in all the dimensions. At the bottom of the main loop these indices are incremented, advancing the last one until it reaches its maximum and then starting it over and incrementing the next to last one, and so on. At the top of the main loop is a call to the function find_index, which finds the one-dimensional array index that corresponds to these different multi-dimensional indices. More importantly, the function find_index also finds the one-dimensional index corresponding to the negative of the frequencies indicated by . From there the implementation is almost identical to fftr1. The inner loop runs over the last dimension and applies the exact same formulas as fftr1 to calculate , the desired transform, from , the output of fftcn. I leave it as an exercise for the reader to show that these formulas are still valid provided the frequencies in all but the last dimension are replaced by their corresponding negative frequencies as shown.
The and (Nyquist) frequencies present a problem. In one dimension these both had to be real because . In multiple dimensions, however, it is no longer true that all values of for which the last index are or are real. To illustrate what I mean let me use to indicate the frequency in the last dimension and to indicate the frequency components in all other dimensions. Now the symmetry is simply . This means that applying equation  to the case will no longer yield a real number as it did in the one-dimensional case and so can no longer be packaged with into one complex array index. So fftrn takes as input an array that it uses for the values and uses both the real and imaginary parts of the first elements of the array for . In practice this means that the output is redundant because it will contain and , which are complex conjugates of each other (and likewise for and in ).
As with fftr1 the algebra for the inverse transform turns out to be virtually identical to the forward case, and once again I leave the details to the reader to verify. Note that it is up to the user to make sure that the symmetries noted above are obeyed by the input to the inverse transform. In section [4.3] I provide sample code for generating data with all the correct symmetries in three dimensions. You can simply cut and paste this sample code into your own program and substitute whatever function you want to generate in frequency space.