Adventuring - The Bear Creek Commentaries

--- DANSYSX: Complex cubic equations ---

2016

Both LibreOffice Calc and DANSYS are able to do basic complex arithmetic, but I want DANSYSX to be able to do advanced complex mathematics. There are a number of algorithms on Jean Pierre Moreau's excellent website (http://jean-pierre.moreau.pagesperso-orange.fr/index.html) that I want to adapt to this purpose.

The first is an algorithm to determine the roots of a cubic equation with complex coefficients. It uses the cubic formula (like the quadratic formula, but for cubic functions):

z(1...3) = cube root (-q/2 + square root(q^2/4 + p^3/27)) + cube root(-q/2 - square root(q^2/4 + p^3/27))

where: p = (c/a) - (b^2/3a^2)
           q = (2b^3/27a^3) + (d/a) - (bc)/3a^2)

I want the complex mathematics programs to be in their own folder so I use Tools>Macros>Organize Macros>LibreOffice Basic>Organizer to create the new folder, Complex under DANSYSX.ods>Standard. To begin programming, I open the LibreOffice Basic Macros dialog and select DANSYSX.ods>Standard>Complex and click the Edit button.

I am more comfortable programming in an environment where indexes run from 1 to whatever (instead of beginning at 0), so I add the Option Base 1 to the top of the editor Window, and change the name of the blank macro from Main to CompCube. I will also want the Macro to be a function so I change the header from Sub to Function and the result looks like this:

REM  *****  BASIC  *****

Option Base 1
Function CompCube

End Function

Now, everything is ready for me to start programming. I have to decide what I want my program to look like - the inputs and outputs. And I want to add a descriptive comment to the top of the code.

Looking at the Moreau algorithm, I see that I will need to input four complex numbers and I want to stick to the LibreOffice convention of designating complex numbers as strings of the form a±bi, so I want to pass four strings to the function. I will also want to output three complex numbers because a cubic equation will have up to three unique roots. So I'll start with the following header:

REM  *****  BASIC  *****

Option Base 1
Function CompCube(Z1 AS STRING, Z2 AS STRING, Z3 AS STRING, Z4 AS STRING)
'The following code determines the roots of a cubic equation
'defined by four complex coefficients:
'a z^3 + b z^2 + c z + d = 0
'The algorithm is an adaptation of one (CRoot3) from the Jean Pierre Moreau
'website: http://jean-pierre.moreau.pagesperso-orange.fr/Basic/croot3_bas.txt
'Several subroutines are used to perform elementary complex
'calculations and are shared with other DANSYSX programs.
'The global variables CMPN 1, 2, ..., 6 are also used to
'transfer complex variables between subroutines.

End Function

Most of the code will be in upper caps. This is not required by LibreOffice Basic, but I have found that I catch errors a lot quicker if everything is capitalized.

The Function statement names the program CompCube. The "Comp" part is a convention I established in DANSYS. Functions that deal with complex numbers begin with "Comp". The variables passed to the function (also called "parameters" or "arguments") are defined in the parentheses following the function name.

One of the flaws in LibreOffice Basic (all programming languages have weak spots) is that you can't transfer matrix values to and from functions and subroutines. (When I say "can't", I mean that I looked in all the resources I could find and tried to figure it out and couldn't - there may be a way but it remains a mystery to me.) I will have to transfer two value matrices (representing complex numbers), so I will store each value in a 1x2 array called CMPN1, CMPN2, to CMPN6. These arrays are defined as global variables (variables that can be used by all routines in DANSYSX and that won't lose their values between calls) which as designated in the Structures folder as, for instance:

Global CMPN1(1,2)

This generates a global array that has 1 row and 2 columns.

Global variables have to be defined outside of any function or subroutine. I collect them in the Structures folder.

I will also have to output three strings representing complex numbers so, to call this function from a spreadsheet, the user will select a 3x1 range of cells and type "=CompCube()' Inside the parentheses will either be four strings or the addresses of four cells containing strings.

I have written programs before that take complex numbers and have saved the code to do it in a KeyNote file of other handy code snippets that I use over and over. I will do that in the initialization section of the program.

Now, I go through the original code and figure out all the variables that I will need and declare them. At first guess, I come up with:

DIM A(2), B(2), C(2), D(2), U(2), V(2), W(2)
DIM DET(2), P(2), Q(2), SDET(2), U1(2), U2(2)
DIM V1(2), V2(2), ZT1(2), ZT2(2), Z(9,2)

