Add quotient fields

[fomichev-d]

For rational functions, for instance.

It seems that would mainly involve:

• copying the RationalField<I> implementation;
• extending it from I::Type: IntegerRing to at least I::Type: EuclideanRing;
• replacing signed_gcd with gcd and signed_lcm with lcm.

[FeanorTheElf]

Yes, that's a project I have been thinking about for a while. Your proposal is one reasonable way on how to implement this (the bounds on the base ring should be PrincipalIdealRing + Domain however). There are also a few variants that might be worth thinking about:

• Weaken the bounds to DivisibilityRing + Domain, which would mean we cannot fully reduce fractions (but there is DivisibilityRing::balance_factor())
• Specialized implementations for the case of rational functions etc. - this might allow for much more efficient implementations at some points (?)
• Should it have a corresponding trait, so as to allow multiple different implementations?

I didn't do it so far, since I didn't need it, and without a use case, I find it very hard to answer these kind of questions. It would be very helpful if you could shortly explain what you need it for, and what kind of design would fit your use case. Generally speaking, fraction fields certainly should be part of feanor-math.

[fomichev-d]

the bounds on the base ring should be PrincipalIdealRing + Domain however

Then implementations of gcd and lcm supporting this will be needed. Currently it requires EuclideanRing.

Specialized implementations for the case of rational functions etc. - this might allow for much more efficient implementations at some points (?)

Specialised implementations might be a good idea. I guess there could be a more efficient QuotientFieldImplBase<I> with stronger trait bounds on I and some GeneralQuotientFieldImplBase<I> with weaker ones. Having a trait QuotientField to unite them and making RationalFieldBase<I> another special case would also make sense.

It would be very helpful if you could shortly explain what you need it for, and what kind of design would fit your use case.

Sure! The brief answer is that I want to solve linear systems of equations over rational functions. My main types (QuotientField is quickly hacked together from RationalField as I described above) look like this:

pub type ZZType = MPZ;
pub const ZZ: MPZ = MPZ::RING;

pub type QQType = RationalField<ZZType>;
pub type QQEl = RationalFieldEl<ZZType>;
pub const QQ: RationalField<MPZ> = RationalField::new(ZZ);

pub type Poly1Ring = DensePolyRing<&'static QQType>;
pub type Poly1RingEl = DensePolyRingEl<&'static QQType>;
pub type RatPoly1Ring = QuotientField<Poly1Ring>;
pub type RatPoly1RingEl = QuotientFieldEl<Poly1Ring>;

A little more detailed and contextualised description of my use case is:

I have a csv table with precomputed values of a conjectured function from graphs to polynomials in one variable with rational coefficients. I want to find graph relations holding for these values. In my code, I define graph templates, graph relations and their evaluation. A graph relation has a coefficient for every term - to find them, I use a separate variable for each, evaluate the expression as a multivariate polynomial and add it to the system I'm solving for the coefficients. There's a large number of possible evaluations so I have an iterator of them and I want to abort early if it's clear there're no solutions.

Previously I was using SageMath, but I've rewritten it in Rust because I needed to make it much faster than it was, including keeping all CPU cores busy (the current speedup for the same subtask is ~24h (maybe more) -> ~40s now). In Sage, I was reducing everything using ideals and Groebner bases and would then call solve() in the end; not the most efficient approach, but it was unavoidable because I was doing something more general there and I actually had three variables for every term. Naturally, I found this crate as one that actually implements Groebner bases in Rust and my initial rewrite used the same approach and called Maxima to do the final solve(). Recently I was optimising it more and I started solving a linear system for the coefficients directly - hence the need for rational functions as a field.

I solve the linear system by (currently straightforward, maybe I should consider the Bareiss algorithm) Gaussian elimination and maintaining the reduced row echelon form. I don't think I'd use a default implementation of this solver even if there were one: the coefficients (rational functions in this case; my implementation is over any field) explode quite a bit in the process (until maybe ~60-70% of the system is solved) and I need to use rayon's parallelism to speed up all the row reduction operations. One typical example I've just seen is an intermediate relation resulting as a solution to a partially solved system for a large relation: that intermediate relation takes >110 full-screen screens just to print.

Oh, and it'd be nice if the printing were better than the one in my implementation (where I just replaced the / separator with |; still no parentheses around the polynomials where needed) and things like 1|-1/2 were simplified properly.

[FeanorTheElf]

Very interesting, thank you for the information!

I think I will go ahead and implement quotient fields, using a trait and a default implementation (as currently done e.g. for polynomial rings).

Then implementations of gcd and lcm supporting this will be needed. Currently it requires EuclideanRing.

I guess this not clear at all from the documentation, but actually PrincipalIdealRing::ideal_gen() is the function we want here, and it already exists. For polynomials over some rings (like the rationals), it will also be much faster than gcd(), since the latter just uses the Euclidean algorithm, while the former does a reduce-mod-p-and-lift approach. I will adjust the documentation accordingly!

