# Approximating the roots of a cubic equation in one variable through simultaneous geometric structures in the plane of complex numbers

Daniel Alba-Cuellar

► To cite this version:

Daniel Alba-Cuellar. Approximating the roots of a cubic equation in one variable through simultaneous geometric structures in the plane of complex numbers. 2025. hal-05062888

HAL Id: hal-05062888

https://hal.science/hal-05062888v1

Preprint submitted on 10 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Approximating the roots of a cubic equation in one variable through simultaneous geometric structures in the plane of complex numbers

Daniel Alba Cuellar

albacd@cimat.mx

## Abstract

This paper describes how to approximate the roots $r_1, r_2, r_3$ of a cubic polynomial equation in one complex variable of the form $x^3 + c_1x^2 + c_2x + c_3 = 0$, by means of simultaneous geometric structures $L_1L_T C(\theta)$ in the plane of complex numbers ($\theta$ is the quantity that varies among these simultaneous structures). Each structure $L_1L_T C(\theta)$ consists of the following three interrelated geometric elements: **1:** A line $L_1(\theta): p_1 + te^{i\theta}$ ($t \in \mathbb{R}$) with fixed point $p_1 = -c_1/2$ and inclination angle $\theta$; $L_1(\theta)$ seeks to contain one of the roots of the cubic equation, say $r_1$.

**2:** A circumference $C(\theta): -c_3/L_1(\theta)$, which seeks to contain the product of roots $r_2r_3$. **3:** A semi-line $L_T(\theta): c_2 - L_1(\theta)(p_1 - te^{i\theta})$, which also seeks to contain $r_2r_3$. The possible intersections $I_1, I_2$ between semi-line $L_T(\theta)$ and circumference $C(\theta)$ provide information on how close line $L_1(\theta)$ is to root $r_1$; this information is encoded into quantities $-1 < e_A(\theta) \le 1$ and $-1 < e_B(\theta) \le 1$: the closer $e_A(\theta)$ or $e_B(\theta)$ is to zero, the closer $L_1(\theta)$ is to $r_1$. In the end, discrete versions of proximity mappings $\theta \to e_A(\theta)$, $\theta \to e_B(\theta)$, with $\theta \in [-\pi, \pi]$, help us approximate the roots of the cubic equation. This procedure, called cubic LC method, is illustrated through a numerical example, whose results can be reproduced with the help of the two R programs listed in the appendices of this document; these two R programs can also be used to analyze arbitrary cubic equations. This paper concludes by briefly outlining how to extend the cubic LC method to univariate polynomial equations of degree $n \ge 4$.

**Keywords:** Cubic Equation; Roots of an Equation; Complex Numbers; Geometric Structure; Line; Terminal Semi-line; Circumference; Proximity Map; Numerical Method; Geometrical Method; Möbius Transformation; Parallel Computing; Vectorization of Operations; Polynomial Root-finding Algorithm; cubic LC method.

## 1. Description of the cubic LC method

Let's consider the cubic equation

$$x^3 + c_1x^2 + c_2x + c_3 = 0 \tag{1}$$

with roots $r_1, r_2, r_3$ such that (1) can be rewritten in factored form as

$$(x - r_1)(x - r_2)(x - r_3) = 0. \tag{2}$$

If we carry out the multiplications on the left side of equation (2), we obtain the following:

$$ \begin{aligned} [(x - r_1)(x - r_2)](x - r_3) &= [x^2 - (r_1 + r_2)x + r_1r_2](x - r_3) \\ &= x^3 - (r_1 + r_2 + r_3)x^2 + (r_1r_2 + r_1r_3 + r_2r_3)x - r_1r_2r_3 \end{aligned} \tag{3} $$

From equations (1) and (3), we see that the coefficients $c_1, c_2, c_3$ can be expressed in terms of the roots $r_1, r_2, r_3$ in the following way:
$$c_1 = -(r_1 + r_2 + r_3)$$

$$c_2 = r_1r_2 + r_1r_3 + r_2r_3$$

$$c_3 = -r_1r_2r_3.$$

(4)

Expressions (4) are called Vieta's formulas for cubic equations; these relations between the roots and the coefficients of equation (1) are central in extending the quadratic LC method [1] to cubic polynomial equations in one complex variable. Our goal in this paper is to describe and illustrate a basic implementation of such an extension, which is called cubic LC method, and whose purpose is to approximate the unknown roots $r_1, r_2, r_3$ of equation (1) from its known coefficients $c_1, c_2, c_3$. In the following description, we will assume that the roots $r_1, r_2, r_3$ are complex numbers, all different from zero and all different from each other.

To begin with, let's consider the following interrelated geometric elements—depicted in Fig. 1, and based on Vieta's formulas (4)—in the plane of complex numbers $C$:

1. A line $L_1$ that in principle does not pass through the origin. This line has the form $L_1: p_1 + tv_{\theta_1}^*,$ where $p_1 := -c_1/2$ is the fixed point of $L_1, t \in \mathbb{R}$ is a parameter used to locate each point on $L_1$, and $v_{\theta_1}^* := e^{i\theta_1^*} = \cos\theta_1^* + i \sin\theta_1^*$ is a unit direction vector for $L_1$, so that $\theta_1^*$ is the inclination angle of $L_1$ that causes $L_1$ to contain one of the roots of equation (1), say $r_1$; notice that the translation of $v_{\theta_1}^*$, when anchored to fixed point $p_1$ (instead of being anchored to the origin), points to root $r_1$. Under these circumstances, $L_1$ also contains the sum of roots $r_2 + r_3$; furthermore, $|r_1 - p_1| = |(r_2 + r_3) - p_1| > 0$; i.e., $p_1$ is the midpoint of rectilinear segment $\overline{r_1(r_2 + r_3)} \subset L_1$.

2. A circumference $C: -c_3/L_1$. Notice that $C$ contains the fixed point $-c_3/p_1$, and the product of roots $r_2r_3$. In [1], it is shown that a complex constant (in this case $-c_3$), divided by a line $L$ that does not pass through the origin, is a circumference that passes through the origin, without containing it.

3. A terminal semi-line $L_T: c_2 - (p_1 + tv_{\theta_1^*})(p_1 - tv_{\theta_1^*})$, derived from line $L_1$, which also contains the product of roots $r_2r_3$.

These geometric elements form a geometric structure $L_1L_T C(c_1, c_2, c_3, \theta_1^*)$ which depends on the coefficients $c_1, c_2, c_3$ of equation (1), and on the inclination angle $\theta_1^*$ of line $L_1$. In this structure $L_1L_T C(c_1, c_2, c_3, \theta_1^*)$, we have one intersection $I$ between $L_T$ and $C$ which coincides with the product $r_2r_3$; this intersection $I$ can be used to find the location of $r_1$ inside line $L_1$, by simply computing

$$r_1 = -c_3/I$$

(5)

(here, we are projecting one point in circumference $C$ onto line $L_1$, through a Möbius transformation). In [1], a general algorithm to compute intersections between an arbitrary line $L$ and an arbitrary circumference $C$ on the complex plane $\mathbb{C}$ is described; this algorithm of course can be adapted to the situation depicted in Fig. 1, in order to compute $I$. Now, to construct $L_T$ from $L_1$, we start by constructing derived semi-line $L_d: (p_1 + tv_{\theta_1^*})(p_1 - tv_{\theta_1^*}) = (p_1 + te^{i\theta_1^*})(p_1 - te^{i\theta_1^*}) = p_1^2 - t^2e^{i2\theta_1^*} = p_1^2 + t^2e^{i(2\theta_1^*+\pi)}$, which has fixed point $p_1^2$ and unit direction vector $e^{i(2\theta_1^*+\pi)} = -v_{\theta_1^*}^2$, with inclination angle $2\theta_1^* + \pi$. Notice that $L_d$ contains the point $r_1r_2 + r_1r_3$. From here, we can define $L_T: c_2 - L_d$. Notice that $L_T: c_2 - [p_1^2 + t^2e^{i(2\theta_1^*+\pi)}] = (c_2 - p_1^2) - t^2e^{i(2\theta_1^*+\pi)} = (c_2 - p_1^2) + t^2e^{i2\theta_1^*}$; thus, $L_T$ is a
semi-line with fixed point $(c_2 - p_1^2)$ and unit direction vector $e^{i2\theta_1^*} = v_{\theta_1^*}^2$, so $L_T$ has an inclination angle equal to $2\theta_1^*$, as depicted in Fig. 1. Notice also that $L_T$ contains the point $c_2 - (r_1r_2 + r_1r_3)$; from Vieta's formulas (4), we see that $c_2 - (r_1r_2 + r_1r_3) = r_2r_3$.

**NOTE:** Another way to find the location of root $r_1$ inside line $L_1$ is by equating intersection $I$ to terminal semi-line $L_T$:

$$ 
\begin{aligned}
I &= L_T \\
I &= (c_2 - p_1^2) + t^2 v_{\theta_1^*}^2 \\
t^2 &= [I - (c_2 - p_1^2)]/v_{\theta_1^*}^2
\end{aligned} 
\qquad (6) $$

The numerator and denominator at the right side of equation (6) are both vectors with the same argument $2\theta_1^*$, so

$$ 
\begin{aligned}
t^2 &= |I - (c_2 - p_1^2)| \\
t &= \pm\sqrt{|I - (c_2 - p_1^2)|}.
\end{aligned} 
\qquad (7) $$

Parameter $t$ in (7) allows us to locate $r_1$ inside line $L_1$:

$$ r_1 = p_1 + \sqrt{|I - (c_2 - p_1^2)|} v_{\theta_1^*}. \qquad (8) $$

Fig. 1 schematically shows this geometric structure $L_1L_T C(c_1, c_2, c_3, \theta_1^*)$.

**Geometric structures** $L_1L_T C(c_1, c_2, c_3, \theta)$. Unlike the quadratic LC method described in [1], this time there is not a direct way to deduce the inclination angle $\theta_1^*$ of $L_1$ solely from the polynomial coefficients of equation (1): terminal semi-line $L_T$ does not have two points that can be directly constructed from the coefficients $c_1, c_2, c_3$ of equation (1); instead, it has only one constructible point (fixed point $c_2 - p_1^2$). The same can be said of line $L_1$, which has only one known point (fixed point $p_1$). As for circumference $C: -c_3/L_1$, it has only two constructible points, $-c_3/p_1$ and 0. We need to consider a more general structure $L_1L_T C(c_1, c_2, c_3, \theta)$, with inclination angle $\theta \in [-\pi, \pi)$ for $L_1$ arbitrarily close to $\theta_1^*$, in order to get additional information that allows us to approximate the numerical value of theta root $\theta_1^*$. Fig. 2 schematically shows such geometric structure $L_1L_T C(c_1, c_2, c_3, \theta)$.

From Figs. 1 and 2, we see that angle $\theta$ of Fig. 2 is somewhat different from angle $\theta_1^*$ of Fig. 1, so line $L_1$ is close to root $r_1$ in Fig. 2, while line $L_1$ contains root $r_1$ in Fig. 1.

Geometric elements $L_1(c_1, \theta): -c_1/2 + te^{i\theta}$, $C(c_1, c_3, \theta): -c_3/L_1(c_1, \theta)$, and $L_T(c_1, c_2, \theta): (c_2 - c_1^2/4) + t^2e^{i2\theta}$ from Fig. 2 are directly obtainable from known values $c_1, c_2, c_3, \theta$. The intersection $I'$ between terminal semi-line $L_T(c_1, c_2, \theta)$ and circumference $C(c_1, c_3, \theta)$ can be obtained by using an algorithm similar to the one described in [1] (in this paper we will adapt that algorithm to the case of cubic polynomial equations). Now, after computing intersection $I'$, we can project it onto line $L_1(c_1, \theta)$ to obtain points $c_0, i_0, i_1$. Projection $c_0$ is based on equation (5):

$$ c_0 = -c_3/I'. \qquad (9) $$
Expression (9) tells us that $c_0 \in L_1(c_1, \theta)$ is obtained through circumference $C(c_1, c_3, \theta)$. Similarly, projections $i_0, i_1 \in L_1(c_1, \theta)$ are obtained through semi-line $L_T(c_1, c_2, \theta)$, and are based on equation (8):

$$i_0 = -c_1/2 + \sqrt{|I' - (c_2 - c_1^2/4)|}v_\theta \quad (10)$$

$$i_1 = -c_1/2 - \sqrt{|I' - (c_2 - c_1^2/4)|}v_\theta, \quad (11)$$

where $v_\theta := e^{i\theta} = \cos\theta + i\sin\theta$.

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Geometric structure $L_1L_T C(c_1, c_2, c_3, \theta_1^*)$, which contains root $r_1$ at line $L_1: p_1 + tv_{\theta_1^*}$, and product of roots $r_2r_3$ at the intersection $I$ between circumference $C: -c_3/L_1$ and terminal semi-line $L_T: c_2 - L_d$, where $L_d: (p_1 + tv_{\theta_1^*})(p_1 - tv_{\theta_1^*})$ is a semi-line directly derived from $L_1$ ($L_d$ is not shown in this figure).</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Geometric structure $L_1L_T C(c_1, c_2, c_3, \theta)$, with the same coefficients $c_1, c_2, c_3$ used in Fig. 1, and inclination angle $\theta$ for line $L_1$ close to angle $\theta_1^*$ from Fig. 1.</figcaption></figure>

To prove the correctness of expressions (10) and (11), let's verify from them that $c_2 - i_0 i_1$ is equal to $I'$:

$$ 
\begin{aligned}
c_2 - i_0 i_1 &= c_2 - \left(-\frac{c_1}{2} + \sqrt{|I' - (c_2 - c_1^2/4)|} v_\theta\right) \left(-\frac{c_1}{2} - \sqrt{|I' - (c_2 - c_1^2/4)|} v_\theta\right) \\
&= c_2 - \left(c_1^2/4 - \left[\sqrt{|I' - (c_2 - c_1^2/4)|} v_\theta\right]^2\right) \\
&= (c_2 - c_1^2/4) + |I' - (c_2 - c_1^2/4)| v_\theta
\end{aligned} 
\quad (12) $$
$$ 
\begin{aligned}
&= (c_2 - c_1^2/4) + [I' - (c_2 - c_1^2/4)] \\
&= I'.
\end{aligned} 
\qquad (13) $$

Expressions (12) and (13) tell us that $|I' - (c_2 - c_1^2/4)|v_\theta^2 = I' - (c_2 - c_1^2/4)$; this is true because $I' - (c_2 - c_1^2/4)$ can be viewed as a linear segment in $L_T(c_1, c_2, \theta)$ with length $|I' - (c_2 - c_1^2/4)|$ and unit direction vector $v_\theta^2$ (see Fig. 2).

Points $c_0$ in (9) and $i_0$ in (10) can both be viewed as approximations to $r_1$ if we recognize their similarities with expressions (5) and (8); their distance on $L_1(c_1, \theta)$ can be used to measure the proximity between line $L_1(c_1, \theta)$ and root $r_1$. Let's look at the following definition:

$$ e(c_1, c_2, c_3, \theta) := \operatorname{sgn} \left( \frac{i_0 - c_0}{v_\theta} \right) \frac{|i_0 - c_0|}{|i_0 - p_1| + |p_1 - c_0|'} \qquad (14) $$

$$ \text{with } \operatorname{sgn}(x) = \begin{cases} 1 & \text{if } x > 0 \\ 0 & \text{if } x = 0 \\ -1 & \text{if } x < 0 \end{cases} \text{ being the sign function.} $$

In expression (14), the fraction $\frac{i_0 - c_0}{v_\theta}$ inside the sign function is a real number, because complex numbers $i_0 - c_0$ and $v_\theta$, when viewed as vectors in $\mathbb{R}^2$, both lie on the same line $L_1$; that is to say, they are collinear, so one vector is a scalar multiple of the other.

Expression (14) is called the **weighted error** associated with geometric structure $L_1L_T C(c_1, c_2, c_3, \theta)$ from Fig. 2.

By definition, $e(c_1, c_2, c_3, \theta)$ is a real quantity that is always between -1 and 1 if there exists the intersection $I'$. To see this, refer to Fig. 2 and try to visualize the following scenarios on the line $L_1$:

a) If $c_0$ is on the same side as $i_1$ with respect to fixed point $p_1$, or if $c_0 = p_1$, then $e(c_1, c_2, c_3, \theta) = 1$.

b) If $c_0$ is on the same side as $i_0$ with respect to fixed point $p_1$, but $|c_0 - p_1| < |i_0 - p_1|$, then $0 < e(c_1, c_2, c_3, \theta) < 1$.

c) If $c_0$ is on the same side as $i_0$ with respect to fixed point $p_1$, but $|c_0 - p_1| > |i_0 - p_1|$, then $-1 < e(c_1, c_2, c_3, \theta) < 0$.

d) If $c_0 = i_0$, then we have the same circumstances as those shown in Fig. 1, and $e(c_1, c_2, c_3, \theta) = 0$, so $c_0 = i_0 = r_1$ and $\theta = \theta_1^*$.

Schematic Figs. 1 & 2, in which there is only one intersection between terminal semi-line $L_T$ and circumference $C$, are useful for visualizing the concepts and ideas we have just developed and summarized into definition (14). In some other situations, however, an arbitrary geometric structure $L_1L_T C(c_1, c_2, c_3, \theta)$ can have two intersections, or no intersections between terminal semi-line $L_T(c_1, c_2, \theta)$ and circumference $C(c_1, c_3, \theta)$. Fig. 3 shows some of the possible scenarios that arise when trying to compute intersections between a terminal semi-line $L_T$ and a circumference $C$.

As we can see from Figs. 3 (a) & 3 (c), it is possible to have a non-tangential intersection or a tangential intersection between $L_T(c_1, c_2, \theta)$ and $C(c_1, c_3, \theta)$; it is also possible to have two different intersections $I_1, I_2$, as in Fig. 3 (b), or no intersections at all, as in Fig. 3 (d).
<figure><img src="image_4.png" /><figcaption>Figure 3. Some possible scenarios for the intersections between terminal semi-line $L_T$ and circumference $C$.</figcaption></figure>

An arbitrary geometric structure $L_1L_T C(c_1, c_2, c_3, \theta)$ can be useful to quantify the proximity of line $L_1(c_1, \theta)$ to one of the roots of equation (1), if we keep track of the weighted errors associated with intersections $I_1, I_2$ between terminal semi-line $L_T(c_1, c_2, \theta)$ and circumference $C(c_1, c_3, \theta)$. This means we must extend our definition (14):

$$ e_A(c_1, c_2, c_3, \theta) := \begin{cases} \operatorname{sgn}\left(\frac{i_{0I_1}-c_{0I_1}}{v_\theta}\right) \frac{|i_{0I_1}-c_{0I_1}|}{|i_{0I_1}-p_1|+|p_1-c_{0I_1}|} & \text{if there exists } I_1 \\ \operatorname{undefined} & \text{if there does not exist } I_1 \end{cases} \tag{15} $$

$$ e_B(c_1, c_2, c_3, \theta) := \begin{cases} \operatorname{sgn}\left(\frac{i_{0I_2}-c_{0I_2}}{v_\theta}\right) \frac{|i_{0I_2}-c_{0I_2}|}{|i_{0I_2}-p_1|+|p_1-c_{0I_2}|} & \text{if there exists } I_2 \\ \operatorname{undefined} & \text{if there does not exist } I_2 \end{cases} $$

In definition (15), $c_{0I_1}$ and $i_{0I_1}$ are the projections of intersection $I_1$ onto line $L_1(c_1, \theta)$; likewise, $c_{0I_2}$ and $i_{0I_2}$ are the projections of intersection $I_2$ onto line $L_1(c_1, \theta)$. These projections $c_{0I_1}, i_{0I_1}, c_{0I_2}, i_{0I_2}$ are computed by means of expressions (9) and (10), using $I_1$ or $I_2$ instead of $I'$. At this point, we need to establish a criterion to decide which intersection is $I_1$ and which is $I_2$. In [1], equations (12) through (15), it is shown that the intersections $I_1, I_2$ between an arbitrary line $L: p + \tau v_\theta (p \in \mathbb{C})$ is any point in $L$,
$v_\theta := e^{i\theta}$ is L's unit direction vector, and $\tau \in \mathbb{R}$ is a parameter for locating points in L) and an arbitrary circumference C with center $c \in \mathbb{C}$ and radius $r \in \mathbb{R}$, can be computed as $I_1 = p + \tau_1 v_\theta$, $I_2 = p + \tau_2 v_\theta$, where

$$ \tau_1 = [(c - p) \cdot v_\theta] - \sqrt{\Delta} \qquad (16) $$

$$ \tau_2 = [(c - p) \cdot v_\theta] + \sqrt{\Delta} $$

with

$$ \Delta = [(c - p) \cdot v_\theta]^2 - (c - p) \cdot (c - p) + r^2. \qquad (17) $$

In expressions (16) and (17), $(c - p)$ and $v_\theta$ can be viewed as vectors in $\mathbb{R}^2$, and the operator $\cdot$ is the dot product between two vectors in $\mathbb{R}^2$. From here, we can use expressions (16) and (17) to compute the intersections $I_1, I_2$ between a terminal semi-line $L_T(c_1, c_2, \theta): (c_2 - p_1^2) + t^2 e^{i2\theta}$ and a circumference $C(c_1, c_3, \theta): -c_3/(p_1 + te^{i\theta})$ by setting

$$ p = c_2 - p_1^2, \qquad (18) $$

$$ v_{\theta} = e^{i2\theta}, \qquad (19) $$

$$ c = \frac{-y_1(x_2^2+y_2^2)+y_2(x_1^2+y_1^2)}{2x_1y_2-2x_2y_1} + i \frac{x_1(x_2^2+y_2^2)-x_2(x_1^2+y_1^2)}{2x_1y_2-2x_2y_1} \quad (\text{see expression (8) in [1]}), \qquad (20) $$

where $x_1 = \text{Re}(-c_3/p_1)$, $y_1 = \text{Im}(-c_3/p_1)$, $x_2 = \text{Re}[-c_3/(p_1 + e^{i\theta})]$, $y_2 = \text{Im}[-c_3/(p_1 + e^{i\theta})]$, and

$$ r = |c|. \qquad (21) $$

Thus, we define the intersections $I_1, I_2$ between terminal semi-line $L_T(c_1, c_2, \theta): (c_2 - p_1^2) + t^2 e^{i2\theta}$ and circumference $C(c_1, c_3, \theta): -c_3/(p_1 + te^{i\theta})$ as

$$ I_1 := (c_2 - p_1^2) + \tau_1 e^{i2\theta} \qquad (22) $$

$$ I_2 := (c_2 - p_1^2) + \tau_2 e^{i2\theta}, $$

where $\tau_1, \tau_2$ are the parametric values computed by means of expressions (16) through (21). For the intersection $I_1$ to exist, it is necessary that $\Delta \ge 0$ and $\tau_1 \ge 0$. Similarly, for the intersection $I_2$ to exist, it is necessary that $\Delta \ge 0$ and $\tau_2 \ge 0$. From expressions (16), we can see that if $\Delta > 0$, then we always have $\tau_1 < \tau_2$. Now, notice that, if for example, $\tau_1 < 0$ and $\tau_2 > 0$, then there are two intersections $I_1, I_2$ between the complete line associated with $L_T$ and the circumference C, but only one intersection, $I_2$, is "relevant" in the context of intersections between terminal semi-line $L_T(c_1, c_2, \theta)$ and circumference $C(c_1, c_3, \theta)$, since the direction vector of $L_T(c_1, c_2, \theta)$, when anchored to $(c_2 - p_1^2)$, points to $I_2$, but not to $I_1$. This is illustrated in Fig. 3 (a), where $I = I_2$, according to criterion (22). In general, we say that an intersection $I_m$ ($m = 1,2$) between terminal semi-line $L_T(c_1, c_2, \theta)$ and circumference $C(c_1, c_3, \theta)$, defined by expressions (22), is **relevant** if $\tau_m \ge 0$.

**A note on computational implementation:** In appendix 1 of this paper, R function `intersect_semiLine_circle` computes possible intersections $I_1, I_2$ between arbitrary semi-lines and circumferences, using general expressions (16) and (17), and incorporating the concept of relevant intersections between semi-lines and circumferences. Also in appendix 1, R function `L1LtC` builds geometric structures $L_1L_T C(c_1, c_2, c_3, \theta)$, then calls `intersect_semiLine_circle` using the parameters
given by expressions (18) through (21) to obtain intersections $I_1, I_2$ given by (22), and finally computes weighted errors $e_A(\theta), e_B(\theta)$ given by (15). These two R functions are central in generating smooth and continuous proximity maps $\theta \to e_A(\theta), \theta \to e_B(\theta)$, with $\theta \in [-\pi, \pi]$. We can use $e_A(\theta), e_B(\theta)$ as a shorthand notation for $e_A(c_1, c_2, c_3, \theta), e_B(c_1, c_2, c_3, \theta)$, as long as the coefficients of equation (1) remain unchanged.

**Summarizing up to this point**, our cubic LC method allows us—through a geometric structure $L_1L_T C(c_1, c_2, c_3, \theta)$ with associated weighted errors $e_A(\theta), e_B(\theta)$ defined by (15)—to know how close a line $L_1(c_1, \theta): -c_1/2 + te^{i\theta}$ is to one of the roots of equation (1): select any angle $\theta \in [-\pi, \pi)$, and the corresponding geometric structure $L_1L_T C(c_1, c_2, c_3, \theta)$ will produce up to two proximity measures $-1 < e_A(\theta) \le 1, -1 < e_B(\theta) \le 1$; the closer $e_A(\theta)$ or $e_B(\theta)$ is to 0, the closer $L_1(c_1, \theta)$ is to a root.

**Parallel strategy to generate proximity maps.** Now imagine that $N$ processing units $u_0, u_1, ..., u_k, ...$, $u_{N-1}$ with parallel computing capabilities are assigned the task of computing the $N$ pairs of weighted errors $\langle e_A(\theta_0), e_B(\theta_0) \rangle, \langle e_A(\theta_1), e_B(\theta_1) \rangle, ..., \langle e_A(\theta_k), e_B(\theta_k) \rangle, ..., \langle e_A(\theta_{N-1}), e_B(\theta_{N-1}) \rangle$, where $\theta_k = -\pi + 2\pi k/N, k = 0, 1, ..., N-1$. A time-efficient way to accomplish this task is to use—for each $k$ at the same time—processing unit $u_k$ to build geometric structure $L_1L_T C(c_1, c_2, c_3, \theta_k)$ and compute weighted errors $\langle e_A(\theta_k), e_B(\theta_k) \rangle$. In this way we can quickly obtain a pair of discrete proximity maps $\theta_k \to e_A(\theta_k), \theta_k \to e_B(\theta_k)$, which will help us approximate (estimate) theta roots $\theta_1^*, \theta_2^*, \theta_3^*$ via linear interpolation. From estimated theta roots $\hat{\theta}_1^*, \hat{\theta}_2^*, \hat{\theta}_3^*$, we can also approximate the roots $r_1, r_2, r_3$ of equation (1) via geometric structures $L_1L_T C(c_1, c_2, c_3, \hat{\theta}_1^*), L_1L_T C(c_1, c_2, c_3, \hat{\theta}_2^*), L_1L_T C(c_1, c_2, c_3, \hat{\theta}_3^*)$. It is to be expected that, to obtain good approximations, this parallel strategy should use large values of $N$ on the order of thousands, tens of thousands or even more.

The feasibility of this parallel strategy can be verified by using the R statistical computing software installed on a typical personal computer, without the need for thousands of physical parallel processing cores. This is possible because R uses vectorization of operations; a vectorized operation is an operation applied to complete arrays of operands, like numerical vectors or matrices. This contrasts with normal operations, which are applied to individual operands. Vectorization of operations can be thought of as a basic and accessible form of parallel computing. Appendix 1 lists a complete R program that implements a basic version of the cubic LC method; each of its user-defined functions uses vectorization of operations, including functions L1LtC and intersect_semiLine_circle mentioned above. The R program from appendix 1 is used to generate most of the results shown in the numerical example from section 3 below. Notice that the R program from appendix 1 neither implements a general-purpose version of the cubic LC method, nor uses advanced parallel computing techniques; it only helps us verify and illustrate the concepts and ideas described in this paper. In other words, the R program from appendix 1 is a **proof of concept** for the cubic LC method.

**A few words about linear interpolation:** Now suppose we have two consecutive discrete values $\theta_k, \theta_{k+1}$ such that $e(\theta_k)e(\theta_{k+1}) < 0$; i.e., $e(\theta_k)$ and $e(\theta_{k+1})$ are of opposite signs, as shown in Fig. 4 (here, function $e(\theta)$ can refer to any of the functions $e_A(\theta)$ or $e_B(\theta)$ associated with intersections $I_1, I_2$). The linear segment joining points $(\theta_k, e(\theta_k)), (\theta_{k+1}, e(\theta_{k+1}))$ in $\mathbb{R}^2$ crosses the horizontal axis $y = 0$ at $(\hat{\theta}^*, 0)$. Abscissa $\hat{\theta}^*$ can be viewed as an approximation to any theta root $\theta^*$. To compute $\hat{\theta}^*$ in terms of coordinates $(\theta_k, e(\theta_k)), (\theta_{k+1}, e(\theta_{k+1}))$, we can use the point-slope formula for a line:

$$y - e(\theta_k) = \left[ \frac{e(\theta_{k+1}) - e(\theta_k)}{\theta_{k+1} - \theta_k} \right] (x - \theta_k). \qquad (23)$$
If we plug values $y = 0$ and $x = \hat{\theta}^*$ into equation (23), we obtain

$$ -e(\theta_k) = \left[ \frac{e(\theta_{k+1}) - e(\theta_k)}{\theta_{k+1} - \theta_k} \right] (\hat{\theta}^* - \theta_k). \quad (24) $$

Now, we solve for $\hat{\theta}^*$ in equation (24):

$$ \hat{\theta}^* = \theta_k + e(\theta_k) \left[ \frac{\theta_{k+1} - \theta_k}{e(\theta_k) - e(\theta_{k+1})} \right]. \quad (25) $$

<figure><img src="image_5.png" /><figcaption>**Figure 4.** A linear segment joining points $(\theta_k, e(\theta_k))$, $(\theta_{k+1}, e(\theta_{k+1}))$ in $\mathbb{R}^2$. Since $e(\theta_k)$ and $e(\theta_{k+1})$ are of opposite signs, this linear segment crosses axis $y = 0$ at point $(\hat{\theta}^*, 0)$. Abscissa $\hat{\theta}^*$ can be viewed as an approximation to any theta root $\theta^*$.</figcaption></figure>

In appendix 1, using vectorization of operations, R function `approximate_theta_root` computes values $\hat{\theta}^*$ by applying formula (25) to input arrays $\{\theta_k\}$, $\{e(\theta_k)\}$, $k = 0,1, ..., N-1$, whenever it finds consecutive elements $\theta_k$, $\theta_{k'+1}$ such that $e(\theta_k)e(\theta_{k'+1}) \le 0$. When approximating theta roots, it is convenient to assume that absolute ordinate differences $|e(\theta_{k'+1}) - e(\theta_k)|$ are not too large relative to abscissa differences $\theta_{k'+1} - \theta_k$ (i.e., function $e(\theta)$ gradually and smoothly crosses horizontal axis $y = 0$ in the vicinity of any theta root $\theta^*$); as a rudimentary attempt to incorporate this assumption, `approximate_theta_root` considers only absolute ordinate differences $|e(\theta_{k'+1}) - e(\theta_k)|$ no greater than an input tolerance parameter TOL in the computation of values $\hat{\theta}^*$.

**How to approximate the roots of equation (1) from theta roots approximations.** To conclude the description of the cubic LC method, let's see how to approximate a root, say $r_1$, from a theta root approximation, say $\hat{\theta}_1^*$. For this task, the first things to do are to build geometric structure $L_1L_T C(c_1, c_2, c_3, \hat{\theta}_1^*)$, compute intersection(s) between $L_T(c_1, c_2, \hat{\theta}_1^*)$ and $C(c_1, c_3, \hat{\theta}_1^*)$, and obtain projections $c_0, i_0 \in L_1(c_1, \hat{\theta}_1^*)$. Now, since projections $c_0, i_0$ defined in (9) and (10) can themselves be viewed as individual approximations to $r_1$, we can simply opt to compute the average of $c_0$ and $i_0$; in this way, we have

$$ \hat{r}_1 := [c_{0I'}(\hat{\theta}_1^*) + i_{0I'}(\hat{\theta}_1^*)] / 2, \quad (26) $$
where $c_{0I'}(\hat{\theta}_1^*)$ and $i_{0I'}(\hat{\theta}_1^*)$ are the projections of the intersection $I'$ between semi-line $L_T(c_1, c_2, \hat{\theta}_1^*)$ and circumference $C(c_1, c_3, \hat{\theta}_1^*)$ onto line $L_1(c_1, \hat{\theta}_1^*)$. **NOTE:** When there is only one non-tangential intersection $I'$ between semi-line $L_T$ and circumference $C$, as in Fig. 2, then $I' = I_2$, according to the concept of relevant intersections related to definitions (22).

Observe that the average $\hat{r}_1$ given by (26) is also a point in $L_1(c_1, \hat{\theta}_1^*)$. Now, in some situations, it is also possible to obtain two such averages, namely, when there are two relevant intersections $I_1, I_2$ between $L_T(c_1, c_2, \hat{\theta}_1^*)$ and $C(c_1, c_3, \hat{\theta}_1^*)$; in those cases, it is best to select the average associated with the smallest absolute weighted error. In this way, adapting expression (26) to the case of two intersections $I_1(\hat{\theta}_1^*)$, $I_2(\hat{\theta}_1^*)$ between $L_T(c_1, c_2, \hat{\theta}_1^*)$ and $C(c_1, c_3, \hat{\theta}_1^*)$, we have

$$ \hat{r}_1 := \begin{cases} \left[ c_{0I_1}(\hat{\theta}_1^*) + i_{0I_1}(\hat{\theta}_1^*) \right] / 2 & \text{if } |e_A(\hat{\theta}_1^*)| < |e_B(\hat{\theta}_1^*)| \\ \left[ c_{0I_2}(\hat{\theta}_1^*) + i_{0I_2}(\hat{\theta}_1^*) \right] / 2 & \text{if } |e_A(\hat{\theta}_1^*)| \ge |e_B(\hat{\theta}_1^*)| \end{cases} \quad (27) $$

R function approxP3, listed in appendix 1, coordinates all other R functions listed in appendix 1, and then generates roots approximations $\hat{r}_1, \hat{r}_2, \hat{r}_3$ for equation (1) by using definitions (26) and (27).

## 2. Algorithm for the cubic LC method

This section lists an algorithm for the cubic LC method, based on the concepts and ideas described in section 1. Before listing the algorithm, next there are some notes on how the R program in appendix 1 implements the cubic LC method, by grouping the steps of the algorithm into R functions:

* **Steps a), b), c), d)** are implemented with the help of R function L1LtC, which is called within R function approxP3. L1LtC receives as inputs the coefficients $\{c_1, c_2, c_3\}$ from equation (1), and an array $\{\theta_j\}$ of inclination angles for lines $L_1$. With these inputs, L1LtC simultaneously builds geometric structures $L_1L_T C(c_1, c_2, c_3, \theta_j)$, and outputs their associated discrete proximity maps $\theta_j \to e_A(\theta_j)$, $\theta_j \to e_B(\theta_j)$, together with projected intersections $c_0, i_0$ used in the computation of weighted errors $e_A(\theta_j), e_B(\theta_j)$. For this, L1LtC uses three R functions: center_zcircle (finds centers of circumferences $C$), intersect_semiLine_circle (finds intersections between semi-lines $L_T$ and circumferences $C$), and weighted_error (computes weighted errors $e(\theta)$).

* **Steps e), f)** are implemented with the help of R function approximate_theta_rootS, which is called within R function approxP3. approximate_theta_rootS receives as input a discrete proximity map $\theta_j \to e(\theta_j)$, and outputs theta roots approximations $\hat{\theta}^*$, via linear interpolation.

* **Steps g), h)** are implemented with the help of approxP3, which acts as a coordinating R function. First, approxP3 calls L1LtC to generate discrete proximity maps $\theta_k \to e_A(\theta_k)$, $\theta_k \to e_B(\theta_k)$ with $\theta_k = -\pi + 2\pi k/N, k = 0,1, ..., N - 1$. Then, approxP3 calls approximate_theta_rootS a couple of times, giving it the discrete proximity maps as inputs, to generate theta roots approximations $\hat{\theta}_l^*$ ($l = 1,2,3$), which serve as inputs for a second call to L1LtC; this second call to L1LtC generates projections $c_0, i_0$ associated with theta roots approximations $\hat{\theta}_l^*$. In the end, approxP3 averages these projections $c_0, i_0$ to generate roots approximations $\hat{r}_l$ for equation (1). approxP3 outputs estimations $\hat{\theta}_l^*, \hat{r}_l$, and discrete mappings $\theta_k \to e_A(\theta_k), \theta_k \to e_B(\theta_k)$.
LC method for approximating the roots $r_1, r_2, r_3 \in \mathbb{C}$ of cubic equation $x^3 + c_1x^2 + c_2x + c_3 = 0$
(this method generally works if $r_1, r_2, r_3$ are all nonzero and different among them)

**Input:** complex coefficients $c_1, c_2, c_3$ from equation $x^3 + c_1x^2 + c_2x + c_3 = 0$

**Algorithm steps:**

a) **Generate simultaneous geometric structures $L_1L_T C(c_1, c_2, c_3, \theta_k)$, with $\theta_k = -\pi + 2\pi k/N$, $k = 0, 1, ..., N-1$.** Notice that interval $[-\pi, \pi)$ is evenly segmented by points $\theta_k$. Structures $L_1L_T C(c_1, c_2, c_3, \theta_k)$ contain the following elements:

- Lines $L_1(\theta_k): p_1 + te^{i\theta_k}$, with $p_1 = -c_1/2, t \in \mathbb{R}$

- Circumferences $C(\theta_k): -c_3/L_1(\theta_k)$

- Terminal semi-lines $L_T(\theta_k): (c_2 - p_1^2) + t^2e^{i2\theta_k}$, with $t \in \mathbb{R}$

b) **Simultaneously, compute relevant intersections $I_1(\theta_k), I_2(\theta_k)$ between terminal semi-lines $L_T(\theta_k)$ and circumferences $C(\theta_k)$.** For this, use the following expressions:

$$I_1(\theta_k) = p + \tau_1(\theta_k)e^{i2\theta_k}$$

$$I_2(\theta_k) = p + \tau_2(\theta_k)e^{i2\theta_k},$$

where

$$p = (c_2 - p_1^2) \text{ is the anchor point of terminal semi-line } L_T(\theta_k),$$

$$\tau_1(\theta_k) = [c(\theta_k) - p] \cdot e^{i2\theta_k} - \sqrt{\Delta(\theta_k)}, \text{ and}$$

$$\tau_2(\theta_k) = [c(\theta_k) - p] \cdot e^{i2\theta_k} + \sqrt{\Delta(\theta_k)},$$

with $\cdot$ being the dot product between two complex numbers viewed as vectors in $\mathbb{R}^2$, $c(\theta_k)$ being the center of circle $C(\theta_k)$, and

$$\Delta(\theta_k) = \left( [c(\theta_k) - p] \cdot e^{i2\theta_k} \right)^2 - [c(\theta_k) - p] \cdot [c(\theta_k) - p] + [r(\theta_k)]^2;$$

$r(\theta_k)$ is the radius of circle $C(\theta_k)$. The center of circle $C(\theta_k)$ is computed as

$$c(\theta_k) = \frac{-y_1(x_{2,k}^2+y_{2,k}^2)+y_{2,k}(x_1^2+y_1^2)}{2x_1y_{2,k}-2x_{2,k}y_1} + i \frac{x_1(x_{2,k}^2+y_{2,k}^2)-x_{2,k}(x_1^2+y_1^2)}{2x_1y_{2,k}-2x_{2,k}y_1}, \text{ where } x_1 = \text{Re}(-c_3/p_1),$$

$$y_1 = \text{Im}(-c_3/p_1), \quad x_{2,k} = \text{Re}[-c_3/(p_1 + e^{i\theta_k})], \quad y_{2,k} = \text{Im}[-c_3/(p_1 + e^{i\theta_k})].$$

The radius of circle $C(\theta_k)$ is computed as $r(\theta_k) = |c(\theta_k)|$.

**NOTE:** If $\Delta(\theta_k) < 0$, there are no intersections between semi-line $L_T(\theta_k)$ and circumference $C(\theta_k)$. Intersection $I_m(\theta_k)$ ($m = 1,2$) is **relevant** only if $\tau_m(\theta_k) \ge 0$.

c) **Simultaneously, compute projections of relevant intersections $I_1(\theta_k), I_2(\theta_k)$ onto lines $L_1(\theta_k)$.** For this, use the following expressions:

$$c_{0I_1}(\theta_k) = -c_3/I_1(\theta_k),$$

$$c_{0I_2}(\theta_k) = -c_3/I_2(\theta_k), \qquad \text{(projections via circumference } C(\theta_k)\text{)}$$

$$i_{0I_1}(\theta_k) = p_1 + \sqrt{|I_1(\theta_k) - (c_2 - p_1^2)|}e^{i\theta_k},$$

$$i_{0I_2}(\theta_k) = p_1 + \sqrt{|I_2(\theta_k) - (c_2 - p_1^2)|}e^{i\theta_k} \quad \text{(projections via terminal semi-line } L_T(\theta_k)\text{)}$$

d) **Simultaneously, compute weighted errors $e_A(\theta_k), e_B(\theta_k)$ associated to geometric structures $L_1L_T C(c_1, c_2, c_3, \theta_k)$.** For this, use the following expressions:

$$e_A(\theta_k) = \text{sgn} \left( \frac{i_{0I_1}(\theta_k) - c_{0I_1}(\theta_k)}{e^{i\theta_k}} \right) \frac{|i_{0I_1}(\theta_k) - c_{0I_1}(\theta_k)|}{|i_{0I_1}(\theta_k) - p_1| + |p_1 - c_{0I_1}(\theta_k)|} \quad \text{if } I_1(\theta_k) \text{ is relevant}$$

$$e_B(\theta_k) = \text{sgn} \left( \frac{i_{0I_2}(\theta_k) - c_{0I_2}(\theta_k)}{e^{i\theta_k}} \right) \frac{|i_{0I_2}(\theta_k) - c_{0I_2}(\theta_k)|}{|i_{0I_2}(\theta_k) - p_1| + |p_1 - c_{0I_2}(\theta_k)|} \quad \text{if } I_2(\theta_k) \text{ is relevant}$$
### Cubic LC method (continued)

e) Obtain continuous functions $\theta \to \hat{e}_A(\theta)$, $\theta \to \hat{e}_B(\theta)$, $\theta \in [-\pi, \pi)$, from discrete proximity maps $\theta_k \to e_A(\theta_k)$, $\theta_k \to e_B(\theta_k)$. For this, join consecutive pairs of points $(\theta_k, e_A(\theta_k))$, $(\theta_{k+1}, e_A(\theta_{k+1}))$ in $\mathbb{R}^2$ with line segments, for all $k$, to obtain function $\theta \to \hat{e}_A(\theta)$. Function $\theta \to \hat{e}_B(\theta)$ is obtained in an analogous way; just replace $e_A$ with $e_B$. **NOTE:** This is a virtual step, based on linear interpolation.

f) Find crossings of continuous functions $\theta \to \hat{e}_A(\theta)$, $\theta \to \hat{e}_B(\theta)$ with the x-axis in $\mathbb{R}^2$. In this way, we obtain numerical approximations $\hat{\theta}_1^*$, $\hat{\theta}_2^*$, $\hat{\theta}_3^*$ to theta roots $\theta_1^*$, $\theta_2^*$, $\theta_3^*$.

To compute such intersections, use formula $\hat{\theta}_l^* = \theta_{k'_l} + \hat{e}(\theta_{k'_l}) \left[ \frac{\theta_{k'_l+1} - \theta_{k'_l}}{\hat{e}(\theta_{k'_l}) - \hat{e}(\theta_{k'_l+1})} \right]$, $l = 1, 2, 3$. $k'_l \in \{0, 1, ..., N-1\}$ are indices such that $\hat{e}(\theta_{k'_l}) \hat{e}(\theta_{k'_l+1}) \le 0$ and $|\hat{e}(\theta_{k'_l+1}) - \hat{e}(\theta_{k'_l})| \le 1$. Function $\hat{e}(\cdot)$ can be either $\hat{e}_A(\cdot)$ or $\hat{e}_B(\cdot)$.

g) Generate geometric structures $L_1L_T C(c_1, c_2, c_3, \hat{\theta}_l^*)$, $l = 1, 2, 3$. For these structures, compute relevant intersections $I_1(\hat{\theta}_l^*)$, $I_2(\hat{\theta}_l^*)$, projections $c_{0I_1(\hat{\theta}_l^*)}$, $c_{0I_2(\hat{\theta}_l^*)}$, $i_{0I_1(\hat{\theta}_l^*)}$, $i_{0I_2(\hat{\theta}_l^*)}$, and weighted errors $e_A(\hat{\theta}_l^*)$, $e_B(\hat{\theta}_l^*)$, in an analogous way to steps a), b), c) and d); instead of using angles $\theta_k$, this time use angles $\hat{\theta}_l^*$.

h) Compute numerical approximations $\hat{r}_1, \hat{r}_2, \hat{r}_3$ for the roots of $x^3 + c_1x^2 + c_2x + c_3 = 0$ by means of projection averages. For this, consider the following situations:

- If geometric structure $L_1L_T C(c_1, c_2, c_3, \hat{\theta}_l^*)$ contains two relevant intersections $I_1(\hat{\theta}_l^*)$, $I_2(\hat{\theta}_l^*)$ between $L_T(\hat{\theta}_l^*)$ and $C(\hat{\theta}_l^*)$, then compute

$$ \hat{r}_l = \left\{ \begin{array}{ll} \frac{[c_{0I_1}(\hat{\theta}_l^*) + i_{0I_1}(\hat{\theta}_l^*)]}{2} & \text{if } |e_A(\hat{\theta}_l^*)| < |e_B(\hat{\theta}_l^*)| \\ \frac{[c_{0I_2}(\hat{\theta}_l^*) + i_{0I_2}(\hat{\theta}_l^*)]}{2} & \text{if } |e_A(\hat{\theta}_l^*)| \ge |e_B(\hat{\theta}_l^*)| \end{array} \right. $$

- If geometric structure $L_1L_T C(c_1, c_2, c_3, \hat{\theta}_l^*)$ contains only one relevant intersection $I_2(\hat{\theta}_l^*)$ between $L_T(\hat{\theta}_l^*)$ and $C(\hat{\theta}_l^*)$, then compute

$$ \hat{r}_l = \frac{[c_{0I_2}(\hat{\theta}_l^*) + i_{0I_2}(\hat{\theta}_l^*)]}{2} $$

**Output:** numerical approximations $\hat{r}_1, \hat{r}_2, \hat{r}_3$ (approximations to roots) and $\hat{\theta}_1^*, \hat{\theta}_2^*, \hat{\theta}_3^*$ (approximations to theta roots); discrete proximity maps $\theta_k \to e_A(\theta_k)$, $\theta_k \to e_B(\theta_k)$ (graphical output)

### 3. Numerical example

Approximate the roots of equation

$$ x^3 + (2 - i)x^2 + (1 - 5i)x + (2 - 6i) = 0, \qquad (28) $$

by using the cubic LC method described above.
## Solution

According to the form of equation (1), the coefficients in equation (28) are $c_1 = 2 - i$, $c_2 = 1 - 5i$, $c_3 = 2 - 6i$. For comparison purposes, the true roots $r_l$ and theta roots $\theta_l^*$ ($l = 1,2,3$) of equation (28) are:

$$r_1 = -1 - i \qquad \theta_1^* = -1.5707963$$

$$r_2 = -2 \qquad \theta_2^* = -2.6779450 \qquad (29)$$

$$r_3 = 1 + 2i \qquad \theta_3^* = 0.6435011$$

Remember that a theta root $\theta_l^*$ is the argument of vector $v_{\theta_l^*} = e^{i\theta_l^*}$ with initial point $p_1 = -c_1/2$ that points to root $r_l$ (see Fig. 1). In this case, $p_1 = -1 + 0.5i$.

If we run the R program from appendix 1 (implementation of the algorithm for the cubic LC method listed in section 2) with the coefficients $c_1, c_2, c_3$ of equation (28) as input, we obtain the numerical approximations listed in Table 1, as well as the corresponding proximity map shown in Fig. 5. To generate these results, we use $N = 100,000$ regularly sampled points $\theta_k \in [-\pi, \pi)$ of the form $\theta_k = -\pi + 2\pi k/N, k = 0,1,2,\dots, N-1$.

<table><caption>**Table 1.** Numerical approximations to roots $r_1, r_2, r_3$, and theta roots $\theta_1^*, \theta_2^*, \theta_3^*$ of equation (28). Values in column $|\hat{e}(\theta_{k_l'+1}) - \hat{e}(\theta_{k_l'})|$ are the vertical distances between consecutive sample points $(\theta_{k_l'}, \hat{e}(\theta_{k_l'}))$, $(\theta_{k_l'+1}, \hat{e}(\theta_{k_l'+1}))$ in $\mathbb{R}^2$ such that $\hat{e}(\theta_{k_l'}) \hat{e}(\theta_{k_l'+1}) \le 0$ and $|\hat{e}(\theta_{k_l'}) - \hat{e}(\theta_{k_l'})| \le 1$ (see step f) in section 2). These numerical approximations were obtained with $N = 100,000$ geometric structures $L_1 L_T C(c_1, c_2, c_3, \theta_k)$, where $\theta_k = -\pi + 2\pi k/N, k = 0,1,2,\dots, N-1$.</caption><thead><tr><th>l</th><th>hat_l</th><th>theta_l^*</th><th>|hat_l - r_l|</th></tr></thead><tbody><tr><td>1</td><td>-1 - 1.000000 &times; 10<sup>0</sup>i</td><td>-1.5707963</td><td>4.537628 &times; 10<sup>-5</sup></td></tr><tr><td>2</td><td>-2 - 1.068727 &times; 10<sup>-9</sup>i</td><td>-2.6779450</td><td>5.834024 &times; 10<sup>-5</sup></td></tr><tr><td>3</td><td>1 + 2.000000 &times; 10<sup>0</sup>i</td><td>0.6435011</td><td>2.041936 &times; 10<sup>-4</sup></td></tr></tbody></table>

Comparing values in Table 1 with those in (29), we see that in this case, the cubic LC method generates excellent numerical approximations $\hat{r}_l, \hat{\theta}_l^*$; they are correct to at least six decimal places, both for their real and imaginary parts. Additionally, roots approximations $\hat{r}_l$ from Table 1 have average absolute relative error $\frac{1}{3}\sum_{l=1}^3 \left|\frac{\hat{r}_l - r_l}{r_l}\right| = 2.527864 \times 10^{-9}$, while their corresponding theta roots approximations $\hat{\theta}_l^*$ have average absolute relative error $\frac{1}{3}\sum_{l=1}^3 \left|\frac{\hat{\theta}_l^* - \theta_l^*}{\theta_l^*}\right| = 1.438584 \times 10^{-9}$.

We expect approximation errors to be larger if we use a smaller number of sample points $\theta_k$. Table 2 shows the numerical approximations obtained when using $N = 1,000$ sample points $\theta_k \in [-\pi, \pi)$ of the form $\theta_k = -\pi + 2\pi k/N, k = 0,1,2,\dots, N-1$.
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Proximity map for cubic equation $x^3 + (2 - i)x^2 + (1 - 5i)x + (2 - 6i) = 0$. Theta roots $\theta_l^*$ ($l = 1,2,3$), indicated by little diamond symbols ◇, are located at the smooth crossings of functions $e_A(\theta)$, $e_B(\theta)$ with horizontal axis $y = e(\theta) = 0$. This proximity map was obtained with $N = 100,000$ geometric structures $L_1L_T C(c_1, c_2, c_3, \theta_k)$, where $\theta_k = -\pi + 2\pi k/N$, $k = 0,1,2,\dots,N-1$.</figcaption></figure>

<table><caption>**Table 2.** Numerical approximations to roots $r_1, r_2, r_3$, and theta roots $\theta_1^*, \theta_2^*, \theta_3^*$ of equation (28). These numerical approximations were obtained with $N = 1,000$ geometric structures $L_1L_T C(c_1, c_2, c_3, \theta_k)$, where $\theta_k = -\pi + 2\pi k/N$, $k = 0,1,2,\dots,N-1$.</caption><thead><tr><th>l</th><th>r&#x0302;<sub>l</sub></th><th>&theta;&#x0302;<sub>l</sub><sup>*</sup></th><th>|ê(&theta;&#x0302;<sub>k'<sub>l+1</sub></sub>) - ê(&theta;&#x0302;<sub>k'</sub>)|</th></tr></thead><tbody><tr><td>1</td><td>-1.0000000 - 1.0000000000i</td><td>-1.5707963</td><td>4.515002 &times; 10<sup>-3</sup></td></tr><tr><td>2</td><td>-1.9999986 - 0.0000118642i</td><td>-2.6779350</td><td>5.867725 &times; 10<sup>-3</sup></td></tr><tr><td>3</td><td>0.9999193 + 1.9998553087i</td><td>0.6434742</td><td>2.048206 &times; 10<sup>-2</sup></td></tr></tbody></table>

The numerical approximations $\hat{r}_l, \hat{\theta}_l^*$ from Table 2 are correct to at least three decimal places, both for their real and imaginary parts, and they have respective average absolute relative errors

$\frac{1}{3}\sum_{l=1}^{3} \left| \frac{\hat{r}_l - r_l}{r_l} \right| = 2.669275 \times 10^{-5}$, $\frac{1}{3}\sum_{l=1}^{3} \left| \frac{\hat{\theta}_l^* - \theta_l^*}{\theta_l^*} \right| = 1.519911 \times 10^{-5}$; these average absolute relative errors are indeed larger (nearly by a factor of $10^4$) compared with the average absolute relative errors associated with Table 1. Even so, the numerical approximations from Table 2 are very good, considering that they were obtained using only a small fraction (1%) of the simultaneous geometrical structures $L_1L_T C(c_1, c_2, c_3, \theta_k)$ used to obtain the results from Table 1.
**Some features of proximity maps.** From Fig. 5, we can notice that function $e_B(\theta)$ has two “abrupt” crossings with horizontal axis $y = e(\theta) = 0$; they occur at $\theta_\alpha = -2.324985$, and $\theta_\beta = 0.8166076$. The horizontal separation $\theta_\beta - \theta_\alpha$ between these two abrupt crossings is equal to $\pi$ units. Fundamentally, this phenomenon is associated with a mathematical singularity: when the inclination angle of line $L_1$ is equal to $\theta_\alpha$ (or, equivalently, to $\theta_\beta$), in principle terminal semi-line $L_T$ would intersect circumference $C$ at $I_2 = 0$. However, since $0 \notin C$, there does not exist an intersection $I_2$ between $L_T$ and $C$ at 0, and therefore, projection $c_{0I_2(\theta_\alpha)} = -c_3/I_2$ is not defined. Now, very close to this singularity, say, at $\theta_\alpha - \epsilon$ (with $\epsilon > 0$ infinitesimally small), intersection $I_2(\theta_\alpha - \epsilon)$ is extremely close to 0, so $c_{0I_2(\theta_\alpha - \epsilon)}$ is a complex number within line $L_1$ very far away from $p_1 = -c_1/2$. At $\theta_\alpha + \epsilon$, $c_{0I_2(\theta_\alpha + \epsilon)}$ is also very far away from $p_1$, but at an opposite extreme location with respect to $c_{0I_2(\theta_\alpha - \epsilon)}$. This causes $e_B(\theta)$ to jump abruptly from one extreme value to the other in its codomain, thus causing an abrupt crossing with horizontal axis $y = e(\theta) = 0$ at $\theta = \theta_\alpha$. The analysis is completely analogous for $\theta_\beta$. For a somewhat more detailed explanation of why abrupt crossings occur in proximity maps associated with cubic polynomial equations, see numerical example 3.1 of chapter 3, in [2].

Also from Fig. 5, we can see that, for example, at $\theta = 3$, functions $e_A(\theta)$, $e_B(\theta)$ are not defined; this occurs simply because there are no intersections between terminal semi-line $L_T$ and circumference $C$ when line $L_1$ has an inclination angle $\theta = 3$. It would be desirable to have an interactive tool to observe the evolution of the geometric structure $L_1L_T C(c_1, c_2, c_3, \theta)$ associated with equation (28) as we vary the inclination angle $\theta$ of line $L_1$ with a slider control, thus confirming the behavior and features of the corresponding proximity map shown in Fig. 5. The R program listed in appendix 2 is an implementation of such interactive tool; Figs. 6 and 7 show a couple of images generated from a run of the R program listed in Appendix 2; in these images we see two states of the geometric structure $L_1L_T C(c_1, c_2, c_3, \theta)$ associated with equation (28), corresponding to two different inclination angles for line $L_1$, namely, $\theta = 0.676$ and $\theta = 1.571$.
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Geometric structure $L_1L_Tc(\theta)$ associated with polynomial $x^3 + (2-i)x^2 + (1-5i)x + (2-6i)$, for $\theta = 0.676$ (and for $\theta = 0.676 - \pi \approx -2.466$). Notice that in this situation, $l_1$ is relatively close to $r_1r_3$, so its corresponding projections $i$ and $c$ onto $L_1$ are close to $r_2$. Similarly, $l_2$ is relatively close to $r_1r_2$, so its corresponding projections $i$ and $c$ onto $L_1$ are close to $r_3$. NOTE: The subscripts of projections $i$ and $c$ onto $L_1$ are dropped in this figure (projections $i$ are associated with $L_T$, while projections $c$ are associated with $C$).</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Geometric structure $L_1L_TC(\theta)$ associated with polynomial $x^3 + (2-i)x^2 + (1-5i)x + (2-6i)$, for $\theta = 1.571$ (and for $\theta = 1.571 - \pi \approx -1.571$). Notice that in this situation, $I_2$ coincides with $r_2r_3$, so its corresponding projections $i$ and $c$ onto $L_1$ coincide with $r_1$. NOTE: The subscripts of projections $i$ and $c$ onto $L_1$ are dropped in this figure (projections $i$ are associated with $L_T$, while projections $c$ are associated with $C$).</figcaption></figure>
**4. Final thoughts and conclusions**

The cubic LC method described in this paper generally works when the roots of equation (1) are non-repeating complex numbers. However, this method also works when equation (1) has two identical complex roots. For instance, if we run the program from appendix 1 with $N = 100,000$ regularly sampled points $\theta_k = -\pi + 2\pi k/N, k = 0,1,2, ..., N-1$, but now using as input the coefficients from equation $x^3 + (0-i)x^2 + (2+8i)x + (-6-4i) = 0$, we will obtain numerical approximations $\hat{r}_1, \hat{r}_2$ to the double root $r_1 = r_2 = 1-i$ that are correct to four decimal places, both for their real and imaginary parts, and a numerical approximation $\hat{r}_3$ to the single root $r_3 = -2+3i$ that is correct to six decimal places, both for its real and imaginary part. Thus, the numerical approximations to the double root in this example are less accurate compared with the numerical approximations to the single roots obtained in section 3 with $N = 100,000$ regularly sampled points $\theta_k$; nevertheless, the numerical approximations in this double root example can still be considered good. **NOTE:** In this double root example, empirical evidence suggests that both weighted error functions $e_A(\theta), e_B(\theta)$ smoothly cross the horizontal axis $y=0$ at $\theta_1^* = \theta_2^* = -0.9827937$, while $e_B(\theta)$ smoothly crosses the horizontal axis $y=0$ at $\theta_3^* = 2.2455373$.

