reduce-historical/r38/log/normform.rlg

Tue Apr 15 00:34:38 2008 run on win32

*** co already defined as operator 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                      %
% Examples of calculations of matrix normal forms using the REDUCE     %
% NORMFORM package.                                                    %
%                                                                      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% load_package normform;

on errcont;

 % So that computation continues after an error.

%
% If using xr, the X interface for REDUCE, then turn on looking_good to 
% improve the appearance of the output.
%
fluid '(options!*);



lisp if memq('fmprint ,options!*) then on looking_good;



procedure test(tmp,A);
  % 
  % Checks that P * N * P^-1 = A where tmp is the output {P,N,P^-1} 
  % of the Normal form calculation on A.
  %
  begin
    if second tmp * first tmp * third tmp = A then
    write "Seems O.K." else rederr "something isn't working.";
  end;


test



%%%%%%%%%%%%%%%%%%%%%%%%%%%% Smithex %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  A := mat((3*x,x^2+x),(0,x^2));


     [3*x  x*(x + 1)]
     [              ]
a := [         2    ]
     [ 0      x     ]


  answer := smithex(A,x);


answer := {

           [x  0 ]
           [     ]
           [    2]
           [0  x ]

           ,


           [1  0]
           [    ]
           [x  1]

           ,


           [3   x + 1]
           [         ]
           [-3   - x ]

           }

  test(answer,A);


Seems O.K.


  %
  % Extend algebraic field to include sqrt2.
  %
  load_package arnum;


  defpoly sqrt2**2-2;


  A := mat((sqrt2*y^2,y+1),(3*sqrt2,y^3+y*sqrt2));


     [       2                ]
     [sqrt2*y       y + 1     ]
a := [                        ]
     [              2         ]
     [3*sqrt2   y*(y  + sqrt2)]


  answer := smithex(A,y);


answer := {

           [1             0           ]
           [                          ]
           [    5          3          ]
           [0  y  + sqrt2*y  - 3*y - 3]

           ,


           [       2   1       ]
           [sqrt2*y   ---*sqrt2]
           [           6       ]
           [                   ]
           [3*sqrt2       0    ]

           ,


           [    1            2         ]
           [1  ---*sqrt2*y*(y  + sqrt2)]
           [    6                      ]
           [                           ]
           [0           - sqrt2        ]

           }

  test(answer,A);


Seems O.K.

  off arnum;



  %
  % smithex will compute the Smith normal form of matrices containing 
  % only integer entries but the integers are regarded as univariate 
  % polynomials in x over a field F (the rationals unless the field has 
  % been extended). For calculations over the integers use smithex_int.
  %
  A := mat((9,-36,30),(-36,192,-180),(30,-180,180));


     [ 9   -36    30 ]
     [               ]
a := [-36  192   -180]
     [               ]
     [30   -180  180 ]


  answer := smithex(A,x);


*** WARNING: all matrix entries are integers.
    If calculations in Z(the integers) are required, use smithex_int.

answer := {

           [1  0  0]
           [       ]
           [0  1  0]
           [       ]
           [0  0  1]

           ,


           [            1  ]
           [ 9   18   -----]
           [           720 ]
           [               ]
           [-36  -24    0  ]
           [               ]
           [30    0     0  ]

           ,


           [1  -6    6   ]
           [             ]
           [         - 3 ]
           [0  1   ------]
           [         2   ]
           [             ]
           [0  0    2160 ]

           }

  test(answer,A);


Seems O.K.


%%%%%%%%%%%%%%%%%%%%%%%%%%%% Smithex_int %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  A := mat((1,2,3),(4,5,6),(7,8,x));


     [1  2  3]
     [       ]
a := [4  5  6]
     [       ]
     [7  8  x]


  answer := smithex_int(A);


***** ERROR: matrix contains non_integer entries. Try smithex.  


  A := mat((9,-36,30),(-36,192,-180),(30,-180,180));


     [ 9   -36    30 ]
     [               ]
a := [-36  192   -180]
     [               ]
     [30   -180  180 ]


  answer := smithex_int(A);


answer := {

           [3  0   0 ]
           [         ]
           [0  12  0 ]
           [         ]
           [0  0   60]

           ,


           [-17  -5   -4 ]
           [             ]
           [64   19   15 ]
           [             ]
           [-50  -15  -12]

           ,


           [1   -24  30 ]
           [            ]
           [-1  25   -30]
           [            ]
           [0   -1    1 ]

           }

  test(answer,A);


Seems O.K.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Frobenius %%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  A := mat(((y+y^2-x^2)/y,(x-x^2-y+y^2)/y,(x^2-y^2)/y),(1+x+y,
            (x-y+y^2+x*y)/y,-x-y),((y-x+y^2-x^2)/y,(x-y+y^2-x^2)/y,
            (x+x^2-y^2)/y));


     [       2    2              2        2           2    2   ]
     [    - x  + y  + y       - x  + x + y  - y      x  - y    ]
     [  ----------------    --------------------    ---------  ]
     [         y                     y                  y      ]
     [                                                         ]
     [                                   2                     ]
     [                        x*y + x + y  - y                 ]
a := [     x + y + 1         ------------------     - (x + y)  ]
     [                               y                         ]
     [                                                         ]
     [     2        2            2        2         2        2 ]
     [  - x  - x + y  + y     - x  + x + y  - y    x  + x - y  ]
     [--------------------  --------------------  -------------]
     [         y                     y                  y      ]


  answer := frobenius(A);


answer := {

           [ x                    ]
           [---  0        0       ]
           [ y                    ]
           [                      ]
           [          - x*(x + y) ]
           [ 0   0  --------------]
           [              y       ]
           [                      ]
           [                    2 ]
           [         x*y + x + y  ]
           [ 0   1  --------------]
           [              y       ]

           ,


                     3      2      2      2            2           2    2
                  - x  - 2*x *y - x  - x*y  + x*y + 2*y  + y      x  - y  - y
           mat((---------------------------------------------,-1,-------------),
                                y*(x + y + 1)                          y

               (x + y + 1,0, - (x + y + 1)),

                     2        2            2        2
                  - x  - x + y  + 2*y     x  + x - y  - y
               (----------------------,0,-----------------))
                          y                      y

,


[                     x - y                                        ]
[0                   -------                            1          ]
[                       y                                          ]
[                                                                  ]
[         3    2      2      2    3    2            2            2 ]
[      - x  - x *y - x  + x*y  + y  + y  + y     - x  - 2*x*y - y  ]
[-1  ----------------------------------------  --------------------]
[                 y*(x + y + 1)                     x + y + 1      ]
[                                                                  ]
[                2        2                                        ]
[               x  + x - y  - 2*y                                  ]
[0             -------------------                      1          ]
[                 y*(x + y + 1)                                    ]

}

  test(answer,A);


Seems O.K.


  %
  % Extend algebraic field to include i.
  %
%   load_package arnum;
    defpoly i^2+1;


    A := mat((-3-i,1,2+i,7-9*i),(-2,1,1,5-i),(-2-2*i,1,2+2*i,4-2*i),
             (2,0,-1,-2+8*i));


     [  - (i + 3)   1   i + 2    - (9*i - 7)]
     [                                      ]
     [     -2       1     1       - (i - 5) ]
a := [                                      ]
     [ - (2*i + 2)  1  2*i + 2   - (2*i - 4)]
     [                                      ]
     [     2        0    -1       8*i - 2   ]


    answer := frobenius(A);


answer := {

           [i + 1  0  0       0      ]
           [                         ]
           [  0    0  0    7*i - 3   ]
           [                         ]
           [  0    1  0   - (8*i - 9)]
           [                         ]
           [  0    0  1    8*i - 3   ]

           ,


           [ 425       189                             ]
           [-----*i + -----  -1   i + 3     18*i - 18  ]
           [ 106       106                             ]
           [                                           ]
           [ 634       258                             ]
           [-----*i + -----  0      2       2*i - 12   ]
           [ 53        53                              ]
           [                                           ]
           [ 150       588                             ]
           [-----*i - -----  0   2*i + 2    4*i - 10   ]
           [ 53        53                              ]
           [                                           ]
           [ 108       7                               ]
           [-----*i + ----   0     -2      - (16*i - 8)]
           [ 53        53                              ]

           ,


           mat((0, - i,1,1),

                        143       268    263       152   491       155
               (-1, - (-----*i - -----),-----*i + -----,-----*i + -----),
                        53        53     53        53    106       106

                       339       368        392       383        370       189
               (0, - (-----*i + -----), - (-----*i - -----), - (-----*i - -----)
                       106       53         53        106        53        53

                ),

                       101        9          7        54
               (0, - (-----*i + -----), - (-----*i - ----),1))
                       106       106        106       53

}

    off arnum;



  A := mat((10,-5,-5,8,3,0),(-4,2,-10,-7,-5,-5),(-8,2,7,3,7,5),
           (-6,-7,-7,-7,10,7),(-4,-3,-3,-6,8,-9),(-2,5,-5,9,7,-4));


     [10  -5  -5   8   3   0 ]
     [                       ]
     [-4  2   -10  -7  -5  -5]
     [                       ]
     [-8  2    7   3   7   5 ]
a := [                       ]
     [-6  -7  -7   -7  10  7 ]
     [                       ]
     [-4  -3  -3   -6  8   -9]
     [                       ]
     [-2  5   -5   9   7   -4]


  F := first frobenius(A);


     [0  0  0  0  0  -867960]
     [                      ]
     [1  0  0  0  0  -466370]
     [                      ]
     [0  1  0  0  0   47845 ]
f := [                      ]
     [0  0  1  0  0   -712  ]
     [                      ]
     [0  0  0  1  0    -95  ]
     [                      ]
     [0  0  0  0  1    16   ]



  %
  % Calculate in Z\23Z...
  %
    on modular;


    setmod 23;


1

    F_mod := first frobenius(A);


         [0  17  0  0  0  0 ]
         [                  ]
         [1  19  0  0  0  0 ]
         [                  ]
         [0  0   0  0  0  10]
f_mod := [                  ]
         [0  0   1  0  0  5 ]
         [                  ]
         [0  0   0  1  0  15]
         [                  ]
         [0  0   0  0  1  20]


 
  %
  % ...and with a balanced modular representation.
  %
    on balanced_mod;


    F_bal_mod := first frobenius(A);


             [0   - 6  0  0  0   0  ]
             [                      ]
             [1   - 4  0  0  0   0  ]
             [                      ]
             [0   0    0  0  0   10 ]
f_bal_mod := [                      ]
             [0   0    1  0  0   5  ]
             [                      ]
             [0   0    0  1  0   - 8]
             [                      ]
             [0   0    0  0  1   - 3]


    off balanced_mod;


    off modular;




%%%%%%%%%%%%%%%%%%%%%%%%%%% Ratjordan %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  A := mat(((y+y^2-x^2)/y,(x-x^2-y+y^2)/y,(x^2-y^2)/y),(1+x+y,
            (x-y+y^2+x*y)/y,-x-y),((y-x+y^2-x^2)/y,(x-y+y^2-x^2)/y,
            (x+x^2-y^2)/y));


     [       2    2              2        2           2    2   ]
     [    - x  + y  + y       - x  + x + y  - y      x  - y    ]
     [  ----------------    --------------------    ---------  ]
     [         y                     y                  y      ]
     [                                                         ]
     [                                   2                     ]
     [                        x*y + x + y  - y                 ]
a := [     x + y + 1         ------------------     - (x + y)  ]
     [                               y                         ]
     [                                                         ]
     [     2        2            2        2         2        2 ]
     [  - x  - x + y  + y     - x  + x + y  - y    x  + x - y  ]
     [--------------------  --------------------  -------------]
     [         y                     y                  y      ]


  answer := ratjordan(A);


answer := {

           [ x             ]
           [---   0     0  ]
           [ y             ]
           [               ]
           [      x        ]
           [ 0   ---    0  ]
           [      y        ]
           [               ]
           [ 0    0   x + y]

           ,


                     3      2      2      2            2           2
                  - x  - 2*x *y - x  - x*y  + x*y + 2*y  + y    - x  - x*y + y
           mat((---------------------------------------------,-----------------,
                                y*(x + y + 1)                              2
                                                                x*y - x + y

                  2        2
                 x  + x - y  - y
                -----------------),
                             2
                  x*y - x + y

                           y*(x + y + 1)    - y*(x + y + 1)
               (x + y + 1,---------------,------------------),
                                      2                 2
                           x*y - x + y       x*y - x + y

                     2        2             2        2        2        2
                  - x  - x + y  + 2*y    - x  - x + y  + y   x  + x - y  - y
               (----------------------,--------------------,-----------------))
                          y                           2                  2
                                           x*y - x + y        x*y - x + y

,


        x - y
mat((0,-------,1),
          y

             3      3    2  2    2      2      3      2            4    3    2
          - x *y + x  - x *y  - x *y + x  + x*y  - x*y  - 2*x*y + y  + y  + y
    (-1,-----------------------------------------------------------------------,
                                     2
                                    y *(x + y + 1)

          2      2        2              3
       - x *y + x  - 2*x*y  + x*y + x - y
     --------------------------------------),
                 y*(x + y + 1)

          - x - y + 1     x + y
    (-1,--------------,-----------))
          x + y + 1     x + y + 1

}

  test(answer,A);


Seems O.K.


  %
  % Extend algebraic field to include sqrt(2).
  %
  %  load_package arnum;
    defpoly sqrt2**2-2;


    A:= mat((4*sqrt2-6,-4*sqrt2+7,-3*sqrt2+6),(3*sqrt2-6,-3*sqrt2+7,
            -3*sqrt2+6),(3*sqrt2,1-3sqrt2,-2*sqrt2));


     [4*sqrt2 - 6   - (4*sqrt2 - 7)   - (3*sqrt2 - 6)]
     [                                               ]
a := [3*sqrt2 - 6   - (3*sqrt2 - 7)   - (3*sqrt2 - 6)]
     [                                               ]
     [  3*sqrt2     - (3*sqrt2 - 1)      - 2*sqrt2   ]


    answer := ratjordan(A);


answer := {

           [sqrt2    0           0        ]
           [                              ]
           [  0    sqrt2         0        ]
           [                              ]
           [  0      0     - (3*sqrt2 - 1)]

           ,


           [                21           49           21           18  ]
           [7*sqrt2 - 6    ----*sqrt2 - ----      - (----*sqrt2 - ----)]
           [                31           31           31           31  ]
           [                                                           ]
           [                21           18           21           18  ]
           [3*sqrt2 - 6    ----*sqrt2 - ----      - (----*sqrt2 - ----)]
           [                31           31           31           31  ]
           [                                                           ]
           [                  3            24       3            24    ]
           [3*sqrt2 + 1   - (----*sqrt2 + ----)    ----*sqrt2 + ----   ]
           [                  31           31       31           31    ]

           ,


           [0       sqrt2 + 1          1   ]
           [                               ]
           [-1     4*sqrt2 + 9      4*sqrt2]
           [                               ]
           [         1                     ]
           [-1   - (---*sqrt2 - 1)     1   ]
           [         6                     ]

           }

    test(answer,A);


Seems O.K.

    off arnum;



  A := mat((-12752,-6285,-9457,-7065,-4939,-5865,-3769),(13028,6430,
           9656, 7213,5041,5984,3841),(16425,8080,12192,9108,6370,7569,
           4871), (-6065,-2979,-4508,-3364,-2354,-2801,-1803),(2968,
           1424,2231, 1664,1171,1404,919),(-22762,-11189,-16902,-12627,
           -8833, -10498,-6760),(23112,11400,17135,12799,8946,10622,
           6821));


     [-12752  -6285   -9457   -7065   -4939  -5865   -3769]
     [                                                    ]
     [13028    6430    9656    7213   5041    5984   3841 ]
     [                                                    ]
     [16425    8080   12192    9108   6370    7569   4871 ]
     [                                                    ]
a := [-6065   -2979   -4508   -3364   -2354  -2801   -1803]
     [                                                    ]
     [ 2968    1424    2231    1664   1171    1404    919 ]
     [                                                    ]
     [-22762  -11189  -16902  -12627  -8833  -10498  -6760]
     [                                                    ]
     [23112   11400   17135   12799   8946   10622   6821 ]


  R := first ratjordan(A);


     [0  2  0  0  0  0  0 ]
     [                    ]
     [1  0  0  0  0  0  0 ]
     [                    ]
     [0  0  0  0  0  0  5 ]
     [                    ]
r := [0  0  1  0  0  0  0 ]
     [                    ]
     [0  0  0  1  0  0  -2]
     [                    ]
     [0  0  0  0  1  0  3 ]
     [                    ]
     [0  0  0  0  0  1  0 ]



  %
  % Calculate in Z/23Z...
  %
    on modular;


    setmod 23;


23

    R_mod := first ratjordan(A);


         [19  0   0   0  0  0  0 ]
         [                       ]
         [0   18  0   0  0  0  0 ]
         [                       ]
         [0   0   17  0  0  0  0 ]
         [                       ]
r_mod := [0   0   0   5  0  0  0 ]
         [                       ]
         [0   0   0   0  0  0  5 ]
         [                       ]
         [0   0   0   0  1  0  19]
         [                       ]
         [0   0   0   0  0  1  10]



  %
  % ...and with a balanced modular representation.
  %
    on balanced_mod;


    R_bal_mod := first ratjordan(A);


             [5   0     0     0    0  0   0  ]
             [                               ]
             [0   - 4   0     0    0  0   0  ]
             [                               ]
             [0   0     - 5   0    0  0   0  ]
             [                               ]
r_bal_mod := [0   0     0     - 6  0  0   0  ]
             [                               ]
             [0   0     0     0    0  0   5  ]
             [                               ]
             [0   0     0     0    1  0   - 4]
             [                               ]
             [0   0     0     0    0  1   10 ]


    off balanced_mod;


    off modular;




%%%%%%%%%%%%%%%%%%%%%%%%%%% jordansymbolic %%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  A := mat(((y+y^2-x^2)/y,(x-x^2-y+y^2)/y,(x^2-y^2)/y),(1+x+y,
            (x-y+y^2+x*y)/y,-x-y),((y-x+y^2-x^2)/y,(x-y+y^2-x^2)/y,
            (x+x^2-y^2)/y));


     [       2    2              2        2           2    2   ]
     [    - x  + y  + y       - x  + x + y  - y      x  - y    ]
     [  ----------------    --------------------    ---------  ]
     [         y                     y                  y      ]
     [                                                         ]
     [                                   2                     ]
     [                        x*y + x + y  - y                 ]
a := [     x + y + 1         ------------------     - (x + y)  ]
     [                               y                         ]
     [                                                         ]
     [     2        2            2        2         2        2 ]
     [  - x  - x + y  + y     - x  + x + y  - y    x  + x - y  ]
     [--------------------  --------------------  -------------]
     [         y                     y                  y      ]


  answer := jordansymbolic(A);


answer := {

           [ x             ]
           [---   0     0  ]
           [ y             ]
           [               ]
           [      x        ]
           [ 0   ---    0  ]
           [      y        ]
           [               ]
           [ 0    0   x + y]

           ,

              lambda*y - x
           {{--------------,lambda - x - y},
                   y

            lambda},


                     3      2      2      2            2           2
                  - x  - 2*x *y - x  - x*y  + x*y + 2*y  + y    - x  - x*y + y
           mat((---------------------------------------------,-----------------,
                                y*(x + y + 1)                              2
                                                                x*y - x + y

                  2        2
                 x  + x - y  - y
                -----------------),
                             2
                  x*y - x + y

                           y*(x + y + 1)    - y*(x + y + 1)
               (x + y + 1,---------------,------------------),
                                      2                 2
                           x*y - x + y       x*y - x + y

                     2        2             2        2        2        2
                  - x  - x + y  + 2*y    - x  - x + y  + y   x  + x - y  - y
               (----------------------,--------------------,-----------------))
                          y                           2                  2
                                           x*y - x + y        x*y - x + y

,


        x - y
mat((0,-------,1),
          y

             3      3    2  2    2      2      3      2            4    3    2
          - x *y + x  - x *y  - x *y + x  + x*y  - x*y  - 2*x*y + y  + y  + y
    (-1,-----------------------------------------------------------------------,
                                     2
                                    y *(x + y + 1)

          2      2        2              3
       - x *y + x  - 2*x*y  + x*y + x - y
     --------------------------------------),
                 y*(x + y + 1)

          - x - y + 1     x + y
    (-1,--------------,-----------))
          x + y + 1     x + y + 1

}



  %
  % Extend algebraic field.
  %
  % load_package arnum;
    defpoly b^3-2*b+b-5;


    A := mat((1-b,2+b^2),(3+b-2*b^2,3));


     [                    2    ]
     [    - (b - 1)      b  + 2]
a := [                         ]
     [       2                 ]
     [ - (2*b  - b - 3)    3   ]


    answer := jordansymbolic(A);


answer := {

           [lambda11     0    ]
           [                  ]
           [   0      lambda12]

           ,

                   2                       2
           {{lambda  + (b - 4)*lambda + 3*b  + 4*b - 8},lambda},


           [  lambda11 - 3       lambda12 - 3   ]
           [                                    ]
           [       2                  2         ]
           [ - (2*b  - b - 3)   - (2*b  - b - 3)]

           ,


                      1966    2     3514         1054                 1
           mat(( - (--------*b  + --------*b - --------)*(lambda11 + ---*b - 2),
                     239891        239891       239891                2

                   127472    2     236383         82923
                (----------*b  + ----------*b + ---------)
                  29986375        29986375       5997275

                              26   2    107       45
                *(lambda11 + ----*b  - -----*b + ----)),
                              11        11        11

                      1966    2     3514         1054                 1
               ( - (--------*b  + --------*b - --------)*(lambda12 + ---*b - 2),
                     239891        239891       239891                2

                   127472    2     236383         82923
                (----------*b  + ----------*b + ---------)
                  29986375        29986375       5997275

                              26   2    107       45
                *(lambda12 + ----*b  - -----*b + ----)))
                              11        11        11

}

    off arnum;



  A := mat((-9,21,-15,4,2,0),(-10,21,-14,4,2,0),(-8,16,-11,4,2,0),
           (-6,12,-9,3,3,0),(-4,8,-6,0,5,0),(-2,4,-3,0,1,3));


     [-9   21  -15  4  2  0]
     [                     ]
     [-10  21  -14  4  2  0]
     [                     ]
     [-8   16  -11  4  2  0]
a := [                     ]
     [-6   12  -9   3  3  0]
     [                     ]
     [-4   8   -6   0  5  0]
     [                     ]
     [-2   4   -3   0  1  3]


  answer := jordansymbolic(A);


answer := {

           [3  0  0  0     0         0    ]
           [                              ]
           [0  3  0  0     0         0    ]
           [                              ]
           [0  0  1  1     0         0    ]
           [                              ]
           [0  0  0  1     0         0    ]
           [                              ]
           [0  0  0  0  lambda31     0    ]
           [                              ]
           [0  0  0  0     0      lambda32]

           ,

                                         2
           {{lambda - 3,lambda - 1,lambda  - 4*lambda + 5},lambda},


           [     - 3         1    6*lambda31 - 17     6*lambda32 - 17  ]
           [3  ------   1   ---  -----------------   ----------------- ]
           [     8           4           2                   2         ]
           [                                                           ]
           [     - 3         1    5*(lambda31 - 3)    5*(lambda32 - 3) ]
           [3  ------   1   ---  ------------------  ------------------]
           [     8           4           2                   2         ]
           [                                                           ]
           [     - 3         1                                         ]
           [3  ------   1   ---   2*(lambda31 - 3)    2*(lambda32 - 3) ]
           [     8           4                                         ]
           [                                                           ]
           [     - 3    3    3    3*(lambda31 - 3)    3*(lambda32 - 3) ]
           [3  ------  ---  ---  ------------------  ------------------]
           [     8      4    8           2                   2         ]
           [                                                           ]
           [     - 3    1    1                                         ]
           [3  ------  ---  ---     lambda31 - 3        lambda32 - 3   ]
           [     8      2    4                                         ]
           [                                                           ]
           [     - 1    1    1      lambda31 - 3        lambda32 - 3   ]
           [2  ------  ---  ---    --------------      --------------  ]
           [     8      4    8           2                   2         ]

           ,


           [                                     - 1        ]
           [       0              0        0   ------  0   1]
           [                                     3          ]
           [                                                ]
           [                                     8          ]
           [       0              0        0    ---    -8  8]
           [                                     3          ]
           [                                                ]
           [       0              -4       6     0     -2  0]
           [                                                ]
           [       0              0        -4    8     -4  0]
           [                                                ]
           [ - lambda31 + 3  lambda31 - 4  1     0     0   0]
           [                                                ]
           [ - lambda32 + 3  lambda32 - 4  1     0     0   0]

           }

  
  % Check to see if looking_good (*) is on as the choice of using 
  % either lambda or xi is dependent upon this.
  % (* -> the use of looking_good is described in the manual.).
  if not lisp !*looking_good then
  <<
    %
    % NB: we use lambda_ in solve (instead of lambda) as lambda is used
    % for other purposes in REDUCE which mean it cannot be used with
    % solve.
    %
    solve(lambda_^2-4*lambda_+5,lambda_);
    J := sub({lambda31=i + 2,lambda32= - i + 2},first answer);
    P := sub({lambda31=i + 2,lambda32= - i + 2},third answer);
    Pinv :=sub({lambda31=i + 2,lambda32= - i + 2},third rest answer);
  >>
  else 
  <<
    solve(xi^2-4*xi+5,xi);
    J := sub({xi(3,1)=i + 2,xi(3,2)= - i + 2},first answer);
    P := sub({xi(3,1)=i + 2,xi(3,2)= - i + 2},third answer);
    Pinv := sub({xi(3,1)=i + 2,xi(3,2)= - i + 2},third rest answer);
  >>;


  test({J,P,Pinv},A);


Seems O.K.


  %
  % Calculate in Z/23Z...
  %
    on modular;


    setmod 23;


23

    answer := jordansymbolic(A)$


    J_mod := {first answer, second answer};


j_mod := {

          [3  0  0  0     0         0    ]
          [                              ]
          [0  3  0  0     0         0    ]
          [                              ]
          [0  0  1  1     0         0    ]
          [                              ]
          [0  0  0  1     0         0    ]
          [                              ]
          [0  0  0  0  lambda31     0    ]
          [                              ]
          [0  0  0  0     0      lambda32]

          ,

                                          2
          {{lambda + 20,lambda + 22,lambda  + 19*lambda + 5},lambda}}


  %
  % ...and with a balanced modular representation.
  %
    on balanced_mod;


    answer := jordansymbolic(A)$


    J_bal_mod := {first answer, second answer};


j_bal_mod := {

              [3  0  0  0     0         0    ]
              [                              ]
              [0  3  0  0     0         0    ]
              [                              ]
              [0  0  1  1     0         0    ]
              [                              ]
              [0  0  0  1     0         0    ]
              [                              ]
              [0  0  0  0  lambda31     0    ]
              [                              ]
              [0  0  0  0     0      lambda32]

              ,

                                            2
              {{lambda - 3,lambda - 1,lambda  - 4*lambda + 5},lambda}}

    off balanced_mod;


    off modular;


 

%%%%%%%%%%%%%%%%%%%%%%%%%%%% jordan %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  A := mat((1,y),(y^2,3));


     [1   y]
     [     ]
a := [ 2   ]
     [y   3]


  answer := jordan(A);


answer := {

           [      3                              ]
           [sqrt(y  + 1) + 2           0         ]
           [                                     ]
           [                           3         ]
           [       0           - sqrt(y  + 1) + 2]

           ,


           [      3                     3          ]
           [sqrt(y  + 1) - 1   - (sqrt(y  + 1) + 1)]
           [                                       ]
           [        2                   2          ]
           [       y                   y           ]

           ,


           [        3                  3         3       ]
           [  sqrt(y  + 1)       sqrt(y  + 1) + y  + 1   ]
           [ --------------     -----------------------  ]
           [       3                    2   3            ]
           [   2*(y  + 1)            2*y *(y  + 1)       ]
           [                                             ]
           [          3                  3         3     ]
           [  - sqrt(y  + 1)     - sqrt(y  + 1) + y  + 1 ]
           [-----------------  --------------------------]
           [       3                    2   3            ]
           [   2*(y  + 1)            2*y *(y  + 1)       ]

           }

  test(answer,A);


Seems O.K.


  A := mat((-12752,-6285,-9457,-7065,-4939,-5865,-3769),(13028,6430,
          9656, 7213,5041,5984,3841),(16425,8080,12192,9108,6370,7569,
          4871), (-6065,-2979,-4508,-3364,-2354,-2801,-1803),(2968,
          1424,2231, 1664,1171,1404,919),(-22762,-11189,-16902,-12627,
          -8833, -10498,-6760),(23112,11400,17135,12799,8946,10622,
          6821));


     [-12752  -6285   -9457   -7065   -4939  -5865   -3769]
     [                                                    ]
     [13028    6430    9656    7213   5041    5984   3841 ]
     [                                                    ]
     [16425    8080   12192    9108   6370    7569   4871 ]
     [                                                    ]
a := [-6065   -2979   -4508   -3364   -2354  -2801   -1803]
     [                                                    ]
     [ 2968    1424    2231    1664   1171    1404    919 ]
     [                                                    ]
     [-22762  -11189  -16902  -12627  -8833  -10498  -6760]
     [                                                    ]
     [23112   11400   17135   12799   8946   10622   6821 ]


  on rounded;


  J := first jordan(A);


*** Domain mode rounded changed to rational 

*** Domain mode rational changed to complex-rational 

*** Domain mode complex-rational changed to rational 

*** Domain mode rational changed to rounded 

j := mat((1.41421356237,0,0,0,0,0,0),

         (0, - 1.41421356237,0,0,0,0,0),

         (0,0, - 1.80492,0,0,0,0),

         (0,0,0, - 1.12491,0,0,0),

         (0,0,0,0,1.03589*i + 0.620319,0,0),

         (0,0,0,0,0, - 1.03589*i + 0.620319,0),

         (0,0,0,0,0,0,1.68919))


  off rounded;



  %
  % Extend algebraic field.
  %
  % load_package arnum;
    defpoly b^3-2*b+b-5;


    A := mat((1-b,2+b^2),(3+b-2*b^2,3));


     [                    2    ]
     [    - (b - 1)      b  + 2]
a := [                         ]
     [       2                 ]
     [ - (2*b  - b - 3)    3   ]


    J := first jordan(A);


           1            2
j := mat((---*(sqrt(11*b  + 24*b - 48)*i - (b - 4)),0),
           2

                1            2
         (0, - ---*(sqrt(11*b  + 24*b - 48)*i + b - 4)))
                2


    off arnum;



end;


Time for test: 1919 ms, plus GC time: 109 ms