PolynomialsImpl.mesa
Last Edited by: Arnon, June 10, 1985 4:19:22 pm PDT
DIRECTORY
Rope,
Basics,
Atom,
Ascii,
Convert,
IO,
Points,
AlgebraClasses,
Variables,
DistribPolys,
Matrices,
MathExpr,
MathConstructors,
Polynomials;
PolynomialsImpl: CEDAR PROGRAM
IMPORTS Rope, Convert, IO, Atom, Points, AlgebraClasses, Variables, DistribPolys, Matrices, MathConstructors
EXPORTS Polynomials =
BEGIN OPEN AC: AlgebraClasses, VARS: Variables, DP: DistribPolys, MAT: Matrices, Polynomials;
TypeError: PUBLIC ERROR [message: ATOM ← $Unspecified] = CODE;
BadElementStructure: PUBLIC ERROR [message: ATOM ← $Unspecified] = CODE;
OperationError: PUBLIC ERROR [message: ATOM ← $Unspecified] = CODE;
Classes for Polynomial Rings
ClassPrintName: AC.PrintNameProc = {
data: PolynomialRingData ← NARROW[structure.instanceData];
RETURN[Rope.Cat[
"Polynomials in ",
VARS.VariableSeqToRope[data.variable],
" over ",
data.coeffRing.class.printName[data.coeffRing]
] ];
};
ClassCharacteristic: AC.StructureRankOp = {
data: PolynomialRingData ← NARROW[structure.instanceData];
RETURN[ data.coeffRing.class.characteristic[data.coeffRing] ]
};
ClassLegalFirstChar: AC.LegalFirstCharOp = {
data: PolynomialRingData ← NARROW[structure.instanceData];
IF data.baseCoeffRing.class.legalFirstChar[char, data.baseCoeffRing] THEN RETURN[TRUE];
RETURN[VARS.VariableFirstChar[char, data.allVariables] ];
Legal first char of polynomial is either legal first char of coefficient, or first char of some variable (this all with reference to distributed rep)
};
ClassZero: AC.NullaryOp = {
RETURN[ NEW[ AC.ObjectRec ← [structure: structure, data: NIL]] ]
};
ClassOne: AC.NullaryOp = {
data: PolynomialRingData ← NARROW[structure.instanceData];
RETURN[ Monom[data.coeffRing.class.one[data.coeffRing], NEW[CARDINAL ← 0], structure] ];
};
polynomialOps: PolynomialOps ← NEW[PolynomialOpsRec ← [
monomial: Monom,
differentiate: Diff,
leadingCoefficient: LeadingCoeff,
degree: Deg,
mainVarEval: MainVarEv,
allVarEval: AllVarEv,
subst: Sub,
sylvesterMatrix: SylvMatrix,
resultant: Res
] ];
polynomialProp: Atom.DottedPair ← NEW[Atom.DottedPairNode← [$PolynomialRing, polynomialOps]];
polynomialsOverCommutativeOrderedRingClass: PUBLIC AC.StructureClass ← NEW[AC.StructureClassRec ← [
category: ring,
flavor: "Polynomials",
printName: ClassPrintName,
structureEqual: AC.defaultStructureEqualityTest,
characteristic: ClassCharacteristic,
isElementOf: AC.defaultElementOfProc,
legalFirstChar: ClassLegalFirstChar,
read: Read,
fromRope: FromRope,
toRope: ToRope,
write: Write,
toExpr: MakePolyExpr,
add: Add,
negate: Negate,
subtract: Subtract,
zero: ClassZero,
multiply: Multiply,
commutative: TRUE,
one: ClassOne,
equal: Equal,
ordered: TRUE,
sign: Sign,
abs: Abs,
compare: Compare,
integralDomain: TRUE, -- not necessarily accurate; need separate classrecs
gcdDomain: FALSE, -- not necessarily accurate; need separate classrecs
gcd: NIL, -- should have a poly gcd proc here when appropriate
euclideanDomain: FALSE,
propList: LIST[polynomialProp]
] ];
polynomialsOverCommutativeUnOrderedRingClass: PUBLIC AC.StructureClass ← NEW[AC.StructureClassRec ← [
category: ring,
flavor: "Polynomials",
printName: ClassPrintName,
structureEqual: AC.defaultStructureEqualityTest,
characteristic: ClassCharacteristic,
isElementOf: AC.defaultElementOfProc,
legalFirstChar: ClassLegalFirstChar,
read: Read,
fromRope: FromRope,
toRope: ToRope,
write: Write,
toExpr: MakePolyExpr,
add: Add,
negate: Negate,
subtract: Subtract,
zero: ClassZero,
multiply: Multiply,
commutative: TRUE,
one: ClassOne,
equal: Equal,
ordered: FALSE,
integralDomain: TRUE, -- not necessarily accurate; need separate classrecs
gcdDomain: FALSE, -- not necessarily accurate; need separate classrecs
gcd: NIL, -- should have a poly gcd proc here when appropriate
euclideanDomain: FALSE,
propList: LIST[polynomialProp]
] ];
polynomialsOverNonCommutativeOrderedRingClass: PUBLIC AC.StructureClass ← NEW[AC.StructureClassRec ← [
category: ring,
flavor: "Polynomials",
printName: ClassPrintName,
structureEqual: AC.defaultStructureEqualityTest,
characteristic: ClassCharacteristic,
isElementOf: AC.defaultElementOfProc,
legalFirstChar: ClassLegalFirstChar,
read: Read,
fromRope: FromRope,
toRope: ToRope,
write: Write,
toExpr: MakePolyExpr,
add: Add,
negate: Negate,
subtract: Subtract,
zero: ClassZero,
multiply: Multiply,
commutative: FALSE,
one: ClassOne,
equal: Equal,
ordered: TRUE,
sign: Sign,
abs: Abs,
compare: Compare,
integralDomain: FALSE, -- integralDomain's commutative by definition
gcdDomain: FALSE, -- not necessarily accurate; need separate classrecs
gcd: NIL, -- should have a poly gcd proc here when appropriate
euclideanDomain: FALSE,
propList: LIST[polynomialProp]
] ];
polynomialsOverNonCommutativeUnOrderedRingClass: PUBLIC AC.StructureClass ← NEW[AC.StructureClassRec ← [
category: ring,
flavor: "Polynomials",
printName: ClassPrintName,
structureEqual: AC.defaultStructureEqualityTest,
characteristic: ClassCharacteristic,
isElementOf: AC.defaultElementOfProc,
legalFirstChar: ClassLegalFirstChar,
read: Read,
fromRope: FromRope,
toRope: ToRope,
write: Write,
toExpr: MakePolyExpr,
add: Add,
negate: Negate,
subtract: Subtract,
zero: ClassZero,
multiply: Multiply,
commutative: FALSE,
one: ClassOne,
equal: Equal,
ordered: FALSE,
integralDomain: FALSE, -- integralDomain's commutative by definition
gcdDomain: FALSE, -- not necessarily accurate; need separate classrecs
gcd: NIL, -- should have a poly gcd proc here when appropriate
euclideanDomain: FALSE,
propList: LIST[polynomialProp]
] ];
polynomialsOverCommutativeOrderedFieldClass: PUBLIC AC.StructureClass ← NEW[AC.StructureClassRec ← [
category: ring,
flavor: "Polynomials",
printName: ClassPrintName,
structureEqual: AC.defaultStructureEqualityTest,
characteristic: ClassCharacteristic,
isElementOf: AC.defaultElementOfProc,
legalFirstChar: ClassLegalFirstChar,
read: Read,
fromRope: FromRope,
toRope: ToRope,
write: Write,
toExpr: MakePolyExpr,
add: Add,
negate: Negate,
subtract: Subtract,
zero: ClassZero,
multiply: Multiply,
commutative: TRUE,
one: ClassOne,
equal: Equal,
ordered: TRUE,
sign: Sign,
abs: Abs,
compare: Compare,
integralDomain: TRUE, -- integralDomain's commutative by definition
gcdDomain: FALSE, -- not necessarily accurate; need separate classrecs
gcd: NIL, -- should have a poly gcd proc here when appropriate
euclideanDomain: TRUE,
remainder: Remainder, -- euclideanDomains only
propList: LIST[polynomialProp]
] ];
polynomialsOverCommutativeUnOrderedFieldClass: PUBLIC AC.StructureClass ← NEW[AC.StructureClassRec ← [
category: ring,
flavor: "Polynomials",
printName: ClassPrintName,
structureEqual: AC.defaultStructureEqualityTest,
characteristic: ClassCharacteristic,
isElementOf: AC.defaultElementOfProc,
legalFirstChar: ClassLegalFirstChar,
read: Read,
fromRope: FromRope,
toRope: ToRope,
write: Write,
toExpr: MakePolyExpr,
add: Add,
negate: Negate,
subtract: Subtract,
zero: ClassZero,
multiply: Multiply,
commutative: TRUE,
one: ClassOne,
equal: Equal,
ordered: FALSE,
integralDomain: TRUE, -- integralDomain's commutative by definition
gcdDomain: FALSE, -- not necessarily accurate; need separate classrecs
gcd: NIL, -- should have a poly gcd proc here when appropriate
euclideanDomain: TRUE,
remainder: Remainder, -- euclideanDomains only
propList: LIST[polynomialProp]
] ];
polynomialsOverNonCommutativeOrderedFieldClass: PUBLIC AC.StructureClass ← NEW[AC.StructureClassRec ← [
category: ring,
flavor: "Polynomials",
printName: ClassPrintName,
structureEqual: AC.defaultStructureEqualityTest,
characteristic: ClassCharacteristic,
isElementOf: AC.defaultElementOfProc,
legalFirstChar: ClassLegalFirstChar,
read: Read,
fromRope: FromRope,
toRope: ToRope,
write: Write,
toExpr: MakePolyExpr,
add: Add,
negate: Negate,
subtract: Subtract,
zero: ClassZero,
multiply: Multiply,
commutative: FALSE,
one: ClassOne,
equal: Equal,
ordered: TRUE,
sign: Sign,
abs: Abs,
compare: Compare,
integralDomain: FALSE, -- integralDomain's commutative by definition
gcdDomain: FALSE, -- not necessarily accurate; need separate classrecs
gcd: NIL, -- should have a poly gcd proc here when appropriate
euclideanDomain: TRUE,
remainder: Remainder, -- euclideanDomains only
propList: LIST[polynomialProp]
] ];
polynomialsOverNonCommutativeUnOrderedFieldClass: PUBLIC AC.StructureClass ← NEW[AC.StructureClassRec ← [
category: ring,
flavor: "Polynomials",
printName: ClassPrintName,
structureEqual: AC.defaultStructureEqualityTest,
characteristic: ClassCharacteristic,
isElementOf: AC.defaultElementOfProc,
legalFirstChar: ClassLegalFirstChar,
read: Read,
fromRope: FromRope,
toRope: ToRope,
write: Write,
toExpr: MakePolyExpr,
add: Add,
negate: Negate,
subtract: Subtract,
zero: ClassZero,
multiply: Multiply,
commutative: FALSE,
one: ClassOne,
equal: Equal,
ordered: FALSE,
integralDomain: FALSE, -- integralDomain's commutative by definition
gcdDomain: FALSE, -- not necessarily accurate; need separate classrecs
gcd: NIL, -- should have a poly gcd proc here when appropriate
euclideanDomain: TRUE,
remainder: Remainder, -- euclideanDomains only
propList: LIST[polynomialProp]
] ];
Polynomial Ring Constructors
MakePolynomialStructure: PUBLIC PROC [coeffRing: AC.Structure, V: VARS.VariableSeq] RETURNS [polynomialRing: AC.Structure] ~ {
baseCoeffRing, strictCoeffRing: AC.Structure;
allVariables, strictVariable: VARS.VariableSeq;
polynomialRingData: PolynomialRingData;
IF V.lengthPlus1 - 1 > 1 THEN {
strictCoeffRing ← MakePolynomialStructure[coeffRing, VARS.RemoveMainVariable[V] ];
strictVariable ← VARS.MainVariable[V];
}
ELSE {
strictCoeffRing ← coeffRing;
strictVariable ← V;
};
IF IsPolynomialRing[strictCoeffRing] THEN { -- polynomials over polynomials, or not?
data: PolynomialRingData ← NARROW[strictCoeffRing.instanceData];
baseCoeffRing ← data.baseCoeffRing;
allVariables ← VARS.AddMainVariable[data.allVariables, strictVariable]
}
ELSE {
baseCoeffRing ← strictCoeffRing;
allVariables ← strictVariable;
};
polynomialRingData ← NEW[PolynomialRingDataRec ← [
coeffRing: strictCoeffRing,
variable: strictVariable,
baseCoeffRing: baseCoeffRing,
allVariables: allVariables
] ];
SELECT strictCoeffRing.class.category FROM
ring, algebra =>
IF strictCoeffRing.class.commutative THEN {
IF strictCoeffRing.class.ordered THEN
RETURN[ NEW[AC.StructureRec ← [
class: polynomialsOverCommutativeOrderedRingClass,
instanceData: polynomialRingData
] ] ]
ELSE
RETURN[ NEW[AC.StructureRec ← [
class: polynomialsOverCommutativeUnOrderedRingClass,
instanceData: polynomialRingData
] ] ]
}
ELSE {
IF strictCoeffRing.class.ordered THEN
RETURN[ NEW[AC.StructureRec ← [
class: polynomialsOverNonCommutativeOrderedRingClass,
instanceData: polynomialRingData
] ] ]
ELSE
RETURN[ NEW[AC.StructureRec ← [
class: polynomialsOverNonCommutativeUnOrderedRingClass,
instanceData: polynomialRingData
] ] ]
};
field, divisionAlgebra =>
IF strictCoeffRing.class.commutative THEN {
IF strictCoeffRing.class.ordered THEN
RETURN[ NEW[AC.StructureRec ← [
class: polynomialsOverCommutativeOrderedFieldClass,
instanceData: polynomialRingData
] ] ]
ELSE
RETURN[ NEW[AC.StructureRec ← [
class: polynomialsOverCommutativeUnOrderedFieldClass,
instanceData: polynomialRingData
] ] ]
}
ELSE {
IF strictCoeffRing.class.ordered THEN
RETURN[ NEW[AC.StructureRec ← [
class: polynomialsOverNonCommutativeOrderedFieldClass,
instanceData: polynomialRingData
] ] ]
ELSE
RETURN[ NEW[AC.StructureRec ← [
class: polynomialsOverNonCommutativeUnOrderedFieldClass,
instanceData: polynomialRingData
] ] ]
};
ENDCASE => ERROR BadElementStructure[];
};
Extract Polynomial Operations from Class Property Lists
IsPolynomialRing: PUBLIC PROC [structure: AC.Structure] RETURNS [BOOL] ~ {
IF Atom.GetPropFromList[structure.class.propList, $PolynomialRing] # NIL THEN RETURN[TRUE] ELSE RETURN[FALSE];
};
Monomial: PUBLIC PROC [polynomialRing: AC.Structure] RETURNS [AC.BinaryImbedOp] ~ {
polynomialOps: PolynomialOps ← NARROW[ Atom.GetPropFromList[polynomialRing.class.propList, $PolynomialRing] ];
IF polynomialOps = NIL THEN ERROR;
RETURN[polynomialOps.monomial];
};
Differentiate: PUBLIC PROC [polynomialRing: AC.Structure] RETURNS [AC.UnaryOp] ~ {
polynomialOps: PolynomialOps ← NARROW[ Atom.GetPropFromList[polynomialRing.class.propList, $PolynomialRing] ];
IF polynomialOps = NIL THEN ERROR;
RETURN[polynomialOps.differentiate];
};
LeadingCoefficient: PUBLIC PROC [polynomialRing: AC.Structure] RETURNS [AC.UnaryOp] ~ {
polynomialOps: PolynomialOps ← NARROW[ Atom.GetPropFromList[polynomialRing.class.propList, $PolynomialRing] ];
IF polynomialOps = NIL THEN ERROR;
RETURN[polynomialOps.leadingCoefficient];
};
Degree: PUBLIC PROC [polynomialRing: AC.Structure] RETURNS [AC.ElementRankOp] ~ {
polynomialOps: PolynomialOps ← NARROW[ Atom.GetPropFromList[polynomialRing.class.propList, $PolynomialRing] ];
IF polynomialOps = NIL THEN ERROR;
RETURN[polynomialOps.degree];
};
MainVarEval: PUBLIC PROC [polynomialRing: AC.Structure] RETURNS [AC.BinaryOp] ~ {
polynomialOps: PolynomialOps ← NARROW[ Atom.GetPropFromList[polynomialRing.class.propList, $PolynomialRing] ];
IF polynomialOps = NIL THEN ERROR;
RETURN[polynomialOps.mainVarEval];
};
AllVarEval: PUBLIC PROC [polynomialRing: AC.Structure] RETURNS [AC.BinaryOp] ~ {
polynomialOps: PolynomialOps ← NARROW[ Atom.GetPropFromList[polynomialRing.class.propList, $PolynomialRing] ];
IF polynomialOps = NIL THEN ERROR;
RETURN[polynomialOps.allVarEval];
};
Subst: PUBLIC PROC [polynomialRing: AC.Structure] RETURNS [AC.BinaryOp] ~ {
polynomialOps: PolynomialOps ← NARROW[ Atom.GetPropFromList[polynomialRing.class.propList, $PolynomialRing] ];
IF polynomialOps = NIL THEN ERROR;
RETURN[polynomialOps.subst];
};
SylvesterMatrix: PUBLIC PROC [polynomialRing: AC.Structure] RETURNS [AC.BinaryOp] ~ {
polynomialOps: PolynomialOps ← NARROW[ Atom.GetPropFromList[polynomialRing.class.propList, $PolynomialRing] ];
IF polynomialOps = NIL THEN ERROR;
RETURN[polynomialOps.sylvesterMatrix];
};
Resultant: PUBLIC PROC [polynomialRing: AC.Structure] RETURNS [AC.BinaryOp] ~ {
polynomialOps: PolynomialOps ← NARROW[ Atom.GetPropFromList[polynomialRing.class.propList, $PolynomialRing] ];
IF polynomialOps = NIL THEN ERROR;
RETURN[polynomialOps.resultant];
};
Constructors
Monom: PUBLIC AC.BinaryImbedOp = {
structureData: PolynomialRingData ← NARROW[structure.instanceData];
refExp: REF CARDINAL ← NARROW[data2];
exp: CARDINAL ← refExp^;
IF NOT structureData.coeffRing.class.isElementOf[data1, structureData.coeffRing] THEN TypeError[]; -- check that we really are constructing an element of polynomialRing
IF structureData.coeffRing.class.equal[data1, structureData.coeffRing.class.zero[structureData.coeffRing]] THEN -- check for zero coeff
RETURN[ NEW[AC.ObjectRec ← [
structure: structure,
data: NIL
] ] ]
ELSE
RETURN[ NEW[AC.ObjectRec ← [
structure: structure,
data: LIST[ NEW[TermRec ← [exponent: exp, coefficient: data1] ] ]
] ] ];
};
Conversion and IO
Read: PUBLIC AC.ReadOp = {
data: PolynomialRingData ← NARROW[structure.instanceData];
dPoly: DP.DPolynomial;
termChar: CHAR;
[dPoly, termChar] ← DP.ReadDPoly[in, data.allVariables, data.baseCoeffRing];
RETURN[ PolyFromDPoly[dPoly, structure] ];
};
FromRope: PUBLIC AC.FromRopeOp = {
stream: IO.STREAM ← IO.RIS[in];
out ← Read[stream, structure];
};
ToRope: PUBLIC AC.ToRopeOp ~ {
data: PolynomialRingData ← NARROW[in.structure.instanceData];
out ← DP.DPolyToRope[DPolyFromPoly[in], data.allVariables, data.baseCoeffRing, NIL];
};
Write: PUBLIC AC.WriteOp ~ {
IO.PutRope[ stream, ToRope[in] ]
};
MakePolyExpr: PUBLIC AC.ToExprOp = {
data: PolynomialRingData ← NARROW[in.structure.instanceData];
dIn: DP.DPolynomial ← DPolyFromPoly[in];
addend, sum: MathExpr.EXPR;
V: VARS.VariableSeq ← data.allVariables;
coeffRing: AC.Structure ← data.baseCoeffRing;
one: AC.Object ← coeffRing.class.one[coeffRing];
firstTerm: BOOL ← TRUE;
trivialMonomial, firstOccurringVar: BOOL;
thisVar: MathExpr.EXPR;
coeff, coeffAbs: AC.Object;
degreeVec: DP.DegreeVector;
coeffSign: Basics.Comparison;
exponent, index: CARDINAL;
CreateDTerm: PROC[firstTerm: BOOL] ~ {
[coeff, degreeVec] ← dIn.first;
IF coeffRing.class.ordered THEN {
coeffSign ← coeffRing.class.sign[coeff];
coeffAbs ← coeffRing.class.abs[coeff]
}
ELSE { -- for unordered coeffRing, act as though coeff is positive
coeffSign ← greater;
coeffAbs ← coeff;
};
trivialMonomial ← TRUE;
firstOccurringVar ← TRUE;
degreeVec ← DP.DVReverse[degreeVec];
WHILE degreeVec#NIL DO
trivialMonomial ← FALSE;
exponent ← degreeVec.first; degreeVec ← degreeVec.rest;
index ← degreeVec.first; degreeVec ← degreeVec.rest;
IF exponent>1 THEN
thisVar ← MathConstructors.MakePow[
MathConstructors.MakeVariable[V[index] ],
MathConstructors.MakeInt[Convert.RopeFromCard[exponent]]
]
ELSE
thisVar ← MathConstructors.MakeVariable[V[index] ];
IF firstOccurringVar THEN addend ← thisVar ELSE
addend ← MathConstructors.MakeProduct[addend, thisVar];
firstOccurringVar ← FALSE;
ENDLOOP;
IF NOT trivialMonomial THEN {
IF NOT coeffRing.class.equal[coeffAbs, one] THEN
addend ← MathConstructors.MakeProduct[coeffRing.class.toExpr[coeffAbs], addend]
}
ELSE addend ← coeffRing.class.toExpr[coeffAbs];
IF firstTerm AND coeffSign = less THEN
addend ← MathConstructors.MakeNegation[addend];
};
IF dIn = DP.ZeroDPoly THEN RETURN[MathConstructors.MakeInt["0"] ];
CreateDTerm[TRUE]; -- first term
sum ← addend;
dIn ← dIn.rest;
WHILE dIn # NIL DO
CreateDTerm[FALSE]; -- not first term
IF coeffSign = greater THEN
sum ← MathConstructors.MakeSum[sum, addend]
ELSE
sum ← MathConstructors.MakeDifference[sum, addend];
dIn ← dIn.rest;
ENDLOOP;
RETURN[sum];
};
PolyToRope: PUBLIC PROC [in: Polynomial, termRope: Rope.ROPE ← NIL] RETURNS [out: Rope.ROPE] = {
data: PolynomialRingData ← NARROW[in.structure.instanceData];
out ← DP.DPolyToRope[DPolyFromPoly[in], data.allVariables, data.baseCoeffRing, termRope];
};
WritePoly: PUBLIC PROC [in: Polynomial, out: IO.STREAM, termRope: Rope.ROPE ← NIL] = {
out.PutF["\n %g \n", IO.rope[PolyToRope[in, termRope]] ];
};
ReadPolySeq: PUBLIC PROC [in: IO.STREAM, polynomialRing: AC.Structure] RETURNS [seq: PolynomialSeq] ~ {
puncChar: CHAR;
nextP: Polynomial;
length: NAT ← 0;
ReadPSFail: PUBLIC ERROR [subclass: ATOM ← $Unspecified] = CODE;
pList, pListTail: LIST OF Polynomial ← NIL;
[]← in.SkipWhitespace[];
puncChar ← in.GetChar[];
IF puncChar # '( THEN ReadPSFail[$LeftParenExpected];
[]← in.SkipWhitespace[];
IF in.PeekChar[] = ') THEN {
puncChar ← in.GetChar[]; -- toss it
[]← in.SkipWhitespace[]; -- move past all trailing white space
RETURN[NEW[PolynomialSeqRec[0] ] ];
};
WHILE puncChar # ') DO
nextP ← Read[in, polynomialRing];
length ← length + 1;
IF pList=NIL THEN pList ← pListTail ←LIST[nextP] ELSE
{ pListTail.rest ← LIST[nextP]; pListTail ← pListTail.rest };
[]← in.SkipWhitespace[];
puncChar ← in.GetChar[]; -- toss it
IF puncChar # ', AND puncChar # ') THEN ReadPSFail[$CommaOrLPExpected];
ENDLOOP;
[]← in.SkipWhitespace[]; -- move past all trailing white space
seq ← NEW[PolynomialSeqRec[length]];
FOR i:NAT IN [1..length+1) DO
seq[i] ← pList.first;
pList ← pList.rest;
ENDLOOP;
};
PolySeqFromRope: PUBLIC PROC [in: Rope.ROPE, polynomialRing: AC.Structure] RETURNS [out: PolynomialSeq] ~ {
PSStream: IO.STREAM ← IO.RIS[in];
RETURN[ ReadPolySeq[ PSStream, polynomialRing ] ];
};
PolySeqToRope: PUBLIC PROC [in: PolynomialSeq] RETURNS [out: Rope.ROPE] ~ {
data: PolynomialRingData;
IF in.lengthPlus1 - 1 = 0 THEN { out ← "(\n)\n"; RETURN };
data ← NARROW[in[1].structure.instanceData];
out ← "(\n";
FOR i:NAT IN [1..in.lengthPlus1) DO
out ← Rope.Concat[out, DP.DPolyToRope[DPolyFromPoly[in[i]], data.allVariables, data.baseCoeffRing] ];
IF i < in.lengthPlus1 - 1 THEN
out ← Rope.Concat[ out, ",\n"]
ELSE
out ← Rope.Concat[ out, "\n"]
ENDLOOP;
out ← Rope.Concat[ out, ")\n" ];
};
WritePolySeq: PUBLIC PROC [in: PolynomialSeq, out: IO.STREAM] ~ {
PSRope: Rope.ROPE ← PolySeqToRope[in];
out.PutF["%g", IO.rope[PSRope] ];
};
PolyFromDPoly: PUBLIC PROC [in: DP.DPolynomial, polynomialRing: AC.Structure] RETURNS [out: Polynomial] = {
data: PolynomialRingData ← NARROW[polynomialRing.instanceData];
numVars: CARDINAL = data.allVariables.lengthPlus1 - 1;
termDegree: CARDINAL; -- degree in main variable of current (output) term
outTermDCoefficient: DP.DPolynomial; -- output term coefficient, still in DP rep
outTermCoefficient: AC.Object; -- output term coefficient, converted to P rep
outTerm: Term; -- completed output term
outTerms, outTermsPointer: LIST OF Term ← NIL;
IF in = DP.ZeroDPoly THEN RETURN[NARROW[polynomialRing.class.zero[polynomialRing],Polynomial]];
WHILE in#NIL DO
termDegree ← DP.DVDegree[in.first.degreeVector, numVars]; -- degree in main variable
IF numVars > 1 THEN {
outTermDCoefficient ← LIST[ [in.first.coefficient, DP.DVRemoveMainVariablePower[in.first.degreeVector, numVars] ] ];
in ← in.rest;
WHILE in#NIL AND DP.DVDegree[in.first.degreeVector, numVars] = termDegree DO
outTermDCoefficient ← CONS[ [in.first.coefficient, DP.DVRemoveMainVariablePower[in.first.degreeVector, numVars] ], outTermDCoefficient];
in ← in.rest;
ENDLOOP;
outTermDCoefficient ← DP.DPReverse[ outTermDCoefficient ];
outTermCoefficient ← PolyFromDPoly[ outTermDCoefficient, data.coeffRing ];
}
ELSE {
outTermCoefficient ← in.first.coefficient;
in ← in.rest;
};
outTerm ← NEW[TermRec ← [
exponent: termDegree,
coefficient: outTermCoefficient
] ];
IF outTerms # NIL THEN outTermsPointer ← outTermsPointer.rest ← LIST[outTerm]
ELSE outTerms ← outTermsPointer ← LIST[outTerm];
ENDLOOP;
RETURN[ NEW[ AC.ObjectRec ← [structure: polynomialRing, data: outTerms]] ]
};
DPolyFromPoly: PUBLIC PROC [in: Polynomial] RETURNS [out: DP.DPolynomial] ~ {
data: PolynomialRingData ← NARROW[in.structure.instanceData];
numVars: CARDINAL = data.allVariables.lengthPlus1 - 1;
inTerms: LIST OF Term ← NARROW[in.data];
inTermDegree: CARDINAL; -- degree in main variable of current (input) term
inTermDPolynomial: DP.DPolynomial;
singleVariableIndex: CARDINAL = 1;
ok: BOOL;
degreeVec: DP.DegreeVector;
IF ZeroPoly[in] THEN RETURN[DP.ZeroDPoly];
out ← DP.ZeroDPoly;
WHILE inTerms # NIL DO
inTermDegree ← inTerms.first.exponent;
IF numVars > 1 THEN {
inTermDPolynomial ← DPolyFromPoly[NARROW[inTerms.first.coefficient, Polynomial]];
WHILE inTermDPolynomial#NIL DO
out ← CONS[ [ inTermDPolynomial.first.coefficient, DP.DVAddMainVariablePower[inTermDPolynomial.first.degreeVector, numVars, inTermDegree] ], out];
inTermDPolynomial ← inTermDPolynomial.rest;
ENDLOOP;
}
ELSE {
[ok, degreeVec] ← DP.DVInsertVariablePower[singleVariableIndex, inTermDegree, NIL];
IF NOT ok THEN ERROR;
out ← CONS[ [inTerms.first.coefficient, degreeVec], out];
};
inTerms ← inTerms.rest;
ENDLOOP;
out ← DP.DPReverse[out];
};
Arithmetic
Add: PUBLIC AC.BinaryOp ~ {
data: PolynomialRingData ← NARROW[firstArg.structure.instanceData];
inTerms1: LIST OF Term ← NARROW[firstArg.data];
inTerms2: LIST OF Term ← NARROW[secondArg.data];
outTerms, outTermsPointer, tail: LIST OF Term ← NIL;
outTerm: Term;
termToAdd: BOOL;
IF ZeroPoly[firstArg] THEN RETURN[secondArg];
IF ZeroPoly[secondArg] THEN RETURN[firstArg];
WHILE inTerms1 # NIL AND inTerms2 # NIL DO
termToAdd ← TRUE;
SELECT inTerms1.first.exponent FROM
< inTerms2.first.exponent => {
outTerm ← inTerms2.first;
inTerms2 ← inTerms2.rest;
};
> inTerms2.first.exponent => {
outTerm ← inTerms1.first;
inTerms1 ← inTerms1.rest;
};
= inTerms2.first.exponent => {
coeff: AC.Object ← data.coeffRing.class.add[inTerms1.first.coefficient, inTerms2.first.coefficient];
IF NOT data.coeffRing.class.equal[coeff, data.coeffRing.class.zero[data.coeffRing] ] THEN
outTerm ← NEW[TermRec ← [
exponent: inTerms2.first.exponent,
coefficient: coeff
] ]
ELSE
termToAdd ← FALSE;
inTerms1 ← inTerms1.rest;
inTerms2 ← inTerms2.rest;
};
ENDCASE;
IF termToAdd THEN {
IF outTerms # NIL THEN outTermsPointer ← outTermsPointer.rest ← LIST[outTerm]
ELSE outTerms ← outTermsPointer ← LIST[outTerm];
};
ENDLOOP;
IF inTerms1 = NIL THEN tail ← inTerms2 ELSE tail ← inTerms1;
IF outTerms # NIL THEN outTermsPointer.rest ← tail
ELSE outTerms ← tail;
RETURN[ NEW[AC.ObjectRec ← [structure: firstArg.structure, data: outTerms] ] ];
};
Negate: PUBLIC AC.UnaryOp ~ {
data: PolynomialRingData ← NARROW[arg.structure.instanceData];
inTerms: LIST OF Term ← NARROW[arg.data];
outTerms, outTermsPointer: LIST OF Term ← NIL;
outTerm: Term;
IF ZeroPoly[arg] THEN RETURN[arg];
WHILE inTerms # NIL DO
outTerm ← NEW[TermRec ← [
exponent: inTerms.first.exponent,
coefficient: data.coeffRing.class.negate[inTerms.first.coefficient]
] ];
IF outTerms # NIL THEN outTermsPointer ← outTermsPointer.rest ← LIST[outTerm]
ELSE outTerms ← outTermsPointer ← LIST[outTerm];
inTerms ← inTerms.rest;
ENDLOOP;
RETURN[ NEW[AC.ObjectRec ← [structure: arg.structure, data: outTerms] ] ];
};
Subtract: PUBLIC AC.BinaryOp ~ {
RETURN[ Add[ firstArg, Negate[ secondArg] ] ];
};
Multiply: PUBLIC AC.BinaryOp ~ {
data: PolynomialRingData ← NARROW[firstArg.structure.instanceData];
inTerms1: LIST OF Term ← NARROW[firstArg.data];
inTerms2: LIST OF Term ← NARROW[secondArg.data];
outSummand: Polynomial;
outSummandTerms, outSummandTermsPointer: LIST OF Term;
outSummandTerm: Term;
scratchInTerms1: LIST OF Term;
result ← NARROW[firstArg.structure.class.zero[firstArg.structure], Polynomial];
IF ZeroPoly[firstArg] OR ZeroPoly[secondArg] THEN RETURN[result];
WHILE inTerms2 # NIL DO
outSummandTerms ← outSummandTermsPointer ← NIL;
scratchInTerms1 ← inTerms1;
WHILE scratchInTerms1 # NIL DO
coeff: AC.Object ← data.coeffRing.class.multiply[
scratchInTerms1.first.coefficient,
inTerms2.first.coefficient
];
outSummandTerm ← NEW[TermRec ← [
exponent: scratchInTerms1.first.exponent + inTerms2.first.exponent,
coefficient: coeff
] ];
IF outSummandTerms # NIL THEN outSummandTermsPointer ← outSummandTermsPointer.rest ← LIST[outSummandTerm]
ELSE outSummandTerms ← outSummandTermsPointer ← LIST[outSummandTerm];
scratchInTerms1 ← scratchInTerms1.rest;
ENDLOOP;
outSummand ← NEW[AC.ObjectRec ← [structure: firstArg.structure, data: outSummandTerms]];
result ← Add[result, outSummand];
inTerms2 ← inTerms2.rest;
ENDLOOP;
};
Remainder: PUBLIC AC.BinaryOp ~ {
newDividend: Polynomial ← firstArg;
data: PolynomialRingData ← NARROW[firstArg.structure.instanceData];
IF data.coeffRing.class.category # field AND data.coeffRing.class.category # divisionAlgebra THEN TypeError[];
WHILE Deg[newDividend] >= Deg[secondArg] DO
coeff: AC.Object ← data.coeffRing.class.divide[LeadingCoeff[newDividend], LeadingCoeff[secondArg] ];
degreeDelta: CARDINAL ← Deg[newDividend] - Deg[secondArg];
multiplier: Polynomial ← Monom[coeff, NEW[CARDINAL ← degreeDelta], firstArg.structure];
product: Polynomial ← Multiply[multiplier, secondArg];
newDividend ← Subtract[newDividend, product];
ENDLOOP;
RETURN[newDividend];
};
Diff: PUBLIC AC.UnaryOp ~ {
Assumes that CARDINAL exponent values can be imbedded in coeffRing, in particular, that coeffRing.class.fromRope[Convert.RopeFromCard[arg.data.first.exponent]] works.
data: PolynomialRingData ← NARROW[arg.structure.instanceData];
inTerms: LIST OF Term ← NARROW[arg.data];
outTerms, outTermsPointer: LIST OF Term ← NIL;
outTerm: Term;
IF ZeroPoly[arg] THEN RETURN[arg];
WHILE inTerms # NIL AND inTerms.first.exponent > 0 DO
imbeddedExp: AC.Object ← data.coeffRing.class.fromRope[Convert.RopeFromCard[inTerms.first.exponent,10, FALSE], data.coeffRing];
newCoeff: AC.Object ← data.coeffRing.class.multiply[imbeddedExp, inTerms.first.coefficient];
outTerm ← NEW[TermRec ← [exponent: inTerms.first.exponent-1, coefficient: newCoeff] ];
IF outTerms # NIL THEN outTermsPointer ← outTermsPointer.rest ← LIST[outTerm]
ELSE outTerms ← outTermsPointer ← LIST[outTerm];
inTerms ← inTerms.rest;
ENDLOOP;
RETURN[ NEW[AC.ObjectRec ← [structure: arg.structure, data: outTerms] ] ];
};
Exp: PROC [val: AC.Object, power: CARDINAL, structure: AC.Structure] RETURNS [result: AC.Object] ~ {
IF power = 0 THEN RETURN[structure.class.one[structure] ];
result ← val;
FOR i:NAT IN [2..power] DO
result ← structure.class.multiply[val, result];
ENDLOOP;
};
MainVarEv: PUBLIC AC.BinaryOp ~ {
data: PolynomialRingData ← NARROW[firstArg.structure.instanceData];
coeffRing: AC.Structure ← data.coeffRing;
inTerms: LIST OF Term ← NARROW[firstArg.data];
IF NOT coeffRing.class.isElementOf[secondArg, coeffRing] THEN TypeError[];
IF ZeroPoly[firstArg] THEN RETURN[coeffRing.class.zero[coeffRing] ];
result ← inTerms.first.coefficient;
WHILE inTerms.rest # NIL DO
degreeDelta: CARDINAL ← DegreeDelta[inTerms];
newAddend: AC.Object ← coeffRing.class.multiply[result, Exp[secondArg, degreeDelta, coeffRing] ];
inTerms ← inTerms.rest;
result ← coeffRing.class.add[newAddend, inTerms.first.coefficient];
ENDLOOP;
result ← coeffRing.class.multiply[result, Exp[secondArg, inTerms.first.exponent, coeffRing] ];
RETURN[ result ];
};
AllVarEv: PUBLIC AC.BinaryOp ~ {
data: PolynomialRingData ← NARROW[firstArg.structure.instanceData];
baseCoeffRing: AC.Structure ← data.baseCoeffRing;
removeMain, mainCood: Points.Point;
pointData: Points.PointData ← NARROW[secondArg.data];
inTerms: LIST OF Term ← NARROW[firstArg.data];
coeff: Polynomial;
at: AC.Object;
IF NOT pointData.dimensionPlus1 = data.allVariables.lengthPlus1 THEN OperationError[$BadPointLength];
IF pointData.dimensionPlus1 - 1 = 1 THEN {
at ← NARROW[secondArg.data, Points.PointData][1];
RETURN[MainVarEv[firstArg, at] ];
};
removeMain ← Points.RemoveMainCood[secondArg];
mainCood ← Points.MainCood[secondArg];
at ← NARROW[mainCood.data, Points.PointData][1];
IF NOT baseCoeffRing.class.isElementOf[at, baseCoeffRing] THEN TypeError[];
IF ZeroPoly[firstArg] THEN RETURN[baseCoeffRing.class.zero[baseCoeffRing] ];
coeff ← NARROW[inTerms.first.coefficient];
result ← AllVarEv[coeff, removeMain]; -- element of baseCoeffRing
WHILE inTerms.rest # NIL DO
degreeDelta: CARDINAL ← DegreeDelta[inTerms];
newAddend: AC.Object ← baseCoeffRing.class.multiply[result, Exp[at, degreeDelta, baseCoeffRing] ];
inTerms ← inTerms.rest;
coeff ← NARROW[inTerms.first.coefficient];
result ← baseCoeffRing.class.add[newAddend, AllVarEv[coeff, removeMain] ];
ENDLOOP;
result ← baseCoeffRing.class.multiply[result, Exp[at, inTerms.first.exponent, baseCoeffRing] ];
RETURN[ result ];
};
Sub: PUBLIC AC.BinaryOp ~ {
polyRing: AC.Structure ← firstArg.structure;
data: PolynomialRingData ← NARROW[polyRing.instanceData];
inTerms: LIST OF Term ← NARROW[firstArg.data];
IF NOT polyRing.class.isElementOf[secondArg, polyRing] THEN TypeError[];
IF ZeroPoly[firstArg] THEN RETURN[firstArg];
result ← Monom[inTerms.first.coefficient, NEW[CARDINAL ← 0], polyRing]; -- lift coefficient
WHILE inTerms.rest # NIL DO
degreeDelta: CARDINAL ← DegreeDelta[inTerms];
trailCoeff: Polynomial;
newAddend: Polynomial ← NARROW[polyRing.class.multiply[result, Exp[secondArg, degreeDelta, polyRing] ] ];
inTerms ← inTerms.rest;
trailCoeff ← Monom[inTerms.first.coefficient, NEW[CARDINAL ← 0], polyRing]; -- lift
result ← NARROW[ polyRing.class.add[newAddend, trailCoeff] ];
ENDLOOP;
result ← NARROW[ polyRing.class.multiply[result, Exp[secondArg, inTerms.first.exponent, polyRing] ] ];
RETURN[ result ];
};
DegreeDelta: PUBLIC PROC [terms: LIST OF Term] RETURNS [CARDINAL] ~ {
IF terms = NIL THEN RETURN[0];
IF terms.rest = NIL THEN RETURN[terms.first.exponent];
RETURN[terms.first.exponent - terms.rest.first.exponent];
};
SylvMatrix: PUBLIC AC.BinaryOp ~ {
data: PolynomialRingData ← NARROW[firstArg.structure.instanceData];
coeffRing: AC.Structure ← data.coeffRing;
degree1: CARDINAL ← Deg[firstArg];
degree2: CARDINAL ← Deg[secondArg];
size: CARDINAL ← degree1 + degree2;
matrixStructure: AC.Structure;
rows: Matrices.RowSeq ← NEW[Matrices.RowSeqRec[size]];
zero: AC.Object ← coeffRing.class.zero[coeffRing];
IF degree1 = 0 AND degree2 = 0 THEN RETURN[NIL]; -- undefined in this case
matrixStructure ← Matrices.MakeMatrixStructure[coeffRing, size];
IF degree1 = 0 THEN {
IF ZeroPoly[firstArg] THEN
RETURN[NARROW[MAT.DiagonalMatrix[matrixStructure][zero, matrixStructure] ] ]
ELSE {
polyData: PolynomialData ← NARROW[firstArg.data];
RETURN[NARROW[MAT.DiagonalMatrix[matrixStructure][polyData.first.coefficient, matrixStructure] ] ];
};
};
IF degree2 = 0 THEN {
IF ZeroPoly[secondArg] THEN
RETURN[NARROW[MAT.DiagonalMatrix[matrixStructure][zero, matrixStructure] ] ]
ELSE {
polyData: PolynomialData ← NARROW[secondArg.data];
RETURN[NARROW[MAT.DiagonalMatrix[matrixStructure][polyData.first.coefficient, matrixStructure] ] ];
};
};
FOR i:NAT IN [1..degree2] DO
row: Matrices.Row ← NEW[Matrices.RowRec[size]];
in1Terms: LIST OF Term ← NARROW[firstArg.data];
FOR j:NAT IN [1..i-1] DO
row[j] ← zero;
ENDLOOP;
row[i] ← in1Terms.first.coefficient;
WHILE in1Terms.rest # NIL DO
degreeDelta: CARDINAL ← DegreeDelta[in1Terms];
FOR j: NAT IN [1..degreeDelta-1] DO
row[i+degree1-in1Terms.first.exponent+j] ← zero;
ENDLOOP;
in1Terms ← in1Terms.rest;
row[i+degree1-in1Terms.first.exponent] ← in1Terms.first.coefficient;
ENDLOOP;
FOR j: NAT IN [1..in1Terms.first.exponent+(degree2-i)] DO
row[size+1-j] ← zero;
ENDLOOP;
rows[i] ← row;
ENDLOOP;
FOR i:NAT IN [1..degree1] DO
row: Matrices.Row ← NEW[Matrices.RowRec[size]];
in2Terms: LIST OF Term ← NARROW[secondArg.data];
FOR j:NAT IN [1..i-1] DO
row[j] ← zero;
ENDLOOP;
row[i] ← in2Terms.first.coefficient;
WHILE in2Terms.rest # NIL DO
degreeDelta: CARDINAL ← DegreeDelta[in2Terms];
FOR j: NAT IN [1..degreeDelta-1] DO
row[i+degree2-in2Terms.first.exponent+j] ← zero;
ENDLOOP;
in2Terms ← in2Terms.rest;
row[i+degree2-in2Terms.first.exponent] ← in2Terms.first.coefficient;
ENDLOOP;
FOR j: NAT IN [1..in2Terms.first.exponent+(degree1-i)] DO
row[size+1-j] ← zero;
ENDLOOP;
rows[degree2+i] ← row;
ENDLOOP;
result ← NEW[AC.ObjectRec ← [structure: matrixStructure, data: rows] ];
};
Res: PUBLIC AC.BinaryOp ~ {
sylvesterMatrix: Matrices.Matrix ← SylvMatrix[firstArg, secondArg];
RETURN[Matrices.Determinant[sylvesterMatrix.structure][sylvesterMatrix] ];
};
Comparison
ZeroPoly: PUBLIC PROC [in: Polynomial] RETURNS [BOOL] ~ {
RETURN[Rope.Equal[in.structure.class.flavor, "Polynomials"] AND in.data = NIL] };
Equal: PUBLIC AC.EqualityOp ~ {
RETURN[ZeroPoly[Subtract[firstArg, secondArg] ] ];
};
Sign: PUBLIC AC.CompareToZeroOp ~ {
data: PolynomialRingData ← NARROW[arg.structure.instanceData];
polyData: PolynomialData;
IF NOT arg.structure.class.ordered THEN TypeError[$UnorderedStructure];
IF ZeroPoly[arg] THEN RETURN[equal]; -- needed since arg.data = NIL for zero poly
polyData ← NARROW[arg.data];
RETURN[data.coeffRing.class.sign[polyData.first.coefficient] ];
};
Abs: PUBLIC AC.UnaryOp ~ {
IF NOT arg.structure.class.ordered THEN TypeError[$UnorderedStructure];
IF Sign[arg] = less THEN RETURN[Negate[arg]] ELSE RETURN[arg];
};
Compare: PUBLIC AC.BinaryCompareOp ~ {
IF NOT firstArg.structure.class.ordered THEN TypeError[$UnorderedStructure];
RETURN[Sign[ Subtract[firstArg, secondArg] ]];
};
Selection
LeadingCoeff: PUBLIC AC.UnaryOp ~ {
data: PolynomialRingData ← NARROW[arg.structure.instanceData];
polyData: PolynomialData;
IF ZeroPoly[arg] THEN RETURN[data.coeffRing.class.zero[data.coeffRing] ];
polyData ← NARROW[arg.data];
RETURN[polyData.first.coefficient];
};
Deg: PUBLIC AC.ElementRankOp ~ {
polyData: PolynomialData;
IF ZeroPoly[arg] THEN RETURN[0];
polyData ← NARROW[arg.data];
RETURN[polyData.first.exponent];
};
Reductum: PUBLIC AC.UnaryOp ~ {
polyData: PolynomialData;
IF ZeroPoly[arg] THEN RETURN[arg];
polyData ← NARROW[arg.data];
RETURN[ NEW[AC.ObjectRec ← [structure: arg.structure, data: polyData.rest] ] ];
};
END.