Solving systems of polynomials with complex coefficients

[alexfikl]

I have a system of polynomials of some moderate size that needs solving, i.e. find all roots. It would be in 2 or 3 variables (i.e. $x, y, z$) with orders of 2, 4, 16 (probably not higher).

Can this be done with this library? Is it in scope? There seems to be some support for complex numbers and Gröbner bases, but I'm not very knowledgeable about how to tie it all together.

From what I can tell, flint also only has support for single variable polynomials in this case, so I imagine it's quite difficult to get right. https://python-flint.readthedocs.io/en/latest/acb_poly.html Annoyingly, Mathematica seems to be doing a very good job at it.

Any thoughts / suggestions would be very welcome! Thank you and congrats for a very nice looking and sorely needed library :grin:

[FeanorTheElf]

Groebner basis are certainly the right tool, and with 2 or 3 variables, performance should not be an issue.

However, the main problem are indeed the complex numbers. Currently there only is an implementation of the complex numbers based on floating point numbers, but that won't work with Groebner basis, as these algorithms are usually very numerically unstable (i.e. rounding errors explode). Instead, the standard way is to solve the system not over C but over the algebraic numbers (assuming your coefficients are algebraic numbers), since these can be represented exactly.

In other words, the question is, in what field are the coefficients of your polynomials? If they are in Q or a known number field, it is (with a little bit of additional code) possible to solve the system in the current version of feanor-math. If they are general algebraic numbers, you would first have to construct the number field that contains all of them, which is a nontrivial task in itself. If they are non-algebraic numbers, things get difficult, although it is often possible to introduce additional indeterminates for transcendent elements.

If you tell me which case we are in, I am happy to give you more details/code examples. Maybe I will also include one in the main Readme.

PS: I am currently refactoring the code to get to version 3.0.0 (there were some breaking changes I wanted to make for quite a while). This includes quite some changes to the interface of multivariate polynomial rings and Groebner basis, but the functionality remains the same (and it hopefully gets faster as well).

[alexfikl]

Thank you for the quick reply!

However, the main problem are indeed the complex numbers. Currently there only is an implementation of the complex numbers based on floating point numbers, but that won't work with Groebner basis, as these algorithms are usually very numerically unstable (i.e. rounding errors explode). Instead, the standard way is to solve the system not over C but over the algebraic numbers (assuming your coefficients are algebraic numbers), since these can be represented exactly.

In other words, the question is, in what field are the coefficients of your polynomials? If they are in Q or a known number field, it is (with a little bit of additional code) possible to solve the system in the current version of feanor-math. If they are general algebraic numbers, you would first have to construct the number field that contains all of them, which is a nontrivial task in itself. If they are non-algebraic numbers, things get difficult, although it is often possible to introduce additional indeterminates for transcendent elements.

If you tell me which case we are in, I am happy to give you more details/code examples. Maybe I will also include one in the main Readme.

I'm afraid my coefficients are in $\mathbb{C}$, but the real and imaginary parts could be assumed rational, so maybe there's some hope there?

The underlying problem that I'm trying to solve is finding the fixed points of (compositions of) a map $f: (A z) * (A z) + c$, where $A$ is a real (even rational) $d \times d$ matrix and $c \in \mathbb{C}$. The multiplication there is supposed to be pointwise. Unfortunately, if you compute iterates $f \circ f \circ f \cdots \circ f$ the complex $c$ constant gets into the coefficients, so I can't assume they're rational..

PS: I am currently refactoring the code to get to version 3.0.0 (there were some breaking changes I wanted to make for quite a while). This includes quite some changes to the interface of multivariate polynomial rings and Groebner basis, but the functionality remains the same (and it hopefully gets faster as well).

That's great to hear! I'll keep an eye on it :grin:

[FeanorTheElf]

I'm afraid my coefficients are in C , but the real and imaginary parts could be assumed rational, so maybe there's some hope there?

That is perfect, since "rational numbers with i" form a very simple number field! In other words, what you want can more or less be done with feanor-math, although the fact that it has very bad floating-point support makes it a little bit ugly. I am assuming here that you want the solution as x + yi with x, y floating point numbers. The alternative would be to again represent them as elements of a number field, which is nice because the representation is exact, but again building the number field is quite a lot of effort - I will consider adding this to the library.

