> restart;
# The following examples is Jan Verschelde's  examples:
# The purpose of this worksheet is to generate instances of multivariate
# polynomials with complex coefficients.  Executing this worksheet
# generates the files exf1, exf2, exf3, and exf4.
# 
# The Jacobian for a Stewart-Gough Platform
# We derive the Jacobian for a platform, first symbolically.
# 1. The Symbolic Jacobian Matrix.
# A rotation is represented by a quaternion q, as an element in
# projective 3-space:
> q := array(0..3):
> Q := sum(q['i']^2,'i'=0..3):
> Rhat := Matrix(3,3):
> R := Rhat/Q:
> Rhat[1,1] := q[0]^2 + q[1]^2 - q[2]^2 - q[3]^2:
> Rhat[2,2] := q[0]^2 - q[1]^2 + q[2]^2 - q[3]^2:
> Rhat[3,3] := q[0]^2 - q[1]^2 - q[2]^2 + q[3]^2:
> Rhat[1,2] := 2*(q[1]*q[2] + q[0]*q[3]):
> Rhat[2,1] := 2*(q[1]*q[2] - q[0]*q[3]):
> Rhat[1,3] := 2*(q[1]*q[3] - q[0]*q[2]):# 
> Rhat[3,1] := 2*(q[1]*q[3] + q[0]*q[2]):
> Rhat[2,3] := 2*(q[2]*q[3] + q[0]*q[1]):
> Rhat[3,2] := 2*(q[2]*q[3] - q[0]*q[1]):
> R;
        [    2       2       2       2
        [q[0]  + q[1]  - q[2]  - q[3]  , 2 q[1] q[2] + 2 q[0] q[3] ,
                                 ]
        2 q[1] q[3] - 2 q[0] q[2]]
        [                                2       2       2       2
        [2 q[1] q[2] - 2 q[0] q[3] , q[0]  - q[1]  + q[2]  - q[3]  ,
                                 ]
        2 q[2] q[3] + 2 q[0] q[1]]
        [
        [2 q[1] q[3] + 2 q[0] q[2] , 2 q[2] q[3] - 2 q[0] q[1] ,
            2       2       2       2]   /
        q[0]  - q[1]  - q[2]  + q[3] ]  /  (
                                       /
            2       2       2       2
        q[0]  + q[1]  + q[2]  + q[3] )
# Next we define the position vector pt:
> pt := 
:
# and we define the locations at endplate and the baseplate joints:
> 
> for i from 1 to 6 do
>   aa||i := :
>   bb||i := :
> end do:
# Now we are ready to define the Jacobian matrix:
> 
> jp := Matrix(6,6):
> for i from 1 to 6 do
>   jpA||i :=
> LinearAlgebra:-ScalarMultiply(pt-bb||i,Q)+LinearAlgebra:-Multiply(Rhat
> ,aa||i);
>   for j from 1 to 3 do
>     jp[i,j] := jpA||i[j];
>   end do;
>   jpB||i :=
> LinearAlgebra:-CrossProduct(LinearAlgebra:-Multiply(Rhat,aa||i),pt-bb|
> |i);
>   for j from 1 to 3 do
>     jp[i,j+3] := jpB||i[j];
>   end do;
> end do:
# We better do not ask here to compute the determinant...
> 
# 2. Numerical Instances
# We have three vectors to specify: the position vector and the
# locations at endplate and the baseplate joints.
#  For the purposes of applying homotopy continuation to factor the
# polynomials, we generate random complex number of modulus one, using
# 20 digit
> 
# precision
> Digits:=20;
                             Digits := 20
> randcomp := proc()
>   description `returns a random complex number of modulus one`;
>   local angle;
> angle := stats[random,uniform[0,2*Pi]](1):
>   evalc(cos(angle)+sin(angle)*I):
> end proc:
> randcomp();
          -0.89780540910193948634 + 0.44039237889329905448 I
# Example 1
# In this first instance, we pick the locations and position vector at
# random:
> randomize:
> for i from 1 to 3 do
>   p[i] := randcomp(): 
> end do:
> for i from 1 to 6 do
>   for j from 1 to 3 do
>     a||i[j] := randcomp();
>     b||i[j] := randcomp();
>   end do:
> end do:
> f1 := LinearAlgebra:-Determinant(jp):
> save f1, "verschelde/svw_exf1";
# Example 2
# Here we set the third element of all locations equal to zero:
> Digits:=20;
                             Digits := 20
> for i from 1 to 6 do
>   a||i[3] := 0:
>   b||i[3] := 0:
> end do:
> f2 := LinearAlgebra:-Determinant(jp):
> save f2,"verschelde/svw_exf2";
# Example 3
# In addition we set q[1] and q[2] equal to zero:
> Digits:=20;
                             Digits := 20
> q[1] := 0: q[2] := 0:
> f3 := LinearAlgebra:-Determinant(jp):
> f3:
> ff:=f3:
> save ff, "verschelde/svw_exf3";
# Example 4
# Finally, we let the point pt have arbitrary coordinates:
> Digits:=20;
                             Digits := 20
> for i from 1 to 3 do
>   p[i] := evaln(p[i]);
> end do:
> pt := ;
                                   [p[1]]
                                   [    ]
                             pt := [p[2]]
                                   [    ]
                                   [p[3]]
> jp := Matrix(6,6);
> for i from 1 to 6 do
>   jpA||i :=
> LinearAlgebra:-ScalarMultiply(pt-bb||i,Q)+LinearAlgebra:-Multiply(Rhat
> ,aa||i);
>   for j from 1 to 3 do
>     jp[i,j] := jpA||i[j];
>   end do;
>   jpB||i :=
> LinearAlgebra:-CrossProduct(LinearAlgebra:-Multiply(Rhat,aa||i),pt-bb|
> |i);
>   for j from 1 to 3 do
>     jp[i,j+3] := jpB||i[j];
>   end do;
> end do:
                        [0    0    0    0    0    0]
                        [                          ]
                        [0    0    0    0    0    0]
                        [                          ]
                        [0    0    0    0    0    0]
                  jp := [                          ]
                        [0    0    0    0    0    0]
                        [                          ]
                        [0    0    0    0    0    0]
                        [                          ]
                        [0    0    0    0    0    0]
> f4 :=LinearAlgebra:- Determinant(jp):
> save f4, "verschelde/svw_exf4";
> 
> 
>