Working sympy into my package (or my package into sympy)?

[Matthew Samuel]

Hello sympy community,

I have a PhD in math and as a result of my research some unique knowledge that has allowed me to implement fast computations with Schubert polynomials. I compiled a python package in 2023 that started out as just console scripts, then it was suggested at a math conference that I try to integrate it into something like, say, Sage.

I did this, but Sage is monstrous and can only be used on a linux-based operating system, which is not ideal. Granted, the "market" for a Schubert polynomials package is quite limited (I say "market" because obviously no one would pay for this, it's a figure of speech), and the depth of computer algebra capabilities that are needed for this is quite low. No offense to sympy, but no one would argue that compared to Sage it is quite lightweight, takes seconds to install, etc.

Anyway my training is in math, not programming, and I do mainly program for a living, but I'm not a software engineer per se, and the subject of my job has little to do with computer algebra. Anyway I have created subclasses of Expr that deal with Schubert polynomials, double Schubert polynomials, quantum Schubert polynomials, quantum double Schubert polynomials, parabolic quantum (double) Schubert polynomials, and interplays between them.

I wanted to give a link to my github project for anyone who is interested at looking at it. I suspect that from your perspective it is pretty much a dumpster fire of bad practice, and any advice you may have is welcome.

https://github.com/matthematics/schubmult

I'm reaching out because I'm running into issues that should be pretty elementary if I were doing it right. For example, the best representation of linear combinations of Schubert polynomials is as a dict of coefficients of basis elements, but there are some times when they need to be left uncombined, but then you should be able to add to the uncombined terms and combine with the cominable ones. sympy does automatically do that, for example,

>>> DSx([3,4,1,2]) + DSx([4,1,3,2],"z")
DSx((3, 4, 1, 2), y) + DSx((4, 1, 3, 2), z)

This is sympy.Add, not a DoubleSchubertAlgebraElement, because the sets of coefficient variables, y and z, are different. This is what I want. When we do this, however,

>>> DSx([3,4,1,2]) + DSx([4,1,3,2],"z") + DSx([3,4,1,2])
2*DSx((3, 4, 1, 2), y) + DSx((4, 1, 3, 2), z)

This looks correct, however that 2*DSx((3, 4, 1, 2), y) is a sympy.Mul object, when what I want is for it to be internally represented as {(3,4,1,2): 2}, which it is not.

I feel like this is similar to the Poly class being a subclass of Basic instead of Expr, but I don't know if I want to be that rigid. And combining terms as Expr alone is inefficient. These can have hundreds or thousands of terms after multiplying, and the hashmap of keys is the most efficient way to combine coefficients.

I guess this is a question? What's the most elegant way to get that 2 into that dict? As well as an introduction and salutations.

Thanks,

Matt

[Richard Katz]

Neat question Matt.  I recall seeing this before - the question of whether to combine or not combine factors.

In fact, someone seems to have answered it in calculus - also saying that SymPy actually has a way to do this using the sympy subs() method.  If you have an expression involving x and you substitute a numeric value  for it sympy will combine the factors. But if instead you substitute a variable expression, sympy will leave the factors uncombined.

Amazingly, we can ask Google and Google has actually already catalogued this issue and its solution using subs().

 

Google writes:

This substitution can control whether factors are combined or left separate. For example:

• x = symbols('x')
• expr = (x + 1) * (x + 2)
• expr.subs(x, 3) will result in 12 (the factors are combined), while expr.subs(x, 'y') will result in (y + 1) * (y + 2) (the factors are left uncombined). 

 I do I remember this issue - of combining or not combining coming up a while back. HTH -- Rich

--
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 sympy+un...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/sympy/fec1e0d2-a271-4cea-8cd0-ab4c216ef721n%40googlegroups.com.

--

Richard Katz

pytho...@gmail.com

(415) 336-9559

https://www.linkedin.com/in/richardkatz/

rich...@berkeley.edu

[Oscar Benjamin]

On Fri, 25 Apr 2025 at 02:31, Matthew Samuel <matthe...@gmail.com> wrote:
>
> I'm reaching out because I'm running into issues that should be pretty elementary if I were doing it right. For example, the best representation of linear combinations of Schubert polynomials is as a dict of coefficients of basis elements, but there are some times when they need to be left uncombined, but then you should be able to add to the uncombined terms and combine with the cominable ones. sympy does automatically do that, for example,
>
> >>> DSx([3,4,1,2]) + DSx([4,1,3,2],"z")
> DSx((3, 4, 1, 2), y) + DSx((4, 1, 3, 2), z)
>
> This is sympy.Add, not a DoubleSchubertAlgebraElement, because the sets of coefficient variables, y and z, are different. This is what I want. When we do this, however,
>
> >>> DSx([3,4,1,2]) + DSx([4,1,3,2],"z") + DSx([3,4,1,2])
> 2*DSx((3, 4, 1, 2), y) + DSx((4, 1, 3, 2), z)
>
> This looks correct, however that 2*DSx((3, 4, 1, 2), y) is a sympy.Mul object, when what I want is for it to be internally represented as {(3,4,1,2): 2}, which it is not.
>
> I feel like this is similar to the Poly class being a subclass of Basic instead of Expr, but I don't know if I want to be that rigid. And combining terms as Expr alone is inefficient. These can have hundreds or thousands of terms after multiplying, and the hashmap of keys is the most efficient way to combine coefficients.

My first thought is that you should probably use something more like
Poly rather than Expr here. Combining as Expr alone will be
inefficient.

Does a typical calculation here use a closed finite set of DSx elements?

If so you could use e.g. a sparse matrix to represent the linear combination.

> I guess this is a question? What's the most elegant way to get that 2 into that dict?

If you are working with Expr then you probably just want to use
something like Expr.replace.

--
Oscar

[Matthew J. Samuel]

Hi Oscar,

Thanks for your reply. So the (double) Schubert polynomials just form a different basis for the regular multivariable polynomial ring, so in quantity they are infinite (indexed by permutations, which I created my own class for layered on top of the sympy Permutation class because it's not quite suitable for what I need) but like the ordinary monomials you only use finitely many at a time, and they're just regular polynomials. So I should be able to just add a normal algebraic combination of variables to it. What Poly would do, as far as I understand, is convert anything to a poly that you add to it. This can be done with the Schubert polynomials, but it's an extremely expensive operation, and you would only want to do it on purpose.

To be a little clearer, let me show an.example with my implementation.

from IPython.display import display
from sympy import init_session, init_printing
from schubmult import *
init_session(order='old',use_latex=True,pretty_print=True)
display(DSx([3,4,1,2]).expand())

$\displaystyle 2 x_{1} x_{2} y_{1} y_{2} + x_{1}^{2} x_{2}^{2} + y_{1}^{2} y_{2}^{2} + x_{1} x_{2} y_{1}^{2} + x_{1} x_{2} y_{2}^{2} + y_{1} y_{2} x_{1}^{2} + y_{1} y_{2} x_{2}^{2} - x_{1} y_{1} x_{2}^{2} - x_{1} y_{1} y_{2}^{2} - x_{1} y_{2} x_{2}^{2} - x_{1} y_{2} y_{1}^{2} - x_{2} y_{1} x_{1}^{2} - x_{2} y_{1} y_{2}^{2} - x_{2} y_{2} x_{1}^{2} - x_{2} y_{2} y_{1}^{2}$

This is substitutible, so if I plug in "z" rather than the default y variable I get

$\displaystyle 2 x_{1} x_{2} z_{1} z_{2} + x_{1}^{2} x_{2}^{2} + z_{1}^{2} z_{2}^{2} + x_{1} x_{2} z_{1}^{2} + x_{1} x_{2} z_{2}^{2} + z_{1} z_{2} x_{1}^{2} + z_{1} z_{2} x_{2}^{2} - x_{1} z_{1} x_{2}^{2} - x_{1} z_{1} z_{2}^{2} - x_{1} z_{2} x_{2}^{2} - x_{1} z_{2} z_{1}^{2} - x_{2} z_{1} x_{1}^{2} - x_{2} z_{1} z_{2}^{2} - x_{2} z_{2} x_{1}^{2} - x_{2} z_{2} z_{1}^{2}$

The polynomial ring is in the indexed x variables, and another set of indexed symbols that are considered constant and part of the "double" Schubert polynomials (DSx), which are arbitrary. But only two sets of variables are required to express everything, so if I include both y and z, I don't have a basis anymore and it's difficult to test equality, but it's a valid operation, which is why it makes sense to leave it as an uncombined sympy.Add.

The coefficients of the double Schubert polynomials can also themselves be polynomials, and I want to leave those as Expr because sometimes I want to express them in a certain way that it's not obvious they can be expressed in if you expand them into monomials, and the expression is not in terms of a basis, and it's extremely expensive to compute. For example,

display(DSx([3,4,1,2],"y")*DSx([4,1,3,2],"z"))

$\displaystyle \left(y_{3} - z_{3}\right) \left(\left(y_{3} - z_{1}\right) \left(\left(y_{1} - z_{1}\right) \left(y_{4} - z_{1}\right) + \left(y_{3} - z_{2}\right) \left(y_{1} + y_{4} - z_{1} - z_{2}\right)\right) - \left(y_{1} - z_{1}\right) \left(y_{3} - z_{1}\right) \left(y_{4} - z_{1}\right)\right) \mathfrak{S}_{\left[ 3, \  4, \  1, \  2\right]}{\left(x; y \right)} + \left(y_{3} - z_{3}\right) \left(\left(y_{3} - z_{1}\right) \left(y_{3} + y_{4} - z_{1} - z_{2}\right) - \left(y_{3} - z_{1}\right) \left(y_{4} - z_{1}\right)\right) \mathfrak{S}_{\left[ 3, \  4, \  2, \  1\right]}{\left(x; y \right)} + \left(y_{3} - z_{3}\right) \left(\left(y_{3} - z_{1}\right) \left(y_{1} + y_{3} - z_{1} - z_{2}\right) - \left(y_{1} - z_{1}\right) \left(y_{3} - z_{1}\right)\right) \mathfrak{S}_{\left[ 3, \  5, \  1, \  2, \  4\right]}{\left(x; y \right)} + \left(\left(y_{3} - z_{1}\right) \left(\left(y_{1} - z_{1}\right) \left(y_{4} - z_{1}\right) + \left(y_{3} - z_{2}\right) \left(y_{1} + y_{4} - z_{1} - z_{2}\right)\right) + \left(y_{4} - z_{3}\right) \left(\left(y_{1} - z_{1}\right) \left(y_{4} - z_{1}\right) + \left(y_{3} - z_{2}\right) \left(y_{1} + y_{4} - z_{1} - z_{2}\right)\right) - \left(y_{1} - z_{1}\right) \left(y_{3} - z_{1}\right) \left(y_{4} - z_{1}\right)\right) \mathfrak{S}_{\left[ 4, \  3, \  1, \  2\right]}{\left(x; y \right)} + \left(\left(y_{3} - z_{1}\right) \left(y_{3} + y_{4} - z_{1} - z_{2}\right) + \left(y_{4} - z_{3}\right) \left(y_{3} + y_{4} - z_{1} - z_{2}\right) - \left(y_{3} - z_{1}\right) \left(y_{4} - z_{1}\right)\right) \mathfrak{S}_{\left[ 4, \  3, \  2, \  1\right]}{\left(x; y \right)} + \left(\left(y_{3} - z_{1}\right) \left(y_{1} + y_{3} - z_{1} - z_{2}\right) + \left(y_{4} - z_{3}\right) \left(y_{3} + y_{4} - z_{1} - z_{2}\right) - \left(y_{1} - z_{1}\right) \left(y_{3} - z_{1}\right)\right) \mathfrak{S}_{\left[ 4, \  5, \  1, \  2, \  3\right]}{\left(x; y \right)} + \left(\left(y_{1} - z_{1}\right) \left(y_{4} - z_{1}\right) + \left(y_{3} - z_{1}\right) \left(y_{1} + y_{3} - z_{1} - z_{2}\right) + \left(y_{3} - z_{2}\right) \left(y_{1} + y_{4} - z_{1} - z_{2}\right) + \left(y_{5} - z_{3}\right) \left(y_{1} + y_{3} - z_{1} - z_{2}\right) - \left(y_{1} - z_{1}\right) \left(y_{3} - z_{1}\right)\right) \mathfrak{S}_{\left[ 5, \  3, \  1, \  2, \  4\right]}{\left(x; y \right)} + \left(y_{3} + y_{4} + y_{5} - z_{1} - z_{2} - z_{3}\right) \mathfrak{S}_{\left[ 5, \  3, \  2, \  1, \  4\right]}{\left(x; y \right)} + \left(y_{3} + y_{4} + y_{5} - z_{1} - z_{2} - z_{3}\right) \mathfrak{S}_{\left[ 5, \  4, \  1, \  2, \  3\right]}{\left(x; y \right)} + \left(y_{1} + y_{3} - z_{1} - z_{2}\right) \mathfrak{S}_{\left[ 6, \  3, \  1, \  2, \  4, \  5\right]}{\left(x; y \right)} + \mathfrak{S}_{\left[ 6, \  3, \  2, \  1, \  4, \  5\right]}{\left(x; y \right)} + \mathfrak{S}_{\left[ 6, \  4, \  1, \  2, \  3, \  5\right]}{\left(x; y \right)}$

This is a tiny example and it's already degree 6 and very long.

So anyway, Poly style makes sense (these are just polynomials after all), but I don't want to be strict about the generating set (there are supposed to be infinitely many variables, I just fix it at 100 because if you're multiplying Schubert polynomials in more than 20 variables you had better have a supercomputer or a lot of time on your hands).

I want to avoid automatically coercing things into the basis (which is what original versions did) because it's expensive, and for example here's what DSx([3,4,1,2],"z") looks like in the y basis:

$\displaystyle \left(2 y_{1} y_{2} z_{1} z_{2} + y_{1}^{2} y_{2}^{2} + z_{1}^{2} z_{2}^{2} + y_{1} y_{2} z_{1}^{2} + y_{1} y_{2} z_{2}^{2} + z_{1} z_{2} y_{1}^{2} + z_{1} z_{2} y_{2}^{2} - y_{1} z_{1} y_{2}^{2} - y_{1} z_{1} z_{2}^{2} - y_{1} z_{2} y_{2}^{2} - y_{1} z_{2} z_{1}^{2} - y_{2} z_{1} y_{1}^{2} - y_{2} z_{1} z_{2}^{2} - y_{2} z_{2} y_{1}^{2} - y_{2} z_{2} z_{1}^{2}\right) + \left(z_{1} \left(- y_{1} y_{2} - y_{1} y_{3}\right) + z_{2} \left(- y_{1} y_{2} - y_{1} y_{3}\right) + z_{1} z_{2} \left(y_{2} + y_{3}\right) + 2 y_{1} z_{1} z_{2} + y_{1} z_{1}^{2} + y_{1} z_{2}^{2} + y_{2} y_{1}^{2} + y_{3} y_{1}^{2} - z_{1} y_{1}^{2} - z_{1} z_{2}^{2} - z_{2} y_{1}^{2} - z_{2} z_{1}^{2}\right) \mathfrak{S}_{\left[ 1, \  3, \  2\right]}{\left(x; y \right)} + \left(z_{1} z_{2} - y_{1} z_{1} - y_{1} z_{2} + y_{1}^{2}\right) \mathfrak{S}_{\left[ 1, \  4, \  2, \  3\right]}{\left(x; y \right)} + \left(y_{1} \left(y_{1} + y_{2}\right) + y_{3} \left(y_{1} + y_{2}\right) + z_{1} z_{2} + z_{1} \left(- y_{1} - y_{2}\right) + z_{2} \left(- y_{1} - y_{2}\right) - y_{3} z_{1} - y_{3} z_{2} + z_{1}^{2} + z_{2}^{2} - y_{1}^{2}\right) \mathfrak{S}_{\left[ 2, \  3, \  1\right]}{\left(x; y \right)} + \left(y_{1} + y_{2} - z_{1} - z_{2}\right) \mathfrak{S}_{\left[ 2, \  4, \  1, \  3\right]}{\left(x; y \right)} + \mathfrak{S}_{\left[ 3, \  4, \  1, \  2\right]}{\left(x; y \right)}$

This is really not an efficient way to express what could be a single term.

The internal representation is

class BasisSchubertAlgebraElement(Expr):

    _op_priority = 1e200

    _kind = NumberKind

    is_commutative = False

    # precedence = 40

    do_parallel = False

    def __new__(cls, _dict, basis):

        obj = Expr.__new__(cls)

        obj._dict = {k: sympify(v) for k, v in _dict.items() if expand(v) != S.Zero}

        if len(obj._dict.keys()) == 1 and next(iter(obj._dict.values())) == S.One:

            obj.precedence = 1000

        else:

            obj.precedence = 40

        obj._basis = basis

        return obj

    @property

    def args(self):

        return (sympy.Dict(self._dict), self._basis)

This is the base class for all the different types of Schubert polynomials. I threw the "is_commutative = False" in there for printing purposes, they are commutative. Anyway I avoid storing the dict and basis directly as args because it's actually not at all worth it to attempt caching non-atomic expressions because of the polynomial coefficients. I also use symengine for the internal algebraic manipulation of the coefficients, and converting back and forth between sympy and symengine causes significant overhead. They are only generated when asked for.

Anyway any advice you may have is welcome.

Thanks,

Matt

--
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 sympy+un...@googlegroups.com.

To view this discussion visit https://groups.google.com/d/msgid/sympy/CAHVvXxSYsJ_YNEn1O38DcJbFx4g7ojtmg%3D_xp6qHGhCbWyeGPQ%40mail.gmail.com.

[Oscar Benjamin]

Hi Matt,

When I say "something more like Poly rather than Expr" I don't necessarily mean that you should use Poly. The difference between these is that:

- Expr is a general representation of expressions as trees.

- Poly is a specialised data structure that represents polynomials in a canonical form over some domain.

Generally if you are doing any significant computational work it is better to use a specialised data structure and/or canonical representation rather than expression trees. I don't understand exactly what sort of calculations you want to do with the Schubert polynomials but I imagine that if I was doing it I would make a special (non-Expr) data structure for the main calculations and then separately have an .as_expr() method (like Poly does) for converting to Expr if there are situations where that is useful.

It is possible that some existing things in SymPy would provide the data structure that you need but I don't know because I'm a bit unclear what it is you need to do with it. One possibility for example is that it could make sense to wrap your Schubert polynomials up as a "domain" so that you could have a Poly whose coefficients are Schubert polynomials. Another possibility is that you could choose a finite basis for a calculation and use matrices to represent e.g. linear combinations of Schubert polynomials.

Oscar

To view this discussion visit https://groups.google.com/d/msgid/sympy/CAF%2B7YfG91iSNELEQmDctbxZdA_3%3DP%2Bgrxifb3njW7GOfBSPV4A%40mail.gmail.com.

[Matthew J. Samuel]

On Sat, Apr 26, 2025 at 8:45 AM Oscar Benjamin <oscar.j....@gmail.com> wrote:

Hi Matt,

When I say "something more like Poly rather than Expr" I don't necessarily mean that you should use Poly. The difference between these is that:

- Expr is a general representation of expressions as trees.

- Poly is a specialised data structure that represents polynomials in a canonical form over some domain.

My BasisSchubertAlgebraElement class is just this, with the canonical form being as a dict with Permutation keys representing the basis elements, with the values being the coefficients. Except that it's also a subclass of Expr. I'm trying to have my cake and eat it too I guess.

My original package was simply to compute coefficients of products of double Schubert polynomials. The output was the coefficients line by line. My advisor wrote the original C program to do this that can actually be installed with apt-get in Ubuntu, but it was only for single Schubert polynomials. My original program mimics his, but was in python and used symengine. It was to multiply pairs of double Schubert polynomials in different sets of coefficient variables. Here is my paper on the topic. You don't need to look at it to understand this, I'm just providing the link on the off chance that you are interested.

https://arxiv.org/abs/2401.11060

This is my original script where symengine was sufficient:

 Then I wanted to test Conjecture 1 in the paper computationally, which took me months to figure out how to do, but it involved some small use of sympy, and ended up with the output re-expressed like this:

It expresses the result as a polynomial in y_i-z_j with nonnegative integer coefficients. In the worst case this involves integer programming, and actually this is one of the worst cases, but it's very fast because it's very small. I have tricks that can do it without integer programming maybe 75% of the time, but for the script it falls back to integer programming if none of those tricks work.

The primary goal of this is personal research, but I'm happy to share it with others who could also make use of it, which is why I'm publishing it as I go. There were research problems I avoided approaching because they seemed difficult to do that are now easier to approach since I can manipulate linear combinations of Schubert polynomials with a computer algebra system.

For example, I came across a paper that said that a Schubert polynomial can be factored as a standard elementary monomial if and only if the permutation avoids 312 and 1432. This is highly useful for my research in that if this is the case, I don't have to do any work to represent the product positively, it just automatically comes out positive with my method. But my result is positive more generally if it can be expressed as an elementary monomial that is not necessarily standard. I was able to use this symbolic library to figure out exactly when that happens. It is if the permutations avoids 1432, 1423, 3142, and 4132.

Then I wanted to know more generally whether it could just be factored with at least one of the factors being an elementary symmetric polynomial. This can be tested using sympy.simplify(). This permutation I keep using as an example, 4132, can be factored in this way but not completely.

That factor on the right makes my life easier; it is automatically positive, so even though this is a fundamentally non-positive case, I only have to worry about the first factor. This has practical value. If I have to resort to integer programming, it reduces the degree by 1, which reduces the dimension by a factor of 2, which is great. This is still new stuff and I have not worked it into the base library as an optimization yet. I never would have figured this out though, or it would have been much more difficult, without this symbolic library.

What am I trying to do? That's a good question, and I'm not 100% sure. The ultimate goal is basically unachievable, which is to prove that every coefficient is positive with a formula. But for example beyond that, this library does things that I'm pretty sure no one else even knows how to approach. There are also quantum double Schubert polynomials, and for example the same product I've been using as an example in quantum form can be expressed positively with one of my newer scripts, only developed last year:

With my method, which turned out to be sort of a lucky accident, quantum products can be computed in the same way, which is a breakthrough. As far as I know, while in theory people know how to compute quantum products, it was computationally infeasible except in the simplest cases. While similar to the base case, many of the most fundamental hooks into the theory of multiplying Schubert polynomials break in quantum. The attempt to find formulas was essentially abandoned (but not entirely of course, the field still exists and people are working on related aspects, just not really finding positive formulas). Quantum, in particular, is relevant in physics. My basic method of multiplying them just lifts and shifts, and I have made a lot of progress beyond that and even found some positive formulas. I have implemented quite a few conjectural results in this library that will be very hard to prove, especially since I have a very hard time writing papers in the first place (which is why I am not an academic, even though I managed to achieve the PhD with a dissertation that required countless revisions, and then this paper that took me years to write that the referee said was poorly written but worth fixing because the results were great).

I have quantum Schubert polynomials implemented as an Expr class as well.

Anyway there is huge potential in this, and I have set a basic goal of making my code 1) clean, 2) usable by others, 3) efficient and using good practices. Which is how I came to this simple question that I asked. The real question is "How can I write this Schubert polynomial class so that it seems like it could've been written by an expert in computer algebra?" Which seems to be the question you are trying to answer.