I have prepared some sample code, using the current version of the library on github (i.e. doesn't work with 2.3.0). It works, with the following caveats:

• Assumes that the solution space consists of isolated points, and is not a curve. However, it shouldn't be too hard to also support curves.
• Assumes that the Lex-GB only has three elements. This is likely to be the case for most "simple" systems. Generalizing should be easy as well, except for the fact that it requires gcd computations. This would be supported over any number field (again requiring to build the correct one), but is not a good idea over C (again because of numerical stability).

Side note: If you know that your system has a solution, and is somewhat well-behaved, you might try just using the multivariate Newton-Raphson method. This is not available in feanor-math, but should be only a few lines to implement in the 3-variable case, with the most difficult part being the inversion of a 3x3 matrix. Also, the Symbolica Rust library has similar functionality as feanor-math when it comes to multivariate polynomials, but might have better floating-point support.

use feanor_math::{field::FieldStore, homomorphism::{Homomorphism, LambdaHom}, integer::BigIntRing, rings::{extension::{extension_impl::FreeAlgebraImpl, FreeAlgebraStore}, field::AsField, float_complex::Complex64, multivariate::{multivariate_impl::MultivariatePolyRingImpl, Lex, MultivariatePolyRingStore}, poly::{dense_poly::DensePolyRing, derive_poly, PolyRingStore}, rational::RationalField}, seq::{VectorFn, VectorView}};
use feanor_math::ring::*;

// the rationals
type QQTy = RationalField<BigIntRing>;
// degree-2 number field, will be "rationals + i"
type QQiTy = AsField<FreeAlgebraImpl<QQTy, [El<QQTy>; 1]>>;
// 3-variate polynomial ring over `QQi`
type QQiXYZTy = MultivariatePolyRingImpl<QQiTy>;
// univariate polynomial ring over floating-point based complex numbers
type CCXTy = DensePolyRing<Complex64>;

///
/// Solves a univariate complex polynomial using the Newton-Raphson method.
/// This of course only gives approximate results
/// 
fn solve_univariate_poly(CCX: &CCXTy, f: &El<CCXTy>) -> El<Complex64> {
    println!("Finding root of univariate poly {}", CCX.format(f));
    let CC = CCX.base_ring();
    let f_prime = derive_poly(&CCX, f);
    let mut current = CC.add(CC.one(), Complex64::I);
    while CC.abs(CCX.evaluate(f, &current, &CC.identity())) > 0.00000001 {
        current = CC.sub(current, CC.div(
            &CCX.evaluate(f, &current, &CC.identity()),
            &CCX.evaluate(&f_prime, &current, &CC.identity())
        ));
    }
    println!("Found {}", CC.format(&current));
    return current;
}

fn find_any_solution_if_exists(QQiXYZ: &QQiXYZTy, system: &[El<QQiXYZTy>]) -> Option<[El<Complex64>; 3]> {
    let mut gb = feanor_math::algorithms::buchberger::buchberger_simple::<_, _, false>(QQiXYZ, system.iter().map(|f| QQiXYZ.clone_el(f)).collect(), Lex);

    // if the GB contains `1`, the system is unsolvable
    if gb.len() == 1 && QQiXYZ.appearing_indeterminates(&gb[0]).into_iter().all(|(_, deg)| deg == 0) {
        return None;
    }

    // the Lex-Groebner basis always has the form `f1(Z), f2(Y, Z), f3(X, Y, Z)`;
    // this allows us to solve the system by only solving a univariate polynomial at a time;
    // in particular, we first find `z` as root of `f1(Z)`, then plug it in in `f2`, 
    // and then find `y` as root of `f2(Y, z)`

    // order them so that `f1` comes last, i.e. we take polynomials from the end, step-by-step
    gb.sort_unstable_by_key(|f| *QQiXYZ.appearing_indeterminates(f).iter().map(|(i, _)| i).min().unwrap() as i64);

    let CC = Complex64::RING;
    let ZZ_to_CC = CC.can_hom(QQiXYZ.base_ring().base_ring().base_ring()).unwrap();
    let QQi_to_CC = |x: &El<QQiTy>| {
        let x_wrt_basis = QQiXYZ.base_ring().wrt_canonical_basis(x);
        CC.add(
            CC.div(&ZZ_to_CC.map(x_wrt_basis.at(0).0), &ZZ_to_CC.map(x_wrt_basis.at(0).1)),
            CC.mul(CC.div(&ZZ_to_CC.map(x_wrt_basis.at(1).0), &ZZ_to_CC.map(x_wrt_basis.at(1).1)), Complex64::I)
        )
    };

    let CCX = DensePolyRing::new(CC, "T");
    let QQi_to_CCX = LambdaHom::new(QQiXYZ.base_ring(), &CCX, |_, _, c| CCX.inclusion().map(QQi_to_CC(c)));

    let f1 = gb.pop().unwrap();
    let variables = QQiXYZ.appearing_indeterminates(&f1);
    // there is a special case that we currently don't handle here: If the solution space of the system is not some
    // isolated points, but a curve, then `f1` will be missing - we can instead pick any `z` we want. However, the
    // "isolated point" case is most common
    assert!(variables.len() == 1);
    assert!(variables.iter().next().unwrap().0 == 2);

    let z = solve_univariate_poly(&CCX, &QQiXYZ.evaluate(&f1, [CCX.zero(), CCX.zero(), CCX.indeterminate()].into_ring_el_fn(&CCX), &QQi_to_CCX));

    let f2 = gb.pop().unwrap();
    assert!(QQiXYZ.appearing_indeterminates(&f2).iter().next().unwrap().0 == 1);

    let f2 = QQiXYZ.evaluate(&f2, [CCX.zero(), CCX.indeterminate(), CCX.inclusion().map(z)].into_ring_el_fn(&CCX), &QQi_to_CCX);
    let y = solve_univariate_poly(&CCX, &f2);

    // now do exactly the same for f3
    let f3 = gb.pop().unwrap();
    assert!(QQiXYZ.appearing_indeterminates(&f3).iter().next().unwrap().0 == 0);

    let mut f3 = QQiXYZ.evaluate(&f3, [CCX.indeterminate(), CCX.inclusion().map(y), CCX.inclusion().map(z)].into_ring_el_fn(&CCX), &QQi_to_CCX);
    let x = solve_univariate_poly(&CCX, &f3);

    return Some([x, y, z]);
}

fn main() {

    let QQ = RationalField::new(BigIntRing::RING);
    let QQi = FreeAlgebraImpl::new(QQ.clone(), 2, [QQ.neg_one()]).as_field().ok().unwrap();
    let QQiXYZ = MultivariatePolyRingImpl::new(QQi.clone(), 3);

    let CC = Complex64::RING;
    let ZZ_to_CC = CC.can_hom(QQiXYZ.base_ring().base_ring().base_ring()).unwrap();
    let QQi_to_CC = LambdaHom::new(QQi, CC, |_, _, x: &El<QQiTy>| {
        let x_wrt_basis = QQiXYZ.base_ring().wrt_canonical_basis(x);
        CC.add(
            CC.div(&ZZ_to_CC.map(x_wrt_basis.at(0).0), &ZZ_to_CC.map(x_wrt_basis.at(0).1)),
            CC.mul(CC.div(&ZZ_to_CC.map(x_wrt_basis.at(1).0), &ZZ_to_CC.map(x_wrt_basis.at(1).1)), Complex64::I)
        )
    });

    // for simplicity, let's just take rational coefficients, as building the polynomials gets much easier.
    // However, `QQi.canonical_gen()` is `i`, so we can use complex coefficients
    let system = QQiXYZ.with_wrapped_indeterminates(|[X, Y, Z]| {
        [
            X * X + 10,           // X^2 = 10
            X * Y * Z * Z - 2,    // XYZ^2 = 1
            Z * Z * Z - X,        // Z^3 = Y
        ]
    });

    let sol = find_any_solution_if_exists(&QQiXYZ, &system).unwrap();

    for f in &system {
        Complex64::RING.println(&QQiXYZ.evaluate(f, sol.into_ring_el_fn(Complex64::RING), &QQi_to_CC));
    }
}

[FeanorTheElf]

Ahh, having a look at your profile I saw that you might be interested in finding all roots. This makes things slightly different of course...

The approach above mainly still works, but solve_univariate_poly() would have to find all roots. Since indeed Newton-Raphson is badly suited for this, we would do the following:

• Factor the polynomial over the current number field (already supported in feanor-math)
• Using the irreducible factors, construct a number field K containing all solutions
• Find a single complex x root using Newton-Raphson of the generating polynomial of K
• Now map the algebraic representations of all roots to the complex numbers, using x

Step 2 here is quite some effort, I think I will implement a function for it in feanor-math. Stay tuned :)

