Jens Taprogge
email: jens.taprogge@csem.ch

Felix Beyeler
email: felix.beyeler@femtotools.com

Alexander Steinecker
email: alexander.steinecker@csem.ch

Bradley J Nelson
email: bnelson@ethz.ch

Parallel packaging of Micro Electro Mechanical Systems (MEMS) using Self-Alignment

Keywords: flip-chip, C-4, self-alignment, MEMS, MOEMS, no clean, force sensor, solder paste, bismuth, FemtoTools, wafer level packaging, WLP, chip scale packaging, CSP, direct chip attach, DCA, DRIE, SOI

Packaging is one of the major cost drivers for MEMS devices. Currently wire bonding is the dominant method for electrically connecting MEMS chips to substrate. Using self-alignment a method for packaging multiple MEMS at the same time has been developed. The presented process achieves high throughput and precise alignment at low cost. The Controlled Collapse Reflow Chip Joining (C-4) process has been adapted to the specific requirements of MEMS. The combination of coarse robotics and liquid solder self-alignment guarantees precise positioning and alignment of the individual MEMS chips to the respective substrates. The new method has been implemented in a case study. In the study force sensors have been packaged. Precise angular alignment of the sensors is critical for receiving accurate measurements. Results of the case study are presented.

Introduction

A large part of MEMS manufacturing costs (60% to 80%) originates from packaging [1, p. 6]. This work deals with reducing these costs. In particular a novel bonding method is presented.

The bonding methods used in MEMS are largely adopted from the integrated circuit (IC) industry. This comes as no surprise since these methods have been highly optimized, are well understood by chip foundries and, furthermore, have been certified for various uses. In fact, MEMS devices are often produced on depreciated CMOS manufacturing equipment. While Controlled Collapse Reflow Chip Joining (C-4) or flip-chip [START_REF] Miller | Controlled collapse reflow chip joining[END_REF] has seen widespread adoption in IC packaging, this is not the case in MEMS packaging. Here wire bonding is still largely preferred. This is due to several factors: Wire-boding is a flexible technology that is easy to adapt to new designs. It is easier to reduce shear stress caused by thermal expansion mismatch (TEM) when using wire bonding. Finally, flip-chip cost savings increase with pin count. [START_REF] Lau | Cost analysis: solder bumped flip chip versus wire bonding[END_REF] However, most MEMS devices typically have a low pin count. By using wirebonding, MEMS devices cannot take advantage of C-4 characteristics such as: ─ Free placement of pads resulting in excellent space utilization (up to chip scale packages (CSP)). ─ Shorter connections leading to better electrical performance (low inductance and resistance of connections) and excellent thermal performance. ─ No separate mechanical fixation.

This work presents a way to use C-4 in MEMS packaging. The C-4 process is adapted to the specific properties of MEMS devices. By doing so it becomes possible to package multiple MEMS devices simultaneously with a single machine. The number of devices that can be packaged in parallel is only limited by the ratio of substrate space to device space and the size of the device wafer. The novel method also reduces the precision requirements of the pick and place robotic manipulator. The above improvements are achieved by relying on liquid solder surface forces causing chip and substrate to self-align.

Attempts to assemble individual micro devices in parallel generally rely on one of two approaches: ─ a deterministic one requiring multiple manipulators and sensory feedback (thereby being very similar to the multiplication of serial processes), ─ or a stochastic one relying on self-assembly which offers low equipment costs but severely limits the size and complexity of the parts that can be assembled [START_REF] Yang | Automated Microassembly[END_REF][START_REF] Böhringer | Microassembly[END_REF].

The approach presented in this paper combines the advantages of both: It has the deterministic properties of the robotic placement of the chips while not requiring multiple manipulators. It uses classical robotics for coarse positioning and self-alignment for the precise positioning that is required for many MEMS and MOEMS.

Parallel C-4 Bonding

In the C-4 process substrate and chip have a matching array of solder bumped contact pads. During the assembly the chip is flipped upside-down and placed at the desired position on the substrate. The solder is liquefied and bridges the facing pads. When the solder cools down, lasting mechanical and electrical connections are formed. Due to the flipping operation this process is often referred to as a flip-chip process. [START_REF] Miller | Controlled collapse reflow chip joining[END_REF] Currently a parallel C-4 process is not feasible since each chip needs to be precisely aligned to a substrate. The combined use of deterministic and stochastic assembly methods presented in this work lowers that requirement.

Process Overview

Our process also relies on solder bump joining. One of its essential features is that the chips are not flipped, however. Instead they are connected to the substrate when they are already separated from the wafer (by means of dicing or etching) but still located at their original position. A robotic manipulator places several substrates above the wafer. The substrates are spaced out at an interval matching the spacing of chips in the wafer or a multiple thereof. As the robotic manipulator lowers the set of substrates, the corresponding pads of substrate and chip are connected through the liquefied solder. The substrate is then lifted pulling the chips along. Once the chips are free hanging the surface forces of the molten solder will work towards minimizing the distance between matching pads and align the chips (see fig. 1).

Maximum alignment precision can be achieved by adding differently sized (not necessarily electrically connected) pads on substrate and chip as suggested by [START_REF] Bache | Bond design and alignment in flip chip solder bonding[END_REF][START_REF] Patra | Modeling of Self-Alignment Mechanism in Flip-Chip Soldering-Part II[END_REF].

Required robotic precision

Assuming that the solder bumps form a spherical shape when molten, the maximum placement precision required is half the size of the solder pad diameter. Once the solder wets, surface forces finely align the pads and, thereby, the MEMS chips.

Space saving

MEMS chips typically require very few electrical connections. The individual pads can be larger and more spread out than the ones found on ICs. Since they can be placed freely on the chip, little chip real estate is lost. Because the pads on the sub-Fig. 1. On the bottom the wafer table is depicted. Above that identical substrate holding devices that are attached to a robotic actuator are located. The arrows are pointing onto the measuring tips of three MEMS chips that have just been separated from the wafer. The separation of the MEMS chips from the wafer is assisted by suction needles that are reaching through the substrate.

strate are placed underneath the chip rather than around the chip, the total package size can be reduced compared to wire bonding in most cases.

Case Study

The invented method was verified in a case study. A demonstrator was built to bond FemtoTools FT-S260 force sensors (see fig. 2) to substrates. The chip is directly attached (DCA) to an organic substrate. The substrate is an FR-4 printed circuit board (PCB). Apart from the MEMS chip the PCB additionally contains readout electronics.

For the final device to produce accurate measurements it is essential that the sensor is exactly orientated with regard to the substrate. The sensor is manufactured from a Silicon on Insulator (SOI) wafer using reactive ion etching (RIE) and deep reactive ion etching (DRIE) processes. Because the force sensor must directly contact its environment, it does not receive protective packaging. The transducer fabrication process does not require dicing of the wafer. Instead the device wafer is placed on a support wafer. During the etching process the individual chips are completely separated from the device wafer. The chips can move freely in a 200 µm gap created by the etching process, and their precise position is no longer defined. The chip as well as the substrate feature square pads with a side length of 380 µm. The maximum diagonal distance between pads is 3.6 mm.

The force sensor chip converts the force to be measured into a deflection using elastic flexures. The deflection is detected by two capacitor electrode arrays (differential comb drives). Liquid coming into contact with the comb drive electrodes will cause the electrodes to permanently stick together. Flux cleaning solutions can, therefore, not be used. [START_REF] Sun | Characterizing fruit fly flight behavior using a microforce sensor with a new comb-drive configuration[END_REF] Prior to this work the devices have been manually wire bonded. The pad layout of the chip and substrate were slightly modified for the new process to maximize selfalignment torque and to balance chip support.

Before being assembled the individual devices are in-wafer tested for defects using the resonance method described in [START_REF] Yang | Automated optical inspection method for MEMS fabrication[END_REF]. Only known good dies (KGD) are assembled.

Fig. 2. FemtoTools FT-S260 force sensor

In the demonstrator implementing the described method, batches of only three devices 1 are bonded at a time. This keeps machine and tooling costs low while still proving the viability of the approach.

Under Bump Metallurgy (UBM)

The pads and traces on the substrate are made of copper. The pads are plated with 3-6 µm of nickel and 50-150 nm of gold. The contact pads of the MEMS chip are also gold plated (200 nm) to prevent oxidation and guarantee good solder wettability. The adhesion layer consists of 10 nm of titanium.

Solder

The solder used in this flip chip application has to fulfill several requirements:

1. It has to be soft enough to limit the stress exposed on the contact pads and thereby the MEMS chip. 2. The melting point has to be lower than that of the solder used to assemble the other parts on the substrate. The other components are assembled to the PCB first, since the sensor is the most fragile part. A low melting point also reduces the stress that results from the mismatch of coefficients of thermal expansion (CTE) during the cool-down after the assembly process. 3. Since the precise alignment of chip and substrate is achieved by self-alignment, the solder has to provide large surface forced in the liquid state.

It has to wet the contact pads well.

A solder paste containing the eutectic Bismuth (46 wt.%), Tin (34 wt.%), Lead (20 wt.%) alloy (Indium Technology Indalloy 42) was selected for this work. It features a melting point of 96°C. According to Indium Technology, the maximum operating temperature of the solder is 90% of the absolute melting point temperature [START_REF] Scalzo | Private correspondence with Indium Corporation employee Mario Scalzo on August 18[END_REF]. In the case of the alloy used this is 59°C. Flux contained in the solder paste guarantees good wetting of the pads and removes the need for cleaning of the pads or a controlled assembly environment.

