MultiPolynomialImpl.mesa
Copyright Ó 1986, 1989, 1992 by Xerox Corporation. All rights reserved.
James Rauen August 14, 1986 6:13:01 pm PDT
Doug Wyatt, September 13, 1989 11:48:38 am PDT
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];
MultiPolynomialImpl: CEDAR PROGRAM
IMPORTS Convert, IO, Polynomial, Rope
EXPORTS MultiPolynomial
~ BEGIN
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;
Conversions to and from other representations
UnivariateFromRef: PUBLIC PROC [a: Ref, var: ATOM] RETURNS [result: Polynomial.Ref] ~ BEGIN
Things needed to assemble the final result.
x: Polynomial.Ref ¬ Polynomial.Linear[[0, 1]];
intermediate: Polynomial.Ref;
degree: NAT ¬ Degree[a];
Make a sequence for the coefficients.
CoeffSeq: TYPE ~ RECORD[coeffs: SEQUENCE length: NAT OF REAL];
coefficients: REF CoeffSeq ¬ NEW[CoeffSeq[degree + 1]];
FOR i: NAT IN [0..coefficients.length) DO coefficients[i] ¬ 0; ENDLOOP;
Add each term of a into the coefficient sequence.
FOR i: NAT IN [0..a.length) DO
SELECT a[i].vars.numVars FROM
=0 =>
coefficients[0] ¬ coefficients[0] + a[i].coefficient;
=1 =>
IF a[i].vars[0].variable= var THEN
coefficients[a[i].vars[0].degree] ¬ coefficients[a[i].vars[0].degree] + a[i].coefficient
ELSE ERROR Error[NotUnivariate];
>1 =>
ERROR Error[NotUnivariate];
ENDCASE;
ENDLOOP;
Convert the coefficient sequence into a Polynomial.Ref.
intermediate ¬ Polynomial.Constant[[0]];
FOR i: NAT IN [0..degree] DO
newConstantTerm: Polynomial.Ref ¬ Polynomial.Constant[[coefficients[degree - i]]];
intermediate ¬ Polynomial.Product[intermediate, x];
intermediate ¬ Polynomial.Sum[intermediate, newConstantTerm];
ENDLOOP;
result ¬ intermediate;
END;
RefFromUnivariate: PUBLIC PROC [pol: Polynomial.Ref, var: ATOM] RETURNS [result: Ref] ~ BEGIN
degree: NAT ¬ Polynomial.Degree[pol];
intermediate: Ref ¬ NEW[MultiPolRec[degree + 1]];
IF ~LegalVariableName[var] THEN ERROR Error[BadVariableName];
FOR i: NAT IN [0..degree] DO
vars: REF PoweredVariableSequence ¬ NEW[PoweredVariableSequence[1]];
vars[0] ¬ [variable: var, degree: i];
intermediate[i] ¬ [coefficient: PolynomialCoefficient[pol, i], vars: vars];
ENDLOOP;
result ¬ Compact[intermediate];
END;
RealFromRef: PUBLIC PROC [a: Ref] RETURNS [constantTerm: REAL] ~ BEGIN
SELECT a.length FROM
0 => RETURN[0];
1 => IF a[0].vars.numVars = 0 THEN RETURN[a[0].coefficient];
ENDCASE;
ERROR Error[NotConstant];
END;
RopeFromRef: PUBLIC PROC [a: Ref] RETURNS [r: Rope.ROPE] ~ BEGIN
result: Rope.ROPE ¬ "";
FOR i: NAT IN [0..a.length) DO
thisMonomial: Monomial ¬ a[i];
Do the sign. Treat the first term specially.
IF i = 0 THEN
{IF a[i].coefficient < 0 THEN result ¬ Rope.Concat[result, "-"]}
ELSE
{IF a[i].coefficient < 0 THEN result ¬ Rope.Concat[result, " - "]
ELSE result ¬ Rope.Concat[result, " + "]};
result ¬ Rope.Concat[result, Convert.RopeFromReal[ABS[thisMonomial.coefficient]]];
FOR j: NAT IN [0..thisMonomial.vars.numVars) DO
thisVar: PoweredVariable ¬ thisMonomial.vars[j];
SELECT thisVar.degree FROM
=0 => {};
=1 => result ¬ Rope.Concat[result, Convert.RopeFromAtom[thisVar.variable, FALSE]];
>1 => {result ¬ Rope.Concat[result, Convert.RopeFromAtom[thisVar.variable, FALSE]];
result ¬ Rope.Concat[result, "^"];
result ¬ Rope.Concat[result, Convert.RopeFromInt[thisVar.degree]];};
ENDCASE;
ENDLOOP;
ENDLOOP;
IF a.length = 0 THEN result ¬ "0.0";
r ¬ result;
END;
RefFromRope: PUBLIC PROC [rope: Rope.ROPE] RETURNS [result: Ref] ~ BEGIN
Declarations.
intermRef: Ref;
Convert the rope into a stream.
inStream: IO.STREAM ¬ IO.RIS[rope];
Define the break procedure used to separate the stream into tokens.
BreakProc: IO.BreakProc ~ BEGIN
charClass: IO.CharClass;
SELECT char FROM
IN [IO.NUL .. IO.SP]  => charClass ¬ sepr;
'+, '-, '^  => charClass ¬ break;
IN ['a..'z], IN ['A..'Z] => charClass ¬ break;
IN ['0..'9], '.  => charClass ¬ other;
ENDCASE  => ERROR Error[BadFormat];
RETURN [charClass];
END;
Define a procedure to get the next token from the stream.
NextToken: PROC [stream: IO.STREAM ¬ inStream] RETURNS [Rope.ROPE] ~ BEGIN
token: Rope.ROPE;
[token: token] ¬ IO.GetTokenRope[stream, BreakProc];
RETURN[token];
END;
These are the intermediate multipolynomial construction procedures.
TemporaryMultiPolRec: TYPE ~ RECORD [
numberOfMonomials: NAT,
monomials: LIST OF TemporaryMonomial];
TemporaryMonomial: TYPE ~ RECORD [
sign: Rope.ROPE,
coefficient: REAL,
numberOfVariables: NAT,
poweredVariables: LIST OF PoweredVariable];
temporaryMultiPol: TemporaryMultiPolRec;
currentMonomial: TemporaryMonomial;
FlushTemporaryMultiPol: PROC[] ~ BEGIN
temporaryMultiPol ¬ [0, NIL];
END;
ClearCurrentMonomial: PROC[] ~ BEGIN
currentMonomial ¬ ["+", 1.0, 0, NIL];
END;
StoreCurrentMonomial: PROC[] ~ BEGIN
temporaryMultiPol.monomials ¬ CONS[currentMonomial, temporaryMultiPol.monomials];
END;
AddNewTerm: PROC[sign: Rope.ROPE] ~ BEGIN
IF temporaryMultiPol.numberOfMonomials > 0 THEN StoreCurrentMonomial[];
temporaryMultiPol.numberOfMonomials ¬ temporaryMultiPol.numberOfMonomials + 1;
ClearCurrentMonomial[];
currentMonomial.sign ¬ sign;
END;
ChangeCoefficient: PROC[coeff: Rope.ROPE] ~ BEGIN
currentMonomial.coefficient ¬ Convert.RealFromRope[coeff !
Convert.Error => {ERROR Error[BadFormat]}];
END;
AddNewVariable: PROC[var: Rope.ROPE] ~ BEGIN
currentMonomial.numberOfVariables ¬ currentMonomial.numberOfVariables + 1;
currentMonomial.poweredVariables ¬ CONS [
[variable: Convert.AtomFromRope[var], degree: 1],
currentMonomial.poweredVariables];
END;
ChangePower: PROC[power: Rope.ROPE] ~ BEGIN
currentMonomial.poweredVariables.first.degree ¬ Convert.IntFromRope[power !
Convert.Error => {ERROR Error[BadFormat]}];
END;
These convert the intermediate representations to actual representations.
MonomialFromTemporary: PROC[temp: TemporaryMonomial] RETURNS[Monomial] ~ BEGIN
vars: REF PoweredVariableSequence ¬ NEW[PoweredVariableSequence[temp.numberOfVariables]];
coefficient: REAL ¬ temp.coefficient;
poweredVariableList: LIST OF PoweredVariable ¬ temp.poweredVariables;
index: NAT ¬ temp.numberOfVariables;
UNTIL poweredVariableList = NIL DO
index ¬ index - 1;
vars[index] ¬ poweredVariableList.first;
poweredVariableList ¬ poweredVariableList.rest;
ENDLOOP;
IF Rope.Compare[temp.sign, "-"] = equal THEN coefficient ¬ - coefficient;
RETURN[FixUpMonomial[[coefficient, vars]]];
END;
RefFromTemporary: PROC RETURNS[Ref] ~ BEGIN
temp: TemporaryMultiPolRec ¬ temporaryMultiPol;
result: Ref ¬ NEW[MultiPolRec[temp.numberOfMonomials]];
monomialList: LIST OF TemporaryMonomial ¬ temp.monomials;
index: NAT ¬ temp.numberOfMonomials;
UNTIL monomialList = NIL DO
index ¬ index - 1;
result[index] ¬ MonomialFromTemporary[monomialList.first];
monomialList ¬ monomialList.rest;
ENDLOOP;
RETURN[result];
END;
These are the kinds of tokens that might appear.
TokenType: TYPE ~ {sign, coefficient, variable, caret, exponent, begin, end, error};
Parse the stream, one token at a time.
previousTokenType: TokenType ¬ begin;
FlushTemporaryMultiPol[];
ClearCurrentMonomial[];
UNTIL previousTokenType = end DO
Declarations for the next token.
nextToken: Rope.ROPE;
firstCharOfNextToken: CHAR;
nextTokenType: TokenType;
Get the next token and decide its type. If the end of the stream is encountered, the token "&EOS" is used, with TokenType end.
nextToken ¬ NextToken[! IO.EndOfStream => {nextToken ¬ "&EOS"; CONTINUE}];
firstCharOfNextToken ¬ Rope.Fetch[nextToken, 0];
nextTokenType ¬ SELECT firstCharOfNextToken FROM
IN ['a..'z], IN ['A..'Z] => variable,
IN ['0..'9], '.  => IF previousTokenType = caret THEN exponent ELSE coefficient,
'+, '-  => sign,
'^  => caret,
'&  => end,
ENDCASE  => error;
Take action, depending on the type of next token and previous token.
SELECT previousTokenType FROM
begin => SELECT nextTokenType FROM
sign => AddNewTerm[sign: nextToken];
coefficient => {AddNewTerm[sign: "+"]; ChangeCoefficient[nextToken]};
variable => {AddNewTerm[sign: "+"]; AddNewVariable[nextToken]};
ENDCASE => ERROR Error[BadFormat];
sign => SELECT nextTokenType FROM
coefficient => ChangeCoefficient[nextToken];
variable => AddNewVariable[nextToken];
ENDCASE => ERROR Error[BadFormat];
coefficient => SELECT nextTokenType FROM
sign => AddNewTerm[sign: nextToken];
variable => AddNewVariable[nextToken];
end => {};
ENDCASE => ERROR Error[BadFormat];
variable => SELECT nextTokenType FROM
sign => AddNewTerm[sign: nextToken];
variable => AddNewVariable[nextToken];
caret => {};
end => {};
ENDCASE => ERROR Error[BadFormat];
caret => SELECT nextTokenType FROM
exponent => ChangePower[nextToken];
ENDCASE => ERROR Error[BadFormat];
exponent => SELECT nextTokenType FROM
sign => AddNewTerm[sign: nextToken];
variable => AddNewVariable[nextToken];
end => {};
ENDCASE => ERROR Error[BadFormat];
ENDCASE => ERROR Error[InternalError];
Now the nextTokenType becomes the previousTokenType.
previousTokenType ¬ nextTokenType;
ENDLOOP;
Convert from the temporary representation to a Ref.
StoreCurrentMonomial[];
intermRef ¬ RefFromTemporary[];
result ¬ Compact[RefFromTemporary[]];
END;
Operations on MultiPols
Sum: PUBLIC PROC [a, b: Ref] RETURNS [r: Ref] ~ BEGIN
r ¬ Compact[AddRefsWithoutCompaction[a, b]];
END;
Difference: PUBLIC PROC [a, b: Ref] RETURNS [r: Ref] ~ BEGIN
r ¬ Sum[a, Scale[b, -1]];
END;
Product: PUBLIC PROC [a, b: Ref] RETURNS [r: Ref] ~ BEGIN
r ¬ Compact[MultiplyRefsWithoutCompaction[a, b]];
END;
Power: PUBLIC PROC [a: Ref, n: NAT] RETURNS [r: Ref] ~ BEGIN
intermediate: Ref ¬ RefFromRope["1"];
FOR i: NAT IN [1..n] DO
intermediate ¬ Product[intermediate, a];
ENDLOOP;
r ¬ intermediate;
END;
Scale: PUBLIC PROC [p: Ref, c: REAL] RETURNS [r: Ref] ~ BEGIN
intermediate: Ref ¬ NEW[MultiPolRec[p.length]];
FOR i: NAT IN [0..p.length) DO
intermediate[i] ¬ [coefficient: p[i].coefficient * c, vars: p[i].vars];
ENDLOOP;
r ¬ Compact[intermediate];
END;
Differentiate: PUBLIC PROC [a: Ref, var: ATOM] RETURNS [r: Ref] ~ BEGIN
intermediate: Ref ¬ NEW[MultiPolRec[a.length]];
FOR i: NAT IN [0..a.length) DO
intermediate[i] ¬ DifferentiateMonomial[a[i], var];
ENDLOOP;
r ¬ Compact[intermediate];
END;
Degree: PUBLIC PROC [a: Ref] RETURNS [d: NAT] ~ BEGIN
degree: NAT ¬ 0;
FOR i: NAT IN [0..a.length) DO
degree ¬ MAX[degree, MonomialDegree[a[i]]];
ENDLOOP;
d ¬ degree;
END;
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];
END;
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 => IF reason=NotConstant THEN GOTO UVError];
result ¬ value;
EXITS UVError => ERROR Error[UnsubstitutedVariables];
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 => IF reason=NotUnivariate THEN GOTO UVError];
EXITS UVError => 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;
END;
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;
END;
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;
END.