The DIM statement is simple. It is just DIM followed by a list of variables separated by commas.

All the variables but the last are arrays with two places. Those are complex numbers.

Using my tried and true complex number decoders in the initialization section to load A, B, C, and D from Z1, Z2, Z3, and Z4, I wrote:

DIM A(2), B(2), C(2), D(2), U(2), V(2), W(2)
DIM DET(2), P(2), Q(2), SDET(2), U1(2), U2(2)
DIM V1(2), V2(2) ZT1(2), ZT2(2), Z(9,2)
DIM I AS INTEGER, ST AS STRING
DIM OutMat(3,1) AS STRING

OutMat is a convention I use to output arrays. I load all the data I want to display into an array and I set the function equal to that array. In this case, I'll be displaying three complex number strings - 3 rows and 1 column.

Notice that the last two variables are explicitly specified to be integers. LibreOffice allows you to specify variables as INTEGER (between values of -32768 and 32767), LONG (integers between -2147483648 and 2147483647), SINGLE (floating-point values between 3.402823E38 and 1.401298R-45), DOUBLE (floating-point values between 1.79769313486232E308 and 4.94065645841247E-324), CURRENCY (64 bit numbers displayed as fixed decimal numbers with 15 non-whole and 4 decimal places ranging from -922337203685477.5808 to 922337203685477.5807), STRING (character strings holding up to 65,535 characters), BOOLEAN (stores TRUE/1 or FALSE/0), and DATE (date and time values stores in an internal format.) Constants can also be declared using the keyword CONST (instead of DIM).

Variable types can also be specified by tagging a declaration symbol to the end of a variable name as follows:

Single !
Double #
Currency @
String $
For instance, DIM Money@ would declare the variable Money as the type Currency.

ST=""
FOR I=1 TO LEN(Z1)
IF MID(Z1,I,1)="+" OR MID(Z1,I,1)="-" THEN
A(1)=VAL(ST)
ST=""
END IF
ST=ST & MID(Z1,I,1)
NEXT I
A(2)=VAL(ST)

ST=""
FOR I=1 TO LEN(Z2)
IF MID(Z2,I,1)="+" OR MID(Z2,I,1)="-" THEN
B(1)=VAL(ST)
ST=""
END IF
ST=ST & MID(Z2,I,1)
NEXT I
B(2)=VAL(ST)

ST=""
FOR I=1 TO LEN(Z3)
IF MID(Z3,I,1)="+" OR MID(Z3,I,1)="-" THEN
C(1)=VAL(ST)
ST=""
END IF
ST=ST & MID(Z3,I,1)
NEXT I
C(2)=VAL(ST)

ST=""
FOR I=1 TO LEN(Z4)
IF MID(Z4,I,1)="+" OR MID(Z4,I,1)="-" THEN
D(1)=VAL(ST)
ST=""
END IF
ST=ST & MID(Z4,I,1)
NEXT I
D(2)=VAL(ST)

MSGBOX A(1)
MSGBOX A(2)
MSGBOX B(1)
MSGBOX B(2)
MSGBOX C(1)
MSGBOX C(2)
MSGBOX D(1)
MSGBOX D(2)

Each block of code beginning with ST works like this: First a blank string is generated in the string variable ST. Then, the following loop looks at each character of an input string Z. If the character is a + or a -, the first cell of a two place array is loaded with the numerical value of whatever is in ST and ST is emptied. If the character is anything else, it is tagged onto the right end of the string in ST and the next character comes up for review.  This goes on until each character of the input string (LEN(Z) is the number of characters in string Z) has been considered, then the numerical value of whatever is left in ST is loaded into the second cell of the two place variable.

The MID(Z,I,1) extracts 1 character from the ith position of string Z.

MSGBOX is a command that causes the following to be displayed in a message box. In this case, I use it to make sure that the input values are being transferred to the complex value arrays. They work so I'm happy. The MSGBOX commands will be deleted.

The next step requires that I add a subroutine that performs complex division and one that performs complex multiplication. I design these so that they work on the complex numbers stored in CMPN1 and CMPN2 and return an answer in CMPN3. You notice that the CMPN global variables act like a stack of registers for complex numbers. Here are CDIV and CMUL:

SUB CDIV()
'Subroutine to divide two complex numbers.
DIM CC(2), D AS DOUBLE, Z(2)

'Determine the value of the denominator
D=CMPN2(1,1)*CMPN2(1,1)+CMPN2(1,2)*CMPN2(1,2)
IF D<1e-15 THEN
MSGBOX "Complex divisiion by 0"
EXIT SUB
ELSE
CC(1)=CMPN2(1,1)
CC(2)=-CMPN2(1,2)
CMPN4(1,1)=CMPN1(1,1)
CMPN4(1,2)=CMPN1(1,2)
CMPN5(1,1)=CMPN2(1,1)
CMPN5(1,2)=CMPN2(1,2)
CMPN2(1,1)=CC(1)
CMPN2(1,2)=CC(2)
CMUL
Z(1)=CMPN3(1,1)
Z(2)=CMPN3(1,2)
CMPN2(1,1)=CMPN5(1,1)
CMPN2(1,2)=CMPN5(1,2)
CMPN3(1,1)=CMPN3(1,1)/D
CMPN3(1,2)=CMPN3(1,2)/D
END IF
END SUB

SUB CMUL()
'DETERMINE THE PRODUCT OF CMPN1 AND CMPN2
'SAVE THE RESULT IN CMPN3
CMPN3(1,1)=CMPN1(1,1)*CMPN2(1,1)-CMPN1(1,2)*CMPN2(1,2)
CMPN3(1,2)=CMPN1(1,1)*CMPN2(1,2)+CMPN1(1,2)*CMPN2(1,1)
END SUB

You will recognize most of the arithmetic operators. * is LibreOffice's multiplication operator (instead of x) and an exponent is denoted by the symbol ^. For instance, x squared would look like: x^2.

The difference between a function and a sub (or subroutine) is that the function returns a value and a sub does something but does not output a value. For instance, CMUL calculates the product of two complex numbers and loads them into a global array, but does not return the product to the program that called the sub. That program has to get the value from the global array. Here's the part of CompCube that uses CDIV and retrieves the result:

CMPN1(1,1)=C(1)
CMPN1(1,2)=C(2)
CMPN2(1,1)=A(1)
CMPN2(1,2)=A(2)
CDIV
U(1)=CMPN3(1,1)
U(2)=CMPN3(1,2)
MSGBOX U(1)
MSGBOX U(2)

This divides the complex number in C by the one in A. In order for CDIV to work C has to be loaded into CMPN1 and A has to be loaded into CMPN2. To call CDIV, I simply type CDIV. Then to get the result back as U, I load CMPN3 into U. Again, the MSGBOX commands display the result and, after I'm satisfied that CDIV works, I'll delete them. That's called "debugging". I prefer to debug as I go.

Now, to complete the calculation of P:

CMPN1(1,1)=B(1)
CMPN1(1,2)=B(2)
CMPN2(1,1)=B(1)
CMPN2(1,2)=B(2)
CMUL
V(1)=CMPN3(1,1)
V(2)=CMPN3(1,2)

CMPN1(1,1)=A(1)
CMPN1(1,2)=A(2)
CMPN2(1,1)=A(1)
CMPN2(1,2)=A(2)
CMUL
W(1)=CMPN3(1,1)
W(2)=CMPN3(1,2)

W(1)=3*W(1)
W(2)=3*W(2)

CMPN1(1,1)=V(1)
CMPN1(1,2)=V(2)
CMPN2(1,1)=W(1)
CMPN2(1,2)=W(2)
CDIV
ZT1(1)=CMPN3(1,1)
ZT1(2)=CMPN3(1,2)
P(1)=U(1)-ZT1(1)
P(2)=U(2)-ZT1(2)

You should be able to follow that. It's the same process I used to calculate c/a. It checks out, so I can continue.

The next piece of the puzzle is Q. It's a straight forward but long calculation (you see why they just teach the quadratic formula in high school and leave the cubic and quartic monsters for folks that are too curious for their own good.) That section of code looks like:

'Calculate q=2b^3/27a^3 + d/a - bc/3a^2
CMPN1(1,1)=B(1)
CMPN1(1,2)=B(2)
CMPN2(1,1)=V(1)
CMPN2(1,2)=V(2)
CMUL
W(1)=CMPN3(1,1)
W(2)=CMPN3(1,2)

W(1)=2*W(1)
W(2)=2*W(2)

CMPN1(1,1)=A(1)
CMPN1(1,2)=A(2)
CMPN2(1,1)=A(1)
CMPN2(1,2)=A(2)
CMUL
U(1)=CMPN3(1,1)
U(2)=CMPN3(1,2)

