These are chat archives for symengine/symengine

21st
Aug 2015
Ondřej Čertík
@certik
Aug 21 2015 16:00
@shivamvats I am free now.
Shivam Vats
@shivamvats
Aug 21 2015 16:00
Hi @certik
Ondřej Čertík
@certik
Aug 21 2015 16:01
So I was going to ask you about your thoughts how polynomials should be implemented in SymEngine.
Shivam Vats
@shivamvats
Aug 21 2015 16:02
I think Sympy's polys is a very good choice. There are broadly three levels - domain, ring and finally polynomial.
In SymEngine we have started with the top most level, ie, polynomial.
We'll definitely be needing different domains (in the form of different implementations of hash_set)
So, we need to decide if we need a ring level.
Ondřej Čertík
@certik
Aug 21 2015 16:07
Right.
So for example, we can have 3 different implementations of hash_set, one for integers, another for rationals, and finally for general expressions.
We can probably use a template for that.
Shivam Vats
@shivamvats
Aug 21 2015 16:09
Yes
Ondřej Čertík
@certik
Aug 21 2015 16:09
So there could be a class Polynomial, that accepts a template (integer, rational, expression)
Would that do everything that we need?
Shivam Vats
@shivamvats
Aug 21 2015 16:10
That would work in a lot of cases, but not all.
To be able to create a totally generic polynomial, we need to have a ring layer of abstraction.
For example, we can create a ring of integers over some other ring.
It will also be needed if we want to implement a symbolic ring.
As in a symbolic ring, you want to differentiate between generators and coefficients.
Ondřej Čertík
@certik
Aug 21 2015 16:13
I am still confused by that.
Aren't the generators just x, y, z, ...?
Shivam Vats
@shivamvats
Aug 21 2015 16:16
Yes. Consider cos(a)*x**2 + sin(a)*x. Here, you might want to have x as a generator and sin(a) and cos(a) as symbolic coefficients.
When you have 4*x, you know that x is a generator (i.e the polynomial is in terms of x). What if you have a*x ?
Ondřej Čertík
@certik
Aug 21 2015 16:31
@shivamvats so if you have the 3 polynomial classes, as I envisioned above, then you still need to tell the class what the symbols are (e.g. x, y, z, ...). Then if you have a*x, you need to use the ExpressionPolynomial, and then you tell it that "x" is a symbol. Then it knows that "a" is the expression coefficient.
Shivam Vats
@shivamvats
Aug 21 2015 16:33
Well, that, effectively is specifying generators to the ring.
Currently, in ring series, we simply go on adding the symbolic elements we need, as symbol/generator to the ring. That way, we avoid using the ExpressionPolynomial.
Ondřej Čertík
@certik
Aug 21 2015 16:42
Ah, I was going to ask about it.
So "a" is also a generator?
Shivam Vats
@shivamvats
Aug 21 2015 16:43
In ring series, yes.
Ondřej Čertík
@certik
Aug 21 2015 16:43
So what kind of polynomial is cos(a)*x**2 + sin(a)*x ?
Shivam Vats
@shivamvats
Aug 21 2015 16:44
In [53]: sring(cos(a)*x**2 + sin(a)*x)
Out[53]: 
(Polynomial ring in x, cos(a), sin(a) over ZZ with lex order,
 x**2*cos(a) + x*sin(a))
So the generators are x, cos(a) and sin(a)
Ondřej Čertík
@certik
Aug 21 2015 16:46
Ah, so the cos(a) is treated as a variable in the polynomial, so it's a multivariate polynomial of 3 variables, x, cos(a) and sin(a)?
Shivam Vats
@shivamvats
Aug 21 2015 16:46
Contrast this with
In [57]: R, x = ring('x', EX)

In [58]: R(cos(a)*x**2 + sin(a)*x)
Out[58]: EX(cos(a))*x**2 + EX(sin(a))*x

