<> <> DIRECTORY Rope, Basics, Atom, Ascii, Convert, IO, AlgebraClasses, Matrices; MatricesImpl: CEDAR PROGRAM IMPORTS Rope, Atom, Convert, IO EXPORTS Matrices = BEGIN OPEN AC: AlgebraClasses, Matrices; <> SyntaxError: PUBLIC ERROR [reason: ATOM] = CODE; BadElementStructure: PUBLIC ERROR [elementStructure: AC.Structure] = CODE; CantInvertError: PUBLIC ERROR [elementStructure: AC.Structure] = CODE; <> ClassPrintName: AC.PrintNameProc = { data: MatrixStructureData _ NARROW[structure.instanceData]; RETURN[Rope.Cat[ Convert.RopeFromCard[data.size], " x ", Convert.RopeFromCard[data.size], " Matrices over ", data.elementStructure.class.printName[data.elementStructure] ] ]; }; ClassCharacteristic: AC.StructureRankOp = { data: MatrixStructureData _ NARROW[structure.instanceData]; RETURN[ data.elementStructure.class.characteristic[data.elementStructure] ] }; ClassDimension: AC.StructureRankOp = { data: MatrixStructureData _ NARROW[structure.instanceData]; RETURN[data.size * data.size] }; ClassLegalFirstChar: AC.LegalFirstCharOp = { SELECT char FROM '[ => RETURN[TRUE]; ENDCASE; RETURN[FALSE]; }; ClassRead: AC.ReadOp = { data: MatrixStructureData _ NARROW[structure.instanceData]; RETURN[ReadMatrix[in, data.size, data.elementStructure] ]; }; ClassFromRope: AC.FromRopeOp = { stream: IO.STREAM _ IO.RIS[in]; RETURN[ ClassRead[stream, structure] ]; }; ClassToRope: AC.ToRopeOp = { matrix: Matrix _ NARROW[in]; data: MatrixStructureData _ NARROW[structure.instanceData]; RETURN[ MatrixToRope[matrix, data.elementStructure] ] -- use defaults }; ClassWrite: AC.WriteOp = { matrix: Matrix _ NARROW[in]; IO.PutRope[stream, ClassToRope[matrix, structure] ] }; ClassZero: AC.NullaryOp = { data: MatrixStructureData _ NARROW[structure.instanceData]; RETURN[ DiagMatrix[data.elementStructure.class.zero[data.elementStructure], data.size, data.elementStructure] ] }; ClassOne: AC.NullaryOp = { data: MatrixStructureData _ NARROW[structure.instanceData]; RETURN[ DiagMatrix[data.elementStructure.class.one[data.elementStructure], data.size, data.elementStructure] ] }; ClassAdd: AC.BinaryOp = { firstMatrix: Matrix _ NARROW[firstArg]; secondMatrix: Matrix _ NARROW[secondArg]; data: MatrixStructureData _ NARROW[structure.instanceData]; RETURN[ Add[firstMatrix, secondMatrix, data.elementStructure] ] }; ClassNegate: AC.UnaryOp = { matrix: Matrix _ NARROW[arg]; data: MatrixStructureData _ NARROW[structure.instanceData]; RETURN[ Negate[matrix, data.elementStructure] ] }; ClassSubtract: AC.BinaryOp = { firstMatrix: Matrix _ NARROW[firstArg]; secondMatrix: Matrix _ NARROW[secondArg]; data: MatrixStructureData _ NARROW[structure.instanceData]; RETURN[ Subtract[firstMatrix, secondMatrix, data.elementStructure] ] }; ClassMultiply: AC.BinaryOp = { firstMatrix: Matrix _ NARROW[firstArg]; secondMatrix: Matrix _ NARROW[secondArg]; data: MatrixStructureData _ NARROW[structure.instanceData]; RETURN[ Multiply[firstMatrix, secondMatrix, data.elementStructure] ] }; ClassEqual: AC.EqualityOp = { firstMatrix: Matrix _ NARROW[firstArg]; secondMatrix: Matrix _ NARROW[secondArg]; data: MatrixStructureData _ NARROW[structure.instanceData]; RETURN[ Equal[firstMatrix, secondMatrix, data.elementStructure] ] }; ClassTranspose: AC.UnaryOp = { matrix: Matrix _ NARROW[arg]; RETURN[ Transp[matrix] ] }; ClassDeterminant: AC.StructuredToGroundOp = { matrix: Matrix _ NARROW[structuredElt]; data: MatrixStructureData _ NARROW[structure.instanceData]; RETURN[ Det[matrix, data.elementStructure] ] }; matrixRingOps: MatrixOps _ NEW[MatrixOpsRec _ [ diagonalMatrix: DiagMatrix, transpose: ClassTranspose, determinant: ClassDeterminant ] ]; matrixProp: Atom.DottedPair _ NEW[Atom.DottedPairNode_ [$MatrixRing, matrixRingOps]]; -- having the $MatrixRing property actually means being a matrix structure but not an algebra <<>> matrixRingClass: PUBLIC AC.StructureClass _ NEW[AC.StructureClassRec _ [ flavor: ring, printName: ClassPrintName, characteristic: ClassCharacteristic, legalFirstChar: ClassLegalFirstChar, read: ClassRead, fromRope: ClassFromRope, toRope: ClassToRope, write: ClassWrite, add: ClassAdd, negate: ClassNegate, subtract: ClassSubtract, zero: ClassZero, multiply: ClassMultiply, commutative: FALSE, one: ClassOne, equal: ClassEqual, integralDomain: FALSE, gcdDomain: FALSE, euclideanDomain: FALSE, propList: LIST[matrixProp] ] ]; ClassInvert: AC.UnaryOp = { matrix: Matrix _ NARROW[arg]; data: MatrixStructureData _ NARROW[structure.instanceData]; RETURN[ Invert[matrix, data.elementStructure] ] }; ClassDivide: AC.BinaryOp = { firstMatrix: Matrix _ NARROW[firstArg]; secondMatrix: Matrix _ NARROW[secondArg]; data: MatrixStructureData _ NARROW[structure.instanceData]; RETURN[ Divide[firstMatrix, secondMatrix, data.elementStructure] ] }; ClassScalarMultiply: AC.BinaryOp = { scalar: REF _ firstArg; inMatrix: Matrix _ NARROW[secondArg]; data: MatrixStructureData _ NARROW[structure.instanceData]; RETURN[ ScalarMultiply[scalar, inMatrix, data.elementStructure] ] }; matrixAlgebraClass: PUBLIC AC.StructureClass _ NEW[AC.StructureClassRec _ [ flavor: algebra, printName: ClassPrintName, characteristic: ClassCharacteristic, dimension: ClassDimension, legalFirstChar: ClassLegalFirstChar, read: ClassRead, fromRope: ClassFromRope, toRope: ClassToRope, write: ClassWrite, add: ClassAdd, negate: ClassNegate, subtract: ClassSubtract, zero: ClassZero, multiply: ClassMultiply, commutative: FALSE, invert: ClassInvert, -- it is up to the user not to try to invert a singular matrix divide: ClassDivide, one: ClassOne, scalarMultiply: ClassScalarMultiply, equal: ClassEqual, integralDomain: FALSE, gcdDomain: FALSE, euclideanDomain: FALSE, propList: LIST[matrixProp] ] ]; <> MakeMatrixRing: PUBLIC PROC [elementStructure: AC.Structure, size: NAT] RETURNS [matrixRing: AC.Structure] ~ { matrixRingData: MatrixStructureData _ NEW[MatrixStructureDataRec _ [ elementStructure: elementStructure, size: size ] ]; IF elementStructure.class.flavor#ring AND elementStructure.class.flavor#algebra THEN ERROR BadElementStructure[elementStructure]; RETURN[ NEW[AC.StructureRec _ [ class: matrixRingClass, instanceData: matrixRingData ] ] ]; }; MakeMatrixAlgebra: PUBLIC PROC [elementField: AC.Structure, size: NAT] RETURNS [matrixAlgebra: AC.Structure] ~ { matrixAlgebraData: MatrixStructureData _ NEW[MatrixStructureDataRec _ [ elementStructure: elementField, size: size ] ]; IF elementField.class.flavor#field THEN ERROR BadElementStructure[elementField]; RETURN[ NEW[AC.StructureRec _ [ class: matrixAlgebraClass, instanceData: matrixAlgebraData ] ] ]; }; <> IsMatrixRing: PUBLIC PROC [structure: AC.Structure] RETURNS [BOOL] ~ { IF Atom.GetPropFromList[structure.class.propList, $MatrixRing] # NIL THEN RETURN[TRUE] ELSE RETURN[FALSE]; }; Transpose: PUBLIC PROC [structure: AC.Structure] RETURNS [AC.UnaryOp] ~ { IF IsMatrixRing[structure] THEN { matrixOps: MatrixOps _ NARROW[ Atom.GetPropFromList[structure.class.propList, $MatrixRing] ]; RETURN[matrixOps.transpose]; } ELSE ERROR; }; Determinant: PUBLIC PROC [structure: AC.Structure] RETURNS [AC.StructuredToGroundOp] ~ { IF IsMatrixRing[structure] THEN { matrixOps: MatrixOps _ NARROW[ Atom.GetPropFromList[structure.class.propList, $MatrixRing] ]; RETURN[matrixOps.determinant]; } ELSE ERROR; }; <> DiagMatrix: PUBLIC PROC [element: REF, size: NAT, elementStructure: AC.Structure] RETURNS [out: Matrix] ~ { rows: RowSeq _ NEW[RowSeqRec[size]]; FOR i:NAT IN [1..size] DO rows[i] _ MostlyZeroRow[element, i, size, elementStructure]; ENDLOOP; out _ NEW[MatrixRec _ [size: size, rows: rows] ]; }; MostlyZeroRow: PROC [element: REF, position, size: NAT, elementStructure: AC.Structure] RETURNS [row: Row] ~ { <> zero: REF _ elementStructure.class.zero[elementStructure]; row _ NEW[RowRec[size]]; FOR i:NAT IN [1..size] DO row[i] _ IF i = position THEN element ELSE zero; ENDLOOP; }; <> ReadMatrix: PUBLIC PROC [in: IO.STREAM, size: NAT, elementStructure: AC.Structure] RETURNS [matrix: Matrix] = { puncChar: CHAR; rows: RowSeq _ NEW[RowSeqRec[size]]; []_ in.SkipWhitespace[]; puncChar _ in.GetChar[ ]; IF puncChar # '[ THEN ERROR SyntaxError[$LeftBracketExpected]; FOR i:NAT IN [1..size] DO row: Row _ NEW[RowRec[size]]; []_ in.SkipWhitespace[]; puncChar _ in.GetChar[ ]; IF puncChar # '[ THEN ERROR SyntaxError[$LeftBracketExpected]; FOR j:NAT IN [1..size] DO row[j] _ elementStructure.class.read[in, elementStructure]; []_ in.SkipWhitespace[]; puncChar _ in.GetChar[]; []_ in.SkipWhitespace[]; IF j < size THEN IF puncChar # ', THEN ERROR SyntaxError[$CommaExpected]; ENDLOOP; IF puncChar # '] THEN ERROR SyntaxError[$RightBracketExpected]; rows[i] _ row; puncChar _ in.GetChar[]; []_ in.SkipWhitespace[]; IF i < size THEN IF puncChar # ', THEN SyntaxError[$CommaExpected]; ENDLOOP; IF puncChar # '] THEN ERROR SyntaxError[$RightBracketExpected]; matrix _ NEW[MatrixRec _ [size: size, rows: rows] ]; }; MatrixFromRope: PUBLIC PROC [in: Rope.ROPE, size: NAT, elementStructure: AC.Structure] RETURNS [out: Matrix] = { stream: IO.STREAM _ IO.RIS[in]; out _ ReadMatrix[stream, size, elementStructure]; }; MatrixToRope: PUBLIC PROC [in: Matrix, elementStructure: AC.Structure] RETURNS [out: Rope.ROPE] = { out _ "[ "; FOR i:NAT IN [1..in.size] DO out _ out.Concat["[ "]; FOR j:NAT IN [1..in.size] DO out _ out.Concat[elementStructure.class.toRope[in.rows[i][j], elementStructure] ]; IF j < in.size THEN out _ out.Concat[", "]; ENDLOOP; out _ out.Concat["] "]; IF i < in.size THEN out _ out.Concat[", "]; ENDLOOP; out _ out.Concat["] "]; }; WriteMatrix: PUBLIC PROC [in: Matrix, elementStructure: AC.Structure, out: IO.STREAM] = { out.PutF["\n %g \n", IO.rope[MatrixToRope[in, elementStructure]] ]; }; <> Create: PROC [size: NAT] RETURNS [a: Matrix] = { a _ NEW[MatrixRec _ [ size: size, rows: NEW[RowSeqRec[size]] ] ]; FOR i: NAT IN [1..size] DO a.rows[i] _ NEW[RowRec[size]]; ENDLOOP; }; Add: PUBLIC PROC [in1, in2: Matrix, elementStructure: AC.Structure] RETURNS [out: Matrix] ~ { out _ Create[in1.size]; FOR i:NAT IN [1..in1.size] DO FOR j:NAT IN [1..in1.size] DO out.rows[i][j] _ elementStructure.class.add[in1.rows[i][j], in2.rows[i][j], elementStructure]; ENDLOOP; ENDLOOP; }; Negate: PUBLIC PROC [in: Matrix, elementStructure: AC.Structure] RETURNS [out: Matrix] ~ { out _ Create[in.size]; FOR i:NAT IN [1..in.size] DO FOR j:NAT IN [1..in.size] DO out.rows[i][j] _ elementStructure.class.negate[in.rows[i][j], elementStructure]; ENDLOOP; ENDLOOP; }; Subtract: PUBLIC PROC [in1, in2: Matrix, elementStructure: AC.Structure] RETURNS [Matrix] ~ { RETURN[ Add[ in1, Negate[ in2, elementStructure], elementStructure ] ]; }; Multiply: PUBLIC PROC [in1, in2: Matrix, elementStructure: AC.Structure] RETURNS [out: Matrix] ~ { zero: REF _ elementStructure.class.zero[elementStructure]; out _ Create[in1.size]; FOR i:NAT IN [1..in1.size] DO FOR j:NAT IN [1..in1.size] DO out.rows[i][j] _ zero; FOR k:NAT IN [1..in1.size] DO out.rows[i][j] _ elementStructure.class.add[ out.rows[i][j], elementStructure.class.multiply[in1.rows[i][k], in2.rows[k][j], elementStructure], elementStructure] ENDLOOP; ENDLOOP; ENDLOOP; }; Invert: PUBLIC PROC [in: Matrix, elementField: AC.Structure] RETURNS [out: Matrix] ~ { det: REF _ Det[in, elementField]; tempVal: REF; Aij: Matrix _ Create[in.size-1]; IF elementField.class.flavor # field THEN ERROR CantInvertError[elementField]; out _ Create[in.size]; FOR i: NAT IN [1..in.size] DO FOR j: NAT IN [1..in.size] DO MakeAij[in, Aij, i, j]; tempVal _ elementField.class.divide[Det[Aij, elementField], det, elementField]; -- error will be generated here if Det[Aij] = 0 IF NOT Even[i + j] THEN tempVal _ elementField.class.negate[tempVal, elementField]; out.rows[j][i] _ tempVal; ENDLOOP; ENDLOOP; RETURN[out]; }; Divide: PUBLIC PROC [in1, in2: Matrix, elementField: AC.Structure] RETURNS [out: Matrix]~ { RETURN[ Multiply[ in1, Invert[ in2, elementField], elementField ] ]; }; ScalarMultiply: PUBLIC PROC [scalar: REF, inMatrix: Matrix, elementField: AC.Structure] RETURNS [outMatrix: Matrix]~ { outMatrix _ Create[inMatrix.size]; FOR i:NAT IN [1..inMatrix.size] DO FOR j:NAT IN [1..inMatrix.size] DO outMatrix.rows[i][j] _ elementField.class.multiply[scalar, inMatrix.rows[i][j], elementField]; ENDLOOP; ENDLOOP; }; <> Equal: PUBLIC PROC [in1, in2: Matrix, elementStructure: AC.Structure] RETURNS [BOOL _ TRUE] ~ { diff: Matrix _ Subtract[in1, in2, elementStructure]; FOR i: INTEGER IN [1..diff.size] DO FOR j: INTEGER IN [1..diff.size] DO IF NOT elementStructure.class.equal[diff.rows[i][j], elementStructure.class.zero[elementStructure], elementStructure] THEN RETURN[FALSE]; ENDLOOP; ENDLOOP; }; Transp: PUBLIC PROC [in: Matrix] RETURNS [out: Matrix] ~ { out _ Create[in.size]; FOR i: INTEGER IN [1..in.size] DO FOR j: INTEGER IN [1..in.size] DO out.rows[j][i] _ in.rows[i][j]; ENDLOOP; ENDLOOP; RETURN[out]; }; MakeAij: PROC[a, Aij: Matrix, i,j: INTEGER] = { --assume Aij is well formed n, m: NAT _ 1; --row index, column index for new matrix FOR row: NAT IN [1..a.size] DO IF row = i THEN LOOP; --row index for original matrix m _ 1; FOR col: NAT IN [1..a.size] DO --column index for original matrix IF col = j THEN LOOP; Aij.rows[n][m] _ a.rows[row][col]; m _ m+1; ENDLOOP; n _ n+1; ENDLOOP; }; Even: PROC[I: NAT] RETURNS[BOOL] = {RETURN[I MOD 2 = 0]}; Det: PUBLIC PROC [in: Matrix, elementStructure: AC.Structure] RETURNS [value: REF] ~ { <> IF in.size = 1 THEN value _ in.rows[1][1] ELSE IF in.size=2 THEN value _ elementStructure.class.subtract[ elementStructure.class.multiply[in.rows[1][1], in.rows[2][2], elementStructure], elementStructure.class.multiply[in.rows[1][2], in.rows[2][1], elementStructure], elementStructure] ELSE { i, j: NAT; zero: REF _ elementStructure.class.zero[elementStructure]; Aij: Matrix _ Create[in.size-1]; value _ zero; j _ 1; --always use column 1 for now FOR i IN [1..in.size] DO tempVal: REF; MakeAij[in, Aij, i, j]; tempVal _ elementStructure.class.multiply[in.rows[i][j], Det[Aij, elementStructure], elementStructure]; IF NOT Even[i + j] THEN tempVal _ elementStructure.class.negate[tempVal, elementStructure]; value _ elementStructure.class.add[value, tempVal, elementStructure]; ENDLOOP; }; RETURN[value]; }; END.