Ilona Bluemke
email: i.bluemke@ii.pw.edu.pl

Marcin Tarka

Learning algorithms in the detection of unused functionalities in SOA systems

The objective of this paper is to present an application of learning algorithms to the detection of anomalies in SOA system. As it was not possible to inject errors into the "real" SOA system and to analyze the effect of these errors, a special model of SOA system was designed and implemented. In this system several anomalies were introduced and the effectiveness of algorithms in detecting them were measured. The results of experiments can be used to select efficient algorithm for anomaly detection. Two algorithms: K-means clustering and Kohonen networks were used to detect the unused functionalities and the results of this experiment are discussed.

Introduction

With the growth of computer networking, electronic commerce, and web services, security of networking systems has become very important. Many companies now rely on web services as a major source of revenue. Computer hacking poses significant problems to these companies, as distributed attacks can make their systems or services inoperable for some period of time. As this happens often, an entire area of research, called Intrusion Detection, is devoted to detect these activities.

Nowadays many system are based on Service Oriented Architecture (SOA) [1,2] idea. A system based on SOA provides functionalities as a suite of interoperable services that can be used within multiple, separate systems from several business domains. SOA also provides a way for consumers of services, such as web-based applications, to be aware of available SOA-based services. Service-orientation requires loose coupling of services with operating systems, and other technologies that underly applications. SOA separates functionality into distinct units, or services, which developers make accessible over a network in order to allow users to combine and reuse them in the production of applications.

The objective of this paper is to present the detection of anomalies in SOA systems by learning algorithms. Related work is presented in section 2 and in section 3, a special model of SOA system which was used in experiments, is presented. In this systems several anomalies were introduced. Four algorithms: Chi-Square statistics, k-means clustering, emerging patterns and Kohonen networks were used to detect anomalies. Detection of anomalies by k-means and Kohonen networks is presented in section 4 and some conclusions are given in section 5.

Related work

Many anomaly detection algorithms have been proposed in the literature. They differ according to the information used for analysis and according to techniques that are used to detect deviations from normal behavior. Lim and Jones in [START_REF] Lim | Network Anomaly Detection System: The State of Art of Network Behavior Analysis[END_REF] proposed two types of anomaly detection techniques based on employed techniques: the learning model method and the specification model. The learning approach is based on the application of machine learning techniques, to automatically obtain a representation of normal behaviors from the analysis of system activities. The specification-based approach requires that someone manually provides specifications of correct behavior. Approaches that concern the model construction are presented in Fig. 1.

Fig. 1. Taxonomy of anomaly detection behavioral model (based on [3])

The specification approach depends more on human observation and expertise than on mathematical models. It was first proposed by C. Ko et. al. [START_REF] Ko | Execution monitoring of security-critical programs in distributed systems: a specification-based approach[END_REF] and uses a logic based description of expected behavior to construct a base model. This specificationbased anomaly detector monitors multiple system elements, ranging from application to network traffic.

In the protocol based approach [START_REF] Lemonnier | Protocol Anomaly Detection in Network-based IDSs[END_REF] a normal use model is built from the protocol specification e.g. TCP/IP. Lemonnier [START_REF] Lemonnier | Protocol Anomaly Detection in Network-based IDSs[END_REF] proposed a protocol anomaly filter able to specifically analyze a protocol and model the normal usage of a specific protocol. This technique can be seen as a filter looking for protocol misuse. The protocol could be interpreted as any official set of rules describing the interaction between the elements of a computer system.

Many protocol anomaly detectors are built as state machines [START_REF] Sekar | Specificationbased anomaly detection: A New Approach for Detecting Network Intrusions[END_REF]. Each state corresponds to a part of the connection, such as a server waiting for a response from client. The transitions between the states describe the legal and expected changes between states. Besides, Z. Shan et. al. [START_REF] Shan | A Network State Based Intrusion Detection Model[END_REF] uses network state based model approach to describe intrusion attacks.

In the transaction based approach the "positive" behavior is formally described. The desired actions and sequence of actions are specified by the definition of transactions. Such explicit definition makes the transaction an integral part of security policy. In the research proposed by R. Buschkes et. al. [START_REF] Buschkes | Transaction-based Anomaly Detection[END_REF] the detection of anomalies is based on the definition of correct transactional behavior.

The learning model must be trained on the specific network. In the training phase, the behavior of the system is observed and logged and machine learning techniques are used to create a profile of normal behaviors. In the process of creating an effective anomaly detection model: rule-based, model-based, and statistical-based approaches have been adopted to create the baseline profile.

Rule-based systems used in anomaly detection describe the normal behavior of users, networks and/or computer systems by a set of rules. These predefined rules typically look for the high-level state change patterns observed in the audit data. SRI International's Next-generation Intrusion Detection Expert System (NIDES) [START_REF] Anderson | Valdes: A Next-generation Intrusion Detection Expert System (NIDES) Summary[END_REF] is an innovative statistical algorithm for anomaly detection, and an expert system that encodes known intrusion scenarios. The features considered by NIDES are related to user activity, for instance ,CPU and file utilization. The information collected by NIDES is compared with the long -term profile of analyzed system and used for the learning process based on Chi-squared statistic. Owens et. al. present in [START_REF] Owens | An adaptive expert system approach for intrusion detection[END_REF] an adaptive expert system for intrusion detection based on fuzzy sets.

In the model based anomaly detector the intrusions are analyzed at a higher level of abstraction than the audit records of the rule based approach. Execution is restricted to pre-computed patterns of expected behavior. Different types of models are used to characterize the normal behavior of the monitored system like data mining, neural networks, pattern matching, etc.

Data mining extracts the behavioral models from a large amount of data. Intrusion detection systems require frequent adaptation to resist new attacks, so data mining techniques that can adaptively build new detection models are very useful. An example of data mining system was proposed by Lee et. al. and is presented in [START_REF] Lee | Data mining approaches for intrusion detection[END_REF]. The key idea is to mine network audit data, then use the patterns to compute inductively learned classifiers that can recognize anomalies and known intrusions.

A neural network learns the user profile in the training process. Training may use data more abstract than the audit records. The fraction of incorrectly predicted events constitutes the variance of the user behavior. [START_REF] Bivens | Network-based intrusion detection using neural networks[END_REF] describes anomalies detection in information system based on Kohonen networks. There are several libraries supporting building neural networks e.g. [13 -15]. Other neural networks based systems for intrusion detection are also described in [START_REF] Ryan | Intrusion Detection with Neural Networks[END_REF][START_REF] Ghosh | A Study in Using Neural Networks for Anomaly and Misuse Detection[END_REF][START_REF] Han | Evolutionary Neural Networks for Anomaly Detection Based on the Behaviour of a Program[END_REF][START_REF] Bivens | Network-based intrusion detection using neural networks[END_REF].

In the pattern matching approach, learning is used to build a traffic profile for a given network. Traffic profiles are built using features such as packet loss, link utilization, number of collisions. Normal behavior is captured as templates and tolerance limits are set based on different levels of standard deviation. The usage of emerging patterns in anomaly detection is described in [START_REF] Ceci | Discovering Emerging Patterns for Anomaly Detection in Network Connection Data[END_REF].

First statistical based anomaly detection was proposed by Denning and Neumann [START_REF] Denning | Requirements and Model for IDES-A Real-Time Intrusion-Detection Expert System[END_REF] in 1985. The anomaly detector observes subjects and generates profiles for them that represent their behavior. The widely known techniques in statistics can often be applied; e.g. data points that lie beyond a multiple of the standard deviation on either side of the mean might be considered anomalous. There are many statistical techniques like Bayesian statistics, covariance matrices and Chi-square statistics [START_REF] Masum | Chi-square statistical profiling for anomaly detection[END_REF] for the profiling of anomaly detection. Example of Chi-square statistic in anomaly detection is described in [START_REF] Ye | An anomaly detection technique based on a chi-square statistic for detecting intrusions into information systems[END_REF]. Statistical approaches disadvantage is the insensitivity to the order of occurrence of events. Sequential interrelationships among events should be considered for more accurate detection. It is also difficult to determine a threshold above which an anomaly should be considered.