Regarding the specific suggestions you are offering,

 

Generally if you are doing any significant computational work it is better to use a specialised data structure and/or canonical representation rather than expression trees. I don't 

understand exactly what sort of calculations you want to do with the Schubert polynomials but I imagine that if I was doing it I would make a special (non-Expr) data structure for the main calculations and then separately have an .as_expr() method (like Poly does) for converting to Expr if there are situations where that is useful.

I accept this. I have taken the route of canonical expression already because it's the only way my method works. You can think of the Schubert polynomials as basically just really weird monomials. I do not know how familiar you are with abstract algebra, but these behave exactly like monomials that are products of variables and this class should in principle be identical to Poly, just with really weird monomials that have a special multiplication method that is extremely complicated and expensive. While this is the case, it can be a black box to the class. It might make sense to make it a subclass of Poly, because all that needs to change is 1) the keys in the dict, 2) the way of multiplying the basis elements. In a perfect world it would just be plug and play, and actually this is how it works in Sage. It was just a few lines of code to get this into Sage because for mathematicians this is an extremely common situation. Like I said, though, Sage is monstrous, written primarily by amateurs (no offense to the wonderful coders, it's just that it's difficult and error prone to install, and even I who grew up immersed in programming and can make almost anything work had a hard time installing it). However, Poly limits to a finite generating set, which actually is OK in practice and would work perfectly well for this, I just dislike the setup and the expressions that explicitly display the generating set doesn't look great to me. However, thinking about it, I could get over this.

So I looked at the source code for Poly and there.s a lot going on. But as one of the primary coders of sympy this is a great question to ask: is the Poly class set up in such a way that I can just override the keys of the monomials and the multiplication method? We would be good to go then, except for the fact that I want to keep the coefficients of the monomials (which are themselves free polynomials) as Expr because I want the freedom of expressing them in different ways, and I also found the internal multiplication methods of sympy to be disappointingly slow for multiplying and adding these coefficients, which is why I use symengine. The bare bones symengine wrapper does a much better job for the simple things I need to do with the coefficients.

It is possible that some existing things in SymPy would provide the data structure that you need but I don't know because I'm a bit unclear what it is you need to do with it. One possibility for example is that it could make sense to wrap your Schubert polynomials up as a "domain" so that you could have a Poly whose coefficients are Schubert polynomials. Another possibility is that you could choose a finite basis for a calculation and use matrices to represent e.g. linear combinations of Schubert polynomials.

Sorry for writing a novel on this but your answer helped me think of what the question could be. I will look into its feasibility myself and again any advice you may have is welcome. 

Thanks,

Matt

To view this discussion visit https://groups.google.com/d/msgid/sympy/CAHVvXxTThndRoBbUZ7hi_16RiAKKc18BoaRY%2Bte7Merjr0GVYg%40mail.gmail.com.