[alexfikl]

Ahh, having a look at your profile I saw that you might be interested in finding all roots. This makes things slightly different of course...

Yes, I unfortunately need all of them. Thank you for the code sample though! It's very useful to get an idea of how to build the necessary things.

The approach above mainly still works, but solve_univariate_poly() would have to find all roots. Since indeed Newton-Raphson is badly suited for this, we would do the following:

* Factor the polynomial over the current number field (already supported in feanor-math)

* Using the irreducible factors, construct a number field `K` containing all solutions

* Find a single complex `x` root using Newton-Raphson of the generating polynomial of `K`

* Now map the algebraic representations of all roots to the complex numbers, using `x`

I was thinking of modifying the code to use something like https://github.com/ickk/aberth/ to find all the roots for the univariate case and the move upwards. Would that work?

I've tried using Newton-Raphson to find all the roots and it mostly just worked up to 4th degree polynomials (to find 16 roots), but couldn't find all of them for higher degrees anymore, which is fair..

Also, the Symbolica Rust library has similar functionality as feanor-math when it comes to multivariate polynomials, but might have better floating-point support.

Yeah, I've taken a look at it following your mention in the README, but their pricing / usage model doesn't really work for me :(

Step 2 here is quite some effort, I think I will implement a function for it in feanor-math. Stay tuned :)

