DIRECTORY
Convert USING [AtomFromRope, Error, IntFromRope, RealFromRope, RopeFromAtom, RopeFromInt, RopeFromReal],
IO USING [BreakProc, CharClass, EndOfStream, GetTokenRope, NUL, RIS, SP, STREAM],
MultiPolynomial,
Polynomial USING [Constant, Degree, Linear, Product, Ref, Sum],
Rope USING [Compare, Concat, Fetch, Length, ROPE];
Implementation notes
Term invariant: An invariant exists on all Monomial records. The variables in the PoweredVariableSequence of the Monomial must (1) be distinct, (2) be sorted in ascending alphabetical order, as determined by Rope.Compare, and (3) have positive degrees. Currently, there is only case where a Monomial not satisfying this invariant exists. This is when a monomial expression is being constructed by RefFromRope. After intermediate Monomials are created, they are converted into a form that satisfies the invariant by FixUpMonomial. Several procedures depend on the term invariant being satisfied.
There is one "unclean" procedure, PolynomialCoefficient. If/when Polynomial gets coefficient selectors, this should be changed. Also, if Polynomial gets a RefFromSequence constructor, UnivariateFromRef can be made a bit smoother. Additionally, it would be nice to have a better IntegerPower procedure.
Type declarations
Ref: TYPE ~ MultiPolynomial.Ref;
MultiPolRec: TYPE ~ MultiPolynomial.MultiPolRec;
Monomial: TYPE ~ MultiPolynomial.Monomial;
PoweredVariable: TYPE ~ MultiPolynomial.PoweredVariable;
PoweredVariableSequence: TYPE ~ MultiPolynomial.PoweredVariableSequence;
EvaluationBinding: TYPE ~ MultiPolynomial.EvaluationBinding;
EvaluationBindings: TYPE ~ MultiPolynomial.EvaluationBindings;
SubstitutionBinding: TYPE ~ MultiPolynomial.SubstitutionBinding;
SubstitutionBindings: TYPE ~ MultiPolynomial.SubstitutionBindings;
TSUBinding: TYPE ~ MultiPolynomial.TSUBinding;
TSUBindings: TYPE ~ MultiPolynomial.TSUBindings;
Error: PUBLIC ERROR[reason: ErrorCode] = CODE;
ErrorCode: TYPE ~ MultiPolynomial.ErrorCode;
Substitution Mechanisms
Evaluate:
PUBLIC
PROC [a: Ref, variable:
ATOM, value:
REAL]
RETURNS [result: Ref] ~
BEGIN
break: SIGNAL = CODE;
intermediate: Ref ← NEW[MultiPolRec[a.length]];
Iterate over each monomial term
FOR i:
NAT
IN [0..a.length)
DO
thisMonomial: Monomial ← a[i];
newCoefficient: REAL;
newVars: REF PoweredVariableSequence;
See if the monomial contains the variable.
variablePresent: BOOLEAN;
variablePosition: NAT;
[variablePresent, variablePosition] ← VariablePositionInMonomial[thisMonomial, variable];
If it does, compute the evaluated monomial.
IF variablePresent
THEN
BEGIN
degree: NAT ← thisMonomial.vars[variablePosition].degree;
newCoefficient ← thisMonomial.coefficient * IntegerPower[value, degree];
newVars ← NEW[PoweredVariableSequence[thisMonomial.vars.numVars - 1]];
FOR j:
NAT
IN [0..variablePosition)
DO
newVars[j] ← thisMonomial.vars[j];
ENDLOOP;
FOR j:
NAT
IN (variablePosition..thisMonomial.vars.numVars)
DO
newVars[j-1] ← thisMonomial.vars[j];
ENDLOOP;
END
If it doesn't, the new monomial is just a copy of the old one.
ELSE
BEGIN
newCoefficient ← thisMonomial.coefficient;
newVars ← NEW[PoweredVariableSequence[thisMonomial.vars.numVars]];
FOR j:
NAT
IN [0..thisMonomial.vars.numVars)
DO
newVars[j] ← thisMonomial.vars[j];
ENDLOOP;
END;
Store the new monomial in the intermediate polynomial.
intermediate[i] ← [coefficient: newCoefficient, vars: newVars];
ENDLOOP;
Compact the intermediate polynomial and return it.
result ← Compact[intermediate];
EvaluateList:
PUBLIC
PROC [a: Ref, bindings: EvaluationBindings]
RETURNS [result: Ref] ~
BEGIN
bindingList: EvaluationBindings ← bindings;
intermediate: Ref ← a;
UNTIL bindingList =
NIL
DO
intermediate ← Evaluate[intermediate, bindingList.first.variable, bindingList.first.value];
bindingList ← bindingList.rest;
ENDLOOP;
result ← intermediate;
END;
TotallyEvaluate:
PUBLIC
PROC [a: Ref, bindings: EvaluationBindings]
RETURNS [result:
REAL] ~
BEGIN
intermediate: Ref ← EvaluateList[a, bindings];
value: REAL;
value ← RealFromRef[intermediate
! Error[NotConstant] => ERROR Error[UnsubstitutedVariables]];
result ← value;
END;
Substitute:
PUBLIC
PROC [a: Ref, variable:
ATOM, replacement: Ref]
RETURNS [result: Ref] ~
BEGIN
intermediate: Ref ← NEW[MultiPolRec[0]];
Iterate over each monomial term
FOR i:
NAT
IN [0..a.length)
DO
thisMonomial: Monomial ← a[i];
newMonomial: Monomial;
newPolynomial: Ref;
See if the monomial contains the variable.
variablePresent: BOOL ← FALSE;
variablePosition: NAT;
FOR j:
NAT
IN [0..thisMonomial.vars.numVars)
DO
IF thisMonomial.vars[j].variable = variable
THEN {
variablePresent ← TRUE;
variablePosition ← j}
ENDLOOP;
If it does, compute newMonomial, the monomial without the variable, and newPolynomial, the product of newMonomial and the proper power of replacement.
IF variablePresent
THEN
BEGIN
degree: NAT ← thisMonomial.vars[variablePosition].degree;
thisMonomialConvertedToPoly: Ref ← NEW[MultiPolRec[1]];
thisMonomialConvertedToPoly[0] ← thisMonomial;
newMonomial ← Evaluate[thisMonomialConvertedToPoly, variable, 1][0];
newPolynomial ← MultiplyMonomialByRef[newMonomial, Power[replacement, degree]];
END
If it doesn't, the new polynomial is just a copy of the old monomial.
ELSE
BEGIN
newPolynomial ← NEW[MultiPolRec[1]];
newPolynomial[0] ← thisMonomial;
END;
Add the new polynomial to the intermediate polynomial.
intermediate ← AddRefsWithoutCompaction[intermediate, newPolynomial];
ENDLOOP;
Compact the intermediate polynomial and return it.
result ← Compact[intermediate];
END;
SubstituteList:
PUBLIC
PROC [a: Ref, bindings: SubstitutionBindings]
RETURNS [result: Ref] ~
BEGIN
bindingList: SubstitutionBindings ← bindings;
intermediate: Ref ← a;
UNTIL bindingList =
NIL
DO
intermediate ← Substitute[
intermediate,
bindingList.first.variable,
bindingList.first.replacement];
bindingList ← bindingList.rest;
ENDLOOP;
result ← intermediate;
END;
TotallySubstituteUnivariates:
PUBLIC
PROC [a: Ref, bindings: TSUBindings]
RETURNS [result: Polynomial.Ref] ~
BEGIN
bindingList: TSUBindings ← bindings;
intermediate: Ref ← a;
UNTIL bindingList =
NIL
DO
intermediate ← Substitute[
intermediate,
bindingList.first.variable,
RefFromUnivariate[bindingList.first.replacement, $t]
];
bindingList ← bindingList.rest;
ENDLOOP;
result ← UnivariateFromRef[intermediate, $t
! Error[NotUnivariate] => ERROR Error[UnsubstitutedVariables]];
END;
Internal Procedures
AddRefsWithoutCompaction:
PROC [a, b: Ref]
RETURNS [result: Ref] ~
BEGIN
result ← NEW[MultiPolRec[a.length + b.length]];
FOR i:
NAT
IN [0..a.length)
DO
result[i] ← a[i];
ENDLOOP;
FOR i:
NAT
IN [0..b.length)
DO
result[i+a.length] ← b[i];
ENDLOOP;
END;
Compact:
PROC [a: Ref]
RETURNS [result: Ref] ~
BEGIN
This procedure compacts a polynomial by combining all equivalent terms. It makes an intermediate polynomial, the same size as the original, with equivalent terms combined. Another pass removes all the zero terms. Then the intermediate polynomial is copied to a possibly smaller and more compact one.
intermediate: Ref ← NEW[MultiPolRec[a.length]];
monomialsCopied: NAT ← 0;
finalIndex: NAT;
Iterate over the monomials in a.
FOR i:
NAT
IN [0..a.length)
DO
thisMonomial: Monomial ← a[i];
copied: BOOLEAN ← FALSE;
If an equivalent monomial has already been copied, add thisMonomial to it.
FOR j:
NAT
IN [0..monomialsCopied)
DO
IF Equivalent[thisMonomial.vars, intermediate[j].vars]
THEN
BEGIN
intermediate[j].coefficient ← intermediate[j].coefficient + thisMonomial.coefficient;
copied ← TRUE;
END;
ENDLOOP;
Otherwise, copy thisMonomial to the intermediate polynomial.
IF copied =
FALSE
AND thisMonomial.coefficient # 0
THEN
BEGIN
intermediate[monomialsCopied] ← thisMonomial;
monomialsCopied ← monomialsCopied + 1;
END;
ENDLOOP;
Squeeze out all the zero terms.
finalIndex ← 0;
FOR i:
NAT
IN [0..monomialsCopied)
DO
IF intermediate[i].coefficient # 0
THEN
BEGIN
intermediate[finalIndex] ← intermediate[i];
finalIndex ← finalIndex + 1;
END;
ENDLOOP;
Make a more compact polynomial from the intermediate.
result ← NEW[MultiPolRec[finalIndex]];
FOR i:
NAT
IN [0..finalIndex)
DO
result[i] ← intermediate[i]
ENDLOOP;
CopyMonomial:
PROC [m: Monomial]
RETURNS [result: Monomial] ~
BEGIN
newVars: REF PoweredVariableSequence ← NEW[PoweredVariableSequence[m.vars.numVars]];
FOR i:
NAT
IN [0..m.vars.numVars)
DO
newVars[i] ← [variable: m.vars[i].variable, degree: m.vars[i].degree];
ENDLOOP;
result ← [coefficient: m.coefficient, vars: newVars];
END;
DifferentiateMonomial:
PROC [m: Monomial, variable:
ATOM]
RETURNS [result: Monomial] ~
BEGIN
newCoefficient: REAL;
newVars: REF PoweredVariableSequence;
See if the variable appears in the monomial.
variablePresent: BOOL ← FALSE;
variablePosition: NAT;
FOR j:
NAT
IN [0..m.vars.numVars)
DO
IF m.vars[j].variable = variable
THEN
BEGIN
variablePresent ← TRUE;
variablePosition ← j;
END;
ENDLOOP;
If it does, compute the derivative.
IF variablePresent
THEN
BEGIN
degree: NAT ← m.vars[variablePosition].degree;
newCoefficient ← m.coefficient * degree;
newVars ← NEW[PoweredVariableSequence[m.vars.numVars]];
FOR j:
NAT
IN [0..m.vars.numVars)
DO
newVars[j] ← m.vars[j];
ENDLOOP;
newVars[variablePosition].degree ← degree - 1; -- if this is zero, it should be eliminated.
END
If it doesn't, the derivative is zero.
ELSE
BEGIN
newCoefficient ← 0;
newVars ← NEW[PoweredVariableSequence[0]];
END;
Create and return the new monomial.
result ← [coefficient: newCoefficient, vars: newVars];
END;
Equivalent:
PROC [p, q:
REF PoweredVariableSequence]
RETURNS [result:
BOOLEAN] ~
BEGIN
Returns true if p and q are equivalent (same powers of same variables). Note that the implementation depends on the term invariant (variables are sorted).
IF p.numVars # q.numVars THEN {result ← FALSE; RETURN};
FOR i:
NAT
IN [0..p.numVars)
DO
IF (p[i].variable # q[i].variable)
OR (p[i].degree # q[i].degree)
THEN
BEGIN
result ← FALSE;
RETURN;
END;
ENDLOOP;
result ← TRUE;
END;
FixUpMonomial:
PROC [m: Monomial]
RETURNS [result: Monomial] ~
BEGIN
The input, m, is a monomial that does not necessarily satisfy the term invariant. Its variables might be unsorted, or it may have several occurrences of the same variable in its vars. The result is a mathematically equivalent monomial that satisfies the term invariant.
First, the powered variables in m are sorted, adding the exponents of identical variables, into an intermediate PoweredVariableSequence. Then this sequence is used to create the result monomial, whose PoweredVariableSequence might be shorter.
Declarations
intermediate: REF PoweredVariableSequence ← NEW[PoweredVariableSequence[m.vars.numVars]];
lastVariableCopied: Rope.ROPE ← "";
lastVariableToCopy: Rope.ROPE ← "";
iIndex: NAT ← 0;
Figure out the last variable, alphabetically.
FOR i:
NAT
IN [0..m.vars.numVars)
DO
thisVar: Rope.ROPE ← Convert.RopeFromAtom[m.vars[i].variable, FALSE];
IF Rope.Compare[thisVar, lastVariableToCopy] = greater THEN lastVariableToCopy ← thisVar;
ENDLOOP;
Copy each variable into the intermediate sequence.
UNTIL Rope.Compare[lastVariableCopied, lastVariableToCopy] = equal
DO
Declarations
varToCopy: ATOM;
variableToCopy: Rope.ROPE ← lastVariableToCopy;
exponent: NAT ← 0;
Figure out the next variable to copy.
FOR i:
NAT
IN [0..m.vars.numVars)
DO
thisVar: Rope.ROPE ← Convert.RopeFromAtom[m.vars[i].variable, FALSE];
IF Rope.Compare[lastVariableCopied, thisVar] = less
AND Rope.Compare[thisVar, variableToCopy] = less
THEN
BEGIN
variableToCopy ← thisVar;
END;
ENDLOOP;
varToCopy ← Convert.AtomFromRope[variableToCopy];
Accumulate that variable's exponents.
FOR i:
NAT
IN [0..m.vars.numVars)
DO
thisVar: Rope.ROPE ← Convert.RopeFromAtom[m.vars[i].variable, FALSE];
IF Rope.Compare[thisVar, variableToCopy] = equal
THEN
exponent ← exponent + m.vars[i].degree;
ENDLOOP;
Put the powered variable into the intermediate sequence, if its exponent is positive.
SELECT exponent
FROM
<0 => ERROR Error[InternalError];
=0 => NULL;
>0 => {
intermediate[iIndex] ← [variable: varToCopy, degree: exponent];
iIndex ← iIndex + 1};
ENDCASE;
Save the variable as the last variable copied, and do the next variable.
lastVariableCopied ← variableToCopy;
ENDLOOP;
Construct the final monomial.
result ← [
coefficient: m.coefficient,
vars: NEW[PoweredVariableSequence[iIndex]]];
FOR i:
NAT
IN [0..iIndex)
DO
result.vars[i] ← intermediate[i];
ENDLOOP;
LegalVariableName:
PROC [name:
ATOM]
RETURNS [legal:
BOOLEAN] ~
BEGIN
A variable name is legal if it is a single alphabetic character.
RETURN[Rope.Length[Convert.RopeFromAtom[name, FALSE]] = 1];
END;
MonomialDegree:
PROC [m: Monomial]
RETURNS [degree:
NAT] ~
BEGIN
d: NAT ← 0;
FOR i:
NAT
IN [0..m.vars.numVars)
DO
d ← d + m.vars[i].degree;
ENDLOOP;
degree ← d;
END;
MultiplyMonomials:
PROC [m, n: Monomial]
RETURNS [result: Monomial] ~
BEGIN
Multiplies two monomials, preserving the term invariant.
Selection: TYPE ~ {fromM, fromN, fromBoth};
intermediate: REF PoweredVariableSequence ← NEW[PoweredVariableSequence[m.vars.numVars + n.vars.numVars]];
resultTerm: REF PoweredVariableSequence;
iIndex: NAT ← 0;
mIndex: NAT ← 0;
nIndex: NAT ← 0;
IF m.coefficient = 0
OR n.coefficient = 0
THEN
BEGIN
resultTerm ← NEW[PoweredVariableSequence[0]];
result ← [coefficient: 0, vars: resultTerm];
RETURN;
END;
This is effectively an iteration over iIndex. Each pass through the loop adds a PoweredVariable to the intermediate term. At the end of each pass, iIndex equals the number of variables so far copied to the intermediate term.
UNTIL (mIndex >= m.vars.numVars)
AND (nIndex >= n.vars.numVars)
DO
selection: Selection;
Determine if the next variable will be taken from m, n, or both.
IF mIndex >= m.vars.numVars THEN selection ← fromN
ELSE IF nIndex >= n.vars.numVars THEN selection ← fromM
ELSE
BEGIN
nextMvariable: Rope.ROPE ← Convert.RopeFromAtom[m.vars[mIndex].variable];
nextNvariable: Rope.ROPE ← Convert.RopeFromAtom[n.vars[nIndex].variable];
SELECT Rope.Compare[nextMvariable, nextNvariable]
FROM
less => selection ← fromM;
equal => selection ← fromBoth;
greater => selection ← fromN;
ENDCASE;
END;
SELECT selection
FROM
fromM =>
BEGIN
intermediate[iIndex] ← m.vars[mIndex];
mIndex ← mIndex + 1;
END;
fromN =>
BEGIN
intermediate[iIndex] ← n.vars[nIndex];
nIndex ← nIndex + 1;
END;
fromBoth =>
BEGIN
intermediate[iIndex].variable ← m.vars[mIndex].variable;
intermediate[iIndex].degree ← m.vars[mIndex].degree + n.vars[nIndex].degree;
mIndex ← mIndex + 1;
nIndex ← nIndex + 1;
END;
ENDCASE;
iIndex ← iIndex + 1;
ENDLOOP;
Compress and construct the resulting monomial.
resultTerm ← NEW[PoweredVariableSequence[iIndex]];
FOR i:
NAT
IN [0..iIndex)
DO
resultTerm[i] ← intermediate[i];
ENDLOOP;
result ← [coefficient: m.coefficient * n.coefficient, vars: resultTerm];
END;
MultiplyMonomialByRef:
PROC [m: Monomial, a: Ref]
RETURNS [result: Ref] ~
BEGIN
If either parameter is zero, return a zero polynomial.
IF m.coefficient = 0
OR a.length = 0
THEN
BEGIN
result ← NEW[MultiPolRec[0]];
RETURN;
END;
Otherwise, multiply each monomial in a by m and store them in the result polynomial. Note that no compaction is needed.
result ← NEW[MultiPolRec[a.length]];
FOR i:
NAT
IN [0..a.length)
DO
result[i] ← MultiplyMonomials[m, a[i]];
ENDLOOP;
END;
MultiplyRefsWithoutCompaction:
PROC [a, b: Ref]
RETURNS [result: Ref] ~
BEGIN
Uses the distributive law to multiply two Refs. Each monomial in a is multiplied by b; the results are accumulated in the intermediate polynomial.
intermediate: Ref ← NEW[MultiPolRec[0]];
FOR i:
NAT
IN [0..a.length)
DO
distributedProduct: Ref ← MultiplyMonomialByRef[a[i], b];
intermediate ← AddRefsWithoutCompaction[intermediate, distributedProduct];
ENDLOOP;
result ← intermediate;
END;
PolynomialCoefficient:
PROC [poly: Polynomial.Ref, degree:
NAT]
RETURNS [coeff:
REAL] ~
BEGIN
THIS BELONGS IN Polynomial MODULE
coeff ← poly[degree];
END;
VariablePositionInMonomial:
PROC [m: Monomial, var:
ATOM]
RETURNS [variablePresent:
BOOLEAN, variablePosition:
NAT] ~
BEGIN
j: NAT ← 0;
variablePresent ← FALSE;
UNTIL variablePresent
OR j = m.vars.numVars
DO
IF m.vars[j].variable = var
THEN {
variablePresent ← TRUE;
variablePosition ← j};
j ← j + 1;
ENDLOOP;
END;
IntegerPower:
PROC [base:
REAL, exponent:
INTEGER]
RETURNS [result:
REAL] ~
BEGIN
SELECT exponent
FROM
0 => result ← 1;
1 => result ← base;
2 => result ← base*base;
3 => result ← base*base*base;
4 => {foo: REAL ← base*base; result ← foo*foo};
5 => {foo: REAL ← base*base; result ← foo*foo*base};
6 => {foo: REAL ← base*base; result ← foo*foo*foo};
>6 => {foo: REAL ← base*base; result ← foo*foo*foo*IntegerPower[base, exponent-6]};
ENDCASE;
END;