CMPN1(1,1)=U(1)
CMPN1(1,2)=U(2)
CMPN2(1,1)=A(1)
CMPN2(1,2)=A(2)
CMUL
V(1)=CMPN3(1,1)
V(2)=CMPN3(1,2)

V(1)=27*V(1)
V(2)=27*V(2)

CMPN1(1,1)=W(1)
CMPN1(1,2)=W(2)
CMPN2(1,1)=V(1)
CMPN2(1,2)=V(2)
CDIV
Q(1)=CMPN3(1,1)
Q(2)=CMPN3(1,2)

CMPN1(1,1)=D(1)
CMPN1(1,2)=D(2)
CMPN2(1,1)=A(1)
CMPN2(1,2)=A(2)
CDIV
U(1)=CMPN3(1,1)
U(2)=CMPN3(1,2)

Q(1)=Q(1)+U(1)
Q(2)=Q(2)+U(2)

CMPN1(1,1)=B(1)
CMPN1(1,2)=B(2)
CMPN2(1,1)=C(1)
CMPN2(1,2)=C(2)
CMUL
U(1)=CMPN3(1,1)
U(2)=CMPN3(1,2)

CMPN1(1,1)=A(1)
CMPN1(1,2)=A(2)
CMPN2(1,1)=A(1)
CMPN2(1,2)=A(2)
CMUL
V(1)=CMPN3(1,1)
V(2)=CMPN3(1,2)

V(1)=3*V(1)
V(2)=3*V(2)

CMPN1(1,1)=U(1)
CMPN1(1,2)=U(2)
CMPN2(1,1)=V(1)
CMPN2(1,2)=V(2)
CDIV
W(1)=CMPN3(1,1)
W(2)=CMPN3(1,2)

Q(1)=Q(1)-W(1)
Q(2)=Q(2)-W(2)

It's very repetitive and allows me to copy, paste,and modify most of it from earlier code. Then there is det:


'Calculate DET = q^2/4 + p^3/27

CMPN1(1,1)=Q(1)
CMPN1(1,2)=Q(2)
CMPN2(1,1)=Q(1)
CMPN2(1,2)=Q(2)
CMUL
U(1)=CMPN3(1,1)
U(2)=CMPN3(1,2)
U(1)=U(1)/4
U(2)=U(2)/4

CMPN1(1,1)=P(1)
CMPN1(1,2)=P(2)
CMPN2(1,1)=P(1)
CMPN2(1,2)=P(2)
CMUL
V(1)=CMPN3(1,1)
V(2)=CMPN3(1,2)

CMPN1(1,1)=P(1)
CMPN1(1,2)=P(2)
CMPN2(1,1)=V(1)
CMPN2(1,2)=V(2)
CMUL
W(1)=CMPN3(1,1)
W(2)=CMPN3(1,2)
W(1)=W(1)/27
W(2)=W(2)/27

DET(1)=U(1)+W(1)
DET(2)=U(2)+W(2)

Now, I need the square root of DET, and for that, I need the following subroutine that calculates the square root of the complex number stored in CMPN1.

SUB CSQRT()
'Determine the square root of CMPN1
'Save the result in CMPN3.
DIM R AS DOUBLE
R=SQR(CMPN1(1,1)*CMPN1(1,1)+CMPN1(1,2)*CMPN1(1,2))
CMPN3(1,1)=SQR((R+CMPN1(1,1))/2)
CMPN3(1,2)=SQR((R-CMPN1(1,1))/2)
IF CMPN1(1,2)<0 THEN CMPN3(1,2)=-CMPN3(1,2)
END SUB

The result is stored in CMPN3. SQR is the LibreOffice command to take the square root of the value in parentheses following it.

Now, the section of CompCube that calculates the square root of DET is:

'Calculate the square root of DET

CMPN1(1,1)=DET(1)
CMPN2(1,2)=DET(2)
CSQRT
SDET(1)=CMPN3(1,1)
SDET(2)=CMPN3(1,2)

V(1)=-Q(1)/2+SDET(1)
V(2)=-Q(2)/2+SDET(2)

If you look back up at the cubic formula, you will see where those last two statements came from.

So, I need the three cube roots of V and to get that, I need a subroutine to calculate the cube root of CMPN1 and store the results in CMPN2, CMPN3, and CMPN4, and that subroutine will need a function to determine the phase of CMPN1.