Commercial products applying rule based, statistical based, model based and neural networks approach to detected anomalies are briefly described in [START_REF] Lim | Network Anomaly Detection System: The State of Art of Network Behavior Analysis[END_REF].

Research model

As we were not allowed to use real SOA system and inject anomalies into it, a special system was implemented. The business idea of this system VTV (Virtual TV) is presented in Fig. 2. The exemplary company is a virtual TV provider. The company does not have its own technical infrastructure and is associating TV digital provider with the telecommunication operator. The receiving equipment is delivered to the client by one of courier companies. The company is also using two applications: Customer Relationship Management (CRM) containing all clients data and a storage management system.

The VTV system simulates a real SOA system and enables to inject typical anomalies into its regular operation. Configurable are frequencies of : single services calls, group of services calls, processes calls, and services in a context. More details can be found in [START_REF] Tarka | Anomaly detection in SOA systems[END_REF].

Fig. 2. The business relations of an exemplary company

The architecture of VTV system is presented in Fig. 3. The request generator simulates activities of clients. It generates different types of request (e.g. create, modify, deactivate) for available services (e.g. TV or hardware). Depending of the value of request some requests can be identified as very important. The generated address of a client determines the currier group. The configurable parameters of the request generator enable to simulate real operation and to inject some anomalies: e.g. request for courier, hardware, TV services can be set, ratio of create, modify, deactivate requests in generated requests can be chosen.

Fig. 3. Architecture of VTV system

The generated request is transferred to the business processes engine (Fig. 3) which composes processes from services available on Enterprise Service Bus (ESB). Outputs of processes are logs. Logs of model contain information similar to logs from monitoring real SOA systems.

Environment of experiment

All experiments were conducted on PC with Intel Core 2 Duo T7200, 512 Mb of memory under Fedora Core operating system. The examined algorithms are using text input files prepared from logs of the VTV system (Fig. 3). In this log file information like number of requests, name of processes, name of services called, execution time are written. These logs files are then transformed by scripts, implemented in R [START_REF]The R Project for Statistical Computing[END_REF] environment, into transactions or summarized reports which are input to detection algorithms. The flow of data from logs of the model is shown in Fig. 4.

Fig. 4. Flow of data from logs

The implementation of algorithms for anomaly detection k-means, Chi-Square and Kohonen networks algorithms was made in R environment while the emerging patterns algorithm, which is more complicated, was implemented in C++.

k-means algorithm was prepared based on [START_REF] Munz | Traffic Anomaly Detection Using K-Means Clustering[END_REF], for emerging patterns algorithm information from [START_REF] Ceci | Discovering Emerging Patterns for Anomaly Detection in Network Connection Data[END_REF][START_REF] Guozhu | Efficient Mining of Emerging Patterns: Discovering Trends and Differences[END_REF] were used, Chi-Square detector was taken from [START_REF] Ye | An anomaly detection technique based on a chi-square statistic for detecting intrusions into information systems[END_REF].

Experiment

The research system presented in section 3 was used to explore four cases typical for SOA systems i.e.: change in the frequency of service calls, change in the frequency of a group of services, change of the context of services calls and not used functionalities.

For each of the above listed cases it was expected, that anomalies detector provides information useful for the maintenance team. Four learning algorithms (section 2) were used in the anomalies detector: Chi-square statistic, Kohonen network, Emerging patterns and k-means clustering.

Each of the above algorithms represents different approach to anomalies detection. The goal of the experiment was to examine advantages and disadvantages of each of these algorithms. Anomalies detection is a kind of clustering with two types of clusters grouping normal and abnormal behaviors. Correctly identified anomaly is an abnormal event which was introduced by purpose by case' scenario and was assigned to the abnormal cluster. Identified as anomalies other events or not recognized anomalies are treated as errors. The values measured in experiments are:  FP (false positive) -number of incorrectly identified anomalies,  FN (false negative) -number of not recognized anomalies,  TP (true positive) -number of correctly identified normal behaviors,  TN (true negative) -number of incorrectly identified normal behaviors.

Good anomalies detector should generate small number of errors. To compare the quality of detectors often sensitivity and specificity measures are used. The sensitivity and specificity are defined accordingly as:

(1) and [START_REF]SOA manifesto[END_REF] The relation between specificity and sensitivity -ROC (Receiver Operating Characteristics) curve [START_REF] Hanley | Receiver operating characteristic (ROC) methodology: the state of the art[END_REF] can be used to compare the quality of models. ROC as a technique to analyze data was introduced during the second world war to identify if signal seen on a radar were coming from an enemy, an alliance or if it were noise. Currently ROC curves are used to analyze different types of data e.g. radiological data.

The construction of ROC curves during experiments was as follows:

FN TP TP y sensitivit   FP TN TN y specificit  
1. create entities with algorithm' parameters 2. for each entity :

 perform experiment  calculate specificity and sensitivity  mark the point on the diagram 3. draw the line connecting points. For each examined algorithm also the learning time was calculated.

Plan for experiment

The examination of anomaly detection algorithms was based on four test cases: changes in frequency of service and groups of services calls, change of the context of services calls, and lacking functionalities. Each of these cases simulates one type of anomaly typical for SOA.

The experiment was conducted in following steps: Below the results for the scenario "not used functionalities" are presented. The detection of this kind of anomalies is very important in distributed systems especially in SOA systems. Some functions may be not used as a result of errors in routing algorithm or in business logic. If not used services were dedicated to some cluster this cluster could be used for other purposes. In ours research' environment this anomaly was obtained by forcing the request generator not to generate courier delivery services by setting parameter Courier_share to zero [START_REF] Tarka | Anomaly detection in SOA systems[END_REF]. In [START_REF] Bluemke | Detection of anomalies in SOA system by learning algorithms[END_REF] the detection of anomalies in the frequency of service calls by K-means clustering algorithm and emerging patterns are described.

Results of experiments

The goal of examination was to find high level of detection with minimal number of false alarms. If the ideal detection was not possible preferred were the results with no false alarms. In practice, if system is generating many false alarms the user will neglect any alarm. [START_REF] Munz | Traffic Anomaly Detection Using K-Means Clustering[END_REF] is a clustering analysis algorithm that groups objects based on their feature values into k disjoint clusters. Objects that are classified into the same cluster have similar feature values. k is a positive integer number specifying the num-ber of clusters, and has to be given in advance. All objects are assigned to their closest cluster according to the ordinary Euclidean distance metric.

k-means algorithm k-means clustering

At the beginning the summarized reports used by the algorithm were created for all logs in the system. The results are given in Table 1. The maximal distance from centroid of cluster for a cluster member is denoted as maxL. Initially the experiment was conducted for 150 logs in training set. The anomaly wasn't satisfactorily detected, the numbers of not identified anomalies were high (e.g. 23, first row in Table 1), also the numbers of incorrectly identified anomalies were high (e.g. 14 first row in Table 1). In the training data there were many summarization reports in which the number of services not called in regular data, was close to zero, so the lack of calls of these services was treated as normal behavior. Next, the number of logs was increased to 350. The results improved (rows 8-12 Table 1).

In Fig. 5 the ROC curve is shown for 350 training logs. The high change in sensitivity is seenfrom zero to one. Point A was obtained by decreasing the value of maxL parameter from the value in point B, till the level in which zero was not assigned to the group with the minimal centroid.

The k-means clustering algorithm can be used in the detection of not used functionalities if the number of training data is high. In [START_REF] Bluemke | Detection of anomalies in SOA system by learning algorithms[END_REF] we have shown that this algorithm was not suitable in detecting the change in frequencies of single services calls.

