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

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;


UnivariateFromRef: PUBLIC PROC [a: Ref, var: ATOM] RETURNS [result: Polynomial.Ref] ~ BEGIN

x: Polynomial.Ref _ Polynomial.Linear[[0, 1]];
intermediate: Polynomial.Ref;
degree: NAT _ Degree[a];

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;

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;

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];
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

intermRef: Ref;

inStream: IO.STREAM _ IO.RIS[rope];

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;

NextToken: PROC [stream: IO.STREAM _ inStream] RETURNS [Rope.ROPE] ~ BEGIN
token: Rope.ROPE;
[token: token] _ IO.GetTokenRope[stream, BreakProc];
RETURN[token];
END;

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;

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;

TokenType: TYPE ~ {sign, coefficient, variable, caret, exponent, begin, end, error};

previousTokenType: TokenType _ begin;
FlushTemporaryMultiPol[];
ClearCurrentMonomial[];
UNTIL previousTokenType = end DO

nextToken: Rope.ROPE;
firstCharOfNextToken: CHAR;
nextTokenType: TokenType;

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;

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];

previousTokenType _ nextTokenType;
ENDLOOP;

StoreCurrentMonomial[];
intermRef _ RefFromTemporary[];
result _ Compact[RefFromTemporary[]];

END;


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;
Evaluate: PUBLIC PROC [a: Ref, variable: ATOM, value: REAL] RETURNS [result: Ref] ~ BEGIN
break: SIGNAL = CODE;
intermediate: Ref _ NEW[MultiPolRec[a.length]];

FOR i: NAT IN [0..a.length) DO
thisMonomial: Monomial _ a[i];
newCoefficient: REAL;
newVars: REF PoweredVariableSequence;

variablePresent: BOOLEAN;
variablePosition: NAT;
[variablePresent, variablePosition] _ VariablePositionInMonomial[thisMonomial, variable];

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

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;

intermediate[i] _ [coefficient: newCoefficient, vars: newVars];
ENDLOOP;

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[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]];

FOR i: NAT IN [0..a.length) DO
thisMonomial: Monomial _ a[i];
newMonomial: Monomial;
newPolynomial: Ref;

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 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
ELSE BEGIN
newPolynomial _ NEW[MultiPolRec[1]];
newPolynomial[0] _ thisMonomial;
END;

intermediate _ AddRefsWithoutCompaction[intermediate, newPolynomial];
ENDLOOP;

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;


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

intermediate: Ref _ NEW[MultiPolRec[a.length]];
monomialsCopied: NAT _ 0;
finalIndex: NAT;

FOR i: NAT IN [0..a.length) DO
thisMonomial: Monomial _ a[i];
copied: BOOLEAN _ FALSE;

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;

IF copied = FALSE AND thisMonomial.coefficient # 0 THEN BEGIN
intermediate[monomialsCopied] _ thisMonomial;
monomialsCopied _ monomialsCopied + 1;
END;
ENDLOOP;

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;

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;

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 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

ELSE BEGIN
newCoefficient _ 0;
newVars _ NEW[PoweredVariableSequence[0]];
END;

result _ [coefficient: newCoefficient, vars: newVars];
END;

Equivalent: PROC [p, q: REF PoweredVariableSequence] RETURNS [result: BOOLEAN] ~ BEGIN

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

intermediate: REF PoweredVariableSequence _ NEW[PoweredVariableSequence[m.vars.numVars]];
lastVariableCopied: Rope.ROPE _ "";
lastVariableToCopy: Rope.ROPE _ "";
iIndex: NAT _ 0;

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;

UNTIL Rope.Compare[lastVariableCopied, lastVariableToCopy] = equal DO

varToCopy: ATOM;
variableToCopy: Rope.ROPE _ lastVariableToCopy;
exponent: NAT _ 0;

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];

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;

SELECT exponent FROM
<0 => ERROR Error[InternalError];
=0 => NULL;
>0 => {
intermediate[iIndex] _ [variable: varToCopy, degree: exponent];
iIndex _ iIndex + 1};
ENDCASE;

lastVariableCopied _ variableToCopy;
ENDLOOP;

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

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

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;

