BUG: special: algorithmic Error in ratevl in cephes/polevl.h

[ZhibingSun]

I am a PhD candidate in Applied Math from the University of Waterloo. I was investigating some numerical issues in my research when I found an algorithmic error for the function ratevl in scipy / scipy / special / cephes / polevl.h. It appears on line 158, which should be i = M - N; instead of i = N - M;. Here I give a brief explanation of what the function does and why there is a bug.
The function ratevl on line 118 evaluates a rational function $\cfrac{P(x)}{Q(x)}$, where $P(x)$ and $Q(x)$ are polynomials of order $M$ and $N$ respectively. For $|x| > 1$, $P(x)$ and $Q(x)$ are evaluated by $\cfrac{1}{x}$ instead of $x$, using Horner's method as referenced on line 55.
Here I give a simple example to illustrate. Let $f_2(x) = a_2 x^2 + a_1 x + a_0$ be a second order polynomial. If $|x| > 1$, then the polynomial is evaluated in reverse order as polynomial in $\cfrac{1}{x}$.
That is, $\cfrac{f_2(x)}{x^2} = a_2 + a_1 \cdot\cfrac{1}{x} + a_0 \cdot\cfrac{1}{x^2} = (a_0 \cdot\cfrac{1}{x} + a_1)\cdot \cfrac{1}{x} + a_2$, thus $f_2(x) = ((a_0\cdot \cfrac{1}{x} + a_1)\cdot\cfrac{1}{x} + a_2) \cdot x^2$. In the code, $P(x) = num\textunderscore ans \cdot x^M$, $Q(x) = denom\textunderscore ans \cdot x^N$, thus $\cfrac{P(x)}{Q(x)} = \cfrac{num\textunderscore ans}{denom\textunderscore ans}\cdot x^{M - N}$, and $i$ on line 158 should be $M - N$, not $N - M$.
In the boost implementation as referenced on line 55, it pre-processed to require $P(x)$ and $Q(x)$ to have the same order so that it's not an issue there. As on line 60, I believe the piece of code for the function ratevl is rewritten instead of copying from the reference directly, which is likely to be a reason for this error.

[lucascolley]

Hi @ZhibingSun , as a small ask, could you edit your post to enclose equations in backticks (this symbol `)? That should make it easier to read! (Or you can use $ signs for LaTeX too).

cc @steppi

[ZhibingSun]

Hi, I have updated my post as required, and thanks for letting me know that I can enclose equations and use LaTex. Please let me know if I need to provide more information.

[steppi]

Thanks for the bug report @ZhibingSun! It looks like this error goes back to when this function was first introduced. It seems no one noticed before you because ratevl is only used in one place in SciPy and in that case the numerator and denominator have the same degree.

scipy/scipy/special/special/cephes/lanczos.h

Line 83 in 3559245

SPECFUN_HOST_DEVICE double lanczos_sum(double x) { return ratevl(x, lanczos_num, 12, lanczos_denom, 12); }

Would you like submit a pull request with the fix?

[ZhibingSun]

Yes, I would love to contribute. Could you give me some guidance of how to submit a pull request?

[steppi]

Yes, I would love to contribute. Could you give me some guidance of how to submit a pull request?

Great! Start with the Contributor quickstart in the Developer documentation, and then check out the Development workflow section. That should be enough to get started.

[ZhibingSun]

I followed Contributor quickstart guide and run command upto python3 dev.py test, but it gives me error meson.build:84:0: ERROR: Compiler gfortran cannot compile programs. I checked that gfortran exists on my system M1 mac by running gfortran --version, which gives me NU Fortran (GCC) 12.2.0 Copyright (C) 2022 Free Software Foundation, Inc. This is free software; see the source for copying conditions. There is NO warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
Could you tell me how I shall proceed?

[andyfaff]

Are you on windows or Linux? I suggesting trying to compile a really simple test program (c.f. hello world) to see if the compiler actually works.

[ZhibingSun]

I am on M1 mac Ventura 13.6.6.

[andyfaff]

How did you install gfortran? I've found homebrew to work.

[ZhibingSun]

I installed it from here. https://github.com/fxcoudert/gfortran-for-macOS/releases. As I mentioned in my previous post, I ran gfortran --version and it works.

[andyfaff]

Running that command doesn't necessarily mean that the compiler works. Try compiling and running a test program.

[fancidev]

I followed Contributor quickstart guide and run command upto python3 dev.py test, but it gives me error meson.build:84:0: ERROR: Compiler gfortran cannot compile programs.

I use Mac M1 and had the same error message. I worked around it by executing the following before python dev.py build:

export LIBRARY_PATH="$LIBRARY_PATH:/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/usr/lib"

This works for me because after I upgraded XCode, some of the environment variables are not correctly set.

Btw, after the config error message, there should be somewhere that shows you a detailed error log. In that log it explains exactly why it “couldn’t compile Fortran program”. For me it was unable to “-lSystem”.

@ZhibingSun

[lucascolley]

These two lines from http://scipy.github.io/devdocs/building/index.html#system-level-dependencies were enough for me.