Kohonen networks

A neural network is a set of simple, highly interconnected units called neurons. Neurons transform a set of inputs to a set of desired outputs. The result of the transformation is determined by the characteristics of the elements and the weights associated with the interconnections among them, therefore by adjusting the weight the output can be controlled. The process of updating the weights and thresholds is called learning.

Training usually takes the form of presenting the network with set of typical inputs vectors. In unsupervised learning mode (self -organizing map [START_REF] Kohonen | The self-organizing map[END_REF] -Kohonen networks), the inputs are presented for training; the adjustments are made so that the network is able to recognize inputs as belonging to a particular profile.

The training utilizes competitive learning. When a training example is fed to the network, its Euclidean distance to all weight vectors is computed. The neuron whose weight vector is most similar to the input is called the best matching unit (BMU). The weights of the BMU and neurons close to it in the lattice are adjusted towards the input vector. The magnitude of the change decreases with time and with distance (within the lattice) from the BMU. This process is repeated for each input vector for a number of cycles. The network winds up associating output nodes with groups or patterns in the input data set.

The results for Kohonen network algorithm are shown in Table 2. Column Lambdatoleration ration, shows the restraint to the distance from BMU. In rows 1-5 the numbers of not recognized (FN) and of incorrectly identified anomalies (FP) are high. After increasing the number of elements in the training phase the results of detection significantly improved. In rows 6-11 only few anomalies were not recognized. In Fig. 6 the ROC curve for 350 training logs is shown. Point A is optimal, with the sensitivity= 0.92 and specificity= 0.94 and the costs of false qualification of anomaly and false qualification of normal behavior are equal. These values may be improved by increasing the training data. Point B represents the case in which all anomalies were detected, the specificity was equal to 0.84.

Fig. 6. ROC curve for Kohonen networks algorithm with 350 training logs

The above presented data show that Kohonen networks can be used to detect the unused functionalities but large training data must be provided.

Conclusions

In this paper some results of the detection of one anomalyunused functionalities in SOA system are presented. An experiment was conducted in an environment designed to introduce several types of anomalies. This environment, described in section 3, enables the detection of anomalies by four learning algorithms, from different types (section 2): emerging patterns, k-means clustering, Kohonen networks and statistical Chi-Square . In this paper the results of only two: k-means clustering and Kohonen networks are presented. Both algorithms were able to satisfactorily detect unused functionalities while two others i.e. Chi-Square and emerging patterns appeared to be inappropriate in the detection of unused functionalities.

In [START_REF] Bluemke | Detection of anomalies in SOA system by learning algorithms[END_REF] the results of detecting other anomaly -the change in frequencies of service calls was described. The least accurate in the detection of this kind of anomaly was the k-means clustering algorithm and the best was emerging pattern algorithm. Kohonen algorithm also produced quite good results.

The results presented in this paper and in [START_REF] Bluemke | Detection of anomalies in SOA system by learning algorithms[END_REF] show that in anomalies detection in SOA systems different algorithms may appear most suitable for different type of anomaly so further research should be conducted. The exemplary SOA system (section 2) enables to conduct other experiments examining the suitability of learning algorithms in the detection of other anomalies. The results of these experiments will be available soon.

Fig. 5 .

 5 Fig. 5. ROC curve for k-means algorithm with 350 training logs

Table 1 .

 1 Results for clustering algorithm.

	row k	maxL	logSize TP TN FP FN	Sensitivity Specificity
	1	4	1	150	0	1	14 23	1	0.6
	2	4	2	150	0	3	12 23	0	0.2
	3	4	2.3	150	23 9	6	0	1	0.6
	4	4	2.5	150	23 3	12 0	1	0.2
	5	5	2	150	23 5	10 0	1	0.33
	6	5	2.3	150	23 6	9	0	1	0.4
	7	5	2.5	150	23 8	7	0	1	0.53
	8	4	1	350	10 4	3	0	1	0.57
	9	4	1.5	350	10 5	2	0	1	0.71
	10 4	2	350	0	5	2	10	0	0.71
	11 4	2.5	350	0	5	2	10	0	0.71
	12 4	3	350	0	6	1	10	0	0.86

Table 2 .

 2 Results for Kohonen network algorithm.

	row N Lambda Logsize TP	TN FP FN	Sensitivity Specificity
	1	25 0.5	150	358 331 137 398	0.47	0.707
	2	25 1	150	324 366 102 432	0.43	0.782

Acknowledgments. We are very grateful to the reviewers for many valuable remarks.


Joris Van Der Hoeven

Grégoire Lecerf
email: lecerf@lix.polytechnique.fr

Modular composition via complex roots

Modular composition is the problem to compute the composition of two univariate polynomials modulo a third one. For polynomials with coefficients in a finite field, Kedlaya and Umans proved in 2008 that the theoretical complexity for performing this task could be made arbitrarily close to linear. Unfortunately, beyond its major theoretical impact, this result has not led to practically faster implementations yet. In this paper, we explore the particular case when the ground field is the field of computable complex numbers. Ultimately, when the precision becomes sufficiently large, we show that modular compositions may be performed in softly linear time.

INTRODUCTION

Let 𝕂 be an effective field, and let f , g, h be polynomials in 𝕂 [x]. The problem of modular composition is to compute g ∘ f modulo h. Modular composition is an important problem in complexity theory because of its applications to polynomial factorization [START_REF] Kaltofen | Subquadratic-time factoring of polynomials over finite fields[END_REF][START_REF] Kedlaya | Fast modular composition in any characteristic[END_REF][START_REF] Kedlaya | Fast polynomial factorization and modular composition[END_REF]. It also occurs very naturally whenever one wishes to perform polynomial computations over 𝕂 inside an algebraic extension of 𝕂. In addition, given two different representations 𝕂[x] / (h(x)) ≅ 𝕂[x ˜]/(h ˜(x ˜)) of an algebraic extension of 𝕂, the implementation of an explicit isomorphism actually boils down to modular composition.

Denote by M(n) the number of operations in 𝕂 required to multiply two polynomials of degree <n in 𝕂 [x]. Let f , g and h be polynomials in 𝕂[x] of respective degrees <n, <n and n. The naive modular composition algorithm takes O(n M(n)) operations in 𝕂. In 1978, Brent and Kung [START_REF] Brent | Fast algorithms for manipulating formal power series[END_REF] gave an algorithm with cost O( n √ M(n) + n 2 ). It uses the baby-step giant-step technique due to Paterson and Stockmeyer [START_REF] Paterson | On the number of nonscalar multiplications necessary to evaluate polynomials[END_REF], and even yields a sub-quadratic cost O(n 𝜛 + n √ M 𝕂 (n)) when using fast linear algebra (see [13, p. 185]). The constant 𝜛 > 1.5 is such that a n √ × n √ matrix over 𝕂 can be multiplied with another n √ × n rectangular matrix in time O(n 𝜛 ). The best current bound 𝜛 < 1.6667 is due to Huang and Pan [START_REF] Huang | Fast rectangular matrix multiplication and applications[END_REF]Theorem 10.1].

A major breakthrough has been achieved by Kedlaya and Umans [START_REF] Kedlaya | Fast modular composition in any characteristic[END_REF][START_REF] Kedlaya | Fast polynomial factorization and modular composition[END_REF] in the case when 𝕂 is the finite field 𝔽 q . For any positive 𝜀 > 0, they showed that the composition g ∘ f modulo h could be computed with bit complexity O((n log q) 1+𝜀 ). Unfortunately, it remains a major open problem to turn this theoretical complexity bound into practically useful implementations.

