|\^/|     Maple 16 (X86 64 LINUX)
._|\|   |/|_. Copyright (c) Maplesoft, a division of Waterloo Maple Inc. 2012
 \  MAPLE  /  All rights reserved. Maple is a trademark of
 <____ ____>  Waterloo Maple Inc.
      |       Type ? for help.
# ======================================================================== #
#
# 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];
                     h1 := [0, 0, 1 - epsilon, 1 + epsilon]

> h2 := [1, -1 + e,     0, 1 + e];
                    h2 := [1, -1 + epsilon, 0, 1 + epsilon]

> 
# 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]);
P := [[-1, 1 - epsilon, -1 + epsilon, -2 - 2 epsilon],

    [-1 - 2 epsilon, 1, -epsilon, 2 - epsilon]]

> Pb := map(Cj, P);       # conjugate quaternion
Pb := [[-1, -1 + epsilon, 1 - epsilon, 2 + 2 epsilon],

    [-1 - 2 epsilon, -1, epsilon, -2 + epsilon]]

> PPb := PolMult(P, Pb);  # a real polynomial
       PPb := [[-2, 0, 0, 0], [5, 0, 0, 0], [-4, 0, 0, 0], [6, 0, 0, 0]]

> f := t^4:
> for i from 1 to nops(PPb) do
>   f := f + t^(4-i) * op(i, PPb)[1]:
> od:
> f := factor(f);
                                 2        2
                          f := (t  + 2) (t  - 2 t + 3)

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

> M[2] := MinPol(h2);
                     M[2] := [[-2, 0, 0, 0], [3, 0, 0, 0]]

> 
# first factorization
> Ra := PolRem(AddLead(P), M[1]);
Ra := [[-1, 1 - epsilon, -1 + epsilon, -2 - 2 epsilon],

    [-3 - 2 epsilon, 1, -epsilon, 2 - epsilon]]

> h2a := Quotl(-Ra[1], Ra[2]);
                         30 epsilon          3 epsilon        13 epsilon
      h2a := [0, - 4/7 + ----------, - 1/7 - ---------, 9/7 + ----------]
                             49                 49                49

> h1a := -PolQuot(P, [-h2a])[1];
                          19 epsilon        46 epsilon        85 epsilon
       h1a := [1, - 3/7 + ----------, 8/7 - ----------, 5/7 + ----------]
                              49                49                49

> 
# second factorization
> Rb := PolRem(AddLead(P), M[2]);
Rb := [[1, 1 - epsilon, -1 + epsilon, -2 - 2 epsilon],

    [-4 - 2 epsilon, 1, -epsilon, 2 - epsilon]]

> h2b := Quotl(-Rb[1], Rb[2]);
                    h2b := [1, -1 + epsilon, 0, 1 + epsilon]

> h1b := -PolQuot(P, [-h2b])[1];
                    h1b := [0, 0, 1 - epsilon, 1 + epsilon]

> 
# check results
> PolMult([-h1a], [-h2a]) - P;
                          [[0, 0, 0, 0], [0, 0, 0, 0]]

> PolMult([-h1b], [-h2b]) - P;
                          [[0, 0, 0, 0], [0, 0, 0, 0]]

> PlugIn(P, h2a);
                                  [0, 0, 0, 0]

> PlugIn(M[1], h2a);
                                  [0, 0, 0, 0]

> PlugIn(P, h2b);
                                  [0, 0, 0, 0]

> PlugIn(M[2], h2b);
                                  [0, 0, 0, 0]

> 
> quit
memory used=2.3MB, alloc=4.0MB, time=0.05