Define NR_Mat(A, L)
  Return Mat([[NR(Poly(X),L) | X In Riga] | Riga In A]);
EndDefine; -- NR_Mat

----------------------------------------------------------------------
-- esercizio 1 (prima matrice)
A := Mat([
	 [1,     2/5, -4/3],
	 [2/5, 29/25, -8/15],
	 [-4/3,-8/15, 25/9]]);  A;

--il minore principale di ordine 1 e' 1 > 0
--il minore principale di ordine 2 e'
Det(Submat(A, [1,2],[1,2])); -- 1 > 0
--il minore principale di ordine 3 e'
Det(A); -- 1 > 0
-- per il criterio di Sylvester, A e' definita positiva
-- quindi esiste la decomposizione di Cholesky.

-- riduzione gaussiana:
E1 := Identity(3);  E1[2,1] := -2/5; E1*A;
E2 := Identity(3);  E2[3,1] := 4/3;  E2*E1*A;
P := Transposed(E2*E1);  Transposed(P)*A*P;
D := Transposed(P)*A*P;  D;
-- quindi D e' una forma diagonale di A
--Mat([
--  [1, 0, 0],
--  [0, 1, 0],
--  [0, 0, 1]
--])

-- calcolo l'inversa di P
InvE1 := Identity(3);  InvE1[2,1]:=2/5;  InvE1*E1; --verificato InvE1*E1=I
InvE2 := Identity(3);  InvE2[3,1]:=-4/3; InvE2*E2; --verificato InvE2*E2=I
-- quindi
InvP := Transposed(InvE2)*Transposed(InvE1);  InvP*P; --verificato InvP*P=I
-- verifico la decomposizione di Cholesky parziale
A = Transposed(InvP)*D*InvP; --> TRUE

-- la matrice diagonale B tale che B^2=D e' B=I, quindi U = B*InvP = InvP
U := InvP;  U;
--   [1, 2/5, -4/3],
--   [0, 1, 0],
--   [0, 0, 1]

-- verifico la decomposizione di Cholesky
A = Transposed(U)*U;  -->  TRUE

-------------------------------
-- esercizio 1 (seconda matrice)
A := Mat([
	 [1,     2/5, -4/3],
	 [2/5, 29/25, -8/15],
	 [-4/3,-8/15,  -2]]);  A;

-- A e' non definita, infatti:
MvE := Transposed(Mat([[1,0,0]]));
MuE := Transposed(Mat([[0,0,1]]));
Transposed(MvE) * A * MvE; --> 1 > 0
Transposed(MuE) * A * MuE; --> -2 < 0
-- quindi non esiste la decomposizione di Cholesky


----------------------------------------------------------------------
-- esercizio 2
A := Mat([
	 [-2, -2,  1],
	 [-2,  1, -1],
	 [-1, -1, -1]]);  A;
Det(A);  --> 9
-- A ha rango massimo, 3, perche' ha un minore di ordine 3 non nullo

-- Metodo di Cholesky
AA := Transposed(A)*A;  AA;
--  [9, 3, 1],
--  [3, 6, -2],
--  [1, -2, 3]

-- A di rango massimo, allora AA e' simmetrica e definita positiva
-- inizio il calcolo della decomposizione di Cholesky:
E1 := Identity(3);  E1[2,1] := -1/3;
E2 := Identity(3);  E2[3,1] := -1/9;
P := Transposed(E2*E1);  Transposed(P)*AA*P;
E3 := Identity(3);  E3[3,2] := 7/15;
P := Transposed(E3*E2*E1);  Transposed(P)*AA*P;
D := Transposed(P)*AA*P;  D;
-- ottengo la matrice diagonale D:
--  [9, 0, 0],
--  [0, 5, 0],
--  [0, 0, 9/5]