Quite surprisingly, the existing literature on modular composition does not exploit the simple observation that composition modulo a separable polynomial h ∈ 𝕂[x] that splits over 𝕂 can be reduced to the well known problems of multi-point evaluation and interpolation [START_REF] Zur Gathen | Modern computer algebra[END_REF]Chapter 10]. More precisely, assume that h = (x -𝜎 1 ) ⋯ (x -𝜎 n ) is separable, which means that gcd(h, h') = 1. If f , g ∈ 𝕂[x] are of degree <n, then g ∘ f mod h can be computed by evaluating f at 𝜎 1 , …, 𝜎 n , by evaluating g at f (𝜎 1 ), …, f (𝜎 n ), and by interpolating the evaluations (g ∘ f )(𝜎 1 ), …, (g ∘ f )(𝜎 n ) to yield g ∘ f mod h.

Whenever 𝕂 is algebraically closed and a factorization of h is known, the latter observation leads to a softly-optimal algorithm for composition modulo h. More generally, if the computation of a factorization of h has a negligible or acceptable cost, then this approach leads to an efficient method for modular composition. In this paper, we prove a precise complexity result in the case when 𝕂 is the field of computable complex numbers. In a separate paper [START_REF] Van Der Hoeven | Modular composition via factorization[END_REF], we also consider the case when 𝕂 is a finite field and h has composite degree; in that case, h can be factored over suitable field extensions, and similar ideas lead to improved complexity bounds.

In the special case of power series composition (i.e. composition modulo h= x n ), our approach is similar in spirit to the analytic algorithm designed by Ritzmann [START_REF] Ritzmann | A fast numerical algorithm for the composition of power series with complex coefficients[END_REF]; see also [START_REF] Van Der Hoeven | Fast composition of numeric power series[END_REF]. In order to keep the exposition as simple as possible in this paper, we only study composition modulo separable polynomials. By handling multiplicities with Ritzmann's algorithm, we expect our algorithm to extend to the general case.

The organization of the present paper is as follows. In section 2, we specify the complexity model to be used, and various standard notations. In section 3, we give a detailed version of the modular composition algorithm that we sketched above for a separable modulus that splits over 𝕂. In order to instantiate this algorithm for the field ℂ com of computable complex numbers, we need additional concepts. In section 4, we recall basic results about ball arithmetic [START_REF] Van Der Hoeven | Ball arithmetic[END_REF]. In section 5, we recall the computation model of straight-line programs [START_REF] Bürgisser | Algebraic complexity theory[END_REF]. In section 6, we introduce a new ultimate complexity model that is convenient for proving complexity results at a "sufficiently large precision". This model has the advantage that complexity results over an abstract effective field 𝕂 can naturally be turned into ultimate complexity results over ℂ com . In section 7, we apply this transfer principle to the modular composition algorithm from section 3-we expect our framework to be useful in many other situations.

One disadvantage of ultimate complexity analysis is that it does not provide us with any information about the precision from which the ultimate complexity is reached. In practical applications, the input polynomials f , g and h often admit integer or rational coefficients. In these cases, the required bit precision is expected to be of order n (l + n) in the worst case, where n = deg h and l is the largest bit size of the coefficients: in fact, this precision allows to compute all the complex roots of h efficiently using algorithms from [START_REF] Neff | An efficient algorithm for the complex roots problem[END_REF][START_REF] Pan | Univariate polynomials: nearly optimal algorithms for numerical factorization and root-finding[END_REF][START_REF] Schönhage | The fundamental theorem of algebra in terms of computational complexity[END_REF]. This precision should also be sufficient to perform the multi-point polynomial evaluations of g and f by asymptotically fast algorithms. We intend to work out more such detailed bit complexity bounds for this situation in a forthcoming paper.

PRELIMINARIES

In the sequel, we consider both the algebraic and bit complexity models for analyzing the costs of our algorithms. The algebraic complexity model expresses the running time in terms of the number of operations in some abstract ground ring or field [START_REF] Bürgisser | Algebraic complexity theory[END_REF]Chapter 4]. The bit complexity model relies on Turing machines with a sufficient number of tapes [START_REF] Papadimitriou | Computational Complexity[END_REF].

Fundamental algebraic complexity bounds. Let 𝕂 be an effective field. We write M: ℕ → ℝ > for a function that bounds the total cost of a polynomial product algorithm in terms of the number operations in 𝕂. In other words, two polynomials of degrees n in 𝕂[x] can be multiplied using at most M(n) arithmetic operations in 𝕂. The schoolbook algorithm allows us to take M(n) = O(n 2 ). The fastest currently known algorithm [START_REF] Cantor | On fast multiplication of polynomials over arbitrary algebras[END_REF] yields 1) g(n) (we refer the reader to [6, Chapter 25, Section 7] for technical details). In order to simplify the cost analysis of our algorithms we make the customary assumption that M(n 1 ) / n 1 ⩽ M(n 2 ) / n 2 for all 0 < n 1 ⩽ n 2 . Notice that this assumption implies the super-additivity of M, namely M(n 1 ) + M(n 2 ) ⩽ M(n 1 + n 2 ) for all n 1 ⩾ 0 and n 2 ⩾ 0.

M(n) = O(n log n log log n) = O ˜(n). Here, the soft-Oh notation f (n) ∈ O ˜(g(n)) means that f (n) = g(n) log O(
Fundamental bit complexity bounds. For bit complexity analyses, we consider Turing machines with sufficiently many tapes. We write I(n) for a function that bounds the bit-cost of an algorithm which multiplies two integers of bit sizes at most n, for the usual binary representation. The best known bound [START_REF] Harvey | Even faster integer multiplication[END_REF] for I(n) is O(n log n 8 log * n ) = O ˜(n). Again, we make the customary assumption that I(n)/n is nondecreasing.

Multipoint evaluation and interpolation. Let 𝕂 again be an effective field. The remainder (resp. quotient) of the Euclidean division of g by h in 𝕂[x] is denoted by g rem h (resp. by g quo h). It may be computed using O(M(n)) operations in 𝕂, if g and h have degrees ⩽n. We recall that the gcd of two polynomials of degrees at most n over 𝕂 can be computed using Algorithm 11.4]. Given polynomials f and g 1 , …, g l over 𝕂 with deg f = n and deg g 1 + ⋯ + deg g l = O(n), all the remainders f rem g i may be computed simultaneously in cost O(M(n) log l) using a subproduct tree [START_REF] Zur Gathen | Modern computer algebra[END_REF]Chapter 10]. The inverse problem, called Chinese remaindering, can be solved with a similar cost O(M(n) log l), assuming that the g i are pairwise coprime. The fastest known algorithms for these tasks can be found in [1, 2, 10].

