File: Algebra3dImpl.mesa
Last edited by Bier on February 18, 1987 6:19:57 pm PST
Copyright © 1984 by Xerox Corporation. All rights reserved.
Contents: The quadratic formula, cubic formula and quartic formula

DIRECTORY
Algebra3d, Complex, RealFns, Vector2;

Algebra3dImpl: CEDAR PROGRAM
IMPORTS Complex, RealFns
EXPORTS Algebra3d =

BEGIN
VEC: TYPE = Vector2.VEC;
GLOBALS
root3Over2: REAL; -- a loadtime constant

LinearFormula: PROC [a, b: REAL] RETURNS [root: REAL, rootCount: NAT] = {
The solution to the equation "ax + b = 0".
IF a = 0 THEN {
IF b = 0 THEN {root ← 0.0; rootCount ← 0; RETURN}
ELSE {root ← 0.0; rootCount ← 0; RETURN} -- inconsistent case
}
ELSE {
IF b = 0 THEN {root ← 0.0; rootCount ← 1; RETURN}
ELSE {root ← -b/a; rootCount ← 1; RETURN};
};
};
QuadraticFormula: PUBLIC PROC [a, b, c: REAL] RETURNS [roots: ARRAY [1..2] OF REAL, rootCount: NAT] = {
The solution to the equation "ax2+bx+c=0". If a=0 this is just a linear equation. If
c = 0, one root is zero and we solve a linear equation. Otherwise, we use the
quadratic formula in either the form [-b+-(b2-4ac)1/2]/2a or (2c)/[-b-+(b2-4ac)1/2]
depending on the sign of b. We will arrange the roots so that roots[1] < roots[2].
discriminant, temp: REAL;
IF a = 0 THEN {
[roots[1], rootCount] ← LinearFormula [b, c];
RETURN;
};
IF c = 0 THEN {
roots[1] ← 0.0;
[roots[2], rootCount] ← LinearFormula[a, b];
rootCount ← rootCount + 1;
IF roots[1] > roots[2] THEN { -- swap them
temp ← roots[1];
roots[1] ← roots[2];
roots[2] ← temp;
};
RETURN};
discriminant ← b*b - 4*a*c; -- 3 mult, 1 add
SELECT discriminant FROM
=0 => {rootCount ← 1;
roots[1] ← -b/(2*a);
RETURN}; -- 1 mult, 1 div (4 mult, 1 div, 1 add total)
<0 => {rootCount ← 0;
RETURN}; -- (3 mult, 1 add total)
>0 => {
sqrtRadical: REAL ← RealFns.SqRt[discriminant];
term: REAL;
rootCount ← 2;
term ← IF b < 0 THEN sqrtRadical - b ELSE -sqrtRadical - b;
roots[1] ← term/(2.0*a);
roots[2] ← (2.0*c)/term;
IF roots[1] > roots[2] THEN { -- swap them
temp ← roots[1];
roots[1] ← roots[2];
roots[2] ← temp;
};
RETURN}; -- 1 SqRt, 2 mult, 2 div, 2 add (1 SqRt, 5 mult, 2 div, 3 add total)
ENDCASE => ERROR;
}; -- end of QuadraticFormula

ToPolar: PROC [vec: VEC] RETURNS [r: REAL, radians: REAL] = {
r ← Complex.Abs[vec];
radians ← Complex.Arg[vec];
};
ComplexRoot: PROC [index: REAL, vec: VEC] RETURNS [rootVec: VEC] = {
r, radians, resultR, resultRadians: REAL;
[r, radians] ← ToPolar[vec];
resultR ← RealFns.Root[index, r];
resultRadians ← radians/index;
rootVec ← Complex.FromPolar[resultR, resultRadians];
};
Init: PROC = {
root3Over2 ← RealFns.SqRt[3]/2.0;
};
Init[];
END.