[Matthew J. Samuel]

Hmm actually the Domain framework looks like pretty much what I want, or at least something related will be, except for the part about coefficients being Exprs. Anyway I will not rapid fire emails, just preempting you saying that this what I should've been doing all along.

[Matthew J. Samuel]

OK. 

ExpressionDomain seems to be basically what I want for coefficients.

The only other issue is I need a basis PolyRing class, where instead of monomials we have Schubert polynomials (double, quantum, quantum double, etc.) with a special multiplication method.

In sage, the latter functionality was achieved with CombinatorialFreeModule that has product_on_basis. Then everything else is handled internally.

Best case, sympy has something like this. If not, we can fix this.

Then we will worry about speed.

[Aaron Meurer]

I'm not sure if it supports what you need yet, but there is some code for free modules in sympy/polys/agca/

Aaron Meurer

To view this discussion visit https://groups.google.com/d/msgid/sympy/CAF%2B7YfGi%3DJg9ChjJ1KqKh76_qT12S6LxXg5NtRH1QvheWwiX7w%40mail.gmail.com.

[Oscar Benjamin]

Note that there are two versions of ExpressionDomain: EX and EXRAW. The EXRAW domain is exactly the same as Expr. The EX domain also uses Expr but calls cancel after each operation to put it into a canonical rational form:

In [4]: x1 = Poly(x, y)

In [5]: x2 = Poly(x, y, domain=EX)

In [6]: x3 = Poly(x, y, domain=EXRAW)

In [7]: (1 + x1)**2
Out[7]: Poly(x**2 + 2*x + 1, y, domain='ZZ[x]')

In [8]: (1 + x2)**2
Out[8]: Poly(x**2 + 2*x + 1, y, domain='EX')

In [9]: (1 + x3)**2
Out[9]: Poly((x + 1)**2, y, domain='EXRAW')

To view this discussion visit https://groups.google.com/d/msgid/sympy/CAF%2B7YfGi%3DJg9ChjJ1KqKh76_qT12S6LxXg5NtRH1QvheWwiX7w%40mail.gmail.com.

[Matthew J. Samuel]

This is great, thanks to both of you.

The good news is that basically I already did the right thing, it's just not integrated into sympy. Integrating it is a relatively minor change. My class is structured just like PolyRing/PolyElement, but less sophisticated and having the problems I mentioned, and the coercion problems that I'm having basically just require subclassing Ring and Element, using EXRAW (a caveat I probably would've missed), and we should be good. I will report back. Thanks a lot.

To view this discussion visit https://groups.google.com/d/msgid/sympy/CAHVvXxR%2B7hCdGesGJ%3DFu3rpTDEQX6JA4BPb3kY4077jO6MpFAw%40mail.gmail.com.

[Matthew J. Samuel]

OK, I thought about this. What if it's both a Ring and an Expr? Is that a thing?

For example, the way I "resolved" (not quite) my original issue is I have EXRAW as the base ring, and radd converts it to an Expr, and I have "as_expr" have it be a sum of Expr objects that I use to stringify the object. Can I use those Exprs just like the ring? So basically there would be two representations, the Expr representation and the Ring representation, and if you have a need to sum them together, one of them is a ring, one is an Expr, and you're sacrificing some speed to do this weird thing just to have consistency.

What it looks like now is,

But then you multiply on the right,

which is expected because that DSx((3,1,4,2), y) only superficially looks the same, it's a shell Expr with no __mul__ method.

[Matthew J. Samuel]

I have a possible solution for which the base class could be independent of Schubert polynomials entirely and actually could quite likely be useful in sympy itself!

We could take the tensor product of domains. This unifies the coefficient rings and joins the generators, which is something that is already there, but also allows us to construct a new basis from the old ones so equality testing is still easy. I will get back to you if I have anything useful to contribute to sympy. Thanks. 

[Aaron Meurer]

You're better off making conversion to Expr for this type of object be explicit. An issue with Expr is that it's slow compared to bespoke classes like Poly, because it is designed to represent arbitrary symbolic expressions, so it can't have any of the optimizations that make the specialized classes faster. You don't want to have an implicit conversion into a slower type.

Aaron Meurer

To view this discussion visit https://groups.google.com/d/msgid/sympy/CAF%2B7YfGp04tt_-RwDu9BmfUF%3DUz-etj58j5yhU8FeoWdJPh66g%40mail.gmail.com.