O(M(n) log n) operations in 𝕂 [6,

ABSTRACT MODULAR COMPOSITION IN THE SEPARABLE CASE

For any field 𝕂 and n ∈ ℕ, we denote

𝕂[x] <n ≔ {P ∈ 𝕂[x]: deg P < n}.
In this section, 𝕂 represents an abstract algebraically closed field of constants. Let h = x n + h n-1 x n-1 + ⋯ + h 0 ∈ 𝕂[x] be a separable monic polynomial, so h admits n pairwise distinct roots 𝜎 1 , …, 𝜎 n in 𝕂. Then we may use the following algorithm for composition modulo h:

Algorithm 1 Input. Polynomials f , g ∈ 𝕂[x] <n and pairwise distinct 𝜎 1 , …, 𝜎 n ∈ 𝕂. Output. f ∘ g rem h, where h = (x -𝜎 1 ) ⋯ (x -𝜎 n ). 1. Compute v 1 = f (𝜎 1 ), …, v n = f (𝜎 n ) using fast multi-point evaluation. 2. Compute w 1 = g(v 1 ), …, w n = g(v n ) using fast multi-point evaluation. 3. Retrieve 𝜚 ∈ 𝕂[x] <n with 𝜚(𝜎 1 ) = v 1 , …, 𝜚(𝜎 n ) = v n using fast interpolation. 4. Return 𝜚. THEOREM 1. Algorithm 1 is correct and requires O(M(n) log n) operations in 𝕂. Proof. By construction, 𝜚(𝜎 i ) = (g ∘ f )(𝜎 i ) = (g ∘ f rem h)(𝜎 i ) for i = 1, …, n. Since deg 𝜚 < n
and the 𝜎 i are pairwise distinct, it follows that 𝜚 = g ∘ f rem h. This proves the correctness of the algorithm. The complexity bound follows from the fact that steps 1, 2 and 3 take

O(M(n) log n) operations in 𝕂. □
We wish to apply the theorem in the case when 𝕂 = ℂ. Of course, on a Turing machine, we can only approximate complex numbers with arbitrarily high precision, and likewise for the field operations in ℂ. For given numbers x and y, approximations at precision p for x + y, xy, x × y and x/y (whenever y ≠ 0) can all be computed in time O(I(p)). In view of Theorem 1, it is therefore natural to ask whether p-bit approximations of the coefficients of g ∘ f rem h may be computed in time

O(I(p) M(n) log n).
In the remainder of this paper we give a positive answer to a carefully formulated version of this question. Our first task is to make the concept of "approximations at precision p" more precise and to understand the way errors accumulate when performing a sequence of computations at precision p. We rely on "fixed point ball arithmetic" for this matter, as described in the next subsection. At a second stage, we prove a complexity bound for modular composition that holds for a fixed modulus h with known roots 𝜎 1 , …, 𝜎 n and for sufficiently large working precisions p.

The assumption that the roots 𝜎 1 , …, 𝜎 n of h are known is actually quite harmless in this context for the following reason: as soon as approximations for 𝜎 1 , …, 𝜎 n are known at a sufficiently high precision, the computation of even better approximations can be done fast using Newton's method combined with multi-point evaluation. Since we are only interested in the complexity for "sufficiently large working precisions", the computation of the initial approximations of 𝜎 1 , …,𝜎 n can therefore be regarded as a precomputation of negligible cost.

BALL ARITHMETIC AND STRAIGHT-LINE PROGRAMS

Fixed point numbers

Let a be a real number, we write ⌊a⌋ for the largest integer less or equal to a and ⌊a⌉ ≔ ⌊a + / 1 2 ⌋ for the closest integer to a.

Given a precision p ∈ ℕ, we denote by 𝔻 p = ℤ 2 -p the set of fixed point numbers with p binary digits after the dot. This set 𝔻 p is clearly stable under addition and subtraction. We can also define approximate multiplication × p on 𝔻 p using x × p y = ⌊2 p x y⌉ 2 -p , so |x × p yx y| ⩽ 2 -p-1 for all x, y ∈ 𝔻 p .

For any fixed constant K > 0 and x, y

∈ 𝔻 p ∩ [-K, K],

Fixed point ball arithmetic

Ball arithmetic is used for providing reliable error bounds for approximate computations. A ball is a set ℬ(c, r) = {z ∈ ℝ: |z -c| ⩽ r} with c ∈ ℝ and r ∈ ℝ ⩾ . From the computational point of view, we represent such balls by their centers c and radii r. We denote by 𝔹 p the set of balls with centers in 𝔻 p and radii in 𝔻 p ⩾ . Given vectors x = (x 1 , …, x n ) ∈ ℝ n and 𝒙 = (𝒙 1 , …, 𝒙 n ) = (ℬ(c 1 , r 1 ), …, ℬ(c n , r n )) ∈ 𝔹 p n we write x ∈ 𝒙 to mean x 1 ∈ 𝒙 1 ∧ ⋯ ∧ x n ∈ 𝒙 n , and we also set rad(𝒙) ≔ max (r 1 , …, r n ).

Let D be an open subset of ℝ n . We say that 𝑫 p ⊆ 𝔹 p n is a domain lift at precision p if 𝒙 ⊆ D for all 𝒙 ∈ 𝑫 p . The maximal such lift is given by 𝑫 p = {𝒙 ∈ 𝔹 p n : 𝒙 ⊆ D}. Given a function f : D → ℝ m , a ball lift of f at precision p is a function 𝒇 p : 𝑫 p → 𝔹 p m , where 𝑫 p = dom 𝒇 p is a domain lift of D at precision p, that satisfies the inclusion property: for any 𝒙 = (𝒙 1 , …, 𝒙 n ) ∈ 𝑫 p n and x = (x 1 , …, x n ) ∈ ℝ n , we have

x ∈ 𝒙 ⟹ f (x) ∈ 𝒇 p (𝒙).
A ball lift 𝒇 of f is a computable sequence (𝒇 p ) p∈ℕ of ball lifts at every precision such that for any sequence (𝒙 p ) p∈ℕ with 𝒙 p ∈ dom 𝒇 p , we have

lim p→∞ rad(𝒙 p ) = 0 ∧ ∩ p∈ℕ 𝒙 p ≠ ∅ ⟹ lim p→∞ rad(𝒇 p (𝒙 p )) = 0.
This condition implies the following:

lim p→∞ rad(𝒙 p ) = 0 ∧ ∩ p∈ℕ 𝒙 p = {x} ⟹ ∩ p∈ℕ 𝒇 p (𝒙 p ) = { f (x)}.
We say that 𝒇 is maximal if dom 𝒇 p is the maximal domain lift for each p. Notice that a function f must be continuous in order to admit a maximal ball lift.

The following formulas define maximal ball lifts ⊕ p , ⊖ p and ⊗ p at precision p for the ring operations +, -and ×:

ℬ(a, r) ⊕ p ℬ(b, s) ≔ ℬ(a + b, r + s) ℬ(a, r) ⊖ p ℬ(b, s) ≔ ℬ(a -b, r + s) ℬ(a, r) ⊗ p ℬ(b, s) ≔ ℬ(a × p b, (|a| + r) × p s + |b| × p s + 2 1-p ).
The extra 2 1-p in the formula for multiplication is needed in order to counter the effect of rounding errors that might occur in the three multiplications a× p b, (|a| + r)× p s and |b| × p s. For ℬ(a, r) ∈ 𝔹 p with r < |a|, the following formula also defines a maximal ball lift 𝜾 p at precision p for the inversion:

𝜾 p (ℬ(a, r)) ≔ ℬ(𝜄 p (a), 𝜄 p (|a| -r) -𝜄 p (|a|) + 2 1-p ).
For any fixed constant K > 0 and a, r, b, s ∈ 𝔻 p ∩ [-K, K], we notice that ℬ(a, r) ⊕ p ℬ(b, s), ℬ(a, r) ⊖ p ℬ(b, s), ℬ(a, r) ⊗ p ℬ(b, s) and 𝜾 p (ℬ(a, r)) can be computed in time O(I(p)).

Let 𝒇 be the ball lift of a function f : D → ℝ m with D ⊆ ℝ n . Consider a second ball lift 𝒈 of a function g: E → ℝ l with f (D) ⊆ E ⊆ ℝ m . Then we may define a ball lift 𝒈 ∘ 𝒇 of the composition g ∘ f : D → ℝ l as follows. For each precision p, we take (𝒈 ∘ 𝒇 ) p = 𝒈 p ∘ (𝒇 p ) |D p , where (𝒇 p ) |D P is the restriction of 𝒇 p to the set D p = {𝒙 ∈ dom 𝒇 p : 𝒇 p (𝒙) ∈ dom 𝒈 p }.

We shall use ball arithmetic for the computation of complex functions ℂ n → ℂ m simply through the consideration of real and imaginary parts. This point of view is sufficient for the asymptotic complexity point of view of the present paper. Of course, it would be more efficient to directly compute with complex balls (i.e. balls with a complex center and a real radius), but this would involve approximate square roots and ensuing technicalities.

The Lipschitz property

Assume that we are given the ball lift 𝒇 of a function f :D→ℝ m with D⊆ℝ n . Given a subset U ⊆D and constants 𝜆 ⩾ 0, 𝜇 ⩾ 0, we say that the ball lift 𝒇 is (𝜆, 𝜇)-Lipschitz on U if

∃p 0 ∈ ℕ, ∃𝜚 > 0, ∀p ⩾ p 0 , ∀𝒙 ∈ 𝔹 p n , 𝒙 ⊆ U ∧ rad(𝒙) ⩽ 𝜚 ⟹ 𝒙 ∈ dom 𝒇 p ∧ rad(𝒇 p (𝒙)) ⩽ 𝜆 rad(𝒙) + 𝜇 2 -p .
For instance, the ball lifts ⊕ and ⊖ of addition and subtraction are (2, 0)-Lipschitz on ℝ 2 . Similarly, the ball lift ⊗ of multiplication is (3 𝜆, 3)-Lipschitz on U = {(x, y) ∈ ℝ 2 : |x| ⩽ 𝜆, |y| ⩽ 𝜆} (by taking 𝜌 = 𝜆), whereas the ball lift 𝜾 of 𝜄 is (𝜆, 3)-Lipschitz on U = {x ∈ ℝ: 𝜆 -1/2 ⩽ |x|}.

Given 𝒇 and 𝜆 > 0, 𝜇 ⩾ 0 as above, we say that 𝒇 is locally (𝜆, 𝜇)-Lipschitz on U if 𝒇 is (𝜆, 𝜇)-Lipschitz on each compact subset of U. We define 𝒇 to be 𝜆-Lipschitz (resp. locally 𝜆-Lipschitz) on U if there exists a constant 𝜇 > 0 for which 𝒇 is (𝜆, 𝜇)-Lipschitz (resp. locally (𝜆, 𝜇)-Lipschitz). If 𝒇 is locally 𝜆-Lipschitz on U, then it is not hard to see that f is necessarily locally Lipschitz on U, with Lipschitz constant 𝜆. That is,

∀x ∈ U, ∃𝜂 > 0, ∀a, b ∈ ℬ(x, 𝜂) ∩ U, ‖ f (b) -f (a)‖ ∞ ⩽ 𝜆 ‖b -a‖ ∞ .
In fact, the requirement that a computable ball lift 𝒇 is 𝜆-Lipschitz implies that we have a means to compute high quality error bounds. We finally define 𝒇 to be Lipschitz (resp. locally Lipschitz) on U if there exists a constant 𝜆 > 0 for which 𝒇 is 𝜆-Lipschitz (resp. locally 𝜆-Lipschitz). Proof. Consider a compact subset C ⊆ U. Since this implies f (C) to be a compact subset of f (U) ⊆ V, it follows that there exists an 𝜀 > 0 such that f (C) + ℬ(0, 𝜀) ⊆ V. Let p 0 ∈ ℕ, 0 < 𝜚 < (𝜀 -𝜇 2 -p )/𝜆 and 0 < 𝜚' be such that for any p ⩾ p 0 , 𝒙 ∈ 𝔹 p n and 𝒚 ∈ 𝔹 p m , we have

𝒙 ⊆ C ∧ rad(𝒙) ⩽ 𝜚 ⟹ 𝒙 ∈ dom 𝒇 p ∧ rad(𝒇 p (𝒙)) ⩽ 𝜆 rad(𝒙) + 𝜇 2 -p < 𝜀 (𝒚 ⊆ f (C) + ℬ(0, 𝜀)) ∧ rad(𝒚) ⩽ 𝜚' ⟹ 𝒚 ∈ dom 𝒈 p ∧ rad(𝒈 p (𝒚)) ⩽ 𝜆' rad(𝒚) + 𝜇' 2 -p .
Given x ∈ 𝔹 p n with 𝒙 ⊆ C and rad(𝒙) ⩽ 𝜚, it follows that 𝒚 ≔ 𝒇 p (𝒙) satisfies rad(𝒚) < 𝜀, whence 𝒚 ⊆ f (C) + ℬ(0, 𝜀). If we also assume that rad(𝒙) ⩽ (𝜚' -𝜇 2 -p ) / 𝜆, then it also follows that rad(𝒚) ⩽ 𝜚', whence 𝒚 ∈ dom 𝒈 p and rad(𝒈 p (𝒚)) ⩽ 𝜆' (𝜆 rad(𝒙) + 𝜇 2 -p ) + 𝜇' 2 -p = 𝜆 𝜆' rad(𝒙) + (𝜇 𝜆' + 𝜇') 2 -p . In other words, if 𝒙⊆C and rad(𝒙)⩽min(𝜚,(𝜚' -𝜇 2 -p )/ 𝜆), then 𝒙∈dom(𝒈 p ∘𝒇 p ) and rad((𝒈 p ∘ 𝒇 p )(𝒙)) ⩽ 𝜆 𝜆' rad(𝒙) + (𝜇 𝜆' + 𝜇') 2 -p . □