Here's the function CPHASE:

FUNCTION CPHASE(NUM,DENOM)
'Returns the phase of a complex number between -pi and pi
DIM PPI AS DOUBLE, CP AS DOUBLE
PPI=4*ATN(1)
IF ABS(DENOM)<1e-15 THEN
IF NUM<0 THEN
CPHASE=-PP1/2
ELSE
CPHASE=PPI/2
END IF
ELSE
CP=ATN(NUM/DENOM)
IF DENOM<0 THEN
CPHASE=CPHASE+PPI
ELSE
CPHASE=CP
END IF
END IF
END FUNCTION

Notice that the header says "FUNCTION" instead of "SUB". This code takes two arguments: NUM and DENOM, and returns a value. The value is determined by the "CPHASE=" statements. When you set a function's name equal to a value, you are outputing the value. To call the function, you set some variable equal to the function name (You'll see that below). ATN is a LibreOffice Basic function that calculates the arctangent of some number. ABS is the function that returns the absolute value of a number. This function is basically a decisiion tree that determines what value CPHASE will take.

If DENOM is less than 1E-15 (or 1 x 10-15), then if NUM is negative CPHASE is -pi/2, else CPHASE is pi/2. If DENOM is more than 1E-15, then if DENOM is negative, CPHASE is the arctangent of NUM/DENOM + pi, otherwise, it is just the arctangent of NUM/DENOM.

Here's the subroutine CCUBRT, which calculates the cube roots of a complex number (there will be three cube roots).

SUB CCUBRT()
'Calculates the three cube roots of CMPN1.
'The results are stored in CMPN2, CMPN3, and CMPN4
DIM PPI AS DOUBLE, R AS DOUBLE, T AS DOUBLE, TT AS DOUBLE

PPI=4*ATN(1)
R=SQR(CMPN1(1,1)*CMPN1(1,1)+CMPN1(1,2)*CMPN1(1,2))
R=R^(1/3)
T=CPHASE(CMPN1(1,2),CMPN1(1,1))
T=T/3
CMPN2(1,1)=R*COS(T)
CMPN2(1,2)=R*SIN(T)
TT=T-(2*PPI/3)
CMPN3(1,1)=R*COS(TT)
CMPN3(1,2)=R*SIN(TT)
TT=T+(2*PPI/3)
CMPN4(1,1)=R*COS(TT)
CMPN4(1,2)=R*SIN(TT)
END SUB

The statement that calls CPHASE is:

T=CPHASE(CMPN1(1,2),CMPN1(1,1))

It feeds the values stored in CMPN1 to the function and transfers the result to the variable T.

The part of CompCube that calculates the cubed roots of V is:

CMPN1(1,1)=V(1)
CMPN1(1,2)=V(2)
CCUBRT
U(1)=CMPN2(1,1)
U(2)=CMPN2(1,2)
U1(1)=CMPN3(1,1)
U1(2)=CMPN3(1,2)
U2(1)=CMPN4(1,1)
U2(2)=CMPN4(1,2)

U, U1, and U2 stores the resulting cubed roots.

And the part that calculates the cube roots of W is almost exactly the same except for the variables used:

W(1)=-Q(1)/2-SDET(1)
W(2)=-Q(2)/2-SDET(2)

CMPN1(1,1)=W(1)
CMPN1(1,2)=W(2)
CCUBRT
V(1)=CMPN2(1,1)
V(2)=CMPN2(1,2)
V1(1)=CMPN3(1,1)
V1(2)=CMPN3(1,2)
V2(1)=CMPN4(1,1)
V2(2)=CMPN4(1,2)

The next section collects all the possible roots of the equation. Only three will be appropriate, but the nine possible roots will all have to be tested.

'Find all possible solutions to the formula.

Z(1,1)=U(1)+V(1)
Z(1,2)=U(2)+V(2)
Z(2,1)=U(1)+V1(1)
Z(2,2)=U(2)+V1(2)
Z(3,1)=U(1)+V2(1)
Z(3,2)=U(2)+V2(2)
Z(4,1)=U1(1)+V(1)
Z(4,2)=U1(2)+V(2)
Z(5,1)=U1(1)+V1(1)
Z(5,2)=U1(2)+V1(2)
Z(6,1)=U1(1)+V2(1)
Z(6,2)=U1(2)+V2(2)
Z(7,1)=U2(1)+V(1)
Z(7,2)=U2(2)+V(2)
Z(8,1)=U2(1)+V1(1)
Z(8,2)=U2(2)+V1(2)
Z(9,1)=U2(1)+V2(1)
Z(9,2)=U2(2)+V2(2)

