MAINT:Translate the entirety of cephes into C++ #20390
Conversation
|
This is rebased on the branch for #20089 which is in turn rebased on the branch for #19954. If someone merges those, the diff on this one will shrink a little. cc @rlucas7, @tirthasheshpatel, @rgommers |
|
I’ll investigate more closely later, but it seems like the failures in stats on Linuxes 32 bit and aarch64 would be ok to resolve with tolerance bumps. Some of them are exact equality checks, though maybe there’s a reason for that. |
|
I’ll take a look this weekend.
|
Thanks. I appreciate it. |
|
I've completely excised the old C cephes and all tests pass locally. What remains is
The trickiest thing in this PR was probably getting the double double arithmetic to work. I opted to write it in a more modern style using operator overloading. Fortunately there was already direct testing for the double double arithmetic, because these tests failed at first despite all tests for functions passing. Phew. That is the largest diff I've ever worked on. |
735e90b to
be39bea
Compare
|
I've tried to rebase this on #20405, but things won't compile due to unsuccessful resolution of conflicts with #20403, which I merged today. @izaid, can you submit a PR to this branch fixing these conflicts? I could probably figure it out, but I'm running out of steam. The issue seems to be with using old cephes functions which no longer exist in the new ufunc infrastructure. Here's a sample of the errors I get when I try to compile. |
Think we fixed this now. |
Yes, fixed thanks to your help. I missed that this morning when editing out comments that were no longer relevant. |
|
One detail I am starting to see is that we are using |
|
I think it is a bit too late but If you also use the git's rename option then it would be possible to get a smaller diff and also a chance to still see the git history. |
Thanks @ilayn! Good point, but it actually does have a use in C++. It allows a function to appear in two different translation units, see the first answer here for instance https://stackoverflow.com/questions/29796264/is-there-still-a-use-for-inline |
I think making everything inline is necessary to avoid symbol clashes in a header only library. I was getting errors in some cases when leaving them out, in particular for functions which are used in ufuncs but also called in C code through wrappers. Things need to be header only so they can work on CUDA, but I've been trying to think of ways to mitigate increase in binary size. Update: @izaid beat me to it. |
I believe there are ways to have it instantiate only once. But if we put all the ufuncs in a single Python extension like |
I don't think it's too late actually. While I was working on this, I had to keep the old cephes around, and didn't think of the now obvious idea of backing up each cephes file before I translated it, doing git mv, and then replacing the cephes C file with its backup. Even before you mentioned it, I've thought about backing up the new C++ files, starting over, and doing git mv, and then immediately replacing each with its updated C++ version. It wouldn't be anywhere near as tedious as translating 80+ files from cephes, so I'd be happy to do it if it will make review easier. Before I do that though, I want to try to fix the failing tests on 32 bit Linux, but just haven't had time to debug yet. |
| #include "cephes/dd_idefs.h" | ||
|
|
||
|
|
||
| /* Computes fl(a+b) and err(a+b). */ |
There was a problem hiding this comment.
I don't think we should change anything here, but I'm curious if we really need volatile.
There was a problem hiding this comment.
I did some research to find why the volatile was there, and apparently it is to prevent aggressive optimizations from altering the results, which may have been more of a problem back then. Now it's only an issue if someone compiles with --ffastmath I think.
| /* multiply by power of 2 */ | ||
| x = std::ldexp(x, n); | ||
|
|
||
| return (x); |
There was a problem hiding this comment.
Don't know why there are brackets here, but also there are a lot of historical things that got carried over from cephes.
There was a problem hiding this comment.
It's just the coding style used in cephes. It was probably popular at the time. When I first started, I may have got rid of them, but after a while I just starting leaving all of the return statements unchanged.
99e61e4 to
820ed56
Compare
|
The failures look unrelated. I think this is good to merge. cc @rgommers |
There was a problem hiding this comment.
Thanks @steppi and @izaid, this looks like a large amount of very nice work. I went through the whole diff (that alone was a bit of work:) ), resolved comments that were addressed, and checked the build time impact (which seems fine).
There was a single unaddressed comments about static variables left - please have a look at that. Aside from that it seems good to go, so I'll hit the approve button here.
Related: when testing I found a build warning (with GCC 12) that came from a previous PR - it'd be nice to address that in a follow-up:
In file included from ../scipy/special/_special_ufuncs.cpp:15:
../scipy/special/special/gamma.h: In function 'std::complex<_Tp> special::gamma(std::complex<_Tp>) [with T = float]':
../scipy/special/special/gamma.h:19:44: warning: infinite recursion detected [-Winfinite-recursion]
19 | SPECFUN_HOST_DEVICE inline std::complex<T> gamma(std::complex<T> z) {
| ^~~~~
../scipy/special/special/gamma.h:20:17: note: recursive call
20 | return gamma(z);
| ~~~~~^~~
scipy/special/special/cephes/jv.h
Outdated
| double ak, bk, akl, bkl; | ||
| int sign, doa, dob, nflg, k, s, tk, tkp1, m; | ||
| static double u[8]; | ||
| static double ai, aip, bi, bip; |
There was a problem hiding this comment.
@steppi this is still unaddressed
Ah yes, in std::complex<double> gamma(std::complex<double> z)that lives in the @izaid's intention in scipy/scipy/special/special/gamma.h Lines 18 to 21 in b4f3037 is to write a templated version, (also in the special namespace), which can handle 32 bit float input in the ufuncs, by passing the inputs to the double version, but that the resolution works correctly seems like it might be based on undefined behavior. @izaid, what do you think we should do here? Just explicitly write a std::complex<float> gamma(std::complex<float> z)maybe? |
|
Thanks @rgommers! I appreciate that you went to the length of checking the full diff. Good catch that I missed the comment about the static variables. |
Yeah, that should just be calling our complex gamma function. I can fix it in another PR, it's unrelated to anything here. |
Update: I've rewritten this PR message since the original I wrote when this was a work-in-progress is no longer relevant.
This PR translates the entirety of Cephes into C++ capable of running on CUDA through conditional compilation. This PR is rebased on #20405 which is rebased on #20089 which is rebased on #19954, so though still gargantuan, the true diff is a little smaller than what is seen now.
Most of this was fairly straightforward, there was just a lot of it. One challenge was translating double-double arithmetic from
cephes/dd_real.cand friends (actually from theqdextended precision package) into C++. I chose to do this in a more modern style with operator overloading.The old cephes has been removed entirely. I wrote C wrappers for use in Cython and C files which have not been translated to C++.
special_c_wrappers.cppfor special functions, anddd_real_wrappers.cppfor double double arithmetic, which is used in_cunity.pxdand which has test code written in Cython.You have your work cut out for you @izaid. Thanks for agreeing to review this in your time between terms.