# ======================================================================== #
#
# Copyright 2012, Gabor Hegedus, Josef Schicho, and Hans-Peter Schroecker
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
# ======================================================================== #

# Maple implementation of Example 1 of Hegedus, Schicho, and
# Schroecker: Factorization of Rational Curves in the Study Quadric.

read "qf.mpl":

# Choose two rotation quaternions.
h1 := [0,      0, 1 - e, 1 + e];
h2 := [1, -1 + e,     0, 1 + e];

# Compute P = (t-h1)*(t-h2), its conjugate polynomial Pb, and multiply
# the two. The result is a real polynomial f of degree four with two
# irreducible quadratic factors M1, M2 that can be computed as minimal
# polynomials of h1 and h2, respectively.
P := PolMult([-h1], [-h2]);
Pb := map(Cj, P);       # conjugate quaternion
PPb := PolMult(P, Pb);  # a real polynomial
f := t^4:
for i from 1 to nops(PPb) do
  f := f + t^(4-i) * op(i, PPb)[1]:
od:
f := factor(f);

# Computing f is actually not necessary. We can compute the factors
# directly.
M[1] := MinPol(h1);
M[2] := MinPol(h2);

# first factorization
Ra := PolRem(AddLead(P), M[1]);
h2a := Quotl(-Ra[1], Ra[2]);
h1a := -PolQuot(P, [-h2a])[1];

# second factorization
Rb := PolRem(AddLead(P), M[2]);
h2b := Quotl(-Rb[1], Rb[2]);
h1b := -PolQuot(P, [-h2b])[1];

# check results
PolMult([-h1a], [-h2a]) - P;
PolMult([-h1b], [-h2b]) - P;
PlugIn(P, h2a);
PlugIn(M[1], h2a);
PlugIn(P, h2b);
PlugIn(M[2], h2b);