# # 1 Maple Code # This chapter introduces Quaternion*Package; # 1.1 In This Chapter Quaternion; # Declare some exported names interface('imaginaryunit' = ii); Q := module () local QMultiplyScalar, QMatrixMultiply, QScalarMultiply; export Qdef, Qreal, Qimag, Qconj, sortcollect, M, `*`, `^`, eval, Qrand, QPrand, Mrand, PMrand, Mconj, Qnorm, Qinv, AXXB, ZeroP, QMFNorm, PRotate, LagrangeP, AXBc, NewtonP, QAdjointMatrix, LFI, NoSim1, NoSim2, Meval, AXB, SemiSim, CSSim, LN1, LN2; option package; # 1.2 Define*a*Quaternion*with*four*coordinates; real, i, j and k; Qdef := proc (a, b, c, d) a+I*b+c*J+d*K end proc; # 1.3 Calculates*the*real*part*of*a*quaternion; Qreal := proc (a) subs(I = 0, J = 0, K = 0, a) end proc; # 1.4 Calculates*the*imaginary*part*of*a*quaternion; Qimag := proc (a) a-Qreal(a) end proc; # 1.5 returns A as a polynomial in standard form of variable a.; sortcollect := proc (A, a) return sort(collect(A, a), a) end proc; # 1.6 Calculates*the*conjugate*of*a*quaternion or a*quaternion*polynomial; Qconj := proc (f) sortcollect(Qreal(f)-Qimag(f), x) end proc; # 1.7 Calculates the inverse of a quaternion or a quaternion polynomial. Qinv := proc (a) sortcollect(Qreal(a)-Qimag(a), x)*(Qreal(a)^2+coeff(a, I)^2+coeff(a, J)^2+coeff(a, K)^2, x)/sortcollect end proc; # 1.8Generates*a*random*quaternion*with*all*four*integer*coordinates*between*m and n; Qrand := proc (m::integer, n::integer) Qdef((rand(m .. n))(), (rand(m .. n))(), (rand(m .. n))(), (rand(m .. n))()) end proc; # 1.9(Generates*a*random*quaternion*polynomial*of*degee*k*with*coefficients*being*random)*quaternions*Qrand(m, n); QPrand := proc (m, n, k) local j; j := k; return convert([seq(Qrand(m, n)*x^(j-i), i = 0 .. k)], `+`) end proc; # 1.10 `Calculates the 2 norm of a quaternion.`; Qnorm := proc (a) (Qreal(a)^2+coeff(a, I)^2+coeff(a, J)^2+coeff(a, K)^2)^(1/2) end proc; # 1.11 Generates a random p by q Matrix with quaternion entries Qrand(m,n); Mrand := proc (m, n, p, q) local M, i, j; M := Matrix(1 .. p, 1 .. q); for i to p do for j to q do M(i, j) := Qrand(m, n) end do end do end proc; # 1.12 Generates a random p by q Matrix with random quaternion polynomial entries; PMrand := proc (m, n, k1, k2, p, q) local M, i, j; M := Matrix(1 .. p, 1 .. q); for i to p do for j to q do M(i, j) := QPrand(m, n, (rand(k1 .. k2))()) end do end do end proc; # 1.13 The multiplication operator for quaternions and quaternion polynomials; M := proc (a, b) local a1, b1, c1, d1, a2, b2, c2, d2; a1 := Qreal(a); b1 := coeff(a, I); c1 := coeff(a, J); d1 := coeff(a, K); a2 := Qreal(b); b2 := coeff(b, I); c2 := coeff(b, J); d2 := coeff(b, K); return sortcollect*(a1*a2-b1*b2-c1*c2-d1*d2+I*(a1*b2+b1*a2+c1*d2-d1*c2)+(a1*c2-b1*d2+c1*a2+d1*b2)*J+(a1*d2+b1*c2-c1*b2+d1*a2)*K, x) end proc; # 1.14 Quaternions Matrix Multiply by Quaternions Matrix; QMatrixMultiply := proc (A::Matrix, B::Matrix) local nrowsA, ncolsA, ncolsB, AB, i, j; nrowsA := LinearAlgebra:-RowDimension(A); ncolsA := LinearAlgebra:-ColumnDimension(A); ncolsB := LinearAlgebra:-ColumnDimension(B); AB := Matrix(1 .. nrowsA, 1 .. ncolsB); for i to nrowsA do for j to ncolsB do AB(i, j) := simplify(add(M(A(i, k), B(k, j)), k = 1 .. ncolsA)) end do end do end proc; # 1.15 Quaternions*Matrix*Multiply*scalar; QMultiplyScalar := proc (A::Matrix, B) local nrowsA, ncolsA, AB, i, j; nrowsA := LinearAlgebra:-RowDimension(A); ncolsA := LinearAlgebra:-ColumnDimension(A); AB := Matrix(1 .. nrowsA, 1 .. ncolsA); for i to nrowsA do for j to ncolsA do AB(i, j) := simplify(M(A(i, j), B)) end do end do end proc; # 1.16 Quaternions*Matrix*Multiply*scalar; QScalarMultiply := proc (A, B::Matrix) local nrowsB, ncolsB, AB, i, j; nrowsB := LinearAlgebra:-RowDimension(B); ncolsB := LinearAlgebra:-ColumnDimension(B); AB := Matrix(1 .. nrowsB, 1 .. ncolsB); for i to nrowsB do for j to ncolsB do AB(i, j) := simplify(M(A, B(i, j))) end do end do end proc; # 1.17 Newly defined `*` for quaternions related elements.; `*` := proc (m, n) option remember; if type(m, Matrix) and type(n, Matrix) then QMatrixMultiply(m, n) elif type(m, Matrix) then QMultiplyScalar(m, n) elif type(n, Matrix) then QScalarMultiply(m, n) else M(m, n) end if end proc; # 1.18 Newly defined `^` to include quaternion raised to positive integer powers.; `^` := proc (m, n::integer) local i, p; option remember; p := 1; if type(m, freeof({I, J, K})) = false and 1 < n then for i to n do p := M(p, m) end do; return p else m^n end if end proc; # # 1.16 Calculates*the*conjugate*of*a*quaternion or quaternion*polynomial*matrix; Mconj := proc (A) local nrowsA, ncolsA, AC, i, j; nrowsA := LinearAlgebra:-RowDimension(A); ncolsA := LinearAlgebra:-ColumnDimension(A); AC := Matrix(1 .. nrowsA, 1 .. ncolsA); for i to nrowsA do for j to ncolsA do AC(i, j) := Qconj(A(i, j)) end do end do end proc; # 1.18 Calculates N(n) where N is a polynomial of x.; eval := proc (N, n) local i, f, m; option remember; f := 1; m := N; if has(N, x) = false then expand(N) elif is(expand(N), `+`) = true then map(eval, expand(N), n) else for i to degree(N, x) do f := M(f, n) end do; M(coeff(N, x, degree(N, x)), f) end if end proc; # 1.19 Calculates N(n) where N is a polynomial matrix of x.; Meval := proc (N, n) local i, j, m; option remember; m := Matrix(LinearAlgebra:-RowDimension(N), LinearAlgebra:-ColumnDimension(N)); for i to LinearAlgebra:-RowDimension(m) do for j to LinearAlgebra:-ColumnDimension(m) do m(i, j) := eval(N(i, j), n) end do end do end proc; # 1.20 Generates*a*block*matrix; PRotate := proc (m::integer) local M, i; M := Matrix(1 .. 2*m, 1 .. 2*m); for i to m do M(i, i+m) := -1 end do; for i from m+1 to 2*m do M(i, i-m) := 1 end do end proc; # 1.23 Ensures that a list of quaternion is ready for zero polynomial interpolation NoSim1 := proc (A::list) local i, j, k, B, n; option remember; B := A; n := nops(B); for i from 2 to n do for j to i-1 do if Qreal(B[j]) = Qreal(B[i]) and Qnorm(B[j]) = Qnorm(B[i]) then for k from j+1 to i-1 do if Qreal(B[k]) = Qreal(B[i]) and Qnorm(B[k]) = Qnorm(B[i]) then B := subsop(i = z, B); next end if end do end if end do end do; return remove(proc (t) options operator, arrow; t = z end proc, B) end proc; # 1.24 Calculates an zero polynomial for a list of quaternions; ZeroP := proc (B::list) local A, i, n, f; option remember; A := [op({op(NoSim1(B))})]; n := nops(A); f[1] := x-op(1, A); for i from 2 to n do f[i] := M(x-M(M(eval(f[i-1], op(i, A)), op(i, A)), Qinv(eval(f[i-1], op(i, A)))), f[i-1]) end do; return sortcollect(f[n], x) end proc; # 1.25 Ensures that a list of quaternion pairs is ready for Lagrange Interpolation.; NoSim2 := proc (A::listlist) local i, j, k, B, C, n1, n2; option remember; B := A; n1 := nops(B); for i from 2 to n1 do for j to i-1 do if A[j][1] = A[i][1] then B := subsop([i, 1] = z, B) end if end do end do; C := remove(proc (t) options operator, arrow; t[1] = z end proc, B); n2 := nops(C); for i from 3 to n2 do for j to i-2 do if Qreal(C[j][1]) = Qreal(C[i][1]) and Qnorm(C[j][1]) = Qnorm(C[i][1]) then for k from j+1 to i-1 do if Qreal(C[k][1]) = Qreal(C[i][1]) and Qnorm(C[k][1]) = Qnorm(C[i][1]) then C := subsop([i, 1] = z, C); next end if end do end if end do end do; return remove(proc (t) options operator, arrow; t[1] = z end proc, C) end proc; # 1.26 Calculates the Lagrange polynomial for a list of quaternion pairs.; LagrangeP := proc (B::listlist) local A, n, b, c, i, j, f, g; option remember; A := NoSim2(B); n := nops(A); for i to n do b[i] := op(1, op(i, A)); c[i] := op(2, op(i, A)) end do; for i to n do f[i] := ZeroP([seq(b[j], j in {seq(1 .. n)} minus {i})]) end do; for i to n do g[i] := M(Qinv(eval(f[i], b[i])), f[i]) end do; return sortcollect(add(M(c[i], g[i]), i = 1 .. n), x) end proc; # 1.27 Calculates the adjoint matrix of a quaternion matrix.; QAdjointMatrix := proc (A::Matrix) local nrowsA, ncolsA, A1, A1C, A2, A2C, AM, i, j; nrowsA := LinearAlgebra:-RowDimension(A); ncolsA := LinearAlgebra:-ColumnDimension(A); A1 := Matrix(1 .. nrowsA, 1 .. ncolsA); A1C := Matrix(1 .. nrowsA, 1 .. ncolsA); A2 := Matrix(1 .. nrowsA, 1 .. ncolsA); A2C := Matrix(1 .. nrowsA, 1 .. ncolsA); AM := Matrix(1 .. 2*nrowsA, 1 .. 2*ncolsA); for i to nrowsA do for j to ncolsA do A1(i, j) := A(i, j)-coeff(A(i, j), J)*J-coeff(A(i, j), K)*K; A1C(i, j) := Qconj(A1(i, j)); A2(i, j) := coeff(A(i, j), J)+I*coeff(A(i, j), K); A2C(i, j) := Qconj(A2(i, j)); AM(i, j) := A1(i, j) end do end do; for i to nrowsA do for j from ncolsA+1 to 2*ncolsA do AM(i, j) := A2(i, j-ncolsA) end do end do; for i from nrowsA+1 to 2*nrowsA do for j to ncolsA do AM(i, j) := -A2C(i-nrowsA, j) end do end do; for i from nrowsA+1 to 2*nrowsA do for j from ncolsA+1 to 2*ncolsA do AM(i, j) := A1C(i-nrowsA, j-ncolsA) end do end do; for i to 2*nrowsA do for j to 2*ncolsA do AM(i, j) := subs(I = ii, AM(i, j)) end do end do; return AM end proc; # 1.28 Calculates the frobinius norm of a quaternion matrix.; QMFNorm := proc (A::Matrix) local nrowsA, ncolsA, i, j, N; nrowsA := LinearAlgebra:-RowDimension(A); ncolsA := LinearAlgebra:-ColumnDimension(A); N := 0; for i to nrowsA do for j to ncolsA do N := N+coeff(A[i, j], I)^2+coeff(A[i, j], J)^2+coeff(A[i, j], K)^2+Qreal(A[i, j])^2 end do end do; return expand(N)^(1/2) end proc; # 1.29Calculates 1, The Newton polynomials and 2,The divided difference (coefficient of 1); NewtonP := proc (B::listlist) local g, i, j, d, b, c, A, n; option remember; A := NoSim2(B); n := nops(A); for i to n do b[i] := op(1, op(i, A)); c[i] := op(2, op(i, A)) end do; d[0] := c[1]; g[0] := 1; for i to n-1 do g[i] := ZeroP([seq(b[j], j = 1 .. i)]); d[i] := M(c[i+1]-c[1]-add(M(d[j], eval(g[j], b[i+1])), j = 1 .. i-1), Qinv(eval(g[i], b[i+1]))) end do; return [seq(g[i], i = 0 .. n-1)], [seq(d[i], i = 0 .. n-1)] end proc; # 1.30 Calculates the Moore-Penrose Pseudoinverse of a matrix.; LFI := proc (A::Matrix) local nrowsA, H, P, a, B, AA, i, j, k, LFI; nrowsA := LinearAlgebra:-RowDimension(A); H := LinearAlgebra:-Transpose(Mconj(A)); P := QMatrixMultiply(A, H); k := 0; a[0] := 1; B[0] := LinearAlgebra:-IdentityMatrix(nrowsA); for i to nrowsA do AA[i] := QMatrixMultiply(P, B[i-1]); a[i] := -LinearAlgebra:-Trace(AA[i])/i; B[i] := AA[i]+QScalarMultiply(a[i], B[0]) end do; for j to nrowsA do if a[j] <> 0 then k := j end if end do; if k = 0 then LFI := 0 else LFI := QScalarMultiply(-1/a[k], QMatrixMultiply(H, B[k-1])) end if end proc; # 1.31 of*the*x^2*Calculates*solution+Ax+C+xB = 0; AXXB := proc (A, B, C) local a0, a1, a2, a3, b0, b1, b2, b3, c, c0, c1, c2, c3, one, two, three, d, f, Z; a0 := Qreal(A); a1 := coeff(A, I); a2 := coeff(A, J); a3 := coeff(A, K); b0 := Qreal(B); b1 := coeff(B, I); b2 := coeff(B, J); b3 := coeff(B, K); c := ((a0+b0)*(1/2))^2-M(A+B, (a0+b0)*(1/2))+C; c0 := Qreal(c); c1 := coeff(c, I); c2 := coeff(c, J); c3 := coeff(c, K); one := -(4*(a1+b1))*x^3+(4*(a2*b3-a3*b2-c1))*x^2+(b1-a1)*(a1^2+a2^2+a3^2-b1^2-b2^2-b3^2)*x+(2*(c2*(b3-a3)+c3*(a2-b2)))*x+(b1-a1)*(c1*(a1-b1)+c2*(a2-b2)+c3*(a3-b3)); two := -(4*(a2+b2))*x^3+(4*(a3*b1-a1*b3-c2))*x^2+(b2-a2)*(a1^2+a2^2+a3^2-b1^2-b2^2-b3^2)*x+(2*(c1*(a3-b3)+c3*(b1-a1)))*x+(b2-a2)*(c1*(a1-b1)+c2*(a2-b2)+c3*(a3-b3)); three := -(4*(a3+b3))*x^3+(4*(a1*b2-a2*b1-c3))*x^2+(b3-a3)*(a1^2+a2^2+a3^2-b1^2-b2^2-b3^2)*x+(2*(c1*(b2-a2)+c2*(a1-b1)))*x+(b3-a3)*(c1*(a1-b1)+c2*(a2-b2)+c3*(a3-b3)); d := 8*x^3+(2*((a1-b1)^2+(a2-b2)^2+(a3-b3)^2))*x; f := 16*x^3+(8*(a1^2+a2^2+a3^2+b1^2+b2^2+b3^2+2*c0))*x^2+4*c0*((a1-b1)^2+(a2-b2)^2+(a3-b3)^2)*x+(a1^2+a2^2+a3^2-b1^2-b2^2-b3^2)^2*x+((8*(a2*b3-a3*b2))*c1-4*c1^2)*x+((8*(a3*b1-a1*b3))*c2-4*c2^2)*x+((8*(a1*b2-a2*b1))*c3-4*c3^2)*x-(a1*c1+a2*c2+a3*c3-b1*c1-b2*c2-b3*c3)^2; Z := evalf(op(select(proc (t) options operator, arrow; 0 < evalf(t) end proc, map(Im*ii+Re, simplify(evalc([solve(f, x)])))))); return map(proc (t) options operator, arrow; t+I*subs(x = t, one/d)+subs(x = t, two/d)*J+subs(x = t, three/d)*K-(1/2)*a0-(1/2)*b0 end proc, {-sqrt(Z), sqrt(Z)}) end proc; # 1.32 of*the*x^2*Calculates*solution+AxB+C = 0; AXBc := proc (A, B, c) local a1, a2, a3, b1, b2, b3, one, two, three, a, b, d, f, sol; a1 := coeff(A, I); a2 := coeff(A, J); a3 := coeff(A, K); b1 := coeff(B, I); b2 := coeff(B, J); b3 := coeff(B, K); one := x*(a3*b2-a2*b3); two := x*(a1*b3-a3*b1); three := x*(a2*b1-a1*b2); a := (a1^2+a2^2+a3^2)*(b1^2+b2^2+b3^2); b := a1*b1+a2*b2+a3*b3; d := 2*x+b; f := x^4+(c-(3/4)*a)*x^2+(1/4)*b*(4*c-a)*x+(1/4)*c*b^2; sol := fsolve(f, x, real); return map(proc (t) options operator, arrow; t+I*subs(x = t, one/d)+subs(x = t, two/d)*J+subs(x = t, three/d)*K end proc, {sol}) end proc; # 1.33 of*the*x^2*Calculates*solution+AxB = 0; AXB := proc (A, B) local a0, a1, a2, a3, b0, b1, b2, b3, m11, m12, m13, m21, m22, m23, m31, m32, m33, d, n1, n2, n3, one, two, three, ab, h, g, solg, solh; a0 := Qreal(A); a1 := coeff(A, I); a2 := coeff(A, J); a3 := coeff(A, K); b0 := Qreal(B); b1 := coeff(B, I); b2 := coeff(B, J); b3 := coeff(B, K); m11 := a0*b0-a1*b1+a2*b2+a3*b3; m12 := a0*b3-a1*b2-a2*b1-a3*b0; m13 := -a0*b2-a1*b3+a2*b0-a3*b1; m21 := -a0*b3-a1*b2-a2*b1+a3*b0; m22 := a0*b0+a1*b1-a2*b2+a3*b3; m23 := a0*b1-a1*b0-a2*b3-a3*b2; m31 := a0*b2-a1*b3-a2*b0-a3*b1; m32 := -a0*b1+a1*b0-a2*b3-a3*b2; m33 := a0*b0+a1*b1+a2*b2-a3*b3; d := 8*x^3+(4*(m11+m22+m33))*x^2+(2*(m11*m22+m11*m33-m12*m21-m13*m31-m23*m32+m22*m33))*x+m11*m22*m33-m11*m23*m32-m12*m21*m33+m13*m21*m32+m12*m23*m31-m13*m22*m31; n1 := -a0*b1-a1*b0-a2*b3+a3*b2; n2 := -a0*b2+a1*b3-a2*b0-a3*b1; n3 := -a0*b3-a1*b2+a2*b1-a3*b0; one := 4*n1*x^3+(2*(n1*(m22+m33)-n2*m12-n3*m13))*x^2+n1*(m22*m33-m23*m32)*x+n2*(m13*m32-m12*m33)*x+n3*(m12*m23-m13*m22)*x; two := 4*n2*x^3+(2*(n2*(m11+m33)-n1*m21-n3*m23))*x^2+n1*(m23*m31-m21*m33)*x+n2*(m11*m33-m13*m31)*x+n3*(m13*m21-m11*m23)*x; three := 4*n3*x^3+(2*(n3*(m11+m22)-n1*m31-n2*m32))*x^2+n1*(m21*m32-m22*m31)*x+n2*(m12*m31-m11*m32)*x+n3*(m11*m22-m12*m21)*x; ab := (a0^2+a1^2+a2^2+a3^2)*(b0^2+b1^2+b2^2+b3^2); h := x^3+2*a0*b0*x^2+((4/3)*a0^2*b0^2-(1/12*(a0^2-3*a1^2-3*a2^2-3*a3^2))*(b0^2-3*b1^2-3*b2^2-3*b3^2))*x+(1/4)*ab*(a0*b0-a1*b1-a2*b2-a3*b3); g := x^4+2*a0*b0*x^3+(1/2*(4*a0^2*b0^2-(a1^2+a2^2+a3^2-a0^2)*(b1^2+b2^2+b3^2-b0^2)))*x^2+(1/2)*a0*b0*ab*x+(1/16)*ab^2; solg := fsolve(g, x, real); solh := fsolve(h, x, real); if a0 = 0 and b0 = 0 then return map(proc (t) options operator, arrow; t+I*subs(x = t, one/d)+subs(x = t, two/d)*J+subs(x = t, three/d)*K end proc, {solh}) else return map(proc (t) options operator, arrow; t+I*subs(x = t, one/d)+subs(x = t, two/d)*J+subs(x = t, three/d)*K end proc, {solg, solh}) end if end proc; # 1.34 Calculates*x, y*such*that*xBy = A and yAx = B; SemiSim := proc (B, A) local a0, a1, a2, a3, b0, b1, b2, b3, Z; a0 := Qreal(A); a1 := coeff(A, I); a2 := coeff(A, J); a3 := coeff(A, K); b0 := Qreal(B); b1 := coeff(B, I); b2 := coeff(B, J); b3 := coeff(B, K); if Qnorm(A) <> Qnorm(B) then return "No Solution" elif abs(a0) <> abs(b0) then return "No Solution" else Z := LinearAlgebra:-LinearSolve(Matrix([[a0^2-b0^2, -a0*a3-b0*b3, a0*a2+b0*b2, a0*a1-b0*b1, 0], [a0*a3+b0*b3, a0^2-b0^2, -a0*a1-b0*b1, a0*a2-b0*b2, 0], [-a0*a2-b0*b2, a0*a1+b0*b1, a0^2-b0^2, a0*a3-b0*b3, 0], [-a0*a1+b0*b1, -a0*a2+b0*b2, -a0*a3+b0*b3, a0^2-b0^2, 0]]), 'free' = 't') end if; return sort(Z[4]+I*Z[1]+Z[2]*J+Z[3]*K, [I, J, K]) and simplify(M(B, M(Qinv(Z[4]+I*Z[1]+Z[2]*J+Z[3]*K), Qinv(A)))) end proc; # 1.35 CSSim := proc (B, A) local a0, a1, a2, a3, b0, b1, b2, b3, Z; a0 := Qreal(A); a1 := coeff(A, I); a2 := coeff(A, J); a3 := coeff(A, K); b0 := Qreal(B); b1 := coeff(B, I); b2 := coeff(B, J); b3 := coeff(B, K); if Qnorm(A) <> Qnorm(B) then return "No Solution" else Z := LinearAlgebra:-LinearSolve(Matrix([[0, a0*b3+a3*b0, -a0*b2-a2*b0, a2*b3-a3*b2, 0], [-a0*b3-a3*b0, 0, a0*b1+a1*b0, -a1*b3+a3*b1, 0], [a0*b2+a2*b0, -a0*b1-a1*b0, 0, a1*b2-a2*b1, 0], [-a2*b3+a3*b2, a1*b3-a3*b1, -a1*b2+a2*b1, 0, 0]]), 'free' = 't') end if; return simplify(Z[4]+I*Z[1]+Z[2]*J+Z[3]*K) and simplify(M(M(Qinv(B), Qinv(Qconj(Z[4]+I*Z[1]+Z[2]*J+Z[3]*K))), A)) end proc; # 1.36 LN1 := proc (A, B, C) local a0, a1, a2, a3, b0, b1, b2, b3, c0, c1, c2, c3, Z; a0 := Qreal(A); a1 := coeff(A, I); a2 := coeff(A, J); a3 := coeff(A, K); b0 := Qreal(B); b1 := coeff(B, I); b2 := coeff(B, J); b3 := coeff(B, K); c0 := Qreal(C); c1 := coeff(C, I); c2 := coeff(C, J); c3 := coeff(C, K); if Qnorm(A) = Qnorm(B) and Qreal(A) = Qreal(B) then Z := LinearAlgebra:-MatrixMatrixMultiply(LinearAlgebra:-MatrixInverse(Matrix([[a0-b0, -a3-b3, a2+b2, a1-b1], [a3+b3, a0-b0, -a1-b1, a2-b2], [-a2-b2, a1+b1, a0-b0, a3-b3], [-a1+b1, -a2+b2, -a3+b3, a0-b0]]) , 'method' = 'pseudo'), Matrix(4, 1, {(1, 1) = c1, (2, 1) = c2, (3, 1) = c3, (4, 1) = c0})) else Z := LinearAlgebra:-LinearSolve(Matrix([[a0-b0, -a3-b3, a2+b2, a1-b1, c1], [a3+b3, a0-b0, -a1-b1, a2-b2, c2], [-a2-b2, a1+b1, a0-b0, a3-b3, c3], [-a1+b1, -a2+b2, -a3+b3, a0-b0, c0]]), 'free' = 't') end if; return simplify(Z[4]+I*Z[1]+Z[2]*J+Z[3]*K) end proc; # 1.37 LN2 := proc (A, B, C) local a0, a1, a2, a3, b0, b1, b2, b3, c0, c1, c2, c3, Z; a0 := Qreal(A); a1 := coeff(A, I); a2 := coeff(A, J); a3 := coeff(A, K); b0 := Qreal(B); b1 := coeff(B, I); b2 := coeff(B, J); b3 := coeff(B, K); c0 := Qreal(C); c1 := coeff(C, I); c2 := coeff(C, J); c3 := coeff(C, K); if Qnorm(A) = Qnorm(B) and Qreal(A) = Qreal(B) then Z := LinearAlgebra:-MatrixMatrixMultiply(LinearAlgebra:-MatrixInverse(Matrix([[a0+b0, -a3+b3, a2-b2, a1-b1], [a3-b3, a0+b0, -a1+b1, a2-b2], [-a2+b2, a1-b1, a0+b0, a3-b3], [-a1-b1, -a2-b2, -a3-b3, a0-b0]]) , 'method' = 'pseudo') , Matrix(4, 1, {(1, 1) = c1, (2, 1) = c2, (3, 1) = c3, (4, 1) = c0})) else Z := LinearAlgebra:-LinearSolve(Matrix([[a0+b0, -a3+b3, a2-b2, a1-b1, c1], [a3-b3, a0+b0, -a1+b1, a2-b2, c2], [-a2+b2, a1-b1, a0+b0, a3-b3, c3], [-a1-b1, -a2-b2, -a3-b3, a0-b0, c0]]) , 'free' = 't') end if; return simplify(Z[4]+I*Z[1]+Z[2]*J+Z[3]*K) end proc; end module: