Thomas Gautier
email: thomas.gautier@latmos.ipsl.fr

Joseph Serigano

Jérémy Bourgalais

Sarah M Hörst

Melissa G Trainer

Decomposition of electron ionization mass spectra for space application using a Monte-Carlo approach

Keywords: Electron Ionization, Deconvolution, Space Instrumentation, Quadrupole mass spectrometer, Monte-Carlo

ABSTRACT

RATIONALE

Quadrupole mass spectrometers equipped with an Electron Ionization (EI) sources have been widely used in space exploration to investigate the composition of planetary surfaces and atmospheres. However, the complexity of the samples and the minimal calibration for the fragmentation of molecules in the ionization chambers have prevented the deconvolution of the majority of the mass spectra obtained at different targets, thus limiting determination of the exact composition of the samples analyzed. We propose a Monte-Carlo approach to solve this issue mathematically.

METHODS

We have decomposed simulated mass spectra of mixtures acquired with unit resolving power mass spectrometers and EI sources into the sum of the single components fragmentation patterns weighted by their relative concentration using interior-point least-square fitting. To fit compounds with poorly known fragmentation patterns, we use a Monte-Carlo method to vary the intensity of individual fragments. We then decompose the spectrum thousands of times to obtain a statistical distribution.

RESULTS

By performing the deconvolution on a mixture of seven different molecules with interfering fragmentation patterns (H2O, O2, CH4, Ar, N2, C2H4 and C2H6) we show that this approach retrieves the mixing ratio of the individual components more accurately than regular mass spectra decomposition methods that rely on fragmentation pattern from general databases. It also provides the probability density function for each species' mixing ratio.

CONCLUSIONS

By removing the solution degeneracy in the decomposition of mass spectra, the method described herein could significantly increase the scientific retrieval from archived space flight mass spectrometry data, where calibration of the ionization source is no longer an option.

INTRODUCTION

The majority of the mass spectrometers (MS) that have flown in space to date have had a resolving power close to the unit resolution, and have used an electron ionization (EI) source to ionize sample gases prior to their mass analysis. The prevalence of this approach is mainly due to the volume, mass, power, and robustness requirements of this subsystem, all of which are critical for space exploration 1,2 . To form ions through electron ionization, the incoming gas sample interacts with an electron beam at a given mean energy, typically 70 eV. It is well known that upon such ionization most, if not all, molecular compounds will fragment. This yields a mass spectrum for a molecular compound comprised not only of one peak at the nominal m/z but of several peaks at different m/z, resulting in a fragmentation pattern for each molecule analyzed. The position of the peaks in the fragmentation pattern is characteristic of the original molecule and the incident electron energy. This means that the fragmentation pattern observed in a sample can be used to trace back to the original molecule analyzed and lead to its identification (Figure 1). While this approach works extremely well for single compounds, it becomes harder to identify the analyzed material in the case of complex mixtures that could be composed of dozens of different molecular compounds. One approach for analyzing complex mixtures is to separate each individual component, for example using gas chromatography (GC) prior to the MS analysis.

However, this approach is not always viable in space exploration where adding a GC to a mass spectrometer adds several kilograms and power in a finite resource environment. This explains why many of the mass spectrometers previously used for space exploration were not coupled to a chromatographic system. There have also been instances in which a GC-capable MS has been operated in MS-mode only due to sparse samples and the need to maximize sensitivity. This was the case for the COSAC (Cometary Sampling and Composition Experiment) instrument aboard the Philae lander of the Rosetta mission that landed on comet 67P/Churyumov-Gerasimenko in 2014 [START_REF] Goesmann | Organic compounds on comet 67P/Churyumov-Gerasimenko revealed by COSAC mass spectrometry[END_REF] and the GCMS (Gas Chromatograph Mass Spectrometer) instrument aboard the Huygens probe that landed on Titan, the largest satellite of Saturn, in 2005 [START_REF] Niemann | Composition of Titan's lower atmosphere and simple surface volatiles as measured by the Cassini-Huygens probe gas chromatograph mass spectrometer experiment[END_REF][START_REF] Niemann | The abundances of constituents of Titan's atmosphere from the GCMS instrument on the Huygens probe[END_REF] .

When only a mass spectrum of the sample is available, the results tend to be an intricate mix of mass peaks corresponding to the addition of multiple fragmentation patterns of the molecules present in the sample (Figure 1). However, while the peak positions in a molecule fragmentation pattern are predictable, their relative intensities are highly dependent of the ion source geometry [START_REF] Mclafferty | Interpretation of Mass Spectra[END_REF] . This explains the variability in the relative weight of fragments in the mass spectra of compounds even with laboratory grade instruments on Earth databases [START_REF] Ausloos | The critical evaluation of a comprehensive mass spectral library[END_REF][START_REF] Linstrom | NIST Mass Spectrometry Data Center[END_REF][START_REF] Wallace | Mass Spectral Library Quality Assurance by Inter-Library Comparison[END_REF][START_REF] Valdez | Assessing the reliability of the NIST library during routine GC-MS analyses: Structure and spectral data corroboration for 5,5diphenyl-1,3-dioxolan-4-one during a recent OPCW proficiency test[END_REF] . This variability is typically of 20-30% of a fragment intensity but can go much higher in extreme cases. For example, for methane, there is a factor 3.5 in the peak intensity of fragment CH2 + at m/z 14 between the value from the NIST database (20.4 % of the base peak at m/z 16) and the value from pre-flight calibration of the Ion and Neutral Mass Spectrometer aboard Cassini (5.7% of m/z 16). Thus, for space flight instruments, the ideal approach would be pre-flight calibration of all the molecules to be detected by the instrument. However, while a few compounds are usually calibrated, this approach is highly unrealistic on a large scale (1) due to the a-priori unknown composition of the samples to be studied, (2) the limited time available for the instrument calibration during its integration, and (3) the desire to prevent contamination of instrument surfaces that would preclude sensitive detections at the target.

The fact that the fragmentation patterns highly depend on the instrument used and that flight instruments are insufficiently calibrated introduces a strong degeneracy in the mathematical decomposition of the mass spectra, and explains why most mass spectra acquired by space instruments have not been fully deconvolved so far, with only very limited counter examples [START_REF] Goesmann | Organic compounds on comet 67P/Churyumov-Gerasimenko revealed by COSAC mass spectrometry[END_REF][START_REF] Cui | Analysis of Titan's neutral upper atmosphere from Cassini Ion Neutral Mass Spectrometer measurements[END_REF] .

In this paper, we present a novel approach to deconvolution, introducing a Monte-Carlo sampling on the expected molecules fragmentation patterns, to overcome the degeneracy inherent within compositional analysis.

THEORY AND CALCULATION

FRAGMENTATION PATTERNS

Since the pressure in an EI ion source is kept low (typically around 1x10 -6 mbar) ensuring a mean free path sufficiently large (ca. 100 m) to avoid bimolecular and collisional reactions. The only possible interaction remaining is between electrons (typically at 70 eV) and the neutral gas phase molecules from the sample. Since these interactions take place at different locations within the ion source, each resulting ion formed from a neutral molecule (M) could have a variety of internal energies. Ions formed with low internal energy will be directly detected as the parent peak of the molecule (M + ) in the spectrum. However, if a molecule is ionized with a sufficient initial energy to reach a dissociative ionization threshold, it will decompose through different reactions depending on its internal energy to generate fragments of the original molecule. In some cases, the fragment ion will have sufficient remaining internal energy to further decompose into an even smaller ion. The decomposition processes explain the fragmentation pattern of a given compound, while the stability of the formed ions and the geometry of the interaction between the incoming neutral molecule and the electron beam explain the variation in fragment intensities depending on the measuring instrument. For further information on the fragmentation processes, we suggest reading the detailed chapters 4, 7 and 8 of Interpretation of mass spectra by McLafferty and Turecek 6 .

MASS SPECTRA DECOMPOSITION

Deconvolving a mass spectrum to identify and quantify the different molecules present in a sample relies on being able to solve the linear system:

𝑀𝑎𝑠𝑠 𝑆𝑝𝑒𝑐𝑡𝑟𝑢𝑚 = 𝐹𝑟𝑎𝑔𝑚𝑒𝑛𝑡𝑎𝑡𝑖𝑜𝑛 𝐷𝑎𝑡𝑎𝑏𝑎𝑠𝑒 × 𝑆𝑝𝑒𝑐𝑖𝑒𝑠 𝑐𝑜𝑛𝑐𝑒𝑛𝑡𝑟𝑎𝑡𝑖𝑜𝑛𝑠 (1)
which equates to, in an ideal case, solving equation ( 1):

( 𝐼 𝑚/𝑧 1 𝐼 𝑚/𝑧 2 ⋮ 𝐼 𝑚/𝑧𝑗 ) = ( 𝐹 1,1 𝐹 2,1 𝐹 1,2 𝐹 2,2 𝐹 1,𝑖 𝐹 2,𝑖 ⋮ ⋱ ⋮ 𝐹 𝑗,1 ⋯ 𝐹 𝑗,𝑖 ) × [ 𝑁 1 𝑁 2 ⋮ 𝑁 𝑖 ] (2)
where Im/z j is the ion count in the measured mass spectrum for a mass bin j, Ni is the concentration of each species i arriving on the detector of the spectrometer, and Fj,i is the fragment intensity for species i in the mass bin j.

Several methods exist to solve such system. For example, the single value decomposition has been used to treat data from the Cassini Ion and Neutral Mass Spectrometer (INMS) instrument [START_REF] Cui | Analysis of Titan's neutral upper atmosphere from Cassini Ion Neutral Mass Spectrometer measurements[END_REF] , while least-square fitting was used for the data from the COSAC instrument on the Philae lander [START_REF] Goesmann | Organic compounds on comet 67P/Churyumov-Gerasimenko revealed by COSAC mass spectrometry[END_REF] . In our case we choose to use least-squares fitting on multiple linear regressions, which can be more time consuming but allows better constraint of the solution if boundary conditions are known on some species. Equation ( 1) is thus solved by finding the least squares minimum of the Euclidian norm of (2) as described in (3) and using condition (4):

min 𝑁 1 2 ‖𝐹 × 𝑁 -𝑀𝑆‖ 2 2 (3)
with

𝐿 ≤ 𝑁 ≤ 𝑈 ( 4
)
where F is the matrix containing the fragmentation database of the molecules suspected to be present in the sample; N is the output vector containing their relative concentration at the detector of the spectrometer; MS is the measured mass spectrum; and L and U are respectively the vectors containing the lower limit and upper limit of molecules in N, both of which must be positive.

A variety of numerical recipes exist to minimize Equation (3) but we chose to use an interiorpoint algorithm that easily accommodates large matrices such as mass spectra [START_REF] Potra | Interior-point methods[END_REF][START_REF] Woitke | Radiation thermo-chemical models of protoplanetary disks[END_REF][START_REF] Karmarkar | A new polynomial-time algorithm for linear programming[END_REF][START_REF] Vanderbei | Symmetric indefinite systems for interior point methods[END_REF][START_REF] Altman | Regularized symmetric indefinite systems in interior point methods for linear and quadratic optimization[END_REF] . However, mathematically decomposing a mass spectrum only works if the fragmentation patterns are well known. Failure to constrain the exact fragment ratios for all molecules present in the sample would most often result in a degeneracy of the solutions. Such degeneracy would render any relative quantification of the molecules in the sample questionable.

MONTE-CARLO APPROACH

To overcome the degeneracy in mass spectra deconvolution we introduced a Monte-Carlo approach in our algorithm. This approach has been successfully applied in planetary sciences to further refine photochemical models of planetary atmospheres such as Earth, Neptune, Saturn and Titan by randomizing the poorly known reaction rate coefficients [START_REF] Stolarski | Uncertainty propagation in a stratospheric model 2. Monte Carlo analysis of imprecisions due to reaction rates[END_REF][START_REF] Dobrijevic | Effect of chemical kinetics uncertainties on hydrocarbon production in the stratosphere of neptune[END_REF][START_REF] Dobrijevic | Effect of chemical kinetic uncertainties on photochemical modeling results: Application to Saturn's atmosphere[END_REF][START_REF] Hébrard | Photochemical kinetics uncertainties in modeling Titan's atmosphere: First consequences[END_REF] .

In our case the randomization is applied to the nominal fragmentation pattern of each molecule in our database to create n different fragmentation databases F, typically n=10,000 for the work presented here, and then decompose the measured mass spectrum with each fragmentation database individually. Out of these 10,000 runs, the best 5% (i.e. with minimal residuals) are kept to retrieve a statistical solution of the problem. This allows retrieving not only the most probable concentrations for a given species, but also its probability distribution function.

The randomization process is only applied to the intensity of the fragments, not their position in the mass spectrum; i.e. N2 will always have peaks at m/z 7, 14 and 28 (corresponding to N 2+ , N + and N2 + ions), and only there, but their relative intensity will be varied in each Monte-Carlo run. The relative intensity variation is made within a chosen range. Typically for molecules whose fragmentation pattern is only constrained from mass spectral databases, we allow a variation, d, of ±30% for each fragment's intensity around its nominal value. For every run n, each member Fj,i of the fragmentation database matrices F can thus be expressed as :

𝐹 𝑗 𝑛 ,𝑖 𝑛 = 𝐹 𝑗 1 ,𝑖 1 × (1 + 𝜀) (5)
with

-𝑑 ≤ 𝜀 ≤ 𝑑 (6)
where 𝐹 𝑗 1 ,𝑖 1 corresponds to the nominal fragmentation of the molecules and 𝜀 is a pseudorandom number comprised between -1 and + 1 and generated using a uniform pseudorandom number generator based on the Mersenne Twister method. To avoid a priori on the validity of the measurement, we chose to use uniformly pseudorandom distributed numbers for 𝜀. Actually it is very likely that epsilon do not follow a uniform distribution but a normal one, centered around the experimental value determined with the same instrument. As our method intends to correct for the poor knowledge of this value when the measurement is performed with a uncalibrated instrument (and thus where epsilon is taken for example from the NIST database), we chose not to make the assumption of a normal distribution. Instead, since we cannot a priori know the bias of the database toward the actual measurement, we chose to use a uniform distribution. While this ensures that the algorithm does not favor an outcome based on false a priori, it also induces larger uncertainties on the retrieval. In cases where sufficient pre-flight calibrations are available, one would be able to determine more precisely the epsilon distribution, which would in return drastically reduce the uncertainties of our method. These pre-flight calibrations database are essential and should always be favored if available.

IONIZATION CROSS SECTIONS

Another parameter of importance to take into account when retrieving the relative concentration of molecules in a sample from an EI mass spectrum is that the intensity of a peak does not directly correlate with the concentration of a molecule, but is a product of its concentration times its ionization cross section. The molecular ionization cross sections (σ) depend on the neutral molecule ionized and on the mean energy of the electrons in the ion source (typically 70 eV). If unavailable in the literature, such cross sections can be approximated using either empirical formulae [START_REF] Fitch | Calculation of relative electron impact total ionization cross sections for organic molecules[END_REF] or through a Binary-Encounter-Dipole (BED) or Binary-Encounter Bethe (BEB) model for electron ionization [START_REF] Kim | Binary-encounter-dipole model for electron-impact ionization[END_REF][START_REF] Hwang | New model for electron-impact ionization cross[END_REF] . The origin of the electron ionization cross section (experimental, empirical, BED or BEB) is not critical. However, since relative ionization probabilities are required to compare abundances, one must be sure to use the same origin for all the cross sections in their database to avoid introducing biases. Given the lack of availability for experimental data on most of the species of interest for planetary science, we recommend using the BED method for such applications. For each species in the sample, the concentration decomposed by the model, N, can thus be expressed as a function of its actual concentration in the sample, C, as follows:

( 𝑁 1 ⋮ 𝑁 𝑖 ) = 𝐴 . ( 𝜎 1 ⋮ 𝜎 𝑖 ) . ( 𝐶 1 ⋮ 𝐶 𝑖 ) ( 7
)
where A is an instrumental function, assumed constant for all of the molecules. Equation 7can thus be transformed to express the relative concentration of two molecules x and y as:

𝐶 𝑥 𝐶 𝑦 = 𝑁 𝑥 𝑁 𝑦 . 𝜎 𝑦 𝜎 𝑥 . ( 8
)
Numerous steps are thus needed before retrieving the concentration of molecules in a sample from a direct mass spectrum measurement. The general workflow to achieve this using our method is summarized in Figure 2.

RESULTS AND DISCUSSION

In order to evaluate the performance of our algorithm we ran it on a synthetic mass spectrum of a mixture of 7 speciesnamely argon, molecular nitrogen, molecular oxygen, methane, ethane, ethylene, and watereach of which was randomly assigned a mixing ratio in the mixture, as given in Table 1, column 2. These species were chosen for testing as their fragmentation patterns particularly interfere with each other in the m/z 10 to 30 range. To test the limitations of our code, the sample mass spectrum was not reconstructed directly from the NIST fragmentation patterns of the 7 molecules. Instead, using the method presented in equations ( 5) and ( 6), we randomly varied their fragments with d ± 30%, and then reconstructed a mass spectrum from the newly generated fragmentation patterns. The fragmentation database used to build the synthetic mass spectrum is given in Table SI 1 and the synthetic mass spectrum in Table SI 2. The synthetic mass spectrum used to test our approach is presented in Figure 3 (black line) with the calculated relative contribution of each species (colored bars). As visible in Figure 3, our algorithm is able to precisely reconstruct the mass spectrum, with the minor exceptions of m/z 24, 25, 30 and 31, indicating inconsequential mismatches in the fragmentation patterns of ethylene and ethane. In order to fully evaluate the performance of our decomposition algorithm we performed three different deconvolutions using: (1) the nominal fragmentation patterns of these molecules from the constructed database given in Table SI 1, which corresponds to having a custom database on the fragmentation in an instrument's ionization chamber and should, as the most accurate method, always be favored if available;

(2) using the nominal fragmentations patterns of the molecules as given by the NIST database [START_REF] Linstrom | NIST Mass Spectrometry Data Center[END_REF] (simulating a decomposition with poorly known fragmentations); and (3) using the Monte-Carlo approach described here. A comparison of the retrieved concentration is summarized in (1) Compared to the exact mixing ratio introduced in the synthetic mass spectrum, the mathematical decomposition of the mass spectrum using exactly known fragmentations yields, by far, the best output. This method allows for the exact determination of the molecular ratio for every single species. This was expected and, mathematically speaking, would remain true for complex mixtures up to z compounds, as long as the number of peaks with signal in the input mass spectrum is greater than z and the fragmentation patterns of the molecules are perfectly characterized with the exact same ionization chamber used to measure the samples.

(2) While the NIST database approach allows the analyst to get a fairly good idea of the sample composition, using fragmentation patterns directly from the NIST database yields the largest inaccuracies in the absolute mixing ratios. Methane and ethane are well determined, but water is severely overestimated. This comes from the differences in the fragmentation patterns between the database and the actual mass spectrum. This is representative of what would happen when deconvolving a mass spectrum measured in space using Earth databases with no fragmentation corrections.

(3) The Monte-Carlo method we propose uses the same fragmentation pattern database as the previous one (NIST), but attempts to correct for the discrepancies by varying the fragmentation patterns a large number of times and deconvolving the spectrum with each variation of the database to run a statistical analysis on the deconvolution. It can be used in two different ways.

The first one is to determine the molecular mixing ratio by taking the average output of the top 5% of simulations (Table 1, column 5), calculated on the least-square regression residue. Note that this 5% threshold is not a requirement of our approach, but mainly a cut-off we introduced to limit computational time of our method. Using this method improves the information one can have of the sample. Compared to method 2, while it slightly degrades the retrieval of ethane mixing ratio, it significantly refines the retrievals of water and ethylene, and improves the retrievals of argon and oxygen. The second use, results of which populate column 6 of Table 1, is to probe the probability density functions for each species as illustrated in Figure 4. The determination of maximum probability density functions is the method we favor for better accuracy. However, this approach is more time consuming than the top 5% method. The probability density functions reflect the likelihood of a species to be at a given mixing ratio in the sample.

These probability densities allow for a quick determination of the parameter space occupied by each species in the sample as well as the determination of a confidence level in the retrieved mixing ratio (Table 1, column 6). Doing so demonstrates the power of this method, which can then accurately retrieve all species in Table 1 but ethane within the predicted 1σ variation. The probability density functions can also be used to probe the interdependencies between species in their mixing ratio retrievals. In our example, using a mixing ratio for ethylene in the lower side of its distribution, automatically induces a higher mixing ratio for ethane. Should ethylene mixing ratio be determined more accurately in a sample using another instrument, one could then infer the accurate mixing ratio of ethane too.

Given our observations, the proposed Monte-Carlo approach seems then to be an extremely valuable tool in mass spectral decomposition that could be applied to flight mass spectra such as the data recorded by the INMS and GCMS instruments during the Cassini-Huygens mission.

CONCLUSION

While numerous mass spectrometers have been flown in space over the last few decades, most of their data have not be deconvolved and remained under used due to the complexity of the analyzed samples and the lack of essential pre-flight calibration of the unique fragmentation patterns of molecules in the instruments.

We propose a method for the analysis of these data using Earth-based determination of the fragmentation patterns of the molecules, readily available in databases. To do so, we correct for While we presented implications for space exploration mass spectrometry, the method proposed in this paper could also be applied to laboratory measurements performed with Residual Gas Analyzers, typically quadrupole mass spectrometers equipped with an electron ionization source, when direct calibration of the component is not possible.

We also believe the method presented here can significantly improve the science retrieval in archived space flight data. For such applications, we want to point out that in cases with low count ratesfor example the most intense peaks in the mass spectrum recorded by the COSAC instrument at the surface of comet 67P were only at ~20 countsthe counting statistics would introduce an additional source of uncertainty. In such cases, this source of error would need to be combined with the fitting ones to retrieve proper uncertainties on the relative abundances in the sample.

Finally, we also confirm that the most accurate method, by far, to determine a sample composition analyzed by a space mass spectrometer relies on the a priori knowledge of the exact fragmentation pattern of every single suspected compound with the exact same instrument. We thus reaffirm the need for extended pre-launch calibrations of future flight mass spectrometers as practicable.

Figure 1 :

 1 Figure 1: Example of the fragmentation of molecular compounds (nitrogen in blue, methane in red and ethane in orange) in an EI ionization source and how the measured mass spectrum (black, top) is the sum of the fragmentation patterns of all the molecules present in the sample weighted by their relative concentrations and ionization cross-sections.

Figure 2 :

 2 Figure 2: Schematic representation of the Monte-Carlo method used to achieve the decomposition of EI mass spectra without instrumental calibration of fragmentation patterns.

Figure 3 :

 3 Figure 3 : Example of mass spectrum decomposition using our Monte-Carlo approach. The synthetic mass spectrum used for the decomposition is shown with the black lines. The colored bars show the calculated contribution of each of the molecules to the mixture: methane (blue), water (orange), nitrogen (yellow), ethylene (purple), ethane (green), argon (cyan) and oxygen (red).

2 (columns 1 and 2 ) and the retrieved mixing ratios from the different decomposition 3 methods: using exactly known fragmentation patterns (column 3 ); using poorly known 4 fragmentation pattern (column 4 ); using a randomized distribution of poorly known 5 fragmentation pattern using a Monte-Carlo method and calculating the average of the top 6 5% of 10 ,

 2233445610 000 decompositions (column 5); and using the Monte-Carlo method but selecting 7 the maximum of the probability density function and its 1σ standard deviation (column 8

9

 9

Figure 4 :

 4 Figure 4: Probability densities of the mixing ratio of the molecules in the sample as retrieved by the Monte-Carlo algorithm

Table 1 .

 1

				Decomposition using:
	Synthetic mass spectrum	Known fragmentations	NIST fragmentation	Monte-Carlo (mean top 5%)	Monte-Carlo (max pdf)
	Compounds	Mixing ratio (%)	Mixing ratio (%)	Mixing ratio (%)	Mixing ratio (%)	Mixing ratio (%) ±1σ
	CH4	23.5	23.5	24.6	21.2	21.7 ±1.9
	H2O	22.1	22.1	31	24.2	22.4 ±3.2
	N2	7.9	7.9	5.2	5.4	7.2 ±3.6
	C2H4	21.4	21.4	25	21.6	22.5 ±2.8
	C2H6	6.3	6.3	5.9	3.7	3.5 ±1.7
	Ar	5.5	5.5	8.3	7.2	6.1 ±1.2
	O2	13.4	13.4	18.4	16.6	14.5 ±2.4
	1

Table 1 : List of compounds present in the mass spectrum with their relative concentration

 1

Table SI 1

 SI : Fragmentation database used to construct the synthetic mass spectrum Table SI 2: Synthetic mass spectrum created to test the algorithm

	6
	4	m/z	Intensity (a.u.)	m/z		7 Intensity (a.u.)
		1		0			21		0		m/z
		12 13 14 2	0	15	16	22	17	18 0	19 20 24 25	26	27	28	29	30 31	32	36 37 38 39 40
	Methane 2.6 9.1 15.3 100.0 64.6 1.4 3 0 23	0
	Water	4		0		0.6 13.7 100.0 0.1 0.2 24 18.79
	Nitrogen	5	8.7	0			25		87.16				100.0 0.5
	Ethylene 0.5 1.1 2.8 6	0	0.4		26		603.7	2.7 11.7 80.0 55.7 100.0 2.0 0.1
	Ethane 0.5 1.6 3.3 7	5.3 0.01	0.2 27		441.25	0.6 5.6 43.1 42.7 100.0 21.8 26.2 0.8
	Argon	8		0			28		1156.17	12.9						0.2	0.0	100.0
	Oxygene	9		0		15.0 29		50.57							100.0
		10		0			30		43.81
		11		0			31		1.38
		12		28.66		32		635.05
		13		95.35		33		0
		14		196.37		34		0
		15		954.51		35		0
		16		711.12		36		0.6
		17		150.87		37		0
		18	1007.72		38		0.1
		19		0.66		39		0
		20		35.99		40		264.81

ACKNOWLEDGEMENTS

This research was supported by the National Aeronautics and Space Administration under16-CDAP16_2-0087, initially issued through the Cassini Data Analysis Program, and NASA's Planetary Science Division Internal Scientist Funding Program through the Fundamental Laboratory Research (FLaRe) work package to M.G.T.


Thomas R Cameron

Stef Graillat
email: stef.graillat@sorbonne-universite.fr

ON A COMPENSATED EHRLICH-ABERTH METHOD FOR THE ACCURATE COMPUTATION OF ALL POLYNOMIAL ROOTS *

Keywords: polynomial evaluation, error-free transformations, polynomial roots, backward error, forward error, rounding error analysis AMS subject classifications. 65H04, 65Y20, 65-04

In this article, we use the complex compensated Horner's method to derive a compensated Ehrlich-Aberth method for the accurate computation of all roots of a polynomial. In particular, under suitable conditions, we prove that the limiting accuracy for the compensated Ehrlich-Aberth iterations is as accurate as if computed in twice the working precision and then rounded into the working precision. Moreover, we derive a running error bound for the complex compensated Horner's and use it to form robust stopping criteria for the compensated Ehrlich-Aberth iterations. Finally, extensive numerical experiments illustrate that the backward and forward errors of the root approximations computed via the compensated Ehrlich-Aberth method are similar to those obtained with a quadruple precision implementation of the Ehrlich-Aberth method with a significant speed-up in terms of the computation time.

1. Introduction. The use of error-free transformations to produce compensated arithmetic routines has a long and interesting history, which includes the works of Dekker, Gill, Goldberg, Kahan, Knuth, and Møller [START_REF] Dekker | A floating-point technique for extending the available precision[END_REF][START_REF] Gill | A process for the step-by-step integration of differential equations in an automatic digital computing machine[END_REF][START_REF] Goldberg | What every computer scientist should know about floating-point arithmetic[END_REF][START_REF] Kahan | Further remarks on reducing truncation errors[END_REF][START_REF] Knuth | The Art of Computer Programming: Seminumerical Algorithms[END_REF][START_REF] Møller | Quasi double-precision in floating point addition[END_REF]. These works were the first to extend the working precision of a computation without the use of a hardware or software implementation of a high precision format [START_REF] Enge | MPC: A library for multiprecision complex arithmetic with exact rounding[END_REF][START_REF] Fousse | MPFR: A multipleprecision binary floating-point library with correct rounding[END_REF]. More recently, Rump, Ogita, and Oishi have developed algorithms for the summation and dot product computed in k-fold working precision [START_REF] Ogita | Accurate sum and dot product[END_REF], k ≥ 2, as well as the computation of a summation that is faithfully rounded [START_REF] Rump | Ultimately fast accurate summation[END_REF][START_REF] Rump | Accurate floating-point summation part I: faithful rounding[END_REF].

In addition, error-free transformations have been used to develop the compensated Horner's method for the evaluation of a polynomial and its derivatives at a real or complex number [START_REF] Graillat | Compensated Horner scheme[END_REF][START_REF] Graillat | Algorithms for accurate, validated and fast polynomial evaluation[END_REF][START_REF] Graillat | Accurate summation, dot product and polynomial evaluation in complex floating point arithmetic[END_REF][START_REF] Jiang | Accurate evaluation of a polynomial and its derivative in Bernstein form[END_REF][START_REF] Jiang | Accurate evaluation of the k-th derivative of a polynomial and its application[END_REF]. The real compensated Horner's method has been used to improve the accuracy of eigenvalue approximations of a symmetric tridiagonal matrix and real root approximations of a real polynomial [START_REF] Gemignani | Accurate polynomial root-finding methods for symmetric tridiagonal matrix eigenproblems[END_REF][START_REF] Graillat | Accurate simple zeros of polynomials in floating point arithmetic[END_REF]. In [START_REF] Graillat | Accurate simple zeros of polynomials in floating point arithmetic[END_REF], Graillat shows that if the real compensated Horner's method is used to evaluate a polynomial, then the Newton iterations converge to an approximate root as accurately as if computed in twice the working precision and then rounded into the working precision.

The Ehrlich-Aberth (Börsch-Supan) method [START_REF] Aberth | Iteration methods for finding all zeros of a polynomial simultaneously[END_REF][START_REF] Börsch-Supan | A posterior error bounds for the zeros of polynomials[END_REF][START_REF] Ehrlich | A modified Newton method for polynomials[END_REF] combines Newton's method with an implicit deflation strategy, which allows for the computation of all roots of a polynomial. In this article, we use the complex compensated Horner's method to develop a compensated Ehrlich-Aberth method that can compute all roots of a polynomial as accurately as the Newton iterations from [START_REF] Graillat | Accurate simple zeros of polynomials in floating point arithmetic[END_REF]. The outline of this article is as follows: In Section 2, we recall the basic properties of real and complex floating-point arithmetic. Then, in Section 3, we describe the Horner method and compensated Horner method for polynomial evaluation in complex floating-point arithmetic. Moreover, we derive a running error bound for the compensated Horner method, which considers the rounding errors that occur during the computation. Section 4 is devoted to the presentation of the Ehrlich-Aberth method and compensated Ehrlich-Aberth method. In particular, the running error bound for the compensated Horner method is used to form robust stopping criteria for the sequence of root approximations. Moreover, under suitable conditions, we prove that these root approximations have a limiting accuracy as if computed in twice the working precision and then rounded into the working precision. Finally, in Section 5, extensive numerical experiments illustrate that the backward and forward errors of the root approximations computed via the compensated Ehrlich-Aberth are similar to those obtained with a quadruple precision implementation of the Ehrlich-Aberth method but with a significant speed-up in terms of the computation time.

2. Floating-Point Arithmetic. Throughout this article, we assume that the computer arithmetic satisfies the IEEE 754 standard [START_REF]ANSI/IEEE, IEEE Standard for Binary Floating Point Arithmetic[END_REF], and that no underflow nor overflow occurs. We denote by F the set of floating-point numbers and by µ the unit roundoff. Note that for single precision, µ = 2 -24 and for double precision, µ = 2 -53 , where the exponent corresponds to the precision of this floating-point format. Finally, we use the standard notation fl (•) to denote floating-point operations in working precision.

2.1. Real Floating-Point Arithmetic. For operations • ∈ {+, -, •}, the IEEE 754 standard requires the result of fl (a • b) to be correctly rounded, i.e., as accurate as if computed exactly and then rounded to the current precision [START_REF] Goldberg | What every computer scientist should know about floating-point arithmetic[END_REF]. In this article, we assume that all the computations are performed with rounding to nearest (using round to even). Throughout this article, we make use of the quantities γ n , which are defined in the usual way [START_REF] Higham | Accuracy and Stability of Numerical Algorithms[END_REF]:

γ n = nµ 1 -nµ ,
where n ∈ N is assumed to satisfy nµ < 1. Moreover, we make use of the round to nearest mode to perform error-free transformations for a floating-point operation. In particular, for each x = fl (a • b), there exists a y ∈ F such that x + y = a • b. The pair (x, y) is called the error-free transformation of (a, b) for the operation •. For instance, Algorithm 2.1 is attributed to Knuth [START_REF] Knuth | The Art of Computer Programming: Seminumerical Algorithms[END_REF] and returns the error-free transformation of (a, b) for addition. Note that Algorithm 2.1 requires 6 flops to be executed. In addition, we make use of the fused multiply-add operation, which is denoted FMA (a, b, c) and results in the floating point number nearest to a•b+c ∈ R. Note that the fused multiply-add operation was added to the IEEE 754 standard in 2008 and is supported by many modern processors [START_REF] Nievergelt | Scalar fused multiply-add instructions produce floating-point matrix arithmetic provably accurate to the penultimate digit[END_REF]. Moreover, the error-free transformation of (a, b) for multiplication can be performed using the fused multiply-add operation as is done in [START_REF] Ogita | Accurate sum and dot product[END_REF], see Algorithm 2.2. For processors that don't support the fused multiply-add operation, the error-free transformation of (a, b) for multiplication can be computed with the splitting operation introduced by Dekker [START_REF] Dekker | A floating-point technique for extending the available precision[END_REF] but requires 17 flops rather than the 2 flops required by Algorithm 2.2, see Algorithm 3.3 of [START_REF] Ogita | Accurate sum and dot product[END_REF]. We conclude this section with a theorem from [START_REF] Ogita | Accurate sum and dot product[END_REF] that summarizes the properties of Algorithm 2.1 and Algorithm 2.2.

a • b = x + y, x = fl (a • b) , |y| ≤ µ |x| , |y| ≤ µ |a • b| .

Complex

Floating-Point Arithmetic. We define C = F + iF to be the set of complex floating-point numbers, where i = √ -1 is the imaginary unit. Also, we use the operators Re (•) and Im (•) to denote the real and imaginary part of a complex number, respectively. As in the real case, we denote by fl (•) the operations that are done in floating-point working precision. The following holds for all a, b ∈ C and • ∈ {+, -}:

fl (a • b) = (a • b)(1 + ),
where | | ≤ µ. In addition, we have

fl (a • b) = (a • b)(1 + ), where | | ≤ √ 2γ 2 . This further implies that (2.2) |fl (a • b) -a • b| ≤ µ |a • b| and |a • b -fl (a • b)| ≤ µ |fl (a • b)| , for • ∈ {+, -} and (2.3) |a • b -fl (a • b)| ≤ √ 2γ 2 |a • b| .
Finally, throughout this article, we make use of the quantity

γ n = n √ 2γ 2 1 -n √ 2γ 2 .
As in the real case, the error-free transformation of the pair of complex floatingpoint numbers (a, b) for the operation • is a pair (x, y) such that x = fl (a • b) and x + y = a • b. The error-free transformation of (a, b) ∈ C 2 for complex addition is a [g 5 , h 5 ] = TwoSum (g 1 , -g 2 ); [g 6 , h 6 ] = TwoSum (g 3 , g 4 )

w = g 5 + ig 6 ; x = h 1 + ih 3 ; y = -h 2 + ih 4 ; z = h 5 + ih 6
We conclude this section with a theorem from [START_REF] Graillat | Accurate summation, dot product and polynomial evaluation in complex floating point arithmetic[END_REF] that summarizes the properties of Algorithm 2.3 and Algorithm 2.4. 3. Horner's Method. Consider the polynomial of degree m in the variable z defined by

(3.1) p(z) = a m z m + • • • + a 1 z + a 0 ,
where a 0 , a 1 , . . . , a m ∈ C, and a m = 0. Given z ∈ C, we can compute the polynomial evaluation p(z) using Horner's method as shown in Algorithm 3.1.

Algorithm 3.1 Horner's Method. function [h 0 ] = Horner (p, z) : h m = a m for k = m -1 to k = 0 do h k = fl (z • h k+1 + a k ) end for
The result h 0 = Horner (p, z) satisfies the following forward error bound [START_REF] Graillat | Accurate summation, dot product and polynomial evaluation in complex floating point arithmetic[END_REF]:

|p(z) -h 0 | ≤ γ 2m p(|z|), where p(z) = m k=0 |a k | z k . Finally, if p(z) = 0, then we have the relative forward error bound |p(z) -h 0 | |p(z)| ≤ γ 2m cond (p, z) ,
where the condition number of the polynomial evaluation of p at z is defined by

(3.2) cond (p, z) = p(|z|)/ |p(z)| .
3.1. Compensated Horner's Method. To improve the standard error bound of Horner's method, we follow the development in [START_REF] Graillat | Accurate summation, dot product and polynomial evaluation in complex floating point arithmetic[END_REF] to produce a compensated Horner's method in complex floating-point arithmetic. To this end, we record the error at each iteration from both the floating-point product and sum operations. Specifically, we produce four error polynomials: p π , p µ , p ν , p σ , which are collected monomial-by-monomial as outlined in Algorithm 3. To obtain a true error-free transformation, we would need to recursively perform Algorithm 3.2 on the four error polynomials of degree (m -1) until the resulting error polynomials were constant. This recursive process would result in many error terms: Four polynomials of degree (m -1), 16 polynomials of degree (m -2), and so on to include 4 m constant polynomials. Moreover, most of the coefficients in these polynomials would suffer from underflow. For this reason, we only consider the first-order error terms. In particular, Algorithm 3.3 shows how we compute the forward error e(z) = (p π + p µ + p ν + p σ )(z) using Horner's method applied to the polynomial whose coefficients are those of (p π + p µ + p ν + p σ ), computed using the doubly compensated summation method from Priest [START_REF] Priest | On properties of floating point arithmetics: Numerical stability and the cost of accurate computations[END_REF].

h m = a m for k = m -1 to k = 0 do [ ĥk , π k , µ k , ν k ] = TwoProductCmplx (h k+1 , z) [h k , σ k ] = TwoSumCmplx ĥk , a k Set π k , µ k , ν k , σ k ,
In fact, any of the summation methods from [START_REF] Demmel | Accurate and efficient floating point summation[END_REF][START_REF] Higham | Accuracy and Stability of Numerical Algorithms[END_REF][START_REF] Priest | On properties of floating point arithmetics: Numerical stability and the cost of accurate computations[END_REF][START_REF] Rump | Ultimately fast accurate summation[END_REF][START_REF] Rump | Accurate floating-point summation part I: faithful rounding[END_REF] that guarantee a relative error bound of 2µ can be used. For example, in [START_REF] Graillat | Accurate summation, dot product and polynomial evaluation in complex floating point arithmetic[END_REF], the authors advocate for the use of the accurate summation method from [START_REF] Rump | Ultimately fast accurate summation[END_REF][START_REF] Rump | Accurate floating-point summation part I: faithful rounding[END_REF], the result of which is guaranteed to be faithfully rounded, i.e., as accurate as if computed exactly and then rounded to an adjacent floating-point number. Hence, the relative error in the accurate summation method is 2µ; however, since we only need to sum 4 complex floating-point numbers, it is more efficient to use the doubly compensated summation method. Although the result of this method is not guaranteed to be faithfully rounded, it follows from the analysis of Priest in [START_REF] Priest | On properties of floating point arithmetics: Numerical stability and the cost of accurate computations[END_REF]Section 4.1] that it has the same relative error of 2µ. If we ignore the cost of the sort, which for the given 4 floating-point values requires 9 absolute values, 6 comparisons, and 3 swaps, then the execution of Algorithm 3.4 requires 30 flops. Since Algorithm 3.3 uses the complex doubly compensated summation from Algorithm 3.5, it follows that its execution requires 68m -8 flops.

Algorithm 3.3 Horner's method applied to the degree (m -1) polynomial (p + q + r + s).

function

[v 0 ] = HornerSum (p, q, r, s, z) : v m-1 = DbleCompSumCmplx (p m-1 , q m-1 , r m-1 , s m-1 ) for k = m -2 to k = 0 do v k = fl (z • v k+1 + DbleCompSumCmplx (p k , q k , r k , s k )) end for Algorithm 3.4 Doubly compensated summation of a 1 , a 2 , a 3 , a 4 ∈ F. function [s n ] = DbleCompSum (a 1 , a 2 , a 3 , a 4 ) : Sort the a i so that |a 1 | ≥ • • • ≥ |a 4 | s 1 = a 1 , c 1 = 0 for k = 2 to k = 4 do y k = fl (c k-1 + a k ) u k = fl (a k -(y k -c k-1 )) t k = fl (y k + s k-1 ) v k = fl (y k -(t k -s k-1 )) z k = fl (u k + v k ) s k = fl (t k + z k ) c k = fl (z k -(s k -t k )) end for Algorithm 3.5 Doubly compensated summation of a 1 , a 2 , a 3 , a 4 ∈ C. function [x + iy] = DbleCompSumCmplx (a 1 , a 2 , a 3 , a 4 ) : x = DbleCompSum (Re (a 1 ) , Re (a 2 ) , Re (a 3 ) , Re (a 4 )) y = DbleCompSum (Im (a 1 ) , Im (a 2 ) , Im (a 3 ) , Im (a 4 ))
Let p π , p µ , p ν , p σ be the error polynomials from Algorithm 3.2. Then, by [17, Lemma 5.4], the result v 0 = HornerSum (p π , p µ , p ν , p σ , z) satisfies the following forward error bound:

(3.3) |e(z) -v 0 | ≤ γ 2m-1 ( (p π + p µ + p ν ) + p σ )(|z|).
It is worth noting that [START_REF] Graillat | Accurate summation, dot product and polynomial evaluation in complex floating point arithmetic[END_REF]Lemma 5.4] is formally stated for the accurate summation method from [START_REF] Rump | Ultimately fast accurate summation[END_REF][START_REF] Rump | Accurate floating-point summation part I: faithful rounding[END_REF]; however, since the doubly compensated method has the same relative error bound the result follows. Now, by combining the error-free transformation of Horner's method and the method for computing the corresponding forward error, we obtain the compensated Horner's method shown in Algorithm 3.6, which requires 100m -7 flops to be executed.

Algorithm 3.6 Compensated Horner's Method. function [c 0 ] = CompHorner (p, z) : [h 0 , p π , p µ , p ν , p σ ] = EFTHorner (p, z) e 0 = HornerSum (p π , p µ , p ν , p σ , z) c 0 = fl (h 0 + e 0 )
From [17, Theorem 5.5], we know that the result c 0 = CompHorner (p, z) satisfies the following forward error bound:

(3.4) |p(z) -c 0 | ≤ µ |p(z)| + γ 2 2m p(|z|).
Furthermore, if p(z) = 0, then we have the following relative forward error bound

|p(z) -c 0 | |p(z)| ≤ µ + γ 2 2m cond (p, z) .
Therefore, the compensated Horner's method is as accurate as if computed in twice the working precision and then rounded into the working precision.

Running Error Bound.

The error bound in (3.4) is not useful in practice since it contains the quantity p(z). Therefore, we establish a running error bound that can be computed using the error polynomials produced in Algorithm 3.2. We begin with the following lemma. Lemma 3.1. Let p, q, r, s be degree (m -1) polynomials with non-negative coefficients a k , b k , c k , d k ∈ F, and let x ∈ F be non-negative. Then,

0 ≤ m-1 k=0 (a k + b k + c k + d k )x k ≤ (1 + 2µ) m-1 HornerSum (p, q, r, s, x) .
Proof. Consider Algorithm 3.3 and the intermediate values v k . We prove that for k = 0, 1, . . . , m -1, we have

s(k) := k j=0 (a m-1-k+j + b m-1-k+j + c m-1-k+j + d m-1-k+j )x j ≤ (1 + 2µ) 2k+1 v m-1-k . (3.5)
The base case, k = 0, holds since the relative error in the doubly compensated summation is bounded above by 2µ, that is,

(a m-1 + b m-1 + c m-1 + d m-1 ) ≤ (1 + 2µ) DbleCompSum (a m-1 , b m-1 , c m-1 , d m-1 ) = (1 + 2µ)v m-1 .
Suppose that (3.5) holds for some integer k, where 0 ≤ k < m -1, and note that s(k + 1) can be split into the sum of two parts:

  k j=0 (a m-1-k+j + b m-1-k+j + c m-1-k+j + d m-1-k+j )x j   x and a m-k-2 + b m-k-2 + c m-k-2 + d m-k-2 .
By the induction hypothesis, the first part is bounded above by (1+2µ) 2k+1 x•v m-1-k and, since the relative error in the doubly compensated summation is bounded above by 2µ, the second part is bounded above by

(1 + 2µ) DbleCompSum (a m-k-2 , b m-k-2 , c m-k-2 , d m-k-2 ) .
Therefore, s(k + 1) is bounded above by

(1 + 2µ) 2k+1 x • v m-1-k + (1 + 2µ) DbleCompSum (a m-k-2 , b m-k-2 , c m-k-2 , d m-k-2 ) , which, in turn, is bounded above by (1 + 2µ) 2k+1 (1 + µ) 2 times fl (x • v m-1-k + DbleCompSum (a m-k-2 , b m-k-2 , c m-k-2 , d m-k-2 )) .
Thus,

s(k + 1) ≤ (1 + 2µ) 2(k+1)+1 v m-k-2 ,
and it follows that (3.5) holds for k = 0, 1, . . . , m -1.

We are now ready to prove the running error bound for the compensated Horner's method in Algorithm 3.6. Theorem 3.2. Let p be a complex polynomial as defined in (3.1). Then, c 0 = CompHorner (p, z) satisfies

|c 0 -p(z)| ≤ fl µ |c 0 | + ( γ 4m+2 HornerSum (|p π | , |p µ | , |p ν | , |p σ | , |z|) + 2µ 2 |c 0 |) . Proof. Note that c 0 = CompHorner (p, z) satisfies |c 0 -p(z)| = |fl (h 0 + e 0 ) -p(z)| ≤ |fl (h 0 + e 0 ) -(h 0 + e 0 )| + |(h 0 + e 0 ) -p(z)| ≤ µ |fl (h 0 + e 0 )| + |e 0 -e(z)| ,
where e(z) = (p π + p µ + p ν + p σ )(z). Applying the forward error bound in (3.3) and Lemma 3.1, we have

|e 0 -e(z)| ≤ γ 2m-1 ( (p π + p µ + p ν ) + p σ )(|z|) ≤ (1 + 2µ) 2m-1 γ 2m-1 HornerSum (|p π | , |p µ | , |p ν | , |p σ | , |z|) .
Next, we make use of the following inequalities:

(1 + √ 2γ 2 ) γ m ≤ γ m+1 and γ m ≤ (1 - √ 2γ 2 ) γ m+1 .
The first inequality implies that (1 9 γ 4m+2 . We summarize these results in the following inequality

+ √ 2γ 2 ) 2m-1 γ 2m-1 ≤ γ 4m-2 . This combined with 2µ ≤ √ 2γ 2 implies that (1 + 2µ) 2m-1 γ 2m-1 ≤ γ 4m-2 . Furthermore, the second inequality implies that γ 4m-2 ≤ (1 - √ 2γ 2 ) 4 γ 4m+2 . Since (1 - √ 2γ 2 ) ≤ (1 -µ) 9/4 , it follows that γ 4m-2 ≤ (1 -µ)
(3.6) (1 + 2µ) 2m-1 γ 2m-1 ≤ γ 4m-2 ≤ (1 -µ) 9 γ 4m+2 .
Now, we analyze the error in the computation of

γ k = k √ 2γ 2 1 -k √ 2γ 2 = 2kµ √ 2 (1 -2µ) -2kµ √ 2 .
Note that 2µ, 2kµ, and (1 -2µ) are floating-point numbers computed exactly, and the IEEE standard requires the square-root operation be rounded to nearest. Therefore, 2kµ

√ 2 ≤ (1-µ) -2 fl 2kµ √ 2 and (1-2µ)-2kµ √ 2 ≥ (1-µ) 3 fl (1 -2µ) -2kµ √ 2
, and it follows that

(3.7) γ k ≤ (1 -µ) -6 fl ( γ k ) .
Finally, we apply the inequalities in (3.6) and (3.7) to obtain the following:

|c 0 -p(z)| ≤ µ |c 0 | + (1 + 2µ) 2m-1 γ 2m-1 HornerSum (|p π | , |p µ | , |p ν | , |p σ | , |z|) ≤ µ |c 0 | + (1 -µ) 9 γ 4m+2 HornerSum (|p π | , |p µ | , |p ν | , |p σ | , |z|) ≤ µ |c 0 | + (1 -µ) 2 fl ( γ 4m+2 HornerSum (|p π | , |p µ | , |p ν | , |p σ | , |z|)) .
Therefore, the error in the computation c 0 is bounded above by

(1 -µ)µ |c 0 | + (1 -µ) 2 fl ( γ 4m+2 HornerSum (|p π | , |p µ | , |p ν | , |p σ | , |z|)) + µ 2 |c 0 | ,
which, since we can always assume that 2(1 -µ) 2 ≥ 1, is bounded above by

(1 -µ)µ |c 0 | + (1 -µ) 2 fl ( γ 4m+2 HornerSum (|p π | , |p µ | , |p ν | , |p σ | , |z|)) + 2µ 2 |c 0 | .
Furthermore, since we assume that no underflow occurs, both µ |c 0 | and 2µ 2 |c 0 | are floating-point numbers computed exactly and it follows that the error in the computation c 0 is bounded above by

(1 -µ) µ |c 0 | + fl γ 4m+2 HornerSum (|p π | , |p µ | , |p ν | , |p σ | , |z|) + 2µ 2 |c 0 | .
Therefore, we have

|c 0 -p(z)| ≤ fl µ |c 0 | + ( γ 4m+2 HornerSum |p π | , |p µ | , |p ν | , |p σ | , |z|) + 2µ 2 |c 0 | .
4. Ehrlich-Aberth Method. Let p be a degree m complex polynomial as defined in (3.1). Throughout this section, we assume that p has simple roots ζ 1 , . . . , ζ m ∈ C. For 1 ≤ i ≤ m and n ≥ 0, let z n,i ∈ C denote the approximation of ζ i after n iterations. The Ehrlich-Aberth method [START_REF] Aberth | Iteration methods for finding all zeros of a polynomial simultaneously[END_REF][START_REF] Bini | Numerical computation of polynomial zeros by means of Aberth's method[END_REF][START_REF] Börsch-Supan | A posterior error bounds for the zeros of polynomials[END_REF][START_REF] Ehrlich | A modified Newton method for polynomials[END_REF] updates each root approximation as follows (4.1)

z n+1,i = fl z n,i - p(z n,i ) p (z n,i ) -p(z n,i )A n,i (z n,i ) , where A n,i (z) = j =i (z -z n,j ) -1
. Note that we reference the fractional expression on the right side of (4.1) as the Ehrlich-Aberth correction term. The polynomial evaluations in (4.1) can be computed using Horner's method from Algorithm 3.1. Moreover, given z ∈ C and h 0 = Horner (p, z), we can push the error in the computation back onto the coefficients as follows [START_REF] Bini | Numerical computation of polynomial zeros by means of Aberth's method[END_REF]Theorem 7]:

h 0 = m i=0 a i (1 + i ) z i , where | i | ≤ (2 √ 2 + 1)i + 1 µ + O(µ 2
). Note that the backward error of z as a root of the polynomial p is defined by

η(z) = min { : (p + ∆p) (z) = 0, |∆a i | ≤ |e i | , i = 0, 1, . . . , m} ,
where the e i are arbitrary and represent tolerances against which perturbations are measured, and ∆p(z

) = ∆a m z m + • • • + ∆a 1 z + ∆a 0 .
It is well-known, e.g., see [START_REF] Chaitin-Chatelin | Lectures on Finite Precision Computations, Software, Environments and Tools[END_REF], that

η(z) = |p(z)| α(|z|) ,
where

α(z) = m i=0 |e i | z i . Let e i = (( 2
√ 2 + 1)i + 1)a i . Then, the inequality (4.2) η(z) ≤ µ
guarantees that z is a root of (p + ∆p), where |∆a i | is no bigger than the upper bound of |a i | i , for all i = 0, 1, . . . , m. Furthermore, if η(z) > µ, then z being a root of (p + ∆p) implies that |∆a i | is bigger than |a i | i for some i = 0, 1, . . . , m. Hence, in exact arithmetic, (4.2) denotes a stopping criterion that guarantees iterations do no terminate until z is a root of a polynomial whose coefficients are no more perturbed than the floating-point computation of p(z). In practice, we stop updating the root approximation z if the following inequality holds:

(4.3) |h 0 | > µ Horner (α, |z|) ,
where the coefficients of α have been stored as floating-point numbers.

4.1. Compensated Ehrlich-Aberth Method. Let z ∈ C be a root approximation of the polynomial p such that the stopping criterion (4.3) holds. Then, h 0 = Horner (p, z) is no longer a reliable computation for updating the root approximation z. However, c 0 = CompHorner (p, z) may still be reliable, depending on the value of cond (p, z) as defined in (3.2). In this case, we update the root approximation using (4.1) and the compensated Horner's method. If the polynomial has well-conditioned simple roots, then we only need to apply the compensated Horner's method to p. However, the more difficult polynomial equations have multiple or near multiple roots, which may be ill-conditioned. For this reason, we use the compensated Horner's method to evaluate p and p . A method for performing these evaluations is derived from [19, Algorithm 2] and presented in Algorithm 4.1.

It is clear that the value of c 0 from Algorithm 4.1 is the same as its value from Algorithm 3.6; therefore, c 0 satisfies the forward error bound in (3.4) and the running error bound in Theorem 3.2. Moreover, the derivative value c 0 can be shown to satisfy a similar forward error bound by appealing to [START_REF] Jiang | Accurate evaluation of the k-th derivative of a polynomial and its application[END_REF]Theorem 4] and making the necessary changes from real to complex floating-point arithmetic.

Furthermore, the value of s 0 from Algorithm 4.1 satisfies

s 0 = HornerSum (|p π | , |p µ | , |p ν | , |p σ | , |z|) ,
where p π , p µ , p ν , and p σ are the polynomials created by EFTHorner (p, z). Therefore, by Theorem 3.2, it follows that

|c 0 -p(z)| ≤ fl µ |c 0 | + (γ 4m+2 s 0 + 2µ 2 |c 0 |) . Hence, if fl µ |c 0 | + (γ 4m+2 s 0 + 2µ 2 |c 0 |) < |c 0 | ,
s m = 0 for k = m -1 to k = 0 do [ ĥ k , π k , µ k , ν k ] = TwoProductCmplx h k+1 , z [h k , σ k ] = TwoSumCmplx ĥ k , h k+1 e k = fl z • e k+1 + e k+1 + DbleCompSumCmplx (π k , µ k , ν k , σ k ) [ ĥk , π k , µ k , ν k ] = TwoProductCmplx (h k+1 , z) [h k , σ k ] = TwoSumCmplx ĥk , a k e k = fl (z • e k+1 + DbleCompSumCmplx (π k , µ k , ν k , σ k )) s k = fl (|z| • s k+1 + DbleCompSum (|π k | , |µ k | , |ν k | , |σ k |)) end for c 0 = fl (h 0 + e 0 ), c 0 = fl (h 0 + e 0 )
then the error in the computation of p(z) is smaller than the size of |c 0 | and is still a reliable computation. Thus, we obtain the following stopping criterion for the compensated Ehrlich-Aberth method (4.4)

|c 0 | ≤ fl µ |c 0 | + (γ 4m+2 s 0 + 2µ 2 |c 0 |) .
Finally, let z ∈ C be an approximate root of the polynomial p, and consider the computation [c 0 , c 0 , s 0 ] = CompHornerDer (p, z). Denote by EA(p, z) the computed Ehrlich-Aberth correction term using the values of c 0 and c 0 . Then, even if (4.4) does not hold, it may happen that the following does hold:

(4.5) |EA(p, z)| ≤ µ |z| ,
which implies that the relative change made to z by the Ehrlich-Aberth correction term is insignificant. In summary, we don't stop updating the root approximations from the compensated Ehrlich-Aberth method until either (4.4) or (4.5) holds, which guarantees that the updates don't cease until the relative error in the compensated Horner method is too big, or the relative size the Ehrlich-Aberth correct term is too small. Before proceeding, note that Horner and compensated Horner methods are prone to overflow; for instance, when a large degree polynomial with positive coefficients is evaluated at z, where |z| > 1. For this reason, the reversal polynomial p R (z) := z m p(1/z), which for ρ = 1/z satisfies

p(z) = z m p R (ρ) p (z) = mz m-1 p R (ρ) -z m-2 p R (ρ),
is used to compute the Ehrlich-Aberth correction term when |z| > 1 [START_REF] Bini | Numerical computation of polynomial zeros by means of Aberth's method[END_REF].

Limiting Accuracy.

In [START_REF] Tisseur | Newton's method in floating point arithmetic and iterative refinement of generalized eigenvalue problems[END_REF], Tisseur examined the multi-variable Newton's method, allowing for extended precision in the computation of the residual, and its application to iterative refinement for the generalized eigenvalue problem. In particular, for a close enough eigenvalue approximation, Newton's method converges to an eigenvalue approximation with limiting accuracy that depends on the relative rounding error µ and the accuracy of the computed residual.

Then, in [START_REF] Graillat | Accurate simple zeros of polynomials in floating point arithmetic[END_REF], Graillat applied this result to real polynomials with real simple zeros. Specifically, if the polynomial evaluation is computed using the real compensated Horner's method, then for a close enough root approximation, Newton's method converges to a root approximation with limiting accuracy as if computed in twice the working precision and then rounded into the working precision. In this section, we show that the compensated Ehrlich-Aberth method has a similar limiting accuracy for all, real or complex, roots of a polynomial.

To this end, we re-write (4.1) as follows

(4.6) z n+1,i = z n,i -J n,i (z n,i ) + e n,i -1 (p(z n,i ) + e n,i ) + n,i ,
where J n,i (z) = p (z) -p(z)A n,i (z), e n,i is the error in computing p(z n,i ), e n,i is the error in computing J n,i (z n,i ) and performing the division, and n,i is the error in the subtraction.

For i = 1, . . . , m, we assume that p(z n,i ) is computed via the compensated Horner's method in Algorithm 3.6. Hence,

|e n,i | ≤ µ |p(z n,i )| + γ 2 2m p(|z n,i |).
We also assume that e n,i ≤ µφ (p, z n,i , m, µ) ,

for some function φ that reflects the error in computing J n,i (z n,i ) and performing the division. For the error n,i , we have

| n,i | ≤ µ |z n,i | + J n,i (z n,i ) + e n,i -1 (p(z n,i ) + e n,i )
Since the roots of p, denoted by ζ 1 , . . . , ζ m , are assumed to be simple, it follows that for i = 1, . . . , m, there exists a closed disk D i in the complex plane centered at ζ i such that ζ j / ∈ D i for all j = i. Furthermore, for close enough root approximations, we can assume that z n,i ∈ D i and z n,j / ∈ D i for all j = i. Under these assumptions, we have the following result. Lemma 4.1. For i = 1, . . . , m, there exists a β i > 0 such that for all z, w ∈ D i , we have

|p(z) -p(w) -J n,i (w)(z -w)| ≤ β i 2 |z -w| 2 + |p(w)| |A n,i (w)| |z -w| .
Proof. Since p is Lipschitz continuous on D i , there exists a

β i > 0 such that |p (z) -p (w)| ≤ β i |z -w| ,
for all z, w ∈ D i . Furthermore, for all z, w ∈ D i , p(z) -p(w) -J n,i (w)(zw) can be written as

1 0 (p (w + t(z -w)) -p (w)) (z -w)dt + p(w)A n,i (w)(z -w).
Therefore, by Holder's inequality and the Lipschitz continuity of p , we have

|p(z) -p(w) -J n,i (w)(z -w)| ≤ β i 2 |z -w| 2 + |p(w)| |A n,i (w)| |z -w| .
From Lemma 4.1 and the continuity of p and A n,i on D i , it follows that there exists M i > 0 such that

|p(z) -p(w) -J n,i (w)(z -w)| ≤ β i 2 |z -w| 2 + M i |z -w| ,
for all z, w ∈ D i . Furthermore, the Lipschitz continuity of J n,i on D i implies that there exists β i > 0 such that

|J n,i (z) -J n,i (w)| ≤ β i |z -w| ,
for all z, w ∈ D i .

Theorem 4.2. For i = 1, . . . , m assume that

(4.7) J n,i (z n,i ) -1 e n,i ≤ ν < 1, and
(4.8) β i J n,i (ζ i ) -1 |z n,i -ζ i | ≤ τ < 1.
Then, z n+1,i in (4.6) is well-defined and satisfies

|z n+1,i -ζ i | ≤ G i |z n,i -ζ i | + g i ,
where

G i = ν 1 -ν + µ(2 + µ) (1 -τ )(1 -ν) + M i + µ + (1 + µ) 2 τ 2(1 -τ )(1 -ν) and g i = µ |ζ i | + γ 2 2m 1 + µ (1 -τ )(1 -ν) p(|z n,i |) |p (ζ i )| .
Proof. By (4.8) and the Lipschitz continuity of J n,i on D i , we have (4.9)

J n,i (ζ i ) -1 |J n,i (z n,i ) -J n,i (ζ i )| ≤ β i J n,i (ζ i ) -1 |z n,i -ζ i | ≤ τ < 1.
Also, note that

J n,i (z n,i ) -1 = 1 + J n,i (z n,i ) -J n,i (ζ i ) J n,i (ζ i ) -1 J n,i (ζ i ) -1 ,
which, combined with (4.9), implies that (4.10)

J n,i (z n,i ) -1 ≤ J n,i (ζ i ) -1 1 -|J n,i (ζ i ) -1 | |J n,i (z n,i ) -J n,i (ζ i )| ≤ J n,i (ζ i ) -1 1 -τ .
Similarly, (4.7) and (4.10) give us

(4.11) J n,i (z n,i ) + e n,i -1 ≤ J n,i (z n,i ) -1 1 -J n,i (z n,i ) -1 e n,i ≤ J n,i (ζ i ) -1 (1 -τ )(1 -ν) .
Therefore, J n,i (z n,i ) + e n,i is non-zero, so z n+1,i in (4.6) is well-defined and

z n+1,i -ζ i = z n,i -ζ i -J n,i (z n,i ) + e n,i -1 (p(z n,i ) + e n,i ) + n,i = 1 - J n,i (z n,i ) J n,i (z n,i ) + e n,i (z n,i -ζ i ) - p(z n,i ) -J n,i (z n,i ) (z n,i -ζ i ) + e n,i J n,i (z n,i ) + e n,i + n,i .
Hence, we have the following bound

|z n+1,i -ζ i | ≤ 1 - J n,i (z n,i ) J n,i (z n,i ) + e n,i |z n,i -ζ i | + |p(z n,i ) -J n,i (z n,i ) (z n,i -ζ i )| + |e n,i | J n,i (z n,i ) + e n,i + | n,i | . (4.12) From 1 - J n,i (z n,i ) J n,i (z n,i ) + e n,i = e n,i J n,i (z n,i ) + e n,i = J n,i (z n,i ) -1 e n,i 1 + J n,i (z n,i ) -1 e n,i
and (4.7), we have

(4.13) 1 - J n,i (z n,i ) J n,i (z n,i ) + e n,i ≤ ν 1 -ν .
From Lemma 4.1, we have

(4.14) |p(z n,i ) -J n,i (z n,i ) (z n,i -ζ i )| ≤ β i 2 |z n,i -ζ i | 2 + M i |z n,i -ζ i | and |p(z n,i ) -J n,i (ζ n,i ) (z n,i -ζ i )| ≤ β i 2 |z n,i -ζ i | 2 .
Hence,

|p(z n,i )| ≤ |p(z n,i ) -J n,i (ζ i ) (z n,i -ζ i )| + |J n,i (ζ i ) (z n,i -ζ i )| ≤ β i 2 |z n,i -ζ i | 2 + |J n,i (ζ i )| |z n,i -ζ i | , (4.15)
which allows us to bound the error term e n,i as follows (4.16) |e

n,i | ≤ µ β i 2 |z n,i -ζ i | 2 + |J n,i (ζ i )| |z n,i -ζ i | + γ 2 2m p(|z n,i |)
and the error term n,i as follows (4.17)

| n,i | ≤ µ (|z n,i -ζ i | + |ζ i | + |d n,i |) ,
where

|d n,i | = J n,i (z n,i ) + e n,i -1 (|p(z n,i )| + |e n,i |) ≤ J n,i (ζ i ) -1 (1 -τ )(1 -ν) (1 + µ) β i 2 |z n,i -ζ i | 2 + |J n,i (ζ i )| |z n,i -ζ i | + γ 2 2m p(|z n,i |) .
Note that conditions (4.7) and (4.8) are necessary for z n+1,i in (4.6) to be defined. Assumption (4.7) is a condition on the accuracy of computing J n,i (z) and assumption (4.8) is a condition on the diameter of the disk D i and the conditioning of the root ζ i . In our final result, we assume that we can make ν, τ and M i small enough, say ν ≤ 1 8 , τ ≤ 1 8 , and M i ≤ 1 8 . Under these assumptions, we have the following result, which implies that under suitable conditions the compensated Ehrlich-Aberth method computes all roots of a polynomial as accurate as if computed in twice the working precision and then rounded into the working precision.

Corollary 4.3. For i = 1, . . . , m, let z 0,i ∈ D i . Then, the compensated Ehrlich-Aberth method generates a sequence of of approximations whose absolute error decreases until the first n for which

(4.18) |z n+1,i -ζ i | ≈ µ |ζ i | + γ 2 2m p(|ζ i |) |p (ζ i )| .
Furthermore, if ζ i = 0, then we have the following relative limiting accuracy

(4.19) |z n+1,i -ζ i | |ζ i | ≈ µ + γ 2 2m cond (p, ζ i ) ,
where

(4.20) cond (p, ζ) = p(|ζ|) |ζ| |p (ζ)|
is the condition number for the computation of the root ζ of p.

Proof. By Theorem 4.2, for i = 1, . . . , m, we have

|z 1,i -ζ i | ≤ G i |z 0,i -ζ i | + g i .
Under the assumption that ν ≤ 1 8 , τ ≤ 1 8 , and M i ≤ 1 8 , we have G i ≤ 1 2 and it follows that the absolute error contracts unless (4.18) already holds. Thus, the result follows via induction.

Numerical Experiments.

In this section, we present the results of several numerical experiments to demonstrate the running error bound of the compensated Horner method in Theorem 3.2 and the limiting accuracy of the compensated Ehrlich-Aberth method in Corollary 4.3. In addition, we compare the computation time of the compensated Horner to the Horner method implemented in double and quadruple precision. Finally, we compare the computation time, backward error, and forward error of the compensated Ehrlich-Aberth method to the Ehrlich-Aberth method implemented in quadruple precision. Note that all higher precision computations, such as those in quadruple precision, are implemented using the GNU MPFR and MPC libraries [START_REF] Enge | MPC: A library for multiprecision complex arithmetic with exact rounding[END_REF][START_REF] Fousse | MPFR: A multipleprecision binary floating-point library with correct rounding[END_REF]. The results that follow are from tests run on an Intel Core i7 CPU running 3.2GHz with 16GB of memory. All code is written in C and compiled using Apple clang version 12.0.5 and is available at https://github.com/trcameron/CompEA. 5.1. Compensated Horner's Method. In this section, we illustrate the running error bound of the compensated Horner method in Theorem 3.2 and the computation time of the compensated Horner method in comparison to the Horner method implemented in double and quadruple precision.

5.1.1. Running Error Bound. The numerical experiment is designed as follows: We evaluate the expanded form of p(z) = (z -(1 + i)) 5 for 2000 z values near the root 1 + i. For each value of z, we compute CompHorner (p, z), the running error bound in Theorem 3.2, and the a priori error bound in (3.4). Then, we compute an accurate evaluation of p(z) in a high precision format, which allows us to measure the forward error in the computation CompHorner (p, z). The results are shown in Figure 1, where Im (z) = 1 and Re (z) ranges from 0.99 to 1.01 with a step size of 10 -5 . Note that as z gets closer to the root 1 + i, the condition number increases and the a priori error bound becomes more pessimistic. The running error bound is more accurate as it takes into account the rounding errors that occur during the computation. In addition to comparing the computation time of the compensated Horner method to the Horner method implemented in double and quadruple precision, we compare the computation time of the compensated Horner Method when the fast accurate summation method from [START_REF] Rump | Accurate floating-point summation part I: faithful rounding[END_REF] is used versus the doubly compensated summation method in Algorithm 3.3. Note that this test illustrates our point made in Section 3.1 that the doubly compensated summation method is more efficient in this case.

The experiment is designed as follows: We form 100 random complex polynomials of each degree from 100 to 2500, incremented by 10. For each polynomial, the computation time required to evaluate the polynomial at a random complex number is recorded. Then, the average time for each degree is displayed in Figure 2. Note that the compensated Horner method is about 2 times faster with the doubly compensated summation method versus with the fast accurate summation method. Moreover, the compensated Horner method with the doubly compensated summation method is about 4 times faster than the Horner method implemented in quadruple precision. In this section, we illustrate the limiting accuracy of the compensated Ehrlich-Aberth method in Corollary 4.3. In addition, we compare the accuracy and computation time of the compensated Ehrlich-Aberth method to the Ehrlich-Aberth method implemented in double and quadruple precision. Finally, the experiments in this section illustrate the robustness of the stopping criteria for the compensated Ehrlich-Aberth method discussed in Section 4.1.

Indeed, in all experiments, the stopping criteria was used to suspend the compensated Ehrlich-Aberth iterations and, as evident by the accuracy of the root approximations, these iterations were not suspended prematurely. 5.2.1. Limiting Accuracy. The experiment is designed as follows: Chebyshev polynomials of the first kind are created for each degree from 5 to 80, incremented by 1. It is known that the roots of the Chebyshev polynomial become more ill-conditioned as the degree of the polynomial increases. In Figure 3, we compare the maximum condition number of the roots of each Chebyshev polynomial to the maximum forward error in the computed roots using the Ehrlich-Aberth method implemented in double precision and the compensated Ehrlich-Aberth method. In addition, we include the limiting accuracy bounds for both methods.

Error.

In this section, we illustrate that the backward errors of the root approximations computed by the compensated Ehrlich-Aberth method are similar to the backward errors of the root approximations computed by the quadruple precision implementation of the Ehrlich-Aberth method. To that end, we investigate polynomials with well-conditioned roots and polynomials with ill-conditioned roots.

Throughout this section, we let Ehrlich-Aberth and Ehrlich-Aberth-Quad denote the Ehrlich-Aberth method implemented in double and quadruple precision, respectively. In addition, we let Ehrlich-Aberth-Comp denote the compensated Ehrlich-Aberth method. Finally, we make use of the following measurement of backward error: Given a set of root approximations z 1 , . . . , z m for the polynomial p, we define the relative backward error of all root approximations by (5.1)

η(z 1 , . . . , z m ) = p -q ∞ p ∞ ,
where • ∞ denotes the infinity vector norm, and q is the polynomial with leading coefficient equal to the leading coefficient of p and whose roots are exactly z 1 , . . . , z m . Note that we compute the coefficients of the polynomial q in a high precision format using the GNU MPC and MPFR libraries [START_REF] Enge | MPC: A library for multiprecision complex arithmetic with exact rounding[END_REF][START_REF] Fousse | MPFR: A multipleprecision binary floating-point library with correct rounding[END_REF]. First, we consider polynomials with well-conditioned roots. For degree m = 10, 20, 40, . . . , 1280, we create 10 polynomials with random complex coefficients that have real and imaginary parts distributed in the interval [-1, 1]. We compute the roots of each polynomial using Ehrlich-Aberth, Ehrlich-Aberth-Quad, and Ehrlich-Aberth-Comp, and the relative backward error and the elapsed computation time is recorded. Then, the average over all 10 polynomials for each degree is recorded on the left of Figure 4. On the right of Figure 4, the same test is repeated for polynomials of the form p(z) = 1 + 2z + Note that slope of the linear regression line is reported for both the backward error and elapsed time in Figure 4. In particular, the dependence of the relative backward error on the degree of the polynomial has been reduced from nearly quadratic for Ehrlich-Aberth to linear for Ehrlich-Aberth-Comp. In addition, the computation time suggests that all three methods have a quadratic cost complexity. However, Ehrlich-Aberth-Quad is approximately 100 -200 times slower than Ehrlich-Aberth, whereas, Ehrlich-Aberth-Comp is only 2 -4 times slower.

