Alexis Fouilloy
email: fouilloy_a@univ-corse.fr

Cyril Voyant
email: cyrilvoyant@gmail.com

Gilles Notton
email: notton@univ-corse.fr

Fabrice Motte
email: motte_f@univ-corse.fr

Christophe Paoli
email: cpaoli@univ-corse.fr

Marie-Laure Nivet
email: nivet@univ-corse.fr

Emmanuel Guillot
email: emmanuel.guillot@promes.cnrs.fr

Duchaud Jean-Laurent
email: duchaud_jl@univ-corse.fr

Solar irradiation prediction with machine learning: forecasting models selection method depending on weather variability

Keywords: Time Series forecasting, machine learning, variability, ARMA, ANN, Regression tree, Gaussian process, SVR

Eleven statistical and machine learning tools are analyzed and applied to hourly solar irradiation forecasting for time horizon from 1 to 6 hours. A methodology is presented to select the best and most reliable forecasting model according to the meteorological variability of the site. To make the conclusions more universal, solar data collected in three sites with low, medium and high meteorological variabilities are used: Ajaccio, Tilos and Odeillo. The datasets variability is evaluated using the mean absolute log return value. The models were compared in term of normalized root mean square error, mean absolute error and skill score. The most efficient models are selected for each variability and temporal horizon: for the weak variability, auto-regressive moving average and multi-layer perceptron are the most efficient, for a medium variability, auto-regressive moving average and bagged regression tree are the best predictors and for a high one, only more complex methods can be used efficiently, bagged regression tree and the random forest approach.

Introduction

The management of an electrical network is a very complex task, particularly on small grids like islands which are generally not interconnected to the mainland grid [START_REF] Diagne | Review of solar irradiance forecasting methods and a proposition for small-scale insular grids[END_REF]. The equality between production and consumption is the most difficult challenge; moreover, the increased rate of intermittent and stochastic renewable energy sources into the energy mix adds a layer of difficulty. For an efficient management of the energy mix and a safest electricity supply, production and consumption must be planned before. There are many forecasting methods for solar radiation, all the methods have advantages and inconvenient and their accuracy depends on the forecasting horizon and on the geographical situation. The three main categories of models are sky imaging, physical models (or numerical weather prediction, e.g. NWP) and machine learning models. The sky imaging models are mainly used for short term forecasting (1 to 60 minutes), the physical models concern numerical weather prediction models often used for long term forecasting (several days). The mediumterm forecasting (from 1 to few hours) is often made by machine learning methods, a large state of the art was made concerning the solar power forecasting (Voyant, Notton, et al., 2017). This study is related to the comparison of eleven machine learning models in order to forecast solar global irradiation and to compare their performances in three different sites to make our results more universal. The main goal is to propose a prioritization methodology to highlight the best forecasting method between the tested ones according to the level of irradiation variability considering the complexity of the method. To distinguish the different weather conditions and to characterize the global behavior of the solar time series in the three sites, a property of time series named variability is calculated, the characterization of this time series property can be made by various statistical parameters [START_REF] Voyant | Statistical parameters as a means to a priori assess the accuracy of solar forecasting models[END_REF].

The solar forecasting is realized with an hourly time granularity for 1 hour to 6 hours' time horizons. The comparison between the models is evaluated by statistical indexes and classical error metrics, the eleven tested models are:

-Persistence (P) and smart persistence (SP) can be considered as naïve models, they are often used in intraday solar radiation forecasting and are usually taken as reference to compare the performances with other models [START_REF] Lauret | A benchmarking of machine learning techniques for solar radiation forecasting in an insular context[END_REF].

-Auto regressive moving average (ARMA) and artificial neural network (ANN) are two types of models allowing to realize respectively linear and non-linear regressions.

ARMA is certainly the most studied tool especially in econometry, but also in global irradiation forecasting [START_REF] Boland | Short term solar radiation forecasting: Island versus continental sites[END_REF]. Artificial neural network (ANN) is taken as well-known global radiation forecasting method (Voyant, Motte, et al., 2017).

-Regression tree based models family is composed by five models with different levels of complexity and optimization. Classical regression tree (RT) and pruned regression tree (RT-pruned) are sometimes used in very short term forecasting and allow good results compared with other machine learning models (Troncoso, Salcedo-Sanz, Casanova-Mateo, Riquelme, & Prieto, 2015). Boosted regression tree (RT-boosted) consists in combine the responses of a collection of weak classifiers which once averaged form a strong classifier to make prediction [START_REF] De | Boosted Trees for Ecological Modeling and Prediction[END_REF]. Then the aggregated versions of regression trees denoted bagged regression tree (RT-bagged) [START_REF] Breiman | Bagging predictors[END_REF] are developed to build an aggregated predictor in order to improve accuracy of the forecast. Finally, we propose also the random forest method (RF) which is another improvement of regression tree methods with more robustness [START_REF] Breiman | Random forests[END_REF].

-Models based on Kernel estimation like Gaussian processes (GP) and support vector regression (SVR). The first one is a nonparametric kernel-based probabilistic model giving very good results in time series prediction [START_REF] Rasmussen | Gaussian Processes in Machine Learning[END_REF]), The second one (SVR) can be considered as a generalization of linear classifiers and can be used to forecast solar irradiation and can be improved by several methods of optimization [START_REF] Jiang | Global horizontal radiation forecast using forward regression on a quadratic kernel support vector machine: Case study of the Tibet Autonomous Region in China[END_REF]).

The study is developed as follow: section 2 describes the materials and methods, the data and the models. Section 3 presents the results and the discussion about the forecasting reliability, the variability of solar data and the link between variability and efficiency of the forecasting models.

Materials and methods

This section presents information related to the preliminary stage before the forecasting phase. 1. A quality control of the solar data: often, mistakes appear in the temporal series of solar data due to problems with the acquisition system; an automatic quality check used in the frame of GEOSS project (Group on Earth Observation System of System) (http://www.earthobservations.org/geoss.php) has been applied to the data; 2. A data preprocessing: the night hours are removed, due to problems occurring at the sunset and the sunrise (mask effects and bad response of pyranometers) [START_REF] Badescu | Modeling solar radiation at the earth's surface: recent advances[END_REF][START_REF] Iqbal | An introduction to solar radiation[END_REF], a filter is applied on the datasets which remove all the data that correspond to a solar elevation angle up to 10°.

3. Stationarity: the solar irradiation time series contains some properties like seasonality and periodicity. Most of the machine learning forecasting methods are applicable only to stationary time series [START_REF] Hornik | Multilayer feedforward networks are universal approximators[END_REF]) describe the necessary steps in order to use data in forecasting with machine learning. Thus, the solar data time series must be made stationary using a solar radiation model by clear sky and computing the clear sky index (paragraph 2.1), this method is described in (Paoli, Voyant, Muselli, & Nivet, 2010). 4. At last, the choice of the number of input data for each forecasting tool is realized by auto mutual information method (Luo, Shi, Zheng, Gang, & Cai, 2017) and described in paragraph 2.2.

The clear sky model (CS) and clear sky index (CSI)

The solar radiation by clear sky computing is the third step of the preprocessing method, it consists in calculate the maximum global horizontal irradiation (GHI) at the ground level in cloudiness conditions (taking into account scattering and absorption by atmosphere). Almost all the machine learning approaches are based on a stationary hypothesis of the time series. they are usually based on the assumption that the data generation mechanism does not change over time. For the global irradiation forecasting in dividing each measured solar irradiation by its corresponding clear sky value, the periodicity of the solar irradiation time series is then removed; thus, a new time series is created with values between 0 and 1 and the forecasting method is then applied to predict only the cloud occurrences and not the fact that GHI is more important at 12h than at 6h, or in summer than in winter. The CS model chosen is the SOLIS model developed by Mueller [START_REF] Mueller | Rethinking satellite-based solar irradiance modelling: The SOLIS clear-sky module[END_REF] and based on the studies of Ineichen et al. [START_REF] Ineichen | Comparison of eight clear sky broadband models against 16 independent data banks[END_REF][START_REF] Ineichen | A broadband simplified version of the Solis clear sky model[END_REF]; it gives excellent results when it was compared with solar measurements in Europe. It is based on radiative transfer models, Beer Lambert function, and their integration on the solar spectra uses the information's of news satellites ERS-2/ENVISAT. The calculation of the clear sky irradiation CS(t) is given by:

𝐶𝑆(𝑡) = 𝐻 0 . 𝑒 -𝜏 𝑠𝑖𝑛 𝑏 (ℎ(𝑡)) . 𝑠𝑖𝑛(ℎ(𝑡)) (1)
h(t) is the solar height in degrees and H0 is the extraterrestrial irradiation (irradiation on the top of the atmosphere before being scattered and absorbed by the atmosphere [START_REF] Iqbal | An introduction to solar radiation[END_REF])); τ is the global total atmospheric depth and the parameter b is a fitting parameter, both of them depending on the meteorological characteristics of the site, Mueller and Ineichen study the different parameters used in modelling of clear sky [START_REF] Ineichen | A broadband simplified version of the Solis clear sky model[END_REF][START_REF] Mueller | Rethinking satellite-based solar irradiance modelling: The SOLIS clear-sky module[END_REF]. τ and b depend on the site and their values can be found in (https://aeronet.gsfc.nasa.gov/).

The clear sky index (CSI) is calculated by:

𝐶𝑆𝐼(𝑡) = 𝐺𝐻𝐼(𝑡) 𝐶𝑆(𝑡) (2)
All the machine learning tools described in the next section can be only applied to stationary time series. CSI is forecasted from previous CSI values and GHI is then obtained using Eq. 2;

the errors metrics are given in term of GHI in absolute and relative values.

Input matrix dimension of the forecasting tool and k-fold sampling

The next step of the preprocessing is to choose the dimension of the input matrix of the forecasting tools. The data are presented under a time series (TS) formalism, a time series is a series of data points indexed in function of the time (𝐶𝑆𝐼(𝑡)). The approach consists in predicting the future clear sky index (at different time horizons) from past observed data.

Mathematically, the formulation is:

𝐶𝑆𝐼(𝑡 + ℎ) = 𝑓(𝐶𝑆𝐼(𝑡), 𝐶𝑆𝐼(𝑡 -1), … , 𝐶𝑆𝐼(𝑡 -𝑛)) + 𝜖(𝑡 + ℎ) (3)
𝜖(𝑡 + ℎ) is a random white noise and f depending on the model. Sometimes, the model is trained on a given percentage of first data and tested on the complementary percentage of last data. But when this data repartition, there exist the risk that the model was trained on a period for which some specific meteorological phenomena occur, conducing to a non-generalizable model. In view to avoid this problem and to make the training less dependent of particular meteorological periods, the k-fold method is used during the validation phase; it consists in dividing the dataset in k samples (here k=10), each sample is used at least one time for the training and one time for the test, and this process is repeated as many time as necessary. Thus, the results are independent of the set of data used for the training because using only one data set (with its own statistical particularities) can reduce the robustness of the conclusions.

Note that the k-fold cross-validation is only used during the retrospective model's comparison but never in operational mode or prediction in real installation. In a retrospective context, the concept of present, future and past data does not make sense…all the data are part of the past.

What is essential is that the data used for training are never used during the test. The k-fold methodology induces that the training set is not defined "in the past" with respect to the testing set. Testing and training data are interleaved.

Meteorological stations and data

The 2

), the station is located in the mountains, at about 100 km from the Mediterranean Sea and presents often a high nebulosity. The second dataset comes from measurements realized in the laboratory at Ajaccio (Corsica, France, 41°55 N, 8°44 E, 4m asl) (blue circle in Fig. 2) at about 100 m from the Mediterranean Sea. The third dataset is constituted by measurements realized in Tilos island (Greece, 36°24 N, 27°22 E, 96 m asl) (green circle in Fig. 2), Tilos is a small Greek island in the Dodecanese archipelago, the tallest mountain is about 650 m high.

The meteorological conditions (variability of solar radiation) in these three stations are very different, the solar radiation is outlined and will be confirmed in paragraph 5.

The periods of the solar data collect and the number of validated measures (after quality check and night data extraction) are given in Table 1.

The forecasting models

Eleven forecasting models are briefly described and more information on each of them are available in literature. These models are classified in three categories: naïve models, classical machine learning models and regression trees-based models. The symbol ̂ indicates that the value is predicted, without this symbol the value is measured.

Naive models