STRAIGHT-LINE PROGRAMS

A signature is a finite or countable set of function symbols ℱ together with an arity r f ∈ ℕ for each f ∈ ℱ. A model for ℱ is a set K together with a function f K :

U f → K with U f ⊆ K r f for each k ∈ ℱ.
If K is a topological space, then U f is required to be an open subset of K r f . Let 𝒱 be a countable and ordered set of variable symbols. A straight-line program Γ with signature ℱ is a sequence Γ 1 ,…,Γ ℓ of instructions of the form

Γ k ≡ X k ≔ f k (Y k,1 , …, Y k,r f k ), where f k ∈ ℱ and X k , Y k,1 , …, Y k,r f k ∈ 𝒱, together with a subset 𝒪 Γ ⊆ {X 1 , …, X ℓ } of output variables.
Variables that appear for the first time in the sequence in the right-hand side of an instruction are called input variables. We denote by ℐ Γ the set of input variables. The number ℓ is called the length of Γ.

There exist unique sequences

I 1 < ⋯ < I n and O 1 < ⋯ < O m with ℐ Γ = {I 1 , …, I n } and 𝒪 Γ = {O 1 , …, O m }. Given a model K of ℱ we can run Γ for inputs in K, provided that the arguments Y k,1 , …, Y k,r f k are always in the domain of f k when executing the instruction Γ k . Let D Γ,K be the set of tuples I = (I 1 , …, I n ) ∈ K n on which Γ can be run. Given I ∈ K n , let Γ K (I) ∈ K m denote the value of (O 1 , …, O m ) at the end of the program. Hence Γ gives rise to a function Γ K : D Γ,K → K m .
Now assume that (ℝ, ( f ℝ ) f ∈ℱ ) is a model for ℱ and that we are given a ball lift 𝒇 of f ℝ for each f ∈ ℱ. Then 𝔹 p is also a model for ℱ at each precision p, by taking f 𝔹 p = 𝒇 p for each f ∈ ℱ. Consequently, any SLP Γ as above gives rise to both a function Γ ℝ : D Γ,ℝ → ℝ m and a ball lift Γ 𝔹 p : D Γ,𝔹 p → 𝔹 p m at each precision p. The sequence (Γ 𝔹 p ) p thus provides us with a ball lift 𝚪 for Γ ℝ .

PROPOSITION 3. If the ball lift 𝒇 of f ℝ is Lipschitz for each f ∈ ℱ, then 𝚪 is again Lipschitz.

Proof. For each model K of ℱ, for each variable v ∈ 𝒱 and each input I = (I 1 , …, I n ) ∈ D Γ,K , let v K,k (I) denote the value of v after step k. We may regard v K,k as a function from D Γ,K to K.

In particular, we obtain a computable sequence of functions v 𝔹 p ,k that give rise to a ball lift 𝒗 (k) of v ℝ,k . Let us show by induction over k that 𝒗 (k) is Lipschitz for every v ∈ 𝒱. This is clear for k = 0, so let k > 0. If v ≠ X k , then we have 𝒗 (k) = 𝒗 (k-1) ; otherwise, we have