That would be lovely!

[FeanorTheElf]

Ok, I have implemented extend_splitting_field() in feanor-math now, although there is a lot left to do...

I was thinking of modifying the code to use something like https://github.com/ickk/aberth/ to find all the roots for the univariate case and the move upwards. Would that work?

That sounds like a great idea! Working over the complex numbers is always preferrable performance-wise, as long as the algorithms have good numerical stability. I am not familiar with the Aberth method, but it looks good, although it might fail if two different complex roots are very close together?.

This is never a problem with a purely algebraic approach, however after implementing it in feanor-math, I ran into huge performance problems. Even very simple cases lead to rational numbers with huge numerator and denominator. There are algorithms designed to counter this, but I have not consistently implemented them in feanor-math. This is planned, but a large project, and will take lots of effort.

Anyway, the following modification of the code above does find all roots of a system. However, its performance is terrible, and I have doubt that it could currently solve even the simple cases you have in mind. This will change once I implement better gcd algorithms for number fields in feanor-math.

#![allow(non_snake_case)]
#![feature(allocator_api)]

use std::alloc::Global;

use feanor_math::{algorithms::{convolution::STANDARD_CONVOLUTION, splitting_field::extend_splitting_field}, field::FieldStore, homomorphism::{Homomorphism, LambdaHom}, integer::BigIntRing, rings::{extension::{extension_impl::FreeAlgebraImpl, FreeAlgebraStore}, field::{AsField, AsFieldBase}, float_complex::Complex64, multivariate::{multivariate_impl::MultivariatePolyRingImpl, Lex, MultivariatePolyRingStore}, poly::{dense_poly::DensePolyRing, derive_poly, PolyRingStore}, rational::*}, seq::{VectorFn, VectorView}};
use feanor_math::ring::*;
use feanor_math::pid::PrincipalIdealRingStore;

// floating-point based complex numbers
type CCTy = Complex64;
// the rationals
type QQTy = RationalField<BigIntRing>;
// degree-2 number field, will be "rationals + i"
type QQiTy<'a> = AsField<FreeAlgebraImpl<&'a QQTy, Vec<El<QQTy>>>>;
// 3-variate polynomial ring over `QQi`
type QQiXYZTy<'a> = MultivariatePolyRingImpl<QQiTy<'a>>;
// univariate polynomial ring over floating-point based complex numbers
type CCXTy = DensePolyRing<Complex64>;

