So I was going to ask you about your thoughts how polynomials should be implemented in SymEngine.

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.
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.

So there could be a class Polynomial, that accepts a template (integer, rational, expression)

Would that do everything that we need?

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

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.

Aren't the generators just x, y, z, ...?

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`

?
@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.

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.

So "a" is also a generator?

```
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)

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)?

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

And in here, we have a polynomial of only one variable (x), and the coefficients are expressions.

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'.
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.
Now the question is how do we do ope rations on polynomial with different generators (or rings)

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?

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.

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).
Right! This is basically composing two rings (in mathematical terms)

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.

And the user just needs to explicitly convert the polynomial, which internally means rewriting the data structures, essentially padding them with zeros.

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.

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.

everything else should raise an exception. And the user needs to explicitly convert it first.

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.

Just for the sake of completeness, for

`a*x`

, we could also have a composite domain.
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?

```
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
```

The polynomial is still univariate, but it allows

`rational*a`

as coefficients.
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?

I haven't gone into the details of its implementation.

But here, it is

But here, it is

`p/q*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
```

It adds multiples of

`a`

into the domain.
Once you have

`p/q*a`

as a coefficient, you can do whatever you could with `p/q`

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]
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?

From what I can see from the code, it is a polynomial in

`a`

with rational coefficients
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'
Generic polynomials as I mentioned earlier.

Hm, that's confusing. Can you go one step further and recurse even more?

How is this used in sympy?

@jksuom might answer better.

It seems this is called a recursive polynomial representation.

If

`Polynomial`

is an `Expression`

then this is just a special case of `Polynomial<Expression>`

right?
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.

@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
```

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.

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.

I think the recursive form is a little different from sparse vs dense representation.

You have dense/sparse, then you have polynomials of one, two, three variables in them.

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.

Yes

So recursive says something about the domain being used, while sparse-dense is about the data structure being used.

So recursive says something about the domain being used, while sparse-dense is about the data structure being used.

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.

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.

That suffices for ring series too.

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?

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.
So when you build a ring series, the

`rs_mul`

etc. functions in SymEngine just accept just a polynomial as an argument?
In SymPy, the arguments of

`rs_mul`

can be any polynomial of any coefficient type?
And a generator with respect to which you cant the order to be calculated.

That's right.

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?
Both need to be in the same ring.

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?

The arguments need to be of the same ring (same generator and same domain)

Mathematically, that's the right way too. Because, on multiplying QQ with say RR, what should the result be?

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)
```

So we use the compiler to ensure that the coefficient type is the same.

mul_poly will do the same thing

And then inside

`mul_poly`

(which presumably `rs_mul`

calls) one checks that the generators are identical, otherwise raise an exception.
Does every polynomial store the list of generators and domain ?

Or do we have a

`ring`

?
In Sympy, ring does so.

Sorry for the confusion.

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?

So (gens + domain) make ring. Every polynomial just stores a

`ring`

I can create polynomials of a ring on the fly by doing

`R(expr)`

where R is a ring.
So in C++, Ring would be a class, and each polynomial would just store an instance of it?

That is a possible approach.

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

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.

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.

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.

Oh got it.

We store a pointer.

Just like this inside poly_mul:

`if (a.ring != b.ring) throw exception`

Or something like that.

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.

Right. We can just store RCP<const Ring> and compare the pointers.

When the last polynomial goes out of scope, it deallocates the Ring.

But that's not a major issue.

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.)

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.

Yes

I started documenting our discussion above.

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.

https://rubygems.org/gems/symengine

I think we will require the libsymengine too. Not sure how to include that with the gem yet.

@abinashmeher999 you should provide step by step instructions how to build it.

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

The main motivation behind having the gem on rubygems.org is that everything is setup only by

`gem install symengine`

@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
```

@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
I didn't have it installed before, let me know if I need to test anything on it.

But the instructions were such. That it is to be pushed after the

`.gem`

file has been generated
Anything I am missing?

@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?

@Sumith1896

Could you try

`$ irb`

an then in the interpreter

```
require 'symengine'
a = SymEngine::Integer.new(1)
```

@abinashmeher999, those instructions are for gems written in pure ruby. For native extensions, you have to compile them when installing

@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>
```

