BUG:special:amos: Fix some mistakes in the AMOS C translation

[inkydragon]

Reference issue

• related pr: ENH:MAINT:special:Rewrite amos F77 code #19587 (comment)

What does this implement/fix?

1. Fixed missing return branch of amos_kscl function in C translation of amos code.
   See description below.
2. Remove extra else branch in amos_kscl.
   See my first comment for detailed description.

C translation:

scipy/scipy/special/_amos.c

Lines 4007 to 4011 in 44b39e1

	if (n == 2) {
	return nz;
	}
	fn = fnu + 1.;
	ck = fn*rz;

ref: amos fortran code:

scipy/scipy/special/amos/zkscl.f

Lines 52 to 56 in e5f2129

	20 CONTINUE
	IF (N.EQ.2) RETURN
	IF (NZ.EQ.0) RETURN
	FN = FNU + 1.0D0
	CKR = FN*RZR

Additional information

• Add some test
  Since this function is not exported, it is difficult to construct a test case to cover this bug.

Some shorter call paths from exported functions to amos_kscl.

amos_airy 
	-> amos_bknu -> amos_kscl

amos_besh, amos_besk
	-> amos_acon
	-> amos_bknu -> amos_kscl

full call graph

It's a bit difficult to test this bug in python, but you can test it in C.
Add the following C code to _amos.c, compile it, and run it.
Then check the printout of the program.

C test code add to _amos.c

run with gcc -g -o amos _amos.c -lm && ./amos

#include <stdio.h>
#define ARRAY_SIZE(x) ((sizeof x) / (sizeof *x))
int main () {
    int n = 5;
    double complex y[5] = { 0.0 };
    int y_size = ARRAY_SIZE(y);
    for (int i=0; i < y_size; i++) {
        y[i] = (i+1) + i*I;
    }
    double complex zr = 1e-9 + 0.0*I;
    double fnu = 2.0;
    double complex rz = 0.5 + 0.5*I;
    double ascle = 1e-30;
    double tol = 1e-15;
    double elim = 20.;

    printf("---- Input: y[%d] =\n", y_size);
    for (int i=0; i < n; i++) {
        printf( "  y[%d] = %g + %g*I\n",
            i, creal(y[i]), cimag(y[i]));
    }
    printf("zr = %g + %g*I;", creal(zr), cimag(zr));
    printf("  fnu=%g;  n=%d\n", fnu, n);
    printf("rz = %g + %g*I;", creal(rz), cimag(rz));
    printf("  ascle=%g;  tol=%g;  elim=%g\n",
        ascle, tol, elim);

    int nz = amos_kscl(zr,fnu,n,y,rz,&ascle,tol,elim);
    puts("---- Output");
    printf("nz=%d;  y[%d]=\n", nz, y_size);
    for (int i=0; i < y_size; i++) {
        printf( "  y[%d] = %g + %g*I\n",
            i, creal(y[i]), cimag(y[i]));
    }

    return 0;
}
#undef ARRAY_SIZE

C test code outout

master output (Without pr)

• nz = 1

$ gcc -O0 -g -o amos _amos.c -lm && ./amos
---- Input: y[5] =
  y[0] = 1 + 0*I
  y[1] = 2 + 1*I
  y[2] = 3 + 2*I
  y[3] = 4 + 3*I
  y[4] = 5 + 4*I
zr = 1e-09 + 0*I;  fnu=2;  n=5
rz = 0.5 + 0.5*I;  ascle=1e-30;  tol=1e-15;  elim=20
---- Output
nz=1;  y[5]=
  y[0] = 0 + 0*I
  y[1] = 2e+15 + 1e+15*I
  y[2] = 2.5e+15 + 4.5e+15*I
  y[3] = 4 + 3*I
  y[4] = 5 + 4*I

with pr

• nz = 0, and y[] outout also changed.

$ gcc -O0 -g -o amos _amos.c -lm && ./amos
---- Input: y[5] =
  y[0] = 1 + 0*I
  y[1] = 2 + 1*I
  y[2] = 3 + 2*I
  y[3] = 4 + 3*I
  y[4] = 5 + 4*I
zr = 1e-09 + 0*I;  fnu=2;  n=5
rz = 0.5 + 0.5*I;  ascle=1e-30;  tol=1e-15;  elim=20
---- Output
nz=0;  y[5]=
  y[0] = 1e+15 + 0*I
  y[1] = 2e+15 + 1e+15*I
  y[2] = 3 + 2*I
  y[3] = 4 + 3*I
  y[4] = 5 + 4*I

[inkydragon]

There appears to be another problem in the amos_kscl function:

Fortran code

scipy/scipy/special/amos/zkscl.f

Lines 109 to 121 in e5f2129

	30 CONTINUE
	NZ = N
	IF(IC.EQ.N) NZ=N-1
	GO TO 45
	40 CONTINUE
	NZ = KK - 2
	45 CONTINUE
	DO 50 I=1,NZ
	YR(I) = ZEROR
	YI(I) = ZEROI
	50 CONTINUE
	RETURN
	END

Regardless of whether the if condition IF(IC.EQ.N) is true or not, we will always jump to 45.
Maybe we don't need else branch here. cc: @ilayn

C translation

scipy/scipy/special/_amos.c

Lines 4052 to 4059 in 44b39e1

	nz = n;
	if (ic == n) {
	nz = n-1;
	} else {
	nz = kk - 2;
	}
	for (int i = 0; i < nz; i++) { y[i] = 0.; }
	return nz;

 if (ic == n) { 
     nz = n-1; 
-} else { 
-    nz = kk - 2; 
 } 

Notice that NZ = KK - 2 is merged into another return path.

scipy/scipy/special/_amos.c

Lines 4035 to 4041 in 44b39e1

	nz -= 1;
	if (ic == kk-1) {
	nz = kk - 2;
	for (int i = 0; i < nz; i++) { y[i] = 0.; }
	return nz;
	}
	ic = kk;

scipy/scipy/special/amos/zkscl.f

Lines 98 to 100 in e5f2129

	NZ = NZ - 1
	IF (IC.EQ.KK-1) GO TO 40
	IC = KK

[ilayn]

Yes it seems like the else branch is not needed. The main issue, probably it seems, is that, I associated the if condition with the jump. So it seems redundant there.

[ilayn]

Ah how nice that you added the call graph. @lucascolley if the failing backend test is not relevant let's get this one in when CI is green please.

[steppi]

Thanks for catching these @inkydragon. The alternative backend test is definitely not relevant. I'll merge when rest of the CI settles.

[lucascolley]

thanks for the ping, that is probably due to the update of array_api_compat that Ralf just merged - will have a look

[lucascolley]

thanks for your first contribution to SciPy @inkydragon !