Hello,
I’ve opened an issue on SciPy recently and I was told that new features should be discussed here.
I implemented several special functions which I thought, you could be interested in:
The Dilogarithm function of a complex argument z:
The dilogarithm Li2(z) is already implemented in SciPy XSF, under spence.h, but I think you should consider my implementation because it’s very precise and fast.
This implementation is almost complete and is already optimized and written in C: https://github.com/Ruda975/LIBkorund/blob/main/dilog.h
It requires NaN and Inf checks but I can do that at your request.
Elliptic, Theta, Modular, Weierstrass functions of complex arguments:
My implementations are available here: https://github.com/Ruda975/elliptix
They are licensed under GPL3, but if you wanted them in SciPy, I would re-license to a compatible license.
I haven’t implemented these in C yet.
Anger-Weber functions of real order v and real argument x:
These require more testing and polishing before I finish the implementations, but testing looks promising.
Barnes-G function:
I’ve implemented this one as well, but as far as I know, it doesn’t have many uses in physics, or other scientific fields.
For precision testing, I’m using MpMath and Mathematica. They have these functions available in arbitrary precision.
For each of these, can you comment on whether “the method is applicable in many fields” (from SciPy Core Developer Guide — SciPy v1.18.0 Manual). There are so many special functions that SciPy should not aim to implement them all but rather focus on those that are most useful.
I’m not a physicist and haven’t used any of these functions myself in my research, but this is what I found:
Theta functions appear in solutions of heat flow PDEs.
Weierstrass P and related functions appear in the theory of elliptic curves and solutions of nonlinear differential equations.
The Dedekind Eta function occurs in bosonic string theory.
The Anger-Weber functions are solutions of the inhomogeneous Bessel equation, and they have applications in plasma physics.
The Barnes-G function doesn’t have many physics applications, so if we’re ranking the functions based on their usefulness in physics, you shouldn’t include this one.
Implementing numerical algorithms is actually one of my hobbies, I mainly write code for numerical integration, but thought special functions would be interesting as well. I decided to ask if anyone wanted these functions implemented, because I thought they would be useful to someone.
I did the theta/modular/Weierstrass functions, because I saw their domain coloring plots and thought they looked cool, but there were no algorithms to compute them very fast and to high precision. I am currently working on an interactive domain coloring plot and I want to include them with it.
I did the dilogarithm, because I saw other open source implementations, and thought I could do better. Also I wanted to practice writing optimized C code.
I implemented the Anger-Weber functions to test out my algorithm that evaluates infinite sums of Bessel-J functions, for example: http://dlmf.nist.gov/11.4.iv, and it worked quite well, so I decided to turn it into an actual implementation.
Thanks for clarifying. My preference would be to wait for a user request rather than add them because you have implementation but thank you putting them forward.
Alright, if new functions aren’t really needed, could I perhaps contribute by improving existing implementations?
I could add the Spence function to my dilogarithm implementation easily:
spence(z) = dilog(1 - z)
I’m confident it would make the spence function more precise and much faster.
It can be hard to know what to prioritise as maintainers when enhancements are suggested without a direct use-case as in this thread. But if you are potentially keen to participate in the process of prioritising different issues and contributions, that’s the kind of thing maintainers are typically quite happy to help with.
Last thing, if I’m coding an improved implementation, do i aim for maximum performance (using optimizations like FMA (fused multiply-add), which make the code way less readable, as in my implementation of the dilogarithm), or do I trade a bit of performance for readability/maintainability?
When making substantial changes, I would focus on ensuring correctness and keeping good readability, at least initially. Feel free to add code comments noting potential performance optimisations that could be added while trading off readability.
I just stumbled across this while looking for python implementations of the Anger-Weber functions. I work in ultrasonics and Scipy implementations of the Anger-Weber functions and the related Lommel functions (s, S, U, and V) would come in handy for certain beam propagation problems. I generally use mpmath when I need the latter, but using it on large arrays is pretty slow.
I am a user of scipy, and I only started to play around with it after I retired. It seems to have ‘everything’! My question: why should scipy only implement the 'most useful ones? Is it a question of maintaining them?
Good question! My phrasing was not the most precise so to quote the SciPy developer guide SciPy Core Developer Guide — SciPy v1.18.0 Manual “The method is applicable in many fields”. This is because every new feature comes with a maintenance burden so we need to consider how much bang for our buck (in this case maintainers time) each new feature gives. Hope this answers your question.
My Anger-Weber functions are written in Python and would need a C/C++ rewrite. I could try to implement the Lommel functions as well.
They currently work for real order v, and real argument x, but if someone finds use for it, I could likely do real order v and complex argument z.
If they don’t get added to scipy, I could open my own repo and share it with you.
The functions would accept numpy arrays. For the Anger-Weber functions, you can expect much higher speeds than mpmath, and similar, maybe faster speeds to scipy.special struve.