To test these values, I have to use them to evaluate the cubic equation. I use another subroutine for that, and that subroutine requires a function to return the absolute value of the complex number. I also realized that I needed to save the complex coefficients for this routine so I went back and did that in MatArray1, a global array I use to transfer matrices around.

'Load equation into MatArray1
REDIM MatArray1(4,2)
MatArray1(1,1)=A(1)
MatArray1(1,2)=A(2)
MatArray1(2,1)=B(1)
MatArray1(2,2)=B(2)
MatArray1(3,1)=C(1)
MatArray1(3,2)=C(2)
MatArray1(4,1)=D(1)
MatArray1(4,2)=D(1)

The sizes of dynamic arrays like the global MatArrays can be changed any time using the REDIM (redimension) command.

Here are the test routine and absolute value function:

FUNCTION CABS(ARG1, ARG2)
'This function returns the absolute value of a complex number
'ARG1+ARG2 i
CABS=SQR(ARG1*ARG1+AGR2*ARG2)
END FUNCTION

SUB CTEST()
'This subroutine tests a solution set of a cubic equation for appropriate
'values.
'The output is z1=a z^3 + b x^2 + c z + d
'In other words, it evaluates the original cubic equation.
'It gets the equation from MatArray1.
'The value of the complex variable is in CMPN4.
DIM U(2), V(2), W(2), Z1(2)

CMPN1(1,1)=MatArray1(1,1)
CMPN1(1,2)=MatArray1(1,2)
CMPN1(2,1)=CMPN4(1,1)
CMPN1(2,2)=CMPN4(1,2)
CMUL
U(1)=CMPN3(1,1)
U(2)=CMPN3(1,2)

CMPN1(1,1)=CMPN4(1,1)
CMPN1(1,2)=CMPN4(1,2)
CMPN1(2,1)=U(1)
CMPN1(2,2)=U(2)
CMUL
V(1)=CMPN3(1,1)
V(2)=CMPN3(1,2)

CMPN1(1,1)=CMPN4(1,1)
CMPN1(1,2)=CMPN4(1,2)
CMPN1(2,1)=V(1)
CMPN1(2,2)=V(2)
CMUL
W(1)=CMPN3(1,1)
W(2)=CMPN3(1,2)

Z1(1)=W(1)
Z1(2)=W(2)

CMPN1(1,1)=MatArray1(2,1)
CMPN1(1,2)=MatArray1(2,2)
CMPN1(2,1)=CMPN4(1,1)
CMPN1(2,2)=CMPN4(1,2)
CMUL
U(1)=CMPN3(1,1)
U(2)=CMPN3(1,2)

CMPN1(1,1)=U(1)
CMPN1(1,2)=U(2)
CMPN1(2,1)=CMPN4(1,1)
CMPN1(2,2)=CMPN4(1,2)
CMUL
V(1)=CMPN3(1,1)
V(2)=CMPN3(1,2)

Z1(1)=Z1(1)+V(1)
Z1(2)=Z1(2)+V(2)

CMPN1(1,1)=MatArray1(3,1)
CMPN1(1,2)=MatArray1(3,2)
CMPN1(2,1)=CMPN4(1,1)
CMPN1(2,2)=CMPN4(1,2)
CMUL
U(1)=CMPN3(1,1)
U(2)=CMPN3(1,2)

Z1(1)=Z1(1)+U(1)+MatArray1(4,1)
Z1(2)=Z1(2)+U(2)+MatArray1(4,2)

CMPN3(1)=Z1(1)
CMPN3(2)=Z1(2)

END SUB

After all that I had some small bugs (small bugs = disaster) and had to do some fine tuning. Debugging got rid of the bugs and my test subroutine was too sensitive for the tiny round off errors I produced so I had to change the test criterion:

IF X<0.001 THEN
to the more liberal:

IF X<10 THEN

Now, it works, and it work within input precision.

I'll be programming a lot more and it will be available at the Timeline. The basics of LibreOffice Basic are here and they should give you enough to figure out what I'm doing in future programs. I have a large, complex project coming up that will give me a chance to share some of my problem solving techniques with you.