Cristina Silvano

Giovanni Agosta

Andrea Bartolini

Andrea R Beccari

Luca Benini

Loïc Besnard

João Bispo

Radim Cmar

João M P Cardoso

Carlo Cavazzoni

Daniele Cesarini

Stefano Cherubin

Federico Ficarelli

Davide Gadioli

Martin Golasowski

Antonio Libri

Jan Martinovič

Gianluca Palermo

Pedro Pinto

Erven Rohou

Kateřina Slaninová

Emanuele Vitali

The ANTAREX domain specific language for high performance computing

Keywords: High Performance Computing, Autotuning, Adaptivity, DSL, Compilers, Energy Efficiency

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

High Performance Computing (HPC) is a strategic asset vigorously pursued by both state actors, supra-national entities, and major industrial players in fields such as finance or oil & gas [START_REF] Curley | HPC and Big Data[END_REF]. A veritable arms race has been waged for decades to build the fastest HPC machine, first in terms of pure floating pooint operations per second (FLOPS), then, with the growing difficulty of supplying power to large HPC centers, in terms of Flops/watt. Currently, the goal is to reach Exascale level (10 18 FLOPS) within the 2023 -24 timeframe -with a ×1000 improvement over Petascale, reached in 2009. The current Top5001 and Green500 2 lists are dominated respectively by IBM's Summit machine, installed at the Oak Ridge National Laboratory (USA), with over 200 PetaFlops, and Japan's Shoubu system B, with a power efficiency of 17.6 GFlops/watt 3 . To reach these levels of performance and power efficiency, heterogeneous computing architectures are critical, in particular through GPGPUs. Traditional NVIDIA GPG-PUs are featured heavily in both lists, though Shoubu's ZettaScaler architecture uses PEZY's PEZY-SC2 instead. The dominance of heterogeneous systems in the Green500 list is expected to continue for the next coming years to reach the target of 20MW Exascale supercomputers set by the DARPA.

When considering the design, development, and deplyoment of HPC applications, each increase in HPC machine size and heterogeneity imposes an increased burden on the developer. Nowadays, few developers have both the domain knowledge and HPC expertise to produce a fully optimised application, and a naïvely parallelised application may easily incur in parallelisation bottlenecks and scalability issues. To address this issue, the current development model assumes that the domain experts do not attempt to optimise their applications on their own, but rely on help from HPC experts, typically part of the HPC center staff, to turn their baseline application into an HPCenabled one. This task requires to handle parallelisation, offloading to heterogeneous resources, as well as performance analysis and tuning. For the latter task, it is worth noting that for large scale applications it is oftern difficult to make predictions on performance based on small-scale runs, so runtime autotuning is desirable.

While many approaches to programming languages and models for HPC attempt to simplify the development process by reducing or removing the need for the HPC expert, the recently completed ANTAREX project [START_REF] Silvano | The ANTAREX Tool Flow for Monitoring and Autotuning Energy Efficient HPC Systems[END_REF][START_REF] Silvano | Antarex: A dsl-based approach to adaptively optimizing and enforcing extra-functional properties in high performance computing[END_REF] fully embraces this split development process, and aims at supporting it. To this end, we introduce the ANTAREX DSL, a Domain Specific Language (DSL) to express the application self-adaptivity and to runtime manage and autotune applications for green heterogeneous HPC systems up to the Exascale level. The ANTAREX DSL allows the introduction of a separation of concerns, where self-adaptivity and energy efficient strategies, managed by the HPC expert, are specified separately from the application functionalities, developed instead by the domain expert. The ANTAREX DSL is based on Aspect Oriented Programming (AOP) principles, embodied in the LARA language [START_REF] Cardoso | LARA: An Aspect-oriented Programming Language for Embedded Systems[END_REF][START_REF] Cardoso | Performance-driven instrumentation and mapping strategies using the LARA aspect-oriented programming approach[END_REF], which effectively support the separation of concern without burdening the domain expert with the need to learn a new language (as they simply develop in C/C++). Since HPC experts already use a variety of scripts and tools to support their work, the ANTAREX DSL also provides a way to unify them under a single interface.

To demonstrate the effectiveness of the ANTAREX DSL and its associated tools, the project employs two use cases taken from highly relevant HPC application scenarios:

1. a biopharmaceutical application for drug discovery developed by one of the leading European companies in the field, Dompé, and deployed on the 1.21 PetaFlops heterogeneous NeXtScale Intel-based IBM system at CINECA; 2. a self-adaptive navigation system for smart cities developed by the top European navigation software company, Sygic, and deployed on the server-side on the 1.46 PetaFlops heterogeneous Intel® Xeon Phi™ based system provided by IT4Innovations National Supercomputing Center.

These scenario also reproduce the aforementioned development model, as each HPC user company is supported by the staff of the HPC center where the application is deployed.

The remaining partners of the ANTAREX Consortium comprise a wealth of expertise in DSL development, code optimisation, energy efficiency, and autotuning. Four top-ranked academic and research partners (Politecnico di Milano, ETHZ Zurich, University of Porto and INRIA) are complemented by the Italian Tier-0 Supercomputing Center (CINECA), the Tier-1 Czech National Supercomputing Center (IT4Innovations) and two industrial application providers, one of the leading biopharmaceutical companies in Europe (Dompé) and the top European navigation software company (Sygic). Politecnico di Milano, the largest Technical University in Italy, played the role of Project Coordinator.

The ANTAREX Approach. The ANTAREX approach and related tool flow, as shown in Figure 1, operate both at design-time and runtime. The application functionality is expressed through C/C++ code (possibly including legacy code), whereas the extrafunctional aspects of the application, including parallelisation, mapping, and adaptivity strategies, are expressed through DSL code (based on LARA) developed in the project. As a result, the expression of such aspects is fully decoupled from the functional code. The Clava tool is the centerpoint of the compile-time phase, performing a refactoring of the application code based on the LARA aspects, and instrumenting it with the necessary calls to other components of the tool flow.

The ANTAREX compilation flow leverages a runtime phase with compilation steps, through the use of partial dynamic compilation techniques enabled by libVC. The application autotuning, performed via the mARGOt tool, is delayed to the runtime phase, where the software knobs (application parameters, code transformations and code variants) are configured according to the runtime information coming from application self-monitoring as well as from system monitoring performed by the ExaMon tool. Finally, the runtime power manager, PowerCapper, is used to control the resource usage for the underlying computing infrastructure given the changing conditions. At runtime, the application control code, thanks to the design-time phase, now contains also runtime monitoring and adaptivity strategy code derived from the DSL extrafunctional specification. Thus, the application is continuously monitored to guarantee the required Service Level Agreement (SLA), while communication with the runtime resource-manager takes place to control the amount of processing resources needed by the application. The application monitoring and autotuning is supported by a runtime layer implementing an application level collect-analyse-decide-act loop.

Organization of the paper. The rest of this paper is organized as follows. In Section 2 we review the technology portfolio provided by the ANTAREX tool flow. In Section 3 we provide an assessment of the impact of the proposed DSL on application specifications, while Sections 4 and 5 show of how the Tool Flow has been applied the two use cases of the project. Finally, in Section 6 we compare the ANTAREX DSL with related works from the recent literature, and in Section 7 we draw some conclusions.

ANTAREX Technology Portfolio

In this Section, we first provide an overview of the ANTAREX DSL itself, and then we show how the DSL is used to apply each of the main components of the Tool Flow, providing examples of the code needed in each case to allow the reader to understand the level of complexity involved.

The ANTAREX DSL

HPC applications might profit from adapting to operational and situational conditions, such as changes in contextual information (e.g., workloads), in requirements (e.g., deadlines, energy), and in availability of resources (e.g., connectivity, number of processor nodes available). A simplistic approach to both adaptation specification and implementation (see, e.g., [START_REF] Floch | Using architecture models for runtime adaptability[END_REF]) employs hard coding of, e.g., conditional expressions and parameterizations. In our approach, the specification of runtime adaptability strategies relies on a DSL implementing key concepts from Aspect-Oriented Programming (AOP) [START_REF] Irwin | Aspect-oriented Programming[END_REF], mainly specifying adaptation concerns, targeting specific execution points, separately from the primary functionality of the application, with minimum or no changes to the application source code.

Our approach is based on the idea that certain application/system requirements (e.g., target-dependent optimizations, adaptivity behavior and concerns) should be specified separately from the source code that defines the main functionality. Those requirements are expressed as DSL aspects that embody strategies. An extra compilation step, performed by a weaver, merges the original source code and aspects into the intended program [START_REF] Elrad | Aspect-oriented Programming: Introduction[END_REF]. Using aspects to separate concerns from the core objective of the program can result in cleaner programs and increased productivity (e.g., reusability of strategies). As the development process of HPC applications typically involves two types of experts (application-domain experts and HPC system architects) that split their responsibilities along the boundary of functional description and extra-functional aspects, our DSL-aided toolflow provides a suitable approach for helping to express their concerns.

The ANTAREX DSL relies on the already existing DSL technology LARA [START_REF] Cardoso | LARA: An Aspect-oriented Programming Language for Embedded Systems[END_REF][START_REF] Cardoso | Performance-driven instrumentation and mapping strategies using the LARA aspect-oriented programming approach[END_REF]. In particular, the LARA technology provides a framework that we adopted to implement the ANTAREX aspects and APIs. Moreover, we developed other LARA-related tools such as the Clava4 weaver to leverage the rest of the ANTAREX tool flow.

LARA is a programming language that allows developers to capture non-functional requirements and concerns in the form of strategies, which are decoupled from the functional description of the application. Compared to other approaches that usually focus on code injection (e.g., [START_REF] Spinczyk | AspectC++: An Aspect-oriented Extension to the C++ Programming Language[END_REF]), LARA provides access to other types of actions, e.g., code refactoring, compiler optimizations, and inclusion of additional information, all of which can guide compilers to generate more efficient implementations. Additional types of actions may be defined in the language specification and associated weaver, such as software/hardware partitioning [START_REF] Cardoso | Controlling a complete hardware synthesis toolchain with LARA aspects[END_REF] or compiler optimization sequences [START_REF] Nobre | Use of Previously Acquired Positioning of Optimizations for Phase Ordering Exploration[END_REF]. One important feature of the LARA-aided source-to-source compiler developed in ANTAREX is the capability to refactor the code of the application in order to expose adaptivity behavior and/or adaptivity design points that can be explored by the ANTAREX autotuning component. In the following sections we show illustrative examples 5 of some of the strategies that can be specified using LARA in the context of a source-to-source compiler and currently used for one of the use cases.

Precision Tuning

Error-tolerating applications are increasingly common in the emerging field of realtime HPC, allowing to trade-off precision for performance and/or energy. Thus, recent works investigated the use of customized precision in HPC as a way to provide a breakthrough in power and performance. We developed a set of LARA aspects enabling mixed precision tuning on C/C++ and OpenCL kernels. In our precision tuning we combine an adaptive selection of floating and fixed point arithmetic, targeting HPC applications.

