// load the package
path_to_package := "...";
if path_to_package eq "..." then
    printf "\n\n************\n* You must edit the variable path_to_package vefore uing this file!\n************\n\n";
    quit;
end if;
    
AttachSpec(path_to_package cat "qpoly.spec");



// EXAMPLE: Bezout coefficients
// define the ring of polynomials over a quaternion algebra
K<t>     := CyclotomicField(7);
A<i,j,k> := QuaternionAlgebra<K | -t+1, -t+3>;
Ax<x>    := PolynomialRing(A);
// take two polynomials
p := (t - 1)*x^6 + (t^2 - 1)*x^5 + ((2*t^2 - 5*t + 3) + (t^2 - 2*t + 1)*i + j)*x^4 + ((-2*t^2 + 2*t) + (2*t - 2)*i + (t + 1)*j - k)*x^3 + (2*t - 3)*j*x^2 + ((t^2 - 2*t - 1)*j + k)*x + (-2*t + 2)*j - 2*t*k;
q := x^5 + ((t + 2) - i)*x^4 + (3*t - 2*i + t*k)*x^3 + (t^2 + (-t + 4)*i + (2*t^2 - 2*t)*j + (t^2 + t)*k)*x^2 + ((4*t - 8) + 4*t*i + (t^3 - t)*j - t*k)*x + -4*t + 4*i + (-2*t^2 + 2*t)*j - 2*t^2*k;
// compute their greatest common RIGHT divisor
// and the Bezout coeffs
r, alpha, beta := XGCRD(p, q); r;
// check the correctness
IsZero(p mod r) and IsZero(q mod r);
alpha*p + beta*q eq r;



// EXAMPLE: ideal arithmetic
// define the ring of polynomials over a quaternion algebra
K<sqrt7> := QuadraticField(7);
A<i,j,k> := QuaternionAlgebra<K | -1, -1>;
Ax<x>    := PolynomialRing(A);
// take two ideals of opposite handedness
I := RightIdeal(Ax, (x^2 + 1)*(x + i + k)); I;
J := LeftIdeal(Ax, (x + j)*(x^2 + 2)); J;
// intersect them
I meet J;
// the reason is:
IdealBound(I) eq IdealBound(J); 




// EXAMPLE: Roots
// define the ring of polynomials over a quaternion algebra
A<i,j,k> := QuaternionAlgebra<Rationals() | -1, -1>;
Ax<x>    := PolynomialRing(A);
// and a polynomial p
p        := Ax![ 6+6*i-36*j+24*k, 18-12*i-6*j+30*k,  
            23+5*i-30*j+20*k, 15-10*i-5*j+25*k, 
            16+i-6*j+4*k, 3-2*i-j+5*k, 3 ];
// find the roots of p
Zp       :=  Roots( p ); Zp;
// verify the roots
[ Evaluate(p, xi[1]) : xi in Zp ];




// EXAMPLE: Factorization
// define the ring of polynomials over a quaternion algebra
K<sqrt7> := QuadraticField(7);
A<i,j,k> := QuaternionAlgebra<K | -1, -1>;
Ax<x>    := PolynomialRing(A);
// take a poly given as a product
r1 := (x^2 - (3+i)*x - k);
r2 := (x - j);
r3 := (x^2 + (1-k)*x + (1-k));
r4 := (k*x - (1+k));
p := r1*r2*r3*r4;
// factor it and see that the factors are different
fp, c := Factorization(p); fp; c;
// yet still, it is a valid factorization, just a different one
// (we are in a non-commutative ring)
c*FactorisationToPolynomial(fp) eq p;
// verify that the factors are similar
test, lambda := IsSimilar(r3, fp[2][1] : Isomorphism := true);
test;



// EXAMPLE: Sum of squares (s.o.s.)
// take a rational polynomial
Qx<x> := PolynomialRing(Rationals());
p := x^6 + 6*x^5 + 8*x^4 - 26*x^3 + 50*x^2 - 22*x + 9; p;
// factor it over (-1,-1 | QQ)
A<i,j,k> := QuaternionAlgebra<Rationals() | -1, -1>;
factors := Factorization( PolynomialRing(A) ! p ); factors;
// the coordinates of any of the two factors form the decomposition into s.o.s.
C := Coordinates(factors[1][1]); C;
// verification
&+ [ q^2 : q in C ] eq p;