Gold embrittlement. The large amount of tin contained in the solder leads to high surface forces in the liquid state. However, it can create an intermetallic bond with the protective gold layers of the contact pads. This in turn tends to accelerate solder joint fatigue. To prevent this from occurring, the thickness of the gold layer has to be limited to approximately 380 nm. [START_REF] Glazer | Effect of Gold on the Reliability of Fine Pitch Surface Mount Solder Joints[END_REF]

Process Steps

Solder dispensing. Due to the relatively large pads used it is not necessary to use expensive solder bumping methods such as vacuum metallization to obtain tinning of the chips' contact pads. Instead, solder paste is dispensed onto the substrate prior to the actual assembly taking place. For batches of significant size mask printing of solder paste is feasible.

Wafer Registration. The wafer is placed onto a holding device. The position of the devices is not necessarily well defined with regard to the wafer. Two alignment marks are located by a custom image processing system [START_REF] Wyss | Tileye-Inspection of Implant Surfaces[END_REF]. Relative to the marks the position of each chip contained in the wafer is known.

Pre-Heating. The wafer holding mechanism contains a closed loop heating system that pre-heats the wafer to a temperature slightly below the melting point of the solder. This reduces the time required to heat the individual chips.

Heating and Chip Pick-up. A robotic manipulator moves an array of substrates to a position above the wafer with the corresponding substrate and chip pads coarsely aligned. The solder on the substrates' pads is liquefied by heating and then brought into contact with the pads of the chips. Once the chips' pads are wetted, the substrate is lifted again pulling the chips along. Depending on the weight of the chips, their adhesion to the wafer, the geometry and number of the solder pads and the surface tension of the solder, the lifting forces exerted by the solder may not be sufficient. This has been observed in the case study. To produce the forces required to achieve separation of the chip from the wafer, retractable vacuum pipes reaching through holes in the substrate have been introduced. These vacuum pipes are retracted as soon as the chip is separated from the wafer and the solder is molten.

Self-Alignment. Substrate, the free hanging chip and solder are heated in a temperature controlled way. Once the solder is molten and the chips are free-hanging, selfalignment takes place. When this is completed the device is cooled down and lasting electrical and mechanical connections are formed.

Results

The demonstrator is operational and has produced several batches of devices (see fig. 3). The devices have been functionally tested and perform as desired. The selected low-temperature solder as well as the relatively large height of the solder connections guarantees a level of shear stress that does not impact the performance of the sensor.

One of the critical points of this assembly method was the flux reaching comb drive, leading to the destruction of the device (see section 3). This has not been observed for any of the manufactured devices.

Shear tests on the devices have shown the connection to withstand a force of up to 10 N. These measurements, the good electrical connections and cross sections (see fig. 4) demonstrate good wetting of the contact pads. The alignment of the chip to the substrate is within the specification of FemtoTools.

Cycle time of the assembly is 18 s per batch of three devices. While this compares favorably with manual wire bonding (15 s per device for wire bonding alone, e.g. excluding mechanical fixation), it is not exceptionally fast. However by increasing the batch size (by means of additional substrate holders), the throughput can be increased substantially. For the devices used in this case study the theoretical minimum assembly time per device is below one second. This could be drastically improved by switching to a larger wafer or smaller substrate. However, this has not yet been experimentally verified.

Conclusions and Outlook

This work presents a method providing high-throughput and precise-placement packaging of MEMS at low cost. The method requires only limited changes to existing MEMS designs. By building on C-4 technology many of the potential advantages become available to the MEMS designer: chip scale packaging (CSP), stacking of MEMS and electronics chips, improved electrical and thermal performance are just some of them. The presented method also features properties that may prove beneficial to new designs. Contrary to regular flip chip approaches the MEMS chip is handled from one side only. This gives the chip designer the option of placing very fragile structures on most of the underside of the chip. Owing to the parallel processing of the assembly, good throughput can be achieved at relatively low handling acceleration. This can be desirable for very sensitive devices.

Due to tooling requirements, the cost of a machine implementing the bonding method will highly depend on the amount of parallelization chosen. The manufacturer can select the appropriate amount of parallelization to match the desired throughput. This significantly reduces the cost for producing small to medium size batches as found in specialized devices.

Fig. 3 .

 3 Fig. 3. Assembled device. The white spot on the chip is an artifact of chip production.

Fig. 4 .

 4 Fig. 4. Micrograph of a cross section through a solder connection The PCB is located at the top, the MEMS chip on the bottom.

When there are not enough devices left in a row, or chips have found to be faulty during in-wafer testing, fewer devices are bonded in one batch.

Acknowledgements. This work was supported by the European Commission via the project HYDROMEL (FP6 NMP2-CT-2006-026622).


Serena Ruggeri
email: serena.ruggeri@ing.unibs.it

Gianmauro Fontana
email: gianmauro.fontana@itia.cnr.it

Claudia Pagano
email: claudia.pagano@itia.cnr.it

Irene Fassi
email: irene.fassi@itia.cnr.it

Giovanni Legnani
email: giovanni.legnani@ing.unibs.it

Handling and Manipulation of Microcomponents: Work-Cell Design and Preliminary Experiments

Keywords:

The paper introduces an experimental setup for the automatic manipulation of microcomponents, based on a 4 dof robot with Shoenflies motion and a two-camera vision system. The general architecture of the workcell is presented. The work-cell functionality was tested via repeatability experiments using a set of vacuum grippers. Due to their intrinsic simplicity, vacuum grippers are very cheap and appear a promising solution for micromanipulation. An innovative nozzle for a vacuum gripper was designed, fabricated and tested, comparing its performance with traditional needles. The design was conceived to reduce the frequency of occlusions of the gripper and handle a wide range of particles. The performed tests evaluate the success and precision of the part release. Indeed, this is a crucial aspect of micromanipulation because microparts tend to stick to the gripper preventing the successful performance of manipulation tasks. The results confirm that adhesive effects prevent the release of components. For this reason different strategies were adopted in order to improve the efficiency in the release of stuck components. This solution increases the percentage of release and, setting appropriately the intensity of the pressure, it does not affect negatively the accuracy nor the repeatability of the positioning.

Introduction

During the past decades several microproducts have been fabricated for a great variety of applications in the traditional fields and in more innovative areas. The majority of microsystems are currently fabricated by semiconductor-based manufacturing techniques, taking advantage of the consolidated experience gained from the fabrication of integrated circuits. More recently, a variety of emergent products (e.g. micro-resonators, radio frequency devices, drug delivery systems, chemical and biochemical sensors) requires components made of non semiconductor material in order to increase the performance and functionalities of the whole system. These components have to be made of metal, plastic and ceramic, which need appropriate fabrication processes and have to be assembled with each other or with semiconductor structures.

These hybrid three dimensional microproducts have still great difficulty in penetrating the market, mainly due to the limits of the fabrication processes that require manipulation and final assembly of microcomponents. However, new market perspectives can be reached automating the assembly phase. The main challenge is due to the new physical scenario when dealing with millimetric and sub-millimetric parts. Indeed, due to the high surface to volume ratio of microcomponents, superficial forces are predominant and the manipulation of microparts significantly differs from that of macroscopic ones. For this reason, the downscaling of standard assembling procedures at the microlevel is infeasible or inefficient and many operations have to be done manually [START_REF] Van Brussel | Assembly of Microsystems[END_REF].

Many issues can affect one or more of the three main manipulation phases (grasp, handling, release) and have to be taken into account during the design of the microgripping systems. Due to the charging effects, the parts move into their energetically most favourable configuration, leading in an uncontrolled grasping or release. Moreover, due to the adhesion forces, objects stick to the gripper preventing the release that is not anymore facilitated by the gravitational force due to the reduced weight.

In the literature, several solutions to these manipulation problems can be found, including hybrid-type grippers where two principles are integrated in order to reduce the adhesive force during the release [START_REF] Chen | A hybrid-type electrostatically driven microgripper with an integrated vacuum tool[END_REF], grippers based on physical principles peculiar of the microscale that control and exploit the adhesive forces [3][4], and handling systems without physical contacts between the gripper and the component, in order to avoid stiction [START_REF] Gosse | Magnetic tweezers: Micromanipulation and force measurement at the molecular level[END_REF]. Among the strategies proposed for the micromanipulation there are: phase transition [START_REF] Kochan | European project develops ice gripper for micro-sized components[END_REF], magnetic [START_REF] Gosse | Magnetic tweezers: Micromanipulation and force measurement at the molecular level[END_REF], van der Waals [START_REF] Feddema | Micro-assembly planning with van der Waals force[END_REF], electrostatic [START_REF] Hesselbach | Centering electrostatic microgripper and magazines for microassembly tasks[END_REF][8], adhesive and capillary interactions [START_REF] Biganzoli | Development of a gripping system based on capillary in[END_REF][9]- [START_REF] Saito | Capillary force with a concave probe-tip for micromanipulation[END_REF] suction [START_REF] Zesch | Vacuum tool for handling microobjects with a NanoRobot[END_REF], and laser [START_REF] Rambin | Micro-assembly with a focused laser beam[END_REF]. Each of these solutions has its own advantages and drawbacks, mainly concerning the cost, accuracy, repeatability, compliancy, versatility and complexity.

