Artifact 292d3bcd45fb86405ed0d6f7d9c9a0f9a5480084bc12670e0ae1e4946b24a8d7:
- Executable file
r38/log/normform.rlg
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 37051) [annotate] [blame] [check-ins using] [more...]
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