[Cyan]<ISL>SolidViews>Algebra3dImpl.mesa!1

File: Algebra3dImpl.mesa

Last edited by Bier on December 18, 1982 1:37 am

Author: Eric Bier on January 4, 1983 12:51 am

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

QuadraticFormula: PUBLIC PROC [a, b, c: REAL] RETURNS [roots: ARRAY [1..2] OF REAL, rootCount: NAT] = {

radical: REAL ← b*b - 4*a*c; -- 3 mult, 1 add

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

a2: REAL ← 2*a;

rootCount ← 2;

roots[1] ← (-b + sqrtRadical)/a2;

roots[2] ← (-b - sqrtRadical)/a2;

RETURN}; -- 1 SqRt, 1 mult, 2 div, 2 add (1 SqRt, 4 mult, 2 div, 3 add total)

ENDCASE => ERROR;

}; -- end of QuadraticFormula

CubicFormula: PUBLIC PROC [a3, a2, a1, a0: REAL] RETURNS [roots: ARRAY [1..3] OF REAL, rootCount: NAT] = {

For polynomial a3x^3 + a2x^2 + a1x + a0 = 0.
Make leading coefficient = 1;
From "Elementary Analytical Methods":
Given z^3 + a2z^2 + a1z + a0 = 0, let
q = 1/3a1-1/9a2^2; r = 1/6(a1a2-3a0)-1/27a2^3
If q^3 + r^2 >0, one real root and a pair of complex conjugate roots,
q^3 + r^2 = 0, all roots real and at least two are equal,
q^3 + r^2<0, all roots real (irreducible case).
Let
s1 = [r + (q^3+r^2)^(1/2)]^(1/3), s2 = [r - (q^3 + r^2)^(1/2)]^(1/3)
then
z1 = (s1 + s2)-a2/3
z2 = -(1/2)(s1+s2)-a2/3 + iRoot(3)/2(s1-s2)
z3 = -(1/2)(s1+s2)-a2/3 - iRoot(3)/2(s1-s2)
[To test answer:]
If z1, z2, z3 are the roots of the cubic equation [then]
z1 + x2 + z3 = -a2
z1z2 + z1z3 + z2z3 = a1
z1z2z3 = -a0

q, r: REAL;

q3PlusR2, a2Over3, a2a2Over9: REAL;

IF a3 = 0 THEN {

quadRoots: ARRAY [1..2] OF REAL;

[quadRoots, rootCount] ← QuadraticFormula [a2, a1, a0];

FOR i: NAT IN [1..rootCount] DO

roots[i] ← quadRoots[i];

ENDLOOP;

RETURN};

a2 ← a2/a3;

a1 ← a1/a3;

a0 ← a0/a3;

a2Over3 ← a2/3;

a2a2Over9 ← a2Over3*a2Over3;

q ← a1/3 - a2a2Over9;

r ← (a1*a2-3*a0)/6 - a2a2Over9*a2Over3;

q3PlusR2 ← q*q*q + r*r;

SELECT q3PlusR2 FROM

>0 => {-- one real root and a pair of complex conjugate roots

s1, s2: REAL;

rootQ3PlusR2, rMinusRootQ3PlusR2, rPlusRootQ3PlusR2: REAL;

rootQ3PlusR2 ← RealFns.SqRt[q3PlusR2];

rMinusRootQ3PlusR2 ← r - rootQ3PlusR2;

rPlusRootQ3PlusR2 ← r + rootQ3PlusR2;

SELECT rPlusRootQ3PlusR2 FROM

<0 => s1 ← -RealFns.Root[3, -rPlusRootQ3PlusR2];

=0 => s1 ← 0;

>0 => s1 ← RealFns.Root[3, rPlusRootQ3PlusR2];

ENDCASE => ERROR;

SELECT rMinusRootQ3PlusR2 FROM

<0 => s2 ← -RealFns.Root[3, -rMinusRootQ3PlusR2];

=0 => s2 ← 0;

>0 => s2 ← RealFns.Root[3, rMinusRootQ3PlusR2];

ENDCASE => ERROR;

roots[1] ← (s1 + s2) - a2Over3;

rootCount ← 1};

=0 => {-- all roots real and at least two are equal

s1 = s2 = r^(1/3); Hence s1 - s2 = 0;

s1: REAL;

s1PlusS2: REAL;

s1 ← RealFns.Root[3, r];

s1PlusS2 ← s1+s1;

roots[1] ← s1PlusS2 - a2Over3;

roots[2] ← -s1PlusS2/2;

roots[3] ← roots[2];

rootCount ← 3};

<0 => {-- all roots real (irreducible case)

The problem here is that s1 and s2 are complex. We represent them as vectors. They are complex conjugates so we only need to calculate one of them.

s1: Vec;

rootQ3PlusR2, s1PlusS2, s1MinusS2: REAL;

rootQ3PlusR2 ← RealFns.SqRt[-q3PlusR2]; -- but remember this is pure complex

s1 ← ComplexRoot[3, [r, rootQ3PlusR2]];

s1PlusS2 ← s1.x + s1.x;

s1MinusS2 ← s1.y + s1.y; -- but remember this is pure complex

roots[1] ← s1PlusS2 - a2Over3;

roots[2] ← -s1PlusS2/2-a2Over3-root3Over2*s1MinusS2;

roots[3] ← -s1PlusS2/2-a2Over3+root3Over2*s1MinusS2;

rootCount ← 3;

};

ENDCASE => ERROR;

};

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

};

