Sun Aug 18 18:21:38 2002 run on Windows
*** 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: 62953 ms, plus GC time: 1061 ms