# Please do not post homework
# Matthew Esaia, 2/9/2025, Homework 4

with(linalg):

# 7.1 Vectors, arrays, and matrices
print(`Section 7.1`);

A := matrix(2, 3, [a, b, c, d, e, f]):

f := (i,j) -> u^(i * j);
A := matrix(4,4,f);
factor(det(A));

# 7.2 Matrix operations
print(`Section 7.2`);

A := matrix(2, 2, [1, 2, 3, 4]):
B := matrix(2, 2, [-2, 3, -5, 1]):
A + B;
evalm(%);

# 7.3 Elementary row operations
print(`Section 7.3`);

A := matrix(3,4, [1,1,3,-3,5,5,13,-7,3,1,7,-11]):
A1 := addrow(A, 1, 2, -5):
A2 := addrow(A1, 1, 3, -3):
A3 := mulrow(A2, 3, -1/2):
A4 := swaprow(A3, 2, 3):
A5 := mulrow(A4, 3, -1/2):
print(A5);

# 7.4 Gaussian elimination
print(`Section 7.4`);
gaussjord(A);

# 7.5 Inverses and determinants
print(`Section 7.5`);

A := matrix(3, 3, [1,1,3,5,5,13,3,1,7]):
det(A);
B := inverse(A);
evalm(B &* A);

# 7.6 Row space, column space, nullspace
# rowspace(A)
# colspace(A)
# nullspace(A)

# 7.7 Eigenvectors and diagonalization
# eigenvals(A)
# eigenvects(A)

# 7.8 Jordan form
print(`Section 7.8`);

C := matrix(4,4,[10,10,-14,15,0,3,0,0,8,1,-13,8,1,-8,-2,-4]):
jordan(C,'Q');
evalm(1/Q&*C&*Q);

# 7.9 Random matrices
# randmatrix(m,n)  // Entries -99 to 99


# 8.1 Conditional statements
print(`Section 8.1`);

x := 1:
if x > 0 then
	y := x + 1
else
	y := x - 1
fi:
y;

# 8.2 Loops
print(`Section 8.2`);

x := 0:
for i from 1 by 3 to 10 do
	x := x + i:
od:
x;

# 8.3 Procedures
print(`Section 8.3`);
print(`f(x)`);

f := proc(x) local z:

	if x >= 0 and x <= 1 then
		z := x^2:
	else
		z := 1 - x:
	fi:
	RETURN(z):
end:

plot(f, -2..2);

## PROBLEM 3

# CT(A): the conjugate transpose of A
CT:=proc(A) local n,i,j:
	n:=nops(A):
	[seq([seq(conjugate(A[j][i]), j=1..n)], i=1..n)]:
end:

# IsUM(A): Is the matrix A unitary
IsUM:=proc(A) local P,i,j,n:
	n:=nops(A);
	P:=multiply(A, CT(A)):
	evalb([seq([seq(P[i,j], j=1..n)], i=1..n)]=[seq([0$(i-1),1,0$(n-i)],i=1..n)]):
end:

# UM() determines if we have a unitary matrix with random alpha, beta, gamma, delta
UM:=proc(alpha,beta,gamma,delta) local A, A1, A2, A3:
	A1 := [[exp(-I * beta / 2), 0], [0, exp(I * beta / 2)]]:
	A2 := [[cos(gamma/2), -sin(gamma/2)], [sin(gamma/2), cos(gamma/2)]]:
	A3 := [[exp(-I * delta / 2), 0], [0, exp(I * delta / 2)]]:
	
	A := exp(I * alpha) * multiply(A1, multiply(A2,A3)):
	
	IsUM(A):
end:

TestUM := proc() local ra:
	ra:=rand(-100..100):

	UM(ra(), ra(), ra(), ra());
end: