Tackling some `sepfir2d` improvements

[rgommers]

From scipy#14367 (comment):

There are a couple of issues (2 open, 2 closed) about sepfir2d, so it does get usage. Now that we've looked into this and found some of the problems, it makes sense to fix them. This could turn into a bit of a project though, so let's not try to do it all in this PR. @Smit-create I suggest the following:

• This PR: incorporate the fix the outptr -> out in memmove (see review comment above), and then we'll merge this.
• Follow-up PRs (all separately, or at least the templating one - to keep diffs readable):
  • deduplicating the four files (C_bspline_util.c, D_bspline_util.c, S_bspline_util.c, Z_bspline_util.c). This can be done with a Tempita template (see scipy/signal/correlate_nd.c.in for an example). At that point bspline_util.c also doesn't need to be a separate file, the one function in it can be included in the same templated C file.
  • add a test for strided input arrays, and handle strides correctly (@peterbell10's comment above)
  • add a test for BUG: signal.sepfir2d: fails with complex input on Windows scipy/scipy#13643 and remove the __GNUC__ guards. I think this requires modernizing the code (__complex__ usage is odd, it's only done in a few signal functions)

This code is pretty arcane, it should be possible to squash quite a few bugs in it here.

EDIT: final thing could be to improve the algorithm to also handle even hrow, hcol, see scipygh-12691.

Next steps

Let's start with the easier topics here:

• deduplicate the four files to get only a single file
• remove the __GNUC__ guards. Instead of using __complex__ float, this could use numpy types like npy_cfloat.

[Smit-create]

Hi Ralf,
I removed the GNUC guards and used npy_cfloat/ npy_cdouble. There were few errors while building like:

scipy/signal/bspline_util.c: In function 'C_IIR_order1':
scipy/signal/bspline_util.c:50:16: error: invalid operands to binary * (have 'npy_cfloat' {aka 'struct <anonymous>'} and 'npy_cfloat' {aka 'struct <anonymous>'})
   50 |  *yvec = *xvec * a1 + *(yvec - stridey) * a2;
      |          ~~~~~ ^
      |          |
      |          npy_cfloat {aka struct <anonymous>}
scipy/signal/bspline_util.c:305:40: error: invalid operands to binary - (have 'npy_cfloat' {aka 'struct <anonymous>'} and 'double')
  305 |     *(y + (N - 1)*stridey) = -c0 / (z1 - 1.0) * yp[N-1];
      |                                        ^
scipy/signal/bspline_util.c: In function 'Z_IIR_forback1':
scipy/signal/bspline_util.c:329:13: error: incompatible types when assigning to type 'npy_cdouble' {aka 'struct <anonymous>'} from type 'double'
  329 |     powz1 = 1.0;
      |             ^~~
scipy/signal/bspline_util.c:334:12: error: invalid operands to binary * (have 'npy_cdouble' {aka 'struct <anonymous>'} and 'npy_cdouble' {aka 'struct <anonymous>'})
  334 |      powz1 *= z1;

They are because of binary operators. Are there any built-in functions such that numpy can handle them? Something like I found here: https://github.com/pydata/numexpr/blob/master/numexpr/complex_functions.hpp. Some of them can be useful from here https://github.com/numpy/numpy/blob/main/numpy/core/include/numpy/npy_math.h but regarding binary operands I didn't find any.

[rgommers]

There really aren't many examples of complex math, and pretty much everything we do have is in Cython. I'd suggest looking at scipy/special/_faddeeva.cxx as the most relevant example here.