One of the strategies downsized from the macroworld is the use of the force generated by the pressure difference between the gripper and the atmosphere. Indeed, vacuum grippers are very common in the assembly of fragile macrocomponents and can be easily miniaturised. They can be very cheap as consist mainly of a micropipette connected to a vacuum pump. However, both the gripper and the part should have smooth surfaces to prevent air leakage. Moreover, like all the contact grippers at the microscale, the adhesion forces significantly affect the release so that the precise positioning is difficult to achieve. This is the reason why some expedients for the release of microcomponents have been proposed. For instance, a short pressure pulse can be applied to assist the release of the microcomponent, even if it can affect the accuracy of the positioning [START_REF] Zesch | Vacuum tool for handling microobjects with a NanoRobot[END_REF]. The adhesion due to the electrostatic force can be reduced coating the glass pipette for the suction with a conductive layer connected to the ground [START_REF] Petrovic | Gripping tools for handling and assembly of microcomponents[END_REF].

In this paper, an experimental setup devoted to the automatic manipulation of microcomponents is introduced and preliminary experiments on the grasping and releasing of parts are described and discussed. The performance of two standard vacuum microgrippers (commercially available needles for dispensing) with respect to an innovative multi-lumen nozzle (Fig. 1a) is critically analysed. The limitations of these grippers, mainly in the releasing phase, are highlighted and some solutions proposed.

Experimental Setup

A suitable experimental setup (Fig. 1) able to move the parts and measure their position in the working area was designed.

The work-cell is equipped with a Mitsubishi Electric RP-1AH robot [START_REF] Van Brussel | Assembly of Microsystems[END_REF]. It presents a 5-joint closed link structure and 4 degrees of freedom with Schoenflies motion: 2 revolute joints for the positioning in the x-y working area, a third revolute joint for rotation and a prismatic joint for the z vertical end-effector motion. The operating limits are 150x105 mm 2 with a vertical stroke of 30 mm. The repeatability is ±5 µm in the x-y plane, ±10 µm for the vertical motion and ±0.02° for the end-effector rotation.

A smart and standard mechanical interface (2) was realized in order to facilitate the tool change. It was directly connected to the bottom part of the hollow screw constituting the third and fourth axis of the robot.

The vacuum generation system is a critical part of the setup, mainly during the releasing phase. Two vacuum generation systems were tested. The former consisted of an air compressor, a FRL (Filter Regulator Lubricator) group, a normally closed solenoid valve mounted on the top of the robot arm and a piINLINE Micro Ti vacuum ejector installed directly on the hose near the suction point. The latter used the same air compressor and FRL group, but was equipped with a piCOMPACT10 vacuum ejector (3). This ejector integrates a vacuum sensor and two normally closed solenoid valves, one for the supply and one for the release. This latter and more complex generation system was chosen to assist the release with a positive pressure. By modifying the throttling, the entity of the blow in the release phase can be set and optimized for the specific component and application.

The measurements of the position of the parts in the focal plane were performed using a suitable vision system, consisting of a camera (4) with field of view (FoV) 16.3x13.5 mm 2 and spatial resolution 6.6 µm and a second camera [START_REF] Gosse | Magnetic tweezers: Micromanipulation and force measurement at the molecular level[END_REF] with FoV 25x18.9 mm 2 and spatial resolution 18.4 µm. The parts lied on a transparent glass substrate (6) so that the camera, fixed on a rigid structure below the robot working area, detects from the bottom their position and orientation. An opportune lighting system is essential for the detection, robust recognition and reliable measurement, thus a diffuse illumination of the scene was adopted, making the disturbance of the environment light negligible. In order to obtain better images, the end-effector was also equipped with a contrast panel [START_REF] Feddema | Micro-assembly planning with van der Waals force[END_REF] on the top of the gripping tool [START_REF] Enikov | Optically transparent gripper for microassembly[END_REF].

Concerning the gripping tools, the nozzle (Fig. 1a) has a central circular hole (radius=55.4µm) and an external circle (radius=140.1µm), circumscribing the internal hole and the other four holes. This geometry was conceived with the purpose of manipulating a wider range of components, avoiding the need of changing the tool for the manipulation of components different in size or geometry. The nozzle should be able to perform like a needle with dimension equal to its internal radius and a needle with dimension equal to its external one at the same time. Keeping that in mind, the performance of the nozzle was compared with the performance of two needles with internal diameters of 100 and 260 µm (Fig. 1b). The camera and the robot were calibrated in order to allow the automatic part detection for the grasp and after the release and for the measure of the part position and orientation. A non-conventional calibration process was conceived and implemented in order to simultaneously calibrate the camera and georeference the camera with respect to the robot, without using external tools. This procedure was based on the use of a virtual grid with the same characteristics of the standard calibration pattern, created by moving a sphere gripped by the robot end-effector in the field of view of the camera. A similar procedure was presented in [START_REF] Tamadazte | Multiscale Calibration of a Photon Videomicroscope for Visual Servo Control: Application to MEMS Micromanipulation and Microassembly Sensors & Transducers Journal[END_REF]. Thus, a suitable vision algorithm was developed in LabVIEW for a complete calibration of the manipulation cell, compensating for perspective, distortion and spatial referencing errors. Moreover, a set of vision functions and image processing algorithms were implemented, including a robust pattern matching algorithm for the detection of the parts in the field of view and the elaboration of the displacement measurements.

Experimental Procedures

Repetitive tests of the grasping and releasing of microcomponents were carried out with the setup previously described. The initial tests showed that the results depended on the boundary conditions, for example the complexity of the trajectory of the robot during the handling of the components. Indeed, the release of the object was achieved more easily when the gripper performed a simple up and down movement, than if the pick and place motion includes a horizontal displacement, as if the gripper motion favoured the adhesion object-gripper. The tests were performed using a pick and place motion. For these preliminary experiments, the performance of each gripper was assessed using spheres with different diameters: 0.8 mm and 1.2 mm and for each sphere, a total of 30 grasp-motion-release cycles were performed. Moreover, three strategies for the release were adopted a) the vacuum pump was switched off; b) the robot moved brusquely 2 mm upwards; c) the controller activated the release valve applying a positive pressure. The feeding pressure was 2 bar, the vacuum flow 0.04 Nl/s and the execution time of the tests approximately 12 seconds. The environmental conditions were kept constant during all the tests: the temperature was 24 °C and the relative humidity 25%. Given the size of the manipulated objects, a clean room environment was not necessary.

The releasing tests were performed approaching the sphere to the planar surface and releasing it from a distance of about 0.15 mm. The precision (accuracy, repeatability) was estimated measuring the position of the sphere just before and after the release. To avoid undesired rolling of the sphere their shape and the horizontality of the plane were verified.

Results and Discussion

For each test cycle the accuracy (AC) and the repeatability (RP) were calculated. The accuracy AC is defined as the distance between the target position and the barycentre of the cloud of the actual achieved positions, while the repeatability RP is defined as the statistical positional deviation from the average, according to the standard ISO 9283. Although this standard is defined for the performance of manipulators, it was assumed significant also in this context.

The smallest needle (internal diameter 100 µm) was initially tested with both the spheres, but its dimension was so small that it grasped the component only if in contact with it. Since the contact could damage the component and the gripper, it was not considered for the present tests. A setup with an appropriate feedback should be developed in the future to allow a deeper investigation. In the following, the term needle will refer to the needle with internal diameter equal to 260 µm. In Table 1 the obtained results are shown. In the experiments exploiting the strategy a described in section 3 the failure of the release was much more probable than the success. Despite of the low releasing percentage, for the sake of completeness, also the accuracy and repeatability of such tests are shown in the table. Regarding the comparison of the performance of the nozzle and the needle using the inertial force (par.3 strategy b) for the release, the results showed that the nozzle achieves lower percentage of release of the biggest sphere. This is probably due to the superficial force contribution of the area of the nozzle between the central hole and the external ones. The opposite is true for the smallest sphere; in this case the nozzle achieves a greater number of releases. This might be due to a smaller surface contact of the nozzle. Indeed, the surface contact between the needle and the sphere is always the same, regardless the dimension of the sphere. Instead, in the case of the nozzle, it is likely that the contact surface is smaller for the smallest sphere than for the biggest one, where all the external circle is probably in contact with the sphere. This could also be the reason why the nozzle showed a higher probability of releasing the smallest sphere compared to the biggest one.

Concerning the precision of the positioning, no great difference can be found for each of the two spheres, regardless the type of gripper. In the case of the needle the accuracy value is high, as well as the repeatability value. Therefore, even when the percentage of the release is very high the accuracy of the positioning is unacceptable for any precise application.

A more positive scenario appears with the nozzle, even if the results are still too inadequate for a microprecision manipulation. It is likely that the procedure can be improved in order to reduce the energy of the sphere after the release. This could reduce the displacement of the sphere after the detachment from the gripper. Nevertheless, all in all, as expected and suggested by the literature, this release strategy seems to be not suitable for precise manipulation.

In order to achieve better results the strategy c was tested. It consists in assisting the release applying a small positive pressure. This strategy increases the percentage of release and, surprisingly it does not affect negatively the accuracy and repeatability of the positioning. However, the parameters that control the blow had to be optimized for each case; indeed the success of the experiments seemed to be significantly influenced by the chosen parameters, including the entity of the depression for the grasp and of the positive pressure for the release. A threshold of values, which allows the achievement of a high percentage of release and a great accuracy of the positioning was in some cases not found, mainly because of the rough regulation of the vacuum pump. In these cases, the experiments were done, trying to improve the precision rather than the success of the release. This strategy was preferred considering the future implementation of a control system able to indicate whether the release occurred or not. In any application such a system should be easily implementable, and preferable than dealing with a warranted but not precise release.

This issue occurred in the test of the nozzle and the small sphere, for which the achieved release percentage was lower than 60%. In fact the release achieved with the nozzle seemed to be more unpredictable than the release using the needle. This might be due to the complex shape of the nozzle, so that uncontrolled air flows through its five holes could occur.

Conclusions