In [59]: R
Out[59]: Polynomial ring in x over EX with lex order
Exactly
Ondřej Čertík
@certik
Aug 21 2015 16:46
And in here, we have a polynomial of only one variable (x), and the coefficients are expressions.
Shivam Vats
@shivamvats
Aug 21 2015 16:47
Right
Ondřej Čertík
@certik
Aug 21 2015 16:47
Let's say we have Polynomial<ZZ>, Polynomial<QQ> and Polynomial<EX> in C++.
Then the first case would be just Polynomial<ZZ> and you just need to tell it what the variables are, so you tell it that the variables are x, sin(a), cos(a).
In the second case, it would be Polynomial<EX> with the variable 'x'.
Shivam Vats
@shivamvats
Aug 21 2015 16:48
Yes
Ondřej Čertík
@certik
Aug 21 2015 16:49
So in the constructor of the polynomial, you need to tell it the list of symbols (generators), and then you just have to specify a list of coefficients and powers.
So it seems that allows everything.
Just like in sympy, we just need to allow the generator to be any expression, don't we?
It would be stored in some kind of a list, e.g. vec_basic in SymEngine.
Shivam Vats
@shivamvats
Aug 21 2015 16:50
Yes, that would be great.
Now the question is how do we do ope rations on polynomial with different generators (or rings)
Ondřej Čertík
@certik
Aug 21 2015 16:51
The polynomial datastructures then only reference the generators by numbers like 1, 2, 3, 4, i.e. the position of the exponent in the exponent list of a given term, i.e. x^2 y^3 z is just (2, 3, 1)
How does sympy do that?
Shivam Vats
@shivamvats
Aug 21 2015 16:52
Yes, since we want to use the sparse representation, we only need a tuple of exponents.
And we need to store the position of each gen at some place.
In Sympy, you can't do arithmetic operations involving different rings.
Ondřej Čertík
@certik
Aug 21 2015 16:53
I would be interested also how Piranha does that. Essentially, imagine you have a polynomial in 'x', let's say you have a term 'x^4', that is represented as a tuple of exponents (4). Then you multiply with a polynomial in 'y', 'z', let's say you have a term 'y^2 z^5', represented by a tuple (2, 5). So you need to create a new polynomial with generators (x, y, z) and then the first tuple becomes (4) -> (4, 0, 0) and the second tuple becomes (2, 5) -> (0, 2, 5) and then you just multiply the tuples (4, 0, 0) and (0, 2, 5), and the result is just adding the exponents, so you get (4, 2, 5).
Shivam Vats
@shivamvats
Aug 21 2015 16:54
Right! This is basically composing two rings (in mathematical terms)
Ondřej Čertík
@certik
Aug 21 2015 16:54
So if SymPy doesn't allow that, we don't need to allow that either. All we need is to just be able to convert a polynomial in 'x' to a new polynomial in 'x, y, z', as I described above.
Shivam Vats
@shivamvats
Aug 21 2015 16:54
Yes
Ondřej Čertík
@certik
Aug 21 2015 16:54
And the user just needs to explicitly convert the polynomial, which internally means rewriting the data structures, essentially padding them with zeros.
Shivam Vats
@shivamvats
Aug 21 2015 16:55
Right
Ondřej Čertík
@certik
Aug 21 2015 16:55
Then once you have the two polynomials in the same ring, as you say, or simply having the same Polynomial<ZZ> (or QQ or EX) and having the same generators, then you can operate on them.
Shivam Vats
@shivamvats
Aug 21 2015 16:56
You want all polynomials to have the same generators and domain, implying they all have the key tuple of the same size and the ith position means the same thing to both of them.
Ondřej Čertík
@certik
Aug 21 2015 16:56
exactly. And the same coefficient types.
Shivam Vats
@shivamvats
Aug 21 2015 16:56
Yes
Ondřej Čertík
@certik
Aug 21 2015 16:57
everything else should raise an exception. And the user needs to explicitly convert it first.
Shivam Vats
@shivamvats
Aug 21 2015 16:57
Yes. In ring series we do that internally.
Ondřej Čertík
@certik
Aug 21 2015 16:58
It's the same philosophy as with the coefficient types. If you know that your final expression needs QQ, you can start in QQ and never need to convert, or you can start in ZZ, which is faster, and then convert later (which takes some time). I can imagine each method to be faster for some application.
It's the same with the generators. If you know that you will need 'x, y, z' at the end, you can either start with 'x, y, z' (and use zero exponents) which is slower, but you don't need to convert later, or you start with just 'x', and then convert later.
Each method might be faster in some circumstances.
Shivam Vats
@shivamvats
Aug 21 2015 16:59
Right!
Just for the sake of completeness, for a*x, we could also have a composite domain.
Ondřej Čertík
@certik
Aug 21 2015 16:59
So the philosophy is to raise an exception and never do any implicit conversion. The user needs to call a method to do the (slow) conversion.
What is a composite domain?
Shivam Vats
@shivamvats
Aug 21 2015 17:00
In [60]: R, x = ring('x', QQ[a])

In [61]: R
Out[61]: Polynomial ring in x over QQ[a] with lex order

In [62]: a*x
Out[62]: a*x
Ondřej Čertík
@certik
Aug 21 2015 17:01
What exactly is QQ[a]?
Shivam Vats
@shivamvats
Aug 21 2015 17:01
The polynomial is still univariate, but it allows rational*a as coefficients.
Ondřej Čertík
@certik
Aug 21 2015 17:02
How exactly does it work, do you know?
if the coefficient is just QQ, then it is essentially two integers p/q.
if it is QQ[a], what kind of type exactly is it?
Shivam Vats
@shivamvats
Aug 21 2015 17:03
I haven't gone into the details of its implementation.
But here, it is p/q*a
Ondřej Čertík
@certik
Aug 21 2015 17:03
apparently it can also be p/q without a.
In [3]: var("a")
Out[3]: a

In [4]: R, x = ring('x', QQ[a])

In [5]: R
Out[5]: Polynomial ring in x over QQ[a] with lex order

