N Prioux
email: nancy.prioux@ensiacet.fr

R Ouaret

G Hetreux

J.-P Belaud

ENVIRONMENTAL ASSESSMENT COUPLED WITH MACHINE LEARNING FOR CIRCULAR ECONOMY

Keywords: Agricultural Waste, Life Cycle Assessment, Artificial Intelligence, Data Science, Clustering, Biorefinery 1

The circular economy and its various recirculation loops have become a major study subject over recent years, particularly in the field of agriculture, which is a significant source of waste production. There have been several studies focused on transforming agricultural lignocellulosic waste with "sustainable" processes: economically viable, socially accepted, and environmentally friendly. Thanks to "life cycle thinking", it is possible to assess these potential environmental impacts. However, these environmental analyses generally require a massive volume of specific data, the collection of which can be time-consuming and tedious, or impossible to practice. On the other hand, scientific articles describing the processes for the valorization of agricultural by-products are intriguing but rarely exploited sources of data. In this paper, a hybridization of data science techniques and environmental analysis proposed to improve Life Cycle Analysis (LCA) thanks to Machine Learning (ML). ML part of the proposed approach is based on unsupervised learning, which is composed by two methods: dimension reduction using the Multidimensional Scaling and clustering technique using k-means. Composed by five steps and dedicated to researchers or R&D engineers, the approach is oriented towards offering decision on technologies and processes for waste to energy in early eco-design step. The case study in the domain of pre-treatment processes for corn stover and rice straw is detailed. The results show that all impacts that concern chemical pollution of soil and water are found in the same cluster. Other impacts are detected in the same cluster which are related to the land use and the land transformation. In the same vein, two purely mechanical pre-treatments has been identified and grouped by Multidimensional Scaling and k-means.

INTRODUCTION AND SCIENTIFIC CONTEXT

The idea of sustainable development is a complex problem applying in almost all areas of life, from social and health to economic development and environmental assessment [START_REF] Wced | World commission on environment and development[END_REF]. In recent years, policymaking increasingly involves action plans for climate, sustainable development or environment like Kyoto Protocol or Paris climate agreement.

In December 2019 at the COP25 Climate Summit in Madrid Spain (COP25), the European Green Deal is presented by the European Commission and is accepted in January 2020 (European Commission 2019). The European Green Deal plan includes several actions aimed at stopping climate change and reducing the level of pollution emitted into the atmosphere. One of the main action plan named "Circular Economy Action Plan" proposes to use the Circular Economy (CE) model to simulate the use of sustainable models. This plan completes the 2015's first Circular Economy Action Plan whose all 54 actions under the plan have been delivered or are being implemented. Ellen Macarthur Foundation (2015) defines the circular economy as "be one that is restorative and regenerative by design and aims to keep products, components, and materials at their highest utility and value at all times, distinguishing between technical and biological cycles". This new concept of sustainable development could reconcile some environmental, economic, and social aspects. A review paper proposed by [START_REF] Ghisellini | A review on circular economy: the expected transition to a balanced interplay of environmental and economic systems[END_REF] describes the origins, principles and the limitations of CE models. Following the first plan in 2015, the French government recommended the SNTEDD (National Ecological Transition Strategy for Sustainable Development) which consists of nine areas, including the CE. According to the French Agency for Environment and Energy Management (ADEME), the CE takes into account three areas of action: (1) consumption through consumer demand and behaviour, (2) supply and economic actors for whom industrial ecology is an accepted and promising path from the initial design of a territorial area and (3) waste management policy (Belaud et al. 2019a). These three areas represent the entire life cycle of a product, service, or a process. To achieve sustainable models, life cycle thinking can help improve environmental performance while maximizing economic and social benefits. Several global methods have emerged to design waste recycling processes that fit into the circular economy [START_REF] Grimaud | Decision Support Methodology for Designing Sustainable Recycling Process Based on ETV Standards[END_REF]. Agriculture is a particular area in which CE and life-cycle thinking have developed over the last decades. [START_REF] Singh | Thematic exploration of sectoral and cross-cutting challenges to circular economy implementation[END_REF] highlights the connections between CE and various major themes and notably the links between CE and LCA, biomass, biorefinery, bioenergy and waste management.

In France, the agri-food industry produces about 2.6 million tons of organic waste per year, a figure that is constantly increasing [START_REF] Barry | La production de déchets non dangereux dans les industries agroalimentaires[END_REF]. This trend is accompanied by the projected increase in the world's population. At the same time, human activities are reducing the amount of land available for agriculture, which inevitably has an impact on farming systems. For some people, new agricultural technologies that facilitate sustainable intensified agriculture seem to be the best solution for the future [START_REF] Garnett | Sustainable Intensification in Agriculture: Premises and Policies[END_REF]. However, such intensification in agriculture could result in more waste of products and resources [START_REF] West | Leverage points for improving global food security and the environment[END_REF]. According to [START_REF] Horton | An integrated theoretical framework to enhance resource efficiency, sustainability and human health in agri-food systems[END_REF], a major challenge in achieving sustainability in agriculture is categorizing wastes. They distinguished them into two categories: wastes generated by inputs, such as fertilizers remains or used water, and wastes generated from following treatment processes. The latter, composed of huge amounts of lignocellulosic by-products, comes mainly from incomplete conversion of biomass or raw materials processing in the supply chain. Lignocellulosic biomass is one of the most abundant and least expensive renewable resources on Earth. The production of biomaterials, biomolecules and bioenergy often relies on the bioconversion of it, where it is enzymatically hydrolysed to produce glucose. It is also possible to produce different biomaterials to replace plastic material, such as composite beads based on olive pomace [START_REF] Lissaneddine | Performance and dynamic modeling of a continuously operated pomace olive packed bed for olive mill wastewater treatment and phenol recovery[END_REF]. Lignocellulosic biomass consists of four main elements: lignin, cellulose, hemicellulose, and phenolic acids. Only cellulose and hemicellulose can be hydrolysed to generate glucose. Although lignocellulosic biomass is a renewable resource, transformation processes must be sustainable to participate in sustainable development. For this reason, more and more agri-food processes have been incorporated with various sustainability assessments (Food SCP Round Table European Commission 2012;[START_REF] Raymond | Spatially and Temporally Explicit Life Cycle Environmental Impacts of Soybean Production in the U.S. Midwest[END_REF]. In order to generate a good yield of glucose, it is important to pre-treat the biomass before hydrolyzing it. Over the past 30 years, numerous pre-treatment processes have been studied and published [START_REF] Davis | Secondary Resources in the Bio-Based Economy: A Computer Assisted Survey of Value Pathways in Academic Literature[END_REF]. Various factors were applied to evaluate and compare the performance, efficiency, or environmental influences of these processes [START_REF] Joglekar | Environmental Impact Study of Bagasse Valorization Routes[END_REF]. Environmental factors, energy consumption and energy efficiency were considered as 3 classical factors [START_REF] Zhu | Pretreatment of woody biomass for biofuel production: energy efficiency, technologies, and recalcitrance[END_REF][START_REF] Barakat | Dry fractionation process as an important step in current and future lignocellulose biorefineries: A review[END_REF]. However, there is a lack of criteria to guide the choice between all these processes. The use of environmental, economic, and social assessment in a CE context is a good way to guide this choice. In this paper, only the environmental dimension is studied, even though it is entirely possible to complete our approach through economic and social manners. Our choice for the environmental analysis is based on the Life Cycle Assessment (LCA) method published by the International Organization for Standardization (ISO 14040:2006(ISO 14040: 2006)).

From the literature on LCA, it appears that the most commonly used approaches rely on knowledge-based approaches and multi-criteria decision-making. The Life Cycle Impact Assessment (LCIA) knowledge acquisition is based on several experts from different backgrounds (ecologist, chemist, lab technician, epidemiologist, ... ). Most previous LCA studies and methods for assessing sustainability are based on environmental, social, and economic dimensions (Svensson et al., 2018;[START_REF] Zhao | Data envelopment analysis for sustainability evaluation in China: Tackling the economic, environmental, and social dimensions[END_REF]. For measuring country sustainability performance, [START_REF] Tan | An adaptive neuro-fuzzy inference system (ANFIS) approach for measuring country sustainability performance[END_REF] used an adaptive neuro-fuzzy inference system (ANFIS) approach. The proposed approach is effective to measure the countries' sustainability performance, however, the results of the ANFIS method are strongly dependent on appropriate training data. In this regard, we have recently seen an increase in the fuzzy-set theory in the sustainability assessment field. Many researchers have applied this approach for country sustainability assessment. The advantage of the fuzzy-set and fuzzy inference system is it can emulate the behaviour of skilled humans and handle the multidimensional complexity of vague situations. Several studies have been conducted using other Machine Learning (ML) algorithms. For example, [START_REF] Kouchaki-Penchah | Environmental management of tea production using joint of life cycle assessment and data envelopment analysis approaches[END_REF] and [START_REF] Zhao | Data envelopment analysis for sustainability evaluation in China: Tackling the economic, environmental, and social dimensions[END_REF] have combined Data Envelope Analysis (DEA) with LCA. Such hybrid systems have been tested for determining the energy efficiency, for measuring the efficiency of a sustainable development system by tracking economic, environmental, and social dimensions.

Some studies have examined the spatiotemporal patterns of Life-Cycle Environmental (LCE) based on multi-year assessment over a large geographic region, and improving the spatial resolution of LCE [START_REF] Lee | Spatially and temporally explicit life cycle global warming, eutrophication, and acidification impacts from corn production in the U.S. Midwest[END_REF]). In the same context, Romeiko et al. (2020) utilized boosted regression tree models to identify the top influential factors among soil, climate, and farming practices, which drive the spatial and temporal heterogeneity of life cycle environmental impacts. The results of this study showed that soil organic content and nitrogen application rate were the top influencing factors for life cycle eutrophication (EU) and acidification (AD). The top influencing factors for life cycle global warming (GWP) impacts were soil texture, nitrogen application rates, and temperature in March. By using Boosted Regression Tree (BRT) model with Gaussian distribution, the same team researchers have conducted another study based on the prediction of LCE impacts of corn production under climate variability and scenarios. To the best of our knowledge, many methods for LCIA or for general sustainability assessment are "supervised", meaning that learning algorithms are performed from inputoutput data [START_REF] Abdella | Sustainability assessment and modeling based on supervised machine learning techniques: The case for food consumption[END_REF]). However, in many cases, one needs to extract latent data structures based on the observations and the mixture nature of the endpoint area protection. For this reason, our research relies on the use of Dimension Reduction (DR) techniques for the assessment of sustainability through a set of impacts (Climate change, Human Health, Ozone depletion, Human toxicity, ... ) and underlying process. Several strategies have been explored to achieve this goal using "Unsupervised" learning which includes dimension reduction and clustering.

The purpose of this paper is to present an approach to aid in the analysis and comparison of different pre-treatment processes for lignocellulosic biomass and different types of biomass in the context of the circular economy. The main contribution of this study is the development of an approach that suggests the enhancement of the traditional LCA method by coupling methods from ML. This enhancement can be found in two main items: (i) the use of experimental opendata extracted from scientific web engines and (ii) the inclusion of unsupervised ML for the interpretation and analysis of environmental impacts. To date, few research studies have been conducted using big data technologies for LCA. In particular, there are studies to complete the background data, necessary for LCA, or to adapt LCA to technological developments [START_REF] Cooper | Big Data in Life Cycle Assessment: Big Data in Life Cycle Assessment[END_REF][START_REF] Bhinge | Data-intensive Life Cycle Assessment (DILCA) for Deteriorating Products[END_REF]. However, in our knowledge, no study uses experimental data from public web, such as data from scientific article or from scientific open database, to complete the foreground database required in a LCA. This point is developed in the step 2 of our approach, and to this end, the proposed approach constitutes a contribution to the classical literature (section 2.2). The use of unsupervised learning for LCA is motivated to detect the environmental fingerprint of lignocellulosic biomass processes. More precisely, this study lies in the characterization of different environmental impacts into clusters to help analysis, interpretation and decision. This point corresponds to the step 4 of our approach (section 2.2 and result section).

This approach is intended for researchers or research and development engineers who would like to make the first screening in pre-treatment processes or biomasses using data contained in the scientific literature. It can be deployed during a preliminary study phase and be used to assist an initial decision for the development of a laboratory or semi-industrial pilot. After an introduction to the description of data availability, section 2 outlines this general approach. Section 3 will use this approach to deal with the study of different pre-treatments processes of rice straw and corn stover and then provides the results and interpretations.

METHODOLOGICAL APPROACH

LCA and ML methods

Data science can be used at different levels of sustainability management. One of the challenges in valorising byproducts in the agricultural supply chain is to design a process that is as sustainable as possible. The supply chain includes several operational steps, from biomass selection to waste disposal, and goes through various processing steps. Each step in the chain can be described with its inputs and outputs, as well as its energy and economic data. All these data are required so that the LCA could be carried out taking into account the data diversity and its heterogeneous sources. Various types of data from heterogeneous sources are required for environmental analysis. ML would provide appropriate methods and technologies for this analysis. The main objective of this approach is to analyse diverse technological processes and to provide decision supports by analysing the results.

LCA method

LCA is a tool for assessing the potential environmental impact of a product or service through its entire life cycle. The life cycle of a product/service can be broken down into several stages which starts from the designing of the product to various stages of processing and utilization, and ends with the disposal or recycling stage of the product. LCA involves four steps: defining the objective and scope of the study, conducting a life cycle inventory, conducting a LCIA, and, at last, analysing and interpreting the results.

The first step is to define the objective, system boundaries, and the functional unit of the life cycle. The functional unit must be well defined so that it will not disturb the outcome of the LCA [START_REF] Burgess | Application of life cycle assessment to chemical processes[END_REF].

The Life Cycle Inventory (LCI), which is the second step, aims at collecting data about the quantities of pollutants emitted and resources extracted throughout the life cycle. This inventory concerns two kinds of data: the background data and the foreground data. The foreground data represents data on processes of interest to decision-makers, while the background data includes all the other processes data that is related to processes of interest [START_REF] Clift | Toward a Coherent Approach to Life Cycle Inventory Analysis[END_REF][START_REF] Elghali | Developing a sustainability framework for the assessment of bioenergy systems[END_REF]. Data for these inventories can be obtained directly, through on-site measurements (primary data), or indirectly, from published scientific papers, models, and databases (secondary data). The foreground system is generally based on primary data, while the background system is based on secondary data sources (Guinée 2002).

The third step is the LCIA, during which numerous data of pollutants and resources collected at the LCI are calculated through specific methods to evaluate their environmental impacts [START_REF] Suh | Methods for Life Cycle Inventory of a product[END_REF].

The last step is the interpretation of the results, which involves identifying important issues based on the results of the LCI and LCIA, assessing their sensitivities, checking for their consistency and completeness, and formulating a report with conclusions, recommendations, as well as limitations. This last stage of LCA is "delicate" for novices and, sometimes, even experts. It may be, therefore, beneficial to enrich the original LCA method with other techniques to assist researchers or engineers during their analysis and interpretation.

In our approach, we use only secondary data. Foreground data are taken from the public web through scientific articles (Web of Science and Science Direct) to constitute an unstructured database, whereas background data are available in structured databases (public or private) such as EcoInvent [START_REF] Frischknecht | The ecoinvent Database: Overview and Methodological Framework (7 pp)[END_REF]). To analyse this data and structured them, engineering and ML techniques for clustering were used.

Multidimensional Scaling method

In this article, we use Multidimensional Scaling (MDS), also called Principal Coordinates Analysis (PCoA) which is a method for visualizing (de)similarity between objects in a reduced dimensional space. It is designed to understand the proximity and opposition structures. Starting from information about the mutual similarities of 𝑛 objects, often with the similarity matrix where ∆ = 𝛿 (𝑖, 𝑗) 0≤𝑖≤𝑛;0≤𝑗≤𝑛 , we look for a configuration of 𝑛 points (𝑅 2 ) which would be like that if the object 𝑖 resembles to 𝑗 more than the object 𝑙 to 𝑘, we have 𝛿 (𝑖, 𝑗) < 𝛿 (𝑙, 𝑘).

For any distance matrix of size 𝑛 × 𝑛, the MDS allows us to find a set of 𝑛 points marked by their coordinates whose similarity matrix is equal or very close according to the data. Let 𝑥 𝑟 (𝑟 = 1, . . . , 𝑛) be the coordinates of 𝑛 points in a 𝑝 dimensional Euclidean space where 𝑥 𝑟 = (𝑥 𝑟1 , 𝑥 𝑟2 , . . . , 𝑥 𝑟𝑝 ) 𝛿 𝑟𝑠 2 , which is deducted from the decomposition 𝑥 𝑟 ⊤ 𝑥 𝑠 , the matrix 𝐵 is obtained using the following relation :

𝐵 = 𝐻𝐴𝐻 (1)
where 𝐻 is the centring matrix: 𝐻 = 𝐼 -𝑛 -1 𝐼 • 𝐼 ⊤ , where 𝐼 , a 𝑛 × 𝑛 identity matrix and𝐼 = (1,1, . . . ,1) ⊤ , a vector of 𝑛 ones. The elements of 𝑏 𝑟𝑠 = 𝑎 𝑟𝑠 -𝑎 𝑟• -𝑎 •𝑠 + 𝑎 •• , where ::

𝑎 𝑟• = 𝑛 -1 ∑ 𝑠 𝑎 𝑟𝑠 , 𝑎 •𝑠 = 𝑛 -1 ∑ 𝑟 𝑎 𝑟𝑠 , 𝑎𝑛𝑑 𝑎 •• = 𝑛 -2 ∑ 𝑟 ∑ 𝑠 𝑎 𝑟𝑠 .
The algorithmic procedure of MDS can be summarized, as illustrated in [START_REF] Cox | Multidimensional scaling[END_REF]) by the following steps:

1.

Obtain the matrix ∆ = 𝛿 𝑟𝑠 of dissimilarities.

2.

Find the matrix 𝐴 = [-

1 2 𝛿 𝑟𝑠 2 ].
3.

Find the matrix 𝐵 = [𝑎 𝑟𝑠 -𝑎 𝑟• -𝑎 •𝑠 + 𝑎 •• ].

4.

Find the eigenvalues 𝜆 1 , 𝜆 2 , . . . , 𝜆 𝑛-1 and the eigenvectors 𝑣 1 , 𝑣 2 , . . . , 𝑣 𝑛-1 . If 𝐵 is semi-defined positive (some eigenvalues are negative), either (𝑖) ignore the negative values and continue, or (𝑖𝑖) add an appropriate constant 𝑐 to the (dis)similarities

5.

Choose an appropriate size number 𝑙, possibly using

∑ 𝑙 1 𝜆 𝑖 ∑ (𝑝𝑜𝑠𝑖𝑡𝑖𝑣𝑒 𝑒𝑖𝑔𝑒𝑛𝑣𝑎𝑙𝑢𝑒𝑠)
.

6.

The coordinates of the n points in Euclidean dimension space 𝑙 are given by 𝑥 𝑟𝑖 = 𝑣 𝑖𝑟 (𝑟 = 1, . . . , 𝑛; 𝑖 = 1, . . . , 𝑙).

Approach

Using a general architecture based on sustainability principles (Belaud et al. 2019b), our approach is composed of five steps (Figure 1): (1) goal and boundaries, (2) data architecture, (3) the environmental assessment, (4) results visualization and analysis, and (5) decision. Each step has sub-steps and the transition from one to the other can be done in either direction through feedback loops. It is recommended to iterate to consolidate the results and, subsequently, the choices resulted from the interpretations.

Goal and scope

In this step, the goal of the study case and the boundaries of the system must be clearly defined. It is possible to come back to this step if a blur or a problem is noticed in the following steps. Life cycle thinking is recommended. This thinking encourages a "cradle to grave" or "cradle to gate" approach. In the CE model, the part of the life cycle in which the product is used is a key element in moving towards the ecological transition. "Cradle-to-Gate" approaches are often preferred because integrating downstream elements into sustainability analyses can be tedious and difficult. Especially, scientists and engineers often find it impossible to consider end-user or consumer behaviour in their models. The system boundary has a significant effect on subsequent evaluations. For example, it needs to be clarified whether the upstream biomass supply chain is considered. Once the objective and scope have been properly defined, the supply chain, technologies and transformation processes should also be described. This description should be as complete as possible concerning the operations of the process, the location of the study, the various inputs and outputs, and the type of energy used. These details ensure the relevancy of the collected data.

Data architecture

The data architecture is directly inspired by the construction of massive data architecture and consists of five substeps: (i) data collection and extraction, (ii) data enrichment and storage, (iii) data processing, (iv) (raw) data analysis, and (v) (raw) data visualization.

Collecting and extracting data from unstructured databases requires special tools, such as data queries (SQL queries) or Online Analytical Processing (OLAP). This sub-step is more complicated for unstructured data collected from the web. The meta-data associated with these web pages can be used for their classification and to provide access to their content. For example, the definition of the MARC (MAchine-Readable Cataloging) format in the early 1960s standardized metadata in documentary resources. With the development of Linked Web, especially the Resource Description Framework (RDF, a W3C standard), it is possible to use SPARQL Protocol and RDF Query Language to request the RDF. The extraction process generates a structured table which differs according to the format of the web pages: API, HTML, or pdf. The extraction process can be automatic, semi-automatic or manual. Data can be extracted automatically from web pages with an API format. However, extracting relevant data from scientific articles requires a reading guide, which is generated from ontologies and pre-processing analysis of experts. This type of extraction is therefore semi-automatic.

For data enrichment and storage, the extracted data is stored in Relational Database Management Systems (RDBMS). Enrichment is the process of adding data to the DBMS. This data comes from experts or pre-processing models. The models are empiric numerical simulations of unit operation type, thermodynamics, and energy for the control of flows, transformations, and transfers.

Data processing involves cleaning, adding, and deleting data for volume and value management. After this step, a second, more accurate and more accessible database can be generated. Keeping this second database saves time for the following sub-steps and avoids misinterpretations during the data analysis sub-step and the environmental analysis step. The data analysis sub-step depends on the objective, data domains and decision-makers. Different types of analysis are possible: descriptive analysis (what happened?), diagnostic analysis (why did it happen?), predictive analysis (what is going to happen?) and prescriptive analysis (how can we make it happen?). Some of these methods are visual, so data visualization can be included in the data analysis sub-step. Data visualization can also be achieved by plotting the raw data in the form of a simple or interactive graph.

2.2.3.

Environmental assessment In this step, the LCI and the LCIA is performed. It is very important to follow the steps recommended by ISO (ISO 14040:2006(ISO 14040: 2006) ) and to use the data from the previous step for the process data. The background data is taken directly from the EcoInvent database. Once the inventory (LCI) is completed, one (or more) impact calculation method(s) must be chosen under Step 1 (The purpose and the boundary of the study). It is possible to choose several impact calculation methods to compare the results, but this may hinder the final analysis and the choice of process or biomass for the researcher in his or her first analysis.

Results visualization and analysis

It is in this stage where we bring the ML method into use to help analyse environmental impacts which is the focus of this paper. At the end of the previous step, the result obtained is a matrix with the processes, biomasses, and environmental impacts, which is difficult to analyse for a non-expert in LCA. From the statistical literature, this step combines traditional techniques for dimension reduction and unsupervised clustering to extract knowledge about life cycle impact assessment (LCIA). More specifically, the hybrid approach is based on Multidimensional Scaling (MDS) using Canberra distance [START_REF] Lance | Computer programs for hierarchical polythetic classification ("similarity analyses[END_REF] and 𝑘-means. MDS is an algorithm that transforms a distance matrix into a set of coordinates such that the distances derived from these coordinates approximate as well as possible the original distances. In this work, we used Canberra distance (a weighted version of the Manhattan distance) for constructing the proximity matrix and then the data is mapped on a lower-dimensional (two or dimensions) spatial representation. These methods reduce the dimensional space (the variable set) while preserving the maximum amount of information. The goal is to seek "hidden" structures in the multidimensional data and to help to interpret the grouped endpoint area of the LCIA assessment matrix. The benefit of such approach is that the data-driven methods require minimal process knowledge to perform this task.

Figure 2 outlines the data-driven processing for LCIA. In the proposed method, first, the process-impacts matrix is used as input (similarity matrix) for learning from data (Part A). Second, DR techniques project the raw process data to a lowerdimensional space (2 or 3) (Part B). There are several statistical methods to achieve this mapping such that uninformative variance in the data is discarded [START_REF] Burges | Dimension Reduction: A Guided Tour[END_REF]. After projection raw data by a DR technique, the clustering approach is then applied to consider similar life cycle impacts as well as in the process within the lower-dimensional space (Part C). Finally, the user (expert) analyses the grouped data in clusters to relate them to meaningful process/impacts (Part D). This last cluster assignment and extracting information step is called knowledge discovery. An advantage of the workflow is that it is relatively simple to use because each DR and clustering combination requires the specification of only one or two parameters and techniques. The idea behind the proposed workflow is motivated by the fact that sustainability assessment requires both qualitative and quantitative criteria.

Decision making

The visualization of data clusters from the previous step can help the researcher with his decision-making. This decision can be made by the researcher himself or by a group composed of different engineers/researchers from different fields.

RESULTS AND DISCUSSION

In this section, the first results of our approach are presented with a case study comparing different pretreatment processes of two biomasses: corn stalk and rice straw. A software programmed on Excel supports the data architecture and the environmental steps as well as on R for results analysis and visualization step. The software part is shown in Figure 3. The internal software has been verified with:

-well-establisted ProSim for the process simulation part. The process simulation is used to raw data processing (step 2 iii), especially for the verification of mass balances. An article whose material balance could not be verified is deleted from the database. -SimaPro for the environmental impact simulation part (step 3) For the origins of data, foreground data come from data collection of our approach -they are called process data -and background data come from. EcoInvent Database. This study concerns mainly the foreground data, the background data are not modified or analysed.

Study case first results

The purpose of the study is to help the researcher to select a process and/or a biomass for the production of glucose. The study has a cradle-to-grave approach i.e. the boundaries of the study range from biomass to the enzymatic hydrolysis. Biomass is here considered as a waste from agriculture that has no impact -the impacts are attributed to the final product of agriculture (corn and rice). Besides, the biorefinery is considered to be relatively close to the site, and therefore the transport stage is negligible. The function of the system is, therefore "producing glucose" and the functional unit is the "production of 1 kg of glucose". All results can therefore be expressed in terms of the functional unit. Figure 4 [START_REF] Liu | Alkaline twin-screw extrusion pretreatment for fermentable sugar production[END_REF] Now, the step 2 of the approach -data architecture -was developed. A more detailed step can be found in [START_REF] Belaud | Intensive Data and Knowledge-Driven Approach for Sustainability Analysis: Application to Lignocellulosic Waste Valorization Processes[END_REF]. The selection of articles describing the processes was made by process experts on corn stove and rice straw thanks to keywords: "rice straw", "corn stover", "treatment","hydrolyze" and "milling" in scientific databases, such as Web of Science or Science Direct. Eighteen papers were initially selected and stored by topics (different type of process). The data from the papers were then extracted and stored in CSV files before being passed into our software (cf. Figure 3). This software developed on Excel did a first cleaning of the data using process simulations to calculate and verify the mass balance. Thirteen articles were selected because they had either too much missing data for verification or too many inconsistencies. The raw data analysis and visualisation are not used in this case study and they will be the subject of future research.

Then the environmental assessment is performed on the remaining process data. The LCA method used here is ReCiPe 2016 [START_REF] Huijbregts | ReCiPe2016: a harmonised life cycle impact assessment method at midpoint and endpoint level[END_REF] , the background database is EcoInvent and the processed data comes from the previous steps. The result of this step is a matrix of 18 mid-point impacts and 13 processes. This matrix is then analyzed by multidimensional scaling.

The interpretation of an MDS result is the same as for any other dimension reduction; objects that are closer together on the scatter plot are more alike than those further apart. That is the projected points are arranged in such a way that the grouped ones (small geometrical distance between them) will reflect original relationships in the data. However, additional information to make the projection more informative. As illustrated in Figure-workflow 1, a clustering algorithm has been applied to the MDS projection to highlight the most similar objects (Impacts -Process). The abbreviation of the impacts is introduced to facilitate the visualization (Table 1, Annex 1). The two-dimensional of MDS results of projected impacts (17 impacts) is shown in Figure 5. It presents 4 sub-figures of the first 4 dimensions with the most significant combinations (dimension 1 vs. dimension 2, ..., dimension 2 vs. dimension 3). For example, in the first figure (top left), we have represented the projection of the 17 impacts on the first two dimensions, which represent a total variance of 45%. The percentage of explained variance for the first four components is 70%. The visualisation of the four dimensions shows the same three groups and we can clearly distinguish three clusters using k-means:

-Group 1: Terrestrial Acidification (TA), Freshwater Eutrophication (FrEu), Terrestrial Ecotoxicity (TecoX) and Freshwater Ecotoxicity (FrEco). TA is most closely associated with FrEu, with respondents considering them almost identical. Other points of Impacts in this cluster are considered similar based on their proximity such TecoX with TA. The FrEco item is furthest in group 1. Almost all impacts that concern chemical pollution of soil and water are found in this group. The exemption is the marine ecotoxicity found in group 2.

-Group 2: Human Toxicity (HT), Particulate Matter formation (PM), Climate Change to Human Health (CCHH), Marine Ecotoxicity (MaEco), Metal Depletion (MeDe), Fossil Depletion (FossDe), Climate Change to Ecosystems (CCE) and Ionizing Radiation (IR). This group forms three sub-clusters with superposed points (from the 2-D perspective). This suggests that these points are highly similar based on the Canberra distance. For example, HT and PM, CCE and IR and, MaEco, MeDe and FossDe. Here, we find a group quite heterogeneous where impacts not presented in groups 1 and 3 are found. The marine ecotoxicity expected rather in group 1 is found in this group.

-Group 3: Urban Land Occupation (UrbLOcc), Photo-Chemical Oxidant formation (Pohto_ChOx), Agricultural Land Occupation (AgLOcc), Ozone Depletion (OD) and Natural Land Transformation (NLTran). This group mainly includes impacts related to land use and land transformation. The results of MDS using process distance matrix is presented in the Figure 6. In this case, the percentage of explained variance for the first four components is 98%, which is an excellent representation in lower-dimensional spaces. That is from a matrix with dimension 15, we lose only 2% to represent the first 4 dimensions. The results are quite similar to those obtained by using MDS on the impacts matrix (Figure 5): three distinguished groups have been identified. The acronyms captioning the dots represent each type of process. These will not be explained here. Three distinct groups of processes can be identified. Very tight and separate clusters appear in the process data, which may suggest that each cluster is a domain or sub-domain that needs to be analysed individually. In group 3, for example, there are two purely mechanical pretreatments (PM-UFM for "pre-milling and ultra-fine milling"). Going back to the impacts, we find that these two pretreatments have a very significant impact on the depletion of fossils compared to the others. For Group 1, the three pretreatments have relatively similar impacts. For all impacts and fossil depletion, the impact costs around $10, whereas group 2 pre-treatments have an impact costs around $1. Through these groups, we can find an analysis that could be done by an LCA expert, whereas the visualization was done using tools from data science field.

Discussion and perspectives

The approach proposed in this article is, to our knowledge, unique. It is intended to be generic and practical for the researcher or the R&D engineer. For instance, before starting a project involving the recovery of energy and valuables from wastes, the researcher will be able to compare different processes existing in the literature. The enhancement of the traditional LCA method by coupling tools from (big) data science and algorithms from artificial intelligence allows to have different discussion on environmental impacts. Indeed, on the one side, tools from data science allow to extract and collect data directly from scientific papers. And on the other, the MDS allows to simplify the discussion on environmental impacts. In this paper and this case study, the MDS method was studied because it was the one that gave the best percentage of explained variance.

Limits of this approach

This case study has shown various limitations of the data, and especially data from scientific databases generated by laboratories. The life cycle analysis is therefore carried out at a research TRL scale (TRL scale 1/2) which can lead to a change of scale if we want to switch to an industrial pilot [START_REF] Bianco | Chapter 4 -Bioreactors for wastewater to energy conversion: from pilot to full scale experiences[END_REF]. Moreover, as the processes are not continuous in the labs, the data had to be processed so that the inputs of a process step correspond to the outputs of the previous step.

The data processing can be improved, for example, adding scale changes made in the process engineering. Another limitation in the data concerning the abundance and quality, which can be very low for some innovative and new processes such as biomass pre-treatment processes. It would be interesting for future studies to change the functional unit and broaden the boundaries of the study to take into account the overall logistical scale or to consider effluent recycling. Indeed, all impacts have been related to the final product (glucose) and no recycling of effluents is taken into account. Recycling could reduce the impacts of some processes more than others. In fact, this approach serves more like aid in deciding on which process or biomass should not be used or which routes could be taken, instead of deciding which one is the best process.

A final limitation comes from the main input, which is biomass, a waste product from agriculture which can be of variable quality depending on time and storage. Furthermore, it is necessary to evaluate its durability and supply, which could also vary over time. These three criteria might be a problem for both the economic and environmental process.

Perspectives of the work

Further improvements to the proposed approach could be obtained with any of the following: (i) including additional specific data sources, methods and visualisations for the economic and social areas. These aspects could enhance the sustainability data inventories and assessment methods. We can cite Social LCA, Life Cycle Cost (LCC), consequential LCA [START_REF] Curran | The international workshop on electricity data for life cycle inventories[END_REF] or dynamic LCA [START_REF] Collinge | Dynamic life cycle assessment: framework and application to an institutional building[END_REF].

(ii) addition of new circular indicators like service-oriented optimization model (Al-Aomar and Alshraideh, 2019) (iii) progress towards the automation in the data extraction step (step 2 in our approach). This would make it possible to save time and to add new sources of data more easily.

(iv) test different ML techniques to analyse and visualise the raw data such as: a) Natural Language Processing (NLP) to analyse the text of the scientific literature b) Other unsupervised learning methods, in particular Dimensionality Reduction Techniques. (v) from our feedback with the Excel-VBA research tool, the development of a complete ergonomic computing framework remains to be achieved. This would encourage stakeholders to adopt this approach and would facilitate decisionmaking through the implementation of collaborative decision-making techniques, such as Delphi-SWOT;

(vi) the design of models for calculating energies, for assessing the impacts on the environment of activities linked to new energy systems;

(vii) the generalisation of this principle and the development of a library of business and domain-specific models from agri-food process engineering. These models could be used to check and validate the data in the data architecture step. Controls could include, for example, an advanced material balance or energy analysis;

(viii) the development of data dispersion propagation and automatic qualitative explanation systems for stakeholders.

ANNEX 1: NOMENCLATURE OF THE ABBREVIATION

⊤

  and [𝐵] 𝑟𝑠 = 𝑏 𝑟𝑠 = 𝑥 𝑟 ⊤ 𝑥 𝑠 . From the Euclidean distance ∆ = 𝛿 𝑟𝑠 , and a matrix 𝐴 of elements [𝐴] 𝑟𝑠 = 𝑎 𝑟𝑠 = -1 2

Figure 1 :

 1 Figure 1: Five main steps of the approach

Figure 2 :

 2 Figure 2: Schematic of data driven processing for LCIA

Figure 3 :

 3 Figure 3: Software tools scheme

  shows a instance pre-treatment process from a Liu et al. (2013)'s paper.

Figure 4 :

 4 Figure 4: Pretreatment process example adapted from[START_REF] Liu | Alkaline twin-screw extrusion pretreatment for fermentable sugar production[END_REF]

Figure 5 :

 5 Figure 5: Scatter plot of MDS projection (two dimensions) and k-means clustering based on Impacts distance matrix. Percentage of explained information for the first four components is 70%.

Figure 6 :

 6 Figure 6: Scatter plot of MDS projection (two dimensions) and k-means clustering based on Process distance matrix. Percentage of explained information for the first four components is 98%.

CONCLUSION

A general approach coupling Multidimensional Scaling, k-means algorithm and environmental analysis was proposed. Composed of five steps, this approach is presented as a decision aid for the researcher in a pre-study. It is designed to save time and money by including no experiments and using public scientific data as a database. This approach has been tested in the example of the valorisation of lignocellulosic biomass into glucose through the comparison of pre-treatment processes and two biomasses: corn stalk and rice straw.

After structuring the data and life cycle analysis steps, the environmental impact-process matrix is analysed using an MDS method. A major result of this article is that it highlights the importance of using a hybridization of LCA and ML methods. In light of the results presented above, it can be concluded that:

• This study demonstrates the importance of data science methods to shed new light on LCIA;

• These results contribute significantly to the very small data set available in the literature on using unsupervised learning for LCA purposes.

DECLARATION -Funding

No funds, grants, or other support was received.

-

Conflicts of interest/Competing interests

The authors have no conflicts of interest to declare that are relevant to the content of this article.

-

Availability of data and code

The datasets and code generated during and/or analysed during the current study are available from the corresponding author on reasonable request.

-

Authors' contributions

All authors contributed to the study conception and design. Material preparation, data collection and analysis were performed by Nancy Prioux, Rachid Ouaret and Jean-Pierre Belaud. The first draft of the manuscript was written by N. Prioux and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript. The work is supervised by Gilles Hetreux and Jean-Pierre Belaud.


Marko Petkovšek

REDUCTION-BASED CREATIVE TELESCOPING FOR DEFINITE SUMMATION OF D-FINITE FUNCTIONS HADRIEN BROCHET AND BRUNO SALVY Dedicated to the memory of

Creative telescoping is an algorithmic method initiated by Zeilberger to compute definite sums by synthesizing summands that telescope, called certificates. We describe a creative telescoping algorithm that computes telescopers for definite sums of D-finite functions as well as the associated certificates in a compact form. The algorithm relies on a discrete analogue of the generalized Hermite reduction, or equivalently, a generalization of the Abramov-Petkovšek reduction. We provide a Maple implementation with good timings on a variety of examples.

Introduction

The algorithmic computation of definite sums originates in Zeilberger's algorithm in the 1990's [START_REF] Zeilberger | A fast algorithm for proving terminating hypergeometric identities[END_REF][START_REF] Zeilberger | The method of creative telescoping[END_REF][START_REF] Wilf | An algorithmic proof theory for hypergeometric (ordinary and "q") multisum/integral identities[END_REF]. Initially designed to deal with hypergeometric sums, his method of creative telescoping has been extended to differential settings [START_REF] Almkvist | The method of differentiating under the integral sign[END_REF][START_REF] Zeilberger | A holonomic systems approach to special functions identities[END_REF][START_REF] Takayama | Gröbner basis, integration and transcendental functions[END_REF][START_REF] Takayama | An algorithm of constructing the integral of a module -an infinite dimensional analog of Gröbner basis[END_REF] and next generalized to the large class of D-finite functions by Chyzak [START_REF] Chyzak | An extension of Zeilberger's fast algorithm to general holonomic functions[END_REF]. In order to compute a definite sum of F (t, x 1 , . . . , x m ) with respect to t, where each x i is a variable with respect to which one can apply a linear operator ∂ i (generally, differentiation or shift or q-shift operator), the creative telescoping algorithm constructs identities of the form [START_REF] Almkvist | The method of differentiating under the integral sign[END_REF] α c α (x 1 , . . . , x m )∂ α (F ) = G(t + 1, x 1 , . . . , x m ) -G(t, x 1 , . . . , x m ).

Here, the sum is over a finite number of multi-indices α and we use the multi-exponent notation

∂ α = ∂ α1 1 • • • ∂ αm m .
In the original version for hypergeometric summation, the monomials ∂ α (F ) are simply successive shifts F (t, n), F (t, n + 1), F (t, n + 2), . . . of a hypergeometric sequence F (t, n). Identities obtained that way can often be summed over t. The right-hand side telescopes by design. Since the coefficients c α do not depend on the variable t, the left-hand side results in an operator applied to the definite sum of F . From there, other algorithms can be applied to compute information on the sum. The operator in the left-hand side of Eq. ( 1) is called a telescoper of F and the function G in the right-hand side is the corresponding certificate. Chyzak's algorithm also deals with the differential analogue of Eq. [START_REF] Almkvist | The method of differentiating under the integral sign[END_REF] where the right-hand side is a derivative; it is used to compute information on definite integrals. Chyzak's algorithm, like Zeilberger's, looks for telescopers with an increasing number of monomials ∂ α with indeterminate coefficients c α and determines c α such that a certificate G exists in the vector space generated by the ∂ β (F ) for β ∈ N m+1 over the field of rational functions. The conditions of being D-finite is that this vector space has finite dimension, which allows for the existence of algorithms based on linear algebra. If no certificate exists, then the support is increased and the process is iterated. This stops either when sufficiently many operators have been found or when a prescribed bound on the orders is reached. (In the original hypergeometric case, no bound on the order is fixed a priori and termination is guaranteed for the family of proper hypergeometric terms [START_REF] Wilf | An algorithmic proof theory for hypergeometric (ordinary and "q") multisum/integral identities[END_REF].) Efficiency issues with this approach have led to the development of heuristics and a very useful Mathematica implementation by Koutschan [START_REF] Koutschan | A fast approach to creative telescoping[END_REF].

The most recent approach to deal with the efficiency issues with creative telescoping was initiated by Bostan, Chen, Chyzak and many co-authors who developed a class of reduction-based algorithms [START_REF] Bostan | Complexity of creative telescoping for bivariate rational functions[END_REF][START_REF] Bostan | Hermite reduction and creative telescoping for hyperexponential functions[END_REF][START_REF] Chen | A modified Abramov-Petkovšek reduction and creative telescoping for hypergeometric terms[END_REF][START_REF] Bostan | Efficient algorithms for mixed creative telescoping[END_REF][START_REF] Chen | Reduction-based creative telescoping for algebraic functions[END_REF][START_REF] Chen | Reduction-based creative telescoping for fuchsian D-finite functions[END_REF]. These algorithms avoid the computation of potentially large certificates. In the differential case, where the right-hand side of Eq. ( 1) is replaced by a derivative ∂ t (G), the principle is to use a variant of Hermite reduction to compute an additive decomposition of each monomial in the form

(2) ∂ α (F ) = R α (t, x 1 , . . . , x m )F + ∂ t (G α ),
where R α is a rational function with a certain minimality property. A telescoper is found by looking for a linear dependency between these rational functions for a family of monomials ∂ α . The computation of the rational function R α by Hermite reduction works by getting rid of multiple poles and isolating a polynomial part. This was first done for the integration of bivariate rational functions [START_REF] Bostan | Complexity of creative telescoping for bivariate rational functions[END_REF], of hyperexponential functions [START_REF] Bostan | Hermite reduction and creative telescoping for hyperexponential functions[END_REF] and of mixed hypergeometric-hyperexponential functions [START_REF] Bostan | Efficient algorithms for mixed creative telescoping[END_REF]. In these three cases, the vector space generated by the functions ∂ β (F ) for β ∈ N 2 has dimension only 1 over the rational functions.

For summation, the analogous problem for bivariate hypergeometric sequences was solved by replacing the Hermite reduction by a modified Abramov-Petkovšek reduction, thereby providing a faster variant of Zeilberger's algorithm [START_REF] Chen | A modified Abramov-Petkovšek reduction and creative telescoping for hypergeometric terms[END_REF]. For bivariate problems of dimension larger than 1, the method was extended to the integration of bivariate algebraic functions [START_REF] Chen | Reduction-based creative telescoping for algebraic functions[END_REF], of Fuchsian functions [START_REF] Chen | Reduction-based creative telescoping for fuchsian D-finite functions[END_REF] and more recently of Precursive sequences [START_REF] Chen | Reduction-based creative telescoping for p-recursive sequences via integral bases[END_REF][START_REF] Du | On the existence of telescopers for p-recursive sequences[END_REF] by means of suitable integral bases. An extension to the integration of purely differential bivariate D-finite functions in arbitrary dimension was first achieved by turning the differential equations satisfied by the function to be integrated into first-order differential systems; then, a variant of Hermite reduction can be designed at the level of vectors of rational functions [START_REF] Van | Constructing reductions for creative telescoping[END_REF]. This approach generalizes to purely differential D-finite functions in more variables [START_REF] Van Der Hoeven | Constructing reductions for creative telescoping: the general differentially finite case[END_REF].

Another method relies on cyclic vectors and allows the integration of arbitrary D-finite functions [START_REF] Bostan | Generalized Hermite reduction, creative telescoping and definite integration of D-finite functions[END_REF]. Without loss of generality, we assume that F is a cyclic vector for ∂ t , which means that all monomials ∂ α (F ) rewrite as M α (F ) with M α a linear operator in ∂ t only. (If F is not a cyclic vector, one finds a cyclic vector G, F = M F (G) for some linear operator M F in ∂ t only and the rest of the reasoning is unchanged.) Next, for any rational function u and any linear operator M in ∂ t , repeated integration by parts implies Lagrange's identity

uM (F ) -M * (u)F = ∂ t (P M (F, u)),
where M * is the adjoint of M and P M is linear in F, ∂ t (F ), . . . and u, ∂ t (u), . . . . Thus, a reduction-based algorithm is obtained by the additive decomposition of Eq. ( 2) with R α a solution of the generalized Hermite reduction

R α = M * α (1) mod Im(L * ).
The tools used in the reduction modulo the image of the linear differential operator L * are classical techniques used when looking for rational solutions of linear differential equations. This method of integration using generalized Hermite reduction extends to other contexts. This has been done in high generality in a preprint by van der Hoeven [START_REF] Van Der Hoeven | Creative telescoping using reductions[END_REF]. Our approach here is different: we focus on the case of summation only and give a simple self-contained presentation of the corresponding algorithm; our algorithm returns operators of minimal order 1 ; we make the choice to avoid algebraic extensions when possible; we present a Maple implementation that performs well in practice. Note that while in terms of complexity, minimal operators cannot be computed in polynomial time in general, in practice this does not seem to be an obstacle.

Example

The multiplication theorem for Bessel functions of the first kind J ν states that [27, 10.23.1] (3)

J ν (λz) = λ ν ∞ n=0 (-1) n (λ 2 -1) n (z/2) n n! J ν+n (z).
This can be proved automatically by showing that the left-hand side and the right-hand side satisfy the same set of mixed differential-difference equations with sufficiently many identical initial conditions.

1 Remark 5.6 in [START_REF] Van Der Hoeven | Creative telescoping using reductions[END_REF] seems to allude to a way of doing this, but the relevant space E may contain rational functions that are not in the image of L. For example take L = 1/z + 1/(z -1)σ -1 , α = 0, and A = {α}, then one can check that 1/z / ∈ Im(L) but 1/z ∈ E.

We write F for the summand in Eq. (3). It is a function of the four variables ν, n, z, λ. Basic properties of the Bessel function give the following four equations:

(λ 2 -1)zS ν (F ) + 2(n + 1)S n (F ) = 0, (4) (λ 2 -1)∂ λ (F ) -2nλF = 0, (5) (1 -λ 2 )z∂ z (F ) + +2(n + 1)S n (F ) + (λ 2 -1)(2n + ν)F = 0, (6) 4(n + 1)(n + 2)S 2 n (F ) + 4(λ 2 -1)(n + 1)(n + ν + 1)S n (F ) + z 2 (λ 2 -1) 2 F = 0, (7)
where S ν denotes the shift with respect to ν: S ν : G(ν) → G(ν + 1) and similarly for S n , while ∂ z and ∂ λ denote partial derivatives. These equations show that any shift or derivative S a ν ∂ b λ ∂ c z S d n F of F with nonnegative integers a, b, c, d rewrites as a Q(ν, λ, z, n)-linear combination of F and S n (F ). In particular, this implies that F is D-finite with respect to these variables. The aim of creative telescoping is to find a similar set of equations, in the variables ν, λ, z only, for the sum in Eq. [START_REF] Barrett | Particular solutions for nonhomogeneous, linear, ordinary difference equations[END_REF].

Let ∆ n be the difference operator

∆ n = S n -1. Any product ϕ(n)S n F with ϕ ∈ Q(ν, λ, z, n) can be rewritten ϕ(n -1)F + ∆ n (ϕ(n -1)F ), i.e.,
∂ α (F ) = R α F + ∆ n (G α ),
with R α a rational function in Q(ν, λ, z, n). This reduction as a sum of a rational function plus a difference is a general phenomenon (see Section 3.2). Reduction-based creative telescoping works by reducing this rational function further by pulling out parts that can be incorporated into the difference ∆ n (G α ). Denote by L the recurrence operator such that Eq. ( 7) is LF = 0. The adjoint of L (see Proposition 2) is

L * = 4(n -1)nS -2 n + 4(λ 2 -1)n(n + ν)S -1 n + z 2 (λ 2 -1) 2 where S -1 n : g(n) → g(n -1). Proposition 3 shows that a rational function R is of the form ∆ n (M (F )) for a recurrence operator M (S n ) if and only if R is in the image L * (Q(ν, λ, z, n))
. This is the basis for the computation of relations of the form [START_REF] Bostan | Creative telescoping for rational functions using the Griffiths-Dwork method[END_REF] where R α is now a reduced rational function (in a sense made precise in Definition 1). The next step is to look for linear combinations of these rational functions that yield telescopers.

The starting point is the monomial 1, which decomposes as

(9) 1 • F = 1 • F + ∆ n (0).
Using Eq. ( 5), the monomial ∂ λ rewrites

(10) ∂ λ (F ) = 2nλ λ 2 -1 F + ∆ n (0)
and the rational function is reduced. Taking the derivative of this equation and using Eq. ( 5) again gives a similar equation for ∂ 2 λ (F ):

∂ 2 λ (F ) = - 2n(λ 2 + 1) (λ 2 -1) 2 F + 2nλ λ 2 -1 ∂ λ (F ) + ∆ n (0), = - 2n(λ 2 + 1) (λ 2 -1) 2 F + (2nλ) 2 (λ 2 -1) 2 F + ∆ n (0).
This time, a reduction is possible. Indeed, Proposition 3 implies that L * (1)F is a difference ∆(A n ) (where A n can be computed explicitly). Since

L * (1) = 4λ 2 n 2 + 4((λ 2 -1)ν -1)n + z 2 (λ 2 -1) 2 ,
we can eliminate the term in n 2 in the expression of ∂ 2 λ (F ) to get

(11) ∂ 2 λ (F ) = - 2n(2ν + 1) (λ 2 -1) F -z 2 F + ∆ n A n (λ 2 -1) 2 .
A simple linear combination of Eqs. ( 9) to [START_REF] Chen | Constructing minimal telescopers for rational functions in three discrete variables[END_REF] then eliminates the term in n, showing that F satisfies the equation

λ∂ 2 λ F + (2ν + 1)∂ λ F + λz 2 F = ∆ n -λA n (λ 2 -1) 2 .
The left-hand side is a telescoper. The right-hand side is a certificate. It can be written more explicitly as

-λA n (λ 2 -1) 2 = - 4(nλ 2 + λ 2 ν -ν -1)nλ (λ 2 -1) 2 F - 4λ(n + 1)n (λ 2 -1) 2 S n (F ).
In general, summation and telescoping of the certificate requires verification. Here, we first observe that the certificate does not have integer poles and thus is well defined at all points over which it is summed. Next, the certificate evaluates to zero at n = 0. Finally, it tends to zero when n tends to infinity, as J ν+n (z) decreases fast as n → ∞ [27, 10.19.1].

In summary, we have obtained that the sum S in the right-hand side of Eq. ( 3) satisfies λ∂ 2 λ (S) + (2ν + 1)∂ λ (S) + λz 2 S = 0. Proceeding similarly with Eqs. ( 4) and ( 6), one gets the equations

zλS ν (S) + ∂ λ (S) = 0, z∂ z (S) -λ∂ λ (S) -νS = 0.
Injecting T = J ν (λz)/λ ν in these equations and using basic equations for J ν shows that it is a solution of this system too. The proof of the multiplication theorem is concluded by checking the equality of the initial conditions for T and for the sum on the right-hand side of Eq. ( 3). As ν is associated to the shift, we need to check initial conditions for any ν satisfying 0 ≤ Re(ν) < 1. Indeed, both term of the identity equal J ν (1) at z = 1, λ = 1, and ν ∈ [0, 1) and both their derivatives with respect to λ equal -J ν+1 (1), which proves the identity.

Background

In this section, we recall the basic framework for reducing creative telescoping to the generalized Abramov-Petkovšek decomposition. Most of this section is identical to the differential case [START_REF] Bostan | Generalized Hermite reduction, creative telescoping and definite integration of D-finite functions[END_REF]Sec. 4], except for the existence and computation of the cyclic vector and the use of the recurrence variant of Lagrange's identity [START_REF] Barrett | Particular solutions for nonhomogeneous, linear, ordinary difference equations[END_REF]. More gentle introductions to Ore algebras, creative telescoping and their applications can be found in the references [START_REF] Chyzak | The ABC of Creative Telescoping[END_REF][START_REF] Chyzak | Non-commutative elimination in Ore algebras proves multivariate holonomic identities[END_REF].

Telescoping ideal.

Ore algebras. Let k be a field of characteristic 0, x 0 , . . . , x m be variables used to form the fields of rational functions K = k(x 1 , . . . , x m ) and K = K(x 0 ). The Ore algebra A = K⟨∂ 0 , . . . , ∂ m ⟩ is a polynomial ring over K, with ∂ i ∂ j = ∂ j ∂ i , and a commutation between the ∂ i s and the elements of K ruled by relations ( 12)

∂ i R = σ i (R)∂ i + δ i (R), R ∈ K,
with σ i a ring morphism of K and δ i a σ i -derivation, which means that δ i (ab) = σ i (a)δ i (b) + δ i (a)b for all a, b in K [START_REF] Bronstein | An introduction to pseudo-linear algebra[END_REF][START_REF] Chyzak | Non-commutative elimination in Ore algebras proves multivariate holonomic identities[END_REF]. The typical cases are when ∂ i is the differentiation d/dx i (then σ i is the identity and δ i = d/dx i ) and the shift operator

x i → x i + 1 (then σ i (a) = a| xi←xi+1 and δ i = 0).
Annihilating and D-finite ideals. For a given function f in a left A-module, the annihilating ideal of f is the left ideal ann f ⊆ A of elements of A that annihilate f . A left ideal I of A is D-finite when the quotient A/I is a finite dimensional K-vector space. A function is called D-finite when its annihilating ideal is D-finite.

Telescoping ideal. As we focus here on summation, from now on, when we use n and S n , they stand for x 0 and the corresponding shift operator ∂ 0 :

x 0 → x 0 + 1.
The telescoping ideal T I of the left ideal I ⊂ A with respect to n is

T I = (I + ∆ n (A)) ∩ K⟨∂ 1 , . . . , ∂ m ⟩, where ∆ n = S n -1.
In other words, if I = ann F , the telescoping ideal T I is the set of operators T ∈ K⟨∂ 1 , . . . , ∂ m ⟩ such that there exists G ∈ A such that T + ∆ n G ∈ I, or equivalently, such that Eq. ( 1) holds (with t = n).

Cyclic vector and Lagrange identity.

Cyclic vector. Let I be a D-finite ideal of A and let r be the dimension of the K-vector space B := A/I.

An element γ ∈ B is called cyclic with respect to ∂ 0 if {γ, . . . , ∂ r-1 0
γ} is a basis of B. In the differential case (∂ 0 = d/dx 0 ), such a vector always exists and can be computed efficiently when I is D-finite [START_REF] Churchill | Cyclic vectors[END_REF]. In the shift case (∂ 0 : x 0 → x 0 + 1), even for a D-finite ideal I, it is not the case that there always exists a cyclic vector: in general, B decomposes as the sum of a vector space where ∂ 0 is nilpotent and a part where it is cyclic [START_REF] Jacobson | Pseudo-linear transformations[END_REF]. However, we have the following.

Proposition 1. [23, Thm. B2] With the notation above, in the case when ∂ 0 is the shift operator x 0 → x 0 +1, let E = (e 1 , . . . , e r ) T be a basis of the vector space B = A/I and A 0 ∈ Kr×r be defined by ∂ 0 E = A 0 E. If A 0 is invertible, then there exists a cyclic vector with respect to ∂ 0 of the form v = a 1 e 1 + • • • + a r e r with polynomial coefficients a i ∈ Z[x 0 ] of degree at most r -1, and coefficients all in {0, . . . , r}.

Sufficient conditions for the matrix A 0 to be invertible are that I = ann f with f in a K[∂ 0 , ∂ -1 0 ]module [START_REF] Hendricks | Solving difference equations in finite terms[END_REF] or that I be a reflexive ideal [START_REF] Van Der Hoeven | Creative telescoping using reductions[END_REF]. In practice, this condition on A 0 can be checked from the input and appears to be always satisfied in the examples we have tried. From this proposition, the computation of a cyclic vector follows the same lines as that of the differential case [START_REF] Churchill | Cyclic vectors[END_REF]. Most often, e 1 = 1 is a cyclic vector, which simplifies the rest of the computation.

Lagrange's identity. For our purpose, the shift version of Lagrange's identity can be viewed as giving an explicit form of the result of the left Euclidean division by the difference operator ∆ n , when applied to a left multiplication by a rational function.

Proposition 2. [3] Let u ∈ K and let L = r i=0 a i S i
n be an operator of order r with a i in K. The adjoint operator L * of L is defined as

L * = r i=0 a i (n -i)S -i
n and it satisfies

(13) uL -L * (u) = ∆ n P L (u)
where

(14) P L (u(n)) = r-1 i=0   r j=i+1 a j (n + i -j)u(n + i -j)   S i n .
Note that the term L * (u) denotes the evaluation of the operator L * at the rational function u, rather than the product of L * by u.

Let γ be a cyclic vector. Then any element of B is of the form Aγ with A ∈ K⟨S n ⟩. Applying Lagrange's identity with u = 1, L = A and multiplying on the right by γ shows that this is a rational multiple of γ up to a difference:

(15) Aγ = A * (1)γ + ∆ n P A (1)γ.
As in the differential case, all computations in B then reduce to K-linear operations on single rational functions, rather than vectors of them, by the following analogue of [6, Prop. 4.2]. Proposition 3. With the notation above, let γ be a cyclic vector of B = A/I and for all i = 0, . . . , m, let

B i ∈ K⟨S n ⟩ be such that ∂ i γ = B i γ. Then for all R ∈ K, ∂ i Rγ = φ i (R)γ + ∆ n Q i (R)γ, ( 16
)
with φ i (R) = B * i (R(x i + 1)), Q i (R) = P Bi (R(x i + 1)) if ∂ i : x i → x i + 1; φ i (R) = B * i (R) + d dxi (R), Q i (R) = P Bi (R(x i )) if ∂ i = d/dx i .
Proof. Multiplying Eq. ( 12) by γ on the right and using the definition of B i gives

∂ i Rγ = σ i (R)B i γ + δ i (R)γ.
The conclusion follows from Lagrange's identity (13) applied with L = B i and u = σ i (R). □

Algorithm 1 Creative Telescoping Algorithm

Input Generators of a D-finite ideal I Output A basis of the telescoping ideal T I of I γ ← a cyclic vector of B := A/I with respect to S n {See Section 3.2} L ← a minimal order operator in S n annihilating γ {See Section 3.3} φ 1 , . . . , φ n the maps described in Proposition 3 Initialize CanonicalForm {See Section 4}

R 0 ← CanonicalForm(A * 1 (1), L * ) {See Section 3.4} L ← [1] {List of monomials in ∂ 1 , . . . , ∂ m to visit} G ← {} {Gröbner basis} Q ← {} {Basis of the quotient} R ← {} {Set of reducible monomials} while L ̸ = ∅ do Remove the first element ∂ α of L if ∂ α is not a multiple of an element of R then if ∂ α ̸ = 1 then Take i such that ∂ α /∂ i ∈ Q; R α ← CanonicalForm(φ j (R ∂ α /∂i ), L * ) if there exists a K-linear relation between R α and {R β | R β ∈ Q} then (λ β ) R β ∈Q ← coefficients of the relation R α = R β ∈Q λ β R β Add ∂ α -R β ∈Q λ β ∂ β to G; Add ∂ α to R else Add ∂ α to Q for j = 1 to m do Append the monomial ∂ j ∂ α to L return G 3.
3. Canonical Form. Proposition 3 shows how, given a cyclic vector γ, all elements of B can be reduced to the product of γ by a rational function, up to a difference in ∆ n B. The starting point of the reduction-based creative telescoping is that one can actually identify those multiples that belong to ∆ n B. Proposition 4. With the same hypotheses as in Proposition 3, let L be a minimal-order operator in K⟨S n ⟩ annihilating γ, ie, the product Lγ is 0 in A/I, L has order r and no operator of order r -1 has that property. Then for all R ∈ K, □

Rγ ∈ ∆ n (B) ⇐⇒ R ∈ L * ( K). Proof. First, if R = L * (R ′ ) with R ′ ∈ K, Lagrange's identity (13) with u = R ′ and Lγ = 0 implies that Rγ = L * (R ′ )γ = ∆ n Gγ for G = -P L (R ′ ). Conversely, if Rγ ∈ ∆ n (B), there exists M ∈ K⟨S n ⟩ such that Rγ = ∆ n M γ. The operator ∆ n M -R
This proposition motivates the following.

Definition 1. [6] A canonical form associated to L * is a K-linear map [•] : K(n) → K(n) such that for all R ∈ K(n), [L * (R)] = 0 and R -[R] ∈ L * (K(n)). A rational function R ∈ K(n) is called reduced when [R] = R.
The computation of canonical forms is the object of Section 4.

3.4.

Creative Telescoping Algorithm via Canonical Forms. With the notation above, the creative telescoping algorithm from [START_REF] Bostan | Generalized Hermite reduction, creative telescoping and definite integration of D-finite functions[END_REF] applies verbatim. It is given in Algorithm 1. Its principle is to iterate on every monomial of the form ∂ α1 1 . . . ∂ αm m by increasing order for some monomial order, e.g., the grevlex order, and to compute the reduced rational functions R α such that ( 17)

∂ α = R α γ + ∆ n G α mod I.
The rational function R α is obtained by

(18) R 0,...,0 = [A * 1 (1)], R α = [φ i (R β )] if ∂ α = ∂ i ∂ β ,
where A * 1 is the adjoint of the operator A 1 verifying 1 = A 1 (γ). When a monomial α is dealt with, two situations are possible. The corresponding R α can be a linear combination of the previous R β . In that case, that linear combination makes the corresponding linear combination of ∂ α and the ∂ β a newly discovered element of the telescoping ideal T I and then it is not necessary to visit the multiples of this monomial. Otherwise, ∂ α is free from the previous ones and thus a new generator of B has been found. The algorithm terminates when there are no more monomials to visit.

The only difference with the differential case lies in the definition of the canonical form [•] associated to the adjoint L * of the minimal-order operator L ∈ K(n)⟨S n ⟩ annihilating γ. By Propositions 3 and 4 and the definition of canonical form, Eq. ( 17) is satisfied and the following equivalence holds:

(19) a 1 R 1 + • • • + a s R s = 0 iff a 1 D 1 + • • • + a s D s ∈ T I .
The following result follows, with the same proof as in [START_REF] Bostan | Generalized Hermite reduction, creative telescoping and definite integration of D-finite functions[END_REF].

Theorem 1. Given as input the generators of a D-finite ideal I and a cyclic vector γ for S n , Algorithm 1 terminates if and only if T I is D-finite. Then, it outputs a Gröbner basis of T I for the grevlex order.

As in the differential situation [START_REF] Bostan | Generalized Hermite reduction, creative telescoping and definite integration of D-finite functions[END_REF], one can modify the algorithm to compute all elements of T I up to a given bound on their degree, or to return as soon as one telescoper is found, thus allowing to recover a generating family of a sub-ideal of T F .

Generalized Abramov-Petkovšek decomposition

The main contribution of this article is an algorithm for the computation of canonical forms as in Definition 1 for the operator

L * = r i=0 p i (n)S -i n with polynomial coefficients p i ∈ K[n].
The modified Abramov-Petkovšek decomposition [START_REF] Chen | A modified Abramov-Petkovšek reduction and creative telescoping for hypergeometric terms[END_REF] is a special case of this reduction when L has order 1 and once the shell [START_REF] Chen | A modified Abramov-Petkovšek reduction and creative telescoping for hypergeometric terms[END_REF] has been removed [START_REF] Bostan | Generalized Hermite reduction, creative telescoping and definite integration of D-finite functions[END_REF]Sec. 3.5.3].

The starting point is a decomposition of any rational function R ∈ K(n) in the form

(20) R(n) = P ∞ (n) + i,h c i,h (n) Q i (n + h) ℓ i,h , with ℓ i,h ∈ N * , polynomials P ∞ , Q i and c i,h in K[n] such that deg c i,h < ℓ i,h deg Q i and gcd(Q i (n + k), Q j (n)) = 1 for all k ∈ Z when i ̸ = j. This is discussed in Section 4.1.
The vector spaces

K Qi (n) def = Vect K n ℓ Q i (n + h) j | h ∈ Z, j ∈ N * , ℓ < j deg(Q i )
are in direct sum for distinct Q i and are left invariant by L * modulo K[n]. This allows the reduction algorithm to operate in each of the K Qi (n) independently. This is described in Sections 4.2 and 4.3, before the reduction of the remaining polynomial part in Sections 4.4 and 4.5.

Notation 1. For two integers a, b with a ≤ b, we write a; b for the set {a, a + 1, . . . , b}.

4.1.

Decomposition of rational functions. Recall that a polynomial Q is square-free when it does not have multiple nontrivial factors. It is shift-free when gcd(Q

(n), Q(n + k)) = 1 for all k ∈ Z * .
A shiftless decomposition of a polynomial Q is a factorization of the form

Q = v i=1 ni j=1 Q i (n + h i,j ) ei,j , Algorithm 2 Weak reduction of poles [•] Q Input R, Q Output a reduced form of R while there exists j < 0 (minimal) such that Q(n -j) | den(R) do R ← R -L * A(n) Q(n-j) ord j (den(R))+ord j (p 0 )
with A as in Eq. ( 22)

while there exists j ≥ r (maximal) such that Q(n -j) | den(R) do R ← R -L * A ′ (n+r) Q(n-j+r) ord j (den(R))+ord j (pr )
with A ′ as in Eq. ( 23) return R with e i,j ∈ N * , h i,j ∈ Z, and

Q i ∈ K[n] are such that each Q i is square-free and gcd(Q i (n + k), Q j (n)) = 1
for all i, j and all k ∈ Z unless i = j and k = 0. Such a factorization can be computed using only gcds, resultants and integer root finding [START_REF] Gerhard | Shiftless decomposition and polynomial-time rational summation[END_REF].

Note that shiftless decompositions are not unique in general. One can be refined when a Q i is not irreducible, by splitting this factor further. In particular, the linear factors of the Q i can be isolated and dealt with more easily.

A polynomial Q is refined with respect to a polynomial P when it is such that for each h ∈ Z, there exists ℓ ∈ N such that gcd(P, Q(n + h) ℓ+1 ) = Q(n + h) ℓ . A shiftless decomposition is called refined with respect to P when each Q i is. This refinement can be computed using gcds only and will be used with P = p 0 and P = p r , the extreme coefficients of L * .

From a shiftless decomposition, the partial fraction decomposition of Eq. ( 20) is then obtained by standard algorithms [34, 5.11].

Weak reduction of the polar part.

Lemma 1. Let Q ∈ K[n] be square-free, shift-free and refined with respect to the coefficients p 0 and p r of L * . Given a rational function

R ∈ K Q (n), Algorithm 2 computes a rational function [R] Q ∈ K Q (n)
with all its poles at zeros of Q(n -j) such that j ∈ 0; r -1 and R -[R] Q = P + L * (T ) for some P ∈ K[n] and T ∈ K Q (n). The algorithm is K-linear.

Proof. Assume that R decomposes as

(21) R = j∈J λ j (n) Q(n -j) sj with deg(λ j (n)) < s j deg(Q).
Let j m = min(J) and ord j (p 0 ) be the largest integer ℓ such that Q(n -j) ℓ | p 0 . Then,

Q(n -j m ) sj m L * 1 Q(n -j m ) sj m +ordj m (p0) = p0 (n) mod Q(n -j m ) sj m ,
where p0 (n) is the remainder in the Euclidean division of p 0 /Q(n -j m ) ordj m (p0) by Q(n -j m ) sj m . The poles of this rational function are at zeros of Q(n -j) with j ∈ Ĵ := J \ {j m } ∪ (j m + 1, r -1 ).

Since Q is reduced with respect to p 0 , the polynomial p0 (n) is relatively prime with Q(n -j m ). Thus, there exist polynomials A and B such that

(22) λ jm (n) = A(n)p 0 (n) + B(n)Q(n -j m ) sj m . Then A(n)p 0 (n) Q(n -j m ) sj m = λ jm (n) Q(n -j m ) sj m -B(n), so that R -L * A(n) Q(n -j m ) sj m +ordj m (p0) is equivalent to R modulo L * (K Q (n))
and with all its poles at zeros of Q(n -j) with j ∈ Ĵ. This operation can be repeated a finite number of times until all poles are at zeros of Q(n -j) with j ≥ 0.

Similarly, let j M = max(J) and ord j M +r (p r ) be the largest integer ℓ such that Q(n-j M ) ℓ divides p r (n-r). Then

Q(n -j M ) sj M L * 1 Q(n -j M + r) sj M +ordj M +r (pr) = pr (n) mod Q(n -j M ) sj M ,
where pr (n) is the remainder in the Euclidean division of p r /Q(n -j M ) ordj M +r (pr) by Q(n -j M ) sj M . The poles of this rational function are at zeros of Q(n -j) with j ∈ Ĵ′ := J \ {j M } ∪ (j M -1, r -1 ). Again, since Q is reduced with respect to p r , the polynomial pr is relatively prime with Q(n -j M ). Thus there exist two polynomials A ′ and B ′ such that ( 23)

λ j M (n) = A ′ (n)p r (n) + B ′ (n)Q(n -j M ) sj M so that R -L * A ′ (n + r) Q(n -j M + r) sj M +ordj M +r (pr) is equivalent to R modulo L * (K Q (n))
and with all its poles at zeros of Q(n -j) with j ∈ Ĵ′ . This operation can be repeated a finite number of times until all poles are at zeros of Q(n -j) with j ∈ 0, r -1 .

Each step being K-linear, so is the algorithm.

□ Example 1. Let R = 8nx + n(8x 2 + x) -16x -1 2(n -1) 2 - 4(x -1)x 2 n - (4x 3 + 8x 2 -31x -32)n + 4x 3 -31x -32 2(n + 1) 2 + 4x(x 2 -x -8) n + 2 + 2x 3 n + 3 and L * = x 2 (n -2)S -3 n -n(4n 2 -x 2 -4n)S -2 n + n(4n 2 -x 2 -4n)S -1 n -x 2 (n + 2
). The poles of R are at {1, 0, -1, -2, -3}. We take Q = n + 1 and follow the steps of the algorithm.

The pole at -3 is easy: from

L * 1 n + 3 = -x 2 + 4n + -x 2 + 8 n + 1 + 2x 2 -16 n + 2 + x 2 n + 3
and the coefficient 2x 3 of (n + 3) -1 in R, the algorithm performs the subtraction

R ← R -2xL * 1 n + 3 = n(8x 2 + x) -16x -1 2(n -1) 2 + 4x 2 n - (8x 2 + x -32)n + x -32 2(n + 1) 2 - 4x 2 n + 2 .
Next, the pole -2 is a simple root of the constant coefficient of L * , leading to the computation of

L * 1 (n + 2) 2 = x 2 (n -2) (n -1) 2 + x 2 n - n(x 2 -4) -4 (n + 1) 2 - x 2 n + 2
so that the pole is removed by

(24) R ← R -4L * 1 (n + 2) 2 = x 2(n -1) - x 2(n + 1)
.

R now has all its poles in {-1, 0, 1} and the weak reduction is finished.

4.3.

Strong reduction of the polar part. By Corollary 1, the weak reduction produces rational functions all whose poles differ from those of Q by an integer in 0, r -1 . The next step of the reduction is to subtract rational functions in L * (K Q (n)) that have this property. It turns out to be possible to focus on a finite-dimensional subspace of L * (K Q (n)) thanks to the following.

Lemma 2. If j < 0, s > ord j (p 0 ) and ℓ < s deg(Q) or if j ≥ 0, s > ord j+r (p r ) and ℓ < s deg(Q) then

L * n ℓ Q(n -j) s Q = 0.

Algorithm 3 Weak reduction of polynomials [•] ∞

Input P and (σ,p) from Eq. ( 25) Output a reduced form of P a ← 0 while deg(P ) ≥ σ do if deg(P ) -σ is a root of p then a ← a + lt(P ); P ← P -lt(P ) else P ← P -lc(P )

p(deg(P )-σ) L * (n deg(P )-σ ) return a + P Proof. Let j, s, ℓ be three integers that satisfy the first assumption. Then L * (n ℓ /Q(n -j) s ) has a denominator that is divisible by Q(n -j) with j < 0 by assumption. No smaller k is such that Q(n -k) divides the denominator. Thus the first pass through the first loop of the weak reduction subtracts L * (n ℓ /Q(n -j) s ) to itself and reduces it to zero. When the second assumption is satisfied, then L * (n ℓ /Q(n -j) s ) has a denominator that is divisible by Q(n -(j + r)) with j + r ≥ r by assumption. No larger k is such that Q(n -k) divides the denominator. Thus again, the second loop reduces that fraction to 0. □

Corollary 1. Let I 0 := {j ∈ Z <0 | gcd(p 0 (n), Q(n-j)) ̸ = 1} and I r := {j ∈ Z ≥0 | gcd(p r (n), Q(n-j)) ̸ = 1}. The K-vector space [L * (K Q (n))] Q is generated by the fractions L * n ℓ Q(n -j) sj Q , with      j ∈ I 0 and 1 ≤ s j ≤ ord j (p 0 ) and 0 ≤ ℓ < s j deg Q or j ∈ I r and 1 ≤ s j ≤ ord j+r (p r ) and 0 ≤ ℓ < s j deg Q.
Corollary 1 gives a generating family of the finite dimensional K-vector space [L * (K Q (n))] Q . These rational functions can be written in the basis (n i /Q(n -j) k ) i,k∈N,j∈Z and one can then compute an echelon basis of this finite-dimensional space. This precomputation step corresponds to the computation of the B Qi 's in Algorithm 4. The strong reduction of a rational function R ∈ K Q then consists in reducing [R] Q with this echelon basis. By this process, we obtain the following.

Proposition 5. Strong reduction reduces every rational function

R ∈ L * (K Q (n)) to a polynomial in K[n].
Example 2. With the same notation as in example 1, Corollary 1 shows that [L * (K n+1 (n))] n+1 is generated by

L * 1 n -2 n+1 = 4n + x 2 n + 1 - x 2 n -1 , L * 1 n -1 n+1 = 4n + x 2 n -1 - x 2 n + 1 .
Thus the strong reduction of the rational function R from Eq. ( 24) is the polynomial

R + L * ((n -2) -1 )/(2x) = -2n/x,
concluding the reduction.

4.4.

Weak reduction of polynomials. The weak reduction of polynomials is a direct adaptation of the differential case [START_REF] Bostan | Generalized Hermite reduction, creative telescoping and definite integration of D-finite functions[END_REF]. The indicial polynomial of L * at infinity is the polynomial p ∈ K[s] defined by

(25) L * (n s ) = n s+σ (p(s) + O(1/n)),
with σ ∈ N. The ensuing weak reduction is presented in Algorithm 3.

Example 3. In Examples 1 and 2, the indicial equation at infinity is

L * (n s ) = n s+2 (8 + 4s + O(1/n)).
The polynomial -2n/x found in Example 2 cannot be reduced further by weak reduction since its degree in n is smaller than 2.

The following properties are proved exactly as those for weak reduction at a pole.

Lemma 3. Algorithm [•] ∞ terminates and is K-linear. For all P ∈ K[n], there exists Q ∈ K[n] such that P -[P ] = L * (Q). If s ∈ N is not a root of p, then [L * (n s )] ∞ = 0.

Certificates

Reduction-based creative telescoping algorithms allow to find a telescoper without having to compute an associated certificate. This has led to faster algorithms as certificates are known to be larger than telescopers [START_REF] Bostan | Creative telescoping for rational functions using the Griffiths-Dwork method[END_REF]. This approach makes sense in the differential case when it is known in advance that the integral of a certificate over a cycle that avoids singularities is equal to zero. The framework is not as favorable for sums. Indeed, it is necessary to detect whether the certificate has poles in the range of summation and it is often unclear whether the certificate becomes 0 at the boundaries of the summation interval.

It is however possible to compute the certificates in a compact way during the execution of our algorithm, with almost no impact on the execution time. The idea is to make the computation and storage of certificates efficient by storing them as directed acyclic graphs (dags) rather than operators with normalized rational function coefficients. These dags have a number of internal nodes of the same order as the number of operations performed when computing the telescoper, so that their computation does not burden the complexity. They can then be evaluated at the endpoints of the range of summation, or expanded in Laurent series there. 5.1. Computation and structure of certificates. Equations ( 17) and [START_REF] Chyzak | A computer-algebra-based formal proof of the irrationality of ζ(3)[END_REF] show how one can compute telescopers without computing certificates. We now show how to compute the certificates simultaneously. These are obtained as sums of monomials in the ∂ i multiplied by rational functions. These rational functions have denominators of the form Q(n -j) s , with integers j, s and polynomials Q that occur either in the shiftless decomposition of p 0 p r or in the denominator of the operator A 1 ∈ K⟨S n ⟩ such that 1 = A 1 γ, or in the denominator of one of the operators B i from Proposition 3. Those sums are not reduced to a common denominator. They share many common coefficients and denominators that are efficiently compacted into dags by sharing common subexpressions (this is how Maple stores them by default).

The starting point is the cyclic vector γ ∈ A/I and the operator A 1 ∈ K⟨S n ⟩ such that 1 = A 1 γ. By Euclidean division by ∆ n on the left, A 1 = R 0 + ∆ n g 0 . In general, the certificate G α in Eq. ( 17) is stored as an (unreduced) element g α of A such that G α = g α γ mod I. The computation is incremental. It is initialized with R 0 and g 0 as above, corresponding to α = 0. In many cases, the vector 1 is cyclic, so that one can take γ = 1 = A 1 = R 0 and g 0 = 0. Otherwise, the denominators of R 0 and g 0 are shifts of the denominator of A 1 .

Next, Eq. ( 18) leads to the computation of the canonical form of the rational function ϕ i (R β ). This computation is performed via a sequence of reductions which consist of subtractions of elements in L * (K(n)). Keeping track of these rational functions (without normalizing them) gives the canonical form R α as

ϕ i (R β ) = R α -L * (c α )
for some rational function c α ∈ K. Both R α and c α have denominators that are shifts of those of R β or of factors of p 0 p r . In view of Eq. ( 16), it follows that

∂ α = ∂ i ∂ β = R α γ + ∆ n (∂ i G β + P Bi (σ i (R β ))γ + P L (c α )γ) mod I.
Thus, the certificate G α is obtained as g α γ with [START_REF] Koutschan | A fast approach to creative telescoping[END_REF]

g α = ∂ i g β + P Bi (σ i (R β )) + P L (c α ).
This proves the claim concerning the structure of the certificates and the factors of their denominators. In our implementation, ∂ i is commuted with the coefficients of the certificate g β only. If desired, one can further use ∂ i = B i mod I so as to write g α as an operator in S n only.

5.2.

Evaluation of the certificates. The output of Algorithm 1 is a set of elements T of the telescoping ideal, which means that

T (F ) = g(F )(n + 1, x 1 , . . . , x m ) -g(F )(n, x 1 , . . . , x m ),
with g a certificate as described above. Summing over n, the right-hand side telescopes and only the values of the certificate g(F ) at the endpoints are needed. It is possible to prove that these evaluations are zero without any evaluation in two important cases. First, if the summand F has finite support (e.g., binomial sums), then the sum of any certificate over Z will be zero provided it has no pole in the summation range. The second case is when one can prove that R(n)γ(F ), . . . , R(n)S r-1 n γ(F ) tend to zero as n tends to ±∞ for any rational function R ∈ K[n] (as in the introductory example). Then again the sum of any certificate over Z will be zero provided it has no pole in the summation range. 5.3. Integer pole detection. By its very nature, the method of creative telescoping requires the certificate not to have poles in the range of summation, so that telescoping can occur. The structure of the certificates described above does not allow the efficient computation of its denominator exactly. However it is possible to compute a multiple of it by taking the least common multiple of the denominators of every rational function in the representation. This can be done efficiently by performing the computation on the dag representation of the certificates.

From this multiple of the denominator of the certificate, one can compute the set of roots that lie in the summation range; this amounts to computing the roots that differ from the endpoints of the summation range by an integer. If that set is not empty, then one can compute a Laurent series expansion of the certificate at any point to check whether it is a pole or not, again by exploiting the dag representation of the rational function coefficients. 2 where J n (x) is the Bessel function of the first kind, Algorithm 1 outputs the telescoper ∂ x and a certificate G in dag form that exhibits poles at n ∈ {-1, 0, 1, 2}. Our implementation produces the polynomial n(n + 1)(n -2)(n -1) containing this information. In this example, the certificate is small enough that it can easily be normalized and one gets its value as

(x) = ∞ n=1 J n (x)
- x 4 (n + 1) S 2 n + n + 1 x S n - 8n 2 -x 2 + 8n 4x (n + 1) ,
showing that the poles at 1 and 2 vanish in the normalization. Without normalizing the certificate, one can still evaluate the series expansions of the certificate at those points to establish that it has a finite limit there, making the summation legitimate. At n = 1, the evaluation is found to be

-J 0 (x)J 1 (x) = 1 2 (J 0 (x) 2 ) ′ .
Thus we have proved

d dx ( 1 2 J 2 0 (x) + J 1 (x) 2 + J 2 (x) 2 + • • • ) = 0
. This shows that the sum is constant and the value is revealed by its value at 0, which follows from J k (0) = 0 for k > 0 and J 0 (0) = 1, so that in the end, we recover the classical identity [27, 10.23

.3] 1 = J 0 (x) 2 + 2 k≥1 J k (x) 2 .
T := (n + 2) 3 S 2 n -(2n + 3) 17n 2 + 51n + 39 S n + (n + 1) 3 .
The singularities of the certificate obtained by our implementation are at k ∈ {n + 1, n + 2}. Indeed, once normalized, the certificate is found to be

C := 4 (3k -4n -8) k 4 (k -n -2) 2 S n - 4k 4 (3k + 4n + 4) (k -n -1) 2 .
Let U n,k denote the product of binomials in the sum. Summing T U n,k from k = 0 to k = n + 2 gives T (A n ). If telescoping is legitimate, then the values of the endpoints are the values of CU n,k at k = 0 and k = n + 3, that are both easily checked to be 0. For this to allow to conclude that T (A n ) = 0, it is then sufficient to check that CU n,k is not singular at k = n + 1 and k = n + 2, even though C is. Indeed, a series expansion of the evaluation of CU n,k at k = n + 1 and k = n + 2 is possible for our implementation, and finds that the sequence has a finite limit there, which concludes the proof that the telescoper Eq. ( 28) cancels the sum Eq. ( 27); see also [START_REF] Chyzak | A computer-algebra-based formal proof of the irrationality of ζ(3)[END_REF] for more on these issues.

5.5.

A larger example. The computation of telescoper and certificate for the sum

∞ n=0 (4n + 1)(2n)! n! 2 2 2n √ x J 2n+1/2 (x)P 2n (u) (29)
takes less than 15 sec. with our current implementation (see Table 1). The telescopers are quite small:

(1 -u 2 )∂ u + xu∂ x , (u 3 -u)∂ 2 u + (1 + u 2 )∂ u -u 3 x 2 .
In this example, not normalizing the certificates during their computation has a cost. The actual certificates, once reduced by the Gröbner basis of the annihilating ideal of the summand, are not very large. They are easily computed by Koutschan's program. Still, the corresponding dags are large. Nonetheless, it takes less than 1 sec. to compute a multiple of the denominators of the certificates and detect that they do not have integer roots. Evaluating the certificates at n = 0 using their dag representation takes less than 2 min. and proves that the telescopers cancel the sum in Eq. ( 29).

Implementation

This algorithm is implemented in Maple 2 . Table 1 gives a comparison of our code with Koutschan's heuristic (HF-FCT) and Chyzak's algorithm (HF-CT) 3 . They are both implemented in Koutschan's Holo-nomicFunctions package in Mathematica [START_REF] Koutschan | Advanced Applications of the Holonomic Systems Approach[END_REF]. The column 'redctsum' corresponds to our algorithm.

These programs have been executed on a list of 21 easy examples that were compiled by Koutschan, as well as more difficult ones given in Eqs. [START_REF] Takayama | Gröbner basis, integration and transcendental functions[END_REF] to (39) below. Eq. ( 30) comes from recent identities involving determinants [START_REF] Amdeberhan | Yay for determinants! Technical Report 2023-18[END_REF], Eqs. [START_REF] Van | Constructing reductions for creative telescoping[END_REF], [START_REF] Van Der Hoeven | Creative telescoping using reductions[END_REF], [START_REF] Zeilberger | A holonomic systems approach to special functions identities[END_REF] and [START_REF] Zeilberger | The method of creative telescoping[END_REF] have been chosen because they looked natural to experiment with, Eq. (34) is a harder example found in Koutschan's list, Eq. (33) as well as Eq. ( 36) and it special case Eq. ( 35) come from the classical book of integral and series by Prudnikov et al. [START_REF] Prudnikov | Integrals and Series[END_REF], and finally Eq. (39) is an example where Koutschan's heuristic does not stop as it does not guess correctly the form of the ansatz to use [START_REF] Chen | Constructing minimal telescopers for rational functions in three discrete variables[END_REF]. 2 The implementation is available at https://github.com/HBrochet/CreativeTelescoping.git, together with sessions of examples. 3 The code was run on a Intel Core i7-1265U with 32 GB of RAM. For any r, our algorithm produces a minimal telescoper of order r and degree r(r -1)/2. The timings are reported in Table 2. It is unclear why the heuristic HF-FCT does not perform well on this family. On most of these examples, the main part of the time of the computation is spent in the reductions in the call to CanonicalForm in Algorithm 1. For the two similar sums of Eqs. [START_REF] Van Der Hoeven | Constructing reductions for creative telescoping: the general differentially finite case[END_REF] and [START_REF] Von | Modern computer algebra[END_REF], almost half of the time is spent in Algorithm 4 performing the reductions needed to compute the bases for the strong reduction. This step is crucial to ensure that the minimal order elements in the telescoping ideal are found.

There are cases, like Eq. ( 36) and the family S r , where the intermediate rational functions R α in Eq. ( 17) become much larger than the telescopers found after linear algebra on them. In such situations, the direct, non-incremental approach taken by HF-CT and HF-FCT can be more efficient, by avoiding an unnecessarily large basis of rational functions.

  annihilates γ. By minimality of L, there exists N ∈ K⟨S n ⟩ such that ∆ n M -R = N L. Taking the adjoint and evaluating at 1 gives R = M * ∆ * n -L * N * and finally R = -L * (N * (1)).

5. 4 .

 4 Examples. 5.4.1. Neumann's Addition Theorem for Bessel functions. On input S

5. 4 . 2 . 2 ,

 422 Apéry's Sequence. The classical sum used by Apéry in his proof of the irrationality of ζ(3), has telescoper[START_REF] Prudnikov | Integrals and Series[END_REF]

P(- 1 )

 1 m -i + j c n,j where c n,j satisfies recurrences of order 2 [2, p. k (4k + 1)J 2k+1/2 (w)P 2k (z) n (x)P n (y)P n (1/2) (37)∞ n=0 P n (x)P n (y)P n (z) (38) y 4x + 2 (45x + 5y + 10z + 47)(45x + 5y + 10z + 2)(63x -5y + 2z + 58)(63x -5y + 2z -5) (39)The family (S r ) is defined by[START_REF] Gillis | On elementary methods in positivity theory[END_REF] k (rn -(r -1)k)!(r!) k(n -k)! r k! .

Table 1 .

 1 Timings. The notation ( * ) means that we could not check whether the telescopers returned by HF-FCT were minimal.

		HF-CT HF-FCT redctsum
	easy examples 6.7s	7s	0.9s
	Eq. (30)	101s	49s	0.8s
	Eq. (31)	52s	4s	1.4s
	Eq. (32)	62s	1.7s	5.7s
	Eq. (33)	4.9s	1.4s	10.3s
	Eq. (34)	4.9s	1.4s	13.5s
	Eq. (35)	1200s	13s	205s
	Eq. (36)	> 6h	108s	3338s
	Eq. (37)	1703s	4.7s	580s
	Eq. (38)	> 1h	3.2s( * )	> 1h
	Eq. (39)	> 1h	> 1h	0.4s

Table 2 .

 2 Timings on the family S r from Eq. (40).

	HF-CT HF-FCT redctsum
	S 6 11s	64s	0.4s
	S 7 32s	331s	0.6s
	S 8 106s	1044s	1.0s
	S 9 325s	3341s	2.5s
	S 10 1035s	>1h	5.7s

Algorithm 4 PrecomputeBases

Input Q 1 , . . . , Q v polynomials that occur in the shiftless decomposition of p 0 p r Output The echelon bases B Q1 , . . . , B Qv , B pol B pol ← {} for i = 1 to v do

| l, j, s j as in Cor.1

Algorithm 5 Reduction of rational functions [•]

Input R and B Q1 , . . . , B Qv , B pol computed by Algorithm 4 Output a reduced form of R by L * (K(n)).

Decompose R as

Here again, a finite number of generators can be obtained thanks to the following.

Lemma 4. Let Q 1 , . . . , Q v be the polynomials that occur in a shiftless decomposition of p 0 p r and let P be a polynomial in L * (K(n)). Then

) is a polynomial, then the poles of R must be cancelled by the zeros of p 0 p r or their shifts. It follows that R decomposes as