Figure 2 presents part of a LARA strategy that changes all declarations of a certain type to a target type (e.g., from double to float) for a given function. We note, however, that a practical and reusable aspect needs to deal with further issues, such as the cloning of functions whose types we want to change but are also called by other unrelated functions in the code, assignments of constants, casts, recursion, changing functions definitions and library functions to the ones related to the type used (e.g., sqrtf vs sqrt in Math.h), etc. In this example, changeType is a function that analyzes and changes compound types, such as double * and double[]. If the type described in $old is found inside the type of the declaration, it is replaced with the type described in $new. To be more specific, if $old is double, $new is float and $decl.type is double * , the type of the declaration will be changed to float * . If the original declaration type does not contain the $old type, it is not changed. A LARA aspect consists of three main steps. Firstly, one captures the points of interest in the code using a select statement, which in this example selects variable declarations. Then, using the apply statement, one acts over the selected program points. In this case, it will define the types of the captured declared variables, using the type attribute. Finally, we can then specify conditions to constrain the execution of the apply (i.e., only if the declared variable has a specific type). This can be done via conditional statements (ifs) as well as via special condition blocks that constrain the entire apply. LARA promotes modularity and aspect reuse, and supports embedding JavaScript code, to specify more sophisticated strategies. As shown in [START_REF] Nobre | Aspect-driven mixed-precision tuning targeting gpus[END_REF], we support exploration of mixed precision OpenCL kernels by using half, single, and double precision floating point data types. We additionally support fixed point representations through a custom C++ template-based implementation for HPC systems, which has already been used in [START_REF] Cherubin | Implications of Reduced-Precision Computations in HPC: Performance, Energy and Error[END_REF]. In both cases the LARA aspects automatically insert code for proper type conversion before and after the critical section that has been converted to exploit a reduced precision data type.

The LARA aspect in Figure 3 shows the generation of different mixed-precision versions to be dynamically evaluated. It is possible to specify -as input of the aspect -the number of mix combination to generate, and a rule set to filter out precision mix combinations which are very likely to lead to useless and/or not efficient results. We exploit programmer's application domain knowledge by relying on them to define test cases to evaluate the different code versions at runtime. LARA automatically inserts code to dynamically perform the exploration over the space of the generated versions with different precision mix.

Code Versioning

One of the strategies supported in the ANTAREX toolflow is the capability to generate versions of a function and to select the one that satisfies certain requirements at runtime. Figure 4 shows an aspect that clones a set of functions and changes the types of the newly generated clones. Each clone has the same name as the original with the addition of a provided suffix. We start with a single user-defined function which is cloned by the aspect CloneFunction (called in line 13). Then, it recursively traverses calls to other functions inside the clone and generates a clone for each of them. Inside the clones, calls to the original functions are changed to calls to the clones instead, building a new call tree with the generated clones. At the end of the aspect CreateFloatVersion (lines 16-17,) we use the previously defined ChangePrecision aspect to change the types of all generate clones. The aspect Multiversion -in Figure 5 -adapts the source code of the application in order to call the original version of a function or a generated cloned version with a different type, according to the value of a parameter given by the autotuner at runtime. The main aspect calls the previously shown aspect, CreateFloatVersion, which clones the target function and every other function it uses, while also changing their variable types from double to float (using the aspects presented in Figure 4 and Figure 2). This is performed in lines 8-9 of the example. From lines 13 to 34, the Multiversion aspect generates and inserts code in the application that is used as switching mechanism between the two versions. It starts by declaring a variable to be used as a knob by the autotuner, then it generates the code for a switch statement and replaces the statement containing the original call with the generated switch code. Finally, in lines 36-38, the aspect surrounds both calls (original and float version) with timing code. An excerpt of the resulting C code can bee seen in Figure 6.

In the ANTAREX toolflow, the capability of providing several versions of the same function is not limited to static features. LIBVERSIONINGCOMPILER [START_REF] Cherubin | libVersioningCompiler: An easy-to-use library for dynamic generation and invocation of multiple code versions[END_REF][START_REF] Festa | Continuous Program Optimization via Advanced Dynamic Compilation Techniques[END_REF] (abbreviated LIBVC) is an open-source C++ library designed to support the dynamic generation and versioning of multiple versions of the same compute kernel in a HPC scenario. It can be used to support continuous optimization, code specialization based on the input data or on workload changes, or to dynamically adjust the application, without the burden of a full just-in-time compiler. LIBVC allows a C/C++ compute kernel to be dynamically compiled multiple times while the program is running, so that different specialized versions of the code can be generated and invoked. Each specialized version can be versioned for later reuse. When the optimal parametrization of the compiler depends on the program workload, the ability to switch at runtime between different versions of the same code can provide significant benefits [START_REF] Chen | Evaluating iterative optimization across 1000 datasets[END_REF][START_REF] Tartara | Continuous learning of compiler heuristics[END_REF]. While such versions can be generated statically in the general case, in HPC execution times can be so long that exhaustive profiling may not be feasible. LIBVC instead enables the exploration and tuning of the parameter space of the compiler at runtime. Figure 7 shows an example of usage of LIBVC through LARA, which demonstrates how to specialize a function. The user provides this aspect with a target function call and a set of compilation options. These include compiler flags and possible compiler definitions, e.g., data discovered at runtime, which is used as a compile-time constant in the new version. Based on the target function call, the aspect finds the function It is worth noting that the combination of LARA and LIBVC can also be used to support compiler flag selection and phase-ordering both statically and dynamically [START_REF] Ashouri | Cobayn: Compiler autotuning framework using bayesian networks[END_REF][START_REF] Ashouri | Micomp: Mitigating the compiler phase-ordering problem using optimization sub-sequences and machine learning[END_REF].

Memoization

Memoization is a technique that trades off computation time for memory space by storing computed values and reusing them instead of recomputing them. It has been shown to be effective in improving performance and reducing energy consumption in large scale applications [START_REF] Agosta | Dynamic memoization for energy efficiency in financial applications[END_REF]. We introduce in this section a memoization technique integrated in the ANTAREX toolflow. Performance can be improved by caching results of pure functions (i.e. deterministic functions without side effects), and retrieving them instead of recomputing a result. We have implemented the work of [START_REF] Suresh | Compile-Time Function Memoization[END_REF] generalized for C++ and aided with extensions regarding user/developer flexibility. We describe here only the principles of this technique and more details can be found in [START_REF] Suresh | Intercepting Functions for Memoization: A Case Study Using Transcendental Functions[END_REF] [START_REF] Suresh | Compile-Time Function Memoization[END_REF]. 2. The substitution of the references to foo by foo_wrapper in the application. This technique has been extended for C++ memoizable methods and takes into account the mangling, the overloading, and the references to the objects. Memoization is proposed in the ANTAREX project by relying on aspects programmed using the DSL. The advantage of these aspects is that the memoization is integrated into the application without requiring user modifications of the source code. The code generated by Clava is then compiled and linked with the associated generated memoization library.

An example of a LARA aspect for memoization is shown in Figure 9. It defines the memoization (lines 1-13) of a method (aMethod) of a class (aClass) with nbArg parameters of same type as the returned type (Type). Note that the inputs nbArg and Type are required to manage the overloading of the object-oriented languages such that C++. Other parameters (from line 4) are provided to improve several memoization approaches. For examples, the user can specify [START_REF] Curley | HPC and Big Data[END_REF] (lines 14-15), the method is searched (lines 17-24). Then, in case of success, the code of the wrapper is added (line 28) to produce the memoization library, and (line 30) the code of the application is modified for calling the created "wrapper", this wrapper is also declared as a new method of the class.

Moreover, some variables are exposed for autotuning in the memoization library. For each function or method, a variable that manages the dynamical "stop/run" of the memoization is exposed, as well as the variable that manages the policy to use in case of conflict to the table. To be complete about the memoization, a LARA aspect is proposed to automatically detect the memoizable functions or methods. Then the user may decide to apply or not the memoization on these selected elements.

Self-Adaptivity & Autotuning

In ANTAREX, we consider each application's function as a parametric function that elaborates input data to produce an output (i.e., o = f (i, k 1 , . . . , k n ) ), with associated extra-functional requirements. In this context, the parameters of the function (k 1 , . . . , k n ) are software-knobs that modify the behavior of the application (e.g., parallelism level or the number of trials in a MonteCarlo simulation). The main goal of mARGOt6 [START_REF] Gadioli | Application autotuning to support runtime adaptivity in multicore architectures[END_REF] is to enhance an application with an adaptive layer, aiming at tuning the software knobs to satisfy the application requirements at runtime. To achieve this goal, the mARGOt dynamic autotuning framework developed in ANTAREX is based on the MAPE-K feedback loop [START_REF] Brun | Engineering Self-Adaptive Systems through Feedback Loops[END_REF]. In particular, it relies on an application knowledge, derived either at deploy time or at runtime, that states the expected behavior of the extra-functional properties of interest. To adapt, on one hand mARGOt uses runtime observations as feedback information for reacting to the evolution of the execution context. On the other hand, it considers features of the actual input to adapt in a more proactive fashion. Moreover, the framework is designed to be flexible, defining the application requirements as a multi-objective constrained optimization problem that might change at runtime.

To hide the complexity of the application enhancement, we use LARA aspects for configuring mARGOt and for instrumenting the code with related API. Figure 10 provides a simple example of a LARA aspect where mARGOt has been configured (lines 5-20) to actuate on a software knob Knob1 and target error and throughput metrics [START_REF] Gadioli | Socrates -a seamless online compiler and system runtime autotuning framework for energy-aware applications[END_REF]. In particular, the optimization problem has been defined as the maximization of the throughput while keeping the error under a certain threshold. The last part of the aspect (lines 23-27) is devoted to the actual code enhancement including the needed mARGOt call for initializing the framework and for updating the application configuration. The declarative nature of the LARA library developed for integrating mARGOt simplifies its usage hiding all the details of the framework.

Monitoring

Today processing elements embed the capability of monitoring their current performance efficiency by inspecting the utilization of the micro-architectural components as well as a set of physical parameters (i.e., power consumption, temperature, etc). These metrics are accessible through hardware performance counters which in x86 systems can be read by privilege users, thus creating practical problems for user-space libraries to access them. Moreover, in addition to sensors which can be read directly from the software running on the core itself, supercomputing machines embed sensors external to the computing elements but relevant to the overall energy-efficiency. These elements include the node and rack cooling components as well as environmental parameters such as the room and ambient temperature. In ANTAREX, we developed ExaMon [START_REF] Beneventi | Continuous learning of hpc infrastructure models using big data analytics and in-memory processing tools[END_REF] (Exascale Monitoring) to virtualise the performance and power monitoring access in a distributed environment. ExaMon decouples the sensor readings from the sensor value usage. Indeed, ExaMon uses a scalable approach were each sensor is associated to a sensing agent which periodically collects the metrics and sends the measured value with a synchronized time-stamp to an external data broker. The data broker organises the incoming data in communication channels with an associated topic. Every new message on a specific topic is then broadcast to the related subscribers, according to a list kept by the broker. The subscriber registers a callback function to the given topic which is called every time a new message is received. To let LARA take advantage of this monitoring mechanism we have designed the Collector API, which allow the initialization of the Collector component associated with a specific topic that keeps an internal state of the remote sensor updated. This internal state can then be queried asynchronously by the Collector API to gather its value. LARA aspects have been designed to embed the Collector API and to make the application code self-aware.

Figure 11 shows a usage example of ExaMon through LARA, which subscribes to a topic on a given broker and inserts a logging message in the application. To define the connection information, the user needs to provide the address to connect to, as well as the name of the topic to subscribe. As for the integration in the original application code the user needs to provide a target function, where the collector will be managed, and a target statement, where the query of the data and logging will be performed.

Power Capping

Today's computing elements and nodes are power limited. For this reason, stateof-the-art processing elements embed the capability of fine-tuning their performance to control dynamically their power consumption. This includes dynamic scaling of voltage and frequency, and power gating for the main architectural blocks of the processing elements, but also some feedback control logic to keep the total power consumption of the processing element within a safe power budget. This logic in x86 systems is named RAPL [START_REF] David | Rapl: Memory power estimation and capping[END_REF]. Demanding the power control of the processing element entirely to RAPL may not be the best choice. Indeed, it has been recently discovered that RAPL is application agnostic and thus tends to waist power in application phases which exhibit IOs or memory slacks. Under these circumstances there are operating points that proved to

ln().text("Time=").double(getTimeExpr(exa))

.text("[s], data=").double(exa.getMean()).ln(); // Add printf after last stmt logger.log($lastStmt); end be more energy efficient than the ones selected by RAPL while still respecting the same power budget [START_REF] Cesarini | Benefits in relaxing the power capping constraint[END_REF]. However, these are only viable if the power capping logic is aware of the application requirements. To do so, we have developed a new power capping run-time based on a set of user space APIs which can be used to define a relative priority for the given task currently in execution on a given core. Thanks to this priority, the run-time is capable of allocating more power to the higher priority process [START_REF] Bartolini | Self-aware thermal management for high performance computing processors[END_REF][START_REF] Cesarini | Prediction horizon vs. efficiency of optimal dynamic thermal control policies in hpc nodes[END_REF]. In ANTAREX, these APIs can be inserted by LARA aspects in the application code.

Evaluation

Tables 1 and2 show static and dynamic metrics collected for the weaving process of the presented examples. In Table 1, we can see the number of logical lines of source code for the LARA strategies, as well as for the input code and generated output code (the SLoC-L columns). In the last two columns we report the difference in SLoC and functions between the input and output code (the delta columns). Note the woven and delta results for the HalfPrecisionOpenCL strategy are the sum of all generated code, totaling 31 versions.

An inspection of columns LARA SLoC-L and Delta SLoC-L reveals that, in most examples, there is a large overhead in terms of LARA SLoC-L over application SLoC-L. While this may seem a problem, we need to consider that a large part of the work being performed by these strategies is code analysis, which does not translate directly to SLoC-L in the final application. Furthermore, the Delta SLoC-L metric does not account for removed application code and for these cases a metric based on the similarity degree among code versions could be of more interest. Also, in real-world applications, the ratio of LARA SLoC-L to application SLoC-L would be definitely more favorable, thanks to aspect reuse. To better understand the impact of analysis, we report in the first two columns of Table 2 the number of code points and of their attributes analysed, which can be compared with the last three column of the same table, which instead report the corresponding effects, in terms of the number of modified points and lines of code inserted. To understand the impact of removed lines of code, we look at the Inserts and Actions columns, which show that circa one half of the actions do not insert code. The end line is that the analysis work exceeds the transformation work by an order of magnitude, and the insertions only underestimate significantly the work performed.

Another benefit for user productivity when using LARA is how the techniques presented in the examples can scale into large-scale applications and scenarios. Most of the presented strategies are parameterized by function, i.e., they receive a function join point or name and act on the corresponding function. This could be performed manually, albeit crudely, using a search function of an IDE. Consider the case where we instead want to target a set of functions, whose names we may not know, based on their function signatures, or based on the characteristics of the variables declared inside their scope. This kind of search and filtering based on syntactic and semantic information available in the program is one of the key features of LARA and it cannot be easily attained with other tools. As the aspects presented here illustrate, LARA strategies can be made reusable and applied over large applications, greatly out scaling the effort needed to develop them.

Case Study: Computer Accelerated Drug Discovery

In this section, we provide an overview of the application of the ANTAREX tool flow to the Drug Discovery application developed in Use Case 1.

Computational discovery of new drugs is a compute-intensive task that is critical to explore the huge space of chemicals with potential applicability as pharmaceutical drugs. Typical problems include the prediction of properties of protein-ligand complexes (such as docking and affinity) and the verification of synthetic feasibility.

In the ANTAREX project it is expected that this application supports the exploration of millions of ligands over several pockets with sizes ranging from a few thousand atoms to several tens of thousands, in a timely manner. Performance is a key factor for this application, and we have used the ANTAREX tool flow to optimize the use case for Exascale systems.

Auto-parallelization and exploration through the ANTAREX DSL

Experiments on using an entire cluster for execution of the Use Case 1 have shown that relying only on MPI to parallelize the code both to distribute parallelism between nodes and to take advantage of parallelism inside a node is highly inadequate. We observed that a few tens of thousand tasks are enough to make MPI fail during execution. To overcome that, we opted to use MPI only for distributing work between nodes, and parallelize the application inside the node with OpenMP.

We have developed an auto-parallelization LARA library that uses static analysis to detect and insert OpenMP pragmas in for loops that can be parallelized. We have identified the function MeasureOverlap in the UC1 application as an important hotspot, and we have extracted it to a test application written in C (the autoparallelization library currently only supports C).

Figure 12 shows the LARA code that parallelizes the application. Line 6 calls the parallelization library, which parallelizes all the for loops in the code that where detected as safe to be parallelized. However, this can be a very inefficient strategy, since nesting parallelism with OpenMP can lead to excessive overhead.

The code in lines 9-28 detects and disables all nested OpenMP pragmas. It starts by iterating over all OpenMP pragmas in the code of the kind parallel for; then, it retrieves all descendants of the loop corresponding to that pragma, that are also OpenMP pragmas of the kind parallel for, and transforms those nested pragmas into inlined comments.

After parallelizing the code, we want to explore two parameters: the number of threads to use, and the size of the pocket to test. Figure 13 shows the LARA code for a strategy that uses the LARA library LAT 7 to explore these parameters. The strategy sets 50.000 atoms. We run the experiments on a computer with an Intel Xeon CPU E5-1650 v4 processor, at 3.60GHz and 64GB of RAM. The code was compiled with gcc 7.3.0.

Figure 14 shows the improvements in execution time (i.e., speedup) and energy (i.e., energy improvement) when the number of threads changes, when considering the biggest pocket size (15.000 atoms). We use as input the sizes of 120k ligands, based on a database provided by the industrial partner.

As expected for this case, the improvement curve of the execution time coincides with the improvement in energy consumed, and increases with the number of threads. The number of atoms in the pocket size modestly changes the speedup, i.e. between 0.7 % and 5.5 %, for the same number of threads. The amount of change increases with the number of threads, but for a given number of threads, an increase in pocket size improves speedup up to a point, and after that the speedup decreases (we consider that this happens due to saturation of memory resources). This kind of exploration can be used to generate data which can be fed to the autotuner. In this section, we provide an overview of the application of the ANTAREX tool flow to the Self-Adaptive Navigation System developed in Use Case 2. The system is designed to process large volumes of data for the global view computation and to handle dynamic loads represented by incoming routing requests from users of the system. Both disciplines require HPC infrastructure in order to operate efficiently while maintaining contracted SLA. Integration of the ANTAREX self-adaptive holistic approach can help the system to meet the mentioned requirements and pave the way to scaling its operation to future Exascale systems. Core of the system is a routing pipeline with several stages which uses our custom algorithm library written in C++. The library provides an API for the individual routing algorithms and for data access layer, which provides abstraction of a graph representation of the road network. The graph is stored in a HDF5 file, which is a well known and convenient storage format for structured data on HPC clusters.

Memoization and Precision Tuning through the ANTAREX DSL

The navigation system in our use case calculates routes simultaneously for possibly all participants in the system and it much relies on precise calculation of the macroscopic traffic flow model. The critical component for our traffic flow calculation is the After applying the strategy, the application has 775 lines of code, and uses a custom memoization library automatically generated by Clava with around 150 lines of C code. The number of generated lines depends on the number of memoized functions (one in this case). So, the total number of lines is 925, considering the memoization library. In order to optimize the memory allocation, we applied another modification to the Betweenness Centrality algorithm. It consists of moving statements outside the body of a loop without affecting the functionality of the program (also called code hoisting). This transformation was applied manually and is not defined in LARA. We also applied the memoization technique to the routing algorithms for one function, with a similar number of generated lines as the Betweenness Centrality algorithm.

Performance and scalability

For testing impact of INRIA tools on Betweenness Centrality algorithm, we chose a graph of the traffic network of Porto city in Portugal and one bigger graph of 4 countries (4EU -Czech Republic, Slovakia, Austria, Hungary). This experiment was focused on comparing the speed of an original and of a modified version of the Betweenness Centrality algorithm. The size of both graphs in terms of vertex and edge count is presented in Table 5.1.1.

The performance was tested on Salomon cluster operated by the IT4Innovations National Supercomputing Center. The cluster consists of 1,008 compute nodes and each of them contains 2x Intel Xeon E5-2680v3 processors clocked at 2.5 GHz and 128 GB DDR4@2133 RAM. When the experiments were performed, the operating system was CentOS 7.6 and the code was compiled using Intel C++ Compiler 17.0.

We tested several versions of the Betweenness Centrality algorithm with various combinations of the ANTAREX tools: Overall, we can see a speedup of 3.7 to 7.8% for the larger graph when adding the hoisting and memoization transformations to the double precision floating point code. Considering also the precision reduction, the performance speedup grows to 14.3 to 20.6% depending on the number of cores employed.

Autotuning for Client-side navigation

The Navigation application is the one which finally sells the whole navigation system, so it requires quite an attention. In a holistic approach the application implementation needs to consider multiple aspects, both in the functional and the performance space as was stipulated in the specification of the system. In simple words:

• Navigation shall not exceed defined data consumption • Navigation shall contribute with car data to build the knowledge of the serverside system

• Navigation secures the required navigation quality

• Navigation optimizes communication with server-system to alleviate it from over loading All these goals have been addressed through the integration of a run-time autotuner. Important from the HPC perspective is that a navigation client acts as an intelligent source of computational demand, able to regulate the number of requests based on the quality of the navigation and the ability to perform local route computation.

Figure 5.2 shows the tool-flow for configuring the autotuner towards statistically the best performance. The input is the data consumption statistics on estimated data transfer per request type for all service modules, while collected from prior test runs (circa 100 hours) of the navigation app under prevailing commute conditions (in Bratislava The mArgot autotuner (green dots) instead takes into account a Navigation Quality Index (NQI) as a constraint. The NQI is defined as a fuction of the remote routing frequency, which saturates at a point dependent on the current traffic level. The constraint in the proposed experiment is for the NQI to stay above a minimum of 6.0. Furthermore, it uses the same data transfer limit (20MB), a minimal data contribution of 20%, and an estimate of 40 monthly driving hours.

As shown in Figure 5.2 the objectives are well met with an exception of few dots about ideal consumption rate, which can be explained with stochasticity of driver behaviour records (data generated only in the case of excessive acceleration, braking, cornering).

The impact of the NQI constraint selection can be understood through the experiments reported in Figure 5.2. In the experiment, the NQI is set to different values in the range [6.0,9.0]. The results show how relaxing the quality level allows to save computational resources, lowering the total cost of routing, e.g., by 12% when dropping from NQI 8.0 to 7.0. The results also allow to identify the best trade-off from a value-at-cost solution in practical cases, at NQI 6.8. Compared to the baseline autotuning, which achieves an NQI of 6.25, the mArgot autotuner allows a saving of 14% in Figure 18: Performance comparison between baseline and mArgot autotuning computational resources while providing a better quality of service, assuming the use of the best trade-off NQI of 6.8.

Related Works

Many Domain Specific Languages have been proposed for the domain of High Performance computing [START_REF] Giles | Trends in high-performance computing for engineering calculations[END_REF]. A survey of tools and approach in generative programming for high performance computing can be found in [START_REF] Cohen | In search of a program generator to implement generic transformations for high-performance computing[END_REF], whereas a collection of techniques for and examples of modern DSL design can be found in [START_REF] Mernik | Formal and Practical Aspects of Domain-Specific Languages: Recent Developments: Recent Developments[END_REF].

Many of them actually focus on a specific class of HPC applications, either in terms of a specific field such as computational biology [START_REF] Starruß | Morpheus: a user-friendly modeling environment for multiscale and multicellular systems biology[END_REF], or a specific mathematical problem, such as partial differential equations [START_REF] Devito | Liszt: a domain specific language for building portable mesh-based pde solvers[END_REF][START_REF] Rathgeber | Pyop2: A high-level framework for performance-portable simulations on unstructured meshes[END_REF][START_REF] Schmitt | Exaslang: a domain-specific language for highly scalable multigrid solvers[END_REF][START_REF] Membarth | Towards domain-specific computing for stencil codes in hpc[END_REF], often using stencil computation, i.e., common functionalities in a given field implemented in an optimised way and offered through the DSL to the user, possibly with some autotuning capabilities. This approach is effective in terms of combining efficiency of the generated code with familiarity of the language for the domain expert. However, a different solution is needed for each domain or problem, and when HPC experts at a supercomputing center are called for to optimise a customer's code, this forces specialisation, which may be undesirable from the point of view of the supercomputing center management.

Indeed, at least some such centers have attempted to develop a coherent set of languages and tools applicable across multiple application domain within the wider context of HPC [START_REF] Duran | Ompss: a proposal for programming heterogeneous multi-core architectures[END_REF], and large scale research programs have been devoted to this goal [START_REF] Weiland | Chapel, fortress and x10: novel languages for hpc[END_REF]. These DSLs are also related to algorithmic skeleton frameworks [START_REF] González-Vélez | A survey of algorithmic skeleton frameworks: high-level structured parallel programming enablers[END_REF], which however do not take into account the presence of the HPC expert. Model-driven approach have been also proposed [START_REF] Palyart | Mde4hpc: an approach for using model-driven engineering in high-performance computing[END_REF][START_REF] Palyart | Hpcml: a modeling language dedicated to highperformance scientific computing[END_REF], with the goal to abstract the application from the concrete and ever-changing HPC infrastructure. A more general-purpose take on the stencil computation is Terra [START_REF] Devito | Terra: a multi-stage language for highperformance computing[END_REF], which employs Lua to support the creation of DSL compilers which can interoperate with the generated code.

All these efforts have the opposite problem -they force the application developer to acquire skills typical of the HPC experts. Furthermore, they do not deal easily with the large amount of legacy code found in HPC, where applications may remain in use for decades [START_REF] Giles | Trends in high-performance computing for engineering calculations[END_REF][START_REF] Palyart | Hpcml: a modeling language dedicated to highperformance scientific computing[END_REF]. The ANTAREX DSL, on the other hand, is designed exactly to cater to these needs. By creating an external set of LARA aspects, the ANTAREX DSL strategies, the HPC expert can manipulate the code developed by the application domain expert. The DSL code can be as small or as large as needed to reach the performance goals, but never intrudes in the ability of the application developers to maintain or extend their applications. At the same time, by providing access to commonly used computation patterns, code transformations, and autotuning tools, the ANTAREX DSL enables the HPC experts to reuse solutions across different applications, removing the need to reinvent the wheel with every new application domain.

Conclusions

In this paper, we have reported the major concrete outcome of the ANTAREX research project, the ANTAREX DSL. The DSL, based on the AOP language LARA, aims at supporting a well-established collaboration model between users and support staff of HPC centers, where the users develop applications focusing on the functional requirements, and then take advantage of the support staff's expertise on HPC to improve the performance of the applications. With the ANTAREX DSL, this process is simplified by removing the need for the HPC support staff to manually modify the user's code, instead relying on aspects which encode available code transformations as well as providing access to libraries for autotuning, system monitoring, and dynamic recompilation.

We demonstrate the use of our tools in a scenario where the HPC center runs remote navigation for a pool of users, considering both urban and transnational scenarios, as well as in a scenario where the HPC center runs a drug discovery appplication.

Figure 1 :

 1 Figure 1: The ANTAREX Tool Flow, comprised of a compile-time part, where the instrumentation of the functional code is performed based on the strategies expressed in the ANTAREX DSL, and a run-time part, where the generated adaptive binary is monitored during its run on the system and optimised by tuning its parameters, including selection of different code versions.

aspectdefFigure 2 :

 2 Figure 2: Example of LARA aspect to change the types of variables declared inside a given function.

Figure 3 :

 3 Figure 3: Example of LARA aspect that generates different precision mix versions of the same OpenCL kernel.

Figure 4 :

 4 Figure 4: Example of LARA aspect to clone an existing function and change the type of the clone.

Figure 5 :

 5 Figure 5: Example of LARA aspect that generates an alternative version of a function and inserts a mechanism in the code to switch between versions.

Figure 6 :

 6 Figure 6: Excerpt of the C code resulting from the generation of alternative code versions.

Figure 7 :

 7 Figure 7: Example of LARA aspect to replace a function call to a kernel with a call to a dynamically generated version of that kernel.

Figure 8 :

 8 Figure 8: A memoizable C function and its wrapper. Consider a memoizable C function foo as shown in Figure 8. The memoization consists in: 1. the insertion of a wrapper function foo_wrapper and an associated table.

Figure 9 :

 9 Figure 9: An example of LARA aspect defined for the memoization.

Figure 10 :

 10 Figure 10: Example of a LARA aspect for autotuner configuration and code enhancement.

  import antarex.examon.Examon; import lara.code.Logger; aspectdef SimpleExamon input name, ip, topic, $manageFunction, $targetStmt end var broker = new ExamonBroker(ip); var exa = new ExamonCollector(name, topic); // manage the collector on the target function select $manageFunction.body end apply exa.init(broker, $body); exa.start($body); exa.end($body); exa.clean($body); end // get the value and use it in the target stmt exa.get($targetStmt); // get the last stmt of the scope of the target stmt var $lastStmt = $targetStmt.ancestor("scope").lastStmt; // Create printf for time and data var logger = new Logger(); logger.

Figure 11 :

 11 Figure 11: Example of a LARA aspect integrate an ExaMon collector into an application.

  import lat.Lat; import lat.LatUtils; import lat.vars.LatVarOmpThreads; import lara.metrics.EnergyMetric; import lara.metrics.ExecutionTimeMetric; aspectdef UC1Exploration var lat = new Lat('uc1_exploration'); // Compiler flags lat.loadCmaker().cmaker.addFlags('-fopenmp', '-O3', '-march=native'); // Scope for variable change select function{"main"}.body.call{"MeasureOverlap"} end apply lat.setMeasure($call); lat.setScope($body); end // Design space pocketSizes = [5000, 7000, 10000, 12000, 50000]; maxThreads = 5; var samples = new LatVarList("num_pocket_atoms", pocketSizes); var threadValues = new LatVarRange("numThreads", 0, maxThreads, 1, function(x){return Math.pow(2, x);}); var numThreads = new LatVarOmpThreads(threadValues); lat.addSearchGroup([samples, numThreads]); // Metrics lat.addMetric(new ExecutionTimeMetric()); lat.addMetric(new EnergyMetric()); // Tune lat.setNumTests(5); var results = lat.tune(); LatUtils.resultToCsv(results, './'); end

Figure 13 :

 13 Figure 13: LARA code for parameter exploration.

Figure 14 :

 14 Figure 14: Results of the DSL-enabled automatic exploration

  Float (F) a version of the algortihm with reduced precision of internal variablesDouble (D) a standard algorithmFloat_Hoisting (FH) a float version with hoistingDouble_Hoisting (DH) a standard version with hoistingFloat_Hoisting_Memoization (FHM) a float version with hoisting and memoizationDouble_Hoisting_Memoization (DHM) a standard version with hoisting and memoizationTable5.1.1 shows run times with the optimizations provided by different ANTAREX tools for different number of nodes.

Figure 17 :

 17 Figure 17: Autotuning framework application to the navigation app

Figure 5 .

 5 2 shows the effectiveness of the mArgot autotuner compared with the baseline autotuner included in the commercial release of Sygic Truck Navigation v13.7.0 (March 2017). The baseline (red dots) is a simple autotuning implementation that supports only a data transfer limit set at 20MB. It does not support constraints on the navigation quality, providing the maximum quality based on its prediction of data usage.

Figure 19 :

 19 Figure 19: Variation of the compute workload (in routing requests per hour) as a function of the Navigation Quality Index for a scenario of medium traffic commute in Bratislava

  the policy in case of conflicts regarding the same table entry (line 11): replacement or not in case of conflict to the same entry of the table for different parameters of the memoized function, and (2) the size of the table (line 15). After some verifications, not detailed here, on the parameters

	aspectdef Memoize_Method_overloading_ARGS
	input
	aClass,	// Name of a class
	aMethod,	// Name of a method of the class aClass
	pType,	// Name of the selected type
	nbArgs,	// Number of parameters of the method
	fileToLoad, // filename for init of the table, or 'none'
	FullOffLine,// yes for a fully offline strategy
	FileToSave, // filemane to save the table, or 'none'
	Replace,	// Always replace in case of collisions
	approx,	// Number of bits to delete for approximation.
	tsize	// Size of the internal table.
	end
	// Control on the parameters of the aspect: nbArgs in [1,3]
	...
	// Searching the method.
	var MethodToMemoize, found=false;
	select class{aClass}.method{aMethod} end
	apply
	if (! found) {
	found = isTheSelectedMethod($method, nbArgs, pType);
	if (found) MethodToMemoize=$method;
	}
	end
	if (!found)
	{ / * message to the user * /}
	else {
	GenCode_CPP_Memoization(aClass, aMethod, pType, nbArgs,
	fileToLoad, FullOffLine, FileToSave, Replace, approx, tsize);
	call CPP_UpdateCallMemoization(aClass, aMethod, pType, nbArgs);
	}
	end

Table 1 :

 1 Static Metrics

	Strategy	LARA	LARA	Input	Input	Woven	Woven	Delta	Delta
		SLoC-L	Aspects	SLoC-L	Func	SLoC-L	Func	SLoC-L	Func
	ChangePrecision	27	1	12	3	13	3	1	0
	SimpleExamon	20	1	12	3	23	5	11	2
	Multiversion	46	2	12	3	43	5	31	2
	CreateFloatVersion	28	2	12	3	24	3	12	0
	SimpleLibVC	12	1	12	3	39	4	27	1
	HalfPrecisionOpenCL*	93	3	9	1	279	31	270	30
	Total	226	10	69	16	421	51	352	35

Table 2 :

 2 Dynamic Metrics

	File	Selects	Attributes	Actions	Inserts	Native
						SLoC
	ChangePrecision	4	109	2	1	0
	SimpleExamon	4	131	18	7	0
	Multiversion	8	477	27	16	9
	HalfPrecisionOpenCL	125	2211	381	159	31
	CreateFloatVersion	2	170	6	3	0
	SimpleLibVC	7	93	13	8	36
	Total	150	3191	447	194	76

Table 3 :

 3 Characterization of the traffic network graphs used in the experiments

	Graph	Vertices	Edges
	Porto city area	7,316	14,793
	4EU countries 3,567,735 8,610,752

Table 4 :

 4 Betweenness Centrality algorithm run times for the Porto graph (in seconds)

	Compute Nodes	F	FH FHM	D	DH DHM
	1	1.20 1.15	1.12 1.34 1.27	1.25
	2	0.66 0.65	0.64 0.75 0.72	0.71
	4	0.41 0.42	0.41 0.47 0.45	0.44

Table 5 :

 5 Betweenness Centrality algorithm run times for the 4EU graph (in 10 3 seconds)

	Compute Nodes	F	FH	FHM	D	DH	DHM
	1 294.575 283.534 282.233 340.523 325.428 315.872
	2 153.296 150.123 147.966 175.940 170.587 167.524
	4	79.484	79.233	75.314	89.905	87.206	86.106
	8	40.811	39.899	38.711	46.031	44.736	44.258
	16	21.305	21.070	20.247	23.725	22.958	22.493
	32	11.122	10.809	10.590	12.228	11.884	11.791
	64	5.922	5.726	5.677	6.489	6.284	6.092

www.top500.org, November

2 www.green500.org, November 2018[START_REF] Silvano | Antarex: A dsl-based approach to adaptively optimizing and enforcing extra-functional properties in high performance computing[END_REF] Summit is also ranked 3rd in the Green500, whereas Shoubu is a comparatively smaller machine

https://github.com/specs-feup/clava

Complete working versions for all examples can be found in https://github.com/specs-feup/specslara/tree/master/2018%20DSD

https://gitlab.com/margot_project/core

LAT is available at https://github.com/specs-feup/LAT-Lara-Autotuning-Tool

Acknowledgements

The ANTAREX project is supported by the EU H2020 FET-HPC program under grant 671623. The IT4Innovations infrastructure is supported by the Large Infrastructures for Research, Experimental Development and Innovations project "IT4Innovations National Supercomputing Center -LM2015070". The CINECA infrastructure is supported by the EU and the Italian Ministry for Education, University and Research (MIUR) under the PRACE project.

the number of threads and the number of atoms of the pocket size to explore, and creates versions of the code for all combinations of these parameters. For each combination, the strategy will compile, run and collect the results, which are aggregated into a CSV file at the end.

In more detail, line 10 creates an instance of the Lat class, which will be used throughout the strategy. Since the code has been instrumented with OpenMP pragmas, line 13 sets some necessary compilation flags. Line 16 select the section of code we want to explore (i.e., the call to the function MeasureOverlap), and lines 17-20 set the call as the place in the code where we want to measure the metrics, and the body of the function where the call appears as the place where the parameters will be changed.

Lines 26-28 create the ranges of values we want to explore: a set of predefined values for the variable num_pocket_atoms, and a set of number of threads that start at 1 and double at each step.

Lines 33-34 set the metrics we want to measure (i.e., execution time and energy), line 37 sets the number of times we should repeat the measures for each version of the code, and line 38 executes the exploration. Finally, line 39 collects the results to a CSV file.

Performance

We have applied the auto-parallelization and exploration strategy to the C version of the MeasureOverlap function, considering pocket sizes from 5.000 atoms to Betweenness Centrality algorithm, which also poses the biggest challenge for computation efficiency. The Betweenness Centrality algorithm can alternatively be used as a tool for a what-if analysis on a real traffic network [START_REF] Hanzelka | Betweenness propagation[END_REF].

We have applied two source-to-source transformations to the Betweenness Centrality algorithm: custom precision and memoization. The original application contains 770 lines of the code (C code + the header files). The LARA aspect in Figure 15 shows the main strategy for the custom precision transformation. It replaces the references to double type by a new symbol and replaces all the syntactic expressions that depend on the double type (for instance, mathematical functions) by new symbols. The custom precision library generates several files that allow to have different definitions for the created symbols: a file that defines the created parameter symbols by the original values (double here); a file that defines the created symbols to have a float version, generated by the call to the CustomPrecisionGenParametersFor('float') aspect; a file that contains a template to have another type representation. The last one must be updated by the user to a specific type representation (for example: half precision, fixed point representation, etc.). After applying the custom precision strategy, the final number of lines of the Betweenness Centrality application is equal to 870. For the memoization strategy, which was used as a proof-of-concept for this technology, the main DSL aspect is shown in Figure 16. It applies the memoization to the method computeMetric of the Betweenness class.


François Boulier

François Lemaire

Marc Moreno Maza

Adrien Poteaux

M Moreno

A Short Contribution to the Theory of Regular Chains

This paper contains short contributions to the theory of regular chains which follow a recent JSC paper by the same authors. These contributions apply to both the nondifferential and the differential context.

Introduction

This paper investigates consequences of the recent [START_REF] Boulier | An Equivalence Theorem for Regular Differential Chains[END_REF] which contains important characterizations of regular chains in the nondifferential case [START_REF] Boulier | An Equivalence Theorem for Regular Differential Chains[END_REF]Theorem 21] and in the differential one [START_REF] Boulier | An Equivalence Theorem for Regular Differential Chains[END_REF]Theorem 37]. This paper: 1) improves an algorithm for computing the normal form of a differential fraction f /g (f and g being two differential polynomials) modulo a regular differential chain A, which was presented in [4, Figure 2]; 2) provides a new proof of the well-known fact that regular chains decide membership to the ideals that they define, which highlights the main theoretical argument underlying this result; 3) extend some results of [START_REF] Boulier | An Equivalence Theorem for Regular Differential Chains[END_REF].

Our new normal form algorithms are better than the former versions for two reasons: 1) they succeed in computing normal forms whenever these ones exist (see Section 2.2.1 for examples); and 2) they are conceptually simpler, though their simplicities hide subtleties (see Section 2.2.2).

The main computational tool involved in this paper is the resultant of a polynomial w.r.t. a regular chain (which is a triangular set of polynomials with good properties). The main recent theoretical argument involved in our proofs is the implication 1 ⇒ 4 of [START_REF] Boulier | An Equivalence Theorem for Regular Differential Chains[END_REF]Theorem 21] (recalled as Theorem 1 in this paper):

Let A be a regular chain, a the ideal that it defines and f be a polynomial of some polynomial ring R. Then, f is regular in R/a if and only if the resultant of f w.r.t. A is not zero. Non expert readers may wonder about the novelty of the above quotation since seemingly similar statements can be found in earlier papers. Indeed, the above statement already appears in [6, Theorem 1]. It was formerly proved, in the zerodimensional case, in [7, Lemma 4]. As far as we know, the first occurrence of a similar statement is [13, Corollary 4, page 150] but its exact relationship with more recent works is difficult to clarify since it is formulated in terms of proper ascending chains. Deciding the regularity of a polynomial w.r.t. a regular chain is also a much studied problem. A careful survey on this question can be found in [5, Section 8].

Non expert readers may also wonder about the relevance of the above quotation. Indeed, the membership test to the ideal a defined by a triangular set A, by means of the pseudoremainder algorithm, seems much more important. Actually, this membership test was established by Ritt in [START_REF] Ritt | Differential Algebra[END_REF], in the easy case of a prime ideal a. In the case of an ideal a which is not prime, this membership test only holds for regular chains -not for general triangular sets (Proposition 11). We believe that any proof for this membership test to a first requires to secure a nontrivial regularity property equivalent to the above quotation. In particular, we believe that many formerly published proofs are -if not flawed -at least implicitly using a nontrivial argument. For more details on this point, we refer to [START_REF] Boulier | An Equivalence Theorem for Regular Differential Chains[END_REF]Section 8]. Our new proof of Proposition 9 aims at highlighting the key role of the above quotation.

The paper is organized as follows. The nondifferential case is addressed in Section 2. Classical results are recalled in Section 2.1. Then the new normal form algorithms are presented in Section 2.2. In particular, Algorithms 1 and 2 as well as Propositions 4 and 7 are new results. The new proof of the membership test is given in Section 2.3. The differential case is addressed in Section 3. Classical results are recalled in Section 3.1. Then the new normal form algorithms are presented in Section 3.2. In particular, Algorithms 3, 4 and 5 as well as Propositions 16 and 18 are new results. Proposition 16 is notably important since it completes [START_REF] Boulier | An Equivalence Theorem for Regular Differential Chains[END_REF]Theorem 37].

The Nondifferential Case

Classical Results

In this section, A = {p 1 , . . . , p n } denotes a subset of the polynomial ring R = K[t 1 , . . . , t m , x 1 , . . . , x n ] which is triangular in the sense that deg(p k , x k ) > 0 and deg(p k , x ) = 0 for all indices 1 ≤ k ≤ n and k < ≤ n. Denote i k the initial of p k i.e. the leading coefficient of p k w.r.t. x k , for 1 ≤ k ≤ n. Let f and g be two polynomials of S[x], where S is a unitary ring of characteristic zero:

f = a m x m + • • • + a 1 x + a 0 , g = b n x n + • • • + b 1 x + b 0 .
If f or g is zero, then the resultant of f and g is taken to be zero. Assume that f and g are nonzero and that at least one of them has positive degree. Then, the resultant of f and g is the determinant of the Sylvester matrix S(f, g) of f and g, which has dimensions (m + n) × (m + n) and rows, from top down x n-1 f, . . . , x f, f, x m-1 g, . . . , x g, g. See [2, 4.2, page 105].

Definition 1.

Let A be a possibly empty triangular set and f ∈ R. The resultant of f by A, denoted res(f, A), is defined as follows:

1. if A = ∅ then res(f, A) = f ; 2. if A = {p 1 , . . . , p n } then res(f, A) = res(res(f, p n , x n ), {p 1 , . . . , p n-1 }).
We could not find any reference for the following Proposition, which is an easy consequence of basic properties of resultants. Observe that this Proposition is algorithmic, thanks to extended versions of algorithms for computing pseudoremainder subresultant sequences. See [START_REF] Ducos | source of the axiom package prs.spad[END_REF].

Proposition 1. Let f be a polynomial and A be a triangular set of R. Then there exist polynomials

u, v 1 , v 2 , . . . , v n of R such that u f = res(f, A) + v 1 p 1 + v 2 p 2 + • • • + v n p n . (1)
Moreover, if f does not depend on x k , . . . , x n for some 1 ≤ k ≤ n, then there exists a formula (1)

such that u, v 1 , . . . , v k-1 do not depend on x k , . . . , x n and v k = • • • = v n = 0.
Proof. Apply n times the fact that the resultant of two polynomials is in the ideal generated by these polynomials. See [2, Proposition 4.18] or [5, Lemma 3]. For the second statement, it is sufficient to notice that, if f does not depend on x k , . . . , x n , then res(f, A) is a power of res(f, {p 1 , . . . , p k-1 }).

To each (non necessarily triangular) set A of polynomials of R \ K, one associates the ideal a = sat (A) = (A) :

(i 1 • • • i n ) ∞ . Definition 2. A triangular set is said to be a regular chain if the initial i k of p k is regular in R/ sat (p 1 , . . . , p k-1 ) for 2 ≤ k ≤ n.
The following characterization is [START_REF] Boulier | An Equivalence Theorem for Regular Differential Chains[END_REF]Theorem 21].

Theorem 1. Let A be a triangular set. The following conditions are equivalent:

1. A is a regular chain; 2. for each 2 ≤ ≤ n and each 1 ≤ k ≤ n we have res(i , {p k , . . . , p n }) regular in R/a; 3. for each 2 ≤ ≤ n we have res(i , A) = 0; 4. for each f ∈ R, f is regular in R/a if and only if res(f, A) = 0;

5. for each f ∈ R, f is regular in R/a for each 1 ≤ k ≤ n, res(f, {p k , . . . , p n }) is regular in R/a.
Normal forms of rational fractions modulo regular chains are defined and studied in [4, Definition 5.1, Proposition 5.2], in the framework of differential algebra. In the nondifferential context, these results can be restated as follows.

Proposition 2. Let A be a regular chain and f /g be a rational fraction with g regular in R/a. The normal form of f /g modulo A exists: it is the unique rational fraction p/q such that 1. deg(p,

x k ) < deg(p k , x k ) for 1 ≤ k ≤ n; 2. q ∈ K[t 1 , . . . , t m ];
3. f /g and p/q are equal in the total ring of fractions of R/a.

An algorithm for computing the normal form of a rational fraction is given in [4, Figure 2]. This algorithm relies on the computation of the inverse of a polynomial modulo a regular chain. The following definition and proposition restate [4, Definition 4.1 and Proposition 4.3]. Definition 3. Let A be a regular chain and g be a nonzero polynomial of R. An inverse of g modulo A is any rational fraction p/q such that q = 0, q ∈ K[t 1 , . . . , t m ] and g p = q in R/a. Proposition 3. Let A be a regular chain and g be a nonzero polynomial of R. The polynomial g is regular in R/a if and only if it admits an inverse modulo A.

A New Normal Form Algorithm

Algorithm 1, respectively 2, computes the normal form of a polynomial, respectively a rational fraction, modulo a regular chain A. These algorithms are new. Their proofs essentially rely on the following Proposition, which actually provides an algorithm for computing inverses, since it relies on Proposition 1, which is itself of algorithmic nature. Proposition 4. Let A be a regular chain and g be a regular element of R/a. Let u be a polynomial such that u g = res(g, A) in R/a. Then u/ res(g, A) is an inverse of g modulo A.

Proof. The existence of u is guaranteed by Proposition 1. The fact that res(g, A) ∈ K[t 1 , . . . , t m ] is an easy property of the iterated resultant. The fact that res(g, A) = 0 follows from the regularity of g and the implication 1 ⇒ 4 of Theorem 1.

Algorithm 1: NFpoly(f, A)

input : a polynomial f , and Algorithm 2: NF(f /g, A) input : a fraction f /g, and a regular chain A = {p 1 , . . . , p n } with g regular in R/a output: the normal form of f /g modulo A 1 compute u/r an inverse of g modulo A ;

A = {p 1 , . . . , p n } a regular chain output: the normal form of f modulo A 1 if A is empty then 2 return f 3 else 4 compute i α n f = f + v n p n ; /* by computing f = prem(f, p n , x n ) */ 5 compute u/
/* Proposition 4 */ 2 return (1/r) × NFpoly(uf, A) 2.2.1. Improvement Compared to the Former Version. As pointed out in [4, Appendix A] the normal form algorithm may fail to compute the normal form because of the failure of the function used for computing inverses modulo regular chains. Consider the regular chain Example 4] and consider the polynomial g = (x 1 -1) x 2 + 1. The polynomial g is invertible in R/a since res(g, {p 2 , p 1 }) = 45. Moreover, the inverse of g in R/a is easily deduced from the following formula:

A = {p 3 := x 3 -x 2 -x 1 , p 2 := x 2 2 -x 3 1 , p 1 := (x 1 -1) (x 1 + 1) (x 2 1 -2)} from [4, Appendix A,
-3 (2 x 1 + 3) (2 x 2 1 + 1) (x 1 x 2 -x 2 -1) g = 45 -(12 x 3 1 + 18 x 2 1 + 6 x 1 + 9) (x 1 -1) 2 p 2 -(12 x 4 1 -6 x 3 1 + 18 x 2 1 -3 x 1 + 18) p 1 .
Thus, using our new algorithm, we have

NF(1/g, A) = - (2 x 1 + 3) (2x 2 1 + 1) (x 1 x 2 -x 2 -1) 15 •
However, calling [START_REF] Boulier | A Normal Form Algorithm for Regular Differential Chains[END_REF]NF](1/g, A) yields an error for the following reason: a call to Algorithm [4, Inverse](g, A) triggers a call to Algorithm [4, AlgebraicInverseNonZero](g, A), which itself triggers a call to Algorithm [4, ExtendedEuclideanAlgorithm](g, p 2 , x 2 , A) which fails since the initial of g (which Proposition 5. Let A be a regular chain and A ⊂ A be nonempty. Then A is a regular chain. Moreover, every f ∈ R which is regular in R/a is regular in R /a .

is x 1 -1) is not invertible in R/a.
It is tempting to try to generalize this proposition as follows: "Let f /g be a fraction with g regular in R/a. Then NF(f /g, A) = NF(NF(f /g, A ), A)". Unfortunately, this generalization is false. A counterexample is suggested by the Remark following [START_REF] Boulier | An Equivalence Theorem for Regular Differential Chains[END_REF]Proposition 17]. Take

p 1 = (x 1 -1)(x 1 -3) , p 2 = x 2 -10 x 1 , p = x 1 + x 2 -31 . Take A = {p 1 , p 2 }, A = {p 1 } and f /g = 1/p. Both A and A are regular chains. We have NF f g , A = 11 x 1 -13 40 and NF f g , A = x 2 -x 1 -27 (x 2 -28) (x 2 -30) •
However, res((x 2 -28) (x 2 -30), A) = 0, proving that this polynomial is a zerodivisor in R/a by the implication 1 ⇒ 4 of Theorem 1. Thus NF(NF(f /g, A ), A) is not defined. In this example, the fact that deg(p 1 , x 1 ) > 1 is important, the next propositions (the second one is new) show.

Proposition 6. Let g be a polynomial and f be a polynomial such that deg(f, x) = 1. Denote n = deg(g, x) and r = prem(g, f, x). Then res(g, f, x) = (-1) n r.

Proof. By [2, Lemma 4.17, page 107] we have res(g, f, x) = (-1) n res(f, r, x). Since deg(f, x) = 1 we have deg(r, x) = 0. Thus the Sylvester matrix defined by f and r has dimension 1 × 1 and involves r as single element. Thus res(g, f, x) = r and the Proposition is proved.

Proposition 7. Let A be a regular chain, p k ∈ A be such that deg(p k , x k ) = 1 and g be a differential polynomial of R. Then g is regular in R/a if and only if res(g,

p k , x k ) is regular in R/a.
Proof. Let r = prem(g, p k , x k ) denote the pseudoremainder of g by p k and n = deg(g, x k ). By Proposition 6, we have res(g, p k , x k ) = (-1) n r.

We claim that g is regular in R/a if and only if r is so. Let p be any associated prime ideal of a. There exists a power h of the initial of p k such that h g = r (mod p). Since a is saturated by the initials of A, we have h / ∈ p hence g ∈ p if and only if r ∈ p. By [14, IV, 7, Corollary 3 to Theorem 10], the claim and the proposition are proved.

The Membership Problem

All propositions stated in this section are well-known: the only difficult one is Proposition 9, which already appears in [1, (i) ⇒ (iii)]. The originality of this section is the proof of that Proposition which involves a single nontrivial argument: the implication 1 ⇒ 4 of Theorem 1. Proposition 8. Let A be a triangular set and f be any polynomial of R. Denote g = prem(f, A). Then

deg(g, x k ) < deg(p k , x k ) (1 ≤ k ≤ n) . (2)
Moreover, there exist polynomials v 1 , v 2 , . . . , v n and a power product h of initials of A and such that

h f = g + v 1 p 1 + v 2 p 2 + • • • + v n p n . ( 3
) Proposition 9. Assume A is a regular chain and f is any polynomial of R. If f ∈ a then prem(f, A) = 0.
Proof. Observe that if a = (0) then the proposition holds so that we only need to consider the case a = (0). Observe also that, if f ∈ a then prem(f, A) ∈ a. Thus, if f ∈ a is such that prem(f, A) = 0 then there exists such a f reduced w.r.t. A i.e. such that deg(f,

x k ) < deg(p k , x k ) for all 1 ≤ k ≤ n.
We thus assume that there exists a nonzero polynomial f ∈ a which is reduced w.r.t. A and seek a contradiction. Since f ∈ a, there exists nonnegative integers α 1 , . . . , α n and polynomials v 1 , . . . , v n such that

i α1 1 • • • i αn n f = v 1 p 1 + v 2 p 2 + • • • + v n p n . (4)
For 1 ≤ k ≤ n, denote r k = res(i k , A) and u k a polynomial such that u k i k = r k (mod a), according to Proposition 1. Multiply both sides of (4) by u α1

1 u α2 2 • • • u αn n and denote h f = r α1 1 • • • r αn n .
Since the initials i k are regular in R/a, their resultants r k are different from zero by the implication 1 ⇒ 4 of Theorem 1. Thus there exists a nonzero element h ∈ K[t 1 , . . . , t m ] (one may take h = h f ) and polynomials w 1 , . . . , w n not all zero, such that

h f = w 1 p 1 + w 2 p 2 + • • • + w n p n F . (5)
To any such formula F , one associates the index j(F ) defined as

j(F ) = max{ ∈ N | ∃ k ∈ [1, n] , deg(w k p k , x ) > 0} . (6)
The index j is well defined since the polynomials w k are not all zero. The definition of j implies that w j+1 = • • • = w n = 0. Moreover, since x j+1 , . . . , x n appear nowhere in the right hand side of (5), they cannot appear either in the left hand side so that f is free of x j+1 , . . . , x n . Last observe that j > 1 since, otherwise, deg(f, x 1 ) would be greater than or equal to deg(p 1 , x 1 ), which would contradict the assumption that f is reduced w.r.t. A. Denote E the set of all possible pairs (f, F ), where f ∈ a is a nonzero polynomial reduced w.r.t. A and there exists some nonzero h ∈ K[t 1 , . . . , t m ] such that (5) holds. Our assumptions imply that E is not empty. Among all possible pairs (f, F ) ∈ E , fix one such that j = j(F ) is minimal and denote d = deg(p j , x j ).

Denote w k = prem(w k , p j , x j ) for 1 ≤ k < j. There exist nonnegative integers β k and polynomials q k such that Denote = max(β 1 , . . . , β j-1 ). Multiply both sides of (5) by i j and recall that w j+1 = • • • = w n = 0. One gets a formula:

i β k j w k = w k + q k p j , (1 ≤ k < j) . (7)
i j h f = i j w 1 p 1 + • • • + i j w j-1 p j-1 + i j w j p j . (8)
Denote k = -β k for 1 ≤ k < j. From ( 7) and ( 8) one gets:

i j h f = i 1 j (i β1 j w 1 ) p 1 + • • • + i j-1 j (i βj-1 j w j-1 ) p j-1 + i j w j p j , (9)
= i 1 j (w

1 + q 1 p j ) p 1 + • • • + i j-1 j (w j-1 + q j-1 p j ) p j-1 + i j w j p j , (10)
= i 1 j w 1

w 1 p 1 + • • • + i j-1 j w j-1 w j-1 p j-1 + (i j w j + i 1 j q 1 p 1 + • • • + i j-1 j q j-1 p j-1 ) w j p j . ( 11
)
We have deg(w k , x j ) < d for 1 ≤ k < j and deg(i j h f ) < d. Thus w j = 0. The initial i j does not depend on x j , . . . , x n . Thus Proposition 1 permits us to express the resultant r j as

u j i j = r j + z 1 p 1 + • • • + z j-1 p j-1 , (12)
where u j , z 1 , . . . , z j-1 are polynomials which do not depend on x j , . . . , x n . Multiplying both sides of ( 11) by u j , one gets

(u j i j ) h f = u j w 1 p 1 + • • • u j w j-1 p j-1 . (13)
Plugging ( 12), one gets

(r j + z 1 p 1 + • • • + z j-1 p j-1 ) h f = u j w 1 p 1 + • • • + u j w j-1 p j-1 . (14)
The polynomials u j , z 1 , . . . , z j-1 are free of x j . The polynomials f, w 1 , . . . , w j-1 have degree in x j less than d. The polynomials r j , h are nonzero elements of K[t 1 , . . . , t m ]. Denote h = r j h. Expanding the left hand side of ( 14) and distributing over the right hand side (there may be many different ways to perform this operation), one eventually gets a formula

h f = w 1 p 1 + • • • + w j-1 p j-1 ( 15
)
where h is a nonzero element of K[t 1 , . . . , t m ] and the w k are polynomials such that deg(w k , x j ) < d for 1 ≤ k < j. All polynomials occurring in (15) are free of x j+1 , . . . , x n .

Let now e be some degree (0 ≤ e < d). Denote f e the coefficient of x e j in f and w k,e the coefficient of x e j in w k , for 1 ≤ k < j. Since the polynomials p 1 , . . . , p j-1 are free of x j , the following formula holds for each degree 0 ≤ e < d. Since f is nonzero, the formula holds for some e such that f e is nonzero.

h f e = w 1,e p 1 + • • • + w j-1,e p j-1 F . (16)
Formula ( 16) implies in particular that h f e ∈ a. Since h is a nonzero element of K[t 1 , . . . , t m ], it is a regular element of R/a, by the implication 1 ⇒ 4 of Theorem 1 and the fact that res(h , A) = h . Thus f e ∈ a. Since f e is a coefficient of a polynomial f reduced w.r.t. A, it is reduced w.r.t. A also.

The pair (f e , F ) thus belongs to E and is such that j(F ) < j(F ). This contradiction with the minimality of j(F ) proves that E must be empty and concludes the proof of the proposition.

The two following Propositions already appear in [1, Theorem 6.1]. Their proofs are elementary.

Proposition 10. Let A be a triangular set and f be a polynomial of R. Assume that prem(f, A) = 0 if f ∈ a. Then A is a regular chain.

Proof. Since prem(f, A) = 0 for each f ∈ a and prem(1, A) = 0, the ideal a is proper. Denote R = K[t 1 , . . . , t m , x 1 , . . . , x ] and a = sat (p 1 , . . . , p ) for each 1 ≤ ≤ n.

We assume the regular chain condition is satisfied up to index k < n but not at index k + 1 i.e. i is regular in R/a -1 for each 2 ≤ ≤ k and i k+1 is a zerodivisor in R/a k . We prove that there exists some f ∈ a such that prem(f, A) = 0.

Denote b

k = a k : i ∞ k+1 = (p 1 , . . . , p k ) : (i 1 • • • i k+1 ) ∞ . Since i k+1 is a zerodivisor in R/a k we have a k b k . Thus there exists some f ∈ b k , f / ∈ a k . Since the defining polynomials of a k and b k belong to R k , we may assume f ∈ R k . Since f / ∈ a k , we have prem(f, A k ) = 0 (Proposition 8). Since f ∈ R k , we have prem(f, A k ) = prem(f, A). Since f ∈ b k ⊂ a, the Proposition is proved.
The following proposition completes [START_REF] Boulier | An Equivalence Theorem for Regular Differential Chains[END_REF]Theorem 21].

Proposition 11. Let A be a triangular set. The following conditions are equivalent:

1. A is a regular chain; 6. for each f ∈ R, f ∈ a if and only if prem(f, A) = 0.

Proof. First observe that, for any set A, if prem(f, A) = 0 then f ∈ a (Proposition 8). The implication 1 ⇒ 6 is then Proposition 9. The converse implication is Proposition 10.

The Differential Case

Classical Results

Reference books are the ones of Ritt and Kolchin [START_REF] Ritt | Differential Algebra[END_REF][START_REF] Kolchin | Differential Algebra and Algebraic Groups[END_REF]. Let R = K{U } be a differential polynomial ring where K is a differential field of characteristic zero, U is a finite set of differential indeterminates u k , endowed with a finite set of derivations. Let Θ denote the multiplicative monoid of derivation operators, generated by the derivations and Θ * ⊂ Θ the set of proper derivation operators. Assume the infinite set of derivatives Θ U is ordered w.r.t. a ranking [9, I, 8, page 75] so that, given any differential polynomial f ∈ R \ K, its leading derivative ld f (called leader by Kolchin), its initial and its separant, which is the partial derivative of f w.r.t. its leading derivative, are well defined.

In the sequel, we will have to consider sets of differential polynomials as particular cases of sets of plain polynomials, in order to apply the results of Section 2. By a triangular set of differential polynomials {p 1 , . . . , p n } we mean a set of differential polynomials of R \ K having pairwise distinct leading derivatives. In order to apply the results of the former sections and fit to their notations, we assume moreover that ld p 1 < ld p 2 < • • • < ld p n . These leading derivatives then correspond to the variables x 1 , x 2 , . . . , x n . In particular, the numbering of the x is imposed by the ranking. The other derivatives occurring in the differential polynomials correspond to t 1 , . . . , t m . This being understood, there is no ambiguity in a statement such as "the triangular set A of differential polynomials is a regular chain". Similarly, if f is any differential polynomial, the differential polynomial res(f, A) is well defined, by means of the iterated resultant.

A differential polynomial f is said to be partially reduced w.r.t. a differential polynomial p / ∈ K if f does not depend on any proper derivative of the leading derivative v of p [9, I, 9, page 77]. It is said to be fully reduced w.r.t. p if, in addition, deg(f, v) < deg(p, v). If f is any element and A is any subset of R \ K, thanks to the pseudo-remainder algorithm, it is easy to compute a differential polynomial g, partially reduced w.r.t. A (i.e. w.r.t. each element of A) and such that h f = g (mod [A]) where [A] denotes the differential ideal generated by A in R and h is some power product of separants of A. The differential polynomial g is called a partial remainder of f by A. It is as well easy to compute a differential polynomial g, fully reduced w.r.t. A and such that h f = g (mod [A]) where h is some power product of initials and separants of A. The differential polynomial g is called a full remainder of f by A. See [9, I, 9] or [11, I, 6].

Consider a triangular set A of differential polynomials of R. Let L denote the set of leading derivatives of A and N = ΘU \ ΘL the possibly infinite set of the elements of ΘU which are not derivatives of any element of L (the derivatives "under the stairs" of A). Then K[N ∪ L] ⊂ R is the ring of the differential polynomials partially reduced w.r.t. A.

Uppercase gothic letters denote differential ideals while lowercase ones denote nondifferential ones. In particular, we denote A the differential ideal [A] : h ∞ of R where h denotes the product of the initials and separants of A. We denote a = sat (A) the nondifferential ideal of R defined by A, viewed as a plain triangular set.

