BUG:special:amos: Fix exit path in amos_asyi

[inkydragon]

Reference issue

xref: #19587 cc @ilayn

What does this implement/fix?

• Fixed an issue where a return branch in amos_asyi was always false

Additional information

C:

scipy/scipy/special/_amos.c

Lines 1078 to 1095 in b882f1b

	dk = ez;
	j = 1;
	for (int j = 1; j < (jl+1); j++)
	{
	ck *= sqk / dk;
	cs2 += ck;
	sgn = -sgn;
	cs1 += ck*sgn;
	dk += ez;
	aa *= fabs(sqk) / bb;
	bb += aez;
	ak += 8.;
	sqk -= ak;
	if (aa <= atol) { break; }
	}
	if ((j == jl) && (aa > atol)) { return -2; }
	/* 50 */

The value of the first expression, (j == jl), on line 1093 is independent of how the for loop exits.
Because the for loop introduces a new variable scope, modifying the loop variable j used by for does not affect the j variable of the same name outside for.

So in fact there is always j==1.

My modification introduces a new flag to check if you're exiting from break, and if so, continue the function.
Otherwise, it returns an error code.

There are a number of ways to fix this, perhaps using aa > atol directly as a judgment condition for exiting the function would work without introducing a new variable, I'm not sure.

Fortran:

scipy/scipy/special/amos/zasyi.f

Lines 99 to 117 in e5f2129

	DO 40 J=1,JL
	CALL ZDIV(CKR, CKI, DKR, DKI, STR, STI)
	CKR = STR*SQK
	CKI = STI*SQK
	CS2R = CS2R + CKR
	CS2I = CS2I + CKI
	SGN = -SGN
	CS1R = CS1R + CKR*SGN
	CS1I = CS1I + CKI*SGN
	DKR = DKR + EZR
	DKI = DKI + EZI
	AA = AA*DABS(SQK)/BB
	BB = BB + AEZ
	AK = AK + 8.0D0
	SQK = SQK - AK
	IF (AA.LE.ATOL) GO TO 50
	40 CONTINUE
	GO TO 110
	50 CONTINUE

scipy/scipy/special/amos/zasyi.f

Lines 162 to 164 in e5f2129

	110 CONTINUE
	NZ=-2
	RETURN

[ilayn]

Actually this is because of that extra in the int word in the loop. I was still fighting with VS Code while writing this part to not to correct things like this.

for (int j = 1; j < (jl+1); j++) 

Without that, j retains its value so you don't need the flag. Also the condition should be

if ((j == jl+1) && (aa > atol)) { return -2; }

because C loop increments j before exiting the loop. In other words the condition here is checking; Loop is exhausted but aa is still larger than atol, return -2.

In Python, loop exhaustion check is wonderful but in these languages, they are left to the developer unfortunately.

[inkydragon]

We might also be missing a return

diff --git a/scipy/special/_amos.c b/scipy/special/_amos.c
index 7629d350b..7bc59cb3d 100644
--- a/scipy/special/_amos.c
+++ b/scipy/special/_amos.c
@@ -1115,7 +1115,10 @@ int amos_asyi(
         if (koded == 0) { return nz; }
         ck = cexp(cz);
         for (int i = 0; i < (nn + 1); i++) { y[i] *= ck; }
+        /* 90 */
+        return nz;
     }
+    /* 100 */
     return -1;
 }

scipy/scipy/special/amos/zasyi.f

Lines 153 to 161 in e5f2129

	DO 90 I=1,NN
	STR = YR(I)*CKR - YI(I)*CKI
	YI(I) = YR(I)*CKI + YI(I)*CKR
	YR(I) = STR
	90 CONTINUE
	RETURN
	100 CONTINUE
	NZ = -1
	RETURN

[ilayn]

We might also be missing a return

Yes I think so too.

By the way, thanks a lot for these @inkydragon ! Keep them coming as you encounter them please. A bit embarrassing on my side but crucial fixes.

[inkydragon]

A bit embarrassing on my side but crucial fixes.

I think this is normal, given the number of lines of code in _amos.c, and the complex function call relationships, it's easy for there to be areas that the tests don't cover.
It's also common for bugs to hide in untested code.

I'm trying to rewrite the amos code in Julia based on your C translation. The way I'm doing it is to start rewriting from the leaf functions that have no dependencies, and then test the new code against the original fortran code and try to improve the test coverage.

So these bugs were found through testing, and I probably wouldn't have been able to find them if I had just stared at the code.

Even with line coverage close to 100%, I was able to find some inconsistencies, mainly in special floating point inputs such as NAN, INF, etc. This could mean that my implementation is still buggy, or that the way fortran handles special floating point numbers is inconsistent with julia.

---

I would recommend my rewrite approach: @ilayn

1. Compile the target fortran library to a dynamic library
2. Compile the rewritten C code into a dynamic library
3. Call both libraries in python using ffi and compare the output of the two libraries on a large number of inputs
4. Furthermore, you can measure the test coverage of the C code.

[ilayn]

Actually, except from item 4, about C coverage, that's how I did it too. But I don't really know how to get the result of gcov triggered from pytest mechanism yet that SciPy uses yet.

Regardless, thanks again for these fixes, let's get this one in.

[j-bowhay]

@ilayn it is fairly painless to get, you do a build with
python dev.py build --gcov and then run the tests as normal then the output ends up in build/**/.gc

[ilayn]

I must be doing something wrong then since I couldn't manage to get things work but thanks for the heads up, I'll try again.

[steppi]

I must be doing something wrong then since I couldn't manage to get things work but thanks for the heads up, I'll try again.

I haven't been able to get it to work either. Let me know if you get it working and if so what you did differently.

[j-bowhay]

What part specifically isn't working? I checked the coverage on one of these recent PRs so pretty sure I had it working

[steppi]

I'm trying to rewrite the amos code in Julia based on your C translation. The way I'm doing it is to start rewriting from the leaf functions that have no dependencies, and then test the new code against the original fortran code and try to improve the test coverage.

I think it's really awesome to see this translation might be helpful for the Julia ecosystem too. It's actually something I'd mentioned in a recent grant proposal asking for funding to rewrite SciPy's F77 code. It's very cool to see that really start taking shape. Good luck in your efforts @inkydragon. Thanks for taking a close look and helping to find these bugs.

[ilayn]

What part specifically isn't working? I checked the coverage on one of these recent PRs so pretty sure I had it working

Two parts; I tried both with the dev.py route that does not generate the gc files and also specifically compiled scipy by modifying meson.build file adding the -lgcov arguments to special libraries. That generates the gcno and other files but can't be used to see the coverage info. Also probably, I am a little bit running out of steam with all kinds of C stuff so probably I screwed up somewhere.

[ilayn]

I think it's really awesome to see this translation might be helpful for the Julia ecosystem too.

Indeed, that's the whole point of all why I got into this and frankly, rejuvenated my motivation to finish things off. I know that there were multiple attempts in the past while I was researching about the appetite such as JuliaMath/SpecialFunctions.jl#52

Thank you indeed.

[j-bowhay]

What part specifically isn't working? I checked the coverage on one of these recent PRs so pretty sure I had it working

Two parts; I tried both with the dev.py route that does not generate the gc files and also specifically compiled scipy by modifying meson.build file adding the -lgcov arguments to special libraries. That generates the gcno and other files but can't be used to see the coverage info. Also probably, I am a little bit running out of steam with all kinds of C stuff so probably I screwed up somewhere.

Feel free to ping me on future PRs and I can just upload the coverage artefacts