UNTIL (mIndex >= m.vars.numVars) AND (nIndex >= n.vars.numVars) DO
selection: Selection;

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;

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 m.coefficient = 0 OR a.length = 0 THEN BEGIN
result _ NEW[MultiPolRec[0]];
RETURN;
END;

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

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
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.
`MultiPolynomialImpl.mesa
Copyright (C) 1986 by Xerox Corporation.  All rights reserved.
James Rauen August 14, 1986 6:13:01 pm PDT


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


Conversions to and from other representations

Things needed to assemble the final result.
Make a sequence for the coefficients.
Add each term of a into the coefficient sequence.
Convert the coefficient sequence into a Polynomial.Ref.
Do the sign.  Treat the first term specially.
Declarations.
Convert the rope into a stream.
Define the break procedure used to separate the stream into tokens.
Define a procedure to get the next token from the stream.
These are the intermediate multipolynomial construction procedures.
These convert the intermediate representations to actual representations.
These are the kinds of tokens that might appear.
Parse the stream, one token at a time.
Declarations for the next token.
Get the next token and decide its type.  If the end of the stream is encountered, the token "&EOS" is used, with TokenType end.
Take action, depending on the type of next token and previous token.
Now the nextTokenType becomes the previousTokenType.
Convert from the temporary representation to a Ref.

Operations on MultiPols



Substitution Mechanisms

Iterate over each monomial term
See if the monomial contains the variable.
If it does, compute the evaluated monomial.
If it doesn't, the new monomial is just a copy of the old one.
Store the new monomial in the intermediate polynomial.
Compact the intermediate polynomial and return it.
Iterate over each monomial term
See if the monomial contains the variable.
If it does, compute newMonomial, the monomial without the variable, and newPolynomial, the product of newMonomial and the proper power of replacement.

If it doesn't, the new polynomial is just a copy of the old monomial.
Add the new polynomial to the intermediate polynomial.
Compact the intermediate polynomial and return it.

Internal Procedures
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.

Iterate over the monomials in a.
If an equivalent monomial has already been copied, add thisMonomial to it.
Otherwise, copy thisMonomial to the intermediate polynomial.
Squeeze out all the zero terms.
Make a more compact polynomial from the intermediate.
See if the variable appears in the monomial.
If it does, compute the derivative.
If it doesn't, the derivative is zero.
Create and return the new monomial.
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).

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
Figure out the last variable, alphabetically.
Copy each variable into the intermediate sequence.
Declarations
Figure out the next variable to copy.
Accumulate that variable's exponents.
Put the powered variable into the intermediate sequence, if its exponent is positive.
Save the variable as the last variable copied, and do the next variable.
Construct the final monomial.
A variable name is legal if it is a single alphabetic character.

Multiplies two monomials, preserving the term invariant.

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.
Determine if the next variable will be taken from m, n, or both.
Compress and construct the resulting monomial.
If either parameter is zero, return a zero polynomial.
Otherwise, multiply each monomial in a by m and store them in the result polynomial.  Note that no compaction is needed.
Uses the distributive law to multiply two Refs.  Each monomial in a is multiplied by b; the results are accumulated in the intermediate polynomial.
THIS BELONGS IN Polynomial MODULE
Κ˜™J™>J™*J™—Icode˜š
Οk	˜	Kšœœ[˜hKšœ
œ3œœœœ˜QKš
œ˜Kšœœ/˜?Kšœœ"œ˜2K˜—K˜šΠlnœœ
˜"K˜Kšœ
œ˜%Kšœ˜K˜Kšœ˜K˜—K™š
ž™K™K™ΩK™K™°K™—K™š
ž™K™Kšœœ˜ Kšœ
œ˜0Kšœ
œ˜*Kšœœ#˜8Kšœœ+˜HKšœœ%˜<Kšœœ&˜>Kšœœ'˜@Kšœœ(˜BKšœœ˜.Kšœ
œ˜0KšΟnœœ
œœ
˜.Kšœœ˜,K˜K˜—K™š
ž-™-K™šŸœœ
œœœ˜[K˜™+K˜.K˜Kšœœ
˜K˜—™%Kš
œ
œœ	œ	œ
œ
œ˜>Kšœœœ˜7Kš
œœ
œœœ
˜GK˜—™1šœœ
œ˜šœ˜˜K˜5—˜šœœ
˜#K˜X—Kšœ
œ˜ —˜Kšœ˜—Kšœ
˜—Kšœ
˜—K˜—™7Kš
œ(˜(šœœ
œ
˜Kš
œR˜RKš
œ3˜3Kš
œ=˜=Kšœ
˜—K˜K˜—Kšœ
˜K˜—šŸœœ
œœœ˜]Kšœœ˜%Kšœœ˜1Kšœœ
œ˜=šœœ
œ
˜Kšœœœ˜DK˜%Jš
œK˜KJšœ
˜—Jš
œ˜Kšœ
˜K˜—šŸœœ
œ
œœ˜Fšœ
˜Kšœœ˜Kšœœœ
œ˜<Kš
˜—Kšœ˜Kšœ
˜K˜—šŸœœ
œ
œ
œ˜@Kšœ
œ˜šœœ
œ˜K˜™-šœœ
˜Kšœ
œœ#˜@—š
˜Kšœ
œœ$˜AKšœ&˜*——Kšœ2œ˜Ršœœ
œ ˜/Kš
œ0˜0šœ˜K˜	KšœJœ˜RKšœKœk˜»
Kšœ
˜—Kšœ
˜—Kšœ
˜—Kšœœ˜$K˜Kšœ
˜K˜—šŸœœ
œ
œœ˜HK˜™
K˜K˜—™Kš	œ
œ
œœ
œ˜#—K˜™CšŸ	œœ
˜Kšœœ˜šœ˜Kš
œœ
œœ
œ˜*Kš
œ!˜!Kšœœ˜.Kšœ$˜&Kšœœ˜#—Kšœ
˜Kšœ
˜—K˜—™9š
Ÿ	œœ
œ
œ
œœ˜JKšœœ
˜Kšœœ!˜4Kšœ˜Kšœ
˜K˜——™Cšœœœ˜%Kšœœ
˜Kšœœ
œ˜&—šœœœ˜"Kšœœ
˜Kšœ
œ
˜Kšœœ
˜Kšœœ
œ˜+—Kš
œ(˜(Kš
œ#˜#šŸœœ˜&Kšœœ˜Kšœ
˜—šŸœœ˜$Kšœ œ˜%Kšœ
˜—šŸœœ˜$Kšœœ/˜QKšœ
˜—šŸ
œœœ˜)Kšœ)œ˜GK˜NJ˜Jš
œ˜Kšœ
˜—šŸœœ
œ˜1˜:Kšœœ˜+—Kšœ
˜—šŸœœœ˜,K˜Jšœ#œ˜)K˜1Kš
œ"˜"—Kšœ
˜—šŸœœ
œ˜+š
œK˜KKšœœ˜+—Kšœ
˜—K˜—™IšŸœœœ
˜NKšœœœ2˜YKšœ
œ˜%Kšœœ
œ)˜EKšœœ˜$šœœ
˜"K˜K˜(K˜/Kšœ
˜—Kšœ&œ˜IKšœ%˜+Kšœ
˜K˜—šŸœœ
œ˜+Kš
œ/˜/Kšœœ&˜7Kšœœ
œ$˜9Kšœœ˜$šœœ
˜K˜Kš
œ:˜:Kš
œ!˜!Kšœ
˜—Kšœ	˜šœ
˜K˜———™0KšœœE˜TK˜—™&K˜%K˜K˜šœ˜ K˜™ Kšœœ
˜Kšœœ
˜K˜—K˜™Kšœœ%œ˜JKš
œ0˜0šœœ˜0Kšœœ˜%Kšœœœ
œ
˜PK˜K˜
K˜Kšœ˜K˜——™Dšœ˜šœ	œ˜"K˜$K˜EK˜?Kšœœ˜"—šœœ˜!K˜,K˜&Kšœœ˜"—šœœ˜(K˜$K˜&K˜
Kšœœ˜"—šœœ˜%K˜$K˜&K˜K˜
Kšœœ˜"—šœ	œ˜"K˜#Kšœœ˜"—šœœ˜%K˜$K˜&K˜
Kšœœ˜"—Kšœœ˜&K˜——™4K˜"Kšœ
˜K˜———™3K˜Kš
œ˜Kš
œ%˜%K˜—Kšœ
˜K˜—K˜—K™š
ž™K™š	Ÿœœ
œ
œ˜5K˜,Kšœ
˜K˜—š	Ÿ
œœ
œ
œ˜<K˜Kšœ
˜K˜—š	Ÿœœ
œ
œ˜9K˜1Kšœ
˜K˜—šŸœœ
œ
œœ˜<K˜%šœœ
œ˜K˜(Kšœ
˜—K˜Kšœ
˜K˜—šŸœœ
œ
œœ˜=Kšœœ˜/šœœ
œ˜K˜GKšœ
˜—K˜Kšœ
˜K˜—šŸ
œœ
œœœ˜GKšœœ˜/šœœ
œ˜K˜3Kšœ
˜—K˜Kšœ
˜K˜—šŸœœ
œ
œœ˜5Kšœœ˜šœœ
œ˜Kšœ	œ˜+Kšœ
˜—K˜Kšœ
˜—K™—K™š
ž™K™š
Ÿœœ
œœ	œœ˜YKšœœœ
˜Kšœœ˜/K˜Kš
œ™šœœ
œ˜K˜Kšœœ
˜Kšœ	œ˜%K˜K™*Kšœœ
˜Kšœœ
˜K˜YK˜K™+šœœ
˜Kšœœ.˜9Kš
œH˜HKšœ
œ9˜Fšœœ
œ˜&K˜"Kšœ
˜—šœœ
œ/˜>K˜$Kšœ
˜—Kš
˜—K˜K™>šœ
˜
K˜*Kšœ
œ5˜Bšœœ
œ œ
˜0Kš
œ"˜"Kšœ
˜—Kšœ
˜—K˜K™6K˜?Kšœ
˜—K˜K™2K˜šœ
˜K˜——š	Ÿœœ
œ(œ˜^K˜+K˜šœœ
˜K˜[K˜Kšœ
˜—K˜Kšœ
˜K˜—šŸœœ
œ(œ
œ˜bK˜.Kšœœ
˜˜ Kšœœ ˜=—K˜Kšœ
˜K˜—šŸ
œœ
œœœ˜`Kšœœ˜)K˜Kš
œ™šœœ
œ˜K˜K˜K˜K˜K™*Kšœœœ
˜Kšœœ
˜šœœ
œ œ
˜0šœ*œ˜2Kšœœ
˜K˜—Kšœ
˜—K˜K™–
šœœ
˜Kšœœ.˜9Kšœ#œ˜7K˜.K˜DK˜OKš
˜—K™K™Ešœ
˜
Kšœœ˜$K˜ Kšœ
˜—K˜K™6K˜EKšœ
˜—K˜K™2K˜Kšœ
˜K˜—š	Ÿœœ
œ*œ˜bKš
œ-˜-K˜šœœ
˜˜K˜K˜K˜—K˜Kšœ
˜—K˜Kšœ
˜K˜—š	Ÿœœ
œ!œ˜rKš
œ$˜$K˜šœœ
˜˜K˜
K˜K˜4K˜—K˜Kšœ
˜—š
œ+˜+Kšœœ ˜?—Kšœ
˜K˜——K™š
ž™K˜šŸœœ
œ˜HKšœ	œ#˜/šœœ
œ˜K˜Kšœ
˜—šœœ
œ˜K˜Kšœ
˜—Kšœ
˜K˜—šŸœœ
œ˜4K˜K™­K™Kšœœ˜/Kšœœ˜Kšœœ
˜K˜K™ šœœ
œ˜K˜Kšœœœ
˜K˜K™Jšœœ
œ˜%šœ5œ
˜AK˜UKšœ	œ
˜Kšœ
˜—Kšœ
˜—K˜K™<š	œ
œ
œœ
˜=K˜-K˜&Kšœ
˜—Kšœ
˜—K˜K™K˜šœœ
œ˜%šœ!œ
˜-K˜+K˜Kšœ
˜—Kšœ
˜K˜—K™5Kšœ	œ˜&šœœ
œ˜ K˜Kšœ
˜K˜—šœ
˜K˜——šŸœœœ˜CK˜Kšœ	œœ*˜Tšœœ
œ˜$K˜FKšœ
˜—K˜5Kšœ
˜K˜—š	Ÿœœœœ˜\K˜Kšœœ
˜Kšœ	œ˜%K˜K™,Kšœœœ
˜Kšœœ
˜šœœ
œ˜$šœœ
˜+Kšœœ
˜K˜Kšœ
˜—Kšœ
˜—K˜K™#šœœ
˜Kšœœ#˜.K˜(Kšœ
œ*˜7šœœ
œ˜$K˜Kšœ
˜—Kšœ/Οc,˜[Kš
˜—K˜K™&šœ
˜
K˜Kšœ
œ˜*Kšœ
˜—K˜K™#K˜6Kšœ
˜K˜—šŸ
œœœœ
œ˜VK˜K™›
K™Kšœœœœ˜7šœœ
œ˜šœ!œœ
˜LKšœ	œ
˜Kšœ
˜Kšœ
˜—Kšœ
˜—Kšœ	œ
˜Kšœ
˜K˜—šŸ
œœœ˜DK˜K™K™K™σ
K™™Kšœœœ*˜YKšœœ˜#Kšœœ˜#Jšœœ˜K˜—™-šœœ
œ˜$Kšœœ,œ˜EKšœ5œ˜YKšœ
˜K˜——™2šœ>˜EK˜™Kšœœ
˜Kšœœ˜/Kšœ
œ˜—K˜™%šœœ
œ˜$Kšœœ,œ˜Ešœ2œ.œ
˜oKš
œ˜Kšœ
˜—Kšœ
˜—š
œ1˜1K˜——™%šœœ
œ˜$Kšœœ,œ˜Ešœ/˜5Kš
œ'˜'—Kšœ
˜K˜——™Ušœ
˜Kšœœ˜!Kšœœ
˜˜Kš
œ?˜?K˜—Kšœ
˜—K˜—™HKš
œ$˜$Kšœ
˜K˜———™˜
K˜Kšœœ#˜,—šœœ
œ
˜K˜!Kšœ
˜—šœ
˜K˜———šŸœœœœ	œ˜EK˜K™@K™Kšœ(œ˜;Kšœ
˜K˜—š	Ÿœœœ
œ˜@Kšœœ˜šœœ
œ˜$K˜Kšœ
˜—K˜Kšœ
˜K˜—šŸœœœ˜KK˜K™8K™Kšœœ˜+Kšœœœ;˜jKšœœ˜(Kšœœ˜Kšœœ˜Kšœœ˜šœœœ
˜4Kšœ
œ˜-K˜,Kšœ
˜Kšœ
˜—K˜Kš
œβ
™β
šœœ˜BK˜K˜K™@Kšœœ˜2Kšœ
œœ˜7šœ
˜
Kšœœ1˜IKšœœ1˜Išœ,˜6K˜K˜K˜Kšœ
˜—Kšœ
˜—K˜šœ˜šœ	˜K˜&K˜Kšœ
˜—šœ	˜K˜&K˜Kšœ
˜—šœ˜K˜8K˜LK˜K˜Kšœ
˜—Kšœ
˜—K˜Kšœ
˜—K˜K™.Kšœ
œ"˜2šœœ
œ
˜K˜ Kšœ
˜—K˜HKšœ
˜K˜—šŸœœœ˜OK˜K™6šœœœ
˜/Kšœ	œ˜Kšœ
˜Kšœ
˜—K˜K™xKšœ	œ˜$šœœ
œ˜K˜'Kšœ
˜—Kšœ
˜K˜—šŸœœ
œ˜MK˜K™“
Kšœœ˜(šœœ
œ˜K˜9K˜JKšœ
˜—K˜Kšœ
˜K˜—šŸœœ œœ	œ˜]K™!K˜Kšœ
˜K˜—š
Ÿœœœœœœ˜{Kšœœ˜Kšœœ
˜šœœ˜.šœœ˜"Kšœœ
˜K˜—K˜
Kšœ
˜—Kšœ
˜K˜—š
Ÿœœœœœ
œ˜Qšœ
˜K˜K˜K˜K˜Kšœœ ˜/Kšœœ%˜4Kšœœ$˜3KšœœC˜SKšœ
˜—Kšœ
˜——K˜Kšœ
˜—…—R|‚ς