BUG: missing complex lqmn derivatives when |z|<1.001

[jorenham]

Closes scipy/scipy#23589

[jorenham]

My VSCode autoformatter was a bit too enthusiastic, but it should be fixed now

[lucascolley]

I promise I haven't hacked your machine to include my changes in your commits

[jorenham]

Erm, I can't seem to find the lqmn tests; are there any?

No xsref tables either 🤔

[steppi]

Thanks @jorenham. That's a good catch. As for missing tests and tables, we're actually missing tests for all gufunc kernels. The process I followed to autogenerate a test suite based on the scipy.special tests only worked for ufunc kernels. I had intended to follow up adding tests for gufunc kernels, but ran out of funding for such work, and until recently, out of bandwidth for unfunded work. I'm trying to triage tasks based on the bandwidth I have, and tests for gufunc kernels is a very important one. I'll try to get that through over the next few weeks. In the meantime, we can just add a few hand coded tests here.

[jorenham]

I wrote smoke tests for lqmn. I based two on the python tests in scipy (https://github.com/scipy/scipy/blob/bdd3b0e77a3813c22c038c908d992b6de23ffcda/scipy/special/tests/test_legendre.py#L693-L708), and the other is the reproducer (real vs complex) for this issue.

[jorenham]

And FWIW, this fix makes it look more like the original fortran code (src):

        SUBROUTINE CLQMN(MM,M,N,X,Y,CQM,CQD)
C
C       =======================================================
C       Purpose: Compute the associated Legendre functions of
C                the second kind, Qmn(z) and Qmn'(z), for a
C                complex argument
C       Input :  x  --- Real part of z
C                y  --- Imaginary part of z
C                m  --- Order of Qmn(z)  ( m = 0,1,2,… )
C                n  --- Degree of Qmn(z) ( n = 0,1,2,… )
C                mm --- Physical dimension of CQM and CQD
C       Output:  CQM(m,n) --- Qmn(z)
C                CQD(m,n) --- Qmn'(z)
C       =======================================================
C
        IMPLICIT DOUBLE PRECISION (X,Y)
        IMPLICIT COMPLEX*16 (C,Z)
        DIMENSION CQM(0:MM,0:N),CQD(0:MM,0:N)
        Z = DCMPLX(X, Y)
        IF (DABS(X).EQ.1.0D0.AND.Y.EQ.0.0D0) THEN
           DO 10 I=0,M
           DO 10 J=0,N
              CQM(I,J)=(1.0D+300,0.0D0)
              CQD(I,J)=(1.0D+300,0.0D0)
10         CONTINUE
           RETURN
        ENDIF
        XC=ABS(Z)
        LS=0
        IF (DIMAG(Z).EQ.0.0D0.OR.XC.LT.1.0D0) LS=1
        IF (XC.GT.1.0D0) LS=-1
        ZQ=SQRT(LS*(1.0D0-Z*Z))
        ZS=LS*(1.0D0-Z*Z)
        CQ0=0.5D0*LOG(LS*(1.0D0+Z)/(1.0D0-Z))
        IF (XC.LT.1.0001D0) THEN
           CQM(0,0)=CQ0
           CQM(0,1)=Z*CQ0-1.0D0
           CQM(1,0)=-1.0D0/ZQ
           CQM(1,1)=-ZQ*(CQ0+Z/(1.0D0-Z*Z))
           DO 15 I=0,1
           DO 15 J=2,N
              CQM(I,J)=((2.0D0*J-1.0D0)*Z*CQM(I,J-1)
     &                -(J+I-1.0D0)*CQM(I,J-2))/(J-I)
15         CONTINUE
           DO 20 J=0,N
           DO 20 I=2,M
              CQM(I,J)=-2.0D0*(I-1.0D0)*Z/ZQ*CQM(I-1,J)-LS*
     &                 (J+I-1.0D0)*(J-I+2.0D0)*CQM(I-2,J)
20         CONTINUE
        ELSE
           IF (XC.GT.1.1) THEN
              KM=40+M+N
           ELSE
              KM=(40+M+N)*INT(-1.0-1.8*LOG(XC-1.0))
           ENDIF
           CQF2=(0.0D0,0.0D0)
           CQF1=(1.0D0,0.0D0)
           DO 25 K=KM,0,-1
              CQF0=((2*K+3.0D0)*Z*CQF1-(K+2.0D0)*CQF2)/(K+1.0D0)
              IF (K.LE.N) CQM(0,K)=CQF0
              CQF2=CQF1
25            CQF1=CQF0
           DO 30 K=0,N
30            CQM(0,K)=CQ0*CQM(0,K)/CQF0
           CQF2=0.0D0
           CQF1=1.0D0
           DO 35 K=KM,0,-1
              CQF0=((2*K+3.0D0)*Z*CQF1-(K+1.0D0)*CQF2)/(K+2.0D0)
              IF (K.LE.N) CQM(1,K)=CQF0
              CQF2=CQF1
35            CQF1=CQF0
           CQ10=-1.0D0/ZQ
           DO 40 K=0,N
40            CQM(1,K)=CQ10*CQM(1,K)/CQF0
           DO 45 J=0,N
              CQ0=CQM(0,J)
              CQ1=CQM(1,J)
              DO 45 I=0,M-2
                 CQF=-2.0D0*(I+1)*Z/ZQ*CQ1+(J-I)*(J+I+1.0D0)*CQ0
                 CQM(I+2,J)=CQF
                 CQ0=CQ1
                 CQ1=CQF
45         CONTINUE
        ENDIF
        CQD(0,0)=LS/ZS
        DO 50 J=1,N
50         CQD(0,J)=LS*J*(CQM(0,J-1)-Z*CQM(0,J))/ZS
        DO 55 J=0,N
        DO 55 I=1,M
           CQD(I,J)=LS*I*Z/ZS*CQM(I,J)+(I+J)*(J-I+1.0D0)
     &              /ZQ*CQM(I-1,J)
55      CONTINUE
        RETURN
        END

[jorenham]

I have no idea what the macos and windows test failures are about though (first time writing non-trivial C++) 😅. Could it be a caching thing, or did I mess up some obvious C pointer thingy or something? Is there a way for me to repro this on my linux (ubuntu) machine?

[steppi]

I have no idea what the macos and windows test failures are about though (first time writing non-trivial C++) 😅. Could it be a caching thing, or did I mess up some obvious C pointer thingy or something? Is there a way for me to repro this on my linux (ubuntu) machine?

Thanks @jorenham for writing the tests. I'll take a look at the failures in the next couple days.

[jorenham]

I have no idea what the macos and windows test failures are about though (first time writing non-trivial C++) 😅. Could it be a caching thing, or did I mess up some obvious C pointer thingy or something? Is there a way for me to repro this on my linux (ubuntu) machine?

I managed to figure it out myself, and it turns out that lqmn apparently indexes the first 2x2 indices of the output matrices, even if it's smaller (m < 1 || n < 1) than that (here and here). On windows this was causing a build failure (stricter compiler or something), and on macos this was causing a runtime crash (without reporting an error code or something). No idea why it was running without problems on linux though.

Either way, I worked around the by increasing the buffer size accordingly if needed, and side-stepped the windows build failure by using 🦆-typed std::mdspan classes.

But as far as I could tell, this lqmn indexing issue is also present in the original Fortran code (here and here). So if this is actually an issue that should be fixed, doing so in a separate PR would have my preference, since it's independent of the one addressed in this PR.

[steppi]

Nice work @jorenham! Rather than fixing lqmn, what I want to do is push forward with rewriting lqmn in the style @izaid used for assoc_legendre_p and assoc_legendre_p_all. I want to get the proper tests set up first though.

[jorenham]

The windows test failures might have something to do with that Windows Server 2025 rollout: actions/runner-images#12677

[jorenham]

It turns out that updating kokkos/mdspan to the latest (0.6) version fixes the windows build errors. I'll open a separate PR for that, and mark this as draft in the meantime

[jorenham]

The windows tests are failing "in the right way" now 😛

[steppi]

Thanks for adding the tests @jorenham. I'm a little curious why the existing SciPy tests use an atol and not an rtol but think it's fine to just copy them over here like this. I'll look into having more systematic tests for legendre functions and other gufuncs.

[steppi]

The windows tests are failing "in the right way" now 😛

Hoping to have a time to investigate fixing the underlying issue soon. It's due to an inconsistency in complex arithmetic on Windows. If I remember right, something to do with them not completely following the spec.