On the other hand, the cubic LC method described in this paper is not guaranteed to work when equation (1) has three identical complex roots. In such circumstances, line $L_1(c_1, \theta) = -c_1/2 + te^{i\theta}$ contains the origin of the complex plane when it touches the triple root $r_1 = r_2 = r_3$; as a consequence, weighted error functions $\theta \to e_A(\theta), \theta \to e_B(\theta)$ in $\mathbb{R}^2$ tend to asymptotically touch the horizontal axis $y=0$ near the triple theta root $\theta_1^* = \theta_2^* = \theta_3^*$; around this point, one of these weighted error functions behaves like a concave upward function, while the other behaves like a concave downward function.
Now, notice that $\mathbb{R}$ function **approximate_theta_root**s in appendix 1 is not designed to detect an intersection between a function $e(\theta)$ and the horizontal axis $y=0$ when $e(\theta)$ does not change sign as it touches the axis $y=0$; **approximate_theta_root**s might still be able to detect zeros of $e(\theta)$ in such circumstances, if by mere chance $e(\theta_{k'})$ is numerically equal to zero (due to numerical inaccuracies given the proximity of $L_1(c_1, \theta)$ to the origin) for a sample value $\theta_{k'}$ sufficiently close to the triple theta root $\theta_1^* = \theta_2^* = \theta_3^*$. As a concrete example, if we run the cubic LC method implemented in appendix 1 with $N = 100,000$ regularly sampled points $\theta_k = -\pi + 2\pi k/N, k = 0,1, ..., N-1$, and using as input the coefficients from cubic equation $x^3 + (9-15i)x^2 + (-48-90i)x + (-198-10i) = 0$, we will obtain very good approximations to the triple root $r_1 = r_2 = r_3 = -3+5i$ (thanks possibly to numerical errors in the computations of $e_A(\theta), e_B(\theta)$ near the triple theta root $\theta_1^* = \theta_2^* = \theta_3^* = -1.030377$), but if we run the program with $N = 1,000$ regularly sampled points $\theta_k$, we will not obtain any numerical approximations. To guarantee the detection of triple theta roots in situations like these, we must add functionality to our basic cubic LC method; one possible way to detect a triple theta root is to obtain the value $\hat{\theta}^* \in [-\pi, \pi)$ where the minimum vertical distance between $e_A(\theta)$ and $e_B(\theta)$ occurs, subject to the constraint $e_A(\theta)e_B(\theta) < 0$. For an illustration of how this added functionality to the cubic LC method works out for the approximation of a triple root associated with a specific cubic equation, see the numerical example 3.5 from chapter 3 in [2].

The cubic LC method described in this paper fails when we want to use it to approximate the roots of a cubic equation with real roots, like $x^3 + 5x^2 - 22x - 56 = 0$ (this equation has roots $4, -2, -7$).
However, if we make a change of variable, say $z = x + 2i$, then when we replace $x$ with $z - 2i$ in the original cubic equation, we get the new equation $z^3 + (5 - 6i)z^2 + (-34 - 20i)z + (-76 + 52i) = 0$, which has roots with the same real parts as the roots from the original equation, and imaginary parts all
equal to 2. If we now apply the cubic LC method to this new equation and ignore the imaginary parts
from the numerical approximations $\hat{r}_l, l = 1,2,3$, we obtain suitable numerical approximations to the
roots of the original cubic equation.

In general, numerical approximations obtained with the cubic LC method can be used as initial approximations to a well-established iterative method, like Newton-Raphson, if we want to obtain accurate numerical approximations to the roots of a cubic equation.

The cubic LC method described here is an extension of the quadratic LC method described in [1]; on the
other hand, the cubic LC method can be extended to univariate polynomial equations of degree $n \ge 4$.
For the case $n = 4$, we seek to construct simultaneous geometric structures $L_1 L_T C(c_1, c_2, c_3, c_4, \theta)$,
$\theta \in [-\pi, \pi)$, associated with the quartic equation $x^4 + c_1 x^3 + c_2 x^2 + c_3 x + c_4 = 0$; in this case, we
want to construct a terminal semi-line $L_T(c_1, c_2, c_3, \theta)$ that approaches root product $r_2 r_3 r_4$, when line
$L_1(c_1, \theta) = -c_1/2 + t e^{i\theta}$ approaches root $r_1$, and circumference $C(c_1, c_4, \theta) = c_4/L_1(c_1, \theta)$ approaches
root product $r_2 r_3 r_4$. It turns out that to construct the desired terminal semi-line $L_T$ through the use of
Vieta's formulas, it is helpful to use the expression $L_T(c_1, c_2, c_3, \theta) = -c_3 - r_x(c_2 - L_d)$, where $r_x$ is the
"best" estimate of $r_1$; to obtain the complex number $r_x$, an option is to define a "terminal curve" of the
form $\mathfrak{C}_T(c_1, c_2, c_3, \theta) = -c_3 - L_1(c_2 - L_d)$, and then compute the parametric value $t^* \in \mathbb{R}$ that
minimizes the squared distance between parametric trajectories $\mathfrak{C}_T(c_1, c_2, c_3, \theta)$ and $C(c_1, c_4, \theta)$ that
also happen to depend on $t \in \mathbb{R}$. In this way, we can obtain the best estimate of $r_1$ as
$r_x = -c_1/2 + t^* e^{i\theta}$. This process can be parallelized for simultaneously varying angles $\theta$, and in this
way, we can build proximity maps $\theta \to e_A(\theta)$, $\theta \to e_B(\theta)$ to obtain theta roots estimates $\hat{\theta}_l^*$, $l = 1,2,3,4$.
Finally, from these theta roots approximations, it is possible to estimate the roots of the quartic
equation, through the formula $\hat{r}_l = -c_1/2 + t^*(\hat{\theta}_l^*) e^{i\hat{\theta}_l^*}$, where $t^*(\hat{\theta}_l^*) \in \mathbb{R}$ is the parametric value that
minimizes the squared distance between parametric trajectories $\mathfrak{C}_T(c_1, c_2, c_3, \hat{\theta}_l^*)$ and $C(c_1, c_4, \hat{\theta}_l^*)$.

The ideas outlined above for the case $n = 4$ can be directly generalized to the case $n > 4$ without
introducing new concepts, and so, the LC method can be applied to the general polynomial equation
$x^n + c_1 x^{n-1} + c_2 x^{n-2} + \cdots + c_{n-2} x^2 + c_{n-1} x + c_n = 0, n \ge 4$. For more details on how to develop
and implement these ideas, consult chapters 4 and 5 in [2].
# Appendix 1

## R code for reproducing the results of numerical example in section 3 through the cubic LC method

The code listed in this appendix can be executed using the R statistical computing software installed on a personal computer with standard hardware resources; in this way, you can reproduce the results shown in section 3 of this paper. This R code implements the cubic LC method described in sections 1 and 2, and uses vectorization of operations, which can be seen as a basic form of parallel computing. With minor modifications, this code can approximate the roots of cubic equations other than equation (28).

```R
rm(list=ls(all=TRUE)) # clears RAM memory
```

```R
 center_zcircle <- function(w1,w2) {
```

```R
    # computes the coordinates h + ki in the complex plane that correspond
    # to the center of a circumference that passes through the origin,
    # given two nonzero points w1 and w2 from that circumference
    x1 <- Re(w1)
    y1 <- Im(w1)
    x2 <- Re(w2)
    y2 <- Im(w2)
    h <- (y2*((x1^2)+(y1^2))-y1*((x2^2)+(y2^2)))/(2*x1*y2-2*y1*x2)
    k <- (x1*((x2^2)+(y2^2))-x2*((x1^2)+(y1^2)))/(2*x1*y2-2*y1*x2)
    return( h + k*1i )
}
```

```R
intersect_semiLine_circle <- function(p,v,cc,r) {
    # computes the intersections, if any,
    # between a semi-line and a circumference in the complex plane
    # INPUTS:      p : anchor point of semi-line
    #             v : unit direction vector of semi-line
    #             cc : circumference's center
    #             r : circumference's radius
    # OUTPUTS:     I1 : intersection I1 (complex value)
    #             I2 : intersection I2 (complex value)
    #             nI : indicates which intersections are relevant
    #                 nI = 0 : I1 and I2 are both irrelevant,
    #                 or there are no intersections
    #                 nI = 1 : only one intersection is relevant,
    #                 nI = 2 : I1 and I2 are both relevant
    # NOTE: Outputs are packed into matrix OUT

    # first two columns in OUT store intersections I1, I2
    # and 3rd column in OUT is the relevance status for intersections I1, I2
    OUT <- matrix(NA,nrow=R length(p),ncol=3)
    nI <- rep(NA,length(p)) # 3rd column in OUT; we store here
                           # the relevance status for I1 and I2

    # we make sure that v is a unit vector
    v <- v/abs(v)

    # vector from p to cc
    c_p <- cc - p
```
```c
# coordinates of vector cc - p
Re_c_p <- Re(c_p)
Im_c_p <- Im(c_p)

# coordinates of semi-line's direction vector
vr <- Re(v)
vi <- Im(v)

# (cc - p) dot v
c_p Dot_v <- Re_c_p * vr + Im_c_p * vi

# squared norm of cc - p
c_p Dot_c_p <- Re_c_p * Re_c_p + Im_c_p * Im_c_p

# discriminant
delta <- c_p Dot_v^2 + r^2 - c_p Dot_c_p

# parametric values tau1, tau2 of the semi-line where
# intersections with the circumference occur
tau1 <- rep(NA, length(p))
tau2 <- rep(NA, length(p))

i <- which delta >= 0
tau1[i] <- c_p Dot_v[i] - sqrt delta[i]
tau2[i] <- c_p Dot_v[i] + sqrt delta[i]

# intersections between semi-line and circumference
I1 <- p + tau1 * v
I2 <- p + tau2 * v

# RELEVANCE CONDITIONS
# I1 and I2 are both irrelevant, or there are no intersections
ii0 <- which(tau1 < 0 & tau2 < 0 | is.na(tau1))

# only one intersection is relevant (I2)
ii1 <- which(tau1 < 0 & tau2 >= 0)

# both intersections are relevant (I1 and I2)
ii2 <- which(tau1 >= 0 & tau2 >= 0)

# assigns relevance values
nI[ii0] <- 0
nI[ii1] <- 1
nI[ii2] <- 2

# prepares output
OUT[,1] <- I1
OUT[,2] <- I2
OUT[,3] <- nI

return(OUT)
}
```
```matlab
weighted_error <- function(p,v,i0,c0) {
    # given an intersection I between a terminal semi-line Lt (derived from
    # line L1) and a circumference C that passes through the origin, this
    # function finds the weighted error e between the projection of I
    # via Lt (i0) and the projection of I via C (c0) onto line L1
    # INPUTS:      p : fixed point of line L1
    #              v : unit direction vector of L1
    #              i0 : projection of I onto L1 via Lt
    #              c0 : projection of I onto L1 via C
    # OUTPUT:       e : weighted error
    e <- sign(Re((i0-c0)/v))*abs(i0-c0)/(abs(i0-p)+abs(p-c0))
    return(e)
}
```

```matlab
L1LtC <- function(CC, theta) {
    # builds structures L1LtC for univariate cubic polynomials.
    # INPUTS:      CC : array with coefficients C1, C2, C3 of a
    #              polynomial of the form
    #              z^3 + C1*z^2 + C2*z + C3
    #              theta : array with inclination angles (with respect
    #              to positive real axis) for lines L1 with
    #              fixed point P1 = -C1/2
    # OUTPUTS:     e_A : array of weighted errors associated with
    #              intersection I1 between Lt and C; an error
    #              for each structure L1LtC
    #              e_B : array of weighted errors associated with
    #              intersection I2 between Lt and C; an error
    #              for each structure L1LtC
    #              i0I1 : array of projections onto L1 of intersections
    #              I1 between Lt theta) and C theta) via Lt
    #              i0I2 : array of projections onto L1 of intersections
    #              I2 between Lt theta) and C theta) via Lt
    #              c0I1 : array of projections onto L1 of intersections
    #              I1 between Lt theta) and C theta) via C
    #              c0I2 : array of projections onto L1 of intersections
    #              I2 between Lt theta) and C theta) via C
    #              ii1 : indicates which structures L1LtC (relative to
    #              array theta) have an intersection I1
    #              ii2 : indicates which structures L1LtC (relative to
    #              array theta) have an intersection I2
    #              iiU : indicates which structures L1LtC (relative to
    #              array theta) have a unique intersection I2
    # NOTE: outputs are packed into a list OUT

    P1 <- -CC[1]/2                      # fixed point of lines L1
    V_th <- cos theta + sin theta)*1i   # unit direction vectors of lines L1
    N <- length theta                   # number of structures L1LtC
    AP <- rep(CC[2]-P1^2,N)            # anchor points of terminal semi-lines
    V_th2 <- V_th^2                     # unit direction vectors of terminal semi-lines
```

```matlab
# centers and radii of circumferences C
P1C <- rep(-CC[3]/P1,N)              # points 1 for z-circumferences
P2C <- -CC[3]/(P1 + V_th)            # points 2 for z-circumferences
Centre_C <- center_zcircle(P1C,P2C)
Radius_C <- abs(Centre_C)
```
```R
# intersections between terminal semi-lines and z-circumferences
I_Lt_C <- intersect_semiLine_circle(AP, V_th2, Centre_C, Radius_C)
I1     <- I_Lt_C[,1]          # intersections I1
I2     <- I_Lt_C[,2]          # intersections I2
ii1    <- which(Re(I_Lt_C[,3])==2) # indices where I1 exists
ii2    <- which(Re(I_Lt_C[,3])>0)  # indices where I2 exists
iiU    <- which(Re(I_Lt_C[,3])==1) # indices where only I2 exists
```

```R
# projections of intersections between terminal semi-lines Lt and
# circumferences C onto line L1 via terminal semi-lines Lt
i0I1 <- P1 + sqrt(Re((I1[ii1]-AP[ii1])/V_th2[ii1])*V_th[ii1])
i0I2 <- P1 + sqrt(Re((I2[ii2]-AP[ii2])/V_th2[ii2])*V_th[ii2])
```

```R
# projections of intersections between terminal semi-lines Lt and
# circumferences C onto line L1 via circumferences C
c0I1 <- -CC[3]/I1[ii1]
c0I2 <- -CC[3]/I2[ii2]
```

```R
# initializes weighted error arrays
e_A <- rep(NA,N)
e_B <- rep(NA,N)
```

```R
# computes discrete proximity maps (fills weighted error arrays)
e_A[ii1] <- weighted_error( rep(P1,length(ii1)), V_th[ii1], i0I1, c0I1 )
e_B[ii2] <- weighted_error( rep(P1,length(ii2)), V_th[ii2], i0I2, c0I2 )
```

```R
# prepares outputs
OUT <- list()
OUT[[length(OUT)+1]] <- e_A
OUT[[length(OUT)+1]] <- e_B
OUT[[length(OUT)+1]] <- i0I1
OUT[[length(OUT)+1]] <- i0I2
OUT[[length(OUT)+1]] <- c0I1
OUT[[length(OUT)+1]] <- c0I2
OUT[[length(OUT)+1]] <- ii1
OUT[[length(OUT)+1]] <- ii2
OUT[[length(OUT)+1]] <- iiU
return(OUT)
}
```

```R
approximate_theta_root <- function( TH, E, TOL ) {
  # Finds the crossings of function E(TH) with horizontal axis y=0
  # in the cartesian plane.
  # Crossings are approximated via linear interpolation.
  # INPUTS : TH : array of discrete values for the independent
  #          variable
  #          E : image of discrete values TH under the real
  #          function of real variable E
  #          TOL : maximum vertical distance allowed between
  #          two discrete consecutive values of E with
  #          opposite signs, involved in a smooth
  #          crossing of E(TH) with horizontal axis y=0
  # OUTPUTS : th : abscissae of the smooth crossings of E(TH)
  #          with axis y=0 obtained by linear interpolation
  #          de : absolute vertical distances from consecutive
```
```matlab
#               elements of E involved in smooth crossings;
#               these distances are not greater than TOL
# NOTE: outputs are packed into a list L

N <- length(TH)

i <- 1:N

ic <- which(E[i]*E[i+1] <= 0 & abs(E[i]-E[i+1]) <= TOL)

th1 <- TH[ic]

th2 <- TH[ic+1]

e1 <- E[ic]

e2 <- E[ic+1]

th <- th1 + e1*(th2-th1)/(e1-e2)

de <- abs(e1-e2)

L <- list()

L[[1]] <- th

L[[2]] <- de

return(L)
}

approxP3 <- function(c1, c2, c3, N, Th_i, Th_f) {
  # approximates the roots of a univariate cubic polynomial
  # given its complex coefficients, by means of the cubic LC method
  # INPUTS :      c1,c2,c3 : coefficients of monic polynomial
  #               z^3 + c1*z^2 + c2*z + c3
  #             N : number of structures L1LtC used in the
  #               construction of the proximity map
  #             Th_i, Th_f : initial and final endpoints for
  #               theta interval of proximity map
  # OUTPUTS:      E_A : discrete values for weighted errors associated
  #               with intersections I1 between Lt and C
  #               E_B : discrete values for weighted errors associated
  #               with intersections I2 between Lt and C
  #               Theta : inclination angles of lines L1 (with respect to
  #                       positive real axis) considered in the discrete
  #                       proximity map; these angles are associated with
  #                       the corresponding components of E_A, E_B;
  #                       together they form the proximity map
  #               APPROX_LC : table of roots approximations
  #                       Rx - approximations to polynomial roots
  #                       Theta - approximations to "theta roots"
  #                       Vert - vertical distances between consecutive
  #                       discrete weighted errors involved
  #                       in the computations of theta roots
  # NOTE: outputs are packed into a list OUT

  i <- 1:N      # indices for the discrete proximity map

  # constructs a regular partition of theta interval
  Theta <- Th_i+(Th_f-Th_i)*(i-1)/N  # interval [Th_i, Th_f]

  CC <- c(c1,c2,c3) # places polynomial coefficients into an array
  OUTPUT1 <- L1LtC(CC,Theta) # generates proximity map

  E_A <- OUTPUT1[[1]]
  E_B <- OUTPUT1[[2]]
```
```c
# obtains numerical approximations of theta roots
# via crossings of functions E_A, E_B with horizontal axis y=0
# (linear interpolation is used to approximate functions E_A, E_B)
if (Th_i==pi & Th_f==pi) {
    # takes advantage of the periodicity of weighted error functions
    O1 <- approximate_theta_rootS(c(Theta,pi), c(E_A,E_A[1]), 1.0)
    O2 <- approximate_theta_rootS(c(Theta,pi), c(E_B,E_B[1]), 1.0)
} else {
    O1 <- approximate_theta_rootS(Theta, E_A, 1.0)
    O2 <- approximate_theta_rootS(Theta, E_B, 1.0)
}

# stores theta roots approximations, and vertical distances between
# associated consecutive discrete weighted errors
theta_rootS <- c( O1[[1]], O2[[1]] )
vertical_dif <- c( O1[[2]], O2[[2]] )

#######################
# obtains roots approximations Rx,
# corresponding to theta roots approximations

# step 1: obtains structures L1LtC
# associated with theta roots approximations
OUTPUT2 <- L1LtC(CC,theta_rootS)

E_Ar    <- OUTPUT2[[1]]
E_Br    <- OUTPUT2[[2]]
i0I1    <- OUTPUT2[[3]]
i0I2    <- OUTPUT2[[4]]
c0I1    <- OUTPUT2[[5]]
c0I2    <- OUTPUT2[[6]]
ii1     <- OUTPUT2[[7]]
iiU     <- OUTPUT2[[9]]

# step 2: averages projections i0 and c0 to obtain Rx
Rx      <- rep(NA,length(Theta_rootS))
Rx[ii1] <- ifelse(abs(E_Ar[ii1])<abs(E_Br[ii1]),
                 (i0I1[ii1]+c0I1[ii1])/2,
                 (i0I2[ii1]+c0I2[ii1])/2)
Rx[iiU] <- (i0I2[iiU]+c0I2[iiU])/2

# listing of roots approximations, theta roots approximations,
# and vertical distances between consecutive discrete
# weighted errors involved in the computation of theta roots
APPROX_LC <- setNames(data.frame(
  Rx[order(Vertical_dif)],
  theta_rootS[order(Vertical_dif)],
  Vertical_dif[order(Vertical_dif)] ),
  c("Rx","Theta","Vert"))

# prepares output
OUT      <- list()
OUT[[length(OUT)+1]] <- E_A
OUT[[length(OUT)+1]] <- E_B
OUT[[length(OUT)+1]] <- Theta
OUT[[length(OUT)+1]] <- APPROX_LC
return(OUT)
}
```
# true polynomial roots (to compare them with approximations from LC method)
R1 <- -1-1i
R2 <- -2+0i
R3 <-  1+2i
R  <- c(R1,R2,R3)

# coefficients of polynomial z^3 + C1*z^2 + C2*z + C3 from numerical example
C1 <- 2-1i
C2 <- 1-5i
C3 <- 2-6i

# true theta roots (to compare them with approximations from LC method)
ThetaR <- atan2(Im(R+C1/2),Re(R+C1/2))

###### approximates polynomial roots by means of cubic LC method
dev.new()
plot(Theta, E_A, 'l', xlim=c(-pi,pi), ylim=c(-1,1), col='red',
lty=1, lwd=2, xlab=expression( theta ), ylab=expression( e( theta )))
abline(h=seq(from=-1,to=1,by=0.25), col='light grey')
abline(h=0, col='black')
abline(v=seq(from=-4,to=4,by=0.25), col='light grey')
lines(Theta, E_A, 'l', col='red', lty=1, lwd=2)
lines(Theta, E_B, 'l', col='blue', lty=1, lwd=2)
points(ThetaR, rep(0,length(ThetaR))), col='black', pch=18)
points(ThetaR, rep(0,length(ThetaR))), col='black', pch=5, cex=1.3)
legend(x='bottomright',
legend=c(expression( e[A](theta))), expression( e[B](theta))),
col=c('red','blue'),
lty=c(1,1),
lwd=c(2,2)
)

