MatricesImpl.mesa
Last Edited by: Arnon, June 10, 1985 4:19:22 pm PDT
DIRECTORY
Rope,
Basics,
Atom,
Ascii,
Convert,
IO,
AlgebraClasses,
Structures,
MathConstructors,
MathExpr,
Bools,
Ints,
Matrices;
MatricesImpl: CEDAR PROGRAM
IMPORTS Rope, IO, Convert, AlgebraClasses, Structures, Ints, MathConstructors
EXPORTS Matrices =
BEGIN OPEN AC: AlgebraClasses, Matrices;
Types
ROPE: TYPE = Rope.ROPE;
STREAM: TYPE = IO.STREAM;
Object: TYPE = AC.Object;
Method: TYPE = AC.Method;
Errors
SyntaxError: PUBLIC ERROR [reason: ATOM] = CODE;
BadElementStructure: PUBLIC ERROR [elementStructure: Object] = CODE;
CantInvertError: PUBLIC ERROR [elementStructure: Object] = CODE;
TypeError: PUBLIC ERROR [message: ATOM ← $Unspecified] = CODE;
Matrix Structure Ops
MakeMatrixStructure: PUBLIC AC.MatrixStructureConstructor ~ {
matrixStructureData: MatrixStructureData ← NEW[MatrixStructureDataRec ← [
elementStructure: elementStructure,
nRows: nRows,
nCols: nCols
] ];
method: Method ← AC.LookupMethodInStructure[$category, elementStructure];
category: REF AC.Category ← NARROW[method.value];
IF nRows#nCols OR category^=group THEN
matrixStructure ← AC.MakeStructure[
name: NIL,
class: matrixClass,
instanceData: matrixStructureData
]
ELSE SELECT category^ FROM -- assumes square matrices
ring, algebra => matrixStructure ← AC.MakeStructure[
name: NIL,
class: squareMatricesOverRingClass,
instanceData: matrixStructureData
];
field, divisionAlgebra => matrixStructure ← AC.MakeStructure[
name: NIL,
class: squareMatricesOverFieldClass,
instanceData: matrixStructureData
];
ENDCASE => ERROR BadElementStructure[elementStructure];
matrixStructure.name ← ShortPrintName[matrixStructure];
IF AC.LookupStructure[matrixStructure.name] = NIL THEN AC.InstallStructure[matrixStructure];
};
PrintName: PUBLIC AC.PrintNameProc = {
data: MatrixStructureData ← NARROW[structure.data];
RETURN[Rope.Cat[
Convert.RopeFromCard[data.size],
" x ",
Convert.RopeFromCard[data.size],
" Matrices over ",
data.elementStructure.class.printName[data.elementStructure]
] ];
};
ShortPrintName: PUBLIC AC.ToRopeOp = {
data: MatrixStructureData ← NARROW[in.data];
shortPrintNameMethod: Method ← AC.LookupMethodInStructure[$shortPrintName, data.elementStructure];
RETURN[Rope.Cat[
Rope.Cat[
"M",
Convert.RopeFromCard[data.nRows],
",",
Convert.RopeFromCard[data.nCols]
],
"(",
NARROW[AC.ApplyNoLkpNoRecastRef[shortPrintNameMethod,LIST[data.elementStructure] ] ],
")"
] ];
};
Characteristic: PUBLIC StructureRankOp = {
data: MatrixStructureData ← NARROW[structure.data];
RETURN[ data.elementStructure.class.characteristic[data.elementStructure] ]
};
Dimension: PUBLIC StructureRankOp = {
data: MatrixStructureData ← NARROW[structure.data];
RETURN[data.size * data.size]
};
IsMatrixStructure: PUBLIC AC.UnaryPredicate = {
IF ~arg.flavor = Structure THEN RETURN[FALSE];
RETURN[ AC.LookupMethodInStructure[$matrixStructure, arg]#NIL ]
};
Conversion and IO
Recast: PUBLIC AC.BinaryOp = {
Args are a StructureElement and a Structure
thisMatrixStructure: Object ← secondArg;
thisMatrixStructureData: MatrixStructureData ← NARROW[thisMatrixStructure.data];
thisStructureElementStructure: Object ← thisMatrixStructureData.elementStructure;
canRecastMethod: Method ← AC.LookupMethodInStructure[$canRecast, thisStructureElementStructure];
recastMethod: Method ← AC.LookupMethodInStructure[$recast, thisStructureElementStructure];
refBOOL: REF BOOL;
Already in this structure
IF AC.StructureEqual[firstArg.class, secondArg] THEN RETURN[firstArg]; -- nothing to do
Try to recast into elementStructure
refBOOL ← NARROW[AC.ApplyNoLkpNoRecastRef[canRecastMethod,LIST[firstArg.class, thisStructureElementStructure] ] ];
IF refBOOL^ THEN {
recastElement: Object ← AC.ApplyNoLkpNoRecastObject[recastMethod, LIST[firstArg, thisStructureElementStructure] ];
RETURN[DiagonalMatrix[recastElement, thisMatrixStructure] ];
};
If a matrix of the right size, try to recast its elements into elementStructure
IF AC.LookupMethodInStructure[$matrixStructure, firstArg.class]#NIL THEN {
inputMatrixStructure: Object ← firstArg.class;
inputMatrixStructureData: MatrixStructureData ← NARROW[inputMatrixStructure.data];
inputStructureElementStructure: Object ← inputMatrixStructureData.elementStructure;
argData: MatrixData ← NARROW[firstArg.data];
resultData: MatrixData;
IF inputMatrixStructureData.nRows # thisMatrixStructureData.nRows THEN RETURN[NIL];
IF inputMatrixStructureData.nCols # thisMatrixStructureData.nCols THEN RETURN[NIL];
refBOOL ← NARROW[AC.ApplyNoLkpNoRecastRef[canRecastMethod,LIST[inputStructureElementStructure, thisStructureElementStructure] ] ];
IF NOT refBOOL^ THEN RETURN[NIL];
result ← Create[thisMatrixStructure];
resultData ← NARROW[result.data];
FOR i:NAT IN [1..thisMatrixStructureData.nRows] DO
FOR j:NAT IN [1..thisMatrixStructureData.nCols] DO
resultData[i][j] ← AC.ApplyNoLkpNoRecastObject[recastMethod, LIST[argData[i][j], thisStructureElementStructure] ];
ENDLOOP;
ENDLOOP;
RETURN[result];
};
Can't do it
RETURN[NIL];
};
CanRecast: PUBLIC AC.BinaryPredicate = {
Args are either [Structure, Structure] or [StructureElement, Structure]
thisMatrixStructure: Object ← secondArg;
thisMatrixStructureData: MatrixStructureData ← NARROW[thisMatrixStructure.data];
thisStructureElementStructure: Object ← thisMatrixStructureData.elementStructure;
canRecastMethod: Method ← AC.LookupMethodInStructure[$canRecast, thisStructureElementStructure];
refBOOL: REF BOOL;
firstArgStructure: Object ← IF firstArg.flavor = StructureElement THEN firstArg.class ELSE IF firstArg.flavor = Structure THEN firstArg ELSE ERROR;
firstArgStructure = thisMatrixStructure
IF AC.StructureEqual[firstArgStructure, thisMatrixStructure] THEN RETURN[TRUE];
CanRecast firstArgStructure into elementStructure?
refBOOL ← NARROW[AC.ApplyNoLkpNoRecastRef[canRecastMethod,LIST[firstArgStructure, thisStructureElementStructure] ] ];
IF refBOOL^ THEN RETURN[TRUE];
If firstArgStructure a matrix Structure of the right size, CanRecast its elementStructure into this elementStructure?
IF AC.LookupMethodInStructure[$matrixStructure, firstArgStructure]#NIL THEN {
inputMatrixStructure: Object ← firstArgStructure;
inputMatrixStructureData: MatrixStructureData ← NARROW[inputMatrixStructure.data];
inputStructureElementStructure: Object ← inputMatrixStructureData.elementStructure;
IF inputMatrixStructureData.nRows # thisMatrixStructureData.nRows THEN RETURN[FALSE];
IF inputMatrixStructureData.nCols # thisMatrixStructureData.nCols THEN RETURN[FALSE];
refBOOL ← NARROW[AC.ApplyNoLkpNoRecastRef[canRecastMethod,LIST[inputStructureElementStructure, thisStructureElementStructure] ] ];
IF NOT refBOOL^ THEN RETURN[FALSE] ELSE RETURN[TRUE];
};
Out of luck
RETURN[FALSE];
};
ToExpr: PUBLIC AC.ToExprOp = {
matrixData: MatrixData ← NARROW[in.data];
matrixStructureData: MatrixStructureData ← NARROW[in.class.data];
outMatrix: LIST OF LIST OF MathExpr.EXPR ← NIL;
method: Method ← AC.LookupMethodInStructure[$toExpr, matrixStructureData.elementStructure];
FOR i:NAT DECREASING IN [1..matrixStructureData.nRows] DO
outRow: LIST OF MathExpr.EXPR ← NIL;
FOR j:NAT DECREASING IN [1..matrixStructureData.nCols] DO
outRow ← CONS[matrixStructureData.elementStructure.class.toExpr[matrixData[i][j] ], outRow];
outRow ← CONS[NARROW[AC.ApplyNoLkpNoRecastRef[method, LIST[matrixData[i][j] ] ] ], outRow];
ENDLOOP;
outMatrix ← CONS[outRow, outMatrix];
ENDLOOP;
out ← MathConstructors.MakeMatrix[matrixStructureData.nRows, matrixStructureData.nCols, outMatrix];
};
LegalFirstChar: PUBLIC AC.LegalFirstCharOp = {
SELECT char FROM
'[ => RETURN[TRUE];
ENDCASE;
RETURN[FALSE];
};
Read: PUBLIC AC.ReadOp = {
data: MatrixStructureData ← NARROW[structure.data];
elementStructure: AC.Object ← data.elementStructure;
puncChar: CHAR;
rows: RowSeq ← NEW[RowSeqRec[data.nRows]];
readMethod: AC.Method ← AC.LookupMethodInStructure[$read, elementStructure];
[]← in.SkipWhitespace[];
puncChar ← in.GetChar[];
IF puncChar # '[ THEN ERROR SyntaxError[$LeftBracketExpected];
FOR i:NAT IN [1..data.nRows] DO
row: Row ← NEW[RowRec[data.nCols]];
[]← in.SkipWhitespace[];
puncChar ← in.GetChar[];
IF puncChar # '[ THEN ERROR SyntaxError[$LeftBracketExpected];
FOR j:NAT IN [1..data.nCols] DO
row[j] ← AC.ApplyReadMethod[readMethod, in, elementStructure];
[]← in.SkipWhitespace[];
puncChar ← in.GetChar[];
[]← in.SkipWhitespace[];
IF j < data.nCols THEN IF puncChar # ', THEN ERROR SyntaxError[$CommaExpected];
ENDLOOP;
IF puncChar # '] THEN ERROR SyntaxError[$RightBracketExpected];
rows[i] ← row;
puncChar ← in.GetChar[];
[]← in.SkipWhitespace[];
IF i < data.nRows THEN IF puncChar # ', THEN SyntaxError[$CommaExpected];
ENDLOOP;
IF puncChar # '] THEN ERROR SyntaxError[$RightBracketExpected];
out ← NEW[AC.ObjectRec ← [flavor: StructureElement, class: structure, data: rows] ];
};
FromRope: PUBLIC AC.FromRopeOp = {
stream: IO.STREAM ← IO.RIS[in];
out ← Read[stream, structure];
};
ToRope: PUBLIC AC.ToRopeOp = {
data: MatrixData ← NARROW[in.data];
matrixStructureData: MatrixStructureData ← NARROW[in.class.data];
elementStructure: AC.Object ← matrixStructureData.elementStructure;
toRopeMethod: AC.Method ← AC.LookupMethodInStructure[$toRope, elementStructure];
out ← "[ ";
FOR i:NAT IN [1..matrixStructureData.nRows] DO
out ← out.Concat["[ "];
FOR j:NAT IN [1..matrixStructureData.nCols] DO
out ← out.Concat[NARROW[AC.ApplyNoLkpNoRecastRef[toRopeMethod, LIST[data[i][j] ] ] ] ];
IF j < matrixStructureData.nCols THEN out ← out.Concat[", "];
ENDLOOP;
out ← out.Concat["] "];
IF i < matrixStructureData.nRows THEN out ← out.Concat[", "];
ENDLOOP;
out ← out.Concat["] "];
};
Write: PUBLIC AC.WriteOp = {
stream.PutF["\n %g \n", IO.rope[ToRope[in]] ];
};
Constructors
DiagonalMatrix: PUBLIC AC.UnaryImbedOp ~ {
doesn't require nRows = nCols, just lays arg down the diagonal
matrixStructureData: MatrixStructureData ← NARROW[structure.data];
rows: RowSeq ← NEW[RowSeqRec[matrixStructureData.nRows]];
IF NOT AC.StructureEqual[in.class, matrixStructureData.elementStructure] THEN TypeError[]; -- check that we really are constructing an element of matrixStructure
FOR i:NAT IN [1..matrixStructureData.nRows] DO
rows[i] ← MostlyZeroRow[in, i, matrixStructureData.nCols];
ENDLOOP;
out ← NEW[AC.ObjectRec ← [class: structure, flavor: StructureElement, data: rows] ];
};
MostlyZeroRow: PROC [element: Object, position, size: NAT] RETURNS [row: Row] ~ {
Zeros everywhere except element in position.
zero: Object ← AC.ApplyLkpRecastObject[$zero, element.class, LIST[element.class] ];
row ← NEW[RowRec[size]];
FOR i:NAT IN [1..size] DO
row[i] ← IF i = position THEN element ELSE zero;
ENDLOOP;
};
MakeMatrix: PUBLIC AC.MatrixImbedOp ~ {
data: MatrixStructureData ← NARROW[structure.data];
nRows: NAT ← data.nRows;
nCols: NAT ← data.nCols;
rows: RowSeq ← NEW[RowSeqRec[nRows]];
FOR i:NAT IN [1..nRows] DO
row: Row ← NEW[RowRec[nCols]];
FOR j:NAT IN [1..nCols] DO
row[j] ← elements.first;
elements ← elements.rest
ENDLOOP;
rows[i] ← row;
ENDLOOP;
out ← NEW[AC.ObjectRec ← [flavor: StructureElement, class: structure, data: rows] ];
};
Arithmetic
Create: PROC [matrixStructure: Object] RETURNS [a: Matrix] = {
data: MatrixData;
matrixStructureData: MatrixStructureData ← NARROW[matrixStructure.data];
a ← NEW[AC.ObjectRec ← [
flavor: StructureElement,
class: matrixStructure,
data: NEW[RowSeqRec[matrixStructureData.nRows]]
] ];
data ← NARROW[a.data];
FOR i: NAT IN [1..matrixStructureData.nRows] DO data[i] ← NEW[RowRec[matrixStructureData.nCols]]; ENDLOOP;
};
Zero: PUBLIC AC.NullaryOp = {
matrixStructureData: MatrixStructureData ← NARROW[structure.data];
elementStructure: Object ← matrixStructureData.elementStructure;
zeroMethod: Method ← AC.LookupMethodInStructure[$zero, elementStructure];
resultData: MatrixData;
result ← Create[structure];
resultData ← NARROW[result.data];
FOR i:NAT IN [1..matrixStructureData.nRows] DO
FOR j:NAT IN [1..matrixStructureData.nCols] DO
resultData[i][j] ← AC.ApplyNoLkpNoRecastObject[zeroMethod, LIST[elementStructure] ];
ENDLOOP;
ENDLOOP;
};
One: PUBLIC AC.NullaryOp = {
data: MatrixStructureData ← NARROW[structure.data];
oneMethod: Method ← AC.LookupMethodInStructure[$one, data.elementStructure];
RETURN[ DiagonalMatrix[
AC.ApplyNoLkpNoRecastObject[oneMethod, LIST[data.elementStructure] ],
structure] ];
};
Add: PUBLIC AC.BinaryOp ~ {
firstData: MatrixData ← NARROW[firstArg.data];
secondData: MatrixData ← NARROW[secondArg.data];
resultData: MatrixData;
matrixStructureData: MatrixStructureData ← NARROW[firstArg.class.data];
result ← Create[firstArg.class];
resultData ← NARROW[result.data];
FOR i:NAT IN [1..matrixStructureData.nRows] DO
FOR j:NAT IN [1..matrixStructureData.nCols] DO
resultData[i][j] ← matrixStructureData.elementStructure.class.add[firstData[i][j], secondData[i][j] ];
resultData[i][j] ← AC.ApplyLkpRecastObject[$sum, matrixStructureData.elementStructure, LIST[firstData[i][j], secondData[i][j] ] ];
ENDLOOP;
ENDLOOP;
};
Negate: PUBLIC AC.UnaryOp ~ {
argData: MatrixData ← NARROW[arg.data];
resultData: MatrixData;
matrixStructureData: MatrixStructureData ← NARROW[arg.class.data];
result ← Create[arg.class];
resultData ← NARROW[result.data];
FOR i:NAT IN [1..matrixStructureData.nRows] DO
FOR j:NAT IN [1..matrixStructureData.nCols] DO
resultData[i][j] ← matrixStructureData.elementStructure.class.negate[argData[i][j]];
resultData[i][j] ← AC.ApplyLkpRecastObject[$negation, matrixStructureData.elementStructure, LIST[argData[i][j]] ];
ENDLOOP;
ENDLOOP;
};
Subtract: PUBLIC AC.BinaryOp ~ {
RETURN[ Add[ firstArg, Negate[ secondArg] ] ];
};
Multiply: PUBLIC AC.BinaryOp ~ {
Assumes matrices from same structure, i.e. square
firstData: MatrixData ← NARROW[firstArg.data];
secondData: MatrixData ← NARROW[secondArg.data];
resultData: MatrixData;
matrixStructureData: MatrixStructureData ← NARROW[firstArg.class.data];
zero: Object ← matrixStructureData.elementStructure.class.zero[matrixStructureData.elementStructure];
zero: Object ← AC.ApplyLkpRecastObject[$zero, matrixStructureData.elementStructure, LIST[matrixStructureData.elementStructure] ];
result ← Create[firstArg.class];
resultData ← NARROW[result.data];
FOR i:NAT IN [1..matrixStructureData.nRows] DO
FOR j:NAT IN [1..matrixStructureData.nCols] DO
resultData[i][j] ← zero;
FOR k:NAT IN [1..matrixStructureData.nCols] DO
resultData[i][j] ← matrixStructureData.elementStructure.class.add[
resultData[i][j],
matrixStructureData.elementStructure.class.multiply[firstData[i][k], secondData[k][j] ] ]
prod: Object ← AC.ApplyLkpRecastObject[$product, matrixStructureData.elementStructure, LIST[firstData[i][k], secondData[k][j]] ];
resultData[i][j] ← AC.ApplyLkpRecastObject[$sum, matrixStructureData.elementStructure, LIST[resultData[i][j], prod ] ];
ENDLOOP;
ENDLOOP;
ENDLOOP;
};
Power: PUBLIC AC.BinaryOp ~ { -- this simple algorithm is Structure independent
power: INT ← Ints.ToINT[secondArg];
structure: Object ← firstArg.class;
one: Object ← AC.ApplyLkpNoRecastObject[$one, structure, LIST[structure] ];
productMethod: Method ← AC.LookupMethodInStructure[$product, structure];
IF power < 0 THEN {
invertMethod: Method ← AC.LookupMethodInStructure[$invert, structure];
temp: Object;
IF invertMethod = NIL THEN ERROR;
temp ← Power[firstArg, Ints.FromINT[ABS[power] ] ];
RETURN[AC.ApplyNoLkpNoRecastObject[invertMethod, LIST[temp] ] ];
};
IF power = 0 THEN RETURN[one];
result ← firstArg;
FOR i:INT IN [2..power] DO
result ← AC.ApplyNoLkpNoRecastObject[productMethod, LIST[firstArg, result] ];
ENDLOOP;
};
Invert: PUBLIC AC.UnaryOp ~ {
resultData: MatrixData;
matrixStructureData: MatrixStructureData ← NARROW[arg.class.data];
elementStructure: Object ← matrixStructureData.elementStructure;
det: Object ← Determinant[arg];
tempVal: Object;
scratchMatrixRing: Object ← MakeMatrixStructure[elementStructure, matrixStructureData.nRows - 1, matrixStructureData.nCols - 1];
Aij: Matrix ← Create[scratchMatrixRing];
IF NOT AC.IsCategory[elementStructure, field] AND NOT AC.IsCategory[elementStructure, divisionAlgebra] THEN ERROR CantInvertError[elementStructure];
result ← Create[arg.class];
resultData ← NARROW[result.data];
FOR i: NAT IN [1..matrixStructureData.nRows] DO
FOR j: NAT IN [1..matrixStructureData.nCols] DO
MakeAij[arg, Aij, i, j];
tempVal ← matrixStructureData.elementStructure.class.divide[Determinant[Aij], det]; -- error will be generated here if Determinant[Aij] = 0
tempVal ← AC.ApplyLkpRecastObject[$fraction, matrixStructureData.elementStructure, LIST[Determinant[Aij], det ] ]; -- error will be generated here if Determinant[Aij] = 0
IF NOT Even[i + j] THEN tempVal ← matrixStructureData.elementStructure.class.negate[tempVal];
IF NOT Even[i + j] THEN tempVal ← AC.ApplyLkpRecastObject[$negation, matrixStructureData.elementStructure, LIST[tempVal ] ];
resultData[j][i] ← tempVal;
ENDLOOP;
ENDLOOP;
RETURN[result];
};
TryInvert: PUBLIC AC.UnaryOp ~ {
Trial inversion routine for squareMatricesOverRingClass
oldMatrixStructure: Object ← arg.class;
matrixStructureData: MatrixStructureData ← NARROW[oldMatrixStructure.data];
elementStructure: Object ← matrixStructureData.elementStructure;
superClass: Object ← elementStructure.class.class;
WHILE superClass#NIL DO
IF superClass.flavor#Structure THEN ERROR CantInvertError[oldMatrixStructure];
IF NOT AC.IsCategory[superClass, field] AND NOT AC.IsCategory[superClass, divisionAlgebra] THEN
{superClass ← superClass.class.class; LOOP}
ELSE {
newMatrixStructure: Object ← MakeMatrixStructure[superClass, matrixStructureData.nRows, matrixStructureData.nCols];
recastArg: Object ← Recast[arg, newMatrixStructure];
RETURN[Invert[recastArg] ];
};
ENDLOOP;
RETURN[ERROR CantInvertError[oldMatrixStructure] ];
};
Divide: PUBLIC AC.BinaryOp~ {
RETURN[ Multiply[ firstArg, Invert[ secondArg] ] ];
};
ScalarMultiply: PUBLIC AC.BinaryOp~ {
firstArg is scalar, secondArg is matrix
does its own recasting, so when installed, supplies DesiredArgStructures ← NIL
secondData: MatrixData ← NARROW[secondArg.data];
resultData: MatrixData;
matrixStructureData: MatrixStructureData ← NARROW[secondArg.class.data];
elementStructure: Object ← matrixStructureData.elementStructure;
recastMethod: Method ← AC.LookupMethodInStructure[$recast, elementStructure];
scalar: Object ← AC.ApplyNoLkpNoRecastObject[recastMethod, LIST[firstArg, elementStructure] ];
productMethod: Method ← AC.LookupMethodInStructure[$product, elementStructure];
IF scalar = NIL THEN TypeError[]; -- recast failed
result ← Create[secondArg.class];
resultData ← NARROW[result.data];
FOR i:NAT IN [1..matrixStructureData.nRows] DO
FOR j:NAT IN [1..matrixStructureData.nCols] DO
resultData[i][j] ← AC.ApplyNoLkpNoRecastObject[productMethod, LIST[scalar, secondData[i][j] ] ];
ENDLOOP;
ENDLOOP;
};
Transpose: PUBLIC AC.UnaryOp ~ {
argData: MatrixData ← NARROW[arg.data];
argStructureData: MatrixStructureData ← NARROW[arg.class.data];
resultData: MatrixData;
nRows: NAT ← argStructureData.nRows;
nCols: NAT ← argStructureData.nCols;
resultStructure: Object ← MakeMatrixStructure[argStructureData.elementStructure, nCols, nRows];
result ← Create[resultStructure];
resultData ← NARROW[result.data];
FOR i: INTEGER IN [1..argStructureData.nRows] DO
FOR j: INTEGER IN [1..argStructureData.nCols] DO
resultData[j][i] ← argData[i][j];
ENDLOOP;
ENDLOOP;
RETURN[result];
};
MakeAij: PROC[a, Aij: Matrix, i,j: INTEGER] = { --assume Aij is well formed
n, m: NAT ← 1; --row index, column index for new matrix
firstData: MatrixData ← NARROW[a.data];
secondData: MatrixData ← NARROW[Aij.data];
matrixStructureData: MatrixStructureData ← NARROW[a.class.data];
FOR row: NAT IN [1..matrixStructureData.nRows] DO
IF row = i THEN LOOP; --row index for original matrix
m ← 1;
FOR col: NAT IN [1..matrixStructureData.nCols] DO --column index for original matrix
IF col = j THEN LOOP;
secondData[n][m] ← firstData[row][col];
m ← m+1;
ENDLOOP;
n ← n+1;
ENDLOOP;
};
Even: PROC[I: NAT] RETURNS[BOOL] = {RETURN[I MOD 2 = 0]};
Determinant: PUBLIC AC.StructuredToGroundOp ~ {
Expansion by minors.
argData: MatrixData ← NARROW[structuredElt.data];
matrixStructureData: MatrixStructureData ← NARROW[structuredElt.class.data];
IF matrixStructureData.nRows = 1 THEN groundElement ← argData[1][1]
ELSE IF matrixStructureData.nRows=2 THEN groundElement ← matrixStructureData.elementStructure.class.subtract[
matrixStructureData.elementStructure.class.multiply[argData[1][1], argData[2][2] ],
matrixStructureData.elementStructure.class.multiply[argData[1][2], argData[2][1] ] ]
ELSE IF matrixStructureData.nRows=2 THEN {
prod1: Object ← AC.ApplyLkpRecastObject[$product, matrixStructureData.elementStructure, LIST[argData[1][1], argData[2][2] ] ];
prod2: Object ← AC.ApplyLkpRecastObject[$product, matrixStructureData.elementStructure, LIST[argData[1][2], argData[2][1] ] ];
groundElement ← AC.ApplyLkpRecastObject[$difference, matrixStructureData.elementStructure, LIST[prod1, prod2 ] ];
}
ELSE {
i, j: NAT;
zero: Object ← AC.ApplyLkpRecastObject[$zero, matrixStructureData.elementStructure, LIST[matrixStructureData.elementStructure] ];
scratchMatrixRing: Object ← MakeMatrixStructure[matrixStructureData.elementStructure, matrixStructureData.nRows - 1, matrixStructureData.nCols - 1];
Aij: Matrix ← Create[scratchMatrixRing];
groundElement ← zero;
j ← 1; --always use column 1 for now
FOR i IN [1..matrixStructureData.nRows] DO
tempVal: Object;
MakeAij[structuredElt, Aij, i, j];
tempVal ← matrixStructureData.elementStructure.class.multiply[ argData[i][j], Determinant[Aij] ];
tempVal ← AC.ApplyLkpRecastObject[$product, matrixStructureData.elementStructure, LIST[argData[i][j], Determinant[Aij] ] ];
IF NOT Even[i + j] THEN tempVal ← matrixStructureData.elementStructure.class.negate[tempVal];
IF NOT Even[i + j] THEN tempVal ← AC.ApplyLkpRecastObject[$negation, matrixStructureData.elementStructure, LIST[tempVal ] ];
groundElement ← matrixStructureData.elementStructure.class.add[groundElement, tempVal];
groundElement ← AC.ApplyLkpRecastObject[$sum, matrixStructureData.elementStructure, LIST[groundElement, tempVal ] ];
ENDLOOP;
};
RETURN[groundElement];
};
Comparison
Equal: PUBLIC AC.BinaryPredicate ~ {
firstData: MatrixData ← NARROW[firstArg.data];
secondData: MatrixData ← NARROW[secondArg.data];
matrixStructureData: MatrixStructureData ← NARROW[firstArg.class.data];
elementEqualsMethod: Method ← AC.LookupMethodInStructure[$eqFormula, matrixStructureData.elementStructure];
FOR i: INTEGER IN [1..matrixStructureData.nRows] DO
FOR j: INTEGER IN [1..matrixStructureData.nCols] DO
IF NOT AC.ApplyPredNoLkpNoRecast[elementEqualsMethod, LIST[firstData[i][j], secondData[i][j] ] ] THEN RETURN[FALSE];
ENDLOOP;
ENDLOOP;
RETURN[TRUE];
};
Standard Desired Arg Structures
ObjectAndIntDesired: PUBLIC AC.UnaryToListOp ~ {
RETURN[ LIST[arg, Ints.Ints] ]; -- arg assumed to be a Matrix Structure
};
MatrixAndElementDesired: PUBLIC AC.UnaryToListOp ~ {
thisMatrixStructureData: MatrixStructureData ← NARROW[arg.data];
elementStructure: Object ← thisMatrixStructureData.elementStructure;
RETURN[ LIST[arg, elementStructure] ]; -- arg assumed to be a Matrix Structure
};
ElementAndMatrixDesired: PUBLIC AC.UnaryToListOp ~ {
thisMatrixStructureData: MatrixStructureData ← NARROW[arg.data];
elementStructure: Object ← thisMatrixStructureData.elementStructure;
RETURN[ LIST[elementStructure, arg] ]; -- arg assumed to be a Matrix Structure
};
MatrixIntAndElementDesired: PUBLIC AC.UnaryToListOp ~ {
thisMatrixStructureData: MatrixStructureData ← NARROW[arg.data];
elementStructure: Object ← thisMatrixStructureData.elementStructure;
RETURN[ LIST[arg, Ints.Ints, elementStructure] ]; -- arg assumed to be a Matrix Structure
};
Start Code
matrixRingClass: PUBLIC ObjectClass ← NEW[ObjectClassRec ← [
category: ring,
multiply: Multiply,
commutative: FALSE,
integralDomain: FALSE,
gcdDomain: FALSE,
euclideanDomain: FALSE,
] ];
matrixAlgebraClass: PUBLIC ObjectClass ← NEW[ObjectClassRec ← [
category: algebra,
multiply: Multiply,
commutative: FALSE,
invert: Invert, -- it is up to the user not to try to invert a singular matrix
divide: Divide,
scalarMultiply: ScalarMultiply,
integralDomain: FALSE,
gcdDomain: FALSE,
euclideanDomain: FALSE,
] ];
matrixClass: AC.Object ← AC.MakeClass["MatrixClass", NIL, NIL];
squareMatricesOverRingClass: AC.Object ← AC.MakeClass["SquareMatricesOverRingClass", matrixClass, NIL]; -- squareMatricesOverRingClass is subclass[matrixClass]
squareMatricesOverFieldClass: AC.Object ← AC.MakeClass["SquareMatricesOverFieldClass", squareMatricesOverRingClass, NIL]; -- squareMatricesOverFieldClass is subclass[squareMatricesOverRingClass]
algebraCategoryMethod: Method ← AC.MakeMethod[Value, FALSE, NEW[AC.Category ← algebra], NIL, NIL];
ringCategoryMethod: Method ← AC.MakeMethod[Value, FALSE, NEW[AC.Category ← ring], NIL, NIL];
setCategoryMethod: Method ← AC.MakeMethod[Value, FALSE, NEW[AC.Category ← set], NIL, NIL];
matrixStructureMethod: Method ← AC.MakeMethod[Value, FALSE, NIL, NIL, "matrixStructure"];
shortPrintNameMethod: Method ← AC.MakeMethod[ToRopeOp, FALSE, NEW[AC.ToRopeOp ← ShortPrintName], NIL, "shortPrintName"];
recastMethod: Method ← AC.MakeMethod[BinaryOp, TRUE, NEW[AC.BinaryOp ← Recast], NIL, "recast"];
canRecastMethod: Method ← AC.MakeMethod[BinaryPredicate, TRUE, NEW[AC.BinaryPredicate ← CanRecast], NIL, "canRecast"];
toExprMethod: Method ← AC.MakeMethod[ToExprOp, FALSE, NEW[AC.ToExprOp ← ToExpr], NEW[AC.UnaryToListOp ← AC.DefaultDesiredArgStructures], "toExpr"];
legalFirstCharMethod: Method ← AC.MakeMethod[LegalFirstCharOp, FALSE, NEW[AC.LegalFirstCharOp ← LegalFirstChar], NIL, "legalFirstChar"];
readMethod: Method ← AC.MakeMethod[ReadOp, FALSE, NEW[AC.ReadOp ← Read], NIL, "read"];
fromRopeMethod: Method ← AC.MakeMethod[FromRopeOp, TRUE, NEW[AC.FromRopeOp ← FromRope], NIL, "fromRope"];
toRopeMethod: Method ← AC.MakeMethod[ToRopeOp, FALSE, NEW[AC.ToRopeOp ← ToRope], NIL, "toRope"];
diagonalMatrixMethod: Method ← AC.MakeMethod[UnaryImbedOp, TRUE, NEW[AC.UnaryImbedOp ← DiagonalMatrix], NIL, "diagonalMatrix"]; -- 4/17/87 recast procs need to be able to handle Structure args; could then have non—NIL DesiredArgStructures
matrixMethod: Method ← AC.MakeMethod[MatrixImbedOp, FALSE, NEW[AC.MatrixImbedOp ← MakeMatrix], NEW[AC.UnaryToListOp ← AC.DefaultDesiredArgStructures], "matrix"];
parenMethod: Method ← AC.MakeMethod[UnaryOp, FALSE, NEW[AC.UnaryOp ← AC.Copy], NIL, "paren"];
zeroMethod: Method ← AC.MakeMethod[NullaryOp, FALSE, NEW[AC.NullaryOp ← Zero], NIL, "zero"];
oneMethod: Method ← AC.MakeMethod[NullaryOp, FALSE, NEW[AC.NullaryOp ← One], NIL, "one"];
sumMethod: Method ← AC.MakeMethod[BinaryOp, TRUE, NEW[AC.BinaryOp ← Add], NEW[AC.UnaryToListOp ← AC.DefaultDesiredArgStructures], "sum"];
negationMethod: Method ← AC.MakeMethod[UnaryOp, TRUE, NEW[AC.UnaryOp ← Negate], NEW[AC.UnaryToListOp ← AC.DefaultDesiredArgStructures], "negation"];
differenceMethod: Method ← AC.MakeMethod[BinaryOp, TRUE, NEW[AC.BinaryOp ← Subtract], NEW[AC.UnaryToListOp ← AC.DefaultDesiredArgStructures], "difference"];
productMethod: Method ← AC.MakeMethod[BinaryOp, TRUE, NEW[AC.BinaryOp ← Multiply], NEW[AC.UnaryToListOp ← AC.DefaultDesiredArgStructures], "product"];
powerMethod: Method ← AC.MakeMethod[BinaryOp, TRUE, NEW[AC.BinaryOp ← Power], NEW[AC.UnaryToListOp ← ObjectAndIntDesired], "power"];
invertMethod: Method ← AC.MakeMethod[UnaryOp, TRUE, NEW[AC.UnaryOp ← Invert], NEW[AC.UnaryToListOp ← AC.DefaultDesiredArgStructures], "invert"];
tryInvertMethod: Method ← AC.MakeMethod[UnaryOp, TRUE, NEW[AC.UnaryOp ← TryInvert], NEW[AC.UnaryToListOp ← AC.DefaultDesiredArgStructures], "invert"];
fractionMethod: Method ← AC.MakeMethod[BinaryOp, TRUE, NEW[AC.BinaryOp ← Divide], NEW[AC.UnaryToListOp ← AC.DefaultDesiredArgStructures], "divide"];
scalarProductMethod: Method ← AC.MakeMethod[BinaryOp, TRUE, NEW[AC.BinaryOp ← ScalarMultiply], NEW[AC.UnaryToListOp ← ElementAndMatrixDesired], "scalarProduct"];
equalMethod: Method ← AC.MakeMethod[BinaryPredicate, TRUE, NEW[AC.BinaryPredicate ← Equal], NEW[AC.UnaryToListOp ← AC.DefaultDesiredArgStructures], "equals"];
transposeMethod: Method ← AC.MakeMethod[UnaryOp, TRUE, NEW[AC.UnaryOp ← Transpose], NEW[AC.UnaryToListOp ← AC.DefaultDesiredArgStructures], "transpose"];
determinantMethod: Method ← AC.MakeMethod[StructuredToGroundOp, TRUE, NEW[AC.StructuredToGroundOp ← Determinant], NEW[AC.UnaryToListOp ← AC.DefaultDesiredArgStructures], "determinant"];
makeMatrixStructureMethod: Method ← AC.MakeMethod[MatrixStructureConstructor, FALSE, NEW[AC.MatrixStructureConstructor ← MakeMatrixStructure], NIL, "makeMatrixStructure"];
AC.AddMethodToClass[$matrixStructure, matrixStructureMethod, matrixClass];
AC.AddMethodToClass[$category, setCategoryMethod, matrixClass];
AC.AddMethodToClass[$shortPrintName, shortPrintNameMethod, matrixClass];
AC.AddMethodToClass[$recast, recastMethod, matrixClass];
AC.AddMethodToClass[$canRecast, canRecastMethod, matrixClass];
AC.AddMethodToClass[$toExpr, toExprMethod, matrixClass];
AC.AddMethodToClass[$legalFirstChar, legalFirstCharMethod, matrixClass];
AC.AddMethodToClass[$read, readMethod, matrixClass];
AC.AddMethodToClass[$fromRope, fromRopeMethod, matrixClass];
AC.AddMethodToClass[$toRope, toRopeMethod, matrixClass];
AC.AddMethodToClass[$matrix, matrixMethod, matrixClass];
AC.AddMethodToClass[$paren, parenMethod, matrixClass];
AC.AddMethodToClass[$zero, zeroMethod, matrixClass];
AC.AddMethodToClass[$sum, sumMethod, matrixClass];
AC.AddMethodToClass[$negation, negationMethod, matrixClass];
AC.AddMethodToClass[$difference, differenceMethod, matrixClass];
AC.AddMethodToClass[$transp, transposeMethod, matrixClass];
AC.AddMethodToClass[$eqFormula, equalMethod, matrixClass];
AC.AddMethodToClass[$category, ringCategoryMethod, squareMatricesOverRingClass];
AC.AddMethodToClass[$one, oneMethod, squareMatricesOverRingClass];
AC.AddMethodToClass[$diagonalMatrix, diagonalMatrixMethod, squareMatricesOverRingClass];
AC.AddMethodToClass[$product, productMethod, squareMatricesOverRingClass];
AC.AddMethodToClass[$pow, powerMethod, squareMatricesOverRingClass];
AC.AddMethodToClass[$det, determinantMethod, squareMatricesOverRingClass];
AC.AddMethodToClass[$invert, tryInvertMethod, squareMatricesOverRingClass];
AC.AddMethodToClass[$scalarProduct, scalarProductMethod, squareMatricesOverRingClass];
AC.AddMethodToClass[$category, algebraCategoryMethod, squareMatricesOverFieldClass];
AC.AddMethodToClass[$invert, invertMethod, squareMatricesOverFieldClass];
AC.AddMethodToClass[$fraction, fractionMethod, squareMatricesOverFieldClass];
AC.AddMethodToClass[$makeMatrixStructure, makeMatrixStructureMethod, Structures.StructuresClass];
END.