𝒗 (k) = 𝒇 k( 𝒀 k,1 (k-1) , …, 𝒀 k,r f k (k-1) ) .
In both cases, it follows from Lemma 2 that 𝒗 (k) is again a Lipschitz ball lift. We conclude by noticing that 𝚪 = (𝑶 1 (ℓ) , …, 𝑶 n (ℓ) ). □

COMPUTABLE NUMBERS AND ULTIMATE COMPLEXITY

A real number x ∈ ℝ is said to be computable if there exists an approximation algorithm x ˇthat takes p ∈ ℕ on input and produces x ˇ(p) ∈ 𝔻 p on output with |xx ˇ(p)| ⩽ 2 -p (we say that x ˇ(p) is a 2 -p -approximation of x). We denote by ℝ com the field of computable real numbers. Let T(p) be a nondecreasing function. We say that a computable real number x ∈ ℝ com has ultimate complexity T(p) if it admits an approximation algorithm x ˇthat computes x ˇ(p) in time T(p + 𝛿) for some fixed constant 𝛿 ∈ ℕ. The fact that we allow x ˇ(p) to be computed in time T(p + 𝛿) and not T(p) is justified by the observation that the position of the "binary dot" is somewhat arbitrary in the approximation process of a computable number.

The notion of approximation algorithm generalizes to vectors with real coefficients: given v ∈ (ℝ com ) n , an approximation algorithm for v as a whole is an algorithm v ˇthat takes p ∈ ℕ on input and returns v ˇ(p) ∈ 𝔻 p r on output with |v ˇ(p) iv i | ⩽ 2 -p for i = 1, …, n. This definition naturally extends to any other mathematical objects that can be encoded by vectors of real numbers: complex numbers (by their real and complex parts), polynomials and matrices (by their vectors of coefficients), etc. The notion of ultimate complexity also extends to any of these objects.