# graphs functions of weighted errors (proximity maps).
# Adds points indicating the true theta roots,
# to compare them with graphical approximations
dev.new()
plot(Theta, E_A, 'l', xlim=c(-pi,pi), ylim=c(-1,1), col='red',
lty=1, lwd=2, xlab=expression( theta ), ylab=expression( e( theta )))
abline(h=seq(from=-1,to=1,by=0.25), col='light grey')
abline(h=0, col='black')
abline(v=seq(from=-4,to=4,by=0.25), col='light grey')
lines(Theta, E_A, 'l', col='red', lty=1, lwd=2)
lines(Theta, E_B, 'l', col='blue', lty=1, lwd=2)
points(ThetaR, rep(0,length(ThetaR))), col='black', pch=18)
points(ThetaR, rep(0,length(ThetaR))), col='black', pch=5, cex=1.3)
legend(x='bottomright',
legend=c(expression( e[A](theta))), expression( e[B](theta))),
col=c('red','blue'),
lty=c(1,1),
lwd=c(2,2)
)

# displays numerical results on R console
# true roots
R

# true theta roots
ThetaR

# table of approximations from cubic LC method
ROOTS_LC3_Table
ťa root of a cubic polynomial; they must correspond with input
# coefficients (red points on the graph, shown only as a reference)
R1 <- -1-1i
R2 <- -2+0i
R3 <-  1+2i

# double products of roots
# (green points on the graph, shown only as a reference)
Q1 <- R1*R2
Q2 <- R1*R3
Q3 <- R2*R3

# coefficients of polynomial z^3 + C1*z^2 + C2*z + C3,
# used to build geometric structures L1-Lt-zC
C1 <- 2-1i
C2 <- 1-5i
C3 <- 2-6i

# fixed points associated with polynomial coefficients
# (black points on the graph)
P1 <- -C1/2              # fixed point for line L1
P2 <- -C3/P1            # fixed point for circumference zC
P3 <- C2 - P1^2         # fixed point for terminal semi-line Lt

ťa root of a cubic polynomial; they must correspond with input
# inclination angle of line L1 through a slider-type control.
# We can observe line L1, terminal semi-line Lt, and z-circumference
# zC; additionally, we can also see the projections of the
# intersections between Lt and zC onto line L1 (red and blue points:
# projections associated with Lt; red and blue diamonds: projections
# associated with zC; color red is associated with intersection I1,
# and color blue is associated with intersection I2; see function
# intersect_semiLine_circle above in this script for more details)

Sliderstart = 0.000
Slidermin   = 0.000
Slidermax   = pi
Sliderstep  = 0.001

tt <- tktplevel()
SliderValue <- tclVar(Sliderstart)

graphic <- function(...) {

  # inclination angle of L1
  THETA1=as.numeric(tclvalue(SliderValue))

  # direction vector of L1
  V <- cos(THETA1) + 1i*sin(THETA1)

  # a moving point on L1
  PL <- P1 + V

Sliderstart = 0.000
Slidermin   = 0.000
Slidermax   = pi
Sliderstep  = 0.001
tt <- tktplevel()
SliderValue <- tclVar(Sliderstart)

graphic <- function(...) {

  # inclination angle of L1
  THETA1=as.numeric(tclvalue(SliderValue))

  # direction vector of L1
  V <- cos(THETA1) + 1i*sin(THETA1)

  # a moving point on L1
  PL <- P1 + V
PLt <- C2 - (P1-100*V)*(P1+100*V)

# a moving point on z-circumference zC
PC <- -C3/PL

# gets coordinates of a discrete set
# of points in the z-circumference
CENTER_ZC <- center_zcircle(P2,PC)
ZC <- points_todraw_z_circle(CENTER_ZC, N=3000)

# prepares graphics area
par(mar=c(1.9,1.9,1.9,1.9))
plot(0, 0, 'n', xlim=c(-5,5), ylim=c(-5,5), asp=1)
abline(h=0) # draws real axis
abline(v=0) # draws imaginary axis

# draws line L1 (red color)
t <- seq(from=-1000, to=1000, by=2000)
lines(P1+t*V,col="red")

# draws z-circumference zC (green color)
lines(ZC, col='green3')

# draws terminal semi-line Lt = C2-LD1 (black color)
segments(Re(P3),Im(P3),Re(PLt),Im(PLt))

# draws polynomial roots (red points)
points(R1,pch=20,col="red")
points(R2,pch=20,col="red")
points(R3,pch=20,col="red")

# draws fixed points of L1, Lt, and zC
# (black points)
points(P1,pch=20,col="black")
points(P2,pch=20,col="black")
points(P3,pch=20,col="black")

# draws double products of roots
# (green points)
points(Q1,pch=20,col="green3")
points(Q2,pch=20,col="green3")
points(Q3,pch=20,col="green3")

# computes intersections between terminal semi-line Lt
# and z-circumference zC
I <- intersect_semiLine_circle(P3, V^2, CENTER_ZC, abs(CENTER_ZC))

# in case there are intersections between
# terminal semi-line Lt and z-circumference zC
if (Re(I[1,3])>0) {
  I1 <- I[1,1]
  I2 <- I[1,2]
# is I1 on terminal semi-line Lt?

if (Re(I[1,3])==2) {

# Computes projections of I1 onto L1, via terminal semi-line Lt

i11 <- P1 - sqrt(abs(I1-P3))*V

i12 <- P1 + sqrt(abs(I1-P3))*V

# Computes projection of I1 onto L1, via z-circumference zC

c1 <- -C3/I1

# draws intersection I1, and their respective projections onto L1

# (red points for projections associated with terminal semi-line Lt,
# and red diamond for projection associated with circumference zC)

points(I1, pch=20, col="red")

points(i11, pch=20, col="red")

points(i12, pch=20, col="red")

points(c1, pch=5, col="red")

# draws text labels (algebraic expressions)

# for intersection I1 and its projections onto line L1

text(I1, pos=2, expression(I[1]))

text(i11, pos=2, expression(i))

text(i12, pos=2, expression(i))

text(c1, pos=1, expression(c))

}

# Computes projections of I2 onto L1, via terminal semi-line Lt

i21 <- P1 - sqrt(abs(I2-P3))*V

i22 <- P1 + sqrt(abs(I2-P3))*V

# Computes projection of I2 onto L1, via z-circumference zC

c2 <- -C3/I2

# draws intersection I2, and their respective projections onto L1

# (blue points for projections associated with terminal semi-line Lt,
# and blue diamond for projection associated with circumference zC)

points(I2, pch=20, col="blue")

points(i21, pch=20, col="blue")

points(i22, pch=20, col="blue")

points(c2, pch=5, col="blue")

# draws text labels (algebraic expressions)

# for intersection I2 and its projections onto line L1

text(I2, pos=2, expression(I[2]))

text(i21, pos=2, expression(i))

text(i22, pos=2, expression(i))

text(c2, pos=1, expression(c))

}

# draws text labels (algebraic expressions)

# for some fixed points

text(P1, pos=4, expression(p[1]))

text(P2, pos=4, expression(-c[3]/p[1]))

text(P3, pos=4, expression(c[2]-p[1]^2))

text(R1, pos=4, expression(r[1]))

text(R2, pos=4, expression(r[2]))
```r
text(R3, pos=4, expression(r[3]))
text(Q1, pos=4, expression(r[1]*r[2]))
text(Q2, pos=4, expression(r[1]*r[3]))
text(Q3, pos=4, expression(r[2]*r[3]))

# draws text labels for geometric elements
text(P1+4.5*V, pos=2, expression(L[1]))
text(C2-(P1-2*V)*(P1+2*V), pos=2, expression(L[T]))
text(-C3/(P1-0.25*V), pos=2, expression(C))

} # ends definition of function graphic

# builds graph of type tcl/tk
img <- tkrplot.tt, graphic, hscale=1.85, vscale=1.85)
showimage <- function(...) tkrreplot(img)
scl = tkscale.tt, from=Slidermin, to=Slidermax,
    showvalue=TRUE, variable=SliderValue, resolution=Sliderstep,
    command=showimage, orient='vertical')
tkpack(scl, side='left')
tkpack(img, side='right')
```

## References

[1] Daniel Alba-Cuellar. Locating the roots of a quadratic equation in one variable through a Line-Circumference (LC) geometric construction in the plane of complex numbers. Maple Trans. 4, 2, Article 17589 (2024), 13 pages. https://doi.org/10.5206/mt.v4i2.17589

[2] Alba-Cuellar, Daniel (2024). The LC Method: A parallelizable numerical method for approximating the roots of single-variable polynomials. arXiv preprint arXiv:2402.15554

[3] Tierney, Luke. (2025). tkrplot: TK Rplot. R package version 0.0-29. https://cran.r-project.org/web/packages/tkrplot/index.html