An experimental setup based on a robot, a vacuum generator system and a vision system was designed and developed. It allowed the automatic manipulation of microcomponents with a variety of vacuum grippers. The experiments were carried out comparing the performance of a standard vacuum gripper (a standard dispending needle) and a multi-lumen nozzle in grasping and releasing microspheres.

Two strategies were tested in order to achieve a correct and precise release, one based on the inertial force of the sphere and the other on the positive pressure. The percentage of release is significantly improved using the second strategy, even though the procedure for the optimization of the parameters is much longer.

In the case of the needle also the accuracy of the release was improved using the positive pressure release strategy, whereas the behaviour of the nozzle was unvaried, probably due to uncontrolled air flows through its five holes. The experiments showed that the performance of a traditional needle was higher than that of the nozzle. However, the nozzle seemed to be more flexible in handling components of various dimensions, whilst the needle size has to be chosen depending on the component dimensions.

Deeper research will be carried out on the nozzle, in order to understand the behaviour showed during the release with positive pressure, developing also some

Fig. 1 .

 1 Fig. 1. Experimental setup.

Table 1 .

 1 Accuracy, repeatability and release percentage of the two spheres by the two grippers.

						Sphere Dimension [mm]
	Release	Gripper type	AC[mm] RP[mm] Release[%] AC[mm] RP[mm] Release[%] 1.2 0.8
	No strategies	Needle Nozzle	0.283 0.181	0.913 0.380	13.33 16.67	0.106 0.382	0.438 1.524	16.67 16.67
	force	Needle	0.175	5.094	96.67	0.775	4.992	43.33
	Inertial	Nozzle	0.402	2.791	40.00	0.383	3.053	70.00
	Positive	pressure	Needle Nozzle	0.026 0.030	0.263 0.330	90.00 96.67	0.027 0.098	0.199 0.589	100.00 56.67

theoretical and numerical models. This will also help to identify if, for suitable dimensions of the spheres, the external holes can work as a mean of auto-centring of the sphere after the release.


domain: info math

text:
S Graillat
email: stef.graillat@lip6.fr

F Jézéquel

M S Ibrahim

Dynamical control of Newton's method for multiple roots of polynomials *

Keywords: Discrete Stochastic Arithmetic, floating-point arithmetic, numerical validation, rounding errors, polynomial, multiple roots, Newton's method, modified Newton's method AMS subject classifications: 65H10, 65G20, 65H05, 65-04

In this article, we show how to perform a dynamical control of Newton's method for the computation of multiple roots of polynomials. Using Discrete Stochastic Arithmetic, root approximations are computed until the difference between two successive approximations is numerical noise. With such a stopping criterion, the optimal number of iterations in Newton's method are performed. Moreover it is possible to estimate in the result obtained which digits are in common with the exact root. Two strategies to estimate the multiplicity of polynomials roots are compared: one requires root approximations computed at different precisions and the other three successive iterates of Newton's method. We show that using such a strategy and then the modified Newton's method, multiple roots can be computed with a requested accuracy.

Introduction

The approximation of polynomial roots is usually based on a sequence computation and it might be difficult to perform the optimal number of iterations. Indeed if the computation stops too early, the result may be still far from the exact root and too many iterations may cause an increase of the rounding error that affects the result. In this article we show how to perform a dynamical control of Newton's method for the approximation of a root α of multiplicity m of a polynomial f , i.e. such that f (i) (α) = 0 for i = 0, . . . , m -1 and f (m) (α) = 0. Thanks to Discrete Stochastic Arithmetic [START_REF] Vignes | Discrete Stochastic Arithmetic for validating results of numerical software[END_REF], a stopping criterion that takes into account rounding errors can be used and the optimal number of iterates can be computed. Moreover, from the root multiplicity m, it is possible to estimate the number of correct digits in the approximation obtained.

In this article, two strategies are compared to estimate the multiplicity of polynomial roots. It is wellknown that deciding whether a univariate polynomial has a multiple root is an ill-posed problem: an arbitrary small perturbation of a polynomial coefficient may change the answer from yes to no. For example a real double root may be changed into two simple roots. As a consequence, it is hardly possible to verify that 2 Estimation of rounding errors using Discrete Stochastic Arithmetic

Principles of the CESTAC method

Based on a probabilistic approach, the CESTAC method [START_REF] Vignes | A stochastic arithmetic for reliable scientific computation[END_REF] allows the estimation of rounding error propagation which occurs with floating-point arithmetic. When no overflow occurs, the exact result, r, of any non exact floating-point arithmetic operation is bounded by two consecutive floating-point values R -and R + . The basic idea of the method is to perform each arithmetic operation N times, randomly rounding each time, with a probability of 0.5, to R -or R + . The computer's deterministic arithmetic, therefore, is replaced by a stochastic arithmetic where each arithmetic operation is performed N times before the next one is executed, thereby propagating the rounding error differently each time. It has been proved [START_REF] Chesneaux | L'arithmétique stochastique et le logiciel CADNA[END_REF] that the computed result R of n elementary arithmetic operations is modelled to the first order in 2 -p as:

R ≈ Z = r + n i=1 g i (d)2 -p z i , ( 1
)
where r is the exact result, g i (d) are coefficients depending exclusively on the data and on the code, p is the number of bits in the mantissa and z i are independent uniformly distributed random variables on [-1, 1].

From Equation 1, we deduce that:

1. the mean value of the random variable Z is the exact result r;

2. the distribution of Z is a quasi-Gaussian distribution.

Then by identifying R and Z, i.e. by neglecting all the second order terms, Student's test can be used to estimate the accuracy of R. From N samples R i , i = 1, 2, ..., N ,

∀β ∈ [0, 1], ∃τ β ∈ R s.t. P R -r ≤ στ β √ N = 1 -β, (2)
where

R = 1 N N i=1 R i and σ 2 = 1 N -1 N i=1 R i -R 2 . ( 3
)
τ β is the value of Student's distribution for N -1 degrees of freedom and a probability level 1 -β. Therefore, if no overflow occurs, the number of decimal significant digits common to R and r can be estimated as

C R = log 10 √ N R στ β . ( 4
)
Thus the implementation of the CESTAC method in a code providing a result R consists in:

• performing N times this code with the random rounding mode, which is obtained by using randomly the upward or downward rounding mode; we then obtain N samples R i of R;

• choosing as the computed result the mean value R of R i , i = 1, ..., N ;

• estimating with Equation 4 the number of exact decimal significant digits of R.

In practice β = 0.05 and N = 3. Indeed, it has been shown [START_REF] Chesneaux | L'arithmétique stochastique et le logiciel CADNA[END_REF][START_REF] Chesneaux | Sur la robustesse de la méthode CESTAC[END_REF] that N = 3 is in some reasonable sense the optimal value. The estimation with N = 3 is more reliable than with N = 2 and increasing the size of the sample does not improve the quality of the estimation. The probability of overestimating the number of exact significant digits of at least one is 0.054% and the probability of underestimating the number of exact significant digits of at least one is 29%. By choosing β = 0.05, we prefer to guarantee a minimal number of exact significant digits with a high probability (99.946%), even if we are often pessimistic by one digit. The complete theory can be found in [START_REF] Chesneaux | L'arithmétique stochastique et le logiciel CADNA[END_REF][START_REF] Vignes | A stochastic arithmetic for reliable scientific computation[END_REF].

Validity of the CESTAC method

Equations 1 and 4 hold if the two following hypotheses are verified:

1. the rounding errors α i are independent, centered uniformly distributed random variables; 2. the approximation to the first order in 2 -p is legitimate.

Concerning the first hypothesis, with the use of the random arithmetic, rounding errors α i are random variables, however, in practice, they are not rigorously centered and in this case Student's test gives a biased estimation of the computed result. It has been proved [START_REF] Chesneaux | L'arithmétique stochastique et le logiciel CADNA[END_REF] that, with a bias of a few σ, the error on the estimation of the number of exact significant digits of R is less than one decimal digit. Therefore even if the first hypothesis is not rigorously satisfied, the estimation obtained with Equation 4 is still correct up to one digit.

Concerning the second hypothesis, the approximation to the first order only concerns multiplications and divisions. Indeed the rounding error generated by an addition or a subtraction does not contain any term of higher order. It has been shown [START_REF] Chesneaux | L'arithmétique stochastique et le logiciel CADNA[END_REF] that, if both operands in a multiplication or the divisor in a division become insignificant, i.e. with no more exact significant digit, then the first order approximation may be not legitimate. In practice, the CESTAC method requires, during the execution of the code, a dynamical control of multiplications and divisions, which is a so-called self-validation of the method.

Principles of DSA

The self-validation of the CESTAC method requires its synchronous implementation. Indeed to enable the estimation of the accuracy, the samples which represent a result must be computed simultaneously. Discrete Stochastic Arithmetic (DSA) [START_REF] Vignes | Discrete Stochastic Arithmetic for validating results of numerical software[END_REF] has been defined from the synchronous implementation of the CESTAC method. With DSA, a real number becomes an N -dimensional set and any operation on these N -dimensional sets is performed element per element using the random rounding mode. The number of exact significant digits of such an N -dimensional set can be estimated from Equation 4. The self-validation of the CESTAC method also leads to the concept of computational zero [START_REF] Vignes | Zéro mathématique et zéro informatique[END_REF] defined below.

Definition 2.1

During the run of a code using the CESTAC method, a result R is a computational zero, denoted by @.0, if ∀i,

R i = 0 or C R ≤ 0.
Any computed result R is a computational zero if either R = 0, R being significant, or R is insignificant. A computational zero is a value that cannot be differentiated from the mathematical zero because of its rounding error. From the concept of computational zero, an equality concept and order relations have been defined for DSA. Definition 2.2 Let X and Y be N -samples provided by the CESTAC method.