///
/// Solves a univariate complex polynomial using the Newton-Raphson method.
/// This of course only gives approximate results
/// 
fn find_one_complex_solution(CCX: &CCXTy, f: &El<CCXTy>) -> El<Complex64> {
    let CC = CCX.base_ring();
    let f_prime = derive_poly(&CCX, f);
    let mut current = CC.add(CC.one(), Complex64::I);
    while CC.abs(CCX.evaluate(f, &current, &CC.identity())) > 0.00000001 {
        current = CC.sub(current, CC.div(
            &CCX.evaluate(f, &current, &CC.identity()),
            &CCX.evaluate(&f_prime, &current, &CC.identity())
        ));
    }
    return current;
}

///
/// Computes a map from some number field K into the complex numbers
///
fn K_to_CC<'a, 'b>(K: &'a AsField<FreeAlgebraImpl<&'b QQTy, Vec<El<QQTy>>>>) -> impl 'a + Fn(&El<AsField<FreeAlgebraImpl<&'b QQTy, Vec<El<QQTy>>>>>) -> El<Complex64> {
    let CC = Complex64::RING;
    let QQ = K.base_ring();
    let CCX = DensePolyRing::new(CC, "X");
    let QQ_to_CC = CC.into_can_hom(QQ).ok().unwrap();
    let gen_poly = K.generating_poly(&CCX, &QQ_to_CC);
    let complex_sol = find_one_complex_solution(&CCX, &gen_poly);
    return move |x| CCX.evaluate(&K.poly_repr(&CCX, x, &QQ_to_CC), &complex_sol, &CC.identity());
}

///
/// Computes a map from `QQ[i]` into some number field K, assuming we have given the value `i` of K
/// 
fn QQi_to_K<'a, 'b>(QQi: &'b QQiTy<'a>, K: &'b AsField<FreeAlgebraImpl<&'a QQTy, Vec<El<QQTy>>>>, i: &'b El<AsField<FreeAlgebraImpl<&'a QQTy, Vec<El<QQTy>>>>>) -> impl 'b + Homomorphism<AsFieldBase<FreeAlgebraImpl<&'a QQTy, Vec<El<QQTy>>>>, AsFieldBase<FreeAlgebraImpl<&'a QQTy, Vec<El<QQTy>>>>> {
    fn map_from_QQi_to_K<'a>(QQi: &QQiTy<'a>, K: &AsField<FreeAlgebraImpl<&'a QQTy, Vec<El<QQTy>>>>, i: &El<AsField<FreeAlgebraImpl<&'a QQTy, Vec<El<QQTy>>>>>, x: &El<QQiTy<'a>>) -> El<AsField<FreeAlgebraImpl<&'a QQTy, Vec<El<QQTy>>>>> {
        let x_wrt_basis = QQi.wrt_canonical_basis(x);
        let result = K.add(K.inclusion().map(x_wrt_basis.at(0)), K.inclusion().mul_ref_map(i, &x_wrt_basis.at(1)));
        return result;
    }
    LambdaHom::new(QQi, K, move |from, to, x| map_from_QQi_to_K(from, to, i, x))
}

///
/// Computes a map `QQ[i]` into the complex numbers.
/// 
/// Technically this is just a special case of `K_to_CC`, but can be implemented much simpler
///
fn QQi_to_CC<'a, 'b, 'c, 'd>(QQi: &'b &QQiTy<'a>, CC: &'c CCTy, x: &'d El<QQiTy<'a>>) -> El<CCTy> {
    let x_wrt_basis = QQi.wrt_canonical_basis(x);
    let QQ = QQi.base_ring();
    let QQ_to_CC = CC.can_hom(&QQ).unwrap();
    CC.add(
        QQ_to_CC.map(x_wrt_basis.at(0)),
        CC.mul(QQ_to_CC.map(x_wrt_basis.at(1)), Complex64::I)
    )
}