Definition 4. A regular chain A is said to be squarefree if, for each 1 ≤ k ≤ n, the separant s k of p k is regular in R/a. Definition 5. A triangular set A of pairwise partially reduced differential polynomials is said to be coherent if, for each pair {p 1 , p 2 } of polynomials of A whose leading derivatives θ 1 u, θ 2 u are the derivatives of some common differential indeterminate u ∈ U , we have ∆(p 1 , p 2 ) ∈ sat (A ) where θ 12 = lcm(θ 1 , θ 2 ), A = {p ∈ ΘA | ld p < θ 12 u} and, denoting s 1 and s 2 the separants of p 1 and p 2 ,

∆(p 1 , p 2 ) = s 2 θ 12 θ 1 p 1 -s 1 θ 12 θ 2 p 2 .
The pairs considered in the above definition may exist only if the number of derivations is greater than or equal to 2. Thus, in the ordinary differential case, every triangular set A of pairwise partially reduced differential polynomials is coherent. Definition 6. A triangular set A of pairwise partially reduced differential polynomials is said to be a regular differential chain if it is a coherent squarefree regular chain.

The following characterization is [START_REF] Boulier | An Equivalence Theorem for Regular Differential Chains[END_REF]Theorem 37].

Theorem 2. Let A be a triangular set of differential polynomials pairwise partially reduced. The following conditions are equivalent:

1. A is a regular differential chain; 2. A is coherent and, for each 1 ≤ ≤ n and each 1 ≤ k ≤ n we have res(i , {p k , . . . , p n }) and res(s , {p k , . . . , p n }) regular in R/a; 3. A is coherent and res(i , A) = 0 and res(s , A) = 0 for each 1 ≤ ≤ n; 4. for each f ∈ R, f is regular in R/A for any triangular finite subset A ⊂ ΘA such that res(f,

A ) ∈ K[N ∪ L], res(f, A ) = 0.
The next proposition is [START_REF] Boulier | An Equivalence Theorem for Regular Differential Chains[END_REF]Proposition 35]. It is an easy corollary to Rosenfeld's Lemma [START_REF] Rosenfeld | Specializations in differential algebra[END_REF].

Proposition 12. Let A be a regular differential chain. Then

A ∩ K[N ∪ L] = a ∩ K[N ∪ L].
The following proposition is well known. Its first statement is a particular case of a result proved in [START_REF] Boulier | Representation for the radical of a finitely generated differential ideal[END_REF][START_REF] Morrison | The Differential Ideal [P ] : M ∞[END_REF]. The second statement goes back to [11, I, 16; and IX, page 166]. Proposition 13. Let A be a regular differential chain. Then A is a radical differential ideal. In particular, it is a finite intersection of differential prime ideals

A = P 1 ∩ • • • ∩ P r , which is unique if it is minimal.
The remaining part of this section is dedicated to normal forms. The following proposition restates [4, Definition 5.1 and Proposition 5.2]. Proposition 14. Let A be a regular differential chain and f /g be a fraction of differential polynomials with g regular in R/A. The normal form of f /g modulo A exists: it is the unique fraction of differential polynomials p/q such that 1. p is fully reduced w.r.t. A; 2. q ∈ K[N ]; 3. f /g and p/q are equal in the total ring of fractions of R/A.

The

Definition 7.

Let A be a regular differential chain and g be a nonzero differential polynomial of R. An inverse of g modulo A is any fraction p/q of differential polynomials such that p ∈ K[N ∪ L], q is a nonzero element of K[N ] and g p = q in R/A.

Proposition 15. Let A be a regular differential chain and g be a differential polynomial of R. Then g admits an inverse modulo A if and only if it is regular in R/A.

A New Normal Form Algorithm

The following proposition is new and completes Theorem 2.

Proposition 16. Let A be a regular differential chain. Then, 5. for each f ∈ R, f is a zerodivisor in R/A for any triangular finite subset A ⊂ ΘA such that res(f, A ) ∈ K[N ], res(f, A ) = 0.

Proof. First observe that there exists a triangular finite subset A of ΘA such that res(f,

A ) ∈ K[N ].
Its finiteness follows from the fact that rankings are well-orderings [9, I, 8].

The bottom-up implication. Assume that, for any triangular finite subset A ⊂ ΘA such that res(f, A ) ∈ K[N ] we have res(f, A ) = 0. Then f is a zerodivisor in R/A by the implication 1 ⇒ 4 of Theorem 2.

The top-down implication. Let A be a triangular finite subset of ΘA such that res(f, A ) ∈ K[N ]. Assume f is a zerodivisor in R/A. Then there exists some differential polynomial g / ∈ A such that f g ∈ A.

Since A is radical, it is a finite, minimal, intersection of differential prime ideals A = P 1 ∩• • •∩P r by Proposition 13. Since f is a zerodivisor in R/A, renumbering the P i if needed, there exists some index 1 ≤ s ≤ r such that f ∈ P i if and only if 1 ≤ i ≤ s. The differential polynomial g thus satisfies: 1) g ∈ P j for all s < j ≤ r; and 2) g / ∈ P i for some 1 ≤ i ≤ s. Let ḡ denote a full remainder of g by A (one may take ḡ = prem(g, A )). For each 1 ≤ i ≤ r, we have g ∈ P i if and only if ḡ ∈ P i . In particular, ḡ / ∈ A hence is nonzero. Consider now the resultant r = res(f, A ). We have r ∈ P i for 1 ≤ i ≤ s by Proposition 1.

The product r ḡ ∈ P i for each 1 ≤ i ≤ r hence it belongs to A. Since ḡ ∈ K[N ∪ L] and r ∈ K[N ], we have r ḡ ∈ A ∩ K[N ∪ L]. Thus, by Proposition 12, r ḡ ∈ a, in some finitely generated polynomial subring of K[N ∪ L]. Since A is a regular chain, we have prem(r ḡ, A) = 0, by Proposition 9. Since r ∈ K[N ] and ḡ is reduced w.r.t. A, the product r ḡ is reduced w.r.t. A i.e. prem(r ḡ, A) = r ḡ. Since ḡ = 0, we must have r = 0.

Proposition 17. Let g be a differential polynomial of R and f be a differential polynomial of R \ K with leading derivative v.

Then for any derivation operator θ ∈ Θ * we have res(g, θf, θv) = ± prem(g, θf, θv). In particular, there exists a power h of the separant of f such that h g = ± res(g, θf, θv) (mod (θf )).

Proof. Since θ is a proper operator, θf has leading derivative θv, degree one in θv and the separant of f as leading coefficient w.r.t. θv. The Proposition then follows from Propositions 6 and 8.

The following proposition is new and provides an algorithm for computing the inverse of a differential polynomial which is regular in R/A, since it relies on Proposition 1, which is itself of algorithmic nature.

Proposition 18. Let A be a regular differential chain and g be a differential polynomial of R. Let A be a triangular subset of ΘA such that res(g, A ) ∈ K[N ]. Then 1. g is regular in R/A if and only if res(g, A ) = 0; 2. there exists a differential polynomial u ∈ K[N ∪ L] such that u g = res(g, A ) in R/A; 3. if g is regular in R/A, then u/ res(g, A ) is an inverse of g modulo A.

Proof. Assume g is regular in R/A. Then res(g, A ) = 0 by the implication 1 ⇒ 4 of Theorem 2. Conversely, assume g is a zerodivisor in R/A. Then res(g, A ) = 0 by Proposition 16. The first item is thus proved.

The existence of a differential polynomial u such that u g = res(g, A ) in R/A follows from Proposition 1. The fact that u ∈ K[N ∪ L] comes from the fact that the iterated resultant by A can be decomposed into elementary resultants, either by elements of A, or by proper derivatives of elements of A. The polynomial u is the product of polynomials (say) u i arising from these elementary resultant computations. In the first case, u i ∈ K[N ∪ L]; in the second case, u i is a separant of some element of A, by Proposition 17, hence belongs also to K[N ∪ L]. The second item is thus proved.

The third item follows from Definition 7 and the two former items.

Algorithms 3 and 4 provide algorithms for computing the normal form of a differential fraction, which are directly derived from the theoretical construction achieved so far. input : a differential fraction f /g, and A a differential regular chain with g regular in R/A output: the normal form of f /g modulo A 1 compute u/r an inverse of g modulo A ;

/* Proposition 18 */ 2 return (1/r) × NFpoly diff(u f, A) Algorithm 5 provides an alternative method for computing the normal form of a differential fraction, which takes advantage of the fact that one can start with a partial reduction step then proceed with Algorithm 2.

Conclusion

In this paper, we have provided new and conceptually simple algorithms for computing normal forms of fractions modulo regular chains, in both the nondifferential and the differential case. We have also provided a new proof of the membership test to ideals defined by regular chains, which highlights the main non trivial argument underlying this well known result. However, one cannot exclude the existence of a completely elementary proof, which would actually lead to a simplified general theory of regular chains. input : a differential fraction f /g, and A a differential regular chain with g regular in R/A output: the normal form of f /g modulo A 1 compute h f f = f (mod [A]) where h f denotes a power product of separants of A ; /* f is a partial remainder of f by A */ 2 compute h g g = ḡ (mod [A]) where h g denotes a power product of separants of A ; /* ḡ is a partial remainder of g by A */

3 return NF h g f h f ḡ , A

2. 2 . 2 .

 22 Caveats. The simplicity of the new algorithms hide subtleties, detailed below. The following proposition is [5, Proposition 23].

6 F.

 6 Boulier, F. Lemaire, M. Moreno Maza and A. Poteaux Math. in Comp. Sc.

F

  . Boulier, F. Lemaire, M. Moreno Maza and A. Poteaux Math. in Comp. Sc.

  next definition and proposition are [4, Definition 4.1, Proposition 4.3 and Proposition 5.1].

Algorithm 3 : 6 if θ = 1 then 7 13 return ( 1 /Algorithm 4 :

 3671314 NFpoly diff (f, A) input : a differential polynomial f and A a differential regular chain output: the normal form of f modulo A 1 if f is fully reduced w.r.t. derivatives of the leading derivatives of the elements of A which actually occur in f , let w be the highest, w.r.t. the ranking ;5 let p k ∈ A, with leading derivative v k be such that w = θv k for some θ ∈ Θ ; compute i α k f = f + q p k ; /* f = prem(f, p k , v k ) */ 8 compute u/r an inverse of i k modulo A ; /* Proposition 4 */ 9 else 10 compute s α k f = f + q (θp k ) ; /* f = prem(f, θp k , θv k ) */11 compute u/r an inverse of s k modulo A ; /* Proposition 4 */ 12 end r α ) × NFpoly diff u α f , A 14 end NF diff (f /g, A)

Algorithm 5 :

 5 NF diff v2 (f /g, A)

return (1/r α ) × NFpoly(u α f , {p 1 , . . . , p n-1 });

This work has been supported by the bilateral project ANR-17-CE40-0036 and DFG-391322026 SYMBIONT..