-- calcolo l'inversa di P
InvE1 := Identity(3);  InvE1[2,1] := 1/3;   InvE1*E1;
InvE2 := Identity(3);  InvE2[3,1] := 1/9;   InvE2*E2;
InvE3 := Identity(3);  InvE3[3,2] := -7/15; InvE3*E3;
InvP  := Transposed(InvE1*InvE2*InvE3); InvP*P;
-- normalizzo le colonne di A*P: mi serve a = rad(5)
Use Q[a];
L := [a^2-5];
B := DiagonalMat([3,a,3a/5]);  B;
-- verifico che B^2 = D
NR_Mat(B*B, L) = D;  --> TRUE
-- calcolo l'inversa di B:
InvB := DiagonalMat([1/3, a/5, a/3]);
-- verifico
NR_Mat(InvB*B,L);

-- definisco Q e R con a=rad(5)
Q := A*P*InvB;  Q;
--   [-2/3, -4/15a, 1/5a],
--   [-2/3, 1/3a, 0],
--   [-1/3, -2/15a, -2/5a]
-- verifico Q ortonormale:
NR_Mat(Transposed(Q)*Q, L) = Identity(3);

R := B*InvP;  R;
--   [3, 1, 1/3],
--   [0, a, -7/15a],
--   [0, 0, 3/5a]
-- verifico la decomposizione QR
NR_Mat(Q*R,L) = A;  --> TRUE

----------------------------------------------------------------------
-- esercizio 3

A := Mat([ [1, 1, 0],
	   [0, 1,-1],
	   [4, 1, 3] ]);
B := Transposed(A);  B; -- vettori in riga
--  [1, 0, 4],
--  [1, 1, 1],
--  [0, -1, 3]
-- prima colonna
E1 := Identity(3); E1[2,1]:=-1; E1*B;
-- seconda colonna
E2 := Identity(3); E2[3,2]:=1; E2*E1*B;
--  [1, 0, 4],
--  [0, -1, 3],
--  [0, 0, 0]
-- quindi B ha rango 2

-- Conclusione:
-- (a) Il sottospazio generato dalle colonne di A ha dimensione 2.
-- (b) una base e' F = ((1,0,4), (0,-1,3))

----------------------------------------------------------------------
-- esercizio 4
A := Mat([
	 [0,3,1],
	 [3,3,4],
	 [1,4,2] ]);  A;
-- calcolo una forma diagonale
E0 := Mat([
	  [0,1,0],
	  [1,0,0],
	  [0,0,1]]);  E0;
P := Transposed(E0);  Transposed(P)*A*P;
E1 := Identity(3);  E1[2,1]:=-1;
P := Transposed(E1*E0);  Transposed(P)*A*P;
P := Transposed(E2*E1*E0);  Transposed(P)*A*P;
E3 := Identity(3);  E3[3,2]:=-1;
P := Transposed(E3*E2*E1*E0);  Transposed(P)*A*P;
MFE := P; MFE;
--   [0, 1, -1],
--   [1, -1, -1/3],
--   [0, 0, 1]
-- quindi la matrice della forma quadratica Q rispetto alla nuova base F e'
MQF := Transposed(MFE)*A*MFE;  MQF;
--   [3, 0, 0],
--   [0, -3, 0],
--   [0, 0, -1/3]
-- quindi Q e' non definita

-- u tale che Q(u) = 0
MuF := Transposed(Mat([[1,1,0]]));
-- verifico
Transposed(MuF) * MQF * MuF; --> 0
MuE := MFE * MuF;  MuE;
-- quindi u = (1,0,0)
-- verifico
Transposed(MuE) * A * MuE; --> 0

-- v tale che Q(v) = 0
MvF := Transposed(Mat([[0,0,1]]));
-- verifico
Transposed(MvF) * MQF * MvF; --> 0
MvE := MFE * MvF;  MvE;
-- quindi v = (-1,-1/3,1)
-- verifico
Transposed(MvE) * A * MvE; --> -1/3

-- Conclusione:
-- Q e' non definita
-- u = (1,0,0)
-- v = (-1,-1/3,1)