I have asked this question some time ago and in the meantime I gained some insights into the problem. Thus I will try to address the questions I had and explain the relationships between continuous problem, discrete problem and FFT algorithm. I would also like to note that I have implemented this in this CUDA program.
First of all, because of the fact that we have just $N \times M$ values of the function $f$ at the gridpoints and we would like to solve continuous equation
$$ \Delta u = f,$$
we would need the mean we understand the objects $f$ and $u$. How $f$ will be defined outside the gridpoints and in which points of the rectangle $[0,A] \times [0,B]$ we actually place gridpoints. I will be using notation, where I identify $N$, $n$ with axis $x$ and $M$, $m$ with axis $y$.
We aim to construct continuous object $f$, which coincides with $f_{nm}$ at the grid points
$$f(n\frac{A}{N}, m \frac{B}{M}) = f_{nm}, \quad n \in \{0...N-1\}, m\in \{0 ... M-1\}.$$
Let's do FFT of $f_{nm}$ to obtain $F_{kl}$, now we have that
$$f_{mn} = \sum_{k=0}^{N-1} \sum_{l=0}^{M-1} F_{kl} e^{2\pi i (k n/N + l m/M)}$$
and thus
$$f(n\frac{A}{N}, m \frac{B}{M}) = \sum_{k=0}^{N-1} \sum_{l=0}^{M-1} F_{kl} e^{2\pi i (k n/N + l m/M)}$$
can be extended to the $f$ defined on the whole domain by the following formula
$$f(x, y) = \sum_{k=0}^{N-1} \sum_{l=0}^{M-1} F_{kl} e^{2\pi i (k x/A + l y/B)}.$$
Partial answer 2:
This definition of continuous object $f$ is possible but probably not the best one. If we subscribe to this definition, the function $f$ will be real at the gridpoints, but outside of them it will be complex function in general. So we use the trick based on the following identity that holds for the discrete IFFT
$$f_{mn} = \sum_{k=0}^{N-1} \sum_{l=0}^{M-1} F_{kl} e^{2\pi i (k n/N + l m/M)}=\sum_{k=-\lfloor \frac{N - 1}{2}\rfloor }^{\lfloor \frac{N}{2}\rfloor} \sum_{l=-\lfloor\frac{M - 1}{2}\rfloor }^{\lfloor \frac{M}{2}\rfloor} F_{kl} e^{2\pi i (k n/N + l m/M)},$$
where we identify $F_{kl} = F_{k+M l+N}$. The identity holds because $e^{2\pi i} = e^{0} = 1$. For us however is more convenient to use the second sum in the identity to define continuous object $f$. When $N$ and $M$ will be odd, then the function
$$f(x, y) = \sum_{k=-\lfloor \frac{N - 1}{2}\rfloor }^{\lfloor \frac{N}{2}\rfloor} \sum_{l=-\lfloor\frac{M - 1}{2}\rfloor }^{\lfloor \frac{M}{2}\rfloor} F_{kl} e^{2\pi i (k x/A + l y/B)}$$
will also be purely real and can also be represented as a expansion by means of trigonometric functions. This interpolation of our grid by the sum of trigonometric functions is actually the object we would like to use to solve our differential equation. Analogously we define the function $u$
$$u(x, y) = \sum_{k=-\lfloor \frac{N - 1}{2}\rfloor }^{\lfloor \frac{N}{2}\rfloor} \sum_{l=-\lfloor\frac{M - 1}{2}\rfloor }^{\lfloor \frac{M}{2}\rfloor} U_{kl} e^{2\pi i (k x/A + l y/B)},$$
and therefore its Laplace will have the form
$$\Delta u(x, y) = -(2 \pi)^2\sum_{k=-\lfloor \frac{N - 1}{2}\rfloor }^{\lfloor \frac{N}{2}\rfloor} \sum_{l=-\lfloor\frac{M - 1}{2}\rfloor }^{\lfloor \frac{M}{2}\rfloor} \left(\left(\frac{k}{A}\right)^2+\left(\frac{l}{B}\right)^2\right) U_{kl} e^{2\pi i (k x/A + l y/B)}.$$
From the theory about Poisson's equation with periodic boundary conditions we know that $f$ and $u$ have to have the zero mean. This can be achieved only by first ensuring that the input data $f_{nm}$ have zero mean, second by fixing $U_{00}=F_{00}=0$.
Partial answer 1: Note what is now our discrete problem. It is to find a values of $u_{mn}$ on the gridpoints such that the Laplace operator on the $u$ defined as the "Fourier extension" the gridpoint values will coincide with the "Fourier extension" of the object $f_{mn}$. We don't discretize the differential operator on the grid. We rather construct Fourier or trigonometric extension of the functions involved and understand the derivatives in the point sense at the gridpoints as the derivatives of interpolating functions $f$ and $u$. From the restriction of the Poisson equation to the gridpoints we know, that we can find $u_{nm}$ by IFFT of $U_{kl}$, where
$$F_{kl} = -(2 \pi)^2 \left(\left(\frac{k}{A}\right)^2+\left(\frac{l}{B}\right)^2\right) U_{kl}.$$
From how we constructed $f$ and $u$ it is evident that the derivatives are understood as a point derivatives at the gridpoints of the interpolations by trigonometric functions. Thus they will not coincide with any finite difference scheme. The referenced article rather states that "when the problem is discretized by the finite difference (FD) method with a five-point stencil on an (n + 1) × (n + 1) equispaced grid, there is a FFT-based algorithm that computes the values of the solution... " I am not aware of such FFT-based algorithm, the referenced book has to be paid for, but for the algorithm I just have outlined the five-point stencil won't coincide with aforementioned definition of Laplacian on trigonometric extension. It does not exclude that it probably could be shown that these two definitions are in some sense close, e.g. by Taylor expansion.
Partial answer 2 - continue: There is also one more argument, that $k$ and $l$ up to $N-1$ and $M-1$ is not a good idea other that $f$ will be complex. The higher values of $|k|$ are there, the higher frequencies will be present in my interpolation of $f$. And the higher differences might be between point derivatives of $f$ and five-point stencil or whatever reasonable discretization of Laplace operator on my grid.
Partial answer 3: We have now the method how to solve this discrete problem by finding functions $f$, $u$ with the support on $[0,A] \times [0,B] $ with given values of $f$ on the gridpoints so that these two functions satisfy Poisson equation. It can be done for any $o_k \in \mathbb{N}$, $o_l \in \mathbb{N}$ and frequencies $k \in \{ o_k ... N - 1 - o_k\}$, $l \in \{ o_l ... M - 1 - o_l\}$. We have seen that by setting $o_k=-\lfloor \frac{M-1}{2} \rfloor$ and $o_l=-\lfloor \frac{N-1}{2} \rfloor$ we obtain real $f$ in odd case. In even case however, there will still be the term correspondig to $k=N/2$
$$ \sum_{l=-\lfloor\frac{M - 1}{2}\rfloor }^{\lfloor \frac{M}{2}\rfloor} F_{(N/2)l} e^{2\pi i (N/2 x/A + l y/B)}$$
which needs to be symmetrized to obtain real function $f$ in the way
$$ \frac{1}{2} \sum_{l=-\lfloor\frac{M - 1}{2}\rfloor }^{\lfloor \frac{M}{2}\rfloor} F_{(-N/2)l} e^{2\pi i (-N/2 x/A + l y/B)} + \frac{1}{2} \sum_{l=-\lfloor\frac{M - 1}{2}\rfloor }^{\lfloor \frac{M}{2}\rfloor} F_{(N/2)l} e^{2\pi i (N/2 x/A + l y/B)}$$
In fact it would need to be symmetrized in $l$ as well. By doing so we have found real $f$ even in unbalanced case. And if moreover we use $k$ and $l$ that are centered around $0$ from the formula
$$F_{kl} = -(2 \pi)^2 \left(\left(\frac{k}{A}\right)^2+\left(\frac{l}{B}\right)^2\right) U_{kl}$$
we also see that $U_{-kl}=\bar{U}_{kl}$ because of $F_{-kl}=\bar{F}_{kl}$ and because of the choice of frequencies I does not break hermitan symmetry.
Partial answer 4:
There is a theory regarding coercive operators, that for the problem with Dirichlet boundary conditions on the part of the boundary there is unique weak solution to the Poisson equation, for sources see this answer. It is natural to ask if we can use spectral methods to solve Poisson equation on square with Neumann or Dirichlet boundary conditions instead of periodic boundary conditions. I am not sure about general answer yet. It seems however that the simplicity of the spectral method stems from the fact, that functions in the Fourier decompositions are eigenfunctions of Laplace operator with periodic boundary conditions. Whenever we lose this, the method gets more ugly.
It is apparently possible to use something called odd or even extension to solve problem with homogeneous Dirichlet or homogeneous Neumann. Key observation here is that any sine function in FFT on $[0,2A]$ is zero at $0$, $A$ and $2A$ and any cosine function in FFT on $[0,2A]$ will have zero derivative at $0$, $A$ and $2A$. And if we construct odd function, it will have zero even components in Fourier transform and vice versa.
So let's say we would like to have homogeneous Dirichlet and we have
$$f(x,y), x \in [0,A], y\in [0, B]$$ and we extend it to $[0,2A] \times [0, 2B]$ such that
$$f(x,y) = \begin{cases}
f(x,y), & x \in [0,A], y \in [0,B],\\
-f(2A-x,y), & x \in [A,2A], y \in [0,B],\\
-f(x,2A-y), & x \in [0,A], y \in [B,2B],\\
f(2A-x,2A-y) & x \in [A,2A], y \in [B,2B].
\end{cases}$$
Now we solve the original problem $\Delta u = f$ with periodic boundary conditions on bigger domain $[0,2A] \times [0, 2B]$. The restriction of the solution to the $[0,A] \times [0, B]$ will be solution to the original problem with Dirichlet boundary conditions. See also Figure 3 of this paper. By solving four times bigger problem, we managed to use this type of spectral method to solve the following boundary value problems
\begin{equation}
\begin{aligned}
u(x,0) &= u(x,B) = 0, \quad x \in [0,A],\\
u(0,y) &= u(A,y)=0, \quad y \in [0,B].
\end{aligned}
\end{equation}
or
\begin{equation}
\begin{aligned}
u(x,0) &= u(x,B) = 0, \quad x \in [0,A],\\
u_x(0,y) &= u_x(A,y) = 0, \quad y \in [0,B].\\
\end{aligned}
\end{equation}
or
\begin{equation}
\begin{aligned}
u_y(x,0) &= u_y(x,B)=0, \quad x \in [0,A],\\
u_x(0,y) &= u_x(A,y)=0, \quad y \in [0,B].\\
\end{aligned}
\end{equation}
Of course that raises further question if there is a way to extend the set of problems even further. And if if it is necessary to solve four times bigger problem to have solution with some kind of normal boundary conditions or if we can computationally exploit the symmetry of the problem.
Apparently when homogeneous Dirichlet conditions are imposed and I need other type of Dirichlet conditions, it is sufficient to compute Laplace equation solution on the square with those Dirichlet conditions and sum the two solutions. This way the method can be generalized for any type of Dirichlet boundary conditions.
Note also when doing the odd or even extension we have to pay attention how FFT and IFFT of sequence [0,n-1] is defined with respect to the domain [0,A]. The last point $n-1$ is mapped not to A but to $A-A/N$ so if we do naïve extension by simple mirroring we obtain lot of aliasing between FFT and FT frequencies. From the perspective of keeping radix-2 algorithm it might be much faster than using different than $2N$ grid.