@isuruf Do we need to bringback

`extconf.rb`

and `mkmf`

then? If specified correctly in `gemspec`

gem should compile the extensions upon install.
@certik When we say list of generators, do we mean

`vec_basic`

?
@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?
@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`

@shivamvats yes, why not?

@shivamvats ah, you mean two different instances? No, then it would raise an exception.

@isuruf . That should do. Where do we keep the source files for

`libsymengine`

? I suppose we need to build that too.
I just don't know if the overhead of checking the generators is worth it.

Anyway, we can check that during implementation.

There are two approaches

- Build both
`libsymengine`

and ruby wrappers - 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?

1 seems easy to setup, but at the same time, won't there be conflicts if there happens to be an existing library?

There won't be any conflicts, if there's no bug in our CMake files

@isuruf What resources would you suggest for a good understanding of cmake, except the official documentation?

@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.

@abinashmeher999 I have no idea, I only refer official docs

@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)?

@certik By default yes. I had to modify the sparse polynomial (

`PolyElement`

) to allow rational exponents. Dense polynomials still follow the maths definition.
By default, do the exponents in sympy.polys allow negative integers?

dense or sparse?

(well I am interested in both, but we only use sparse for ring series)

@abinashmeher999, you can take a look at

`setup.py`

which needs to be done in ruby in a `extconf.rb`

file
I suggest to modify it.

to allow rational or to disable it?

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.

Personally I am fine with it.

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.

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"

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?

Is there any other issue with the ring series? Any other road block to overcome?

I have added the changes I made to polys and their reasons, to the docs.

We now need to extend it to other functions.

I see. That should be straightforward and other people can help out.

I have figured out how to make

`rs_series`

work with puiseux series. I will send a PR soon.
Yes

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`

.
One doubt, I always thought of

`a*x`

as in Ring<EX>(x). Any advantage of doing Ring<ZZ>(x, a)?
@Sumith1896 you can do both.

`Ring<ZZ>(x, a)`

will obviously be a lot faster.
@Sumith1896 How is

`EX`

represented internally in SymEngine?
@Sumith1896 I've extended the Example 4 (https://github.com/sympy/symengine/wiki/En-route-to-Polynomial) with this.

So

`Ring<ZZ>(x, a)`

will be a tuple of numbers with a pointer to the list `[x, a]`

@Sumith1896 yes,

`a*x`

is a multivariate polynomial.
@certik Sorry, I mean a dictionary. I was referring to the key of the dict.

@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.
I don't think it has a Ring. Does it store the generators inside each polynomial instance?

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.
Key represents the monomial type.

Looking here:

https://github.com/bluescarni/piranha/blob/master/tests/fateman1.hpp

https://github.com/bluescarni/piranha/blob/master/tests/fateman1_perf.cpp

It looks like the ring is just a

https://github.com/bluescarni/piranha/blob/master/tests/fateman1.hpp

https://github.com/bluescarni/piranha/blob/master/tests/fateman1_perf.cpp

It looks like the ring is just a

`piranha::integer`

.
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?

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.

That's an overhead isn't it?

The same problem we are facing now.

The same problem we are facing now.

Yes, I am looking into the source code to figure out how it works.

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.
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.

I thought it does hence, ours was the slower one

Only if the generators are different. In our case they are the same, so it avoids

`_normalise_polymul`

The relevant part is here: https://github.com/bluescarni/piranha/blob/03ae86014e1e1b90d7768ea48ffc012ebae4bb29/src/series.hpp#L543

It just calls

`==`

on `std::vector<symbol>`

to check if the generators are the same.
Just like you do.

I was searching for this piece of code in Piranha, good that we found it now

The generators comparison is here: https://github.com/bluescarni/piranha/blob/03ae86014e1e1b90d7768ea48ffc012ebae4bb29/src/symbol_set.hpp#L526

Hi. Question: For the symengine python bindings, is there a way to get the coefficients from a (linear) expression like .coeff(name) in SymPy?