On Fri, Mar 27, 2020, at 6:46 PM, Aaron Meurer wrote:
> On Fri, Mar 27, 2020 at 4:58 PM Ondřej Čertík <[email protected]> wrote:
> >
> >
> >
> > On Fri, Mar 27, 2020, at 4:47 PM, Aaron Meurer wrote:
> > > On Fri, Mar 27, 2020 at 4:19 PM Ondřej Čertík <[email protected]> wrote:
> > > >
> > > > Hi,
> > > >
> > > > I am looking for students again this year for any related ideas:
> > > >
> > > > https://github.com/sympy/sympy/wiki/GSoC-2020-Ideas#sympy---fortran-code-generation-and-jit
> > > >
> > > > The ideas are from last year, and we made good progress on them last
> > > > year. There is plenty still to do. Here are some ideas in random order
> > > > and you are welcome to propose your own idea. If anybody is interested
> > > > in applying to any of those, let me know and I am happy to help
> > > > brainstorm more and review the application.
> > > >
> > > > Ondrej
> > > >
> > > > 1) Implement a module for efficient rational function approximation
> > > > with code generation to Fortran and C. It would accept a SymPy
> > > > expression as a function of one variable and it would output an
> > > > efficient rational function approximation. Here is an example how the
> > > > final function looks like for a modified Bessel function of half
> > > > integer order (this one was written by hand):
> > > >
> > > > https://github.com/certik/hfsolver/blob/b4c50c1979fb7e468b1852b144ba756f5a51788d/src/special_functions.f90#L371
> > >
> > > I have some code that computes the Chebyshev rational approximation of
> > > exp(-t) on [0, infinity) using the Remez algorithm here
> > > https://github.com/ergs/transmutagen/blob/master/transmutagen/cram.py.
> > > It uses SymPy and mpmath to compute an arbitrary number of digits of
> > > the approximation to an arbitrary degree. We also have a
> > > work-in-progress paper (never published) describing it in more detail
> > > which I can share if anyone is interested. Extending this to work with
> > > arbitrary functions would be very difficult, and would be a project
> > > all of its own (nothing to do with Fortran or code generation). One
> > > should also look at the chebfun MatLab package (which is open source).
> >
> > I implemented the basic chebfun algorithms in 1D and 2D (in Fortran), it's
> > very easy. The Remez algorithm is the most optimal, but my understanding is
> > that it is quite fragile sometimes. The Chebyshev approximation is up to a
> > factor of 2 worse in the error (I believe), but very robust / fast / easy
> > to implement. Here is Mathematica code that I used to produce that
> > approximation:
>
> Yes, it's very fragile. You can see in the code some of the things I
> had to do to make it work. I can't imagine what would be needed to
> make things work generally. It's especially bad once you increase the
> number of digits. The primary problem is that the algorithm involves
> root finding, one for finding the roots to a system of equations, and
> one for finding all the roots an in interval. If either one fails the
> algorithm fails. It is a good fit for SymPy because symbolics help a
> lot in the calculations, and also you want to use mpmath to compute
> more digits than machine precision. The nice thing about it though is
> that once you get the coefficients you are done. You can save them for
> later and you never need to run the algorithm again unless you want to
> increase the degree.
The reason it's good to have the Remez algorithm is that it gives you the best
possible approximation, and so it's helpful for comparisons. The Chebyshev
approach seems better suited for regular use. We should have both.
If any student is interested in any of these ideas, please let me know as soon
as possible.
Ondrej
--
You received this message because you are subscribed to the Google Groups
"sympy" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
To view this discussion on the web visit
https://groups.google.com/d/msgid/sympy/dcd674f0-266a-448e-b691-808f4ac2432c%40www.fastmail.com.