Next, we compare the backward error and computation time when working with polynomials with ill-conditioned roots. For degree m = 10, 20, 30, . . . , 150, we create 10 polynomials with random complex roots on the unit circle. Then, we compute their roots using Ehrlich-Aberth, Ehrlich-Aberth-Quad, and Ehrlich-Aberth-Comp, and the relative backward error and the elapsed computation time is recorded. The average over all 10 polynomials for each degree is recorded on the left of Figure 5. On the right of Figure 5, the same test is repeated for the truncated exponential:

p(z) = 1 + x + 1 2 x 2 + 1 3! x 3 + • • • + 1 m! x m ,
maximum relative forward error of the root approximations is recorded on the left of 6. Next, for degree m = 5, 6, . . . , 25, we consider monic mth degree polynomials whose roots have small imaginary parts: j + (-1) j 8µi, for j = 1, 2, . . . , m. For each polynomial, the maximum relative forward error of the root approximations is recorded in the middle of Figure 6. Finally, we consider Wilkinson polynomials of degree m = 5, 6, . . . , 20. For each polynomial, the maximum relative forward error of the root approximations is recorded on the right of Figure 6. Note that the results in Figure 6 illustrate the limiting accuracy of the compensated Ehrlich-Aberth method indicated by 4.3. That is, the root approximations in examples have a relative forward error as if they were computed in twice the working precision and then rounded into the working precision. Moreover, in each of the tests, we can see the point at which the maximum condition number of the polynomial roots became too large and we lost the additional pseudo-precision afforded to us by the compensated Horner's method.

Next, we consider polynomials with multiple or near multiple roots. In particular, the Mandelbrot polynomial, which has roots that lie in a fractal (Mandelbrot set), and is defined recursively by p 0 (z) = 1 and p j (z) = zp 2 j-1 + 1, for j = 1, 2, . . . , k, with degree m = 2 k -1. For degree our purposes, we consider m = 63, i.e., k = 6. Also, we consider the Kameny polynomial, which is defined by

p(z) = (c 2 z 2 -3) 2 + c 2 z 9 .
For c = 10, 10 3 , these polynomials have two close real roots with 3, 10 common decimal digits, respectively, and a complex pair with very small imaginary parts: 10 -4 , 10 -13 , respectively. Finally, we consider the following polynomials with multiple roots: Note that the non-multiple roots of p 1 (z) and p 2 (z) are well-conditioned, whereas the non-multiple roots of p 3 (z) are ill-conditioned.

In Table 1, we display the maximum relative forward error in the root computations via Ehrlich-Aberth, Ehrlich-Aberth-Quad, and Ehrlich-Aberth-Quad. In addition, the worst case limiting accuracy is displayed for each polynomial to illustrate the result in Corollary 4.3. Finally, note that this test illustrates the importance of using the compensated Horner's method to evaluate the polynomial derivative. Indeed, if the Horner method implemented in double precision is used, then the forward error of Ehrlich-Aberth-Comp is 2.42 • 10 -1 for the Mandelbrot polynomial. Hence, in this case, the limiting accuracy is not attained unless the compensated Horner's method is used to evaluate the polynomial and its derivative.

Conclusion.

The compensated Ehrlich-Aberth method is effective for the accurate computation of all roots of a polynomial. In Theorem 3.2, we proved a running error bound for the compensated Horner method. Then, in Section 4.1, we used this error bound to form robust stopping criteria for the compensated Ehrlich-Aberth method that guarantees iterations do not terminate until the relative error in the polynomial evaluation is too large or until the relative size of the Ehrlich-Aberth correction term is too small. Moreover, in Corollary 4.3, we showed that under suitable conditions, all root approximations have a limiting accuracy as if computed in twice the working precision and then rounded into the working precision. Finally, in Section 5, extensive numerical experiments illustrate the accuracy of the compensated Horner and Ehrlich-Aberth methods as well as the speed-up in terms of computation time as compared to the quadruple precision implementations of the Horner and Ehrlich-Aberth methods, respectively.

  As a result, for a, b ∈ F, floating-point operations satisfy fl (a • b) = (a • b)(1 + ), where | | ≤ µ. This further implies that (2.1) |fl (a • b) -a • b| ≤ µ |a • b| and |a • bfl (a • b)| ≤ µ |fl (a • b)| .

Algorithm 2 . 1

 21 Error-free transformation of (a, b) ∈ F 2 for addition. function [x, y] = TwoSum (a, b) : x = fl (a + b) z = fl (xa) y = fl ((a -(xz)) + (bz))

Algorithm 2 . 2

 22 Error-free transformation of (a, b) ∈ F 2 for multiplication. function [x, y] = TwoProduct (a, b) : x = fl (a • b) y = FMA (a, b, -x)

Theorem 2 . 1 .

 21 Let a, b ∈ F. Then, for [x, y] = TwoSum (a, b) we have a + b = x + y, x = fl (a + b) , |y| ≤ µ |x| , |y| ≤ µ |a + b| , and for [x, y] = TwoProduct (a, b) we have

Algorithm 2 . 3

 23 Error-free transformation of (a, b) ∈ C 2 for addition. function [x, y] = TwoSumCmplx (a, b) : [Re (x) , Re (y)] = TwoSum (Re (a) , Re (b)) [Im (x) , Im (y)] = TwoSum (Im (a) , Im (b)) straightforward extension of Algorithm 2.1 and is shown in Algorithm 2.3. Note that Algorithm 2.3 requires 12 flops to be executed. The error-free transformation of (a, b) ∈ C 2 for complex multiplication requires multiple products of the real and imaginary parts of a and b as shown in Algorithm 2.4. Note that Algorithm 2.4 requires 20 flops to be executed. In contrast, if the splitting operation from [6] was used in Algorithm 2.2, then Algorithm 2.4 would require 64 flops to be executed, see Algorithm 3.3 from [17]. Algorithm 2.4 Error-free transformation of (a, b) ∈ C 2 for multiplication. function [w, x, y, z] = TwoProductCmplx (a, b) : [g 1 , h 1 ] = TwoProduct (Re (a) , Re (b)); [g 2 , h 2 ] = TwoProduct (Im (a) , Im (b)) [g 3 , h 3 ] = TwoProduct (Re (a) , Im (b)); [g 4 , h 4 ] = TwoProduct (Re (a) , Im (b))

Theorem 2 . 2 .

 22 Let a, b ∈ C. Then, for [x, y] = TwoSumCmplx (a, b) we have a + b = x + y, x = fl (a + b) , |y| ≤ µ |x| , |y| ≤ µ |a + b| , and for [w, x, y, z] = TwoProductCmplx (a, b) we have a • b = w + x + y + z, w = fl (a • b) , |x + y + z| ≤ √ 2γ 2 |a • b| .

Algorithm 4 . 1

 41 Compensated Horner's Method for p and p . function [c 0 , c 0 , s 0 ] = CompHornerDer (p, z) : h m = a m , h m = 0, h m = 0 e m = 0, e m = 0, e m = 0

Fig. 1 .

 1 Fig. 1. Running Error Bound Test

Fig. 2 .

 2 Fig. 2. Computation Time Test

Fig. 3 .

 3 Fig. 3. Limiting Accuracy Test

Fig. 4 .

 4 Fig. 4. Well-Conditioned Polynomial Roots Test

Fig. 6 .

 6 Fig. 6. Simple Ill-Conditioned Roots Test

p 1 (

 1 z) = (z + 1) 5 (z 50 + z + 1),p 2 (z) = (z -1) 4 (z 2 + z + 5) 3 (3z -1) 2 (z 50 + 1), p 3 (z) = (z -15) 2 15 i=1 (zi).

  respectively, as the coefficient of degree k in p π , p µ , p ν , p σ .

	end for
	It is immediately clear that h 0 = Horner (p, z). Furthermore, by induction, it is
	easy to show that
	p(z) = h 0 + (p π + p µ + p ν + p σ )(z),
	and it follows from [17, Proposition 5.3] that
	( (p π + p µ + p ν ) + p σ )(|z|) ≤ γ 2m p(|z|).

Table 1

 1 PolynomialLimiting Accuracy Ehrlich-Aberth Ehrlich-Aberth-Comp Ehrlich-Aberth-Quad Mandelbrot 2.53 • 10 -5 4.42 • 10 -1 3.04 • 10 -8 1.70 • 10 -11 Kameny, c = 10 1.11 • 10 -16 1.42 • 10 -14 1.77 • 10 -16 1.11 • 10 -16 Kameny, c = 10 3 1.11 • 10 -16 5.11 • 10 -11 1.25 • 10 -16 1.11 • 10 -16 p Multiple and Near Multiple Roots Test

	1 (z) p 2 (z) p 3 (z)	N/A N/A N/A	3.76 • 10 -3 7.19 • 10 -4 1.56 • 10 -2	3.02 • 10 -6 8.40 • 10 -8 7.86 • 10 -8	8.88 • 10 -7 1.90 • 10 -8 1.04 • 10 -8

This work was partly supported by the NuSCAP (ANR-20-CE48-0014) project of the French National Agency for Research (ANR).

for degree m = 10, 20, 30, . . . , 100. Note that the growing condition number of the polynomial roots causes the accuracy of Ehrlich-Aberth to rapidly decline. In contrast, Ehrlich-Aberth-Comp is able reduce the relative backward error of all computed roots to O(µ) for polynomials of degree less than 80. Of course, once the condition number of the polynomial roots is too large, say greater than 10 16 , the additional accuracy afforded by the compensated Horner method is lost.

Forward Error.

In this section, we illustrate the improvements made to the forward error of root approximations using the compensated Ehrlich-Aberth method developed in Section 4.1. In particular, we investigate polynomials with simple ill-conditioned roots and polynomials with multiple and near multiple roots. Throughout this section, we let Ehrlich-Aberth and Ehrlich-Aberth-Quad denote the Ehrlich-Aberth method implemented in double and quadruple precision, respectively. In addition, we let Ehrlich-Aberth-Comp denote the compensated Ehrlich-Aberth method.

First, we consider polynomials with simple ill-conditioned roots. In particular, for degree m = 5, 6, . . . , 20, we consider the monic mth degree polynomial with prescribed roots 2 -m/2 +j -3, for j = 0, 1, . . . , m -1, where • denotes the floor function. For each polynomial, the