• Discrete stochastic equality denoted by ds= is defined as Xds= Y if and only if X -Y = @.0. Stochastic relational operators ensure that in a branching statement the same sequence of instructions is performed for all the samples which represent a variable. DSA enables to estimate the impact of rounding errors on any result of a scientific code and also to check that no anomaly occurred during the run, especially in branching statements.

Implementation of DSA

The CADNA1 software [START_REF] Chesneaux | L'arithmétique stochastique et le logiciel CADNA[END_REF][START_REF] Eberhart | High performance numerical validation using stochastic arithmetic[END_REF][START_REF] Jézéquel | CADNA: a library for estimating round-off error propagation[END_REF][START_REF] Jézéquel | A new version of the CADNA library for estimating round-off error propagation in Fortran programs[END_REF][START_REF] Lamotte | CADNA_C: A version of CADNA for use with C or C++ programs[END_REF] is a library which implements DSA in any code written in C, C++ or Fortran and allows one to use new numerical types: the stochastic types. In essence, classical floating-point variables (in single or in double precision) are replaced by the corresponding stochastic variables, which are composed of three perturbed floating-point values. The library contains the definition of all arithmetic operations and order relations for the stochastic types. The control of the accuracy is performed only on variables of stochastic type. When a stochastic variable is printed, only its exact significant digits appear. For a computational zero, the string "@.0" is printed. In contrast to interval arithmetic, that computes guaranteed results, the CADNA software provides, with the probability 95% the number of exact significant digits of any computed result. CADNA has been successfully used for the numerical validation of academic and industrial simulation codes in various domains such as astrophysics [START_REF] Jézéquel | Reliable computation of a multiple integral involved in the neutron star theory[END_REF], atomic physics [START_REF] Scott | Numerical 'health check' for scientific codes: the CADNA approach[END_REF], chemistry, climate science [START_REF] Brajard | Numerical Validation of Data Assimilation Codes Generated by the YAO Software[END_REF][START_REF] Jézéquel | Parallelization of Discrete Stochastic Arithmetic on multicore architectures[END_REF], fluid dynamics, geophysics [START_REF] Jézéquel | Estimation of numerical reproducibility on CPU and GPU[END_REF].

The SAM library2 [START_REF] Graillat | Stochastic Arithmetic in Multiprecision[END_REF] implements in arbitrary precision the features of DSA: the stochastic types, the concept of computational zero and the stochastic operators. The SAM library is written in C++ and is based on MPFR. The particularity of SAM (compared to CADNA) is the arbitrary precision of stochastic variables. The SAM library with 24-bit (resp. 53-bit) mantissa length is similar to CADNA in single (resp. double) precision, except the range of the exponent is only limited by the machine memory. In SAM, the number of exact significant digits of any stochastic variable is estimated with the probability 95%, whatever its precision. Like in CADNA, the arithmetic and relational operators in SAM take into account rounding error propagation.

In CADNA and SAM all operators are overloaded, therefore their use in a program requires only a few modifications: essentially changes in the declarations of variables and in input/output statements. CADNA and SAM can detect numerical instabilities which occur during the execution of the code. These instabilities are either operations involving a computational zero or cancellations, i.e. subtractions of two very close values which generates a sudden loss of accuracy. At the end of the run, each type of instability together with their occurrences are printed.

Dynamical control of Newton's method

In this section, we will first recall Newton's method and some results on its convergence. We will then state a criterion to stop the iterations and give some information about the accuracy of the computed solution. Numerical experiments will confirm our results.

Newton's method

Newton's method enables one to approximate a root α of a function f by computing the sequence

x n+1 = g(x n ) with g(x) = x -f (x) f (x) .
Newton's method is based on the fixed-point method; indeed g(x) = x ⇐⇒ f (x) = 0. If α is a simple root (f (α) = 0), then there is a local quadratic convergence. If now α is a multiple root of multiplicity m > 1 then the convergence is linear. If one knows the multiplicity, it is possible to get a quadratic convergence via the modified Newton's method defined by the iteration

x n+1 = x n -m f (xn) f (xn)
. One can find more information on Newton's method in the book [START_REF] Dedieu | Points fixes, zéros et la méthode de Newton[END_REF]. For some results on the convergence and the accuracy of Newton's method in floating-point arithmetic, see [START_REF] Graillat | Accurate simple zeros of polynomials in floating point arithmetic[END_REF][START_REF] Jiang | Accurate evaluation of the k-th derivative of a polynomial and its application[END_REF][START_REF] Tisseur | Newton's method in floating point arithmetic and iterative refinement of generalized eigenvalue problems[END_REF].

Optimal stopping criterion

The following definition makes clear the notion of decimal significant digits in common between two numbers.

Definition 3.1 The number of decimal significant digits in common between two real numbers a and b is defined in R by

• for a = b, C a,b = log 10 a + b 2(a -b) ; • for all a ∈ R, C a,a = +∞. Then |a -b| = a+b 2 10 -C a,b . For instance, if C a,b = 3
, the relative difference between a and b is of the order of 10 -3 , which means that a and b have three significant decimal digits in common.

Remark 3.1

The value of C a,b can seem surprising if we consider the decimal notations of a and b. For example, if a = 2.4599976 and b = 2.4600012, then C a,b ≈ 5.8. The difference due to the sequences of "0" or "9" is illusive. The significant decimal digits of R and r become actually different from the sixth position.

The following lemma gives a relation between the derivative of the function g at the root α and the multiplicity.

Lemma 3.1 Let us assume α is a root of multiplicity

m ≥ 2 of a polynomial f and g(x) = x -f (x) f (x) . Then g (α) = 1 -1 m .
Proof: Let us assume that α is a root of f of multiplicity m, then we can find a function h such that

∀x ∈ I, f (x) = (x -α) m h(x) with h(α) = 0, (5)
with I a small open interval around α. As a consequence, it turns out that

f (x) = (x -α) m-1 [mh(x) + (x -α)h (x)] . ( 6
)
Then

f (x) = (x -α) m-1 [(m + 1)h (x) + (x -α)h (x)] + (m -1)(x -α) m-2 [mh(x) + (x -α)h (x)] . ( 7
)
The derivative of g is

g (x) = 1 - (f (x)) 2 -f (x)f (x) (f (x)) 2 = f (x)f (x) (f (x)) 2 . ( 8
)
Therefore we have

g (x) = (x -α)h(x) [(m + 1)h (x) + (x -α)h (x)] [mh(x) + (x -α)h (x)] 2 + (m -1)h(x) mh(x) + (x -α)h (x) , ( 9
)
which implies that

g (α) = 1 - 1 m . ( 10
)
In the next two lemmas, we recall well-known results on finite expansions that can be found in any textbook on mathematical analysis. They will be used in the sequel.

Lemma 3.2 Let us consider (u n ) and (v n ) two real sequences (n ∈ N ). Then we have lim n→∞ u n v n = 1 ⇐⇒ u n ∼ ∞ v n . Lemma 3.3 Let us consider (u n ) and (v n ) two real sequences such that ∀n ∈ N , u n > 0 and v n > 0, lim n→∞ u n = L = 1 and u n ∼ ∞ v n . Then log u n ∼ ∞ log v n .
The following theorem gives a relation between the common significant digits of two successive approximations of the root and the common significant digits of an approximation and the root when the root is simple.

Theorem 3.1 Let x n and x n+1 be two successive approximations computed using Newton's method of a polynomial root α of multiplicity m = 1 (simple root).

Then

C xn,xn+1 ∼ ∞ C xn,α . Proof: If n ∈ N exists such that x n = x n+1 . Then x n+1 = α, C xn,α = ∞, and C xn,xn+1 = ∞. So, in this case, Theorem 3.1 is satisfied. Let us now consider ∀n ∈ N , x n = x n+1 . From Taylor expansion, f (x n ) = f (α) + f (λ n )(x n -α) where λ n ∈ [min(α, x n ), max(α, x n )].
Moreover, with Newton iteration, we have

x n+1 = x n -f (xn) f (xn) . As a consequence, x n+1 = x n + M n (α -x n ) with M n = f (λ n )/f (x n ) and M n → 1 when n → +∞. It turns out that lim n→+∞ (α -x n )(x n+1 + x n ) (x n+1 -x n )(α + x n ) = 1.
Applying Lemma 3.2, we obtain

x n+1 + x n x n+1 -x n ∼ ∞ α + x n α -x n . ( 11
)
Therefore, we have

x n+1 + x n 2(x n+1 -x n ) ∼ ∞ α + x n 2(α -x n ) . ( 12
)
Applying Lemma 3.3, we also obtain

log 10 x n+1 + x n 2(x n+1 -x n ) ∼ ∞ log 10 α + x n 2(α -x n ) . ( 13
)
Finally, from Definition 3.1, it turns out that

C xn,xn+1 ∼ ∞ C xn,α . ( 14
)
From Theorem 3.1, if the convergence zone is reached, then the digits common to two successive approximations x n and x n+1 are also in common with the exact root α. If iterations are performed until the difference between two successive approximations x n and x n+1 is not significant, then the transformation of x n into x n+1 is due to rounding errors. Further iterations are useless: the number of iterations has been optimized. Moreover, if the difference x n -x n+1 is not significant, the digits of x n+1 which are not affected by rounding errors are in common with x n . Therefore, from Theorem 3.1, the decimal digits of the approximation x n which are not affected by rounding errors are also in common with the exact root α.