I don't think I'd use a default implementation of this solver even if there were one:

There is one in feanor-math, via the trait LinSolveRing. However, it might indeed not be fast enough for your task. Generally speaking, there are specialized algorithms to solve linear systems over the rationals (and that avoid the coefficient explosion, a common approach is to compute a solution modulo many primes p and then reconstruct the solution over the integers. Unfortunately, this is currently not implemented in feanor-math for linear systems (another project on my todo list). For polynomials, one can also compute a solution for the system evaluated at many points, and reconstruct the polynomial solution from that.

[fomichev-d]

I guess this not clear at all from the documentation, but actually PrincipalIdealRing::ideal_gen() is the function we want here, and it already exists. For polynomials over some rings (like the rationals), it will also be much faster than gcd(), since the latter just uses the Euclidean algorithm, while the former does a reduce-mod-p-and-lift approach. I will adjust the documentation accordingly!

I see, thank you! It makes some things 2x faster already.

There is one in feanor-math, via the trait LinSolveRing. However, it might indeed not be fast enough for your task.

Oh, I somehow missed it. But yeah, I think I need my own here and I have some ideas on improving it.

[FeanorTheElf]

Version 3.2.0 has now (with feature unstable-enable) an implementation of fraction fields, via FractionFieldImpl. However, there are still some unsolved questions.

• Currently, this implementation supports all rings R: Domain + DivisibilityRing. In particular, this means that it cannot and won't fully reduce fractions. Instead, it currently relies on DivisibilityRing::balance_factor(), which however for polynomials does not actually compute gcds. It would also be valuable to know how much faster it will get if we correctly reduce by computing polynomial gcds.
• Would it make sense to add a special implementation for fraction fields of polynomial rings? This would solve 1. in this case at least.
• Is there anything that is missing in the interface, or that would be nice to have?

If you have any opinions on the current design, I'd be happy to consider them. Otherwise, I will do some benchmarking when I find time and figure out how to proceed.

[fomichev-d]

Thank you for working on this!

• Would it be difficult/undesirable for some reason to compute gcd in DivisibilityRing::balance_factor() for polynomials? It would seem this should be the proper thing to do.
• I imagine one wouldn't need a special implementation for fraction fields for polynomial rings if balance_factor() also multiplied the factor by such a unit that the last non-zero polynomial in the iterator became monic. But it'd require a method like PolyRing::LT() (which would make sense to have in any case) and it'd also probably be only implemented this way over fields probably? Maybe a separate trait FieldPolyRing with an automatic optimised implementation could be added just for this if specialisation is too messy in Rust yet.
• Generally, it'd be good for divisibility rings to be able to account for things such as LT() and so on for various ring extensions where everything is minimised modulo some unit now - if the ring can determine what a normalised form of its element is; it could even work for fraction fields themselves, but I'm not sure how I'd design traits around this. Maybe it could be a trait with a method producing an El<Self> from an El<Self> such that it guarantees that the resulting element is a unit and that dividing by this unit makes the element normalised in some sense specific to the ring? Maybe it could be a provided method in DivisibilityRing itself with the default implementation of self.one(). I think that last option seems to make the most sense.

Unfortunately, I can't test everything extensively right now as I have a series of computations that'll each take maybe days to complete and I don't want to potentially disrupt them by changing something right now. But the implementation looks good to me visually modulo the behaviour of balance_factor().

[FeanorTheElf]

Computing the polynomial gcd in DivisibilityRing::balance_factor() is somewhat difficult, since polynomial rings can be instantiated for rings where the gcd is not defined (i.e. not fields). Hence, we'd need specialization for that. My current manual implementation/workaround for specialization should work, although it is very messy. But as long as builtin specialization in Rust is as unstable as it currently is, I need some amount of workarounds anyway.

The other point is that in general, computing the gcd of two polynomials is a very expensive operation (especially over the rationals), and it might harm performance in some cases instead of helping. However, I haven't benchmarked it so far, hence I'm not sure.

Nevertheless I think you are right, and DivisibilityRing::balance_factor() for polynomials should generally compute the gcd. I'll do some more research once I find time, but I think adding some kind of specialization to achieve this is the best way forward.

Regarding your other point: I don't think we need to normalize the leading term of the denominator. Is there any reason why you need a canonical representation of you fraction field elements? The only thing I can think of is if you intend to hash them, but otherwise reduction is "just" a performance optimization. However, since we need specialization anyway to compute the gcd in balance_factor(), it wouldn't be too difficult to also make the last polynomial monic/primitive/whatever other notion, if you require that.

[fomichev-d]

Nice!

As of normalisation, my main arguments for it are:

• it may help printed results be more uniform and easily interpretable;
• other than fractions, having some sort of universal element normalisation API in DivisibilityRing in general (where it'd mean different things for different rings) would make it easier to analyse and interpret a ring element structurally in code, should one need that, and also to get rid of a constant multiplier if it doesn't matter for a particular use case.