In [6]: a*x
Out[6]: a*x

In [7]: x**2
Out[7]: x**2

In [8]: x**2+3
Out[8]: x**2 + 3

In [9]: x**2+3*a
Out[9]: x**2 + 3*a
Shivam Vats
@shivamvats
Aug 21 2015 17:04
Yes
It adds multiples of a into the domain.
Ondřej Čertík
@certik
Aug 21 2015 17:04
And powers too:
In [14]: x**2+3*a**4
Out[14]: x**2 + 3*a**4
Shivam Vats
@shivamvats
Aug 21 2015 17:05
Once you have p/q*a as a coefficient, you can do whatever you could with p/q
Ondřej Čertík
@certik
Aug 21 2015 17:06
With p/q, you can do any kind of arithmetic operation, since if you do (p/q)^10, it just becomes some other p/q.
If on the other hand you have p/q*a, it's not clear to me how it is represented internally and what operations are allowed. You can also do QQ[a, b]
Shivam Vats
@shivamvats
Aug 21 2015 17:07
Yes
Ondřej Čertík
@certik
Aug 21 2015 17:08
Isn't QQ[a, b] just a polynomial with rational coefficients in a and b?
I.e. aren't 'a' and 'b' just regular polynomial generators?
Shivam Vats
@shivamvats
Aug 21 2015 17:09
I guess so.
Isuru Fernando
@isuruf
Aug 21 2015 17:09
From what I can see from the code, it is a polynomial in a with rational coefficients
Shivam Vats
@shivamvats
Aug 21 2015 17:09
Only, they aren't treated as such
Isuru Fernando
@isuruf
Aug 21 2015 17:09
In [34]: x*a*(a+1)
Out[34]: (a**2 + a)*x
Ondřej Čertík
@certik
Aug 21 2015 17:10
So R, x = ring('x', QQ[a, b]) is some kind of recursive polynomial. The generator is 'x' and the coefficients are polynomials in QQ with generators 'a' and 'b'
Shivam Vats
@shivamvats
Aug 21 2015 17:10
Yes.
Generic polynomials as I mentioned earlier.
Ondřej Čertík
@certik
Aug 21 2015 17:10
Hm, that's confusing. Can you go one step further and recurse even more?
How is this used in sympy?
Shivam Vats
@shivamvats
Aug 21 2015 17:11
I haven't gone beyond one layer, but it should work
@jksuom might answer better.
Ondřej Čertík
@certik
Aug 21 2015 17:12
It seems this is called a recursive polynomial representation.
Isuru Fernando
@isuruf
Aug 21 2015 17:12
If Polynomial is an Expression then this is just a special case of Polynomial<Expression> right?
Ondřej Čertík
@certik
Aug 21 2015 17:12
If you have a polynomial 3xyz^2 + yz+..., you can treat it as a polynomial in 'x' with coefficients of a polynomial in 'y' with coefficients a polynomial in 'z'.
@isuruf I think so.
Shivam Vats
@shivamvats
Aug 21 2015 17:13
@isuruf Yes. But you need to use only a part of the symbolic terms, which makes it faster.
@certik yes
But I don't see much benefit in going beyond one layer for similar terms.
It makes sense if we have very different rings- say polynomial in x over a over matrix
In [81]: R, x = ring('x', QQ[a][b])
In [83]: R((a + 2)*(b + 1)*x + a*b)
Out[83]: ((a + 2)*b + a + 2)*x + a*b
Ondřej Čertík
@certik
Aug 21 2015 17:17
There are many ways how polynomials can be represented. This recursive way is just one of them. We can write a specialized data structures for that. Then you have dense and sparse representations (and various ways to implement those). Given our experience with SymPy, it seems we want to do sparse representation only and concentrate on ZZ, QQ and general Expression (EX) coefficients.
And in particular, do the implementation that Piranha is using.
It's just one of the many ways to represent polynomials, but it seems it is very fast and competitive for a lot of applications/benchmarks.
Shivam Vats
@shivamvats
Aug 21 2015 17:18
Yes, I mentioned it for the sake of completeness.
Ondřej Čertík
@certik
Aug 21 2015 17:18
I think when Mateusz implemented it, he didn't know which way is the best.
So he implemented all of them.
That's why we have dense and sparse, as well as this recursive representation in sympy.
Shivam Vats
@shivamvats
Aug 21 2015 17:19
I think the recursive form is a little different from sparse vs dense representation.
Ondřej Čertík
@certik
Aug 21 2015 17:20
You have dense/sparse, then you have polynomials of one, two, three variables in them.
Shivam Vats
@shivamvats
Aug 21 2015 17:20
Right
Ondřej Čertík
@certik
Aug 21 2015 17:20
Then you have recursive, where you can combine e.g. polynomial of one (or two, three) variables with coefficients of a polynomial in one (or two, three) variables.
And you can switch dense/sparse.
Shivam Vats
@shivamvats
Aug 21 2015 17:21
Yes
So recursive says something about the domain being used, while sparse-dense is about the data structure being used.
Ondřej Čertík
@certik
Aug 21 2015 17:22
It's like a 3D space --- one axis is dense/sparse, another axis is the number of generators, and another axis are the coefficient type (ZZ, QQ or a polynomial, which by itself has a 3D space of options).
Or as you said.
Shivam Vats
@shivamvats
Aug 21 2015 17:22
Exactly!
Ondřej Čertík
@certik
Aug 21 2015 17:23
and as it goes, you can probably figure out an application where each of these options will be the fastest.
However, I think our experience with sympy and piranha suggests that sparse and non-recursive (e.g. just using all the generators that we need, non-recursively) seems to be very very competitive.
Shivam Vats
@shivamvats
Aug 21 2015 17:24
Right.
That suffices for ring series too.
Ondřej Čertík
@certik
Aug 21 2015 17:24
Right.
So it seems the conclusion is to just have Polynomial<ZZ>, Polynomial<QQ>, Polynomial<EX> which are sparse polynomials (Piranha style) and then any number of generators, and those generators can be any SymEngine expressions.
@shivamvats, Is that consistent with your experience as well?
Shivam Vats
@shivamvats
Aug 21 2015 17:26
Yes.
It seems to be the best choice.
Btw, QQ[a][b] is isomorphic to QQ[a, b]. So not having a recursive representation doesn't really mean we are losing out on anything.
Ondřej Čertík
@certik
Aug 21 2015 17:27
So when you build a ring series, the rs_mul etc. functions in SymEngine just accept just a polynomial as an argument?
Shivam Vats
@shivamvats
Aug 21 2015 17:27
Yes
Ondřej Čertík
@certik
Aug 21 2015 17:28
In SymPy, the arguments of rs_mul can be any polynomial of any coefficient type?
Shivam Vats
@shivamvats
Aug 21 2015 17:28
And a generator with respect to which you cant the order to be calculated.
That's right.
Ondřej Čertík
@certik
Aug 21 2015 17:28
Right, and a generator.
Shivam Vats
@shivamvats
Aug 21 2015 17:29
*want
Ondřej Čertík
@certik
Aug 21 2015 17:29
So if we use templates for ZZ, QQ, EX, this would be a bit tricky I think.
What happens if you do rs_mul of polynomials in SymPy, let's say one is ZZ and the other one is QQ?
Shivam Vats
@shivamvats
Aug 21 2015 17:30
That can't be done
Both need to be in the same ring.
Ondřej Čertík
@certik
Aug 21 2015 17:31
In SymPy, you can't multiply polynomials unless the generators are the same. But can you multiply polynomials of the same generators, but one is ZZ and one is QQ, using the polys module? Or will they give you an exception?
Shivam Vats
@shivamvats
Aug 21 2015 17:31
The arguments need to be of the same ring (same generator and same domain)
Ondřej Čertík
@certik
Aug 21 2015 17:32
That's good and clean.
Shivam Vats
@shivamvats
Aug 21 2015 17:34
Mathematically, that's the right way too. Because, on multiplying QQ with say RR, what should the result be?
Ondřej Čertík
@certik
Aug 21 2015 17:34
In that case, there shouldn't be an issue, because we simply declare one template (coefficients) for the whole rs_mul function in SymEngine, and simply declare it something like:
template<C>
rs_mul(const Polynomial<C> &a, const Polynomial<C> &b, Polynomial<C> &c, Generator x)
Shivam Vats
@shivamvats
Aug 21 2015 17:34
Right
Ondřej Čertík
@certik
Aug 21 2015 17:35
So we use the compiler to ensure that the coefficient type is the same.
mul_poly will do the same thing
Shivam Vats
@shivamvats
Aug 21 2015 17:35
Nice
Ondřej Čertík
@certik
Aug 21 2015 17:35
And then inside mul_poly (which presumably rs_mul calls) one checks that the generators are identical, otherwise raise an exception.
Shivam Vats
@shivamvats
Aug 21 2015 17:36
Does every polynomial store the list of generators and domain ?
Or do we have a ring?
Ondřej Čertík
@certik
Aug 21 2015 17:37
In SymPy?
Shivam Vats
@shivamvats
Aug 21 2015 17:38
No, in Symengine
In Sympy, ring does so.
Ondřej Čertík
@certik
Aug 21 2015 17:38
In SymEngine we are still designing the polynomials.
Shivam Vats
@shivamvats
Aug 21 2015 17:38
I mean, how we want it to be.
Sorry for the confusion.
Ondřej Čertík
@certik
Aug 21 2015 17:39
For example, here: https://github.com/sympy/symengine/pull/597/files#diff-37b4e114aeb8a00554803e0174218a8cR344, it automatically padds the exponents with zero. Which, based on our discussion above, we don't want to do --- we want to raise an exception if the generators are not the same.
How exactly does this work in SymPy? I thought each polynomial stores the generators and a domain, doesn't it? Or what does the ring do?
Shivam Vats
@shivamvats
Aug 21 2015 17:41
So (gens + domain) make ring. Every polynomial just stores a ring
Ondřej Čertík
@certik
Aug 21 2015 17:41
Is ring an instance of a class?
Shivam Vats
@shivamvats
Aug 21 2015 17:41
I can create polynomials of a ring on the fly by doing R(expr) where R is a ring.
Ondřej Čertík
@certik
Aug 21 2015 17:42
So in C++, Ring would be a class, and each polynomial would just store an instance of it?
That is a possible approach.
Ondřej Čertík
@certik
Aug 21 2015 17:43
I think that's the way to go.
Because you want a very quick comparison if the ring is the same. Here we are comparing a vector of symbols: https://github.com/sympy/symengine/pull/597/files#diff-37b4e114aeb8a00554803e0174218a8cR348
Shivam Vats
@shivamvats
Aug 21 2015 17:43
This allows us to create polynomials of the same type on the fly. We dont need to specify the list of generators and domain, again.
Ondřej Čertík
@certik
Aug 21 2015 17:44
Each time a polynomial is multiplied. That's slow. There should be a specialized class handling the generators and the coefficient type, called a Ring.
Shivam Vats
@shivamvats
Aug 21 2015 17:44
So, we only compare the hash.
Ondřej Čertík
@certik
Aug 21 2015 17:44
If the hash is unique, which it typically isn't safe to assume.
But we can for example restrict multiplication of polynomials if they are constructed using the exact same instance of a Ring. Then you just compare a pointer.
Shivam Vats
@shivamvats
Aug 21 2015 17:45
How do we check that?
Oh got it.
We store a pointer.
Ondřej Čertík
@certik
Aug 21 2015 17:46
Just like this inside poly_mul:
if (a.ring != b.ring) throw exception
Or something like that.
Shivam Vats
@shivamvats
Aug 21 2015 17:47
I meant if we store copies of the same ring, we would need to compare the list of symbols. But if we store pointers, we just compare the addresses.
Ondřej Čertík
@certik
Aug 21 2015 17:48
Right. We can just store RCP<const Ring> and compare the pointers.
Shivam Vats
@shivamvats
Aug 21 2015 17:48
Yes
Ondřej Čertík
@certik
Aug 21 2015 17:48
When the last polynomial goes out of scope, it deallocates the Ring.
Shivam Vats
@shivamvats
Aug 21 2015 17:49
We might not always want that.
But that's not a major issue.
Ondřej Čertík
@certik
Aug 21 2015 17:49
We have to see what leads to the least amount of repetitive code. We might want to access the type of the coefficients, so perhaps we need something like RCP<const Ring<ZZ>>, so that you can get the type of the coefficients out of it inside the Polynomial and use it in defining the internal data structures.
(You can always store the Ring somewhere else using RCP, then it doesn't get deallocated obviously.)
Shivam Vats
@shivamvats
Aug 21 2015 17:51
Yes
Ondřej Čertík
@certik
Aug 21 2015 17:51
The generators can in principle be templates as well, but there I think a simpler code is simply to store them in a vector.
Actually, generators can't be templates, since you won't know until runtime what generators the user wants to use.
Since they can be any kind of expressions.
Shivam Vats
@shivamvats
Aug 21 2015 17:52
Moreover, there could be a lot of generators.
Yes
Ondřej Čertík
@certik
Aug 21 2015 17:52
Right.
I started documenting our discussion above.
Shivam Vats
@shivamvats
Aug 21 2015 17:56
Great!
Abinash Meher
@abinashmeher999
Aug 21 2015 17:57
Hi everyone, I just pushed the gem to RubyGems.org.
https://rubygems.org/gems/symengine
I think we will require the libsymengine too. Not sure how to include that with the gem yet.
Ondřej Čertík
@certik
Aug 21 2015 17:59
@abinashmeher999 you should provide step by step instructions how to build it.
Abinash Meher
@abinashmeher999
Aug 21 2015 17:59
Can anyone, preferably one who doesn't have the gem already, please confirm if you can install the gem?
gem install symengine
@certik
The main motivation behind having the gem on rubygems.org is that everything is setup only by gem install symengine
Sumith Kulal
@Sumith1896
Aug 21 2015 18:01
@abinashmeher999
sumith@sumith-Lenovo-Z50-70:~$ sudo gem install symengine
[sudo] password for sumith: 
Fetching: symengine-0.0.1.gem (100%)
Successfully installed symengine-0.0.1
Parsing documentation for symengine-0.0.1
Installing ri documentation for symengine-0.0.1
Done installing documentation for symengine after 0 seconds
1 gem installed
Isuru Fernando
@isuruf
Aug 21 2015 18:01
@abinashmeher999, yes, we should have it on rubygems.org, so that you can do install it easily
you are including the compiled symengine.so which will only work for the particular architecture it was compiled for
Abinash Meher
@abinashmeher999
Aug 21 2015 18:01
@Sumith1896 Thanks! :smile:
Sumith Kulal
@Sumith1896
Aug 21 2015 18:01
I didn't have it installed before, let me know if I need to test anything on it.
Shivam Vats
@shivamvats
Aug 21 2015 18:02
Nice work @abinashmeher999 :)
Abinash Meher
@abinashmeher999
Aug 21 2015 18:02
@isuruf I see.
But the instructions were such. That it is to be pushed after the .gem file has been generated
Anything I am missing?
Ondřej Čertík
@certik
Aug 21 2015 18:04
@shivamvats can you read https://github.com/sympy/symengine/wiki/En-route-to-Polynomial#design and let me know if I captured our discussion above?
Did I forget about anything important?
Abinash Meher
@abinashmeher999
Aug 21 2015 18:05