///
/// Computes a number field K that contains all solution to the given system, together with these solutions
/// 
fn find_all_solutions<'a>(QQiXYZ: &QQiXYZTy<'a>, system: &[El<QQiXYZTy<'a>>]) -> (AsField<FreeAlgebraImpl<&'a QQTy, Vec<El<QQTy>>>>, Vec<[El<AsField<FreeAlgebraImpl<&'a QQTy, Vec<El<QQTy>>>>>; 3]>) {

    let QQi = QQiXYZ.base_ring();

    let mut gb = feanor_math::algorithms::buchberger::buchberger_simple::<_, _, false>(QQiXYZ, system.iter().map(|f| QQiXYZ.clone_el(f)).collect(), Lex);

    println!("Found Groebner basis");

    // if the GB contains a unit, the system is unsolvable
    if gb.len() == 1 && QQiXYZ.appearing_indeterminates(&gb[0]).into_iter().all(|(_, deg)| deg == 0) {
        return (QQiXYZ.base_ring().clone(), Vec::new());
    }

    // the Lex-Groebner basis always has the form `f1(Z), f21(Y, Z), ..., f2k(Y, Z), f31(X, Y, Z), ..., f3l(X, Y, Z)`;
    // this allows us to solve the system by only solving a univariate polynomial at a time;
    // in particular, we first find `z` as root of `f1(Z)`, then plug it into `f21(Y, Z), ..., f2k(Y, Z)`
    // and compute `y` as a root of their gcd. Then continue for `X` and `f31, ..., f3l`

    // order them so that `f1` comes last, i.e. we take polynomials from the end, step-by-step
    gb.sort_unstable_by_key(|f| *QQiXYZ.appearing_indeterminates(f).iter().map(|(i, _)| i).min().unwrap() as i64);

    let f1_poly = gb.pop().unwrap();
    println!("f1 = {}", QQiXYZ.format(&f1_poly));
    let variables = QQiXYZ.appearing_indeterminates(&f1_poly);
    // there is a special case that we currently don't handle here: If the solution space of the system is not some
    // isolated points, but a curve, then `f1` will be missing - we can instead pick any `z` we want. However, the
    // "isolated point" case is most common
    assert!(variables.len() == 1);
    assert!(variables.iter().next().unwrap().0 == 2);

    let mut f2_polys = Vec::new();
    let mut f3_polys = Vec::new();
    while let Some(f) = gb.pop() {
        if QQiXYZ.appearing_indeterminates(&f).iter().next().unwrap().0 == 1 {
            f2_polys.push(f);
        } else {
            f3_polys = gb;
            f3_polys.push(f);
            break;
        }
    }
    for (i, f) in f2_polys.iter().enumerate() {
        println!("f2{} = {}", i, QQiXYZ.format(&f));
    }
    for (i, f) in f3_polys.iter().enumerate() {
        println!("f3{} = {}", i, QQiXYZ.format(&f));
    }

    // K is the current number field that contains all solutions we found so far; in order to be able to map `QQ[i]` into
    // K, we also keep track of the value `i` within K. This will of course also change whenever we replace K with an extension field
    let mut K: AsField<FreeAlgebraImpl<&QQTy, Vec<El<QQTy>>>> = QQiXYZ.base_ring().clone();
    let mut i = K.canonical_gen();

    let mut result = Vec::new();

    let KiZ = DensePolyRing::new(RingRef::new(K.get_ring()), "Z");
    let f1 = QQiXYZ.evaluate(&f1_poly, [KiZ.zero(), KiZ.zero(), KiZ.indeterminate()].into_ring_el_fn(&KiZ), &KiZ.inclusion());

    println!("Computing splitting field of {}", KiZ.format(&f1));
    result.push(i);
    let (new_K, all_z, result_mapped) = extend_splitting_field(&KiZ, vec![f1], Vec::new(), result);
    result = result_mapped;
    i = result.pop().unwrap();
    K = new_K;
    assert!(K.is_neg_one(&K.pow(K.clone_el(&i), 2)));

    for mut z in all_z {

        println!("Found root z = {}", K.format(&z));
        let KY = DensePolyRing::new(RingRef::new(K.get_ring()), "Y");
        let f2 = f2_polys.iter().fold(KY.zero(), |current, f| KY.ideal_gen(&current, &QQiXYZ.evaluate(&f, &[KY.zero(), KY.indeterminate(), KY.inclusion().map(K.clone_el(&z))].into_ring_el_fn(&KY), &KY.inclusion().compose(QQi_to_K(QQi, &K, &i)))));

        println!("Computing splitting field of {}", KY.format(&f2));
        result.push(i);
        result.push(z);
        let (new_K, all_y, result_mapped) = extend_splitting_field(&KY, vec![f2], Vec::new(), result);
        result = result_mapped;
        z = result.pop().unwrap();
        i = result.pop().unwrap();
        K = new_K;
        assert!(K.is_neg_one(&K.pow(K.clone_el(&i), 2)));
        assert!(K.is_zero(&QQiXYZ.evaluate(&f1_poly, &[K.zero(), K.zero(), K.clone_el(&z)].into_ring_el_fn(&K), &QQi_to_K(QQi, &K, &i))));

        for mut y in all_y {

            println!("Found root y = {}", K.format(&y));
            let KX = DensePolyRing::new(RingRef::new(K.get_ring()), "X");
            let f1 = f3_polys.iter().fold(KX.zero(), |current, f| KX.ideal_gen(&current, &QQiXYZ.evaluate(&f, &[KX.indeterminate(), KX.inclusion().map(K.clone_el(&y)), KX.inclusion().map(K.clone_el(&z))].into_ring_el_fn(&KX), &KX.inclusion().compose(QQi_to_K(QQi, &K, &i)))));

            println!("Computing splitting field of {}", KX.format(&f1));
            result.push(i);
            result.push(z);
            result.push(y);
            let (new_K, all_x, result_mapped) = extend_splitting_field(&KX, vec![f1], Vec::new(), result);
            result = result_mapped;
            y = result.pop().unwrap();
            z = result.pop().unwrap();
            i = result.pop().unwrap();
            K = new_K;
            assert!(K.is_neg_one(&K.pow(K.clone_el(&i), 2)));
            assert!(K.is_zero(&QQiXYZ.evaluate(&f1_poly, [K.zero(), K.zero(), K.clone_el(&z)].into_ring_el_fn(&K), &QQi_to_K(QQi, &K, &i))));
            assert!(K.is_zero(&QQiXYZ.evaluate(&f2_polys[0], &[K.zero(), K.clone_el(&y), K.clone_el(&z)].into_ring_el_fn(&K), &QQi_to_K(QQi, &K, &i))));

            for x in all_x {
                println!("Found root x = {}", K.format(&x));

                result.push(x);
                result.push(K.clone_el(&y));
                result.push(K.clone_el(&z));
            }
        }
    }

    assert!(result.len() % 3 == 0);
    let mut result_grouped = Vec::new();
    while result.len() > 0 {
        let [z, y, x] = [result.pop().unwrap(), result.pop().unwrap(), result.pop().unwrap()];
        result_grouped.push([x, y, z]);
    }

    return (K, result_grouped);
}

fn main() {
    let QQ = RationalField::new(BigIntRing::RING);
    let QQi: QQiTy = FreeAlgebraImpl::new_with(&QQ, 2, vec![QQ.neg_one()], "i", Global, STANDARD_CONVOLUTION).as_field().ok().unwrap();
    let QQiXYZ: QQiXYZTy = MultivariatePolyRingImpl::new(QQi.clone(), 3);

    // for simplicity, let's just take rational coefficients, as building the polynomials gets much easier.
    // However, `QQi.canonical_gen()` is `i`, so we can use complex coefficients
    let system = QQiXYZ.with_wrapped_indeterminates(|[X, Y, Z]| {
        [
            X * X + 10,
            Z * Y * X - 2, 
            Z - X - 1, 
        ]
    });

    let CC = Complex64::RING;
    let QQi_to_CC = LambdaHom::new(&QQi, CC, QQi_to_CC);

    let (extension_field, sols) = find_all_solutions(&QQiXYZ, &system);
    let K_to_CC = K_to_CC(&extension_field);

    println!();
    println!("Presenting solutions");
    println!();
    for sol in &sols {
        let [x, y, z] = [K_to_CC(&sol[0]), K_to_CC(&sol[1]), K_to_CC(&sol[2])];
        println!("Solution ({}, {}, {})", CC.format(&x), CC.format(&y), CC.format(&z));
        for f in &system {
            Complex64::RING.println(&QQiXYZ.evaluate(f, [x, y, z].into_ring_el_fn(Complex64::RING), &QQi_to_CC));
        }
        println!();
    }
}