The following theorem gives a relation between the common significant digits of two successive approximations of the root and the common significant digits of an approximation and the root when the root is multiple. Theorem 3.2 Let x n and x n+1 be two successive approximations computed using Newton's method of a polynomial root α of multiplicity m ≥ 2.

Then

C xn,xn+1 ∼ ∞ C xn+1,α + log 10 (m -1)
.

Proof: If n ∈ N exists such that x n = x n+1 . Then x n+1 = α, C xn+1,α = ∞ and C xn,xn+1 = ∞. So, in this case, Theorem 3.2 is satisfied. Let us now consider ∀n ∈ N , x n = x n+1 . Applying Taylor-Lagrange formula, it exists λ n ∈ [min(α, x n ), max(α, x n )] such that g(α) -g(x n ) = g (λ n )(α -x n ). ( 15
)
As lim

n→∞ x n = α and lim n→∞ g (λ n ) = g (α) = 1 -1 m with m ≥ 2 from Lemma 3.1, we can assume that g (λ n ) = 0 and g (λ n ) = 1.
By definition, we have

g(α) -g(x n ) = α -x n+1 . ( 16
)
By combining Equations 15 and 16 we can write

α -x n = α -x n+1 g (λ n ) . ( 17
)
Therefore, we have

α -x n+1 = α -x n+1 g (λ n ) + x n -x n+1 . ( 18
)
Because x n = x n+1 , from Equation 18, we can write

α -x n+1 x n -x n+1 = g (λ n ) g (λ n ) -1 . ( 19
)
Because lim n→∞ λ n = α, from Lemma 3.1 and Equation 19, we deduce

lim n→∞ α -x n+1 (1 -m)(x n -x n+1 ) = 1. (20)
Moreover, because lim

n→∞ x n = α, it follows that lim n→∞ (α -x n+1 )(x n + x n+1 ) (1 -m)(x n -x n+1 )(α + x n+1 ) = 1. ( 21
)
Applying Lemma 3.2, we obtain

x n + x n+1 x n -x n+1 ∼ ∞ (1 -m)(α + x n+1 ) α -x n+1 . ( 22
)
and so

x n + x n+1 2(x n -x n+1 ) ∼ ∞ (m -1) α + x n+1 2(α -x n+1 ) . ( 23
)
Applying Lemma 3.3, we obtain

log 10 x n + x n+1 2(x n -x n+1 ) ∼ ∞ log 10 α + x n+1 2(α -x n+1 ) + log 10 (m -1). (24)
Finally, from Definition 3.1, we deduce

C xn,xn+1 ∼ ∞ C xn+1,α + log 10 (m -1). ( 25
)
From Theorem 3.2, if the convergence zone is reached, then the digits common to two successive approximations x n and x n+1 are also in common with the exact root α, up to log 10 (m -1). If iterations are performed until the difference between two sucessive approximations x n and x n+1 is not significant, then the significant digits of the last approximation x n+1 which are not affected by rounding errors are in common with the exact root α, up to δ = log 10 (m -1) , as illustrated in Figure 1.

x n+1

x n α rounding errors δ Figure 1: In case α is a multiple root, representation of the last iterates x n and x n+1 , and the first digits of the exact root α: if the difference x n -x n+1 is not significant, then the digits of x n+1 which are not affected by rounding errors are in common with α, up to δ.

Numerical experiments

In the numerical experiments presented here, each polynomial root is approximated by computing a sequence (x n ) using Newton's method with DSA until the difference x n -x n+1 is not significant, i.e. x n ds= x n+1 . The algorithm used is given below. We recall that in the SAM library all arithmetic and relational operators are redefined for stochastic variables. Therefore, as usual in C implementations of Newton's method, the stopping criterion is actually written as while (x!=y).

Algorithm 1: Newton using DSA Data: x init (initial approximation) Result: x (approximation of a root of a polynomial f )

x = x init ; do y = x; x = y -f (y)/f (y); while x ds = y;
Let us first consider the following polynomial

P (x) = (19x + 5) 5 (19x + 21) 9 (19x + 46) 13 (19x + 67) 25 . ( 26
)
Its four roots with their multiplicities are presented in Table 1. From Theorem 3.2 and Table 1, in the approximations of the roots obtained, the digits which are not affected by rounding errors are in common with the exact roots, up to 1 for α 1 and α 2 , up to 2 for α 3 and α 4 . This is observed in the numerical expriment described here, despite a few exceptions due to the probabilistic aspect of DSA. 13 (19x + 67) 25 with their multiplities.

Figure 2 and Table 2 present the accuracy of the roots computed with the SAM library with respect to the precision used. For each root, 15 different precisions are considered, from 24 to 10,000 bits (equivalent to 3,010 decimal digits). In Figure 2 and Table 2 are reported: the accuracy estimated by SAM, i.e. the number of significant digits which are not affected by rounding errors, and the number of significant digits in common with the exact roots. As a remark, a log scale is used in Figure 2 for the y-axis to improve the readability of the results obtained with the lowest precisions. The accuracy is linear with respect to the decimal precision used and, as expected, it decreases as the multiplicity increases. The digits not affected by rounding errors which are provided by SAM are in common with the exact roots, up to 1 or 2. Actually the difference between the accuracy provided by SAM and the number of exact significant digits is 2 in six results out of sixty:

• one result is an approximation of α 1 , its accuracy estimated by SAM is 12 digits and its exact accuracy, 10 digits, determined from Definition 3.1 is actually C xn+1,α ≈ 10.83 digits;

• one result is an approximation of α 2 , its accuracy estimated by SAM is 6 digits and its exact accuracy, 4 digits, determined from Definition 3.1 is actually C xn+1,α ≈ 4.92 digits;

• in four results, which are approximations of α 4 , this difference is consistent with Theorem 3.2 and Table 1.

A polynomial with a root of multiplicity 100 is also considered. In this case, from Theorem 3.2, in the approximations of the root obtained, the digits which are not affected by rounding errors are in common with precision 25 , number of decimal digits not affected by rounding errors estimated by SAM and number of decimal digits in common with the exact root.

α 1 = -5/
the exact root, up to 2. Table 3 presents for P (x) = (3x -1) 100 the number of decimal digits not affected by rounding errors estimated by SAM and the number of decimal digits in common with the exact root α = 1/3. The digits not affected by rounding errors which are provided by SAM are in common with the exact root up to 2, except in one case. Indeed, if the precision is set to 750 bits (equivalent to 225 decimal digits), the accuracy of the root estimated by SAM is 4 decimal digits and its exact accuracy, reported in Table 3, is 1 decimal digit. However this exact accuracy, determined from Definition 3.1 is actually C xn+1,α ≈ 1.97 digit.

Theorems 3.1 and 3.2 can be applied to other numerical examples presented in Section 4 that describes strategies to approximate polynomial roots along with their multiplicity.

Estimation of the multiplicity of polynomial roots 4.1 Previous work

Multiple roots of polynomials have been studied for a long time. See for example [START_REF] Kyurkchiev | Initial approximations and root finding methods, volume 104 of Mathematical Research[END_REF] and especially chapter 7 in [START_REF] Traub | Iterative methods for the solution of equations[END_REF].

In [START_REF] Bini | Design, analysis, and implementation of a multiprecision polynomial rootfinder[END_REF], Bini and Fiorentino used an adaptive multiprecision version of the Ehrlich-Aberth iteration [START_REF] Aberth | Iteration methods for finding all zeros of a polynomial simultaneously[END_REF] to compute clusters of roots. For that they use a suitable restarting criterion and cluster analysis to overcome the numerical instabilities arising from multiple roots. With their method, it is possible to count the number of roots in a cluster and so to have an approximation of the multiplicity.

More recently, Zeng [START_REF] Zeng | A method computing multiple roots of inexact polynomials[END_REF][START_REF] Zeng | Algorithm 835: MultRoot-a Matlab package for computing polynomial roots and multiplicities[END_REF][START_REF] Zeng | Computing multiple roots of inexact polynomials[END_REF] showed that it is possible to accurately compute a multiple root when one knows its multiplicity. This work is based on computing the pejorative manifold introduced by Kahan [START_REF] Kahan | Conserving confluence curbs ill-condition[END_REF]. The algorithm does not require the use of multiprecision. It first computes the multiplicity structure by developing a numerical GCD-finder using a successive singular value updating and an iterative GCD refinement algorithm using a Gauss-Newton iteration. Table 3: For P (x) = (3x -1) 100 , number of decimal digits not affected by rounding errors estimated by SAM and number of decimal digits in common with the exact root.

In [START_REF] Rump | Ten methods to bound multiple roots of polynomials[END_REF], Rump proposed ten methods to bound multiple roots of polynomials. Given a univariate polynomial with a multiple root of multiplicity k, the paper discusses methods for computing a disc containing exactly k roots.

One can find reviews and presentations of the previous methods in the following books [START_REF] Mcnamee | Numerical methods for roots of polynomials[END_REF][START_REF] Mcnamee | Numerical methods for roots of polynomials[END_REF]. In the sequel, we show how to compute polynomial roots with a requested accuracy using Newton's method and taking into account rounding errors. Then we present methods to compute the multiplicity of roots and therefore benefit from the quadratic convergence of the modified Newton's method.

Computation of polynomial roots with a requested accuracy

Algorithm 2 enables one to compute using DSA an approximation of a polynomial root, taking into account a requested accuracy, i.e. a minimal number of significant digits that will not be affected by rounding errors. Newton's method, implemented as in Algorithm 1, is applied iteratively until the desired root accuracy, computed from Equation 4, is reached. The precision, i.e. the number of digits used for the computation, is initially set to the requested accuracy multiplied by a rate arbitrarily chosen by the user, for instance 1.3. This rate enables one to take into account rounding errors that are necessarily generated; it avoids to perform two iterations in stable cases when one is sufficient. In accordance with the strategy described in [START_REF] Kreinovich | Towards optimal use of multi-precision arithmetic: a remark[END_REF], the precision is doubled at each iteration. If the requested accuracy and consequently the initial precision are too low, the root accuracy is also too low and useless iterations are performed. Therefore if the accuracy of the root computed at the first iteration is very low, a message is left to the user and the computation stops.

Algorithm 2: Newton with a requested accuracy using DSA Data: Requested_accuracy, x init (initial approximation) Result: x (root approximation with a requested accuracy) P recision = Requested_accuracy * Rate ; step = 0; do step = step + 1; x = Newton (x init ); if step = 1 and C x ≤ 2 then print "increase requested accuracy" and stop ; end x init = x; P recision = P recision * 2; while C x ≤ Requested_accuracy;

Estimation of the multiplicity from a linear regression

We now study the accuracy of the root compared to the precision of the computation. We will show that the relation is linear with a slope depending on the multiplicity. Proposition 4.1 Let f be a polynomial and let α be a root of f with multiplicity m. Let x be an approximation of α computed using Newton's method. Let h = α -x and ε

= f (m) (α)h m m!
. It follows that

• if α = 0, then -log 10 |ε| = m -log 10 h α + log 10 m! f (m) (α) α m ; • if α = 0, then -log 10 |ε| = m (-log 10 |h|) + log 10 m! f (m) (α)
.

Proof: Because α is a root of f with multiplicity m, f (α) = 0, ∀i ∈ {1, ..., m -1} f (i) (α) = 0 and f (m) (α) = 0. If ε = f (m) (α)h m m!
, then

h m = ε m! f (m) (α) , (27)
and

m log 10 |h| = log 10 ε m! f (m) (α) . ( 28
)
If α = 0, we deduce

-m log 10 h α = -log 10 ε m! f (m) (α) α m , (29)
and

-log 10 |ε| = m -log 10 h α + log 10 m! f (m) (α) α m . ( 30
)
If α = 0, from Equation 28we deduce

-log 10 |ε| = m (-log 10 |h|) + log 10 m! f (m) (α) . ( 31
)
From Proposition 4.1, the multiplicity m satisfies

P = m C x + K m , (32)
where

• P = -log 10 |ε|. From Taylor expansion of f , f (x) = ε + O(h m+1 ). Therefore ε ≈ f (x),
x being an approximation of the root α. In general floating-point computations are affected by a relative error, the magnitude of which being the precision. This is similar as a perturbation of the given function f so that ε can be seen as the precision of the computation. As a consequence, P is assumed to be very close to the precision (i.e. the number of decimal digits) used for the computation.

• If α = 0, C x = -log 10 x-α α
and if α = 0, C x = -log 10 |x|. In both cases C x is the number of exact significant digits in x.

• If α = 0, K m = log 10 m! f (m) (α) α m and if α = 0, K m = log 10 m! f (m) (α) . Because of the log 10 function, K m is considered not preponderant compared to m C x .
If Algorithm 2 is used, we compute approximations of a polynomial root at several precisions. For each precision P i , we obtain an approximation x i and its accuracy C xi computed from Equation 4. From Equation 32, the different couples (C xi , P i ) are close to a line of slope m. Therefore the multiplicity m can be estimated from a linear regression. To enable the estimation of m even if one iteration in Algorithm 2 is sufficient, the line is assumed to include the origin: in this case both the precision and the accuracy are zero.

Estimation of the multiplicity from successive approximations

Here we assume that ∀n ∈ N , x n = α. If the exact solution α is known, then we can determine the multiplicity m by performing successive euclidean divisions of the polynomial f by (x -α).

To estimate the multiplicity m of a polynomial root, Yakoubsohn [START_REF] Yakoubsohn | Numerical elimination, Newton method and multiple roots[END_REF] suggested a strategy based on the following property. Proposition 4.2 Let (x n ) be the sequence of approximations computed using Newton's method of the root α of multiplicity m of a polynomial f . Then

lim i→∞ x i+2 -x i+1 x i+1 -x i = 1 - 1 m .
Proof: The Taylor expansion of f at α leads to

f (x) = f (m) (α) (x -α) m m! + ϕ(x -α)(x -α) m+1 with lim x→0 ϕ(x) = 0. (33)
Then, it follows that

f (x) = f (m) (α) (x -α) m-1 (m -1)! + [(m + 1)ϕ(x -α) + (x -α)ϕ (x -α)](x -α) m , (34)
and

f (x) = f (m) (α) (x -α) m-1 (m -1)! + ψ(x -α)(x -α) m with lim x→0 ψ(x) = 0. (35)
Therefore, from Equation 33, we have

lim x→α f (x) = f (m) (α) (x -α) m m! , (36)
and, from Equation 35,

lim x→α f (x) = f (m) (α) (x -α) m-1 (m -1)! . ( 37
)
If α = 0, then lim x→α f (x) = 0. By division, we obtain

f (x) f (x) ≈ α x -α m , ( 38
) x -α - f (x) f (x) ≈ α x -α - x -α m , ( 39
) g(x) -α ≈ α 1 - 1 m (x -α), (40)
where

g(x) = x -f (x)/f (x). Finally, it turns out that g(x) ≈ α α m + 1 - 1 m x. ( 41
)
Therefore if approximations x i are close to the root α, the points of coordinates (x i , x i+1 ) are close to a line of slope 1 -1 m and lim

i→∞ xi+2-xi+1 xi+1-xi = 1 -1 m .
From Proposition 4.2, in the convergence phase, the multiplicity m can be estimated from three successive Newton's approximations.

Numerical experiments

In the numerical experiments presented here, several polynomial roots along with their multiplicities are computed using three strategies described below.

1. Newton's method is applied iteratively using Algorithm 2 until the accuracy of the approximated root is satisfactory. To estimate the multiplicity, a linear regression is performed using the different couples (accuracy, precision) obtained, as described in 4.3.

2. Algorithm 2 is used to compute the polynomial root with the requested accuracy. The root multiplicity is estimated from three successive approximations computed using Newton's method, as described in 4.4. For performance reasons, the multiplicity is estimated only at the first iteration of Algorithm 2.

3. The convergence of Newton's method is linear in case of a root multiplicity m > 1. However if the multiplicity is known, a quadratic convergence can be obtained using the modified Newton's method. In Algorithm 3, the root is approximated with a requested accuracy using firstly Newton's method and then the modified Newton's method. The multiplicity m is estimated from three successive approximations computed using Newton's method and then m is used in the modified Newton's method. Algorithm 4 shows how the modified Newton's method is implemented using DSA. Like in Algorithm 1, computations are performed until the difference between two successive iterates is a computational zero. Because of the fast convergence of the modified Newton's method, it is sometimes observed that, for an approximation x, f (x) is a computational zero. In this case, further computations are useless. Such an approximation is returned as a solution.

by rounding errors provided by SAM are always in common with the exact root, up to δ = log 10 (n -1) if n > 1. In the case of a simple root, i.e. when the polynomial considered is 3x -1, in accordance with Theorem 3.1, the accuracy estimated by SAM is the exact accuracy, except if 10 correct decimal digits are requested. However, in this case, the accuracy estimated by SAM is 12 digits and the exact accuracy, 11 digits, determined from Definition 3.1 is actually 11.95 digits.

The first strategy and the second one differ by the estimation of the multiplicity. The estimation performed by the second strategy from three successive root approximations is always correct. It is better than the estimation performed by the first strategy that requires a linear regression. Furthermore with the second strategy, the multiplicity is correctly estimated at the first step of Algorithm 2, whereas the first strategy requires results computed at several steps.

Therefore with the third strategy, the multiplicity is estimated from three successive approximations computed by Newton's method and then the modified Newton's method can be used. One can observe that, with the third strategy, two steps of Algorithm 3 are sufficient; furthermore in the roots provided by SAM the digits not affected by rounding errors are always in common with the exact ones. Because of the quadratic convergence of the modified Newton's method, the third strategy performs better than the others: the execution times reported in Table 4 can be up to five orders of magnitude lower.

Table 5 presents results obtained with the polynomial (3x -1) n for higher values of n. Because of the execution time, only the third strategy is used. Results computed with a sufficiently high requested accuracy are reported. One can observe that two steps of Algorithm 3 are sufficient, the multiplicity is correctly estimated and in the approximations provided by SAM the digits not affected by rounding errors are always in common with the exact root.

Results obtained using the third strategy with the polynomial P defined by Equation 26 are reported in Table 7. Table 6 presents for each root of the polynomial P , its multiplicity, the initial approximation x init and the rate used to initialize the precision in Algorithm 3. Indeed a particular rate is set for each root. If the rate and consequently the initial precision are too low, the computed root may be numerical noise, i.e. a result with no correct digit because of rounding errors. If the rate is too high, the execution time may be unsatisfactory, because the computation may be performed with a precision that would be too high for the requested accuracy. In the same way we need an initial approximation of the root, we also need an approximation of the rate that can depend on the degree of the polynomial and on an approximation of the multiplicity. For each root approximation, the multiplicity is correctly determined. One can observe in Table 7 that, in the approximations provided by SAM, the digits not affected by rounding errors are always in common with the exact root, up to one. As a remark, all the approximations are obtained after two steps of Algorithm 3.

Conclusion and perspectives

In this article we show how to numerically compute polynomial roots along with their multiplicity taking into account a requested accuracy. We present a strategy to perform the optimal number of iterations with Newton's method. This strategy requires to control rounding error propagation using for instance Discrete Stochastic Arithmetic. Furthermore we show how to estimate in the root approximation obtained which digits are in common with the exact root. Two methods are compared to estimate the multiplicity of polynomial roots: one based on a linear regression from results computed with several precisions and one proposed by Yakoubsohn that uses three successive iterates of Newton's method. Yakoubsohn's method is particularly satisfactory and enables one to compute the root multiplicity and then benefit from the quadratic convergence of the modified Newton's method.

Attention should be paid to the initial precision to avoid both insignificant results and costly computations due to a too high precision. From the numerical experiments presented in this article, the initial precision should depend on the requested accuracy, the degree of the polynomial and a rough approximation of the multiplicity of the root. A perspective to this work consists in determining automatically the optimal initial precision. Another prospect would be to extend the theoretical results concerning the dynamical 4: Approximation of the root of (3x -1) n and estimation m of its multiplicity with several strategies: Newton's method and linear regression, Newton's method and estimation of the multiplicity from three successive approximations, modified Newton's method. Are reported: the number of decimal digits not affected by rounding errors estimated by SAM, the number of decimal digits in common with the exact root, the execution time, and the number of steps in Algorithms 2 and 3. 25 , number of decimal digits not affected by rounding errors estimated by SAM, number of decimal digits in common with the exact root, and execution time.

control of Newton's method to the modified Newton's method. The two stopping criteria in Algorithm 4 should be taken into account to estimate the number of correct digits in the solution obtained.

roots α 1 =

 1 -5/19 α 2 = -21/19 α 3 = -46/19 α 4 = -67/19 multiplicity m

1 :Figure 2 :

 12 Figure 2: Accuracy of the roots of P (x) = (19x + 5) 5 (19x + 21) 9 (19x + 46) 13 (19x + 67) 25 computed by SAM with respect to the decimal precision used.

Table 1 :

 1 Roots of P (x) = (19x + 5) 5 (19x + 21)9 (19x + 46)

Table 2 :

 2 19 α 2 = -21/19 α 3 = -46/19 α 4 = -67/19 bits digits SAM exact SAM exact SAM exact SAM exact

	24	7	1	0	0	0	0	0	0	0
	53	15	2	1	1	0	1	0	0	0
	100	30	5	5	2	2	1	0	1	0
	200	60	12	10	6	4	3	3	3	1
	300	90	17	16	9	8	6	5	3	2
	500	150	29	28	15	14	10	9	6	5
	750	225	44	43	24	23	16	15	9	8
	1000	301	59	59	32	31	21	21	12	11
	1500	451	90	89	49	48	33	32	18	17
	2000	602	120	119	66	65	45	44	24	23
	3000	903	180	179	99	98	68	67	36	35
	4000 1204	240	239	132	131	91	90	49	47
	5000 1505	300	300	166	165	114	113	61	59
	7500 2257	450	449	250	249	172	171	90	89
	10000 3010	601	601	333	332	230	229	121	119

For each root of P (x) = (19x + 5) 5 (19x + 21) 9 (19x + 46)

13 (19x + 67)

Table

	n requested	method	m	#digits	time (s) #steps
		accuracy			SAM exact
	1 1	10 10	lin. regr. succ. app.	1.08 1.00	12	11	5. E-3	1
	1 1	25 25	lin. regr. succ. app.	1.04 1.00	31	31	6. E-3	1
	1 1	50 50	lin. regr. succ. app.	1.01 1.00	64	64	6. E-3	1
	1 1	100 100	lin. regr. succ. app.	1.01 129 1.00	129	7. E-3	1
	1 1	500 500	lin. regr. succ. app.	1.00 649 1.00	649	8. E-3	1
	10 10	25 25	lin. regr. succ. app. 10.00 9.91	26	25	1.1E-1	4
	10	25	modif.	10.00	35	35	1.1E-2	2
	10 10	50 50	lin. regr. succ. app. 10.00 9.96	52	51	4.5E-1	4
	10	50	modif.	10.00	70	70	2.4E-2	2
	10 10	100 100	lin. regr. succ. app. 10.00 9.98 104	103	2.2E+0	4
	10	100	modif.	10.00 142	142	5.6E-2	2
	10 10	500 500	lin. regr. 10.00 520 succ. app. 10.00	519	1.5E+2	4
	10	500	modif.	10.00 714	714	1.3E+0	2
	25 25	50 50	lin. regr. 24.66 succ. app. 25.00	83	82	3.4E+1	6
	25	50	modif.	25.00	66	66	3.5E-2	2
	25 25	100 100	lin. regr. 25.32 166 succ. app. 25.00	165	2.3E+2	6
	25	100	modif.	25.00 133	133	1.0E-1	2
	25 25	500 500	lin. regr. 25.00 833 succ. app. 25.00	831	1.8E+4	6
	25	500	modif.	25.00 675	675	2.9E+0	2
	50 50	100 100	lin. regr. 49.43 167 succ. app. 50.00	165	3.8E+3	7
	50	100	modif.	50.00 131	131	1.5E-1	2
	50 50	500 500	lin. regr. 49.99 833 succ. app. 50.00	831	2.1E+5	7
	50	500	modif.	50.00 662	662	5.5E+0	2

Table 5 :

 5 For P (x) = (3x -1) n , estimation m of the root multiplicity and then computation of the root using modified Newton's method. Are reported: the number of decimal digits not affected by rounding errors estimated by SAM, the number of decimal digits in common with the exact root, the execution time, and the number of steps in Algorithm 3. roots α 1 = -5/19 α 2 = -21/19 α 3 = -46/19 α 4 = -67/19

	n	requested	m	#digits	time (s) #steps
		accuracy		SAM exact
	100	100		100.00 130	130	1.1E-1	2
	100	500		100.00 655	655	9.9E+0	2
	100	1000		100.00 1311 1311 6.6E+1	2
	200	500		200.00 653	653	1.6E+1	2
	200	1000		200.00 1305 1305 1.2E+2	2
	500	500		500.00 651	651	1.3E+1	2
	500	1000		500.00 1301 1301 2.7E+2	2
	multiplicity m	5		9		13	25
	x init		0		-1		-2	-3
	Rate		3		5		7	12

Table 6 :

 6 Roots of P (x) = (19x + 5) 5 (19x + 21) 9 (19x + 46)13 (19x + 67)25 with their multiplities, values x init (initial approximations) and Rate used in Algorithm 3.

	root requested	#digits	time (s)
		accuracy SAM exact
	α 1	50	62	61	3.0E-1
	α 1	100	123	123	1.2E+0
	α 1	200	243	243	5.1E+0
	α 1	500	603	603	5.1E+1
	α 2	50	69	68	7.9E-1
	α 2	100	124	123	3.4E+0
	α 2	200	235	235	1.9E+1
	α 2	500	569	569	1.9E+2
	α 3	50	80	79	1.7E+0
	α 3	100	134	133	8.2E+0
	α 3	200	242	241	4.5E+1
	α 3	500	564	563	4.8E+2
	α 4	50	73	72	5.2E+0
	α 4	100	122	122	3.2E+1
	α 4	200	218	218	2.0E+2
	α 4	500	507	507	3.2E+3

Table 7 :

 7 For each root of P (x) = (19x + 5) 5 (19x + 21) 9 (19x + 46)13 (19x + 67)

http://cadna.lip6.fr

http://www-pequan.lip6.fr/~jezequel/SAM

This work was partly supported by the project FastRelax ANR-14-CE25-0018-01.

1

Algorithm 3: Modified Newton with a requested accuracy using DSA Data: Requested_accuracy, x init (initial approximation) Result: x (root approximation with a requested accuracy) P recision = Requested_accuracy * Rate ; step = 0;

print "increase requested accuracy" and stop ; end else x = Modified_Newton(x init ,m); end x init = x; P recision = P recision * 2; while C x ≤ Requested_accuracy; Algorithm 4: Modified Newton using DSA Data: x init (initial approximation), m (multiplicity) Result: x (approximation of a root of a polynomial f ) x = x init ; do if f (x) ds= 0 then return x and stop ; end y = x; x = y -m * f (y)/f (y); while x ds = y;

The numerical experiments have been carried out using the SAM library on an Intel Core 2 Quad Q9550 processor at 2.83 GHz using the gcc version 4.9.2 compiler. In the numerical experiments presented here, the self-validation described in 2.2 is activated: all multiplications and divisions are controlled during the execution.

Tables 4 and5 present results obtained with the polynomial (3x -1) n for different values of the degree n that here equals the root multiplicity. In Algorithms 2 and 3, the initial approximation x init and the rate used to initialize the precision are respectively set to 0.4 and 1.3.

In Table 4, for each degree n, the root is approximated using the three strategies previously described with the following requested accuracies: 10, 25, 50, 100, and 500 decimal digits. However no result is reported if the requested accuracy and consequently the initial precision are too low. In this case, a message is printed and the computation stops. If the degree n increases, the minimal requested accuracy also increases.

Table 4 presents the multiplicity m estimated using each strategy, the accuracy estimated by SAM, i.e. the number of significant digits which are not affected by rounding errors, the number of significant digits in common with the exact root, the execution time, and the number of steps performed in Algorithms 2 and 3. In the case of a simple root, it has been observed that, one step is sufficient; furthermore the modified Newton's method is useless. Therefore if n = 1, the results obtained using the third strategy are not reported: they are the same as those provided by the second one.

The first strategy and the second one provide the same root approximation with similar execution times. In accordance with Theorem 3.2, one can observe that, in this root approximation, the digits not affected