@Sumith1896
Could you try

$ irb

an then in the interpreter

require 'symengine'
a = SymEngine::Integer.new(1)
Isuru Fernando
@isuruf
Aug 21 2015 18:05
@abinashmeher999, those instructions are for gems written in pure ruby. For native extensions, you have to compile them when installing
Shivam Vats
@shivamvats
Aug 21 2015 18:06
@certik Sure
Sumith Kulal
@Sumith1896
Aug 21 2015 18:06
@abinashmeher999 Nope
sumith@sumith-Lenovo-Z50-70:~$  irb
irb(main):001:0> require 'symengine'
LoadError: libsymengine.so: cannot open shared object file: No such file or directory - /usr/local/lib/ruby/gems/2.2.0/gems/symengine-0.0.1/lib/symengine/symengine.so
    from /usr/local/lib/ruby/2.2.0/rubygems/core_ext/kernel_require.rb:54:in `require'
    from /usr/local/lib/ruby/2.2.0/rubygems/core_ext/kernel_require.rb:54:in `require'
    from /usr/local/lib/ruby/gems/2.2.0/gems/symengine-0.0.1/lib/symengine.rb:1:in `<top (required)>'
    from /usr/local/lib/ruby/2.2.0/rubygems/core_ext/kernel_require.rb:128:in `require'
    from /usr/local/lib/ruby/2.2.0/rubygems/core_ext/kernel_require.rb:128:in `rescue in require'
    from /usr/local/lib/ruby/2.2.0/rubygems/core_ext/kernel_require.rb:39:in `require'
    from (irb):1
    from /usr/local/bin/irb:11:in `<main>'
irb(main):002:0>
Abinash Meher
@abinashmeher999
Aug 21 2015 18:07
Thanks. Error, as expected.
@isuruf Do we need to bringback extconf.rb and mkmf then? If specified correctly in gemspec gem should compile the extensions upon install.
Sumith Kulal
@Sumith1896
Aug 21 2015 18:08
@certik When we say list of generators, do we mean vec_basic ?
Shivam Vats
@shivamvats
Aug 21 2015 18:09
@certik A question. Say I create two instances of the same ring(with same gens and domain), will rs_mul allow multiplication of polynomials over them?
Isuru Fernando
@isuruf
Aug 21 2015 18:10
@abinashmeher999, you need to have an extconf.rb and you can call cmake, make and then make install from it, instead of going through mkmf
Ondřej Čertík
@certik
Aug 21 2015 18:10
@Sumith1896 yes.
@shivamvats yes, why not?
Sumith Kulal
@Sumith1896
Aug 21 2015 18:11
Cool, this implementation is much cleaner.
Ondřej Čertík
@certik
Aug 21 2015 18:11
@shivamvats ah, you mean two different instances? No, then it would raise an exception.
Shivam Vats
@shivamvats
Aug 21 2015 18:11
Yes
Ondřej Čertík
@certik
Aug 21 2015 18:11
Well, we can design it either way.
Abinash Meher
@abinashmeher999
Aug 21 2015 18:11
@isuruf . That should do. Where do we keep the source files for libsymengine? I suppose we need to build that too.
Ondřej Čertík
@certik
Aug 21 2015 18:12
I just don't know if the overhead of checking the generators is worth it.
Shivam Vats
@shivamvats
Aug 21 2015 18:12
In most cases, it shouldn't be much.
Anyway, we can check that during implementation.
Isuru Fernando
@isuruf
Aug 21 2015 18:14

There are two approaches

  1. Build both libsymengine and ruby wrappers
  2. Build only the ruby wrappers and specify libsymengine as a library dependency.

1 is the approach for Python wrappers right now, but I like 2 better.
@certik, @abinashmeher999, thoughts?

Abinash Meher
@abinashmeher999
Aug 21 2015 18:16
1 seems easy to setup, but at the same time, won't there be conflicts if there happens to be an existing library?
Isuru Fernando
@isuruf
Aug 21 2015 18:16
Btw, for python, both methods are supported with #596
There won't be any conflicts, if there's no bug in our CMake files
Abinash Meher
@abinashmeher999
Aug 21 2015 18:22
@isuruf What resources would you suggest for a good understanding of cmake, except the official documentation?
Ondřej Čertík
@certik
Aug 21 2015 18:24
@shivamvats I have extended the notes: https://github.com/sympy/symengine/wiki/En-route-to-Polynomial#design, we also need to allow to specify the exponents type. So far we assume positive or negative integers (ZZ), but we also need to allow rationals QQ. See also https://github.com/sympy/symengine/wiki/En-route-to-Polynomial#ring-series which now lists exactly how Taylor, Laurent and Puiseux series should be implemented. Let me know if you agree.
Shivam Vats
@shivamvats
Aug 21 2015 18:25
I read it. You have summarised it pretty well.
Isuru Fernando
@isuruf
Aug 21 2015 18:25
@abinashmeher999 I have no idea, I only refer official docs
Shivam Vats
@shivamvats
Aug 21 2015 18:26
@certik Yes, that is the way to go
Ondřej Čertík
@certik
Aug 21 2015 18:26
@shivamvats thanks for your time, this was very useful for me to understand it. One last question --- in SymPy, the polynomials only allow integer exponents, is that right?
How did you get around to make it work with negative integers (for Laurent), as well as rational exponents (for Puiseux)?
Shivam Vats
@shivamvats
Aug 21 2015 18:28
@certik By default yes. I had to modify the sparse polynomial (PolyElement) to allow rational exponents. Dense polynomials still follow the maths definition.
Ondřej Čertík
@certik
Aug 21 2015 18:28
I see. I suggest to modify it.
By default, do the exponents in sympy.polys allow negative integers?
Shivam Vats
@shivamvats
Aug 21 2015 18:29
I have written about it in detail here
dense or sparse?
Ondřej Čertík
@certik
Aug 21 2015 18:29
sparse
(well I am interested in both, but we only use sparse for ring series)
Isuru Fernando
@isuruf
Aug 21 2015 18:30
@abinashmeher999, you can take a look at setup.py which needs to be done in ruby in a extconf.rb file
Shivam Vats
@shivamvats
Aug 21 2015 18:30

I suggest to modify it.

to allow rational or to disable it?

Abinash Meher
@abinashmeher999
Aug 21 2015 18:30
Thanks! I will have a look at it.
Ondřej Čertík
@certik
Aug 21 2015 18:30
I suggest to do your modification to allow to use both negative integers as well as positive and negative rationals in the sparse polys in sympy.
Unless you see some reason not to.
Shivam Vats
@shivamvats
Aug 21 2015 18:31
@certik I did that already to make ring series work :)
Ondřej Čertík
@certik
Aug 21 2015 18:31
Yes, but then you said you want to revert it.
Shivam Vats
@shivamvats
Aug 21 2015 18:32
jksuom suggested so.
Personally I am fine with it.
Ondřej Čertík
@certik
Aug 21 2015 18:32
What was his argument for reverting it?
Shivam Vats
@shivamvats
Aug 21 2015 18:33
He doesn't want to modify the internals of polynomials for a client module.
And then it also breaks the mathematical definition of a polynomial.
Ondřej Čertík
@certik
Aug 21 2015 18:33
I am totally fine with it --- yes, with rational or negative exponents it is not a polynomial in a strict mathematical sense, but it follows the same structure and it is very useful to be able to represent such "generalized polynomials"
Shivam Vats
@shivamvats
Aug 21 2015 18:34
Great!
Ondřej Čertík
@certik
Aug 21 2015 18:34
Not just for ring series, but for any application that requires fractional exponents. Piranha also allows that.
So as it is now in sympy, is it enabled, or is it reverted?
Shivam Vats
@shivamvats
Aug 21 2015 18:34
It is enabled.
Ondřej Čertík
@certik
Aug 21 2015 18:34
Ok, cool. Let's keep it that way.
Is there any other issue with the ring series? Any other road block to overcome?
Shivam Vats
@shivamvats
Aug 21 2015 18:35
I have added the changes I made to polys and their reasons, to the docs.
Ondřej Čertík
@certik
Aug 21 2015 18:35
Thanks.
Shivam Vats
@shivamvats
Aug 21 2015 18:35
I think all the major ones are done.
We now need to extend it to other functions.
Ondřej Čertík
@certik
Aug 21 2015 18:36
I see. That should be straightforward and other people can help out.
Shivam Vats
@shivamvats
Aug 21 2015 18:36
I have figured out how to make rs_series work with puiseux series. I will send a PR soon.
Yes
Ondřej Čertík
@certik
Aug 21 2015 18:37
Thanks, please do.
I read through http://doc.sagemath.org/html/en/tutorial/tour_polynomial.html, and it looks like we should also eventually implement CC (complex numbers) and any other field if needed.
I.e. if you do complex series expansion, you want to probably work using the complex numbers in SymEngine, which are of the rational form a/b + I*p/q.
Shivam Vats
@shivamvats
Aug 21 2015 18:38
I think, that will only need a complex domain.
Sumith Kulal
@Sumith1896
Aug 21 2015 18:39
One doubt, I always thought of a*xas in Ring<EX>(x). Any advantage of doing Ring<ZZ>(x, a)?
Ondřej Čertík
@certik
Aug 21 2015 18:40
@Sumith1896 you can do both. Ring<ZZ>(x, a) will obviously be a lot faster.
Shivam Vats
@shivamvats
Aug 21 2015 18:41
@Sumith1896 How is EX represented internally in SymEngine?
Sumith Kulal
@Sumith1896
Aug 21 2015 18:42
Value class of Basic is Expression
Ondřej Čertík
@certik
Aug 21 2015 18:42
@Sumith1896 I've extended the Example 4 (https://github.com/sympy/symengine/wiki/En-route-to-Polynomial) with this.
Shivam Vats
@shivamvats
Aug 21 2015 18:42
So Ring<ZZ>(x, a) will be a tuple of numbers with a pointer to the list [x, a]
Sumith Kulal
@Sumith1896
Aug 21 2015 18:43
Cool, a multivariate polynomial basically
Ondřej Čertík
@certik
Aug 21 2015 18:43
What kind of "tuple of numbers"?
@Sumith1896 yes, a*x is a multivariate polynomial.
Shivam Vats
@shivamvats
Aug 21 2015 18:44
@certik Sorry, I mean a dictionary. I was referring to the key of the dict.
Ondřej Čertík
@certik
Aug 21 2015 18:44
@Sumith1896 Ring<ZZ>(x, a) will internally just have a list [x, a] with the generators. The dictionary is inside a Polynomial, not a Ring.
Shivam Vats
@shivamvats
Aug 21 2015 18:45
Yes. Sorry for the confusion.
Ondřej Čertík
@certik
Aug 21 2015 18:45
Yes.
Sumith Kulal
@Sumith1896
Aug 21 2015 18:45
Thanks, things are clear now
Ondřej Čertík
@certik
Aug 21 2015 18:46
@Sumith1896 how does Piranha do that?
I don't think it has a Ring. Does it store the generators inside each polynomial instance?
Sumith Kulal
@Sumith1896
Aug 21 2015 18:47
I don't know how Piranha does it, it doesn't have a Ring though
It seems it accepts the Ring as the Cf template argument.
Sumith Kulal
@Sumith1896
Aug 21 2015 18:53
Oh yes
Key represents the monomial type.
It could also be piranha::rational.
The Key argument can be monomial<short>, but if you use monomial<rational>, then you get rational exponents, i.e. for a Puiseux series.
@bluescarni did I get it right?
Ondřej Čertík
@certik
Aug 21 2015 19:01
But looking here: https://github.com/bluescarni/piranha/blob/master/tests/fateman1.hpp, it looks like the coefficient types (and exponent types and internal data structure) are part of the templates, but the generators themselves seems to be treated at runtime somehow inside the polynomial class.

I.e.

    typedef polynomial<Cf,Key> p_type;
    p_type x("x"), y("y"), z("z"), t("t");
    auto f = x + y + z + t + 1;

the x is a polynomial with only one generator 'x'. The same for y, z and t. Then when you write f = x + y + z + t + 1, Piranha somehow converts the polynomials over different generators into a common base (by padding the zeros in the exponents) and adds them together.

Sumith Kulal
@Sumith1896
Aug 21 2015 19:05
That's an overhead isn't it?
The same problem we are facing now.
Ondřej Čertík
@certik
Aug 21 2015 19:09
Yes, I am looking into the source code to figure out how it works.
Ondřej Čertík
@certik
Aug 21 2015 19:15
I didn't find any operator+ on the polynomial type, so I don't understand how the line f = x + y + z + t + 1; can compile at all.
But polynomial seems to have m_symbol_set with the generators. So they seem to be saved inside the instance itself.
Ondřej Čertík
@certik
Aug 21 2015 20:29
Found it: http://bluescarni.github.io/piranha/doxygen/classpiranha_1_1series__operators.html#details, in particular: "in case two series arguments have different symbol sets, either one or both series will be copied in a new series in which the symbols have been merged, and the operation will be performed on those series instead;", so Piranha automatically adds the zero padding and copies both of the series so that they have the same generators.
Sumith Kulal
@Sumith1896
Aug 21 2015 20:30
So it doesn't avoid _normalise_polymul
I thought it does hence, ours was the slower one
Ondřej Čertík
@certik
Aug 21 2015 20:34
Only if the generators are different. In our case they are the same, so it avoids _normalise_polymul
Sumith Kulal
@Sumith1896
Aug 21 2015 20:35
Yes, the benchmark doesn't have different symbol set.
It just calls == on std::vector<symbol> to check if the generators are the same.
Just like you do.
Sumith Kulal
@Sumith1896
Aug 21 2015 20:41
I was searching for this piece of code in Piranha, good that we found it now
Michael Schubert
@mschubert
Aug 21 2015 23:52
Hi. Question: For the symengine python bindings, is there a way to get the coefficients from a (linear) expression like .coeff(name) in SymPy?