File: Algebra3dImpl.mesa
Last edited by Bier on December 18, 1982 1:37 am
Author: Eric Bier on August 4, 1983 3:36 pm
Contents: The quadratic formula, cubic formula and quartic formula
DIRECTORY
Algebra3d,
Complex,
Matrix3d,
RealFns,
SV2d,
SVVector3d;
Algebra3dImpl:
PROGRAM
IMPORTS Complex, RealFns
EXPORTS Algebra3d =
BEGIN
Point2d: TYPE = Matrix3d.Point2d;
Polygon: TYPE = REF PolygonObj;
PolygonObj: TYPE = SV2d.PolygonObj;
Vec: TYPE = Complex.Vec;
Vector: TYPE = SVVector3d.Vector;
GLOBALS
root3Over2: REAL; -- a loadtime constant
LinearFormula:
PRIVATE
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:
PRIVATE
PROC [vec: Vec]
RETURNS [r:
REAL, radians:
REAL] = {
r ← Complex.Abs[vec];
radians ← Complex.Arg[vec];
};
ComplexRoot:
PRIVATE
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.