These are chat archives for symengine/symengine

17th
Nov 2015
Ralf Stephan
@rwst
Nov 17 2015 07:23
@certik At least this doesn't work directly in Piranha:
        p = (x - 1).invert();

function: pow
where: /usr/local/include/piranha/series.hpp, 2210
what: invalid argument for series exponentiation: negative integral value
Francesco Biscani
@bluescarni
Nov 17 2015 07:46
@rwst yes that function will work only for trivial series consisting of a single constant term
you can represent 1/(x-1) exactly if you use a divisor series, but I doubt that is what you are after
Ralf Stephan
@rwst
Nov 17 2015 07:47
@bluescarni any other piranha function providing this?
Francesco Biscani
@bluescarni
Nov 17 2015 07:48
you mean to provide a binomial-like expansion for 1/(x-1)?
or to represent that exactly?
Ralf Stephan
@rwst
Nov 17 2015 07:52
I think I want what you call binomial-like, and it's certainly not exact because series have finite precision (OK infinite is just polynomial).
Francesco Biscani
@bluescarni
Nov 17 2015 07:54
no there is no direct support for that type of expansions. I mean to write at one point a module for these expansions but so far I have been too busy with the core of the library
Ralf Stephan
@rwst
Nov 17 2015 07:54
NP @certik already gave a Python implementation, and GiNaC has such code too.
Francesco Biscani
@bluescarni
Nov 17 2015 07:54
actually an earlier version of Piranha from some years ago had these things implemented, but I was not really satisfied with the interface and haven't re-implemented them yet in the current version
ok nice :) If I recall correctly, in my old implementation most of this functionality could be implemented on top of a general hypergeometric series
but I was focused mostly on expansions for application to celestial mechanics
Ralf Stephan
@rwst
Nov 17 2015 07:58
Hypergeometrical not needed if the base is a poly. Simple recursion. One could give even a closed form for the nth coeff if the poly has degree 2.
Francesco Biscani
@bluescarni
Nov 17 2015 07:59
I think in general with a polynomial base one can use the multinomial theorem to have a closed form for the coefficients of the expansion?
Ralf Stephan
@rwst
Nov 17 2015 08:00
I'm only versed in the univarate case. There it is a general form of the Binet formula.
But since we need all coeffs up to N a closed form is not needed.
Francesco Biscani
@bluescarni
Nov 17 2015 08:03
it's something I wanted to investigate to see if it's possible to implement a faster exponentiation of polynomials, but I never looked into it much farther than speculation
But I guess the function would have to be generalized for negative exponents.
The Python implementation I posted above is using some kind of Newton iterations.
Ralf Stephan
@rwst
Nov 17 2015 14:45
I think the algorithm achieves the same as if you would have computed the individual coefficients isolated using the underlying recurrence, but by using polynomial arithmetic for chunks of coeffs simultaneously.
For example, for 1/(1-x) the intermediate results are:
1+x
1+x+x**2+x**3
1+x+x**2+x**3+x**4+x**5
1+x+x**2+x**3+x**4+x**5+x**6+x**7
Ralf Stephan
@rwst
Nov 17 2015 16:23
Ths works now with univariates, see https://github.com/rwst/symengine/blob/expansion4/symengine/tests/basic/test_series_expansion.cpp the last test case
Next, I'll probably add more functions.
Ondřej Čertík
@certik
Nov 17 2015 16:58
Nice.
So how would you handle 1/(1-(x^2+x+1))?
Ah I see, you already do.
I like that you are leveraging Piranha as much as you can. I think that's the way to go.
@rwst in Release mode, is it about as fast as your original raw Piranha benchmark?
Ondřej Čertík
@certik
Nov 17 2015 17:11
@rwst I reviewed what you did so far, and that's precisely what I had in mind. Functions like series_invert operate directly with Piranha's types and so should have zero overhead. Then you have _series which converts a Basic expression into a Piranha's polynomial. There the overhead also seems minimal --- looking at it, I don't see how it could be sped up more. Finally, the possible overhead is in the series method, that just calls _series and converts back to a vector of symengine numbers. That will have some overhead. This only happens at the very end, so this is independent of the complexity of the expression, it only depends on the number of terms. If this is an issue, we can figure out a different design. For now this is nice, since only the expansion.cpp file depends on Piranha, so the dependency is nicely isolated.
Let me know if you need any help with testing or benchmarking or developing.
Ralf Stephan
@rwst
Nov 17 2015 20:03
@certik The conversion overhead is negligible compared to, say one multiplication, so I'm not worried atm. The design is a copy of rs_seriesso the blame is not on me :smile: Complexity will come with multivar but I'd like to have univar first.
Ondřej Čertík
@certik
Nov 17 2015 20:11
rs_series is thanks to @shivamvats I think. I am glad you like it. :)