QuarticFormula: PUBLIC PROC [a4, a3, a2, a1, a0: REAL] RETURNS [roots: ARRAY [0..3] OF REAL, rootCount: NAT] = {

For polynomial a4x^4 + a3x^3 + a2x^2 + a1x + a0 = 0.

From "Elementary Analytical Methods":
Given z^4 + a3z^3 + a2z^2 + a1z + a0 = 0, find the real root u1 of the cubic equation:

u^3 -a2u^2 + (a1a3 - 4a0)u - (a1^2 + a0a3^2 - 4aoa2) = 0

and determine the four roots of the quartic as solutions of the two quadratic equations

v^2 + [a3/2 +-(a3^2/4 + u1 - a2)^(1/2)]v + u1/2+-[(u1/2)^2-a0]^(1/2) = 0

If all roots of the cubic equation are real, use the value of u1 which gives real cooefficients in the *quadratic equation and select signs so that if

z^4 + a3z^3 + a2z^2 + a1z + a0 = (z^2 + p1z+ q1)(z^2 + p2z + q2)

then

p1 + p2 = a3, p1p2 + q1 + q2 = a2, p1q2 + p2q1 = a1, q1q2 = a0.
We begin by normalizing the polynomial.

cubicRoots: ARRAY [1..3] OF REAL;

cubicRootCount: NAT;

a3Over2, a3a3Over4, u1, u1Over2: REAL;

quadRoots1, quadRoots2: ARRAY [1..2] OF REAL;

quadRootCount1, quadRootCount2: NAT;

term1, term2: REAL;

IF a4 = 0 THEN {

[cubicRoots, cubicRootCount] ← CubicFormula [a3, a2, a1, a0];

FOR i: NAT IN [1..cubicRootCount] DO

roots[i-1] ← cubicRoots[i];

ENDLOOP;

rootCount ← cubicRootCount;

RETURN};

a3 ← a3/a4; a2 ← a2/a4; a1 ← a1/a4; a0 ← a0/a4;

We continue by finding the roots of the cubic equation:

a3Over2 ← a3/2;

a3a3Over4 ← a3Over2*a3Over2;

[cubicRoots, cubicRootCount] ← CubicFormula[1, -a2, (a1*a3 - 4*a0), -(a1*a1 + a0*a3*a3 - 4*a0*a2)];

SELECT cubicRootCount FROM

1 => u1 ← cubicRoots[1];

2 => ERROR;

3 => {-- the hard case

Of the three roots, choose the first one which satifies both:
u1 >= a2-a3^2/4
(u1/2)^2 >= a0

a2Minusa3a3Over4: REAL ← a2 - a3a3Over4;

SELECT TRUE FROM

cubicRoots[1] >= a3a3Over4 AND cubicRoots[1]*cubicRoots[1]/4 >= a0 =>
u1 ← cubicRoots[1];

cubicRoots[2] >= a3a3Over4 AND cubicRoots[2]*cubicRoots[2]/4 >= a0 =>
u1 ← cubicRoots[2];

cubicRoots[3] >= a3a3Over4 AND cubicRoots[3]*cubicRoots[3]/4 >= a0 =>
u1 ← cubicRoots[3];

ENDCASE => ERROR;

};

ENDCASE => ERROR;

Finally we solve the two quadratic equations:

u1Over2 ← u1/2.0;

term1 ← RealFns.SqRt[a3a3Over4+u1-a2];

term2 ← RealFns.SqRt[u1Over2*u1Over2-a0];

[quadRoots1, quadRootCount1] ←
QuadraticFormula[1, a3Over2 + term1, u1Over2 + term2];

[quadRoots2, quadRootCount2] ←
QuadraticFormula[1, a3Over2 - term1, u1Over2 - term2];

rootCount ← 0;

IF quadRootCount1>0 THEN {roots[rootCount] ← quadRoots1[1];

rootCount ← rootCount + 1};

IF quadRootCount1>1 THEN {roots[rootCount] ← quadRoots1[2];

rootCount ← rootCount + 1};

IF quadRootCount2>0 THEN {roots[rootCount] ← quadRoots2[1];

rootCount ← rootCount + 1};

IF quadRootCount2>1 THEN {roots[rootCount] ← quadRoots2[2];

rootCount ← rootCount + 1};

}; -- end of Quartic Formula

Init: PROC = {

root3Over2 ← RealFns.SqRt[3]/2.0;

};

Init[];

END.