A ball lift 𝒇 is said to be computable if there exists an algorithm for computing 𝒇 p for all p∈ ℕ. A computable ball lift 𝒇 of a function f : D → ℝ m with D ⊆ ℝ n allows us to compute the restriction of f to D ∩ (ℝ com ) n : given x ∈ D ∩ (ℝ com ) n with approximation algorithm x ˇ, by taking 𝒙 p = ℬ(x ˇ(p), 2 -p ) ∈ 𝔹 p n , we have ∩p∈ℕ 𝒙 p = {x}, ∩p∈ℕ 𝒇 p (𝒙 p ) = { f (x)}, and lim p→∞ rad(𝒇 p (𝒙 p )) = 0. Proof. Let x ˇbe an approximation algorithm for x of complexity T(p + 𝛿), where 𝛿 ∈ ℕ. There exist p 0 ∈ ℕ and a compact ball C around x with C ⊆ dom f and such that 𝒙 p = ℬ(x ˇ(p), 2 -p ) ∈ 𝔹 p n is included in C for all p ⩾ p 0 . There also exists a constant 𝛿' ∈ ℕ such that 𝒇 p (𝒙 p ) can be computed in time F(p + 𝛿') for all p ⩾ p 0 . Since 𝒇 is locally Lipschitz, there exists yet another constant 𝛿'' ∈ ℕ such that rad(𝒇 p (𝒙 p )) ⩽ 2 𝛿''-p for p ⩾ p 0 . For q = p -𝛿'' ⩾ max (p 0 -𝛿'', 0) and 𝛿''' = max (𝛿, 𝛿'), this shows that we may compute a 2 -q -approximation of f (x) in time T(q + 𝛿''') + F(q + 𝛿'''

Let

(p). Let Γ = Γ 1 , …, Γ ℓ be an SLP over ℱ whose k-th instruction Γ k writes X k ≔ f k (Y k,1 , …, Y k,r f ). Then, the ball lift 𝚪 of Γ ℝ has ultimate complexity F Γ (p) ≔ F f 1 (p) + ⋯ + F f ℓ (p).
Proof. This is a direct consequence of Proposition 5. □ COROLLARY 7. Let Γ be an SLP of length ℓ over ℱ = {0, 1,+, -,×, 𝜄} (where 0 and 1 are naturally seen as constant fonctions of arity zero). Then, there exists a ball lift 𝚪 of Γ ℝ with ultimate complexity O(I(p) ℓ).

Proof. We use the ball lifts of section 4 for each f ∈ {+, -, ×, 𝜄}: they are locally Lipschitz and computable with ultimate complexity O(I(p)). We may thus apply the Theorem 6 to obtain F Γ (p) = O(I(p) ℓ). □

Let us now analyze the ultimate complexity C(n, p) of our algorithm. For large p ⩾ p 1 , the algorithm essentially performs two multi-point evaluations of ultimate cost 𝜅' I(p) M(n) log n for some constant 𝜅' that does not depend on p, and a recursive call. Consequently,

C(n, p) ⩽ 𝜅' I(p) M(n) log n + C(n, ⌈𝛽 p⌉).

We finally obtain an other constant 𝜅 ⩾ 𝜅' such that C(n, p) = 𝜅 I(p) M(n) log n, by summing up the geometric progression and using the fact that I(p)/ p is nondecreasing. The conclusion now follows from Lemma 4. □ Remark 10. A remarkable feature of the above proof is that the precision p 1 at which the Newton iteration can safely be used does not need to be known in advance. In particular, the proof does not require any a priori knowledge about the Lipschitz constants.

THEOREM 11. There exists a constant 𝜅 > 0 such that the following assertion holds. Let f , g ∈ ℂ com [x] <n and let h ∈ ℂ com [x] be separable and monic of degree n. Assume that ( f , g, h) has ultimate complexity T(n, p). Then 𝜚 = g ∘ f rem h has ultimate complexity T(n, p) + 𝜅 I(p) M(n) log n.

Proof. This is an immediate consequence of the combination of the two above lemmas. □

CONCLUSION AND FINAL REMARKS

With some more work, we expect that all above bounds of the form O(I(p) M(n) log n) can be lowered to O(I(n p) log n). Notice that I(n p) log n = O(I(p) n log n) for p ⩾ n, when taking I(p) = Θ(n log n 8 log * n ) [START_REF] Harvey | Even faster integer multiplication[END_REF]. In order to prove this stronger bound using our framework, one might add an auxiliary operation × [n] for the product of two polynomials of degrees <n to the set of signatures ℱ. Polynomial products of this kind can be implemented for coefficients in 𝔻 p [i] with p ⩾ n using Kronecker substitution. For bounded coefficients, this technique allows for the computation of one such product in time O(I(n p)). By using Theorem 6, a standard complexity analysis should show that multi-point evaluation and interpolation have ultimate complexity O(I(n p) log n). By Theorem 11, the actual bit complexity of modular composition is of the form T(n, p + 𝛿)+ 𝜅 I(p + 𝛿) M(n) log n for some value of 𝛿 that depends on f , g, h (hence of n). An interesting problem is to get a better grip on this value 𝛿, which mainly depends on the geometric proximity of the roots of h.

If f , g, h belong to ℚ[x], then T(n, p) = O(n I(p)) and we may wish to bound 𝛿 as a function of n and the maximum bit size l of the coefficients of f , g and h. This would involve bit complexity results for root isolation [START_REF] Neff | An efficient algorithm for the complex roots problem[END_REF][START_REF] Pan | Univariate polynomials: nearly optimal algorithms for numerical factorization and root-finding[END_REF][START_REF] Schönhage | The fundamental theorem of algebra in terms of computational complexity[END_REF], for multi-point evaluation, and for interpolation. The overal complexity should then be compared with the maximal size of the output, namely g ∘ f rem h, which is in general much larger than the input size.

If h is not separable, but if a separable decomposition is known, then the techniques developed in this paper could be combined with Ritzmann's algorithm for the composition of formal power series [START_REF] Ritzmann | A fast numerical algorithm for the composition of power series with complex coefficients[END_REF]. If such a separable decomposition is not known, then it is an interesting problem to obtain a general algorithm for modular composition with a similar complexity (but this seems far beyond the scope of this paper).

BIBLIOGRAPHY

[1] D. Bernstein. Scaled remainder trees. Available from https://cr.yp.to/arith/scaledmod-20040820.pdf, 2004.

LEMMA 2 .

 2 Let 𝒇 be a locally (𝜆, 𝜇)-Lipschitz ball lift of f : D → ℝ m on an open set U. Let 𝒈 be a locally (𝜆', 𝜇')-Lipschitz ball lift of g: E → ℝ l on an open set V. If f (D) ⊆ E and f (U) ⊆ V, then 𝒈 ∘ 𝒇 is a locally (𝜆 𝜆', 𝜇 𝜆' + 𝜇')-Lipschitz ball lift of g ∘ f on U.

  we notice that x + y and xy can be computed in time O(p), whereas x × p y can be computed in time I(p) + O(p). Similarly, one may define an approximate inversion 𝜄 p on 𝔻 p

≠ ≔ 𝔻 p ∖ {0} by 𝜄 p (x) = ⌊2 p x -1 ⌉ 2 -p . For any fixed constant K > 0 and x ∈ 𝔻 p ≠ ∩ [-K, K],

we may compute 𝜄 p (x) in time O(I(p)).

  F be a nondecreasing function and assume that D is open. We say that 𝒇 has ultimate complexity F(p) if for every compact set C ⊆D, there exist constants p 0 ∈ℕ, 𝜚 >0 and 𝛿 ∈ℕ such that for any p ⩾ p 0 and 𝒙 p ∈ dom 𝒇 p with 𝒙 p ⊆ C and rad(𝒙 p ) ⩽ 𝜚, we can compute 𝒇 p (𝒙 p ) in time F(p + 𝛿). For instance, ⊕ and ⊖ have ultimate complexity O(p), whereas ⊗ and 𝜾 have ultimate complexity O(I(p)). PROPOSITION 4. Assume that 𝒇 is locally Lipschitz. If 𝒇 has ultimate complexity F(p) and x ∈ D ∩ (ℝ com ) n has ultimate complexity T(p), then f (x) has ultimate complexity T(p) + F(p).

ULTIMATE MODULAR COMPOSITION FOR SEPARABLE MODULI

LEMMA 8. There exists a constant 𝜅 > 0 such that the following assertion holds. Let f , g ∈ ℂ com [x] <n , let 𝜎 1 , …, 𝜎 n be pairwise distinct elements of ℂ com , and let h = (x -𝜎 1 ) ⋯ (x -𝜎 n ). Assume that ( f 0 , …, f n-1 , g 0 , …, g n-1 , 𝜎 1 , …, 𝜎 n ) has ultimate complexity T(n, p). Then 𝜚 = g ∘ f rem h has ultimate complexity T(n, p) + 𝜅 I(p) M(n) log n.

Proof. The algorithm for fast multi-point evaluation of a polynomial P = ∑ i=0 n-1 P i x i ∈ 𝕂[x] <n at 𝜉 1 , …, 𝜉 n ∈ 𝕂 can be regarded as an SLP over ℱ = {0, 1, +, -, ×, 𝜄} of length O(M(n) log n) that takes (P 0 , …, P n-1 , 𝜉 1 , …, 𝜉 n ) ∈ 𝕂 2n on input and that produces (P(𝜉 1 ), …, P(𝜉 n )) ∈ 𝕂 n on output. Similarly, the algorithm for interpolation can be regarded as an SLP over ℱ of length 2n on input and that produces (P 0 , …, P n-1 ) ∈ 𝕂 n on output with v 1 = P(𝜉 1 ), …, v n = P(𝜉 n ). Altogether, we may regard the entire Algorithm 1 as an

The conclusion now follows from Proposition 4. □ According to the above lemma, we notice that the time complexity for computing 𝜚 = g ∘ f rem h is T(n, p + 𝛿) for some constant 𝛿 that depends on n, f , g, and the 𝜎 i . LEMMA 9. There exists a constant 𝜅 > 0 such that the following assertion holds. Let h ∈ ℂ com [x] be separable and monic of degree n, and denote the roots of h by 𝜎= (𝜎 1 , …,𝜎 n ). If h has ultimate complexity T(n, p), then 𝜎 has ultimate complexity T(n, p)

Proof. There are many algorithms for the certified computation of the roots of a separable complex polynomial. We may use any of these algorithms as a "fall back" algorithm in the case that we only need a 2 -p -approximation of 𝜎 at a low precision p determined by h only.

For general precisions p, we use the following strategy in order to compute a ball 𝝈 ∈ 𝔹 p n with 𝜎 ∈ 𝝈 and rad(𝝈) ⩽ 2 -𝛼p for some fixed threshold / 1 2 < 𝛼 < 1. For some suitable p 0 ∈ ℕ and p ⩽ p 0 , we use the fall back algorithm. For p > p 0 and for a second fixed constant / 1 2 < 𝛽 < 1, we first compute a ball enclosure 𝝉 ∈ 𝔹 q n at the lower precision q = ⌈𝛽 p⌉ using a recursive application of the method. We next compute 𝝈 using a ball version of the Newton iteration, as explained below. If this yields a ball 𝝈 with acceptable radius rad(𝝈)⩽2 -𝛼p , then we are done. Otherwise, we resort to our fall-back method. Such calls of the fall-back method only occur if the default threshold precision p 0 was chosen too low. Nevertheless, we will show that there exists a threshold p 1 such that the computed 𝝈 by the Newton iteration always satisfies rad(𝝈) ⩽ 2 -𝛼p for p ⩾ p 1 .

Let us detail how we perform our ball version of the Newton iteration. Recall that 𝝉 ∈ 𝔹 q n with 𝜎 ∈ 𝝉 and rad(𝝉) ⩽ 2 -𝛼𝛽p is given. We also assume that we computed once and for all a 2 -p -approximation of h, in the form of a ball polynomial 𝒉 p ∈ 𝔹 p [i][x] of radius 2 -p that contains h. Now we evaluate 𝒉 p and 𝒉 p ' at each of the points 𝝈 1 , …, 𝝈 n using fast multi-point evaluation. Let us denote the results by 𝒗 = 𝒉 p (𝝈) and 𝒘 = 𝒉 p ' (𝝈). Let 𝜏, v and w denote the balls with radius zero and whose centers are the same as for 𝝉, 𝒗 and 𝒘. Using vector notation, the Newton iteration now becomes:

If 𝜎 ∈ 𝝉, then it is well-known [START_REF] Krawczyk | Newton-Algorithmen zur Bestimmung von Nullstellen mit Fehler-schranken[END_REF][START_REF] Rump | Kleine Fehlerschranken bei Matrixproblemen[END_REF] that 𝜎 ∈ 𝝈. Since rad(𝝉) ⩽ 2 -𝛼𝛽p , the fact that multipoint ball evaluation (used for 𝒉 p and 𝒉 p ' ) is locally Lipschitz implies the existence of a constant 𝛿 > 0 with rad(𝝂) ⩽ 2 𝛿-𝛼𝛽p and rad(𝒘) ⩽ 2 𝛿-𝛼𝛽p . Since h'(𝜎 i ) ≠ 0 for i = 1, …, n, there also exists a constant 𝛿' > 0 with 1 -𝜾 p (w) 𝒘 ⊆ ℬ(0, 2 𝛿'-𝛼𝛽p ). Altogether, this means that there exists a constant 𝛿'' > 0 with rad(𝝈) ⩽ 2 𝛿''-2𝛼𝛽p . Let p 1 = ⌈𝛿'' / (𝛼 (2 𝛽 -1))⌉. Then for any p ⩾ p 1 , the Newton iteration provides us with a 𝝈 with rad(𝝈) ⩽ 2 -𝛼p .