The two naive models are generally used as a reference in view to compare it with more sophisticated models. Indeed, if a complex model is not more efficient than these two models which are easy to implement and not requiring historical data, their use is not justified. The first model is the simplest one, the persistence, it is the repetition of the measure at the instant t to the predicted value at the instant t+h (h being the forecasting horizon) [START_REF] Diagne | Review of solar irradiance forecasting methods and a proposition for small-scale insular grids[END_REF]:

𝐺𝐻𝐼 ̂(𝑡 + ℎ) = 𝐺𝐻𝐼(𝑡) (4)
𝐺𝐻𝐼 ̂(𝑡) and 𝐺𝐻𝐼(𝑡) are respectively the predicted and measured hourly global horizontal solar irradiation at time t. This model is very simple to implement but it gives generally results with low accuracy.

The daily profile of the solar radiation can be added in using the Solis clear sky model to give the smart persistence, a simple improvement of the previous model developed for solar forecasting and used as reference on comparison with other models [START_REF] Voyant | Statistical parameters as a means to a priori assess the accuracy of solar forecasting models[END_REF] .

𝐺𝐻𝐼 ̂(𝑡 + ℎ) = 𝐺𝐻𝐼(𝑡). 𝐶𝑆(𝑡+ℎ) 𝐶𝑆(𝑡) (5)
The accuracy of such models decreases rapidly with the time horizon and is generally not adequate for a horizon higher than one hour.

CSI(t+h) being the clear sky index at time t+h, φ and θ the ARMA parameters deduced by a least square method, p and q are the model orders and ε(t) is the error or a noise related to a normal distribution.

Artificial neural network (ANN): MultiLayer Perceptron (MLP)

A feed forward MLP, a type of ANN, with one hidden layer and one output layer is used. The utilization of this technique for application on energy systems and in forecasting and modelling solar radiation [START_REF] Kalogirou | Applications of artificial neural-networks for energy systems[END_REF][START_REF] Mellit | Artificial intelligence techniques for modelling and forecasting of solar radiation data: A review[END_REF]. The equation for a MLP with one hidden layer of m neurons, one output neuron and n input variables is given by:

𝐶𝑆𝐼 ̂(𝑡 + ℎ) = ∑ 𝜔 𝑗 𝑚 𝑗=1
. (𝑔(∑ 𝜔 𝑖,𝑗 . 𝐶𝑆𝐼(𝑡 -𝑗)

𝑛-1 𝑖=0 + 𝑏 𝑗 ) (7)
with CSI the input vector of n clear sky indexes, 𝐶𝑆𝐼 ̂(𝑡 + ℎ) the predicted value, 𝑏 𝑗 the biases of the hidden neuron j and 𝜔 𝑖,𝑗 the weights between the input i and the hidden node j, g is the transfer function,𝜔 𝑗 the weight between the output and the hidden neuron j

Gaussian Process (GP)

The Gaussian process, nonlinear model, is a Gaussian distribution with an infinity of variables, [START_REF] Rasmussen | Gaussian Processes in Machine Learning[END_REF]) details the work and implementation of GP. Every predicted CSI is represented by the sum of a function 𝑓(𝐶𝑆𝐼(𝛕)) and an independent Gaussian noise 𝒩(𝑂, 𝜎 𝑛 2 ) with a variance 𝜎 𝑛 2 and CSI(τ) = (CSI(t), CSI(t -1), ⋯ , CSI(t -n)):

𝐶𝑆𝐼 ̂(𝑡 + ℎ) = 𝑓(𝐶𝑆𝐼(𝛕)) + 𝒩(𝑂, 𝜎 𝑛 2 ) (8)
GP is defined as a mean function 𝑚(𝐶𝑆𝐼(𝛕)) and a covariance function k. k is an exponential squared function relying 𝐶𝑆𝐼 ̂(𝑡 𝑝 + ℎ) with 𝐶𝑆𝐼 ̂(𝑡 𝑞 + ℎ) , tp and tq being two successive instants:

𝑘(𝐶𝑆𝐼 ̂(𝑡 𝑝 + ℎ), 𝐶𝑆𝐼 ̂(𝑡 𝑞 + ℎ)) = 𝜎 𝑓 2 𝑒𝑥𝑝 [ -(𝐶𝑆𝐼(𝑡 𝑝 )-𝐶𝑆𝐼(𝑡 𝑞 )) 2 2𝑙 2 ] + 𝛿 𝑝𝑞 𝜎 𝑛 2 (9)
𝛿 𝑝𝑞 is Kronecker delta, 𝜎 𝑓 2 and 𝜎 𝑛 2 are the hyper parameters of the covariance function and are responsible of the complexity of the model, l is a length parameter.

Support Vector Regression (SVR)

SVR is a Kernel based model, firstly developed for regression problems and applied now for forecasting purposes [START_REF] Lauret | A benchmarking of machine learning techniques for solar radiation forecasting in an insular context[END_REF][START_REF] Vapnik | The Nature of Statistical Learning Theory[END_REF]. With a training dataset 𝐷 = {𝐶𝑆𝐼(), 𝐶𝑆𝐼(𝑡 + ℎ)} ; the predicted CSI can be expressed by:

𝐶𝑆𝐼 ̂(𝑡 + ℎ) = ∑ 𝛼 𝜏 𝑡-1 𝜏=1
. 𝑘 𝑟𝑏𝑓 (𝐶𝑆𝐼(𝑡 + ℎ), 𝐶𝑆𝐼(𝑡 -𝜏) ) + 𝑏 [START_REF] Breiman | Random forests[END_REF] and the Kernel radial basis function is:

𝑘 𝑟𝑏𝑓 (𝐶𝑆𝐼(𝑡 𝑝 ), 𝐶𝑆𝐼(𝑡 𝑞 ) ) = 𝑒𝑥𝑝 [ -(𝐶𝑆𝐼(𝑡 𝑝 )-𝐶𝑆𝐼(𝑡 𝑞 )) 2 2𝜎 𝑓 2 ] (11)
αi is the Lagrange multipliers, solutions of a quadratic problem, b is the bias determined by specific conditions.

Regression Trees based methods

Decision trees based on "If-Then" rules are used for classification and produces understandable models due to their graphical representation, they are adapted to work on forecasting problems, like solar radiation forecasting (Aggarwal & Saini, 2014; Burrows, 1997). They have been then extended for predicting numerical values of attributes and led to regression trees which are a decision tree in which the leaf nodes have been set as regression models, and therefore, continuous numeric values can be predicted. Some results are given on forecasting of photovoltaic production system using regression trees (Persson, Bacher, Shiga, & Madsen, 2017).

A. Standard and pruned regression trees (RT and RT-pruned)

Hastie and Tibshirani [START_REF] Hastie | Generalized additive models[END_REF]) proposed a formalization of the RT models:

𝐶𝑆𝐼 ̂(𝑡 + ℎ) = ∑ 𝑘 𝑖 . 𝐼(𝐶𝑆𝐼(𝑡 -𝑖))

𝑡-1 𝑖=1 [START_REF] Jiang | Global horizontal radiation forecast using forward regression on a quadratic kernel support vector machine: Case study of the Tibet Autonomous Region in China[END_REF] with ki constant factors, I a function which return 1 if input is used and 0 elsewhere; The trees are built by splitting the data based on the values of predictive attributes. A regression model is then computed for each node. The pruning aspect of RT is operated with an elevation of the quadratic error tolerance per node. Splitting nodes stops when the quadratic error per node drops below tolerance. For normal RT, the tolerance is close to zero, while for the pruned RT, a higher value is chosen using a heuristic method based on the minimizing of the global error of prediction (Pedro, Coimbra, David, & Lauret, 2018) .

B. Boosted and bagged regression trees (RT-boosted and RT-bagged)

It is interesting in "ensemble learning" methods to generate many regressors and aggregate their results. Two usual methods are boosting and bagging of RT, these two methods are improvement of the classical regression trees.

The boosting method consists in assembly weak RT classifiers and take the average of predictions in order to improve the efficiency [START_REF] Huang | A semi-empirical approach using gradient boosting and k-nearest neighbors regression for GEFCom2014 probabilistic solar power forecasting[END_REF]. A weak predictor is a simple single split RT and the next trees give more weights to the data badly predicted at the previous point, De 'Ath explains in details the method to use boosted RT in modelling and forecasting applications [START_REF] De | Boosted Trees for Ecological Modeling and Prediction[END_REF]. The function for additive models applied to solar forecasting by boosted RT is:

𝐶𝑆𝐼 ̂(𝑡 + ℎ) = ∑ 𝛽 𝑚 𝑚 𝑏(𝐶𝑆𝐼 ̂(𝑡 + ℎ), 𝛾 𝑚 ) (13)
The function 𝑏 represents the individual trees with 𝛾 𝑚 the split variable and 𝛽 𝑚 the weight at each node. The Bagging method (for bootstrap aggregating) is another improvement of the RT [START_REF] Breiman | Bagging predictors[END_REF]; the model is an aggregation of RT which grows from samples of dataset :

𝐶𝑆𝐼 ̂(𝑡 + ℎ) = 𝑎𝑣 𝑘 𝜑 𝑘 (𝐶𝑆𝐼 ̂(𝑡 + ℎ)) (14)
𝜑 𝑘 are the different predictors before the aggregation and 𝑎𝑣 𝑘 the mean of the different predictors.

C. Random forests

Random forests [START_REF] Breiman | Bagging predictors[END_REF] add an additional layer of randomness to bagging. In a random forest, the dataset is equally divided in samples but each regression tree grows differently, each node is split using the best among a subset of predictors randomly chosen at that node [START_REF] Ibrahim | A novel hybrid model for hourly global solar radiation prediction using random forests technique and firefly algorithm[END_REF]. This improvement by randomness gives robustness to the model and decreases the over-training risks.

Evaluation of models performance

To evaluate the accuracy of the models, three different statistics are used. The mean absolute error (MAE), defines the absolute value of the gap between the observed and the predicted value.

𝑀𝐴𝐸 = ∑ |𝐺𝐻𝐼 ̂(𝑖)-𝐺𝐻𝐼(𝑖)| 𝑁 𝑖=1
𝑁 [START_REF] Hornik | Multilayer feedforward networks are universal approximators[END_REF] with 𝐺𝐻𝐼 ̂(𝑖) the predicted variable, 𝐺𝐻𝐼 (𝑖) the observed variable and N the number of data.

The normalized root mean squared error (nRMSE) is more sensitive to large forecast errors, and hence is suitable for applications where small errors are more tolerable and larger errors cause disproportionately high costs, as for example in the case of utility applications. It is probably the reliability factor that is most appreciated and used: it's a good statistical index to evaluate the accuracy of a models, the aim of an operator is to minimize it in order to improve model performances:

𝑛𝑅𝑀𝑆𝐸 = √ 1 𝑁 ∑ (𝐺𝐻𝐼 ̂(𝑖)-𝐺𝐻𝐼(𝑖)) 2 𝑁 𝑖=1
𝐺𝐻𝐼 ̅̅̅̅̅̅ [START_REF] Paoli | Forecasting of preprocessed daily solar radiation time series using neural networks[END_REF] where 𝐺𝐻𝐼 is the algebraic average of the observed values.

The skill score is an index calculated in order to compare the performance of a given model with a reference model, here the reference model is the smart persistence model (SP) described in 3.1.

𝑆𝑘𝑖𝑙𝑙 𝑆𝑐𝑜𝑟𝑒 = 𝑀𝑒𝑡𝑟𝑖𝑐 𝑓𝑜𝑟𝑒𝑐𝑎𝑠𝑡 -𝑀𝑒𝑡𝑟𝑖𝑐 𝑟𝑒𝑓𝑒𝑟𝑒𝑛𝑐𝑒 𝑀𝑒𝑡𝑟𝑖𝑐 𝑝𝑒𝑟𝑓𝑒𝑐𝑡_𝑓𝑜𝑟𝑒𝑐𝑎𝑠𝑡 -𝑀𝑒𝑡𝑟𝑖𝑐 𝑟𝑒𝑓𝑒𝑟𝑒𝑛𝑐𝑒 = 1 - 𝑛𝑅𝑀𝑆𝐸 𝑓𝑜𝑟𝑒𝑐𝑎𝑠𝑡 𝑛𝑅𝑀𝑆𝐸 𝑟𝑒𝑓𝑒𝑟𝑒𝑛𝑐𝑒 (17)
The skill score is always inferior at 1, negative if the forecaster is less performant than the reference, 0 if the performances are similar, and positive if it is better.

Evaluation of the variability

The evaluation of the variability of different datasets is made thanks to a statistical parameter, Voyant et al. [START_REF] Voyant | Statistical parameters as a means to a priori assess the accuracy of solar forecasting models[END_REF] described several parameters to quantify the variability of a solar irradiation time series and the mean absolute log return appeared to be the most efficient:

𝑚𝑒𝑎𝑛𝑎𝑏𝑠 (𝑙𝑜𝑔𝑟) = ∑ |log(𝐶𝑆𝐼(𝑖))-log(𝐶𝑆𝐼(𝑖-1))| 𝑁 𝑖=1 𝑁 (18)
Other variability metrics are available like the P parameter described by Perez et al. (Perez, Kivalov, Schlemmer, Hemker Jr., & Hoff, 2012) and computed from the standard deviation (std) by 𝑃 = 𝑠𝑡𝑑(𝐶𝑆𝐼(𝑡) -𝐶𝑆𝐼(𝑡 -1)). In order to not overburden the manuscript, the variability is only estimated with the mean absolute log return. The variability of the three solar dataset were computed using Eq [START_REF] Mueller | Rethinking satellite-based solar irradiance modelling: The SOLIS clear-sky module[END_REF] and the Mean absolute log return values for each site are shown in Table 2.

Table 2. Variability for three CSI datasets according to the mean absolute log return.

Site Odeillo Tilos Ajaccio

Mean absolute log return 0.5028 0.3732 0.1961

Type of variability Strong Medium Weak

A large difference between the three datasets is noted, the variability of the solar irradiance measurements differs from one site to another; thus, the performance of the forecasting methods presented here will have a more universal character in nature because they are the result of calculations realized from various weather conditions. Moreover, it will be possible to draw a correlation between the site variability and the ranking of the models in term of performance.

Solar irradiance forecasting and performances of forecasting models

This section presents the results of solar irradiance forecasting for the three sites. The preprocessing is the same for every dataset:

-Calculation of the clear-sky irradiation by the Solis model; -Calculation of clear sky indexes;

-Remove of night hours and hour with the sun elevation up to 10°.

The auto-mutual information which is a part of information theory, is a way to select the number of input given in the models [START_REF] Fu | Uncertainty analysis of an integrated energy system based on information theory[END_REF][START_REF] Parviz | Mutual Information Based Input Variable Selection Algorithm and Wavelet Neural Network for Time Series Prediction[END_REF]. For each time series shows that the optimal number of input data n (first minimum of the auto mutual information criteria, see paragraph 2.2) is n=8 for Ajaccio and Odeillo and n=6 for Tilos. Thus, for estimating the future solar irradiation, the n previous measured solar irradiations are used. The training set is about 80% of the dataset and the testing set about 20% considering a k-fold sampling equal to 10.

Tables 3 and4 show the results of the nRMSE and MAE calculations for the eleven models and a forecast horizon from 1 to 6 hours with an hourly resolution. Conclusions are difficult to drawn for these three locations. More meteorological sites would have been welcome but it is difficult to find reliable and long hourly solar irradiation data set and to propose a synthetic paper with so many forecasting models applied to a larger number of locations. The results generalization is, as often, very complicated and depends also on the error metrics used.

However, some first conclusions can be drawn, the first one the naïve models (persistence and smart persistence) give always bad performances for all the sites and are the worst models. For a high variability (Odeillo), the best performances are obtained with ensemble learning models (bagged regression trees and random forest), the machine learning models are less performant.

These results show that ensemble learning models have the capacity to apprehend complicated phenomena. The high dataset variability explains the bad performances of "classical" models as MLP, ARMA or simple regression trees.

For a medium variability dataset (Tilos), the best results come from the ARMA model, but occasionally smart persistence and random forest gave good results. The poorer performances of machine learning models, for Tilos, is probably due to the small length of the dataset (smaller than for Ajaccio and Odeillo).

For a low variability dataset (Ajaccio) the best models are ARMA and MLP, followed by smart persistence and bagged regression trees; however, the errors are of the same order of magnitude for all the models excepted for persistence and simple regression tree. With such a variability, machine learning models are strong predictors.

The results from the MAE point of view are similar:

-For a high variability dataset, the ensemble learning models present the smallest absolute error.

-For a medium variability and the two first hours of forecasting horizon, the statistical models give the best performances and for deeper horizons, the ensemble learning models are the best predictors -For a low variability dataset, excepted for the naïve model, all the models are similar in term of absolute error, the results are of the same order of magnitude.

The skill score (related to the smart persistence reference model) is computed for each model, The skill score compares the performances of the model with the smart persistence. If the performances of the models are generally better than the smart persistence for Odeillo, this fact is not checked for the two other sites. Thus, as expected, the improvement on the forecasting due to the use of machine learning methods compared with naïve models increases with the variability of the dataset to predict.

Comments concerning prediction error and variability

Depending on the variability of the site, the models inducing the lowest prediction error are different. In Table 5, the main conclusions concerning the benchmarking study of forecasting model are drawn for the three sites: Ajaccio, Odeillo and Tilos.

For a low variability, (Ajaccio), a linear model ARMA and a classical MLP give very good results. With a lower reliability, the smart persistence or a RT with the bagged mode for the horizons inferior to 4 hours can be also used. For a high variability site (here Odeillo), the It is obvious than these results must be confirmed on other meteorological sites with a wide range of variability. Gaussian process

Conclusion

Eleven statistical and machine learning tools for global solar irradiation forecasting were analyzed and compared in term off performances on three sites with different meteorological characteristics. To characterize the solar data time series measured in each location, an evaluation of the variability was realized. A ranking of the forecasting methods according to the prediction horizon and depending on the site was realized; the main conclusions can be drawn:

-For Ajaccio, with a weak variability, ARMA and MLP are the most efficient tools;

-For Tilos, with a medium variability, ARMA and bagged regression tree are recommended; -For Odeillo with a higher variability, the forecasting reliability is lower, but the best results were obtained for the bagged regression tree and the random forest approach.

Even as the results are not totally conclusive, it appears that that higher is the variability, more complex is the forecasting tool to be used. In general, when the cloud occurrences are low, the utilization of an ARMA model is sufficient else the regression tree based on the bagging mode is the most reliable. The conclusions drawn here concern only three sites and to confirm them, other simulations must be carried out on well-chosen sites with particular weather characteristics (desert, monsoon, extreme longitude, etc.).

Fig 1

 1 Fig 1 shows the different steps (detailed afterwards) of this methodology.

Figure 1 .

 1 Figure 1. Flowchart of the model development before forecasting

  The future time step (t+h) is forecasted from a given number of observed data at previous times (t, t-1…, t-n) which must be optimally determined. In other words, the objective is to calculate the value of n (number of previous values) and to obtain 𝜖 as low as possible. The choice of n, i.e. the dimension of the input matrix, is made by an auto mutual information method, in (Luo et al., 2017; Parviz, Nasser, & Motlagh, 2008) some details concerning this method. This auto mutual information is a property of a time series and depends on each dataset. It determinates the degree of statistical dependence of the variables. The lag corresponding to the first minimum of this parameter corresponds to the best n to consider. In machine learning method, the times series is divided into two groups: a training and a testing group. The first one is used to train the model i.e. inputs and output are given to the model and a training algorithm determines the parameter values of the model; then, once these parameters known, the model is tested on a second set of data called testing set (never used for the training) and the error metrics are calculated only on this testing set.

  datasets are time series of horizontal global solar irradiation measurements (GHI) in three different sites (Fig 2).

Figure 2 .

 2 Figure 2. Geographical situation of experimental sites

3. 2 . Classical machine learning models 3 . 2 . 1 .

 2321 Auto Regressive Mobile Average (ARMA) The ARMA model, used in energy consumption forecasting (de Oliveira & Cyrino Oliveira, 2018) includes two parts, an auto regressive one and a mobile average one. This model predicts the future values, as described in Ref (De Gooijer & Hyndman, 2006; Faraday & Chatfield, 1998), from a linear combination of past values and a past residue: 𝐶𝑆𝐼 ̂(𝑡 + ℎ) = 𝜀(𝑡) + ∑ 𝜑 𝑖 𝑝 𝑖=0. 𝐶𝑆𝐼(𝑡 -𝑖) + ∑ 𝜃 𝑖 . 𝜀(𝑡 -𝑖) 𝑞 𝑖=0

  each horizon and site and is plotted in Fig 3.

Figure 3 .

 3 Figure 3. Skill score vs forecast horizon for a) Odeillo, b) Tilos and c) Ajaccio

  complicated, ensemble machine learning methods like random forest or bagged RT are more recommended and linear models or smart persistence must be rejected . For an intermediate variability (here Tilos), the two best models are ARMA and Bagged RT, but all the other machine learning models have about the same performances and the model ranking is very difficult.

Table 1 .

 1 Period of availability and number of validated data for the three sites.

	Location	Period of collected data	Number of validated data

Table 3 .

 3 nRMSE vs forecast horizon, the two best models are highlighted.

		Horizon	h+1	h+2	h+3	h+4	h+5	h+6
	Odeillo (variability: 0.5028)	Persistence 40.91% 61.81% Smart persistence 37.00% 54.59% ARMA 34.42% 47.12% MLP 29.88% 43.51% Regression tree (RT) 36.73% 51.94% Boosted RT 30.20% 43.24% Bagged RT 28.80% 42.00% Pruned RT 29.90% 43.97% Random forest 28.76% 42.75% Gaussian process 28.65% 42.15%	76.18% 66.68% 49.26% 45.86% 54.95% 45.59% 44.66% 46.47% 44.89% 45.37%	84.53% 74.64% 49.76% 47.87% 56.76% 47.71% 46.63% 48.47% 46.56% 46.67%	88.35% 77.77% 49.16% 48.44% 57.86% 48.57% 47.83% 50.01% 47.78% 48.54%	86.25% 77.76% 48.56% 48.72% 57.45% 48.78% 47.52% 49.84% 48.34% 48.42%
		Support vector regression 34.18% 57.71%	57.25%	57.76%	58.27%	57.25%

Table 4 .

 4 MAE values (in Wh/m²) vs forecast horizon, the two best models are highlighted

		Horizon	h+1	h+2	h+3	h+4	h+5	h+6
	Odeillo (variability: 0.5028)	Persistence 147.07 Smart persistence 124.62 ARMA 130.73 MLP 105.87 Regression tree (RT) 119.94 Boosted RT 107.77 Bagged RT 98.94 Pruned RT 105.68	238.13 199.81 179.45 162.23 186.50 161.64 156.32 163.08	297.75 249.00 189.20 172.95 201.09 171.97 166.26 173.14	328.68 341.65 277.70 288.13 190.79 188.41 181.43 185.81 208.70 212.29 180.56 185.47 175.28 179.78 181.84 188.12	331.23 285.29 187.33 186.29 210.24 185.90 179.82 186.80
		Random forest	97.48	156.73	167.45	174.15 180.13	181.94

Table 5 .

 5 Ranking of the prediction methodologies according to the variability

		Weak variability	Medium	Variability	Strong variability
		𝒎𝒆𝒂𝒏𝒂𝒃𝒔 (𝒍𝒐𝒈𝒓) < 0.2	0.2< 𝒎𝒆𝒂𝒏𝒂𝒃𝒔 (𝒍𝒐𝒈𝒓) < 0.4	0.4 < 𝒎𝒆𝒂𝒏𝒂𝒃𝒔 (𝒍𝒐𝒈𝒓)
	Recommended models	ARMA, MLP	ARMA, Bagged RT	Bagged RT, Random forest
	Usable models	SP (horizon < 4h), Bagged RT	Random forest, Gaussian process, MLP, SP (horizon <3h)

Acknowledgment

This work is supported by European Union's Horizon 2020 research and innovation programme under grant agreement No 646529 through the TILOS project (Technology innovation for the Local Scale, Optimum integration of Battery Energy Storage). The data sets were providing by the CNRS UPR PROMES 8521 laboratory in Odeillo (Pyrénées Orientales, France), the CNRS UMR SPE 6134 Laboratory in Ajaccio (Corsica, France) and the local weather station on TILOS Island.


Daphné Trama
email: daphne.trama@cea.fr

Pierre-Emmanuel Clet
email: pierre-emmanuel.clet@cea.fr

Aymen Boudguiga
email: aymen.boudguiga@cea.fr

Renaud Sirdey
email: renaud.sirdey@cea.fr

A Homomorphic AES Evaluation in Less than 30 Seconds by Means of TFHE

Keywords: Security and privacy → Public key (asymmetric) techniques AES, Fully Homomorphic Encryption, Transciphering, TFHE, Programmable Bootstrapping

Since the pioneering work of Gentry, Halevi, and Smart in 2012 [18], the state of the art on transciphering has moved away from work on AES to focus on new symmetric algorithms that are better suited for a homomorphic execution. Yet, with recent advances in homomorphic cryptosystems, the question arises as to where we stand today. Especially since AES execution is the application that may be chosen by NIST in the FHE part of its future call for threshold encryption. In this paper, we propose an AES implementation using TFHE programmable bootstrapping, which runs in less than a minute on an average laptop. We detail the transformations carried out on the original AES code as well as the optimized FHE operators we developed to lead to a more efficient homomorphic evaluation. We also duly give several execution times on different machines, depending on the type of execution (sequential or parallelized). These times vary from 4.5 minutes (resp. 54 secs) for sequential (resp. parallel) execution on a standard laptop down to 28 seconds for a parallelized execution over 16 threads on a multi-core workstation.

INTRODUCTION

With recent advances in FHE, is a homomorphic AES still as impractical as it was ten years ago? The work of Gentry, Halevi, and Smart [START_REF] Gentry | Homomorphic Evaluation of the AES Circuit[END_REF] in 2012 pushed research towards new symmetric cryptosystems designed primarily to be faster to evaluate over FHE. Indeed, they performed an AES-128 homomorphic evaluation with BGV using HElib, with now obsolete parameters that did not allow bootstrapping. They then obtained an execution time of 4.1 minutes, but without allowing further operations on the final ciphertext. With bootstrapping, thus allowing further calculations after the homomorphic execution of the AES, their runtime grew to 17.5 minutes. So neither of these two approaches could be used in practice. Since transciphering (the ability to homomorphically turn low overhead symmetrically encrypted data into homomorphically encrypted ones) is an important issue for FHE practicality, several teams then decided to create new symmetric cryptosystems, whose encryption operations were specifically chosen to be more rapidly executed in the homomorphic domain. As of today, there are many proposals, from block ciphers (LowMC [START_REF] Albrecht | Ciphers for MPC and FHE[END_REF], PRINCE [START_REF] Borghoff | PRINCE -A Low-Latency Block Cipher for Pervasive Computing Applications -Extended Abstract[END_REF], CHAGHRI [START_REF] Ashur | Chaghri -A FHE-Friendly Block Cipher[END_REF]) to stream ciphers (Elisabeth [START_REF] Cosseron | Towards Globally Optimized Hybrid Homomorphic Encryption -Featuring the Elisabeth Stream Cipher[END_REF], PASTA [START_REF] Dobraunig | Pasta: A Case for Hybrid Homomorphic Encryption[END_REF], Kreyvium [START_REF] Canteaut | Stream ciphers: A Practical Solution for Efficient Homomorphic-Ciphertext Compression[END_REF]). Each comes with its pros and cons. For instance, PRINCE [START_REF] Borghoff | PRINCE -A Low-Latency Block Cipher for Pervasive Computing Applications -Extended Abstract[END_REF] is a block cipher especially created to be lightweight and, although it was initially proposed independently of Gentry's breakthrough on FHE, has a number of desirable properties with respect to homomorphic execution: a moderate number of rounds, small depth (for a block-cipher) and a low gate count/footprint of the decryption and encryption functions. PRINCE was one of the first symmetric algorithms for which an FHE execution attempt was done [START_REF] Lepoint | A Comparison of the Homomorphic Encryption Schemes FV and YASHE[END_REF]. LowMC [START_REF] Albrecht | Ciphers for MPC and FHE[END_REF], on the other hand, is the first block cipher explicitly designed with FHE and MPC in mind. Although demonstrating competitive FHE execution performances at the time of proposal, its design was intrinsically bit-oriented while the FHE state of the art has moved away from bit-level FHE operations due to the relative inefficiency of this latter approach. In 2022 Ashur et al. presented CHAGHRI [START_REF] Ashur | Chaghri -A FHE-Friendly Block Cipher[END_REF], an FHE-friendly block cipher enabling efficient transciphering in BGV-like schemes. A complete CHAGHRI circuit can be implemented using 16 multiplications, 48 Frobenius automorphisms, and 32 rotations. The authors implemented it with HElib in order to compare it with Gentry et al. work on AES. Although their implementation is claimed to be 63% faster than [START_REF] Gentry | Homomorphic Evaluation of the AES Circuit[END_REF], an attack on CHAGHRI has recently been proposed [START_REF] Liu | Coefficient Grouping: Breaking Chaghri and More[END_REF]. Also introduced in 2022, Elisabeth is a family of stream ciphers especially designed to be efficient for Hybrid Homomorphic Encryption (HHE). The authors use TFHE and propose a Rust implementation (using the Concrete library) of Elisabeth-4, that is to say, a cryptosystem in which inputs are on 4 bits. So it would take 32 executions of the cipher to obtain a 128-bits ciphertext. Before Elisabeth, the PASTA cryptosystem, implemented with BGV/BFV proposed an optimized cipher for integer HHE use cases. They also benchmark several HHE schemes, using the HElib library. But the use of a non-bootstrapping-able scheme limits the number of operations to be further performed on the ciphertexts. Kreyvium is a stream cipher, which is a variant of Trivium [START_REF] De | Trivium. In New Stream Cipher Designs -The eSTREAM Finalists[END_REF] (a stream cipher belonging to the eSTREAM portfolio). The main motivation for introducing Kreyvium was to propose an FHE-friendly symmetric primitive with 128-bits of security, based on the sound design rationals of Trivium. Additionally, the state of the art also includes homomorphic evaluation of several variants of the Grain-128 stream cipher by means of TFHE either in gate-bootstrapping mode or exploiting its functional bootstrapping capabilities [START_REF] Bendoukha | Optimized stream-cipher-based transciphering by means of functional bootstrapping[END_REF][START_REF] Bendoukha | Revisiting Stream-Cipher-Based Homomorphic Transciphering In The TFHE Era[END_REF]. Often compared to the homomorphic execution times of AES as a guarantee of efficiency, none of these cryptosystems has been standardized (with the notable exception of Grain-128, which was a finalist in the recent NIST competition on lightweight cryptography). Yet, an "efficient-by-FHE-standards" homomorphic AES execution remains interesting for the research community working on transciphering, even if it does not bring any revolution. This is especially so, since AES execution may be the application chosen by NIST in the FHE part of its future call for proposals on threshold encryption 1 .

Contribution-In this paper, we propose an AES implementation using TFHE programmable bootstrapping, which runs in less than a minute on a standard laptop PC. We first detail the modifications carried out on the original AES code as well as the optimized FHE operators we developed to lead to an efficient homomorphic evaluation of the algorithm. Then we give details about the benchmark made to determine which decomposition basis to use to have a faster evaluation of the algorithm. We finally provide experimental execution times on different machines, depending on the type of execution (sequential or parallelized).

Paper Ogranization-This paper is organized as follows: Section 2 reviews the basics of the TFHE cryptosystem and gives the necessary details of the tree-based method for bootstrapping with multi-input ciphertexts and its optimization with multi-value bootstrapping. Section 3 gives a brief reminder on the AES. Section 4 provides a detailed exposition of our approaches to transform the original AES code and implement the most optimized and efficient version of it with TFHE programmable bootstrapping. Section 5 presents the performances and results of our methods.

TFHE PRELIMINARIES 2.1 Notations

Let T = R/Z be the real torus, that is to say, the additive group of real numbers modulo 1 (R mod 1). We will denote by T 𝑁 [𝑋 ] 𝑛 the set of vectors of size 𝑛 whose coefficients are polynomials of T [𝑋 ] mod (𝑋 𝑁 + 1). 𝑁 is usually a power of 2. Let B = {0, 1}. ⟨ , ⟩ denotes the inner product.

1 https://csrc.nist.gov/Projects/threshold-cryptography.

TFHE Scheme

The TFHE scheme is a fully homomorphic encryption scheme introduced in 2016 in [START_REF] Chillotti | TFHE: Fast Fully Homomorphic Encryption Library[END_REF] and implemented as the TFHE library2 . TFHE defines three structures to encrypt plaintexts, which we summarize below as fresh encryptions of 0:

• TLWE sample: A pair (𝑎, 𝑏) ∈ T 𝑛+1 , where 𝑎 is uniformly sampled from T 𝑛 and 𝑏 = ⟨𝑎, 𝑠⟩ + 𝑒. The secret key 𝑠 is uniformly sampled from B 𝑛 , and the error 𝑒 ∈ T is sampled from a Gaussian distribution with mean 0 and standard deviation 𝜎. Let M denote the discrete message space (M ∈ T 𝑁 [𝑋 ] or M ∈ T). To encrypt a message 𝑚 ∈ M, we add what is called a noiseless trivial ciphertext (0, 𝑚) to a fresh encryption of 0. We denote by 𝑐 = (𝑎, 𝑏) + (0, 𝑚) = (𝑎, 𝑏 + 𝑚) ∈ T(R)LWE 𝑠 (𝑚) the T(R)LWE encryption of 𝑚 with key 𝑠. A message 𝑚 ∈ Z[𝑋 ] can also be encrypted in TRGSW samples by adding 𝑚 • 𝐻 to a TRGSW sample of 0, where 𝐻 is a gadget decomposition matrix. As we will not explicitly need such an operation in this paper, more details about TRGSW can be found in [START_REF] Chillotti | TFHE: Fast Fully Homomorphic Encryption Library[END_REF].

To decrypt a ciphertext 𝑐, we first calculate its phase 𝜙 (𝑐) = 𝑏 -⟨𝑎, 𝑠⟩ = 𝑚 + 𝑒. Then, we need to remove the error, which is achieved by rounding the phase to the nearest valid value in M. This procedure fails if the error exceeds half the distance between two elements of M.

TFHE Bootstrapping

Bootstrapping is the operation that reduces the noise of a ciphertext, thus allowing further homomorphic calculations. It relies on three basic operations, which we briefly review in this section.

• BlindRotate: rotates a polynomial encrypted as a TRLWE ciphertext by an encrypted index (under the form of a TLWE encryption). It takes several inputs: a ciphertext 𝑐 ∈ TRLWE 𝑘 (𝑚), a vector (𝑎 It is important to note that, during a BlindRotate operation, an excessive noise level in the input TLWE ciphertext (the encrypted index which we use to rotate the polynomial) can lead to errors in the bootstrapping output resulting in incorrect ciphertexts (i.e., ciphertext which does not decrypt to correct calculation results). This has implications for parameters and data representation choices (number of digits and basis).

Algorithm 1 shows the TFHE Gate Bootstrapping [START_REF] Chillotti | TFHE: Fast Fully Homomorphic Encryption Library[END_REF], which aims to evaluate a binary gate operation homomorphically and reduce the output ciphertext noise at the same time. To that end, 0 and 1 are respectively encoded as 0 and 1 2 over T. The first step of this algorithm consists of selecting a value m ∈ T, which will be used afterward to compute the coefficients of the polynomial, which will rotate during the BlindRotate. We call this polynomial 𝑡𝑒𝑠𝑡𝑣 as seen in Step 3. Note that for any 𝑝 ∈ ⟦0, 2𝑁 ⟧ (where ⟦0, 2𝑁 ⟧ corresponds to the set of integers {0,

• • • , 2𝑁 }), the constant term of 𝑡𝑒𝑠𝑡𝑣 • 𝑋 𝑝 is m if 𝑝 ∈⟧ 𝑁 2 , 3𝑁 2 ⟧ and -m otherwise.
Step 4 returns an accumulator 𝐴𝐶𝐶 ∈ TRLWE 𝑠 ′ (𝑡𝑒𝑠𝑡𝑣 • 𝑋 ⟨𝑎,𝑠 ⟩-𝑏 ). Indeed, the constant term of 𝐴𝐶𝐶 ism if 𝑐 is an encryption of 0 and m if 𝑐 is an encryption of 1 2 . Then step 5 creates a new ciphertext 𝑐 by extracting the constant term in position 0 from 𝐴𝐶𝐶 and adding (0, m). Thus, 𝑐 corresponds to an encryption of 0 if 𝑐 is an encryption of 0 and 𝑚 otherwise. On the other hand, if 𝑐 is an encryption of 1 2 and if we choose 𝑚 = 1 2 , the algorithm returns a fresh ciphertext of 1 2 , that is to say the encoding of 1.

In Fig. 1, we present an example of TFHE gate bootstrapping algorithm with Z 4 = {0, 1, 2, 3} as input space. The outer circle in Fig. 1 corresponds to the plaintext encoding in T as {0, 1 4 , 2 4 , 3 4 }. Meanwhile, the inner circle sets the coefficients of the test polynomial 𝑡𝑒𝑠𝑡𝑣 to 1, i.e., m = 1 4 . Then, we rotate the test polynomial during the bootstrapping by the phase 𝜙 (𝑐 0 ) of the input ciphertext 𝑐 0 . In our example, we obtain as bootstrapping output either an encryption of the encoding of 1 for positive inputs {0, 1 4 }, or an encryption of the encoding of -1 for negative inputs { 2 4 , 3 4 }.

Algorithm 1 TFHE gate bootstrapping [START_REF] Chillotti | TFHE: Fast Fully Homomorphic Encryption Library[END_REF] Require:

a constant 𝑚 ∈ T, a TLWE sample 𝑐 = (𝑎, 𝑏) ∈ TLWE 𝑠 (𝑥 • 1 2 ) with 𝑥 ∈ B, a bootstrapping key 𝐵𝐾 𝑠→𝑠 ′ = (𝐵𝐾 𝑖 ∈ TRGSW 𝑆 ′ (𝑠 𝑖 )) 𝑖 ∈⟦1,
= ⌊2𝑁 𝑎 𝑖 ⌉ ∈ Z, ∀𝑖 ∈ ⟦1, 𝑛⟧ 3: Let 𝑡𝑒𝑠𝑡𝑣 := (1 + 𝑋 + • • • + 𝑋 𝑁 -1 ) • 𝑋 𝑁 2 • m ∈ T 𝑁 [𝑋 ]
4: 𝐴𝐶𝐶 ← BlindRotate( (0, 𝑡𝑒𝑠𝑡 𝑣), (𝑎 1 , . . . , 𝑎 𝑛 , 𝑏 ), (𝐵𝐾 1 , . . . , 𝐵𝐾 𝑛 ) )

5: 𝑐 = (0, m) + SampleExtract(𝐴𝐶𝐶) 6: return KeySwitch 𝑠 ′ →𝑠 (𝑐)

TFHE Functional Bootstrapping

We can use Look-Up Tables to compute functions during the bootstrapping operation. To do so, we replace the coefficients of the test polynomial 𝑡𝑒𝑠𝑡𝑣 with the corresponding values of the LUT. Let us assume that we want to evaluate the function 𝑓 T via a LUT. Then, if we retrieve the 𝑖 𝑡ℎ coefficient of 𝑡𝑒𝑠𝑡𝑣, we actually get 𝑓 T (𝑚 𝑖 ) where 𝑚 𝑖 is the encrypted input to the bootstrapping. We refer to this idea by programmable or functional bootstrapping [START_REF] Chillotti | Programmable Bootstrapping Enables Efficient Homomorphic Inference of Deep Neural Networks[END_REF][START_REF] Chillotti | Improved Programmable Bootstrapping with Larger Precision and Efficient Arithmetic Circuits for TFHE[END_REF][START_REF] Clet | Putting up the swiss army knife of homomorphic calculations by means of TFHE functional bootstrapping[END_REF][START_REF] Kluczniak | FDFB: Full Domain Functional Bootstrapping Towards Practical Fully Homomorphic Encryption[END_REF][START_REF] Yang | TOTA: Fully Homomorphic Encryption with Smaller Parameters and Stronger Security[END_REF].

In Fig. 1, we give an example of functional bootstrapping with Z 4 = {0, 1, 2, 3} as input space. We encode the images of {0, 1 4 } by 𝑓 T as coefficients of the test polynomial (in the inner circle). Meanwhile, we deduce the images of { 2 4 , 3 4 } by negacyclicity. Indeed, in T, we can encode negacyclic functions, i.e., antiperiodic functions with period 1 2 (verifying

𝑓 T (𝑥) = -𝑓 T (𝑥 + 1 2 ))
, where [0, 0.5[ corresponds to positive values and [0.5, 1[ to negative ones. In our example, if we encrypt one of the following values {0, 1 4 , 2 4 , 3 4 } and we give it as input to the functional bootstrapping algorithm, we get {𝑓 T (0), 𝑓 T (1), -𝑓 T (0), -𝑓 T (1)}, respectively.

Almost all of the functional bootstrapping methods from state of the art ( [START_REF] Chillotti | Programmable Bootstrapping Enables Efficient Homomorphic Inference of Deep Neural Networks[END_REF][START_REF] Chillotti | Improved Programmable Bootstrapping with Larger Precision and Efficient Arithmetic Circuits for TFHE[END_REF][START_REF] Clet | Putting up the swiss army knife of homomorphic calculations by means of TFHE functional bootstrapping[END_REF][START_REF] Kluczniak | FDFB: Full Domain Functional Bootstrapping Towards Practical Fully Homomorphic Encryption[END_REF][START_REF] Yang | TOTA: Fully Homomorphic Encryption with Smaller Parameters and Stronger Security[END_REF]) take as input a single ciphertext. In 2021, Guimarães et al., [START_REF] Guimarães | Revisiting the functional bootstrap in TFHE[END_REF] discussed two methods for performing functional bootstrapping with larger plaintexts. They combine several bootstrappings with different encrypted inputs by using a tree or a chain structure. The ciphertexts are encryptions of digits that come from the decomposition of plaintexts in a certain basis 𝐵.

Tree-based Method

Let 𝐵, 𝐵 ′ , 𝑑 ∈ N * and 𝑚 be an integer message. 𝐵 and 𝐵 ′ are the basis on which to decompose the message. We then have 𝑚 = 𝑑 -1 𝑖=0 𝑚 𝑖 𝐵 𝑖 , with 𝑚 𝑖 ∈ ⟦0, 𝐵 -1⟧. From this decomposition, we obtain 𝑑 TLWE encryptions

(𝑐 0 , 𝑐 1 , • • • , 𝑐 𝑑 -1 ) of (𝑚 0 , 𝑚 1 , • • • , 𝑚 𝑑 -1
) on half of the torus T. We denote 𝑓 : ⟦0, 𝐵 -1⟧ 𝑑 → ⟦0, 𝐵 ′ -1⟧ the target function and define 𝑔 as the following bijection:

𝑔 : ⟦0, 𝐵 -1⟧ 𝑑 → ⟦0, 𝐵 𝑑 -1⟧ (𝑎 0 , 𝑎 1 , • • • , 𝑎 𝑑 -1 ) ↦ → 𝑑 -1 𝑖=0 𝑎 𝑖 • 𝐵 𝑖
We then encode a LUT for 𝑓 under the form of 𝐵 𝑑 -1 TRLWE ciphertexts. Each of these ciphertexts encodes a polynomial 𝑃 𝑖 such that:

𝑃 𝑖 (𝑋 ) = 𝐵-1 ∑ 𝑗=0 𝑁 𝐵 -1 ∑ 𝑘=0 𝑓 • 𝑔 -1 ( 𝑗 • 𝐵 𝑑 -1 + 𝑖) • 𝑋 𝑗 • 𝑁 𝐵 +𝑘
Then, we apply the BlindRotateAndExtract (the BlindRotate directly followed by the SampleExtract in position 0) to each test polynomial 𝑡𝑒𝑠𝑡𝑣 = TRLWE(𝑃 𝑖 ) with 𝑐 0 as a selector. We obtain 𝐵 𝑑 -1 TLWE ciphertexts, each corresponding to the encryption of

𝑓 • 𝑔 -1 (𝑚 𝑑 -1 • 𝐵 𝑑 -1 + 𝑖), for 𝑖 ∈ ⟦0, 𝐵 𝑑 -1 -1⟧.
Finally, we use the KeySwitch operation from TLWE to TRLWE to gather them into 𝐵 𝑑 -2 encrypted TRLWE, corresponding to the LUT of ℎ, with:

ℎ : ⟦0, 𝐵 -1⟧ 𝑑 -1 → ⟦0, 𝐵 ′ -1⟧ (𝑎 0 , 𝑎 1 , • • • , 𝑎 𝑑 -2 ) ↦ → 𝑓 (𝑎 0 , 𝑎 1 , • • • , 𝑎 𝑑 -2 , 𝑚 𝑑 -1 )
We then repeat this operation, using the ciphertext 𝑐 𝑖 at step 𝑖, until we obtain a single TLWE ciphertext of 𝑓 (𝑚 0 , 𝑚 1 , • • • , 𝑚 𝑑 -1 ). Note Meanwhile, the inner circles represent the coefficients of the test polynomial 𝑡𝑒𝑠𝑡𝑣. One of these coefficients is extracted as the output of the bootstrapping after the BlindRotate. that the tree-based method must be run independently as many times as the number of digits in the output.

Multi-Value Bootstrapping

Multi-Value Bootstrapping (MVB) [START_REF] Carpov | New techniques for Multi-value input Homomorphic Evaluation and Applications[END_REF] refers to the method for evaluating 𝑘 different LUTs on a single input with a single bootstrapping. MVB factors the test polynomial 𝑃 𝑓 𝑖 associated with the function 𝑓 𝑖 into a product of two polynomials 𝑣 0 and 𝑣 𝑖 , where 𝑣 0 is a common factor to all 𝑃 𝑓 𝑖 . In practice, we have:

(1 + 𝑋 + • • • + 𝑋 𝑁 -1 ) • (1 -𝑋 ) ≡ 2 mod (𝑋 𝑁 + 1)
Now by writing 𝑃 𝑓 𝑖 in the form 𝑃 𝑓 𝑖 = 𝑁 -1 𝑗=0 𝛼 𝑖,𝑗 𝑋 𝑗 with 𝛼 𝑖,𝑗 ∈ Z, we get from the previous equation:

𝑃 𝑓 𝑖 = 1 2 • (1 + 𝑋 + • • • + 𝑋 𝑁 -1 ) • (1 -𝑋 ) • 𝑃 𝑓 𝑖 mod (𝑋 𝑁 + 1) = 𝑣 0 • 𝑣 𝑖 mod (𝑋 𝑁 + 1)
where:

𝑣 0 = 1 2 • (1 + 𝑋 + • • • + 𝑋 𝑁 -1 ) 𝑣 𝑖 = 𝛼 𝑖,0 + 𝛼 𝑖,𝑁 -1 + (𝛼 𝑖,1 -𝛼 𝑖,0 ) • 𝑋 + • • • + (𝛼 𝑖,𝑁 -1 -𝛼 𝑖,𝑁 -2 ) • 𝑋 𝑁 -1
This factorization makes it possible to compute many LUTs using a unique bootstrapping. Indeed, it is enough to initialize the test polynomial 𝑡𝑒𝑠𝑡𝑣 with the value of 𝑣 0 during bootstrapping. Then, after the BlindRotate operation, one has to multiply the obtained 𝐴𝐶𝐶 by each 𝑣 𝑖 corresponding to the LUT of 𝑓 𝑖 to get 𝐴𝐶𝐶 𝑖 . Figure 2 illustrates the advantage of this method. This optimization reduces the number of bootstrapping required for an operation and, thus, the overall computation time.

A SHORT REMINDER ON THE AES

Advanced Encryption Standard (AES) is the name given to the Rijndael algorithm, the winner of the NIST standardization competition in 2000 [START_REF] Daemen | The Design of Rijndael: AES -The Advanced Encryption Standard (Information Security and Cryptography[END_REF]. It is a symmetric block cipher, defined to work with different key sizes. Several rounds are applied to the original message to obtain an encrypted message. Each round consists of the same operations performed in the same order. We chose to work on the 128-bits AES, which uses 10 rounds. The ciphertext is composed of 16 bytes such as 16 and is encoded in what we call a state matrix in the following way:

𝑐 = 𝑐 0 𝑐 1 • • • 𝑐 15 ∈ (F 2 8 )
𝑐 0 𝑐 4 𝑐 8 𝑐 12 𝑐 1 𝑐 5 𝑐 9 𝑐 13 𝑐 2 𝑐 6 𝑐 10 𝑐 14 𝑐 3 𝑐 7 𝑐 11 𝑐 15
The round operations affect this matrix as follows:

• SubBytes: the SubBytes operation is the only non-linear transformation of the cipher. It is a permutation consisting of an S-box applied to the bytes of the state matrix. As it acts on the individual bytes of the state, it can be parallelized for efficient execution. • AddRoundKey: before the encryption, the secret key is "expanded" into several round keys. The encryption process relies only on these round keys and not the initial secret key. In this transformation, the state is modified by combining it with a round key with the bitwise XOR operation. Of course, to do so, the size of each round key is equal to the size of the ciphertext (in our case, 128 bits encoded in a round key matrix to match the state matrix). That is to say, we have at round 𝑖 ∈ {0, A typical execution of the 128-bits AES begins with the first AddRoundKey directly followed by the first iteration of the rounds. Each round proceeds as follows:

(1) SubBytes (2) ShiftRows (3) MixColumns (4) AddRoundKey Except for the last round, which does not require the MixColumns step.

About the Key Expansion-Key expansion is an operation that may be performed once and for all, from the secret key. Indeed, from the 128-bit key are derived eleven 128-bit round keys, which are used in the AddRoundKey operation. As a result, to evaluate the AES encryption or decryption algorithm, a server only needs to know the round keys. This operation, consisting of XOR and 𝐺𝐹 (256) multiplication, is expensive in the homomorphic domain. It is therefore more efficient for a client to generate its own key, derive the round keys, and then homomorphically encrypt them. Sending these eleven homomorphically encrypted round keys is faster than creating the initial key, encrypting it in homomorphic, and sending it to a server to perform the Key Schedule in the homomorphic domain. For these reasons, we remove the key schedule from the encrypted-domain computations.

AES GOES HOMOMORPHIC

To use the full potential of programmable bootstrapping, we aim to transform the AES algorithm into a succession of LUT evaluations. The permutation given by the Sbox is already in LUT form, so all that remains is to modify the multiplication in 𝐺𝐹 (256) and XOR operations. These are binary operations, which are not immediate to put into LUT form. Indeed, evaluating a LUT is like dereferencing an array from an encrypted index. The LUT evaluation operator, therefore, takes only one input: this encrypted index. If possible, these operations should thus be transformed into unary table indirection operations. Since AES works on a byte-by-byte basis, LUTs must take 8-bits encrypted indexes and return an evaluation of this index on 8 bits. In this section, we explain the various steps involved in transforming the original code into the sequence of 8-bits-to-8-bits LUTs required for efficient homomorphic execution.

Optimizing the Original AES Code

The first step towards the homomorphization of a symmetric cryptosystem is to look at the original code for small changes that would ease the transition. AES makes no exception. However, we quickly observed that the proposed C++ implementation (given by the creators of the AES in [START_REF] Daemen | The Design of Rijndael: AES -The Advanced Encryption Standard (Information Security and Cryptography[END_REF]) is only "pseudo" 8-bits, as a carry that requires an additional ninth bit is necessary for specific bytes operations.

Indeed, this occurs during the 𝐺𝐹 (256) multiplication of two integers (see Fig. 3). To compute this operation efficiently, the authors use a generator 𝛼 = 𝑋 + 1 of the message space 𝐺𝐹 (256) F2 [𝑋 ] /⟨𝑋 8 + 𝑋 4 + 𝑋 3 + 𝑋 + 1⟩. From this, they construct two tables of 256 elements. The first one is Logtable, defined such as Logtable[𝛼 𝑖 ] = 𝑖. The second is Alogtable where Alogtable[𝑖] = 𝛼 𝑖 . Then, instead of a naive 𝐺𝐹 (256) multiplication, the result is obtained with a simple sum of logs and two table indirections. Consequently and as seen in Figure 3, when the mul(word8 𝑎, word8 𝑏) function is called, two 8-bits integers representing the logs are added together. And this sum may exceed 255, which is why the authors apply a %255 to the result. We must change the whole function structure to avoid this overhead and have an actual 8-bits implementation. But it also implies some more profound changes in the entire code.

We deal with this issue by observing that there are only a few calls to the mul function throughout the algorithm. During these calls, the parameter 𝑎 only takes six different values. Indeed, a careful reading of the code shows that 𝑎 ∈ {2, 3, 9, 𝑏, 𝑑, 𝑒}. This means that we can separate the mul(𝑎, 𝑏) function into six mul_a(𝑏) functions, with 𝑎 ∈ {2, 3, 9, 𝑏, 𝑑, 𝑒} by removing the variable 𝑎 of the parameters and just encoding the correct values in the corresponding new functions. Then, to work around the problematic addition of two 8-bits integers, we compute and hardcode six tables Tlog_a such that Tlog_a[𝑏] = (Logtable[a]+Logtable[b])%255. This pre-calculation of the tables enables us to eliminate the problem of the "pseudo" 8-bits implementation. At the same time, this allows us to avoid an explicit modulo operation, which is difficult to efficiently perform in the homomorphic domain. To optimize a step further and reduce the number of tables indirections per encrypted index, we finally consider the six new T_a tables such that T_a[𝑏] =Alogtable[Tlog_a [b]] and use them directly into our implementation. We now have six functions of the following form:

word8 mul_a(word8 b) { if (b) return T_a[b];
else return 0; } It is a simple function, but regarding a homomorphic evaluation, we want our function to be as light as possible regarding the number of operations performed. This is why, to avoid an explicit evaluation of the if condition on 𝑏 (which should be performed via a conditional assignment as FHE disallows branching for obvious reasons), we modified the T_a tables so that they consider the case where 𝑏 = 0. That is to say, for every 𝑎 ∈ {2, 3, 9, 𝑏, 𝑑, 𝑒}, T_a[0] = 0. So our final implementations of the 𝐺𝐹 (256) multiplications are straightforward and only require one indirection, as seen in Figure 4. Therefore, this work on the multiplication functions allows us to operate with only one indirection instead of three indirections, an addition, a modulo operation, and an if condition in the original code version. These changes allow us to convert the initial binary 𝐺𝐹 (256) multiplication operator into several unary ones, which is ideal for a LUT transformation.

When looking for other possible optimizations, we realized we achieved an optimal or nearly optimal form regarding LUT factorization and, thus, LUT-based homomorphic evaluation (at least when starting from the standard AES implementation). Indeed, the structure of the MixColumns instructions prevent further factorization. For example, if one considers the following instructions (in

MixColumn): b[0][j] = mul_2(a[0][j]) ^mul_3(a[1][j]) ^a[2][j] ^a[3][j]; b[1][j] = mul_2(a[1][j]) ^mul_3(a[2][j]) ^a[3][j] ^a[0][j]; b[2][j] = mul_2(a[2][j]) ^mul_3(a[3][j]) ^a[0][j] ^a[1][j]; b[3][j] = mul_2(a[3][j]) ^mul_3(a[0][j]) ^a[1][j] ^a[2][j];
where a denotes the state matrix, then, although the mul_2 and mul_3 LUTs could directly embed the Sbox LUT (assuming the indices are made consistent with the effect of ShiftRows) the two last terms cannot (since they do not require any 𝐺𝐹 (256) multiplication, hence there is no LUT to factor the Sbox LUT with). For instance, a[0][j] is used via mul_2 in b[0][j] calculations and then as is for the others

(b[1][j], b[2][j], b[3][j]
). Finally, as an additional optimization, we have merged the AddRoundKey function with the SubBytes one, as the first always precedes the second one, for increased parallelism efficiency.

And Finally Going Homomorphic for Real

To make the most of TFHE programmable bootstrapping, we implement all operations (XOR, Sbox, 𝐺𝐹 (256) multiplication) via LUT evaluations. As we now have a ciphertext-by-cleartext 𝐺𝐹 (256) multiplication function that is performed only by means of a unique indirection, this causes no issue. And by definition, the Sbox is a LUT. So, the actual work here is to transform the XOR into a LUT. The XOR is a binary operation in which, in our case, inputs are bytes. That means we need a 256 × 256 table to encode our XOR. As [START_REF] Daemen | The Design of Rijndael: AES -The Advanced Encryption Standard (Information Security and Cryptography[END_REF]. The sum of Logtable[a] + Logtable[b] may exceed 255. Therefore, this operation requires more than 8 bits in practice, which is unsuitable for a homomorphic evaluation with an 8-bits processor.

it is a binary operation, we must either use the tree-based method (recalled in Sect. 2.5) or the chaining method, both introduced in [START_REF] Guimarães | Revisiting the functional bootstrap in TFHE[END_REF], which allows bootstrapping with several entries. Table 1 compares the parameters needed depending on the method. In this table, 𝐵 𝑔 and 𝑙 denote the basis and levels associated with the gadget decomposition, 𝐵 𝐾𝑆 and 𝑡 denote the decomposition basis and the precision of the decomposition of the KeySwitch. Finally, 𝑞 denotes the size of the used plaintext size (meaning that the torus is discretized on 𝑞 values), and 𝜖 is the error probability of one MVB evaluation or one evaluation of the chaining method. We also give the noises associated with the TRLWE and TLWE ciphertexts. Note that these parameters are not universal; other parameters with another error probability could also work. Given a message space of size 𝐵, the chaining method requires using a plaintext space of size 2𝐵 2 (instead of only 2 × 𝐵 for the tree-based method). As such, the size of the parameters dramatically increases as the basis 𝐵 increases. This growth of parameters jeopardizes the other speed improvements that could come with the chaining method. Relying on this comparison, we chose to work with the tree-based method, which is thus more efficient in our case.

On top of this, the MVB optimization (Sect. 2.6) allows us to evaluate a LUT with two input digits and one output digit in only two bootstrappings. Still, as the parameters needed for calculation in basis 256 are significantly larger (we have to work on the discretized torus on 512 values, using only the positive half), it will not be efficient. Indeed, we need to use the cyclotomic polynomial 𝑋 32768 + 1, which results in a very slow bootstrapping and a non-implementable small variance (still for 135 bits of security according to the lattice-estimator). Recall that the parameter sets given for the chaining method are for illustration purposes as we do not use that method in the paper and just aim at showing that the tree-based method is more practical for basis 16. This makes it interesting to break down the messages into smaller bases and to use the tree-based method for every LUT evaluation. For example, in order to perform an indirection into a 256 8-bits entries LUT, when working in basis 16, we have to split the original table into two 256 4-bits entries. For each of these two tables, we apply a two levels tree-based method where, as illustrated by Figure 5 with the LUT of the identity function, the first level extracts one encrypted 4-bits value from each block of 16 values (assuming, as we do, that the first digit contains the least significant 4 bits) in the original table. This leads to a 16 (encrypted) 4-bits entries table, and the second level extracts the encrypted 4-bits value from the latter.

More generally, the number 𝑁 𝑏𝑜𝑜𝑡 of bootstrapping needed with the tree-based method increases with the number of input digits 𝑑 and the number of output digits 𝑑 ′ as well as the basis 𝐵 chosen for the decomposition. In fact, for a complete evaluation of a LUT via the classic tree-based method, we have 𝑁 𝑏𝑜𝑜𝑡 = 𝑑 ′ × 𝑑 -1 𝑖=0 𝐵 𝑖 and with the MVB, we have

𝑁 𝐵 ′ 𝑏𝑜𝑜𝑡 = 𝑑 ′ × (1 + 𝑑 -2 𝑖=0 𝐵 𝑖 ).
Table 2 provides the number of bootstrapping per operation (256 × 8bits LUT dereferencing and 8-bits-by-8-bits XOR) depending on the decomposition basis for the ciphertexts (using the MVB-optimized tree-based method).

Although the tree-based method requires more operations, the smaller the decomposition basis, the smaller the parameters to be used. Small ciphertext sizes allow small parameters and, therefore, faster operations. A tradeoff must therefore be achieved between the number of bootstrapping operations performed and the parameters' size. For this reason, we produce a benchmark of the execution time of LUT evaluation depending on the ciphertext decomposition. We consider the following different decompositions:

• basis 256: it is not a decomposition per se, but we have to be sure that this basis is not the most advantageous one • basis 16: the message is decomposed into two digits in basis 16 (4 bits per digit) • basis 8: the message is decomposed into three digits in basis 8 (3 bits per digit, but only 2 bits for the most significant digit) • basis 4: the message is decomposed into four digits in basis 4 (2 bits per digit)

For this, we first implemented an efficient homomorphic operator to evaluate any LUT using TFHE programmable bootstrapping and MVB. We then used it with our different decompositions and their associated parameter sets. The results of this experiment are given in Table 3. It is thus clear that basis 16 is the optimal choice regarding the execution time per LUT evaluation as well as for a full AES execution (as the overall number of LUT evaluations is independent of the basis choice). Note that there is no linearity between the timings due to the number of multiplications involved by the MVB during the first bootstrapping. Additionally, Table 4 provides the number of bootstrapping needed with the MVB optimization for a full AES evaluation, depending on the decomposition basis. We now have everything needed to obtain homomorphic versions of the mul_a functions. To do so, as already emphasized, we decompose every table T_a into two tables of 256 basis 16 digits, i.e., one per digit of the decomposition. We do the same with the Sbox table. Finally, we use a similar approach to create the 4-bits-by-4-bits XOR operator. We do so by means of a 16 × 16 table with 4-bit entries, which is much easier to handle than a 256 × 256 (with 8-bits

Details about the Homomorphization of the XOR Operator

Following the previous section, we give details about our homomorphic XOR operator, which is a relevant concrete example of our tree-based method approach. On the one hand, we have to transform the classic 8-bits-to-8-bits XOR operator to fit the chosen decomposition basis, and on the other hand, we have to transform it so we can evaluate it via a LUT. As already emphasized (Table 3), basis 16 is the best one for our AES homomorphic execution attempt. This means that we decompose our messages on the following form 𝑚 = 𝑚 0 +𝑚 1 • 16 and that the corresponding ciphertext is a vector of the form 𝑐 = (𝑐 0 , 𝑐 1 ) = ([𝑚 0 ] , [𝑚 1 ]). So, to compute the XOR of two ciphertexts 𝑐 and 𝑐 ′ in basis 16, we have to compute:

𝑐 ⊕ 𝑐 ′ = (𝑐 0 , 𝑐 1 ) ⊕ (𝑐 ′ 0 , 𝑐 ′ 1 ) = (𝑐 0 ⊕ 𝑐 ′ 0 , 𝑐 1 ⊕ 𝑐 ′ 1 )
We thus have to evaluate a 4-bits-by-4-bits XOR on the ciphertexts. For this, we dereference the double input table of 4-bits XOR, which has size 16 × 16. The tricky part is to choose the correct way to construct the test polynomials from this table for the tree-based method. Indeed, on the first step of the tree-based method, we have 16 test polynomials. Each one must encode the coefficient of a unary XOR. That is to say, the first polynomial 𝑃 0 encodes the values of the unary operation xor_by_0(𝑖) such that for 𝑖 ∈ {0, 1, • • • , 15} xor_by_0(𝑖) = 𝑖 ⊕ 0, the second polynomial 𝑃 1 encodes xor_by_1, etc. So when applying the bootstrapping on these polynomials with selector 𝑐 0 = [𝑚 0 ], we obtain 16 new ciphertexts encoding xor_by_m0 that we put together on a polynomial 𝑃 𝑓 𝑖𝑛𝑎𝑙 . We then apply the bootstrapping on the new polynomial 𝑃 𝑓 𝑖𝑛𝑎𝑙 with selector 𝑐 ′ 0 = 𝑚 ′ 0 , and we obtain the final results, which is an encryption of 𝑐 0 ⊕ 𝑐 ′ 0 . The method is illustrated in Fig. 6. By using the same method, we can easily compute 𝑐 1 ⊕ 𝑐 ′ 1 , and thus obtain the vector 𝑐 ⊕ 𝑐 ′ = (𝑐 0 ⊕ 𝑐 ′ 0 , 𝑐 1 ⊕ 𝑐 ′ 1 ) encoding 𝑚 ⊕ 𝑚 ′ .

Remarks and Perspectives on LUT Evaluated Operations

In this section, we briefly present our various operations in a unified way and illustrate that it allows to express more general operations than needed for our AES execution attempt. A distinction must be made between ciphertext-ciphertext binary operators and cleartextciphertext binary operators. The latter is the case of the 𝐺𝐹 (256) multiplication operator in the AES algorithm, which we turned into several unary operators (the mul_a functions, which could easily have been as many as 256 if it had been required). Regarding the case of ciphertext-ciphertext operators, we can divide them into two types: the bitwise ones and the others.

Bitwise ciphertext-ciphertext operator -The AES implementation requires an XOR operation, which can finally be implemented using the same tool as for the unary functions (Sbox, mul_a). Slightly loosely speaking, let us say that our 8-bits-to-4-bits table dereferencing tool takes the following form LUT(𝑐 0 , 𝑐 1 , 𝑡𝑎𝑏𝑙𝑒) and in basis 16 LUT(𝑐 0 , 𝑐 1 , 𝑡𝑎𝑏𝑙𝑒) = 𝑡𝑎𝑏𝑙𝑒 [𝑐 0 + 𝑐 1 • 16]. To evaluate the unary operator of the Sbox (or mul_a), one must call LUT(𝑐 0 , 𝑐 1 , Sbox_0) and then LUT(𝑐 0 , 𝑐 1 , Sbox_1), where Sbox_0 and Sbox_1 respectively denote the LSB and MSB of Sbox, to obtain the two result digits and reform Sbox[𝑐 0 + 𝑐 1 • 16] (respectively LUT(𝑐 0 , 𝑐 1 , mul_a_0) and LUT(𝑐 0 , 𝑐 1 , mul_a_1)). But for XOR, we use LUT(𝑐 0 , 𝑐 ′ 0 , XOR) and then LUT(𝑐 1 , 𝑐 ′ 1 , XOR) as explained above. In the first case, we use the same ciphertexts on different tables; in the second case, we use different ciphertexts on the same table. This means a binary bitwise homomorphic operator is as easy (and costly) to implement as a unary non-bitwise operator. Furthermore, such a binary operator can easily be adapted to any bitwise operation (AND, NOR, etc.) or unary operation by modifying the initial polynomials accordingly.

Non-bitwise ciphertext-ciphertext operator -Although not needed for AES, the case of non-bitwise ciphertext-ciphertext operators is a bit more involved. Let us take the example of the addition of two encrypted basis 256 messages 𝑐 and 𝑐 ′ decomposed in basis 16. Two different tables with 256 4-bits entries then come into play: ADD, which encodes the addition of two digits, and CAR, which handles the carry. So, using the same dereferencing tool as above, a simple addition over Z 256 can be broken down into 4 invocations of the operator:

LUT (𝑐 0 , 𝑐 ′ 0 , ADD ) = 𝑟𝑒𝑠 0 LUT (𝑐 0 , 𝑐 ′ 0 , CAR ) = 𝑟 LUT (𝑐 1 , 𝑐 ′ 1 , ADD ) = 𝑣 LUT (𝑣, 𝑟, ADD ) = 𝑟𝑒𝑠 1
where 𝑟𝑒𝑠 0 + 𝑟𝑒𝑠 1 • 16 = 𝑐 + 𝑐 ′ . The basic dereferencing tool is the same for all of these cases and is, in essence, universal, although on a case-by-case basis, there might be more optimized methods for certain operations (e.g., by combining the tree-based and chainbased methods, the latter being, for example, more efficient for carry computations). This gives interesting perspectives for building a more general optimized "instruction set".

EXPERIMENTAL RESULTS

TFHE Library Implementation

We worked with TFHElib3 using the parameters given in Table 1 for the tree-based method with basis 16. As a reminder, we take 𝐵 𝑔 = 256 and 𝑙 = 3 as the basis and levels associated with the gadget decomposition. For the KeySwitch, we take 𝐵 𝐾𝑆 = 1024 and 𝑡 = 2 as the decomposition basis and the precision of the decomposition. Finally, as we work with a plaintext space of size 16, we work on the discretized torus on 𝑞 = 32 values. We chose these parameters specifically to enable MVB evaluation with error probability 𝜖 = 2 -23 . According to the lattice estimator 4 , these give us at least 128 bits of security. The first step was to implement the tree-based method and its MVB optimization. Indeed, none of these methods are part of the TFHE library. Furthermore, we had to adapt the code to be able to use vectors of TLWE and TRLWE samples (as our ciphertexts are decomposed into 2 digits in basis 16, we put them in vectors of 2 elements). Then, we modified the MVB operator to turn it into a generic LUT evaluation tool. This new operator takes as input a vector of ciphertexts, the number 𝑑 of input digits, the number 𝑑 ′ of output digits, the decomposition basis 𝐵 and the LUT to be evaluated. Depending on how it is called, this operator can perform our XOR as well as our 𝐺𝐹 (256) multiplications and Sbox.

Parallelization

The purpose of transciphering is to avoid transferring large ciphertexts. Since a server has more computing power than a client, it can efficiently exploit several cores to parallelize computations and optimize execution times. With that respect, AES is naturally parallelizable up to a certain degree: within each round step, operations can be performed simultaneously on each byte of the state matrix (except for the ShiftRows step). This is why, in continuing our work, we use the OpenMP library to parallelize our code and optimize the execution times of our homomorphic AES. We run our tests on two machines. The first one is a 12th Gen Intel(R) Core(TM) i7-12700H CPU (using six cores) laptop with 64 Gib total system memory with an Ubuntu 22.04.2 LTS server. The second is an AMD EPYC 7702P 64-cores Processor server with an Ubuntu 20.04.6 LTS server. In the following sections, we will refer to them respectively as i7-laptop and AMD-server.

Using the OpenMP library, we first parallelize onto the six cores on the i7-laptop. As it is a standard laptop, it is interesting to see how a partially parallelized homomorphic AES can work. Indeed, it gives a reasonable hint of the potential of a larger scale parallelization. That is why we parallelize every round function, except for the ShiftRows one, as it only involves ciphertexts reorganization within the state matrix. Still, thanks to OpenMP, we could use 16 cores for a larger-scale parallelization on the AMD-server (which hosts 64 cores in total). As the server has less powerful cores, it is slower than the i7-laptop in sequential time. But, by running a sequential execution and a full parallelized one, we can measure the speedup factor induced by parallelization. For instance, the execution parallelized on 16 cores is 9.5 times faster than the sequential one on the AMD-server. Results can be found in Table 5. The last column of this table indicates the ratio between the (fastest) "i7-server" time and the other (slower) ones.

A perspective is to go further in the parallelization of the AES evaluation. Indeed, we can also parallelize the tree-based method and execute the computation of each decomposition digit at the same time. This method will need 32 cores, but it can divide the execution time almost by 2.

Computation Times

This section summarizes the results of our implementations in Table 5. Even if our measured speedups are not linear in the number of cores, we improve the state of the art as the so far best-known implementation runs in 4.2 minutes with 16-thread parallelization [START_REF] Stracovsky | Faster evaluation of AES using TFHE[END_REF] using the standard programmable bootstrapping in base 16 as implemented in Zama's Concrete5 . However, note that no implementation details were provided (to the best of our knowledge) beyond a poster presented by Stracovsky et al. at FHE.org 2022 [START_REF] Stracovsky | Faster evaluation of AES using TFHE[END_REF]. In [START_REF] Gentry | Homomorphic Evaluation of the AES Circuit[END_REF], Gentry et al. use a packed implementation with BGV so that each slot can hold a state byte. The round operations are computed as single linear transformations over 𝐺𝐹 (2 8 ) 16 . These linear transformations combine several permutations that are computed via automorphism evaluation and additions. In [START_REF] Mella | On the Homomorphic Computation of Symmetric Cryptographic Primitives[END_REF], Mella and Susella improve these latter results by using an optimized data representation to further factor some of these operations in ciphertext slots. As they relied heavily on batching, these works focused primarily on optimizing the amortized time metric which is relevant only when one has to perform many independent AES executions in parallel. As Table 5 shows, we obtain a sequential time comparable to [START_REF] Stracovsky | Faster evaluation of AES using TFHE[END_REF] on i7-laptop and a very improved one with parallel execution. Indeed, our parallelized version on i7-laptop is 4.6 times more efficient than Stracovsky et al. 's [START_REF] Stracovsky | Faster evaluation of AES using TFHE[END_REF] and almost 20 times more efficient than Gentry et al. 's. [START_REF] Gentry | Homomorphic Evaluation of the AES Circuit[END_REF] in terms of latency. Furthermore, our parallelized version on AMD-server is almost 7 times faster than [START_REF] Stracovsky | Faster evaluation of AES using TFHE[END_REF] and 30 times faster than [START_REF] Gentry | Homomorphic Evaluation of the AES Circuit[END_REF] as it only takes 3% of its execution time.

CONCLUSION

In this paper, we have proposed a homomorphic AES implementation relying on two encrypted-domain "instructions", a unary 256-by-8 bits table indirection, and a binary 8-bit XOR instruction, both running over encryptions of nibbles (the cute name for hex digits) and relying on functional bootstrapping for efficiency. Our work illustrates that, even when starting from the standard AES implementation, this approach significantly improves the state of Table 5: Our execution times compared to the state of the art (the "i7-server" time corresponds to an extrapolation of the speedup observed with 16 threads on the AMD-server from the i7-laptop sequential time as we did not have an i7-based server to perform experiments). execution time time ratio i7-laptop (1 thread) 4.5 mins = 270 secs 9. 4 i7-laptop (6 threads) 54.31 secs 1.9 AMD-server (1 thread) 5.7 mins = 342 secs 11.9 AMD-server (16 threads) 36.39 secs 1.3 "i7-server" (16 threads) 28.73 secs 1 Gentry et al. [START_REF] Gentry | Homomorphic Evaluation of the AES Circuit[END_REF] (1 thread) 18 mins 37.6 Mella and Susella [START_REF] Mella | On the Homomorphic Computation of Symmetric Cryptographic Primitives[END_REF] (1 thread) 22 mins 45.9 Stracovsky et al. [START_REF] Stracovsky | Faster evaluation of AES using TFHE[END_REF] (16 threads) 4.2 mins = 252 secs 8.8 the art of homomorphic AES execution timings. In terms of perspectives, beyond improved parallelism, it would be interesting to consider other non-standard forms for the AES as a starting point in a search for ones that may lead to smaller numbers of homomorphic operations (i.e., in fine, fewer bootstrapping). For example, a non-public implementation optimized for constrained embedded systems by inlining, loop unrolling, and careful instructions reorganization was brought to our attention. In our terms, such an implementation would lead to 68 XOR and 32 LUT per AES round (vs. 64 XOR and 48 LUT when starting from the standard implementation as we did). Although we do not expect to gain one or more orders of magnitude, this could lead to an additional 10% performance improvement, at least in the sequential execution times we report in this paper (as this version appears more difficult to parallelize at first glance). We could also try to exploit the emerging full-Torus functional bootstrapping methods, such as ComBo [START_REF] Clet | ComBo: A Novel Functional Bootstrapping Method for Efficient Evaluation of Nonlinear Functions in the Encrypted Domain[END_REF], which may allow smaller parameters and, thus, faster bootstrapping evaluations. Another approach that could potentially lead to further timing improvements would be to investigate how the techniques introduced in [5] may help to either or both increase the basis size or factor larger portions of the AES algorithm in single WoP-PBS operator evaluations. Yet another interesting perspective would be, as hinted in Sect. 4.4, to extend the embryonic "instruction set" defined in this paper in order to apply this approach to other algorithms more easily.

Figure 1 :

 1 Figure 1: TFHE Bootstrapping examples: the outer circles describe the inputs to the bootstrapping (i.e., ciphertexts over T).Meanwhile, the inner circles represent the coefficients of the test polynomial 𝑡𝑒𝑠𝑡𝑣. One of these coefficients is extracted as the output of the bootstrapping after the BlindRotate.

Figure 2 :

 2 Figure 2: Illustration of the MVB optimization. (a) represents the classic method to process several bootstrapping, while (b) represents the MVB optimization. As seen here, it reduces the number of BlindRotate operations, which is the most expansive one of the bootstrapping.

Figure 4 :

 4 Figure 4: Our final 𝐺𝐹 (256) multiplication operator. The new multiplication functions are very simple and easy to execute in the homomorphic domain, thanks to pre-calculation and table hardcoding.

Figure 3 :

 3 Figure3: The original 𝐺𝐹 (256) multiplication implementation given in[START_REF] Daemen | The Design of Rijndael: AES -The Advanced Encryption Standard (Information Security and Cryptography[END_REF]. The sum of Logtable[a] + Logtable[b] may exceed 255. Therefore, this operation requires more than 8 bits in practice, which is unsuitable for a homomorphic evaluation with an 8-bits processor.

Figure 5 :

 5 Figure 5: Illustration of the tree_based method on the identity function. The message is 𝑚 = 9 = 1 • 4 0 + 2 • 4 1 and its encryption is 𝑐 = (𝑐 0 , 𝑐 1 ) = ([1] , [2]). Red arrows indicate bootstrappings.

Figure 6 :

 6 Figure 6: The principle of the tree-based method applied to the XOR with basis 16. The green color indicates that the content of the box is encrypted.

  • TRLWE sample: A pair (𝑎, 𝑏) ∈ T 𝑁 [𝑋 ] 𝑘+1 , where 𝑎 is uniformly sampled from T 𝑁 [𝑋 ] 𝑘 and 𝑏 = ⟨𝑎, 𝑠⟩ + 𝑒. The secret key 𝑠 is uniformly sampled from B 𝑁 [𝑋 ] 𝑘 , the error 𝑒 ∈ T is a polynomial with random coefficients sampled from a Gaussian distribution with mean 0 and standard deviation 𝜎. One usually chooses 𝑘 = 1; therefore, 𝑎 and 𝑏 are vectors of size 1 whose coefficient is a polynomial.• TRGSW sample: a vector of (𝑘 + 1)𝑙 TRLWE fresh samples.

  TLWE 𝑠 (𝑚 𝑖 ) where 𝑚 𝑖 is the 𝑖 𝑡ℎ coefficient of the polynomial 𝑚. This paper will refer to this algorithm as SampleExtract.• Public Functional Keyswitching: allows the switching of keys and parameters from 𝑝 ciphertexts𝑐 𝑖 ∈ TLWE 𝑘 (𝑚 𝑖 ) to one 𝑐 ′ ∈ T(R)LWE 𝑠 (𝑓 (𝑚 1 , • • • , 𝑚 𝑝 ))where 𝑓 is a public linear morphism between T 𝑝 and T 𝑁[𝑋 ]. That is to say, this operation not only allows the packing of TLWE ciphertexts in a TRLWE ciphertext, but it can also evaluate a linear function 𝑓 over the input TLWEs. This paper will refer to this algorithm as KeySwitch.

	It returns a ciphertext 𝑐 ′ ∈ TRLWE 𝑘 (𝑚 • 𝑋 ⟨𝑎,𝑠 ⟩-𝑏 ). This
	paper will refer to this algorithm as BlindRotate.
	• TLWE Sample Extract: extracts a coefficient of a TRLWE
	ciphertext and converts it into a TLWE ciphertext. It takes
	as inputs both a ciphertext 𝑐 ∈ TRLWE 𝑠 (𝑚) and an index
	𝑝 ∈ {0, • • • , 𝑁 -1}. The result is a TLWE ciphertext 𝑐 ′ ∈

, • • • , 𝑎 𝑝 , 𝑏) where ∀𝑖, 𝑎 𝑖 ∈ Z 2𝑁 , and a TRGSW ciphertext encrypting the secret key 𝑠 = (𝑠 1 , • • • , 𝑠 𝑝 ).

  𝑛⟧ where 𝑆 ′ is the TRLWE interpretation of a secret key 𝑠 ′ .

	Ensure: a TLWE sample 𝑐 ∈ TLWE 𝑠 (𝑥 .𝑚)
	1: Let m = 1 2 𝑚 ∈ T (pick one of the two possible values)
	2: Let 𝑏 = ⌊2𝑁𝑏⌉ and 𝑎 𝑖

  𝑐 0 𝑐 4 𝑐 8 𝑐 12 𝑐 1 𝑐 5 𝑐 9 𝑐 13 𝑐 2 𝑐 6 𝑐 10 𝑐 14 𝑐 3 𝑐 7 𝑐 11 𝑐 15 ⊕ 𝑘 𝑖,0 𝑘 𝑖,4 𝑘 𝑖,8 𝑘 𝑖,12 𝑘 𝑖,1 𝑘 𝑖,5 𝑘 𝑖,9 𝑘 𝑖,13 𝑘 𝑖,2 𝑘 𝑖,6 𝑘 𝑖,10 𝑘 𝑖,14 𝑘 𝑖,3 𝑘 𝑖,7 𝑘 𝑖,11 𝑘 𝑖,15 = 𝑐 0 ⊕ 𝑘 𝑖,0 𝑐 4 ⊕ 𝑘 𝑖,4 𝑐 8 ⊕ 𝑘 𝑖,8 𝑐 12 ⊕ 𝑘 𝑖,12 𝑐 1 ⊕ 𝑘 𝑖,1 𝑐 5 ⊕ 𝑘 𝑖,5 𝑐 9 ⊕ 𝑘 𝑖,9 𝑐 13 ⊕ 𝑘 𝑖,13 𝑐 2 ⊕ 𝑘 𝑖,2 𝑐 6 ⊕ 𝑘 𝑖,6 𝑐 10 ⊕ 𝑘 𝑖,10 𝑐 14 ⊕ 𝑘 𝑖,14 𝑐 3 ⊕ 𝑘 𝑖,3 𝑐 7 ⊕ 𝑘 𝑖,7 𝑐 11 ⊕ 𝑘 𝑖,11 𝑐 15 ⊕ 𝑘 𝑖,15 • ShiftRows: the ShiftRows step is a byte transposition that cyclically shifts the rows of the state over different offsets. For AES-128, row 0 is shifted over 0 bytes, row 1 over 1 byte, row 2 over 2 bytes and row 3 over 3 bytes. As this operation only alters the position of the bytes in the state matrix, it does not require a homomorphic equivalent (instead of shifting regular bytes, we shift homomorphic bytes). • MixColumns: the MixColumns step is operating on the state column by column via matrix multiplication. But in practice, the authors do not implement the naive matrix product but work with each byte of the state matrix individually. They use a mix of scalar 𝐺𝐹 (256) multiplication and XOR, as visible in the original code. Therefore, this operation can be parallelized.

• • • , 9}:

Table 1 :

 1 Parameter sets depending on the LUT evaluation method and the chosen decomposition basis (𝜆 ≈ 128).

		basis 𝑛	𝑁	𝑙 𝐵 𝑔	𝐵 𝐾𝑆 𝑡 𝑞	𝜖	TRLWE noise TLWE noise
		4	850	2048	2 2048	1024 2 32	2 -32 9.6 × 10 -11	1.27 × 10 -6
	chaining method	8	1024 8192	1 268435456	1024 2 128 2 -27 10 -45	5.6 × 10 -8
		16	1100 32768 1 4294967296 8192 2 512 2 -30 1.4 × 10 -8	10 -248
		4	700	1024	5 16	1024 2 8	2 -30 5.6 × 10 -8	1.9 × 10 -5
	tree-based method 8	700	2048	2 2048	1024 2 16	2 -23 9.6 × 10 -11	1.9 × 10 -5
		16	1024 2048	3 256	1024 2 32	2 -23 9.6 × 10 -11	6.5 × 10 -8

Table 2 :

 2 Number of bootstrapping for one XOR or one LUT evaluation (using MVB) depending on the decomposition basis.

	basis # XOR # LUT
	256	2	1
	16	4	4
	8	6	30
	4	8	88
	entries) needed without the decomposition (yet another pro for
	the basis 16 decomposition). This being done, we now have duly
	"translated" our optimized-for-FHE version of the AES code into a
	real homomorphic one.

Table 3 :

 3 Unitary timings for bootstrapping and full LUT evaluation depending on basis and parameter choices (𝜆 ≈ 128).

	basis single boot. complete	𝑛	𝑁
			LUT eval.
	256	1.5s	1.5 s	1024 32768
	16	0.029s	0.3 s	1024 2048
	8	0.015s	1.4 s	700	2048
	4	0.007s	2.0 s	700	1024

Table 4 :

 4

	basis function	# XOR # LUT # Boot
		AddRoundKey 16	0	32
		SubBytes	0	16	16
	256	MixColumns	48	32	128
		Round	64	48	176
		Full AES	608	448	1664
		AddRoundKey 16	0	64
		SubBytes	0	16	64
	16	MixColumns	48	32	320
		Round	64	48	448
		Full AES	608	448	4224
		AddRoundKey 16	0	96
		SubBytes	0	16	480
	8	MixColumns	48	32	1248
		Round	64	48	1824
		Full AES	608	448	17088
		AddRoundKey 16	0	128
		SubBytes	0	16	1408
	4	MixColumns	48	32	3200
		Round	64	48	4736
		Full AES	608	448	44288

Number of operations and the corresponding number of bootstrapping depending on decomposition basis. Readers are reminded that AES-128 consists of a RoundKey, followed by 9 classic rounds and a final round that does not include the MixColumns operation.

https://tfhe.github.io/tfhe/

https://tfhe.github.io/tfhe/

https://github.com/malb/lattice-estimator

https://github.com/zama-ai/concrete

ACKNOWLEDGEMENTS

This work was supported by the France 2030 ANR Project ANR-22-PECY-003 SecureCompute. The authors also would like to thank Nicolas Belleville for sharing with them the non-standard AES implementation referred to in the conclusion.