# OK to post homework
# James Betti, 9 Feb 2025, Assignment 5

RandomState := proc(n,K) local ra,i,A,k;
    ra := rand(-K..K);
    A := [seq(ra()+I*ra(),i=1..n)];
    k := sqrt(add(abs(A[i])^2,i=1..n));
    [seq(A[i]/k,i=1..n)];
end:

GeneralPauli := proc(n) local i;
    if nops(n)<>3 or sqrt(add(abs(n[i])^2,i=1..3))<>1 then
        RETURN(FAIL);
    fi;
    [[n[3],n[1]-I*n[2]],[n[1]+I*n[2],-n[3]]];
end:

# For RHM, check that eigenvectors are orthogonal and eigenvalues are real.
for _ from 1 to 3 do;
    L,E := op(EVC(RHM(3,10)));
    evals := [seq(evalc(Im(L[i])),i=1..3)]=[0,0,0];
    evecs := [IP(E[1],E[2]),IP(E[2],E[3]),IP(E[1],E[3])]=[0,0,0];
    print(evals and evecs);             # should be true
od:

# Check that the probability distribution adds up to 1.
for _ from 1 to 5 do;
    PD := ProbDist(RHM(3,10),RandomState(3,10))[1];
    print(evalf(add(PD[i],i=1..3)));    # should be close to 1
od: