A ( Ragt

C ( Limagrain

Vilmorin & Co

Towards seed protection using biocontrol strategies

The French "Consortium public-privé Recherche-Développement-Innovation sur le Biocontrôle" and the "Groupement d'Intérêt Scientifique Biotechnologies Vertes" (GIS BV), both public-private groups, organized a scientific workshop on "Seeds and biocontrol: solutions for tomorrow's agricultures" on October 2nd, 2019 in Paris. Seed biocontrol is an emerging research topic that could meet economic issues and environmental and societal expectations to reduce the use of chemical farm inputs and adapt agriculture practices to climate change. Seventy-eight scientists, including forty-five from the private sector joined this workshop in order to debate on the most promising research areas to address a successful implementation of biocontrol solutions. They unanimously agreed that effective biological control solutions will result from multidisciplinary approaches and require innovative in plant breeding developments. This article describes the current strengths of French research in an international context around seed and biocontrol and highlights the research priorities between the public and the private sectors, the upmost being transdisciplinarity for the benefit of future agricultures.

Introduction

The French private seed sector represents 3,3 Billion euros of Annual Turnover, almost 12 000 employees and 400 Million Euros invested in R&D in 2018/2019. The seed, as propagating material, is the first input in cropping system. Seed quality, usually defined by germination, sanitary and genetic quality, is therefore a priority for growers. Seed-and soil-borne pathogens, nematodes and insect pests can limit or inhibit germination and seedling emergence, thus impacting crop yield, and can also play a major role in epidemic of bioagressors. To secure seed quality and high yield, plant protection products are commonly used during seed production and on seeds after harvest before sowing. The development of biocontrol products in the last decades can be a key lever in seed and plantlet protection. According to the French legislation, biocontrol refers to a set of methods relying on natural mechanisms to protect plant against bioagressors, in the framework of integrated pest management. These methods are based on the action of macroorganisms, microorganisms, semiochemicals or natural substances (article L253-6 of the Code Rural et de la Pêche Maritime). Mostly defined as minor uses, seed treatments are not widely supported for authorizations. Only three biocontrol seed treatments are currently recorded in the French list of biocontrol products 1 and they are authorized only for a few crops (ANSES, 2020/03/26). 1 According to articles L253-5 and L253-7 of the Code Rural et de la Pêche Maritime (dated February 26, 2020): VOTIVO (Bacillus firmus I-1582) on sugar beet and forage beet seeds, INTEGRAL PRO (Bacillus amyloliquefaciens MBI 600) on OilSeedRape and field mustard or "Navette", against fungi other than Pythiaceae and CERALL (Pseudomonas chlororaphis MA342) on rye and wheat seeds against fungi other than Pythiaceae.

Plant-variety improvement through partial and complete resistance selection exists to cope with many diseases, however very few traits of resistance to seed or soil-borne diseases have been characterized to date. Agronomic practices (i.e. rotations, tillage...), are also used but not always sufficient to protect seeds and seedlings. This bottleneck leads actors from research and technical institutes to explore alternative solutions considering seed as a material to protect and as a vector for innovation.

A number of research areas have been identified in this study to contribute to the expansion of biocontrol use for seed protection. They are summarized here below:

Seed biology and seed's biotic environment to identify new biocontrol methods and to optimize their efficacy within various seed development stages and pedoclimatic conditions.

Role of the mother plant in enhancing seed protection (seed microbiota, seed and seedling defense gene priming) through vertical transmission of biocontrol agents to the seed and role of the treated seed in vectorizing biocontrol strategies to protect mature crops (defense elicitation, endophytes...).

Combination between genotype and biocontrol products, especially microbial inoculants known to establish molecular dialogues, to adapt varieties to biocontrol products and vice versa.

Adaptation of seed treatment formulation to biocontrol products to make them compatible with seed processes.

Development and use of new technologies, like sensors and other phenotyping tools, which could play a key role in the management of seed protection with biocontrol products and foliar applications after sowing.

Evaluation of the protection mechanism sustainability and of biocontrol application efficacy.

At the forefront of science

Seeds allow plant multiplication in most agrosystems. The development and dispersal of seeds as well as their transition to seedlings represent the most critical stages of a plant life cycle. The formation of mature seed starts after the plant pollination and goes through the development of three seed tissues: the embryo, the endosperm and the seed coat. Besides the deposition of storage reserves, seed maturation is characterized by the acquisition of functional traits that allow germination and seedling establishment until favourable conditions.

Seed vigour is a concept used in the seed industry to define the seed performances in the field, namely rate and uniformity of seed germination and seedling growth, emergence ability under unfavourable environmental conditions and retention of these characteristics during storage [START_REF] Rajjou | Seed germination and vigor[END_REF]. These characteristics are acquired during seed development and influenced by genotype x environment interactions. The production of highly vigorous seeds is a key leverage to improve and stabilize yield. It is also admitted that rapid germination limits pathogen damages. However, germination phenology is likely to be largely affected by biotic and abiotic stresses, at any stage of seed life from development to germination. The development of biocontrol solutions as key success factor to increase seed protection, needs in-depth and multidisciplinary knowledge relative to seed biology and seed interactions with the environment.

Innate seed resistance

Seeds have well-established passive physical and chemical defence mechanisms that protect them from decay-inducing organisms and phytophagous pests. For instance, they constitutively produce antimicrobial compounds or can have reinforced cell-walls. In the dry quiescent state, the seed cannot activate defence responses, so the expression of defence-related components is a programmed mechanism activated during seed maturation, while seeds are still attached to the mother-plant. It allows to install passive barriers as a preemptive strategy to protect the embryos from pathogen attack. These mechanisms are also known to be activated in dormant seeds that go through hydration cycles in the soil [START_REF] Bolingue | The MtSNF4b subunit of the sucrose non-fermenting-related kinase complex connects after-ripening and constitutive defense responses in seeds of Medicago truncatula[END_REF][START_REF] Fuerst | Defense Enzyme Responses in Dormant Wild Oat and Wheat Caryopses Challenged with a Seed Decay Pathogen[END_REF]. The activation of these passive defence mechanisms is influenced both by genetic and environmental factors on the mother plant during seed production. Indeed, there exists a large natural variation in the level of phytochemical compounds in seeds throughout species, and elicitation of defence responses in a seed can occur when the mother plant is under attack from pathogens. In addition, the environmental conditions experienced by the mother plants during seed development have a direct effect on seed dormancy and vigour at harvest. A better knowledge of these mechanisms and their activation represent an interesting lever to increase seed innate resistance during seed production or seedling establishment, or for engineering biocontrol products

Immune priming during seed multiplication

Plant defense is triggered by exposure to microorganisms and their derived molecules. In the immune response that is orchestrated by plant regulators, defense genes and antimicrobial molecules are generally induced locally and in whole-plant, commonly referring to the local acquired resistance (LAR), induced systemic resistance (ISR) and systemic acquired resistance (SAR) respectively [START_REF] Pieterse | Induced disease resistance[END_REF]. Interestingly, the sensitivity of their defense systems could also be enhanced avoiding a direct induction of defense genes [START_REF] Pieterse | Induced Systemic Resistance by Beneficial Microbes[END_REF]. This "alert stratus", also called immune priming

Seed microbiota

Dry seeds represent a stressful environment. Microorganisms associated to seeds, also called seed microbiota, have developed biological mechanisms in order to colonize them (N'Guyen et al., 2019). They are considered as the primary inoculum for plants and may have an impact on plant's possibilities to deal with biotic and abiotic stresses [START_REF] Wassermann | Understanding the Indigenous Seed Microbiota to Design Bacterial Seed Treatments[END_REF]. So, a promising lever to improve seed health and seedling establishment is to act on this primary inoculum in order to limit the transmission of seed-borne pathogens and to increase plant protection abilities against soil-borne pathogens. Seed microbiota result from distinct way of contamination, with microorganisms vertically transmitted from the mother plant and horizontally transmitted from the environment. Geographic site, soil type, soil microbiome and plant genotype impact seed microbiota composition (Klaedtke et al., 2016; Pérez-Jaramillo et al., 2019). Seed microbiota optimization could be considered using multi-strains microbial inoculum also called Synthetic communities (SynComs) applied to the mother-plant or at the surface of the seed. As the robustness of beneficial microbial trait expression is increased by redundancy (i.e by multiple strains) and dominance (i.e whenever one strain has the respective trait), communities' efficiency will depend on SynCom's design. The use of native, locally adapted plant associated microbes for plant field protection has been highlighted in several studies [START_REF] Vannier | Microbiota-mediated disease resistance in plants[END_REF]. These results raise research questions about the adequacy of the SynCom composition with the plant and the local environment, and only little is known today. Moreover, due to the influence of plant genotype and environment on seed microbiota, breeding strategies that include microbial beneficial traits and agronomic practices impact also need to be considered.

Crop genotype and biocontrol product

Each plant species is able to recruit specific rhizosphere microbial assemblages and genotypes of the same species may differ in their rhizosphere microbiome composition (Pérez-Jaramillo et al., 2016), which may affect the successful establishment of introduced biocontrol microorganisms. When colonizing roots and/or internal tissues, a biochemical dialog appears between the biocontrol agent and the crop. The plant itself, by its exudates or rootlets architecture, can favor the beneficial agent development. Hence, considering the interaction biocontrol x genotype x cultural practice will be a success key in the development of both biocontrol solutions and crop varieties. In the target environment, objectives will be to develop (i) varieties favouring the attraction and establishment of biocontrol agents as well as their elicitation effect on constitutive and inducible defences, and (ii) biocontrol solutions combining as many beneficial effects as possible: biocide for pathogens, competing for space and resource with pathogens, enhancing the response of the plant towards abiotic stresses -concomitant biostimulation effect, etc. . This will need an indepth knowledge on molecular dialogues and genetic determinisms which promote the interaction of plants with beneficial microbes, on biological mechanisms mediating constitutive and inducible defences activation in response to biocontrol and on the conditions of successful establishment and action of biocontrol agents.

Seed technology

Seed technology gathers techniques that improve seed quality. In particular, seed treatment reduces damages due to seed transmitted pathogens and damping-off diseases that can drastically damage germination and impact final crop yield. While chemicals are often small organic or inorganic molecules for which seed treatment processes and specific formulations have been developed, biocontrol products are biological agents or biomolecules that need specific formulations to guarantee their stability and efficacy [START_REF] Rocha | Seed Coating: A Tool for Delivering Beneficial Microbes to Agricultural Crops[END_REF].

Selection and production of microorganisms and bio-based molecules

Until now biocontrol products applied to seeds were previously homologated for foliar application. These products have been historically hampered by the variation in efficacy of the microbial strains employed, partly explained by the empirical selection of biocontrol agents [START_REF] Barret | Plant microbiota affects seed transmission of phytopathogenic microorganisms[END_REF]. Most of microorganisms used were isolated from soil and rhizosphere of wild and cultivated plants [START_REF] Harman | Changes in perception derived from research on Trichoderma harzianum T-22 (Myths and dogmas of biocontrol)[END_REF][START_REF] Hökeberg | Biological control of cereal seed-borne diseases by seed bacterization with greenhouse-selected bacteria[END_REF] whereas microorganisms sourced from seeds seemed underutilized.

In the current context, it appears crucial to characterize novel treatment solutions adapted to seeds, either to cope with pathogens or to positively influence seedling establishment. Special attention should be brought on the influence of (i) plant genotype, (ii) seed physiological quality during maturation and at harvest, and (iii) environmental factors on seed microbiota structure. On the one hand and interestingly, decrease of microbiota diversity (including soil-borne pathogens) is often observed during germination and seedling establishment [START_REF] Barret | Emergence shapes the structure of the seed microbiota[END_REF]. The best would be to keep only the beneficial ones. On the other hand, biobased molecules (i.e. polyphenols, polysaccharides, proteins/peptides, nucleic acids, organic acids and hormones) coming from seed exudates at imbibition and/or seed coat interfere with the dynamics of microbial communities. Interdisciplinary research and development are essential to optimize microbes (single or communities) and/or biomolecules to display a variety of physico-chemical properties adapted to seeds.

The use of micoorganisms for biocontrol can be oriented towards several kinds of products like:

• Spores (when the microorganisms can sporulate) • Microorganisms in a vegetative form • Molecules (metabolites, peptides, small RNA…) produced by microorganisms • Mix of microorganisms and bio-based molecules.

In the fermentation processes a production of a large biomass is often at the expense of the desired activity. Overcoming the trade-off between the cost of production and the expected product efficiency is a real industrial challenge that can only be solved by mobilizing front science technologies like high volume of data processing and machine learning processes. It must be underlined that production chain optimization need quantitation of the final activity that is not easy to reach with biocontrol products. For this correlation with a marker that can be traced during the fermentation is even more challenging.

Stability and conservation of biocontrol solutions applied on seeds

Several challenges appear when biocontrol solutions are applied on seeds.

First when choosing the biocontrol solution to apply on a given seed lot.

The "Classical" supply chain responds to several constraints at the scale of the seed lot (genetic, sanitary, germination) and at the scale of the whole seed process (seed processing, quality requirements, number and diversity of crops and varieties, requiring several production lines...). At the top, comes the correct association of a compatible biocontrol product and a plant genotype, or concordant supply chains of both seeds and biocontrol products in suitable quantity and quality (stability, concentration, nature of impurities, delivery time...). There are also regulatory constraints, like the country of seed origin, or the diversity of legislation according to countries (authorization of sowing of treated seeds of a given crop).

Second during application on seeds.

In most of cases, seed applications are done with water-based slurries. Extemporaneous incorporation of biocontrol agent in the seed application slurry a few minutes before seed application can limit the physiological activation of microorganisms and increase their survival at seed surface. The ideal situation would be to apply microorganisms without having to dehydrate and rehydrate them, like by: • Extemporaneous mixing of the microorganisms with the seeds at sowing. In addition, biocontrol agents and substances need to be compatible with the other compounds of the seed application (fungicides, insecticides, pigments…) and their associated coformulants (preservatives...). The slurry mixing technology can be considered only for spore seed coating; new mild coating technologies should be developed for surface coating allowing oxygenation to preserve living microorganisms. Application processes contains or are followed by a drying step, which can also lead to microorganisms death (several logs of mortality) and denature some complex substances. All these parameters imply to develop a biocontrol formulation compatible with seed treatment and different from other uses (plant or soil applications).

Third during seed logistics between application on seeds and sowing.

Challenges are adequate storage before and during shipment (including delays to perform Quality Control tests just after application or validate phytosanitary checks to cross countries boarders). Once seeds are by the final client, the conditions of storage (usually between 0 to 18 months) directly impact the biocontrol efficacy.

Evaluation of biocontrol products: from their effectiveness to their environmental impact Evaluation methods and tools allow the characterisation of products and their integration into cropping systems.

Existing evaluation methods have been developed for chemicals in conventional crop production systems based on their spectrum of action, effectiveness and mode of action. Given the diversity, selectivity and the nature of potential biocontrol products, the current evaluation methods are not appropriate. The development and deployment of biocontrol products for seed require new tool development for the assessment of their efficacy, their environmental impact or the quality of the final agricultural product. This assessment is essential to ensure the proper integration of these solutions into cropping systems.

Methods and tools for the determination of the efficacy

Screening and characterisation of the effect of biocontrol products are long and costly processes that are mainly based on trials under optimal and controlled growth conditions. However, effectiveness in the laboratory and/or greenhouse does not always translate into success in the field [START_REF] Parnell | From the Lab to the Farm: An Industrial Perspective of Plant Beneficial Microorganisms[END_REF]. The variability in the effectiveness of these products could result from variations in environmental conditions, the physiological state of the plants, or the crop practices. According to Rocha and colleagues (2019), this inconstancy field performance can be one of the main restraints for the wide application of seed coating with plant beneficial microbes. Contrary to conventional pesticides, efficacy of biocontrol products is not based on the eradication of bioaggressors but rather on the installation of a biological balance in the plant ecosystem. Thus, although they have a very positive image and a major interest in pesticide reduction, biocontrol products do not seem to offer enough guarantees as to their effectiveness, which limits their use by farmers. In order to accelerate the availability of effective biocontrol products, methods and tools for evaluating the effectiveness of these alternatives to pesticides must be developed, adapted or optimised to provide an adequate experimental framework. Currently, robustness and plasticity of biocontrol solutions checked in numerous real representative conditions are highly recommended in biocontrol seed treatment development, with associated statistical analysis that will help concluding. Environmental parameters measurements can help to identify best conditions for expression of biocontrol seed treatment efficacy. Either direct (speed of beneficial agent root colonisation, % of contaminated plantlets, severity of symptoms, grow out level in the field…) or indirect aspects (final crop yield, elicitation markers quantification…) can be quantified to check this success. But research efforts are still necessary for better understanding of the key parameters of field efficacy related to the diversity and nature of these biocontrol products and their mode of action [START_REF] Nicot | Perspectives for future research-and-development projects on biological control of plant pests and diseases[END_REF]. The study of efficacy cannot be carried out without the help of highperformance measurement, calculation and analysis tools, and their development depend on the understanding of the biological mechanisms and ecological interactions, as well as the implementation of the digital and new technologies for relevant decision-support tools, sensors and phenotyping tools. It is thanks to such development that the methodologies used, whether under controlled conditions or in the field, will be better adapted to determine in the end the field efficacy of biocontrol products.

Potential impact of biocontrol seed treatments on the performance (reliability and detection threshold) of standard methods used to detect pathogens on seeds

To a large extent, seeds traded around the world have received treatment, which can be a physical treatment, chemical or biological treatment. All these different treatments can have consequences regarding to the methods of detection of pathogens in seed lots [START_REF] Sérandat | How to adapt methods to evaluate seed health of treated seeds[END_REF]. It is important for seed trade to be able to define if pests can be detected on treated seeds and to know if the result is due to the efficiency of treatment or if the treatment affects the performance of the method for detecting the pest. With most non-chemical seed treatments, pathogens may still be detectable even though they are non-viable. So, additional measures to determine the viability of any pest detected after treatment is required. An evaluation of the effectiveness of these biocontrol treatment methods remains essential.

To reach this goal, the capacity of transmission of the pathogen from seed to plant (pathosystem) is being evaluated and the treatment effectiveness is being assessed. Viability and damage potential of the pathogens can be evaluated using the existing detection methods (grow out, culture on media…) but there is also a need to develop new complementary methods (vital staining, bio test...). Current detection methods have mainly been developed and validated for non-treated seeds. Due to the high diversity of treatments, it could be difficult to validate the methods for each type of treatment and so it is necessary either to review and adapt current methods or develop and validate new ones. In our globalized market context, the presence of pathogens on seeds needs to be detected and their damage potential needs to assessed.

Evaluation of the sustainability of biocontrol product

Evaluation of the risk of release into the environment

Biocontrol products are composed of microbial agents and/or complex mixture of molecules naturally present in the environment. They are therefore generally considered as a safe and environmentally friendly alternative to conventional pesticides. Whether one agrees or not with this assumption, all agree that further studies are needed to provide scientific evidences on the environmental fate and behaviour of biocontrol products in order to ensure environment protection, to support the sustainability of these innovative products or to encourage and accelerate the placing on the market of new biocontrol products. Indeed, placing a Plant Protection Product on the European market, biocontrol or not, requires an environmental risk assessment (Regulation EC 414/1991) by authorities involving extensive testing and a significant amount of data to be provided by the applicants. Better knowledge and understanding of potential adverse effects are crucial for the determination of appropriate data requirements as well as for the development of reliable, relevant and reproducible assessment methods.

Evaluation of the sustainability of the crop protection

According to the official French definition, biocontrol is based on managing the balance of pest populations through natural mechanisms rather than on their eradication. Under the principles of agroecology, biocontrol products must be used in combination with other leverages to ensure effective and sustainable crop protection and not be part of a dynamic of substitution of chemical plant protection products. However, care must be taken not to rule out the risks of circumvention of natural mechanisms by phytopathogenic agents that may occur with the massive and repeated use of biocontrol products. Indeed, it is important to bear in mind that natural does not mean sustainable, depending on the use made of it. It is therefore essential to improve our knowledge of potential resistance to biocontrol products and to develop tools and methods to assess and control this risk.

Conclusion

New biocontrol solutions, crop varieties and their combination are important levers to reach sustainable seed and crop protection strategies. Biology and ecology of three parties have now to be put in balance when developing biocontrol methods and crop varieties : the plant genotype, the biocontrol agent at the strain level for microorganisms and their environment at the scale of soil and rhizosphere for seed protection.

In France, research community on biocontrol started to jointly address those questions by transdisciplinary approaches however a few numbers of initiatives are specifically centered on seed protection. The scientific workshop organized in Paris on October 2019 was a first step to establish close connections between research communities focusing on seeds and biocontrol and will lead to new partnerships.

It is noteworthy that the successful integration of biocontrol strategies in technical routes will rely on the design of new cropping systems combining several levers, such as biocontrol or genetic, to reach an equilibrium between yield and sustainability. Some biostimulation strategies were already successfully implemented in technical routes. For example, commercialization of Myc (Mycorrhizal) and Nod (Nodulation) factors, both involved in the establishment of microbial (fungal and bacterial respectively) root symbiosis, allowed sustainable yield increasing by enhancing plant capacities to fix nitrogen without adding chemical fertilizers. Nod factors are lipochitooligosaccharide (LCO) able to promote germination and plant growth of both legumes and non-legumes when applied as a seed treatment. Seed-applied LCO technology for crops are already available [START_REF] Smith | [END_REF]. Developing optimal associations between biocontrol and biostimulation strategies in innovative cropping systems will lead to an agriculture less and less dependent of chemical inputs. Furthermore, in upcoming years, innovative cropping systems will more and more lean on precision agriculture. Therefore, bridges have also to be built with research communities on digital farming and cropping systems.

  (Martinez-Medina et al., 2016), involves a low fitness costs (van Hulten et al., 2006) and maintains an immunological memory that can be transmitted to the next progeny during seed formation. Noticeably, epigenetic modifications induced by priming stimulation, i.e. chromatin remodeling or DNA methylation levels, have been documented by several studies (Jaskiewicz et al., 2011) in a way that defense gene promoters become easier to activate (Conrath et al., 2015) in case of subsequent stresses (Bruce et al., 2007; van den Burg and Takken, 2009). Epigenetic modifications are one of the mechanisms that enable plants to acquire stress memory that can even be inherited by subsequent generations as illustrated by Molinier et al. (2006). A better understanding of these mechanisms transmitted to the next progeny during seed development will open opportunities to improve pest management of crops during seed production.

Annex 1. List of French academic laboratories involved in biocontrol research community (with potential application for seed protection)


Sylvie Boldo

Diane Gallois-Wong

Thibault Hilaire

A Correctly-Rounded Fixed-Point-Arithmetic Dot-Product Algorithm

Keywords: Dot Product, Sum-of-Products, Correct Rounding, Odd Rounding, Fixed-Point Arithmetic

Dot products (also called sums of products) are ubiquitous in matrix computations, for instance in signal processing. We are especially interested in digital filters, where they are the core operation. We therefore focus on fixed-point arithmetic, used in embedded systems for time and energy efficiency. Common dot product algorithms ensure faithful rounding. For the sake of accuracy and reproducibility, we want to ensure correct rounding. This article describes an algorithm that computes a correctly-rounded sum of products from inputs whose format is known in advance. This algorithm relies on odd rounding (that is easily implemented in hardware) and comes with a careful proof and some cost analysis.

I. INTRODUCTION

Many algorithms used in embedded systems rely on dot product algorithms, for example in linear algebra based systems. Notably, many applications in automotive, communication, aerospace domains use signal processing and control algorithms named filters and controllers that have dot product evaluation as core.

Because of cost, power efficiency, or performance constraints, these low-level algorithms are often implemented using fixed-point arithmetic (eg. based on the integer arithmetic unit) rather than floating-point arithmetic. Moreover, in hardware implementation (ASIC or FPGA), the multiple word-length paradigm [START_REF] Constantinides | The multiple wordlength paradigm[END_REF] (ie. the use of dedicated word-length for each computation) is key to minimize the area and/or the power consumption while preserving a good performance, and bounding the error due to the finite precision by an acceptable level.

In this article, we only focus on the dot product (also called sum-of-products) in fixed-point arithmetic. We also assume that the fixed-point formats involved (the Most Significant Bit and Least Significant Bit positions, as defined in Section II) are fixed and given by the user. We propose an efficient algorithm that ensures correct rounding and we prove its correctness. It is based on odd rounding, an unusual rounding, defined in Section II-A, which is easy to implement (the proof will be useful for the proof of our fixed-point dot product, and then of a complete filter or controller implementation).

A correctly rounded dot product is known to be difficult to achieve, as we cannot dismiss the least significant bits. A This work is supported by a grant from the "Fondation CFM pour la Recherche". solution is to rely on a Kulisch accumulator [START_REF] Bohlender | Comments on fast and exact accumulation of products[END_REF], but this may be costly when the exponent range is large. Indeed, it uses the complete bit-range for all the accumulation in order to perform only exact operations before the last rounding. Our algorithm is more efficient and still ensures correct rounding.

Other practical solutions for the dot product exist, but only provide last-bit accuracy dot product (also denoted faithful rounding: it is not a correct rounding, only a round up or round down of the exact result). The main idea is to use a slightly extended precision for the internal computation of the sum, using some extra guard bits. The number of them depends on the log2 of the number of terms to add [START_REF] Volkova | Towards hardware iir filters computing just right: Direct form i case study[END_REF].

Section II gives a brief overview of fixed-point arithmetic and odd rounding, and introduces related notations. The correctly-rounded dot product is detailed in Section III: a simple version is first presented, then an optimized one illustrated by two examples. The proof of the correct rounding is given in Section IV. Some elements of comparison of the proposed algorithm to classical ones are given in Section V. Section VI features a conclusion.

II. ODD ROUNDING AND FIXED-POINT ARITHMETIC

As for the pre-requisite, we first describe odd rounding and how to implement it, and then fixed-point arithmetic, and the used notations.

A. Odd rounding

Correct rounding is notably difficult to achieve for some functions due to the Table Maker Dilemma [START_REF] Lefèvre | Worst cases for correct rounding of the elementary functions in double precision[END_REF]. In particular, common algorithms for sums or dot products usually only ensure faithful rounding [START_REF] Lefèvre | Worst cases for correct rounding of the elementary functions in double precision[END_REF], [START_REF] Ogita | Accurate sum and dot product[END_REF]. As explained above, we want to ensure correct rounding, that requires more care to ensure.

A solution is to rely on a specific rounding mode, called rounding to odd, and using it for intermediary roundings. This odd rounding is also called von Neumann rounding because it has been first designed in the 40's by Von Neumann for the design of the arithmetic unit of the EDVAC computer [6, section 9.4] [START_REF] Burks | Preliminary discussion of the logical design of an electronic computing instrument[END_REF]. It was lately used by Goldberg when converting binary FP numbers to decimal representations [START_REF] Goldberg | What every computer scientist should know about floating point arithmetic[END_REF] and formally studied later, notably to emulate the FMA operator [START_REF] Boldo | Emulation of FMA and Correctly-Rounded Sums: Proved Algorithms Using Rounding to Odd[END_REF].

The informal definition of odd rounding is the following one. When a real number is not in the set of representable numbers F, it is rounded to the adjacent number in the format with an odd integer significand. More precisely, let be the

m + 1 - w s -2 2 0 2 -1 2 -1 2
Fig. 1: Fixed-point number in format ( , ).

rounding up and the rounding down to F. Then rounding to odd is defined by:

(x) =    x if x ∈ F, (x) if the mantissa of (x) is odd, (x) otherwise.
Note that the result of x rounded to odd can be even only when x is already a representable number.

This rounding belongs neither to the IEEE-754 required roundings for floating-point arithmetic, nor to the IEEE-1666 roundings for fixed-point arithmetic in SystemC [START_REF]IEEE Standard for Standard SystemC Language Reference Manual[END_REF]. Nevertheless, implementations are available. An effective software implementation is given in [START_REF] Boldo | Emulation of FMA and Correctly-Rounded Sums: Proved Algorithms Using Rounding to Odd[END_REF]. As explained, a hardware one is much simpler: round to zero (with the IEEE-754 flags), and perform a logical or between the inexact flag (or the sticky bit) of this first rounding step and the last bit of the mantissa to get the result (or equivalently it forces the last bit of the mantissa to 1 when the deleted bits are not all zero). From an IEEE-compliant architecture, the additional cost is therefore slight.

Numerous properties are given in [START_REF] Boldo | Emulation of FMA and Correctly-Rounded Sums: Proved Algorithms Using Rounding to Odd[END_REF]. They are mostly dedicated to floating-point arithmetic, while we are interested in fixed-point arithmetic as explained above. Fortunately, many theorems are generic enough to be used both in floatingand fixed-point arithmetics. In particular, we will use here Theorem 1, that ensures correct rounding to nearest, provided a previous odd rounding with sufficient precision: Theorem 1 (round N odd, from [START_REF] Boldo | Emulation of FMA and Correctly-Rounded Sums: Proved Algorithms Using Rounding to Odd[END_REF]). Assume an even radix β and two different formats: a working format and an extended format. Assume that the extended format has at least two more digits at all points. If ext denotes the rounding to odd in the extended format, and if denotes a rounding to nearest with any tie-breaking rule in the working format, we have

∀x ∈ R, ( ext (x)) = (x).
We refer to Corollary 2 for an explicit version of this theorem applied to fixed-point arithmetic (see Section IV).

B. Fixed-Point Arithmetic

We consider here numbers and computations based on two's complement fixed-point arithmetic. Let x be a w-bit long fixed-point number, with bits {x i } ≤i≤ :

x = -2 x + -1 i= 2 i x i (1)
where and are the positions of the most significant bit (MSB) and the least significant bit (LSB) of x, respectively. So the word-length w and the Most and Least Significant Bit positions, and are linked with

w = -+ 1. (2)
The couple ( , ) denotes the fixed-point format of x. Figure 1 exhibits a fixed-point number and its binary representation.

Internally, the fixed-point numbers are represented with an integer X (denoted mantissa) formed by the w-bit {x j } and an implicit scaling factor 2 , ie,

X = -2 w-1 x + w-2 i=0 2 i x i+ ∈ Z (3)
or equivalently x = X2 . All the operations are then integer operations on the mantissa X. We denote F the set of fixed-point integers with as LSB (or equivalently F = Z2 ), and we do not consider overflow here.

We will also need the following notations: (x), (x), (x), (x) are the rounding down, rounding up, rounding to odd and rounding to nearest (given a tie-breaking rule) of x at the bit , respectively. Moreover, we denote x + y and x + y the addition of x and y rounded at bit , using rounding down and rounding odd respectively. We assume correct rounding for addition, so we have x + y = (x + y) and x + y = (x + y).

III. ALGORITHM

Let us now describe our algorithm. A first point is that we decide to compute exactly all multiplications. The fixed-point numbers obtained each have an LSB equal to the sum of the LSBs of both corresponding operands. What is left to do is to compute a correct sum of fixed-point numbers with LSBs given in advance. Therefore, we just consider that we have n non-null inputs noted x 0 , . . . , x n-1 with respective LSBs 0 , . . . , n-1 known in advance. We are also given a desired LSB f for the output. The algorithm computes an output noted result which is a correct rounding of the sum of the inputs:

result = f (x 0 + • • • + x n-1
). This summation relies on odd rounding presented in Section II-A. Then, the intuition is that, to obtain a correct rounding to nearest at LSB f , we need an odd rounding at LSB f -2. Indeed, Theorem 1 will then guarantee a correct rounding of the sum. Of course, we do not use a Kulisch accumulator for computing this odd rounding as it would defeat our efficiency wish. Instead, we add the values sequentially, using odd rounding with a well-chosen LSB at

x 0 0 x 1 1 . . . x n-2 n-2
x n-1 n-1 lsb: inputs:

k 0 + k 1 + k n-2 + k n-1 . . . f result t 0 t 1 t n-2 t n-1
Fig. 2: General algorithm for any input LSBs.

each step to ensure that we keep the following invariant: the current computed value is the odd rounding of the sum of all the previously-seen values.

The key ingredient is the choice of the LSB at each step of the algorithm. Too large, and we cannot ensure correct rounding; too small, and the algorithm will be less efficient. For instance, if we used the smallest LSB of all inputs for all computations, we would get a correct result à la Kulisch.

A first version of the algorithm is described in Section III-A. Section III-B presents an optimized version, where inputs are sorted by LSB (the LSBs are known in advance so there is no runtime cost) to allow for overall larger computing LSBs. Examples are given in Section III-C. Note that all the proofs will be presented in Section IV; only the intuition is given here.

A. General algorithm

The general algorithm described by Algorithm 1 computes a correct rounding of the sum for any input LSBs. First of all, the LSBs k i used in intermediary computations are chosen recursively. Then, we sum the inputs x i from left to right using odd rounded additions to F ki as illustrated by Fig. 2. The t i could all be replaced with a single variable t that is updated at each step of the loop. However, using indices to identify the value at each iteration will be useful for the proof in Section IV-B.

Algorithm 1: General algorithm for any input LSBs.

Inputs:

x 0 , . . . , x n-1 , 0 , . . . , n-1 // x i has LSB i Output: result = f (x 0 + • • • + x n-1 ) 1 k n-1 ← f -2 2 for i ← n -1 downto 1 do 3 k i-1 ← min(k i , i ) -1 4 t 0 ← k0 (x 0 ) 5 for i ← 0 to n -2 do 6 t i+1 ← t i + ki+1 x i+1
are maximized when the input LSBs i are sorted in increasing order. The optimized algorithm presented in the next section builds on this remark and similar considerations to maximize the LSB at each step.

B. Optimized algorithm

The optimized algorithm is derived from the general algorithm above through a few successive optimizations. As explained previously, the first one consists in sorting the inputs by LSB, as this maximizes the computing LSBs k i . As the input LSBs are known in advance, sorting by LSB has no runtime cost. Therefore, we can assume that 0

≤ 1 ≤ • • • ≤ n-1 .
Yet to maximize efficiency, the hypothesis we need is

0 < • • • < n-1 .
Therefore, we do a pre-processing where inputs with the same LSB are grouped together as described by Algorithm 2. This may decrease the number of inputs n and modify the values of the x i and i . But it preserves the exact sum x 0 + • • • + x n-1 , as well as the fact that each x i has LSB i . And afterwards, we do have

0 < • • • < n-1 .
Algorithm 2: Pre-processing to make input LSBs distinct. 1 while there is an i such that i = i+1 do 2

x i ← x i + i x i+1 3 for j ← i + 1 to n -2 do 4 x j ← x j+1 5 j ← j+1 6 n ← n -1
The optimized algorithm for strictly increasing input LSBs is described in Algorithm 3 and illustrated by Fig. 3. Inputs are separated into two groups: on the left hand side x 0 , . . . , x m-1 whose LSBs are at most f -2, and on the right hand side x m , . . . , x n-1 whose LSBs are greater than f -2. The left hand side is handled similarly to the previous general algorithm. The main difference is that thanks to the LSBs being strictly increasing, we do not need to define computation LSBs recursively. Instead, each former k i is simply equal to i+1 -1 (except for k m-1 = f -2). On the right hand side, inputs with large LSBs are added together from right to left. At each step, we use the smaller LSB of both operands so that the operation is exact. Finally, the results for both groups of inputs are added together at f -2 (the operation is once again guaranteed to be exact) then rounded to nearest at f .

x 0 0 x 1 1 . . . . . . x m-2 m-2 x m-1 m-1 x m m x m+1 m+1 . . . . . . x n-2 n-2
x n-1 n-1 lsb: inputs:

< < < < ≤ f -2 < ≤ ≤ ≤ ≤ l 1 -1 + 2 -1 + m-1 -1 + f -2 . . . t 0 t 1 t m-3 t m-2 + n-2 + m+1 + m . . . t n-1 t n-2 t m+2 t m+1 + f -2 t m-1 t m f result
Fig. 3: Optimized algorithm.

Note that every time we use a rounded down addition + , we know that the operation is actually exact. This means we could have used any rounding mode. We choose rounding down for its efficiency: it is a simple truncation in two's complement fixed-point arithmetic.

It may be the case that all inputs are either in the set with large LSBs or the one with small LSBs, making the other set empty. Then, m is respectively 0 or n. The same algorithm holds, provided we consider all t i for i outside of the bounds {0, . . . , n -1} to be null. When m = 0, this means that we can ignore Lines 2 to 5 (both included), and we use t m-1 = 0 at Line 9. When m = n, we should ignore Lines 6 to 8, and we use t m = 0 at Line 9. Note that this convention also affects the case m = 1. Indeed, Line 5 becomes

t 0 = 0 + f -2 x 0 = f -2 (x 0 ) which overrides Line 2.
Algorithm 3: Optimized algorithm. Inputs: x 0 , . . . , x n-1 , 0 , . . . , n-1 // x i has LSB i and

0 < • • • < n-1 Output: result = f (x 0 + • • • + x n-1 ) 1 let m such that m-1 ≤ f -2 < m 2 t 0 ← 1 -1 (x 0 ) 3 for i ← 0 to m -3 do 4 t i+1 ← t i + i+2-1 x i+1 5 t m-1 ← t m-2 + f -2 x m-1 6 t n-1 ← x n-1 7 for i ← n -1 downto m + 1 do 8 t i-1 ← t i + i-1 x i-1 // addition is exact 9 result ← f (t m-1 + f -2 t m ) // add. is exact C. Example
Let us now illustrate the optimized algorithm on two examples with explicit LSBs. The first one is a toy example that showcases every part of the algorithm. The second one was generated from an actual digital filter as explained below.

As a first example, we are given inputs with respective LSBs 0, 4, 4, 5, 8, 9, 9, 12, and 20. The desired output LSB is f = 10. After pre-processing, x 0 is the input with LSB 0 = 0, x 1 is the exact sum of the two inputs with LSB 1 = 4 (computed with an addition truncated to F 4 ), and so on: Another point of view is given in Fig. 5. All bits of the inputs are shown. The red line shows the intermediary LSBs used for computations (that is to say the LSBs of the t i ) to give an idea on how many bits are kept and dismissed at each step. More precisely, each line i shows both x i and a vertical red line at the LSB of t i (remember t i ← t i-1 + i+1-1 x i ). In particular, Fig. 5 pictures the dropped bits and we clearly see that the smallest LSBs are quickly dropped.

Let us now consider an example from signal processing. We consider a 7-th order low-pass Butterworth filter with normalized cutoff frequency equal to 0.41. Its evaluation using a Direct Form algorithm1 [START_REF] Hanselmann | Implementation of digital controllers -a survey[END_REF], [START_REF] Porat | A course in digital signal processing[END_REF] requires the use of a dot product that sums 10 terms2 with respective LSBs -34, -30, -27, -25, -23, -22, -22, -22, -21, and -11. Moreover, the output should have LSB f = -7.

Once again, inputs are grouped by LSB then labeled x 0 , . . . , x 7 , where x 5 is the exact sum of the three inputs with LSB -22. But this time, all input LSBs are smaller than f -2. Technically, this means that m = n = 8 and the group of inputs with large LSBs is empty so their sum is zero. In practice, we can remove all the truncated additions from Fig. 3. We obtain Fig. 6.

As before, the bit view is given in Fig. 7 with the red line showing the intermediary LSBs used for computations. Again, the dropped bits are clear, even if we need to consider them all to ensure correct rounding.

IV. PROOF

Now let us prove the correctness of these algorithms. This Section is organized as follows: Section IV-A provides general lemmas about odd rounding. Section IV-B tackles the general algorithm and Section IV-C the optimized one.

A. Useful properties about odd rounding

Odd rounding was defined in Section II-A. We have put here all the lemmas about odd rounding that we need to prove the correctness of the algorithm.

Corollary 2. For all x ∈ R and k, ∈ Z, if k ≤ -2 then ( k (x)) = (x).
This is only an instance of Theorem 1, where the format is fixed-point arithmetic.

Lemma 3. For all x ∈ R and k, ∈ Z, if k ≤ then ( k (x)) = (x).
Proof. When k = , this is trivial:

as k (x) has LSB k = , then ( k (x)) = k (x) = (x). When x ∈ F k , we have k (x) = x, therefore ( k (x)) = (x).
Let us now assume k < and x / ∈ F k . To prove that ( k (x)) is the odd rounding of x to F , it is sufficient to prove that it is a faithful rounding of x to F and it is odd in F . As a double rounding is always a faithful rounding, the former holds. Moreover,

x / ∈ F k implies that k (x) is odd in F k so k (x) / ∈ F then ( k (x)) is odd in F . This means that ( k (x)
) is indeed an odd faithful rounding of x to F , therefore equal to (x).

Lemma 4. For all x ∈ R and y ∈ F +1 , (x) + y = (x + y).

Proof. The idea is quite simple: as y ∈ F +1 , it is even with LSB . Therefore adding it to an odd rounding makes an odd rounding of the sum. Getting a clean proof is slightly more complicated as many subcases have to be tackled.

If x ∈ F , then (x) = x. We also have that x + y ∈ F , so (x + y) = x + y. Therefore, (x) + y = x + y = (x + y).

Let us now assume that x ∈ F , meaning that (x) is odd in F . Then (x)+y is also odd in F . There is left to prove, by uniqueness of odd rounding, that (x) + y is a faithful rounding of x + y. It is sufficient to prove that

(x + y) ≤ (x) + y ≤ (x + y). When (x) = (x), we easily have (x) + y = (x) + y ≤ x + y ≤ (x + y).
The other inequality is slightly more complicated to get: we want to prove that (x + y) ≤ (x) + y, that is (x + y) -y ≤ (x). By the properties of the rounding down of x, it is enough to prove that (x + y) -y ≤ x and (x + y) -y ∈ F , and both easily hold. The case where (x) = (x) is similar.

Corollary 5. For all x, y ∈ R and k, ∈ Z, if y ∈ F and k < then

( k (x) + y) = (x + y).
Proof. As y ∈ F with > k, then y ∈ F k+1 . Therefore:

( k (x) + y) = ( k (x + y) (Lemma 4) = (x + y) (Lemma 3)
Corollary 6. For all x, y ∈ R and k, , y ∈ Z, if y ∈ F y and k < and k < y then

( k (x) + y) = (x + y).
Proof.

• If ≤ y then y ∈ F and we conclude using Corollary 5.

• If > y then:

( k (x) + y) = ( y ( k (x) + y)) (Lemma 3) = ( y (x + y)) (Corollary 5) = (x + y) (Lemma 3)

B. Proof of the general algorithm

We want to prove that the general algorithm presented in Section III-A computes a correct rounding of the sum of the inputs. Variables k i and t i and output result are defined by Algorithm 1 as illustrated by Fig. 2.

Lemma 7. For all 0 ≤ i ≤ n-1, t i = ki (x 0 +x 1 +• • •+x i ).
Proof. We proceed by induction over i.

• For i = 0 we get t 0 = k0 (x 0 ) by construction.

• For 0 ≤ i ≤ n -2 we have:

t i+1 = ki+1 (t i + x i+1 ) = ki+1 ( ki (x 0 + • • • + x i ) + x i+1 ) (I.H.) = ki+1 (x 0 + • • • + x i + x i+1 )
using Corollary 6: indeed x i+1 ∈ F i+1 and by con-

struction k i = min(k i+1 , i+1 ) -1 so k i < k i+1 and k i < i+1 .
Theorem 8. The output of the general algorithm verifies:

result = f (x 0 + • • • + x n-1 ). Proof. Using Lemma 7 for i = n -1 then Corollary 2 with k n-1 = f -2: result = f (t n-1 ) = f ( kn-1 (x 0 + • • • + x n-1 )) = f (x 0 + • • • + x n-1 )

C. Proof of the optimized algorithm

Let us now prove the correctness of the optimized algorithm presented in Section III-B. First of all, we must show that the pre-processing described by Algorithm 2 has the right properties. Its inputs are once again x 0 , . . . , x n-1 with respective LSBs 0 , . . . , n-1 . But here, we require the LSBs to have been sorted: 0 ≤ • • • ≤ n-1 . Pre-processing can modify the x j and j and even the number of inputs n. However, the fact that x j ∈ F j is obviously preserved through each step of the loop. Each step also preserves the sum x 0 +• • •+x n-1 . Indeed, all the x j are preserved (though indices get shifted down for j > i + 1) except for x i and x i+1 , which are replaced by a new x i with value x i + i x i+1 . This is where the LSBs being sorted is important: as i ≤ i+1 , both x i and x i+1 are in F i , so the addition is exact. Finally, exiting the loop means that

0 < • • • < n-1 at the end.
From now on, we just assume that we have inputs x 0 , . . . , x n-1 with respective LSBs 0 < • • • < n-1 . We want to prove that the result of Algorithm 3 (illustrated by Fig. 3) is a correct rounding to f of the sum of the inputs. Recall that m is chosen such that m-1 ≤ f -2 < m . With the convention that -1 = -∞ and n = +∞, this means that m is well defined and verifies 0 ≤ m ≤ n.

Lemma 9. For all 0 ≤ i ≤ m-2, t i = i+1-1 (x 0 +• • •+x i ). Proof.
As in the proof of Lemma 7, by induction over i:

• For i = 0 we get t 0 = 1-1 (x 0 ) by construction. • For 0 ≤ i ≤ m -3 we have:

t i+1 = i+2-1 (t i + x i+1 ) = i+2-1 ( i+1-1 (x 0 + • • • + x i ) + x i+1 ) (I.H.) = i+2-1 (x 0 + • • • + x i + x i+1 ) using Corollary 6 (indeed i+1 -1 < i+2 -1
since the LSBs are strictly increasing; and of course

i+1 -1 < i+1 with x i+1 ∈ F i+1 ). Lemma 10. If m ≥ 1 then t m-1 = f -2 (x 0 + • • • + x m-1 ). Proof. If m ≥ 2, we use Lemma 9 for i = m -2 then Corol- lary 6 with x m-1 ∈ F m-1 and m-1 -1 < m-1 ≤ f -2: t m-1 = f -2 (t m-2 + x m-1 ) = f -2 ( m-1-1 (x 0 + • • • + x m-2 ) + x m-1 ) = f -2 (x 0 + • • • + x m-2 + x m-1 ) If m = 1, recall that Line 5 of Algorithm 3 overrides Line 2 so that t 0 = f -2 (x 0 ). Lemma 11. If m ≤ n -1 then t m = x m + • • • + x n-1 ∈ F m .
Proof. By induction over i decreasing from n -1 to m, it is easy to prove that t i ∈ F i , the additions t i + i-1

x i-1 are all exact and

t i = x i + x i+1 + • • • + x n-1 .
Theorem 12. output of the algorithm verifies:

result = f (x 0 + • • • + x n-1
). Proof. Let us assume for now that 1 ≤ m ≤ n -1, so that we can use both Lemma 10 and Lemma 11 (the special cases m = 0 and m = n are handled further below). A corollary of Lemma 10 is that t m-1 ∈ F f -2 . Lemma 11 states that t m ∈ F m . By construction of m we have f -2 < m so t m ∈ F f -2 and the addition t m-1 + f -2 t m is exact. This also means that f -1 ≤ m so t m-1 ∈ F f -1 , which allows us to use Lemma 4 thereafter.

result = f (t m-1 + f -2 t m ) = f (t m-1 + t m ) = f ( f -2 (x 0 + • • • + x m-1 ) + t m ) (Lemma 10) = f ( f -2 (x 0 + • • • + x m-1 + t m )) (Lemma 4) = f ( f -2 (x 0 + • • • + x m-1 + x m + • • • + x n-1 )) (Lemma 11) = f (x 0 + • • • + x n-1 ) (Corollary 2)
Let us now handle the special cases m = 0 and m = n. As mentioned in Section III-B, any t i for i outside of the range 0, . . . , n -1 is considered to be null.

When m = 0, Lemma 11 still gives us t m ∈ F m so t m ∈ F f -2 then:

result = f (0 + f -2 t m ) = f (t m ) = f (x m + • • • + x n-1 ) (Lemma 11) = f (x 0 + • • • + x n-1 )
When m = n, we still have t m-1 ∈ F f -2 and:

result = f (t m-1 + f -2 0) = f (t m-1 ) = f ( f -2 (x 0 + • • • + x m-1 )) (Lemma 10) = f (x 0 + • • • + x m-1 ) (Corollary 2) = f (x 0 + • • • + x n-1 )

V. COMPARISON WITH OTHER ALGORITHMS

Let us now give some intuition about the cost of this algorithm, compared to the literature. Note this is a very approximated cost, guessing the number of and/or gates involved. A first approximation is that we do not consider multiplications: we only consider additions. One reason is that both Kulisch and our algorithm compute exactly all multiplications. This notably disfavors the faithful algorithm that is allowed to round the multiplication and therefore drop bits.

A second approximation is that we are not taking overflow into account: all algorithms use a MSB sufficient to ensure that no overflow will occur. If the x i have a MSB worth m i , we may use as computation MSB: f = max i i + (n -1). This could be improved if we have prior knowledge, but it will be sufficient for this rough analysis.

The chosen cost model is as follows: for the exact sum of two values x 1 and x 2 with LSBs 1 and 2 with 1 ≤ 2 , the cost is 2 ( f -1 +2). When a subsequent rounding to LSB is done, the additional cost is 0 for rounding down (assuming two's complement arithmetic); it is -1 for rounding to odd (as it only needs an or of the trailing bits, see Section II-A); and it is f -1 for rounding to nearest as there may be a carry.

The Kulisch accumulator [START_REF] Bohlender | Comments on fast and exact accumulation of products[END_REF] is a sum of all x i computed with LSB 0 , followed by a rounding to nearest. (In this section, we assume the input LSBs to be increasing, though they do not have to be stricly increasing.) The cost is then:

C K = (n -1) * 2 ( f -0 + 2) + ( f -0 ).
The faithful algorithm [START_REF] Volkova | Towards hardware iir filters computing just right: Direct form i case study[END_REF] is also a sum of all x i with the same LSB. Instead of 0 , the computation LSB is c = flog(n) -1, followed by a rounding to nearest. The cost is then:

C F = (n -1) * 2 ( f -c + 2) + ( f -c ).
Our Algorithm 3 based on odd rounding has a more complex cost. First, sorting the numbers by their LSB has cost 0 as the LSBs are known in advance. Then, let us consider the addition

t i-1 + i+1-1 x i for 1 ≤ i ≤ n -1 (with n = f -1 by convention). The LSB of t i-1 is i -1 and the LSB of x i is i . Then then cost of this addition is 2( f -( i -1) + 2) + ( i+1 -1 -( i -1)) = 2( f -i + 3) + ( i+1 -i ).
The cost is then:

C O = n-1 i=1 (2( f -i + 3) + ( i+1 -i )) + ( f -( f -2)) = n-1 i=1 2( f -i + 3) + n -0 + f -( f -2) = n-1 i=1 2( f -i + 3) + f -0 + 1.
Note that we have considered the worst case where we only have rounding to odd: the cases of equal LSBs and of large LSBs that involve correct rounding (that will be computed with rounding down) have a smaller cost.

An abstract view of the previous formulas is given by Fig. 8 and Fig. 9: they exhibit the bits involved in each of the three algorithms. The cost of the Kulisch accumulator is roughly the area of the bigger rectangle with width f -0 . The cost of the faithful algorithm is the area of the smaller rectangle with width fc . The cost of our algorithm is the area of the roughly trapezoidal region delimited on the right hand side by the LSBs i+1 -1 used in intermediary computations.

VI. CONCLUSION

We have presented an efficient algorithm for computing correctly-rounding dot products in the context of two's complement fixed-point arithmetic. Note that it easily applies to sums too: it is basically a summation algorithm that takes correct results of multiplications.

As explained, it is clearly more efficient than Kulisch accumulator (even if this effiency depends upon the input LSBs). Moreover, it is more accurate than the classical faithful rounding algorithms. Indeed, correct rounding means more accuracy (the error bound is twice as small). This could be useful in further computations. More importantly, the result is now unique, allowing reproducibility of the computation. An immediate application of this algorithm is the sum of products used inside numerical filters. Relying on this algorithm inside numerical filters iterations will improve the accuracy of the computations. It could be used for the wordlength allocation problem where the word-lengths are minimized while guarantying a final accuracy [START_REF] Hilaire | Optimal Word-Length Allocation for the Fixed-Point Implementation of Linear Filters and Controllers[END_REF].

Another perspective would be the formal proof of this algorithm. Most of the basic ingredients are already available in the Flocq library [START_REF] Boldo | Computer Arithmetic and Formal Proofs: Verifying Floating-point Algorithms with the Coq System[END_REF] (odd rounding, fixed-point algorithm, Theorem 1), but other lemmas are not. Filters have already been formalized in a previous work [START_REF] Gallois-Wong | A Coq formalization of digital filters[END_REF]. Then, the formal description of the algorithm(s) will need to be done with care, to be as readable as possible, in order to be as convincing as possible. And then this proof will allow the formal proof of more larger implementations like filters or controllers.

A more challenging perspective is the handling of overflow. We assumed here that fixed-point numbers are big enough (in the sense of a big enough MSB) in all our operations, but this has a cost. Two directions are possible. As for LSBs, we could provide correct MSBs for each operation to ensure no overflow may happen. We could also try to have smaller MSBs inside the computations by using modulo overflows if we know in advance a maximal MSB for the final result. Both directions would need careful proofs, or even formal ones, as overflow proofs are quite error-prone.

2 = 5 ,Fig. 4 :Fig. 5 :

 2545 Fig. 4: Optimized algorithm: toy example.

Fig. 6 :Fig. 7 :

 67 Fig. 6: Optimized algorithm: example from signal processing.

Fig. 8 :

 8 Fig.8: Toy example: bit view and comparison of the involved bits in the Kulisch accumulator , the faithful algorithm , and our algorithm .

Fig. 9 :

 9 Fig. 9: Signal Processing example: bit view and comparison of the involved bits in the Kulisch accumulator , the faithful algorithm , and our algorithm .

result ← f (t n-1 )The larger the LSBs used in computations, the more efficient the algorithm. Note that for i ≤ n -1, we define k i-1 as the larger LSB that verifies both k i-1 < k i and k i-1 < i . Indeed, these inequalities are required to keep on with odd rounding. For that, we will rely on Corollary 6 in the proof in Section IV. It is easy to notice that all these computation LSBs

Other possibilities for the implementation of this filter exist, implying the use of multiple-but-shorter dot products. The comparison of these different algorithms and their implementation is outside of the scope of this article.

The coefficients of this filter can be obtained with the butter function of the Scipy library (or equivalently with Matlab) with the command butter(7,0.41).

ACKNOWLEDGMENTS

We are thankful to the anonymous referees for their remarks and contructive suggestions.


Yacine El

Jean-Baptiste Pomet

Stéphane Régnier

Laetitia Giraldi

Optimal Actuation of Flagellar Magnetic Micro-Swimmers

We present an automated procedure for the design of optimal actuation for flagellar magnetic microswimmers based on numerical optimization. Using this method, a new magnetic actuation method is provided which allows these devices to swim significantly faster compared to the usual sinusoidal actuation. This leads to a novel swimming strategy which makes the swimmer perform a 3D figure-8 trajectory. This shows that a faster propulsion is obtained when the swimmer is allowed to go out-of-plane. This approach is experimentally validated on a scaled-up flexible swimmer.

I. INTRODUCTION

Untethered robotic micro-swimmers have the potential to make a serious impact on the development of new promising therapeutic techniques [START_REF] Bradley J Nelson | Microrobots for minimally invasive medicine[END_REF][START_REF] Sitti | Biomedical applications of untethered mobile milli/microrobots[END_REF][START_REF] Qiu | Magnetic helical microswimmers functionalized with lipoplexes for targeted gene delivery[END_REF][START_REF] Patra | Intelligent, self-powered, drug delivery systems[END_REF][START_REF] Michael | Minimally invasive and robotic surgery[END_REF][START_REF] Fusco | An integrated microrobotic platform for on-demand, targeted therapeutic interventions[END_REF][START_REF] Gao | Cargo-towing fuel-free magnetic nanoswimmers for targeted drug delivery[END_REF]. However, due to their microscopic size, controlling their navigation through fluids inside the body faces numerous challenges [START_REF] Lauga | The hydrodynamics of swimming microorganisms[END_REF][START_REF] Qiu | Swimming by reciprocal motion at low reynolds number[END_REF]. These devices can take different shapes and designs (helicoidal [START_REF] Tabak | Computationally-validated surrogate models for optimal geometric design of bioinspired swimming robots: Helical swimmers[END_REF][START_REF] Barbot | The rotation of microrobot simplifies 3d control inside microchannels[END_REF][START_REF] Qiu | Magnetic helical microswimmers functionalized with lipoplexes for targeted gene delivery[END_REF], beating flexible tail [START_REF] Dreyfus | Microscopic artificial swimmers[END_REF][START_REF] Islam | Magnetosperm: A microrobot that navigates using weak magnetic fields[END_REF], 2-linked structure [START_REF] Jang | Fabrication and locomotion of flexible nanoswimmers[END_REF], ciliary micro-robots [START_REF] Kim | Fabrication and manipulation of ciliary microrobots with non-reciprocal magnetic actuation[END_REF], etc.), and are actuated using multiple methods like chemical fuel propulsion [START_REF] Ma | Reversed janus micro/nanomotors with internal chemical engine[END_REF][START_REF] Jang | Multiwavelength light-responsive au/b-tio2 janus micromotors[END_REF], acoustic-based actuation [START_REF] Ahmed | Artificial swimmers propelled by acoustically activated flagella[END_REF], or external magnetic fields [START_REF] Bradley J Nelson | Microrobots for minimally invasive medicine[END_REF][START_REF] Xu | Magnetic actuation based motion control for microrobots: An overview[END_REF]. In this paper, we focus on sperm cell-inspired robotic swimmers composed of a magnetic head and an elastic tail and actuated by external magnetic fields [START_REF] Islam | Magnetosperm: A microrobot that navigates using weak magnetic fields[END_REF][START_REF] Oulmas | 3d closed-loop motion control of swimmer with flexible flagella at low reynolds numbers[END_REF]. As far as we know, the most commonly used actuation method for these swimmers is sinusoidal actuation. This consists in applying the superposition of a static orientating magnetic field parallel to the desired swimming direction and a perpendicular sinusoidal field that induces a planar symmetric beating of the tail, allowing a displacement along the swimming direction ( [START_REF] Abbott | How should microrobots swim?[END_REF]). The purpose of this paper is to propose an automatic design method for the actuation of flexible magnetic micro-robots that optimizes the swimming speed. The aim of this work is twofold: firstly, to provide a more efficient alternative to the commonly used sinusoidal actuation, and, secondly, to present an optimal control design method than can be adapted to deal with more demanding tasks, like swimming in a complex environment. Our method is based on the resolution of an optimal control problem under the constraints of an approximate dynamical model of the swimmer's displacement. In particular, we focus on computing a magnetic field shape that maximizes the horizontal swimming speed of these devices.

This actuation pattern is experimentally validated on a scaled-up flexible magnetic swimmer and compared to the classical sinusoidal actuation. The key for this approach is the use of a dynamical model that is sufficiently computationally inexpensive to be used as a constraint for an optimization process but still accurate enough to predict the displacement of a real-life low Reynolds swimmer. For this purpose, we use a 3D simplified dynamic model based on the approximation of the hydrodynamical forces using Resistive Force Theory, [START_REF] Gray | The propulsion of seaurchin spermatozoa[END_REF][START_REF] Friedrich | High-precision tracking of sperm swimming fine structure provides strong test of resistive force theory[END_REF] and the discretization of the curvature and elasticity of the tail of the swimmer, generalizing the planar "N-Link" models of [START_REF] Moreau | The asymptotic coarse-graining formulation of slender-rods, bio-filaments and flagella[END_REF][START_REF] Alouges | Self-propulsion of slender micro-swimmers by curvature control: N-link swimmers[END_REF][START_REF] Alouges | Can magnetic multilayers propel artificial microswimmers mimicking sperm cells?[END_REF]. Using this dynamic model as a constraint, the optimal control problem is numerically solved with a direct method using the software ICLOCS [START_REF] Falugi | Imperial college london optimal control software (iclocs)[END_REF]. This results in a magnetic actuation pattern that induces a significantly faster propulsion than the common sinusoidal actuation. The optimal magnetic field has two time-varying components which makes the trajectory of the swimmer non planar. In particular, the optimal swimming strategy is to perform a 3D 'figure-8' trajectory around the swimming direction. This indicates the importance of going out of plane in order to maximize the propulsion speed of flexible swimmers. The optimal actuation is experimentally validated on a scaled-up flexible magnetic swimmer using a similar setting as [START_REF] Oulmas | 3d closed-loop motion control of swimmer with flexible flagella at low reynolds numbers[END_REF] and compared to the sinusoidal field. The horizontal displacement of the experimental swimmer is accurately predicted by the dynamic model for both fields. These results show the usefulness of the simplified swimmer model for the design of optimal controls that are usable in experiments.

II. EXPERIMENTAL SETTING

The flexible swimmer used in experiments consists of magnetic disk (Neodyminum-Iron-Boron permanent magnet) with 0.3mm in height and 0.77mm in diameter attached to a silicone tail with 7mm in length and 1mm in diameter (see Figure 1). The swimmer is immersed in pure glycerol to ensure low-Reynolds conditions (≈ 10 -2 ). It is placed at the center of three orthogonal Helmholtz coil pairs that generate the actuating magnetic field. Each pair of coils is driven by a servoamplifier (Maxon Motor) which outputs a constant current for a fixed input voltage. Two cameras provide a side view and a top view of the swimmer, and its position in 3D is tracked in real-time using Visp software [START_REF] Marchand | Visp for visual servoing: a generic software platform with a wide class of robot control skills[END_REF]. More information on the experimental setting used can be found in [START_REF] Oulmas | 3d closed-loop motion control of swimmer with flexible flagella at low reynolds numbers[END_REF][START_REF] Xu | Propulsion characteristics and visual servo control of scaled-up helical microswimmers[END_REF]. We consistently use a fixed frame (x, y, z) where the z axis is vertical and the y axis is along the axis of the side camera. In what follows , the desired swimming direction will be along the x-axis.

The propulsion of the experimental swimmer is characterized by measuring its velocity-frequency response (given by the dotted line in Figure (3)) under a magnetic field of the form :

B(t) = B x B y sin(2πf t) 0 T , (1)
where B x = 2.5mT ,B y = 10mT .

FIG. 1. Experimental setup ([21]

). The three orthogonal Helmholtz coils generate a homogeneous magnetic field in the center, where the swimmer has been placed. The swimmer is tracked using two perpendicular cameras. In the corner: Scaled-up flexible magnetic swimmer used in the experiments.

III. MODELING OF THE SWIMMER

A simplified dynamical model of the swimmer is used to simulate the swimmer's displacement under a threedimensional magnetic field. The hydrodynamic effects on the swimmer are approximated by Resistive Force Theory [START_REF] Gray | The propulsion of seaurchin spermatozoa[END_REF] and the shape of the tail is discretized into an articulated chain of N rigid slender rods, as done in the planar swimmer models of [START_REF] Moreau | The asymptotic coarse-graining formulation of slender-rods, bio-filaments and flagella[END_REF][START_REF] Alouges | Self-propulsion of slender micro-swimmers by curvature control: N-link swimmers[END_REF][START_REF] Alouges | Can magnetic multilayers propel artificial microswimmers mimicking sperm cells?[END_REF].

A. Kinematics of the swimmer FIG. 2. Reference and local frames of the discrete-shape model. The swimmer's head frame is oriented relative to the reference frame. For each link i, the corresponding local frame Ri is oriented relative to R h .

We consider a moving frame R head = (O h , e h

x , e h y , e h z ) associated to the head of the swimmer, where O h is the center of the head. The orientation of each link i is represented by the moving frame R i = (O i , e i

x , e i y , e i z ), where O i is the extremity of the i-th link. We call R head ∈ SO(3) the rotation matrix that allows the transformation of coordinates from the fixed reference frame to R head . Similarly, the matrix R i ∈ SO(3), for i = 1 • • • N , denotes the relative rotation matrix that transforms coordinates from R head to R i . We use the angles (θ x , θ y , θ z ) resulting from a (Z -X -Y ) rotation sequence to parametrise R head . Since each link is considered to be slender and axisymmetric, the vectors e i y and e i z are only defined up to one rotation around, hence we parametrize each matrix R i by only two angles (φ i y , φ i z ), resulting from a Z -Y rotation sequence relative to the head frame. With these notations, the swimmer is described by two sets of variables : The 6 Position variables: (X, Θ) where

X = (x h , y h , z h ) ∈ R 3 and Θ = (θ x , θ y , θ z ) ∈ [0, 2π] 3 , (2)
, and the 2N Shape variables, denoted by

Φ = (φ 1 y , φ 1 z , • • • , φ N y , φ N z ). (3)

B. Dynamics of the swimmer

We assume that the swimmer is immersed in an unbounded domain of a viscous fluid. Due to its scale, we consider that the swimmer moves at low Reynolds Number, i.e., that viscosity prevails over inertia, and thus that the fluid is governed by Stokes equations. In what follows, the hydrodynamical interactions between the fluid and the swimmer are further simplified by using the Resistive Force Theory (RFT) framework [START_REF] Gray | The propulsion of seaurchin spermatozoa[END_REF], where the interactions on the global scale between a slender Stokesian swimmer and the surrounding fluid are neglected in favor of the local anisotropic friction of the surface of the slender body with the nearby fluid. This results in explicit expressions of the density of force applied by the fluid to the swimmer that are linear with respect to the velocities, hence to ( Ẋ, Θ, Φ).

From the integration of the local hydrodynamic force densities on each link, we are able to obtain the expression of F h head and F h i , the hydrodynamic drag forces on the head of the swimmer and on each link i ∈ (1, • • • , N ) respectively. The expressions for T h head,P and T h i,P , the moments of these forces about any given point P , are similarly derived from the local drag densities. We consider that the head of the swimmer is magnetized along the e h x axis. Denoting by M the magnetization vector of the head and considering an external homogeneous timevarying field B(t), the following torque is applied to the swimmer:

T mag = M × B(t). (4)
The acceleration terms in the dynamics are neglected due to the Low Reynolds assumption [START_REF] Yates | How microorganisms move through water: The hydrodynamics of ciliary and flagellar propulsion reveal how microorganisms overcome the extreme effect of the viscosity of water[END_REF], thus, the balance of forces and torques applied on the swimmer gives:

             F h head + N i=1 F h i =0 3 , T h head + N i=1 T h i,H = -T mag , (5)
which leads to 6 independent equations. In addition to these equations, the internal contributions of the tail are taken into account by adding the balance of torque on each subsystem consisting of the chain formed by the links i to N for i = (1 • • • N ). These additional 3N equations reduce to 2N non-trivial equations by taking only the components perpendicular to the link k when calculating the sum of the torques from k to N . The elasticity of the tail is discretized by considering a restoring elastic moment T el i at each joint O i that tends to align each pair (i, i + 1) of adjacent links with each other:

T el i = k el e i x × e i-1 x . (6)
Thus, the dynamics of the swimmer are described by the following system of 2N + 6 equations :

                                                       F h head + N i=1 F h i = 0 3 , T h head + N i=1 T h i,H = -T mag , N i=1 T h i,1 .e 1 y = -T el 1 .e 1 y , N i=1 T h i,1 .e 1 z = -T el 1 .e 1 z ,
. . .

T h N,N .e n y = -T el n .e n y , T h N,N .e n z = -T el n .e n z . (7)
Following Resistive Force Theory, the hydrodynamic contributions (left-hand side of the previous system) are linear with respect to the rotational and translational velocities, thus, the previous system can be rewritten matricially in the form :

M h (Θ, Φ)   Ẋ Θ Φ  = B(X, Θ, Φ) , (8)
where M h (Θ, Φ) is a 2N + 6 × 2N + 6 matrix. The left-hand side of the equation represents the hydrodynamic effects on the swimmer and the right-hand side B(X, Θ, Φ) is the magnetic and elastic contributions on the swimmer. We refer the reader to the appendixA for the full derivation of the matrix M h . The previous equation can be rewritten as a control system where the dynamics of the swimmer are affine with respect to the components of the actuating magnetic field viewed as a control:

  Ẋ Θ Φ  = F 0 (Θ, Φ)+ B x (t) B y (t) B z (t)   F 1 (Θ, Φ) F 2 (Θ, Φ) F 3 (Θ, Φ)   , (9)
where the vector fields F 0 , • • • , F 3 are functions of the columns of (M h ) -1 and of the magnetic and elastic constants.

C. Parameter Identification

The hydrodynamic parameters (RFT coefficients) and the elasticity coefficient of the model are identified by non-linear fitting to match the observed velocityfrequency response curve of the experimental swimmer under sinusoidal actuation. Three links were used for the approximation of the tail. Figure 3 shows the agreement between the experimental and simulated frequency responses. The fitted parameter values are given in the appendix I. It is worth noting that using a finer discretization (more than 3) of the tail only marginally improves the fitting error while adding to the computational cost of the magnetic field optimization, as it is shown in figure 9 of the appendix.

IV. MAGNETIC FIELD OPTIMIZATION

In what follows we focus on finding the time-varying magnetic field that maximizes the mean horizontal propulsion speed of the swimmer under the constraint of a periodic deformation. The admissible controls are of the form (B x , B y (t), B z (t)), where the static orientating field B x is fixed and the two dimensional time-varying actuating field (B y (t), B z (t)) is optimized. The bounds on the orientating field and the time-varying fields are taken to be the same as the sinusoidal field used in the experiments (B x = 2.5mT and ||(B y (t), B z (t))|| < 10mT ). The deformation period is chosen to be the same as the observed optimal period for a sinusoidal actuation and it is imposed that the swimmer returns in its initial orientation and position on the y and z-axis at the end of the period. This optimal control problem is numerically solved with the direct solver ICLOCS [START_REF] Falugi | Imperial college london optimal control software (iclocs)[END_REF]. Figure 4 shows the result of the numerical optimization, where we can see the y (a) and z (b) components of the optimal magnetic field. Interestingly enough, this results in a simple magnetic field shape, shown in (c), where the actuating strategy over one deformation period is a partial rotation (about two thirds of a circle) of the magnetic field followed by a full rotation in the opposite direction. Repeating this pattern over time makes the swimmer revolve around the x axis, drawing a "figure-8" in the y -z plane. Figure 5 shows the simulated optimal trajectory of the swimmer compared to the trajectory under the sinusoidal field. The simulated fields perform better (mean propulsion speed of 1.45 10 -3 ms -1 ) than the sinusoidal actuation (mean propulsion speed of 1.2 10 -3 ms -1 ).

V. EXPERIMENTAL RESULTS

The numerical solution of the optimal control problem (Figure 4) is interpolated by its truncated Fourier expansion (first 10 modes) and then implemented in the magnetic generation system to actuate the swimmer. As predicted by the simulations, the optimal field out-performs the sinusoidal field in terms of horizontal propulsion speed, as shown in Figure 6. Figure 7 shows the 3D trajectory of the experimental swimmer during one period. Figure 8 shows the deformation pattern undergone by the swimmer in the two perpendicular planes. Under optimal actuation, the swimmer reaches a mean horizontal propulsion speed of 1.54 ± 0.3 10 -3 ms -1 (N = 6). The mean relative error (∞-norm) between the simulated and observed x-displacement is 0.16(±0.02).

We refer the reader to the supplemental movies [32] for a side-by-side comparison of the displacements of the swimmer actuated by the optimal and sinusoidal field.

FIG. 6. Simulated and experimental horizontal displacements of the swimmer actuated by the optimal magnetic field and the sinusoidal magnetic field for 3 seconds of straight swimming under both actuation patterns. FIG. 7. Trajectory of the experimental swimmer actuated by the optimal magnetic field given in Figure 4 during three periods of the magnetic field.

VI. DISCUSSION

The optimization process exploits two time varying components of the field in order to maximize the horizontal speed of the swimmer, which allows the swimmer to perform a 3D trajectory. This swimming strategy shows the necessity of allowing flagellar swimmers to go out-of plane in order to swim at a maximal propulsion speed, as illustrated in Figure 5 and Figure 7. The effectiveness of non-planar actuation has been corroborated in the literature by studies where non-planar helical waves have been shown to induce a faster propulsion speed for flagellar swimmers. For example, in [START_REF] Islam | Controllable switching between planar and helical flagellar swimming of a soft robotic sperm[END_REF], the sperm-like swimmer's swimming speed increases between 1.2 and 2 times (depending on the viscosity of the fluid) when switching between a planar swimming induced by sinusoidal actuation and helical swimming induced by a conical mag- netic field. This characteristic is also shown in [START_REF] Chwang | A note on the helical movement of micro-organisms[END_REF] for self propelled swimmers. However, our work differs from these approaches as it does not rely on an a priori prescribed actuation pattern or shape deformation but optimizes the 3D driving magnetic field of the swimmer which allows the generation of swimmer-specific optimal actuation. From figure 6, we see that the dynamic model accurately predicts the horizontal displacements of the swimmer under both magnetic fields. It is less accurate in predicting the value of the amplitudes of the oscillations of the swimmer along the y and z axis as can be seen from comparing figures 7 and 5. However, the shapes of the predicted and experimental trajectories match qualitatively. The swimmer's displacements are characterized by the deformations of its tail, as a consequence of the forces and torques balances in equation 5. The experimental deformation patterns of the tail under both magnetic fields are shown in Figure 8, which illustrates how the swimmer breaks the time-reversible symmetry of the Stokes flow to propel in both cases. In particular, the deformation pattern under the optimal field shows a torsion of the tail of the swimmer in addition to a beating pattern. This torsion of the tail causes the swimmer to go out-of-plane. This same effect of the torsion on the trajectory has been shown in the case of flagellated cells that self propel along a 3D chiral path when the torsion of the flagellum is coupled with an oscillating waveform [START_REF] Jikeli | Sperm navigation along helical paths in 3d chemoattractant landscapes[END_REF].

Although the main focus of this work was to optimize the speed of the flexible robot for swimming along a straight line, this swimming strategy is easily implementable for open loop or closed loop path following (see [START_REF] Oulmas | 3d closed-loop motion control of swimmer with flexible flagella at low reynolds numbers[END_REF] for an example) by applying the static component of the magnetic field in the direction tangent to the curve, and the time-varying components in the normal and binormal directions.

VII. CONCLUSION

We have investigated the design of optimal actuation patterns for a flexible low Reynold number swimmer actuated by magnetic fields using a simple computationally inexpensive model that predicts the horizontal displacement of the swimmer. This provides an automated procedure for the optimal control design of flexible magnetic low-Reynolds swimmers. Using this approach, we simulate magnetic fields that maximize the horizontal propulsion speed of the swimmer. From this, we are able to propose a novel magnetic actuation pattern that allows the swimmer to swim significantly faster compared to the usual sinusoidal actuation. This actuation pattern is experimentally validated on a flexible low Reynolds swimmer. The dynamic model used accurately predicts the horizontal displacement of the experimental swimmer under both the optimal and sinusoidal actuation. The simulation and experimental results show that it is necessary to go out-of plane in order to maximize the propulsion speed of flexible magnetic low-Reynolds swimmers. Although the dynamic model used is limited in its accuracy, since it is based on an approximation of the fluid-structure interaction of the swimmer, this study showcases the usefulness of such a simplified model for computationally-inexpensive control design for the actuation of flagellar magnetic swimmers. Current research is focused on finding a solution to allow propulsion of flagellated swimmers in confined environments ( [START_REF] Ishimoto | A study of spermatozoan swimming stability near a surface[END_REF][START_REF] Shum | Hydrodynamic analysis of flagellated bacteria swimming in corners of rectangular channels[END_REF]) or in the presence of obstacles, such as the human vasculature ( [START_REF] Jang | Fabrication and locomotion of flexible nanoswimmers[END_REF]). For this purpose, the technique presented in this paper could be generalized to find magnetic fields that allows the propulsion of the robot in a confined environnement by adapting the model.

The rotation matrix R head represents the coordinates transform between the world frame and R head . It is parametrized by three angles (θ x , θ y , θ z ) resulting from a Z -Y -X rotation sequence :

R head = R x (θ x )R y (θ y )R z (θ x ), ( A1
)
where R x , R y , and R z are the elementary rotation matrices around the x,y, and z axes. Similarly, for each link i, the rotation matrix R i represents the coordinates transform between the head's frame R head and the frame associated with link i R i . It is parametrized by two angles (φ i y , φ i z ) resulting from a Z -Y rotation sequence :

R i = R head Ri R T head for i ∈ (1, • • • , N ). ( A2
)
where Ri = R y (φ i y )R z (φ i z ). For all skew-symmetric matrix A ∈ R 3×3 , we can define a vector Ω verifying :

∀V ∈ R 3 AV = Ω × V . (A3)
Using this property, we define the following angular velocity vectors :

• Ω head is the vector defining the cross product associated with the skew-symmetric matrix Ṙhead R T head . Ω head depends linearly on ( θx , θy , θz ) as follows :

Ω head = L head   θx θy θz   , (A4)
where

L head =   1 0 sin(θ y ) 0 cos(θ x ) -cos(θ y ) sin(θ x ) 0 sin(θ x ) cos(θ x ) cos(θ y )     θx θy θz   (A5)
• Ω i is the vector defining the cross product associated with the skew-symmetric matrix Ṙi Ri

T for i ∈ (1 • • • N )
. Ω i depends linearly on ( φi y , φi z ) as follows :

Ω i = L i   φi x φi y φi z   , (A6)
where

L i =   sin(φ i y ) 0 cos(φ i y ) 0 0 1     φi x φi y φi z   . (A7)
2. Hydrodynamical force and torque on the head of the swimmer

We consider a drag force acting on the head of the swimmer that is proportional to its velocity in each direction of the head frame's R head and a resistive torque proportional to the angular velocity of the head :

F h head = -R head   k H, 0 0 0 k H,⊥ 0 0 0 k H,⊥   R T head Ẋ , T h head = -k R Ω head = -k R L head Θ, ( A8
)
where k H, and k H,⊥ are the parallel and perpendicular hydrodynamic coefficients of the head and k R is a rotationnal drag coefficient.

Expression of the hydrodynamical force density on a point of the tail of the swimmer

For i = (1 • • • N ), we consider a point x i (s) on the i-th link of the tail parametrized by its arclength s such as

x i (s) = O i -se i x . (A9)
Using the rotation matrices defined in the previous paragraph, x i (s) is written as :

x i (s) = O H -rR head e x -l i-1 k=1 R head Rk e x -sR head Ri e x , (
A10) where r is the distance between O H and O i . We differentiate the previous equation to obtain the expression of the velocity of x i (s):

ẋi (s) = Ẋ -r Ṙhead e x -l i-1 k=1 Ṙhead Rk e x -l i-1 k=1 R head Ṙk e x -s Ṙhead Ri e x -sR head Ṙi e x .
(A11) Following Resistive Force Theory, the density of hydrodynamic force f i is linear with respect to the components of ẋi (s):

f i (s) = -k ( ẋi (s).e i x )e i x -k ⊥ ( ẋi (s).e i y )e i y -k ⊥ ( ẋi (s).e i z )e i z =S i ẋi (s) , (A12)
where k and k ⊥ are respectively the parallel and prependicular drag coefficients of the swimmer's tail and, for each link i :

S i = (R head Ri )D(R head Ri ) T . (A13)
where

D = -   k 0 0 0 k ⊥ 0 0 0 k ⊥   (A14)
Using the expression of ẋi (s) ( equation (A11)) in equation (A12), the hydrodynamic force density reads:

f i (s) =S i Ẋ -r(R head Ri D RT i )R T head Ṙhead e x -l(R head Ri D RT i ) i-1 k=1 R T head Ṙhead Rk e x -l(R head Ri D RT i ) i-1 k=1 Rk RT k Ṙk e x -s(R head Ri D RT i )R T head Ṙhead Ri e x -s(R head Ri D RT i ) Ṙi e x . (A15)
Using the definition of the angular velocity vectors (equations A4,A6 ) , we rewrite the expression above as :

f i (s) =S i Ẋ -r(R head Ri D RT i )Ω head × e x -l(R head Ri D RT i ) i-1 k=1 Ω head × ( Rk e x ) -l(R head Ri D RT i ) i-1 k=1 Rk Ω k × e x -s(R head Ri D RT i )Ω head × ( Ri e x ) -s(R head Ri D)Ω i × e x .

(A16)

In order to write cross products in matricial form, we introduce the following notation :

∀V = V 1 V 2 V 3 T ∈ R 3 [V ] × =   0 -V 3 V 2 V 3 0 -V 1 -V 2 V 1 0   .
(A17) Using this notation and the linear dependency of Ω head and Ω i on the angular velocities (equations A4 and A6), the hydrodynamical drag force density on link i is written as a linear function of Ẋ, Ω, and Φ :

f i (s) = A i (s)   Ẋ Θ Φ  , (A18)
where A i (s) ∈ R 3,2N +6 is defined block-wise as :

A i (s) =        A i X ∈ R 3×3 A i θ ∈ R 3×3 A i 1 ∈ R 2×3 . . . A i N ∈ R 2×3        , (A19)
where Using this notation, the hydrodynamical force on link i depends linearly on the position and shape velocities as follows :

                                         A i X (s) =S i , A i θ (s) =(
F h i = B i   Ẋ Θ Φ  . (A22)
The hydrodynamic torque on link i calculated about the point O H is computed from the hydrodynamical force density as follows : (A25) Thus, the matricial form of T h i,H is :

T h i,H = -([rR head e x + i-1 k=1 lR head Rk e x ] × B i )   Ẋ Θ Φ  -([R head Ri ] × C i )   Ẋ Θ Φ  (A26) .
Similarly, for k = (i, • • • , N ), the hydrodynamic torque on link i about O k reads :

T h i,j = -( i-1 k=j [lR head Rk e x ] × B i -[R head Ri ] × C i )   Ẋ Θ Φ  .
(A27)

Equations of motion

Using the expressions of the hydrodynamical forces and torques in the previous section (equations A22,A26, and A27), the left-hand side of the dynamical system (7) is written as a linear function of the state derivatives :

                    N i=1 F h i + F h head T h head + N i=1

FIG. 3 .

 3 FIG. 3. Experimental (N = 6) and simulated horizontal swimming velocities for the frequency range f = (0 • • • 3Hz).

FIG. 4 .

 4 FIG. 4. Actuating magnetic field that maximizes the horizontal speed of the swimmer. (a) and (b) are respectively the y and z components of the magnetic field. (c) : Shape of the actuating optimal magnetic field deriving from the optimization process during one period.

FIG. 8 .

 8 FIG. 8. Deformation pattern of the experimental swimmer actuated by the optimal magnetic field and the sinusoidal magnetic field over one period in two planes. ( (a) : Optimal field, top view, (b) : Sinusoidal field, top view, (c): Optimal field, side view, (d) : Sinusoidal field, side view.). The snapshots of the experimental swimmer are taken at equal time steps over a period of the actuating fields. We observe (as expected) no deformation in the side plane for the sinusoidal field.

4 .

 4 Hydrodynamical forces and torques on each linkWe define the two 3 × 2N + 6 matrices B i and C i as follows :

1 k=1lR head Rk e x ] × l 0 f

 10 (s) -O H ) × f i (s)ds. (A23)We rewrite (x i (s) -O H ) as :(x i (s) -O H ) = -rR head e x -i-1 k=1 lR head Rk e x -sR head Ri e x .(A24)T h i,H is then expressed as :T h i,H = -[rR head e x + ii (s)ds -[R head Ri ] × l 0 sf i (s)ds.

  R head Ri D RT i )[re x + l Rk e x ] × L head + (R head Ri D RT i )[s Ri e x ] × L head , A i 1 (s) =l(R head Ri D Rt i ) R1 [e x ] × L 1 , . . . (s) =l(R head Ri D Rt i ) Ri-1 [e x ] × L i-1 , A i i (s) =s(R head Ri D)[e x ] × L i , A i j =0 3 ∀j > i .

	i-1
	k=1
	A i i-1 (A20)

ACKNOWLEDGMENTS

The authors were partially funded by the CNRS through the project Défi InFIniti-AAP 2018-2019.

Appendix A: Formulation of the hydrodynamics of the swimmer.

1. Parametrization of the rotation matrices and angular velocities.

,

where

) where R * i is the 2 × 3 matrix consisting of the second and third line of the matrix (R head Ri ) T which represents the projection on the plane (e i y , e i z ). Figure 9 shows the result of the parameter fitting (relative l 2 -norm) between the experimental velocityfrequency response curve and the simulated velocityfrequency response curve.

As seen in the figure, the model fails to accurately match the experimental velocity-frequency curve with N = 1 and N = 2 links. The fitting results in a small error for N ≥ 3.