Léleng Kebalo
email: kebalo.leleng@gmail.com.

Effects of oil price shocks on economic sectors of net oil-importing countries: case of Togo

, we find that the impact of world oil price on economic sectors is nonlinear. Thus, positive changes in world oil price do not affect the value-added of economic sectors considered while the negative changes in oil price contribute to improve significantly the value-added of primary and secondary sectors, but not the tertiary sector. Finally, our analysis shows that the value-added of primary and secondary sectors affect respectively the value-added of the tertiary sector. The inverse is not true. This paper recommends that Togo must seek to take benefit from all negative changes in world oil price for boosting the value-added of their economic sectors.

Special thanks to John P. Conley and the two anonymous referees for the comments and suggestions that

Introduction

The relationship between energy prices, in particular, oil prices and macroeconomic performance, has drawn attention in many recent studies. The literature boom on the influence of oil prices on an economy was initiated by [START_REF] Hamilton | Oil and the macroeconomy since World War II[END_REF], who intended to measure the impact of oil prices on US macro-economic aggregates. Treating the oil price as an exogenous variable, [START_REF] Hamilton | Oil and the macroeconomy since World War II[END_REF] found that oil prices had a significant impact on the US economy. [START_REF] Mork | Oil and the macroeconomy when prices go up and down: an extension of Hamilton's results[END_REF] 1 , [START_REF] Lardic | The impact of oil prices on GDP in European countries: an empirical investigation based on asymmetric cointegration[END_REF] and Cunado and de [START_REF] Cunado | Oil prices, economic activity and inflation: evidence for some Asian countries[END_REF] found consistent negative impacts of oil prices on GDP for industrialized, industrializing, oil-importing and oil-exporting economies. Like [START_REF] Hamilton | Oil and the macroeconomy since World War II[END_REF], the negative relationship between oil shocks and growth is supported by a series of empirical works, such as [START_REF] Gisser | Crude oil and the macroeconomy: tests of some popular notions[END_REF], [START_REF] Rotemberg | Imperfect competition and the effects of energy price increases on economic activity[END_REF], [START_REF] Brown | Energy prices and aggregate economic activity: an interpretative survey[END_REF], and [START_REF] Zhang | Oil shock and economic growth in Japan: a nonlinear approach[END_REF]. However, it was shown by [START_REF] Hooker | What Happened to the Oil Price-Macroeconomy Relationship?[END_REF] using data after 1985 to be less convincing. But a consensus was found by [START_REF] Blanchard | The macroeconomic effects of oil price shocks: why are the 2000s so different from the 1970s?[END_REF] and Baumeister andPeersman (2013a, 2013b) that the impact of oil price on economies declined strongly over time.

Oil price shocks affect economies differently depending on whether they are net oil-importers or oil-exporters (Berunment et al., 2010). Oil price increases might be considered bad for net oilimporting countries [START_REF] Yanagisawa | Impact of Rising oil prices on the Macro Economy[END_REF], Berunment et al., 2010) but good news for oil-exporting countries. The reverse might be expected for oil price decreases. High oil prices for net oilimporting countries could lead to high import costs with an adverse effect on GDP, exchange rate, inflation and balance of payment. For this category of countries, the immediate effect of positive oil price shocks is to increase the cost of production and therefore affect negatively growth [START_REF] Sachez | Welfare effects of rising oil prices on oil-importing developing countries[END_REF]. Higher oil prices lower disposable income and this decreases consumption. However, high oil prices for oil-exporters improves the general balance of payment due to the increase in oil revenue. Besides, high oil price volatility increases uncertainty regarding cashflows which can be challenging for government in policy decisions. Many studies have been conducted on oil-exporting countries around the world [START_REF] Alimi | The threshold effect of oil-price shocks on economic growth[END_REF][START_REF] Anyanwu | Regional Comparison of Direct Investment to Africa: Empirical Analysis[END_REF][START_REF] Anyanwu | Does Oil Wealth Affect Democracy in Africa?[END_REF][START_REF] Anyanwu | Oil Wealth, Ethno-Religious-Linguistic Fractionalization and Civil Wars in Africa: Cross-Country Evidence[END_REF][START_REF] Ogwumike | Resource-led Development: An Illustrative Example from Nigeria[END_REF][START_REF] Jimenez-Rodriguez | The Impact of Oil Price Shocks: Evidence from the Industries of Six OECD Countries[END_REF]Sanchez, 2005 and[START_REF] Cunado | Oil prices, economic activity and inflation: evidence for some Asian countries[END_REF][START_REF] Gracia | Do Oil Price Shocks Matter? Evidence for some European countries[END_REF].

Oil price shocks affect oil-importing economies differently depending on whether they are developed or developing countries. Several empirical studies on the impact of oil price shocks on economic components conducted for oil-importing developed economies such as United States, European Union, and Japan as in [START_REF] Kurihara | Oil prices and economic growth in developed countries[END_REF] and [START_REF] Akira | Impact of Rising Oil Prices on the Macroeconomy[END_REF], show that oil price increase causes economic growth. However, studies carried out on oil-importing developing countries such as Kenya, Nicaragua, Tanzania and Thailand reveal that increase in oil price generally has a negative impact in GDP causing an increase in unemployment, higher consumer prices and reduced welfare for citizens [START_REF] Sachez | Welfare effects of rising oil prices on oil-importing developing countries[END_REF].

However, recent studies of this relationship on small oil-importing countries claimed that rising oil prices will stimulate economic growth (Gbatu et al., 2017, Liew and[START_REF] Liew | Oil Price Shocks and Sectoral Outputs: Empirical Evidence from Malaysia[END_REF], which is not consistent with other studies that claimed that rising oil prices have an adverse effect on net oil-importers [START_REF] Lemazoshvili | Effect of Oil Price Shocks on Oil Importing Developing Economies: the Case of Georgia and Armenia[END_REF]Shabhaz et al., 2018). Using an unrestricted vector autoregressive (VAR) model, [START_REF] Wesseh | Exchange rate fluctuations, oil price shocks and economic growth in a small net-importing economy[END_REF] evaluate the real impacts of oil price fluctuations on real GDP in Liberia, a small country in West Africa, and find that, contrary to the bulk of the literature, increase in oil price is sometimes good for net-oil importers and this is, according to [START_REF] Wesseh | Causal independence between energy consumption and economic growth in Liberia: Evidence from a non-parametric bootstrapped causality test[END_REF], probably due to the fact that oil price increase in Liberia is linked to an increased substitution for other factors of production particularly labor. And since the service sector constitutes a large portion of real GDP in Liberia (more than 50% of real GDP comes from services), the substitution of energy for labor (due to higher oil prices) leads to an increase in aggregate output.

Concerning oil-importing countries (developing or developed countries), several studies have investigated the relationship between oil price shocks and economic activities since the first oil crisis of the 1970s, but no study has focused on that of Togo, a small net oil-importing country. For improving its economic growth with the purpose to face his economic and social challenges, Togo has decided to increase its economic production. Therefore, the level of energy consumption and more precisely that of oil has increased (Figure 1). Thus, since 1986, Togo has used more energy to boost their economic production and therefore improve the level of its economic growth. But compared to the huge oil-importer such as South Africa2 for example, Togo remains a small net oil-importer. Being a net oil-importer, Togo's economy depends therefore on the world oil prices. The different sectors (primary, secondary and tertiary) of Togolese economy need oil. Then, it would be possible to see these three sectors be affected by world oil price fluctuations. Indeed, according to Berunment et al. (2010), oil prices are not stable, and the fluctuations in oil price have significant effects on the economic activity of a country. [START_REF] Liew | Oil Price Shocks and Sectoral Outputs: Empirical Evidence from Malaysia[END_REF] 3 give real empirical proof that oil price fluctuations can affect the economic sectors of a country.

However, in the context of Togo, it is often argued that, because of the government subventions on oil-derived goods, world oil price shocks/fluctuations do not affect the different sectors of the Togolese economy. Thus, Togolese economy is perceived as resilient to world oil price fluctuations. de jure, this claim seems true when we analyze the economic environment; but is not supported by studies. However, de facto, it is difficult to find legal published documents that show and explain the subvention mechanisms defended and put in place by Togolese government which allow to conclude that, oil price fluctuations do not affect the Togolese economy. Besides, there is no serious study or research paper carried out on the effect on world oil price shocks on Togolese economy that allows to confirm that the Togolese economy is resilient or not to the world oil price fluctuations. This paper fills these gaps in the literature.

The contribution of our paper is two-fold. There is no study in our knowledge that has investigated the linear and nonlinear impact of world oil price shocks on Togo's economy. The second contribution is the analysis of world oil price shocks on the three sectors of an economy; this for knowing the economic sector that is the most affected by the oil price shocks. By performing this kind of analysis, we would be able to conclude precisely on the effect of world oil price on economic sectors of Togo. This allows to propose policies need to make the Togolese economy more resilient to world oil price shocks.

To fulfil the lack of analysis in the context of Togo in the literature, we intend to analyze the effect of world oil price shocks on the three different sectors (primary, secondary and tertiary) of Togolese economy based on an annual time series from 1970-2017. To achieve our purpose, we use like in [START_REF] Du | The relationship between oil price shocks and China's macroeconomy: An empirical analysis[END_REF], an unrestricted vector autoregressive (VAR) model which is one the leading approach employed in the existing literature. According to Schirber (2009), [START_REF] Mork | Oil and the macroeconomy when prices go up and down: an extension of Hamilton's results[END_REF]), Berunment et al. (2010[START_REF] Hamilton | A neoclassical model of unemployment and the business cycle[END_REF], the increase and decrease of oil price can have a different impact on an economic aggregate. Based on their results, besides the analysis with the linear impact model, we use an unrestricted VAR asymmetric impact model proposed by [START_REF] Mork | Oil and the macroeconomy when prices go up and down: an extension of Hamilton's results[END_REF] to analyze the asymmetric effects of oil price shocks on Togo's economic sectors. We carry out Granger causality tests, impulse-response functions and variance decomposition in both linear and nonlinear models (asymmetric world oil price shock).

By using a linear impact model, the results show on one hand that, the world oil price does not affect the three economic sectors considered. On the other hand, as expected, Togo economic sectors fail to affect the world oil price markets, which confirms that Togo, as a small open economy and also a small net oil-importing country, has no pricing power in the world oil markets. However, by using an unrestricted VAR asymmetric impact model proposed by [START_REF] Mork | Oil and the macroeconomy when prices go up and down: an extension of Hamilton's results[END_REF], we find that the impact of world oil price is nonlinear. The positive changes in world oil price do not affect the value-added of the primary, secondary and tertiary sectors while the negative changes in world oil price contribute to improve significantly the value-added of the primary and secondary sectors. The paper finishes the analysis by showing that the primary and secondary sectors affect respectively the tertiary sector. The inverse is not true. This paper recommends that Togo must seek to take benefit from all decreases of world oil price for boosting the value-added of the primary and the secondary sectors, and therefore improve its economic activity.

The remainder of this paper is organized as follows. Section 2 describes the methodology and justify the choice of the data. Section 3 presents and interprets both the results of the linear and nonlinear impact models. Section 4 concludes.

Methodology and data

VAR model

The VAR model has become since the pioneering work of [START_REF] Sims | Macroeconomics and reality[END_REF], one of the best approaches employed in the analysis of the dynamic economic system, especially in the research of interactions, between commodity and energy price shocks and macroeconomy [START_REF] Du | The relationship between oil price shocks and China's macroeconomy: An empirical analysis[END_REF][START_REF] Brown | Energy prices and aggregate economic activity: an interpretative survey[END_REF]. This study adopts an unrestricted vector autoregressive (VAR) model to analyze the effect of world oil price shocks on the three Togo's economic sectors. Consider the following VAR model of order :

= + ∑ -+ = ,
where = , , … , ′ is a × vector of endogenous variables, while -is the corresponding lag terms of order . is the × matrix of autoregressive coefficients of the vector -for = , , … , . = , , … , ′ is the × intercept vector of the VAR model, = , , … , ′ is the × vector of white noise process.

We consider the annual VAR model including the following four endogenous variables: the real oil price, the value-added of the primary sector in the economy, the value-added of the secondary sector in the economy, and the value-added of the tertiary sector in the economy. The effects of oil price shocks on the three sectors are our main objects of interest.

Our data are annual time series and the sample period is 1970-2017. The VAR model is estimated with a constant and two lags, which is determined by the likelihood ratio test, Akaike information criteria, and Schwartz information criteria.

Asymmetric analysis

After analyzing the linear effect of oil price shock on the three considered-sectors (primary, secondary, and tertiary) of Togolese economy, we analyze the nonlinear effect of oil price shock on the three sectors. This could lead us to achieve our specific objective which consists to know whether the effects of oil price shock on Togolese economic sectors are nonlinear or not. Thus, according to [START_REF] Du | The relationship between oil price shocks and China's macroeconomy: An empirical analysis[END_REF], we consider the nonlinear transformation of oil prices developed by [START_REF] Mork | Oil and the macroeconomy when prices go up and down: an extension of Hamilton's results[END_REF], which are widely employed in the empirical literature. [START_REF] Mork | Oil and the macroeconomy when prices go up and down: an extension of Hamilton's results[END_REF] proposes to consider the asymmetric responses by separating the oil price variable into upward (positive changes) and downward (negative changes) movements. [START_REF] Mork | Oil and the macroeconomy when prices go up and down: an extension of Hamilton's results[END_REF] allows an asymmetric response to oil price changes by specifying that the oil price increases and decreases as separate variables. Their definitions are as follows:

Δ + = { Δ Δ > ℎ and Δ -= { Δ Δ < ℎ ,
where Δ is the rate of changes in the world oil price, while Δ + and Δ -are respectively, the positive and negative rate of changes in the world oil price. As the linear case, it is important to test the stationarity of the transformed series by Augmented Dickey-Fuller (ADF) and Phillips-Perron (PP) tests. If the transformed series are stationary, then, we can estimate the VAR model.

Data description

The four endogenous variables considered are constructed as follows:

1. Real oil price ( ). We choose the West Texas Intermediate brent annual crude oil price as the proxy of world oil price. The data are derived from the website of Fred of St-Louis. But, like in [START_REF] Balcilar | The Impact of Oil Price on South African GDP Growth: A Bayesian Markov Switching-VAR Analysis[END_REF], to obtain the real oil price, nominal values of world oil price are deflated using consumer price index (CPI) from the International Financial Statistics (IFS) of the International Monetary Fund (IMF).

2. Primary sector ( . The primary sector tends to make up a larger portion of the economy in developing countries than it does in developed countries. But due to the hard problem of availability of data related to the primary sector of Togo, we choose to quantify the primary sector by the agriculture, forestry, hunting and fishing value-added to the economy in real value.

Secondary sector (

). This sector includes the industrial activity (including construction, aeronautic, electronic). But, due to the lack of data related to the secondary sector of Togo, we choose to quantify the secondary sector by the industry (including construction) value-added to the economy in real value.

Tertiary sector (

). The tertiary sector is also known as the services sector. The service sector produces intangible goods, more precisely services instead of goods, and it comprises various service industries including warehousing and transportation services, information services, securities and other investment services, professional services, waste management; health care and social assistance; and arts, entertainment, and recreation. We quantify the tertiary sector by the services value-added to the economy in real value.

The data are obtained either from the International Monetary Fund-International Financial Statistics (IFS), World Development Indicators (WDI) and Fred of St Louis databases, depending on the availability.

Figure 2 below presents the evolution of the three sectors considered in this analysis and Table 1 summarizes the descriptive statistics. Overall, the service sector is the largest contributor to the Togolese economy. It is followed by the primary sector (agriculture, forestry, fishing) whose contribution has tended over the last five years to outweigh the contribution of the services sector to the economy. Finally, we have the secondary sector (industry) whose added value has been steadily declining over the years. Table 1 shows that oil price is more volatile compared to the value-added of the three sectors because oil price has the highest standard deviation. The four variables are normally distributed according to the Jarque-Bera statistics. This is because the null hypothesis of normal distribution cannot be rejected at 5% significance level for these variables.

Unit root test

According to [START_REF] Sims | Macroeconomics and reality[END_REF], the VAR system can be transformed into its infinite moving average representation only if the variables used in the model are stationary, and the moving average representation is used to obtain both the forecast error variance decomposition and the impulseresponse functions. We test the stationarity of the time series by performing two-unit root tests. More specifically, we perform the Augmented Dickey-Fuller (ADF) test and Phillips-Perron (PP) test. The null hypothesis of the test is that the series has a unit root. If the null hypothesis is rejected, we can conclude that the series is stationary. The results of these two-unit root tests are summarized in Table 2. The first two columns of Table 2 report the test results of the variables in level (log). These results reveal that all the variables are nonstationary in level. Then, in order to get the stationary series, we take the first difference for all the variable in logarithmic form and we perform again the two-unit root tests. The second two columns of Table 2 reporting the test results of each variable (in first-difference) show that the null hypothesis of having a unit root is rejected for all the variables at 1% significant confidence level, meaning that all the first difference variables are stationary. In the following analysis, we use the first difference form of all the series. Concerning the transformed oil price series, the results are reported in Table 3. The results indicate that the transformed series are stationary. Thus, the VAR models can be estimated with these transformed series directly.

Results

In this section, we analyze the empirical results of the VAR model described in the previous section, including the Granger causality test, impulse-response functions and variance decomposition analysis, etc. We focus on the effect of oil price shocks on the three selected sectors of Togolese economy.

Two VAR model are estimated. The first one analyzes the linear effect of world oil price shocks on Togolese's economic sectors and the second one analyzes the nonlinear effects (asymmetric analysis) of world oil price shocks on the Togolese's economic sectors. For the two models, the number of two lags is retained and determined by the likelihood ratio test, Akaike information criteria, and Schwartz information criteria. Then, after estimating the VAR models with two lags on the variables (transformed in first difference), we perform two validation test such us the White test (residual heteroskedasticity, cf. Table 4) and the Residual serial correlation LM Tests (serial correlation, cf. Table 5). The two tests confirm that our residuals are BLUE.

Granger causality test

To investigate the causal relationships of the variables, we perform the Granger causality tests. Granger causality test investigates whether a variable can help forecast another variable . If it cannot, then we say that does not Granger cause . More formally, fails to Granger cause if for all > the mean squared error (MSE) of a forecast of + based on ( , -, …) is the same as the MSE of a forecast of + that uses both ( , -, … and ( , -, …). For more technical representation, please refer to [START_REF] Hamilton | Time Series Analysis[END_REF]. The results of the Granger causality test are summarized in Table 3. We mainly focus on the causal relationship between oil price (change rate) and the three economic sectors, and also between the three economic sectors.

From Table 6, we can see that oil price does not cause any of the value-added of the three sectors (primary, secondary and tertiary) of Togolese economy. The value-added of the three sectors cannot cause the world oil price too (results not summarized in this paper). Concerning the Togolese economy, the Granger causality tests reveal that the value-added of the primary sector and the secondary sector cause respectively at 1% and 5%, the value-added of the tertiary sector.

The secondary sector causes at 10%, the primary sector. This is also true for VAR models on the transformed oil price series.

Impulse-response functions

We resort to impulse-response functions to investigate the dynamic/direct impact of the oil price shocks on the VAR system. Impulse-response functions are used to plot the effects of a shock to one endogenous variable on the other variables in the VAR system. Figure 3 plots the responses of the value-added of the primary, secondary and tertiary sectors to one-unit innovation of the world oil price with two standard error bands. We can see that a linear oil price shock does not affect the value-added of the three considered-sectors of Togolese economy.

Figure 4 plots the responses of primary, secondary and tertiary sectors value-added to the positive oil price shocks (positive change in oil price) on one hand and the responses of the value-added of those sectors to the negative oil price shocks (negative change in oil price) on the other hand. Figure 4 allows to analyze the asymmetric impact of oil price shocks.

Figure 4 shows that positive oil price shock does not affect the value-added of the primary, secondary and tertiary sectors. This leads us to conclude that Togo's economic sectors are resilient to positive changes in world oil price. Our results do not corroborate those of [START_REF] Sachez | Welfare effects of rising oil prices on oil-importing developing countries[END_REF]Berunment et al. (2010). But these results corroborate those of [START_REF] Du | The relationship between oil price shocks and China's macroeconomy: An empirical analysis[END_REF]. This lack of effect and the divergence of results can be explained by the potential subvention mechanisms put in place de jure by the Togolese government to control the effects of oil price fluctuations on Togolese sectors. de facto, no legal published documents are proving the existence of these subvention mechanisms. Encouraging to improve its governance, Togo is a country that communicates very little about its economic policies, especially those related to petroleum products. Togolese government claims to have subsidized the price of oil at the pump for years when the international price of barrel of oil rises so that this does not adversely affect the purchasing power of the Togolese people and the various sectors of the economy. However, when world prices fall, the drop in the price at the pump in Togo is not immediate and proportional to that on the world market. However, since 2010, Togo has a new adjustment mechanism4 that limits the increase in pump prices to 30% per year, independently of the evolution of oil prices on the international market.

However, negative changes in world oil price contribute to improve the added value of the primary and secondary economic sectors. These results corroborate those of Berunment et al. (2010) but do not corroborate those of [START_REF] Wesseh | Exchange rate fluctuations, oil price shocks and economic growth in a small net-importing economy[END_REF] and [START_REF] Liew | Oil Price Shocks and Sectoral Outputs: Empirical Evidence from Malaysia[END_REF]. Besides, our results reveal that negative changes in world oil price do not affect the added value of the tertiary sector. In full modernization and economic transformation, the primary sector in Togo, like that of developing countries, is becoming more mechanized and motorized over the years. Thus, the whole sector has become over the years more dependent on the world oil price.

As in Berunment et al. (2010) and [START_REF] Sachez | Welfare effects of rising oil prices on oil-importing developing countries[END_REF], falls in energy prices, such as oil prices, are an opportunity (lower costs) for the primary sector, while increases in oil price rise production costs. The secondary sector dominated by industrial activity (because of the availability of raw materials), takes benefit from the decrease in energy prices, and therefore from the decrease in the world oil price to improve its value-added and to support its activity. Hence, the positive effect of the negative oil price shock on the value-added of the primary and secondary sectors. The lack of effect of the negative oil price shock on the tertiary sector can be explained by the fact that this sector in Togo, is more dominated by telecommunications, banks and not by transport (maritime, land and air) which requires oil-derived products.

Table 7 shows that negative oil price shocks or negative changes in oil price have a positive cumulative effect on primary and secondary sectors. More specifically, a 100% decrease of Mork (1989)'s transformation of the oil price will cumulatively increase after one year the value-added of primary and secondary sectors by about 4.445% and 4.625%, respectively.

Variance decomposition

Variance decomposition separates the variation in an endogenous variable into the component shocks to the model and provides information about the relative importance of each random innovation in affecting the variables in the VAR model. In this section, we mainly focus on the variance decomposition of the primary, secondary and tertiary sector to see how many of the unanticipated changes of these three variables are explained by each type of different oil price shock (linear, nonlinear oil price shock). The results are presented in Table 8.

The table shows that linear oil price shocks are not a substantial source of volatility in value-added of the primary, secondary and tertiary sectors. However, the value-added of the primary and secondary sectors are substantial sources of volatility in the value-added of the tertiary sector. The contributions of the primary and secondary sectors in the volatility in the value-added of the tertiary sector are 44% and 25%, respectively. We note that the value-added of the secondary sector is at 13%, source of volatility in the value-added of the primary sector. Concerning the asymmetry analysis, Table 8 shows that the negative world oil price shocks or negative changes in the world oil price are substantial sources of volatility in the value-added of the primary and secondary sectors. The contributions of negative changes in world oil price or negative world oil price shocks to the primary and secondary sectors are 13% and 15%, respectively. Overall, although it is often argued that world oil price shocks do not affect Togolese economic sectors, the analysis carried out shows that some effects are hidden. Considering the linear approach, analysis shows that the direct/linear world oil price shock does not affect Togolese economic sectors. However, the adoption of the asymmetric approach proposed by [START_REF] Mork | Oil and the macroeconomy when prices go up and down: an extension of Hamilton's results[END_REF] reveals that world oil price shocks have a nonlinear effect on Togolese economic sectors. Negative changes in world oil price have a positive effect on primary and secondary sectors but do not affect the tertiary sector. However, positive changes in oil price do not affect the three sectors of activity, due to the subvention/subsidy mechanisms put in place by the Togolese government, but not demonstrated.

Conclusion

The impact of the world oil price on the economic activity of net oil-importing countries is becoming more and more significant because of the increase in oil consumption. In this paper, we analyze the direct and asymmetric effects of the world oil price shocks on Togolese economic sectors based on an annual time series from 1970 to 2017. To do so, we have used an unrestricted Vector autoregressive model which is one of the leading approaches employed in the existing literature. In our analysis, we consider the Granger causality test, the impulse response functions and the variance decomposition. As results, as expected, the world oil price is exogenous to Togo' economic sectors. Togo is a small net oil-importing country that cannot affect the world oil price. The impulse-response functions of the linear impact model show that the world oil price shocks do not affect the value-added of the primary, secondary and tertiary sectors of Togolese economy. However, the results of the asymmetric impact model show that the impact of the world oil price on Togo's economic sector is nonlinear. The impact of the positive changes in oil price is not significant on each on the value-added of the three sectors while the negative world oil price shocks increase significantly the value-added of the primary and secondary sectors, but not the tertiary sector. Cumulatively, a 100% decrease in the world oil price, increases the value-added of the primary and secondary sectors by about 4.45% and 4.63%, respectively. The paper finishes the analysis by showing that the value-added of the primary and secondary sectors affect respectively the value-added of the tertiary sector. The inverse is not true. Togo must seek to take benefit from all decreases of world oil price like remarked in March 2020 to boost its economy, and more precisely improve the value-added of primary and the secondary sectors.

Figure 1 :

 1 Figure 1: Energy use (kg of oil equivalent per capita)

Figure 2 :

 2 Figure 2: Sectors of Togolese's economy (1970-2017) Panel a: Agriculture, forestry, and fishing, value added (% of GDP)

Figure 3 :

 3 Figure 3: Impulse-responses functions to linear oil price shocks

Table 1 :

 1 Descriptive statistics (1970Descriptive statistics ( -2017) )

		Oil price	Agriculture,	Industry, value	Service, value
			value added (%	added (% of	added (% of
			of GDP)	GDP)	GDP)
	Mean	43.160	37.027	19.158	41.925
	Std. Dev.	29.184	3.681	2.475	6.512
	Jarque-Bera	4.435	1.669	3.892	2.129
	Probability	0.109	0.434	0.143	0.345
	Observations	48	48	48	48

Source: Author's calculations.

Table 2 :

 2 Unit root test

	Variables	Level (log.)	First difference
		ADF test	PP test	ADF test	PP test
		0.63	0.64	-5.37***	-5.38***
		3.34	5.53	-6.07***	-7.79***
		2.25	-1.69	-6.92***	-7.31***
		-2.01	-1.83	-8.48***	-8.95***
		*, ** and *** denote significant at 10%, 5% and 1% levels, respectively.

Table 3 :

 3 Unit root test of the transformed oil price series

		ADF test	PP test
	Δ + Δ -	-4.66*** -5.63***	-4.63*** -5.97***

*, ** and *** denote significant at 10%, 5% and 1% levels, respectively.

Table 4 :

 4 White test

	Lags

Null Hypothesis: no heteroskedasticity.

Table 5 :

 5 VAR Residual Serial Correlation LM Test

	Lags	LM-Stat	Prob
		Linear impact model (VAR)
	1	28.47	0.287
	2	21.81	0.647
	3	28.78	0.273
	4	25.56	0.431
	Nonlinear impact model (positive changes in oil price)
	1	23.47	0.550
	2	26.22	0.396
	3	25.51	0.434
	4	29.52	0.243
	Nonlinear impact model (negative changes in oil price)
	1	34.06	0.107
	2	22.18	0.625
	3	29.18	0.256
	4	26.53	0.380

Note: Null Hypothesis: no serial correlation at lag order ℎ.

Table 6 :

 6 Granger causality test

	Null hypothesis

Source: Author's estimates.

Table 7 :

 7 Cumulative responses to nonlinear oil price shocks (negative change)

	1 st year	4.445	4.625	-1.809
	2 nd year	4.503	3.589	-0.438
		Source: Author's estimates.

Table 8 :

 8 Variance decomposition

			Linear impact model (oil price shock)
	Variance Decomposition: value-added of the primary sector
	Period	S.E.	Δ
	1	10.94	4.20	95.80	0.00	0.00
	2	12.25	3.39	83.73	12.69	0.18
	Variance Decomposition: value-added of the secondary sector
	1	11.88	10.93	8.33	80.74	0.00
	2	12.10	11.30	8.79	79.85	0.05
	Variance Decomposition: value-added of the tertiary sector
	1	13.43	0.00	38.93	17.05	44.02
	2	16.03	0.19	44.04	24.68	31.09

VAR model with transformed oil price series (positive change in oil price)

	Variance Decomposition: value-added of the primary sector Period S.E. Δ +
	1	12.70	0.86	99.14	0.00	0.00
	2	13.95	0.70	84.36	14.44	0.50
	Variance Decomposition: value-added of the secondary sector
	1	10.36	2.59	17.72	79.68	0.00
	2	11.65	2.94	17.54	79.43	0.09
	Variance Decomposition: value-added of the tertiary sector
	1	13.13	3.09	38.03	24.23	34.65
	2	13.40	2.34	41.20	30.30	26.16

VAR model with transformed oil price series (negative change in oil price)

	Variance Decomposition of Primary sector Period S.E. -
	1	16.55	15.82	84.18	0.00	0.00
	2	18.10	12.61	76.34	11.03	0.03
	Variance Decomposition of Secondary sector
	1	11.17	15.05	4.00	80.96	0.00
	2	12.52	15.22	5.15	79.59	0.03
	Variance Decomposition of Tertiary sector
	1	13.12	1.81		35.11	15.04	48.03
	2	13.63	1.94	42.59 Source: Author's estimates	22.70	32.77

Indeed, Mork (1989) in his study on the role of oil price shocks on economic activity, finds oil price increases to affect economic growth negatively while a decline in oil price does not have the opposite effect.

South Africa's economy is the most transformed in Sub-Saharan Africa. Source: African Development Group.

[START_REF] Liew | Oil Price Shocks and Sectoral Outputs: Empirical Evidence from Malaysia[END_REF], through the estimation of a nonlinear ARDL model on Malaysia, show that nonlinear oil price-out nexus is revealed for the manufacturing and industrial outputs. They find that rising output is significantly associated to oil price increase, while to falling output is significantly related to oil price decrease.

source: Abdoulaye Tchagawou, Revue de presse, Semaine No.

50, Press Service of the French Embassy in Togo, December, 2010.


Jean-Guillaume Dumas
email: jean-guillaume.dumas@imag.fr.

Laurent Fousse
email: laurent.fousse@imag.fr.

Bruno Salvy
email: bruno.salvy@inria.fr..

Simultaneous Modular Reduction and Kronecker Substitution for Small Finite Fields *

We present algorithms to perform modular polynomial multiplication or modular dot product efficiently in a single machine word. We pack polynomials into integers and perform several modular operations with machine integer or floating point arithmetic. The modular polynomials are converted into integers using Kronecker substitution (evaluation at a sufficiently large integer). With some control on the sizes and degrees, arithmetic operations on the polynomials can be performed directly with machine integers or floating point numbers and the number of conversions can be reduced. We also present efficient ways to recover the modular values of the coefficients. This leads to practical gains of quite large constant factors for polynomial multiplication, prime field linear algebra and small extension field arithmetic.

Introduction

While theoretically well understood, the basic routines for linear algebra or polynomial arithmetic over finite fields are difficult to implement efficiently. Important factors of speed can be gained by exploiting machine integer or floatingpoint arithmetic and taking cache effects into account. This has been demonstrated for instance by [START_REF] Dumas | Finite field linear algebra subroutines[END_REF][START_REF] Dumas | Efficient dot product over finite fields[END_REF] who wrapped cache-aware routines for efficient small finite field linear algebra in their FFLAS/FFPACK project.

The elements of small enough prime fields can be represented as integers fitting in a machine word or half-word, or in exact floating point numbers. For extension fields, the elements can be represented as polynomials over a prime 1 field. A further compression is proposed by [START_REF] Dumas | Finite field linear algebra subroutines[END_REF][START_REF] Dumas | Q-adic transform revisited[END_REF], using Kronecker substitution (Gathen and Gerhard, 1999, §8.4): the polynomials are represented by their value at an integer q larger than the characteristic of the field. We call DQT, for Discrete Q-adic Transform, this simple map:

DQT : Z/pZ[X] → Z k i=0 α i X i → k i=0 α i q i .
(

) 1
With some care, in particular on the size of q, it is possible to map the operations in the extension field into the floating point arithmetic realization of this q-adic representation.

In computational number theory, matrices over F 2 are often compressed by fitting several entries into one via the binary representation of machine integers [START_REF] Coppersmith | Solving linear equations over GF (2): block Lanczos algorithm[END_REF][START_REF] Kaltofen | Distributed matrix-free solution of large sparse linear systems over finite fields[END_REF]. The need for efficient matrix computations over very small finite fields also arises in other areas and in particular graph theory (adjacency matrices), see e.g., [START_REF] May | Efficient matrix rank computation with application to the study of strongly regular graphs[END_REF] or [START_REF] Weng | Pseudo-Paley graphs and skew Hadamard difference sets from presemifields[END_REF]. In these cases, one is thus led to work with polynomials as in (1), using a small q.

Recovering the results is obtained by an inverse DQT. This consists of a radix conversion (Gathen and Gerhard, 1999, Algorithm 9.14), followed by a reduction of the coefficients modulo p. On modern processors, machine division or remaindering, that are used in radix conversion, are quite slow when compared to other arithmetic operations 1 .

In this article, we use two techniques to avoid some of these remainderings: (i) delayed reduction; (ii) simultaneous reduction. Delayed reduction means that after having been transformed into their DQT form, polynomials may undergo several arithmetic operations before being converted back into coefficient form. This is possible as long as the coefficients remain smaller than q, which is prevented by a priori bounds. Simultaneous reduction is the operation of recovering α i mod p from α i , for i = 1, . . . , k, in the context of the inverse DQT. We propose a new algorithm called REDQ performing these k modular reductions by a single division, ⌈ k+1 2 ⌉ additions and multiplications and some table look-up. We also discuss the possibility of replacing the remaining division by floating point operations, taking into account the rounding modes.

We recall in Section 2 the Kronecker substitution and delayed reduction algorithms. Then we present our new simultaneous reduction algorithm and give its complexity in Section 3 and we discuss how to replace the remaining machine division by floating point operations with different rounding modes in Section 4. Our new reduction algorithm has two parts. The first one is a compression, performed by arithmetic operations, which reduces the size of the polynomial entries. The second one is performed only when required and is a correction, which gives to the residues their correct value modulo p when the compression has shifted the result. We also present a time-memory trade-off enabling some table look-ups computing this correction. Then we apply the DQT to different contexts: modular polynomial multiplication in Section 5; linear algebra over small extension fields in Section 6; compressed linear algebra over small prime fields in Section 7. This gives some constraints on the possible choices of q and k. In the three applications anyway, we show that these compression techniques represent a speed-up factor usually of the order of the number k of residues stored in the compressed format.

Preliminary versions of this work have been presented by [START_REF] Dumas | Q-adic transform revisited[END_REF] and [START_REF] Dumas | Compressed modular matrix multiplication[END_REF]. Here, we give an improved version of the simultaneous reduction where the number of operations has been divided by two. We also give a complete study of the behavior of the division of integers by floating point routines, depending on the rounding modes. Finally, we present more experimental results and faster implementations of the applications, namely a Karatsuba version of the polynomial multiplication and a right compressed matrix multiplication.

Q-adic Representation of Polynomials

Kronecker Substitution

The principle of Kronecker substitution is very simple. It consists in evaluating polynomials at a given integer as in Eq. (1). For instance, for k = 2, the substitution is performed by the following compression: double& init3( double& r, const double u, const double v, const double w) { // _dQ is a floating point storage of Q r=u; r*=_dQ; r+=v; r*=_dQ; return r+=w; }

The integer q can be chosen to be a power of 2 in a binary architecture. Then the Horner like evaluation of a polynomial at q is just a left shift. One can then compute this shift with exponent manipulations in floating point arithmetic and use native shift operators (e.g., the < < operator in C) as soon as values are within the 32 (or 64 when available) bit range.

The motivation for this substitution is its use in multiplication.

Example 1. To multiply a = X + 1 by b = X + 2 in Z/3Z[X] one can use the substitution X = q := 100; compute 101 × 102 = 10302; use radix conversion to write 10302 = q 2 + 3q + 2; reduce the coefficients modulo 3 to get

a× b = X 2 + 2. More generally, if p is prime, a = k-1 i=0 a i X i and b = k-1 i=0 b i X i are two polynomials in Z/pZ[X], then one can perform the polynomial multiplication ab via Kronecker substitution. The product of ã = k-1 i=0 a i q i and b = k-1 i=0 b i q i is given by ab = 2k-2 j=0 j i=0 a i b j-i q j .
(2)

Now if q is large enough, the coefficient of q j does not exceed q -1. If moreover k is not too large, the product fits in a machine number (floating point number or integer). Thus in that case, it is possible to evaluate ã and b as machine numbers, compute the product of these evaluations, and convert back to polynomials by radix conversion. There just remains to perform reductions of the coefficients modulo p.

We show in Section 6 that one can also tabulate the evaluations at q, and that one can access directly the required part of the machine words (using e.g., bit fields and unions in C) instead of performing a radix conversion.

Delayed Reduction

With current processors, machine division (as well as modular reduction) is still much slower in general than machine addition and machine multiplication. It is possible to replace the machine division by some other implementation, such as floating point multiplication by the inverse [START_REF] Shoup | NTL 5.4.1: A library for doing number theory[END_REF] or Montgomery reduction [START_REF] Montgomery | Modular multiplication without trial division[END_REF]; see e.g., [START_REF] Dumas | Efficient dot product over finite fields[END_REF] and references therein for more details.

It is also important to reduce the number of machine remainderings when performing modular computations. This is achieved by using Kronecker substitution seldom and perform as many arithmetic operations as possible on the transformed values.

Discrete Q-adic Transform

The idea of the Q-adic transform is to combine Kronecker substitution with delayed reduction.

We call DQT the evaluation of polynomials modulo p at a sufficiently large q. Since this is a ring morphism, products as well as additions can be performed on the transformed values. We call DQT inverse the radix conversion of a q-adic expansion followed by a modular reduction. For a given procedure ⋆(a 1 , . . . , a m ) performing only ring operations, the idea is to compute the DQT of the a i 's, perform ⋆ on these DQT's and compute the inverse DQT at the end. For appropriate choices of the parameters, this procedure recovers the exact value.

As an example, we recall the dot product of [START_REF] Dumas | Finite field linear algebra subroutines[END_REF] in Algorithm 1.

The following bounds generalize those of (Gathen and Gerhard, 1999, §8.4): Theorem 2. [START_REF] Dumas | Finite field linear algebra subroutines[END_REF] Let β be the number of mantissa bits available within the machine numbers. If

q > nk(p -1) 2 and (2k -1) log 2 (q) ≤ β, (3)
Algorithm 1 Polynomial dot product by DQT Input: Two vectors of polynomials v 1 and v 2 in Z/pZ[X] n of degree less than k; Input: a sufficiently large integer q. Output: R ∈ Z/pZ[X], with R = v T 1 • v 2 . Polynomial to q-adic conversion 1: Set v 1 and v 2 to the floating point vectors of the evaluations at q of the elements of v 1 and v 2 . {Using e.g., Horner's formula} One computation

2: Compute r = v 1 T • v 2
Building the solution

3: r = 2k-2 i=0 µ i q i . {Using radix conversion} 4: For each i, set µ i = µ i mod p 5: Return R = 2k-2 i=0 µ i X i
then Algorithm 1 is correct. Dumas et al. (2002, Figures 5 & 6) show that this wrapping is already a pretty good way to obtain high speed linear algebra over small extension fields. They reach high peak performance, quite close to those obtained with prime fields, namely 420 Millions of finite field operations per second (Mop/s) on a Pentium III, 735 MHz, and more than 500 Mop/s on a 64-bit DEC alpha 500 MHz. This is roughly 20% below the pure floating point performance and 15% below the prime field implementation. We show in Section 6 that the new algorithms presented in this article enable to reduce this overhead to less than 4%.

REDQ: Modular Reduction in the DQT Domain

The first improvement we propose to the DQT is to replace the costly modular reduction of the polynomial coefficients (e.g., in Step 4 of Algorithm 1) by a single division by p followed by several shifts. In the next section, we replace this division by a multiplication by an inverse.

Examples

We first illustrate the basic idea on a simple example.

Example 3. Let a = X 2 + 2X + 3 and b = 4X 2 + 5X + 6 be unreduced modulo 5, so that we want to compute a × b = 4X 4 + 3X 3 + 3X 2 + 2X + 3. We are free to choose the integer q at which the evaluation takes place in such a way that shifting by powers of q is cheap. For clarity we take here q = 10000. Thus, we start from r := ã × b = 40013002800270018 for which we need to reduce five coefficients modulo 5. In the direct approach, the coefficients would be recovered by computing 4 = 0 × 5 + 4, 13 = 2 × 5 + 3, 28 = 5 × 5 + 3, 27 = 5 × 5 + 2, 18 = 3 × 5 + 3.

The trick is that we can recover all the quotients at once. Indeed, we compute s := ⌊r/p⌋ = 08002600560054003 that contains all the quotients (in boldface).

The remainders (the coefficients) can then be recovered as u i := ⌊r/q i ⌋-p⌊s/q i ⌋, for i = 0, . . . , 4.

Some more work is needed in order to make this idea correct in general.

Example 4. We now consider the polynomial R = 1234X 3 +5678X 2 +9123X + 4567, the prime p = 23 and use q = 10 6 . Note that this time, p does not divide q.

The polynomial we want to recover is R mod p = 15X 3 + 20X 2 + 15X + 13. We start from r = 1234005678009123004567 and the division gives s := ⌊r/p⌋ = 53652420783005348024. Proceeding as before, we get

u 0 = 15, u 1 = 8, u 2 = 18, u 3 = 15.
These coefficients are small, but they are not the correct ones except for u 3 . Indeed, if C = αp + u i , then Cq = u i q = 0 mod p so that each coefficient needs to be corrected to take into account the values of the preceding ones. We thus consider this first computation as a compression stage and use a correction stage that produces the correct values. The correction is obtained by taking µ 3 = u 3 and µ i = u i -qu i+1 mod p for i = 0, 1, 2 and returning the µ i 's.

Algorithm

Theorem 5. Algorithm REDQ is correct.

The following lemma is probably classical. We state and prove it here for completeness. PROOF. [of Theorem 5] First we prove that 0 ≤ u i < p. Let T i = ⌊r/q i ⌋.

By Lemma 6 ⌊s/q i ⌋ = ⌊T i /p⌋, so that u i = T i -p⌊T i /p⌋ which proves the inequalities.

Algorithm 2 REDQ Input: Two integers p and q satisfying Conditions (3); Input: r = d i=0 µ i q i ∈ Z. Output: ρ ∈ Z, with ρ = d i=0 µ i q i where µ i = µ i mod p.

REDQ COMPRESSION

1: s = r p ; 2: for i = 0 to d do 3:

u i = r q i -p s q i ; 4: end for REDQ CORRECTION {when p ∤ q} 5: µ d = u d 6: for i = 0 to d -1 do 7: µ i = u i -qu i+1 mod p; 8: end for 9: Return ρ = d i=0 µ i q i ;
Next we compute the value of u i by the following sequence of identities modulo p:

u i = T i = r q i = d j=i μj q j-i = d j=i µ j q j-i mod p.
When q = 0 mod p, we thus have u i = µ i mod p and the proof is complete. Otherwise, in view of the previous identity, the correction step on Line 7 produces

u i -qu i+1 = d j=i µ j q j-i -q d j=i+1
µ j q j-i-1 = µ i mod p.

Definition 7. We call REDQ k a simultaneous reduction of k residues performed by Algorithm 2 (in other words k = d + 1 if d is the degree of the Kronecker substitution).

Binary Case

When q is a power of 2 and elements are represented using an integral type, division by q i and flooring are a single operation, a right shift, or direct bit fields extractions when available. Moreover, the remainders can be computed independently and thus the loop of REDQ COMPRESSION can be performed by only half of the required k axpy (combined addition and multiplication, or fused-mac) as shown below:

Proposition 8. Let q be a power of two. Then, a REDQ k COMPRESSION requires ⌈ k+1 2 ⌉ machine word multiplications and additions. PROOF. We use the notations of Algorithm 2. Let b q = log 2 (q) be the number of bits of q = 2 bq and let β be the number of bits in the mantissa of a machine word. If the REDQ k is correct, we have q k ≤ 2 β , or b q k ≤ β. The first value u 0 = r -s × p requires the whole mantissa. Now both ⌊ r q i ⌋ and ⌊ s q i ⌋ can be stored on at most k × b q -i × b q bits. Moreover, by Theorem 5, 0 ≤ ⌊ r q i ⌋ -p × ⌊ s q i ⌋ < ⌊ r q i ⌋ so that the result does not overflow these (k -i)b q bits. This proves that the operations can be done independently on the different parts of the machine words. Now, the total number of bits required to perform the loop of REDQ COMPRESSION is

b q k-1 i=0 k -i = b q k(k + 1) 2 ≤ β k + 1 2 .
This, combined with the non-overlapping of the operations, proves the proposition.

Here is an instance of a REDQ compression with 3 residues in C++, using the syntactic sugar of bit field extraction: inline void REDQ_COMP(UINT32_three& res, const double r, // to be reduced const double p){ // modulo _ULL64_unions rll, tll; // union of 64, 17/34 or 34/17 bits rll._64 = static_cast<UINT64>( r ); tll._64 = static_cast<UINT64>( r/p ); // One division res.high = static_cast<UINT32>(rll._64-tll._64*p); // One axpy rll._17_34.low = rll._34_17.high; // Packing tll._17_34.low = tll._34_17.high; // Packing rll._64 -= tll._64*p; // Two axpy in one res.mid = static_cast<UINT32>(rll._17_34.high); res.low = rll._17_34.low; } In general also the algorithm is efficient because one can precompute 1/p, 1/q, 1/q 2 etc. The computation of each u i and µ i can also be pipelined or vectorized since they are independent. As is, the benefit when compared to direct remaindering by p is that the corrections occur on smaller integers. Thus the remaindering by p can be faster. Actually, another major acceleration can be added: the fact that the µ i are much smaller than the initial μi makes it possible to tabulate the corrections as shown next.

Matrix Version of the Correction

In Example 4, the final correction can be written as a matrix vector product:

µ =     1 17 0 0 0 1 17 0 0 0 1 17 0 0 0 1     u mod p.
More generally, the corrections of lines 5 to 8 of Algorithm 2 are given by a matrix-vector multiplication with an invertible matrix Q k :

Q k =           1 -q 0 . . . 0 0 . . . . . . . . . . . . . . . . . . . . . . . . 0 . . . . . . . . . -q 0 . . . . . . 0 1           .

Tabulations of the Matrix-Vector Product and Time-Memory Trade-off

Tabulating fully the multiplication by Q k requires a table of size at least p k+1 . However, the matrices Q k can be decomposed recursively into smaller similar matrices as follows:

1 0 0 = Q i+j Q i Q j
It is therefore easy to tabulate the product by Q k with an adjustable table of size p j . Proposition 9. Let q be a power of 2. Given a table of size p j (1

≤ j ≤ k + 1), a REDQ k CORRECTION can be performed with ⌊(k -1)/(j -1)⌋ table accesses.
When q is a power of 2, the computation of the u i in the first part of Algorithm 2 requires 1 div and (k + 1)/2 axpy as shown in Proposition 8. Now, the time memory trade-off enables to compute the second part efficiently.

Example 10. We compute the corrections for a degree 6 polynomial. One can tabulate the multiplication by Q 6 , a 7 × 7 matrix, or actually, by only the first 6 rows of Q 6 , denoted by Q 6 , with therefore p 7 entries, each of size 6 log 2 (p). Instead, one can tabulate the multiplication by

Q 2 , a 2 × 3 matrix. To compute [µ 0 , . . . , µ 6 ] T = Q 6 [u 0 . . . , u 6 ] T = [Q 6 [u 0 . . . , u 6 ] T , u 6 ]
it is sufficient to use three multiplications by Q 2 as shown in the following algorithm:

Algorithm 3 Q 6 with an extra memory of size p 3 Input: [u 0 . . . , u 6 ] ∈ (Z/pZ) 7 ; Input: a table Q 2 of the associated 2×3 matrix-vector multiplication over Z/pZ. Output: [µ 0 , . . . , µ 6 ] T = Q 6 [u 0 . . . , u 6 ] T . 1: a 0 , a 1 = Q 2 [u 0 , u 1 , u 2 ] T ; 2: b 0 , b 1 = Q 2 [u 2 , u 3 , u 4 ] T ; 3: c 0 , c 1 = Q 2 [u 4 , u 5 , u 6 ] T ; 4: Return [µ 0 , . . . , µ 6 ] T = [a 0 , a 1 , b 0 , b 1 , c 0 , c 1 , u 6 ] T ;
Note on Indexing. In practice, indexing by a tuple of integers mod p is made by evaluating at p, as u i p i . If more memory is available, one can also directly index in the binary format using u i 2 ⌈log 2 (p)⌉ i . On the one hand all the multiplications by p are replaced by bit fields extraction. On the other hand, this makes the table grow a little bit, from p k to 2 ⌈log 2 (p)⌉k .

In the rest of this article we restrict to the case when q is a power of two.

Euclidean Division by Floating Point Routines

In the implementation of REDQ above, there remains one machine division (Line 1). Since exact division is rather time-consuming on modern processors, it is natural to try and replace it by floating point operations. If r and p are integers and we want to compute r/p, then computing r/p by a floating point division with a rounding to nearest mode followed by flooring produces the expected value (Lefèvre, 2005a, Theorem 1).

Instead of a division, a multiplication by a precomputed inverse of p can be used (this is done e.g., in NTL [START_REF] Shoup | A computational introduction to number theory and algebra[END_REF], Theorem 1.5)). However, in that case, the correctness is not guaranteed for all representable r (see Lefèvre (2005b) for details). Thus, as is done e.g. in NTL, one has to make some additional tests and corrections. In this section we propose bounds on r for which this correctness is guaranteed, taking the rounding modes into account. Outside of these bounds, we show that a difference of (at most) 1 after the flooring is possible for some values of r, p and the rounding modes. This can be detected by two tests on the resulting residue (below 0 or above p) as is done by [START_REF] Shoup | NTL 5.4.1: A library for doing number theory[END_REF]. We show that only one of these tests is mandatory if the rounding modes can be mixed. The inverse of the prime can be precomputed for each mode, which avoids the need for costly changes of modes. This is summarized in Algorithm 5 at the end of this section.

Rounding Modes

We assume that the floating-point arithmetic follows the IEEE 754 standard and we denote by ulp(z) the unit in the last place of z for a floating-point number z such that 2 β-1 ≤ |z| ≤ 2 β -1.

For each operation three rounding modes are possible (up △(•), down ▽(•) and nearest2 ⋄(•)). If two computations take place at different times, each of them can be performed in a different rounding mode, so that we have 9 cases to consider. It is worth investigating all cases since changing rounding modes is a costly operation and the FPU may be in a particular rounding mode due to constraints in the surrounding code, or some rounding modes are simply not available in the particular computing environment. Also, computing the best bounds enables to make the best of delayed reduction in modular computations.

Floating Point Division

Denoting by r = kp+ u the Euclidean division of r by p where 0 ≤ r ≤ 2 β -1 we are interested to know under which conditions on r and p Algorithm 4 returns the quotient k, depending on the rounding modes

• 1 and • 2 . Algorithm 4 FDIV Input: One integer r such that 0 ≤ r ≤ 2 β -1; Input: one integer p such that 1 ≤ p ≤ 2 β -1; Input: two choices of rounding-modes • 1 and • 2 . Output: ⌊ r p ⌋. 1: invp ← • 1 (1/p) 2: x ← • 2 (r • invp) 3: Return ⌊x⌋.

Results

Our results are summarized in Table 1.

The column "Range" gives guaranteed bounds on the result. This shows which tests and corrections may be needed in order to compute the expected value k. The interval given in this column is optimal, in the sense that we have examples where ⌊x⌋ = k in each possible direction.

The column "Bound on r" gives a strict upper bound on r under which ⌊x⌋ ≤ k, in those cases where the result could indeed overflow. We do not know whether these bounds are optimal; for some cases we could find a systematic family of examples indexed by β that reach the bounds asymptotically; other bounds could be approached by exhaustive search on small β. These examples are not included here. We believe that all bounds are optimal except for case 2 where a value closer to 3 8 (instead of 1 3 ) could be reached (take β = 2n + 1 and p = 2 n -1, r = (3 • 2 n-2 + 3)p -1). For our purpose 1 3 is close enough. The column "Lost bits" gives a simplified version of this bound: if r fits in β minus this number of bits, then ⌊x⌋ ≤ k.

Case • 1 • 2 Range Bound on r Lost bits 1 △(•) △(•) k ≤ ⌊x⌋ ≤ k + 1 2 β /(4 + 2 2-β ) 3 2 △(•) ⋄(•) k ≤ ⌊x⌋ ≤ k + 1 2 β /(3 + 2 1-β ) 2 3 △(•) ▽(•) k ≤ ⌊x⌋ ≤ k + 1 2 β /2 1 4 ⋄(•) △(•) k ≤ ⌊x⌋ ≤ k + 1 2 β /(3 + 2 1-β ) 2 5 ⋄(•) ⋄(•) k -1 ≤ ⌊x⌋ ≤ k + 1 2 β-1 /(1 + 2 -1-β ) 2 6 ⋄(•) ▽(•) k -1 ≤ ⌊x⌋ ≤ k - 0 7 ▽(•) △(•) k -1 ≤ ⌊x⌋ ≤ k + 1 2 β /2 1 8 ▽(•) ⋄(•) k -1 ≤ ⌊x⌋ ≤ k - 0 9 ▽(•) ▽(•) k -1 ≤ ⌊x⌋ ≤ k - 0
Table 1: Possible values of ⌊x⌋ and bounds on r such that ⌊x⌋ ≤ k

Proof of the Bounds on ⌊x⌋

We denote by ǫ 1 and ǫ 2 the errors in rounding as follows:

invp = 1 p (1 + ǫ 1 ), x = (r • invp)(1 + ǫ 2 ).
Thus the value that is computed is

x = r p (1 + ǫ 1 )(1 + ǫ 2 ) = k + u p + (ǫ 1 + ǫ 2 + ǫ 1 ǫ 2 ) r p =: k + R (4)
where R is the term we want to bound. For example R < 1 means ⌊x⌋ ≤ k.

Bounds on ǫ 1 and ǫ 2 depend on the rounding mode, but in all cases we have

|ǫ i | ≤ 2 1-β , i ∈ {1, 2}. (5)
Lemma 1. The result of Algorithm 4 is never off by more than one unit, that is

k -1 ≤ ⌊x⌋ ≤ k + 1.
Proof. First we show R < 2, which gives the upper bound. This is obtained by injecting the inequalities (5), in the definition (4) of R:

R ≤ p -1 p + (2 2-β + 2 2-2β ) r p (6) ≤ 1 - 1 p + (2 2-β + 2 2-2β ) 2 β -1 p = 1 - 1 p (3 -2 2-2β ).
Thus R < 2 for p ≥ 3. In the special case p = 2, we have ǫ 1 = 0 so that the result still holds.

Similarly, in the other direction we have

R ≥ -(2 2-β + 2 2-2β ) r p ≥ - 1 p (4 -2 2-2β )
and R ≥ -1 for p ≥ 4. For p = 2, the results follows from ǫ 1 = 0. For the last case, p = 3, we first analyze more precisely the rounding error ǫ 1 . The binary expansion of 1 3 = 0.01010101010 . . . implies that

•( 1 3 ) =      1 3 (1 + 2 -β-1 ) if β is even and • (•) ∈ {△(•), ⋄(•)}, 1 3 (1 + 2 -β-1 ) if β is odd and • (•) ∈ {⋄(•), ▽(•)}, 1 3 (1 + 2 -β ) otherwise.
Thus in all cases |ǫ 1 | ≤ 2 -β . Using this better bound in the computation above completes the proof.

Lemma 2. In cases 1, 2, 3 and 4, ⌊x⌋ ≥ k.

Proof. In the first three cases, invp is rounded up, and thus

r • invp = k(p • invp) + u • invp ≥ k.
This implies that •(r • invp) ≥ k because rounding modes are monotone and k is exactly representable.

In case 4, since

|ǫ 1 | ≤ 2 -β , we have invp > (1 -2 -β ) 1 p . Then r • invp = (kp + u)invp ≥ kp • invp > k(1 -2 -β ).
Again, k is an integer thus exactly representable. Denote by k -the largest floating-point number that is strictly less than k:

k -= k -1 2 ulp(k) if k is a power of 2, k -ulp(k)
otherwise.

If k is a power of 2, then the way of rounding k(1 -2 -β ) is the same as that of 1-2 -β since k only changes the exponent in the result. Since △(1-2 -β +δ) ≥ 1 for any δ > 0 we get x = △(r • invp) ≥ k. If k is not a power of 2, then

k -= k -ulp(k) < k(1 -2 -β ) < r • invp and thus again x = △(r • invp) ≥ k.
Lemma 3. In cases 6, 8 and 9, ⌊x⌋ ≤ k.

Proof. In case 6, we have

|ǫ 1 | ≤ 2 -β . This implies r • invp = r p (1 + ǫ 1 ) ≤ r p (1 + 2 -β ) < r p + 1 p < k + 1 and therefore x = ▽(r • invp) < k + 1, whence ⌊x⌋ ≤ k.
In cases 8 and 9, since invp = ▽( 1 p ), we have

r • invp ≤ r p ≤ k + 1 - 1 p . (7)
In case 9, this is sufficient to conclude. Otherwise, as in case 4 above, denote by (k + 1) -the floating-point number that is just below k + 1 and let m be the midpoint between k + 1 and (k + 1) -. If we show r • invp < m, then ⋄(r • invp) ≤ (k + 1) -and therefore ⌊x⌋ ≤ k.

First assume that k + 1 is a power of 2. In this case ulp

(k + 1) = 2 1-β (k + 1) and m = k + 1 -1 4 ulp(k + 1) = (k + 1)(1 -2 -β-1
). Since by assumption p < 2 β and r < 2 β , so p(k + 1) ≤ r + p < 2 β+1 and it follows that

1 p > (k + 1)2 -β-1 .
Together with (7) this gives

r • invp < (k + 1)(1 -2 -β-1 ) = m.
Assume now that k + 1 is not a power of 2 (and k = 0). Then ulp

(k + 1) = ulp(k) ≤ 2 1-β k, m = k + 1 -1 2 ulp(k) and 1 p > k2 -β ≥ ulp(k)/2. Thus finally, r • invp < k + 1 - 1 2 ulp(k) ≤ m. If k = 0 then 1 2 ulp(k + 1) = 2 -β < 1
p and the result follows.

Bounds on r making the result exact

In cases 6, 8 and 9, the result of Algorithm 4 is smaller than k for all values of r. We now proceed to a case-by-case proof of each remaining row of Table 1, giving bounds on r such that this happens.

Case 1. In this case 0 ≤ ǫ i ≤ 2 1-β for i ∈ {1, 2} and thus (6) becomes

R ≤ p -1 p + (2 2-β + 2 2-2β ) r p so that |R| < 1 is implied by r < 2 β 1 4 + 2 2-β
that is close to 2 β /4 and we lose less than 3 bits compared to the bound r < 2 β required to have no loss of precision on r.

Case 2. In this case |ǫ 2 | ≤ 2 -β and (6) becomes

R ≤ p -1 p + (3 • 2 -β + 2 1-2β ) r p
The condition R < 1 is implied by

3 • 2 -β + 2 1-2β < 1 r that is r < 2 β /(3 + 2 1-β
) which is close to 2 β /3, and less that two bits are lost.

Case 3.

In this case -2 1-β ≤ ǫ 2 ≤ 0, |ǫ 1 + ǫ 2 | ≤ 2 1-β and ǫ 1 ǫ 2 ≤ 0. We get R ≤ 1 - 1 p + r p 2 1-β
and the condition R < 1 is ensured by

r < 1 2 2 β
which means we lose one bit.

Case 4. This is as in case 2 since △(•) and ⋄(•) play a symmetric role in the analysis.

At this stage, we have obtained the following.

Proposition 11. In Cases 1-4, Algorithm 4 is correct for r obeying the bounds of Table 1.

The remaining cases are proved similarly:

Case 5. We have

|ǫ 1 | ≤ 2 -β , |ǫ 2 | ≤ 2 -β and |ǫ 1 + ǫ 2 + ǫ 1 ǫ 2 | ≤ 2 1-β + 2 -2β . Then (6) becomes R ≤ 1 + 1 p (1 -2 -β -2 -2β ). The condition R < 1 is implied by r < 2 β-1 /(1 + 2 -1-β ) which is close to 1 2 2 β .
Case 7. The bound on r follows from case 3 as ǫ 1 and ǫ 2 play a symmetric role in the error analysis of case 3.

Using Algorithm 4

Algorithm 5 Applied FDIV Input: r, integer such that 0 ≤ r ≤ 2 β -1; Input: p, integer such that 1 ≤ p ≤ 2 β -1; Output: ⌊ r p ⌋. Constants end if 14: end if 15: Return y.

1: B △ ← 2 β /(3 + 2 1-β ) 2: B ⋄ ← 2 β /(3 + 2 1-β ) 3: B ▽ ← 2 β-1 Precomputation 4: invp △ ← ⋄(1/p) 5: invp ⋄ ← △(1/p) 6: invp ▽ ← △(1/p) Division 7: x ← •(r • invp • ) 8: y ← ⌊x⌋
Algorithm 5 demonstrates how to apply the results of Table 1 in a program. We precompute 1/p in several rounding modes and use the best version in the subsequent multiplication, depending on the current rounding mode •(•). The strategy used here is to make sure ⌊x⌋ ≥ k after the multiplication so only one test at most is needed for the correction. In addition we choose the version that maximizes the bound B.

It should be noted that there is no strategy in the choice of • 1 (•) that guarantees ⌊x⌋ ≤ k for any choice of • 2 (•) (take • 2 (•) = △(•)), meaning that the described strategy is indeed the only one than minimizes the number of tests needed for the correction.

In a typical application of Algorithm 5 where r is the accumulation of several products modulo p, the bound B can be interpreted in the numbers of operations that can be performed before a reduction is necessary. If this number is not exceeded then the correction is never needed.

5 Application 1: Polynomial Multiplication

Delayed Reduction

A classical technique for modular polynomial multiplication is to use only delayed reductions. The idea is to compute each polynomial coefficient by a delayed dotproduct e.g., as in [START_REF] Dumas | Efficient dot product over finite fields[END_REF]): products of the form i a i b k-i are accumulated without reductions, and the overflow is dealt with in one final pass. Thus, with a centered representation modulo p for instance (integers from (1 -p)/2 to (p -1)/2), it is possible to accumulate at least n d products as long as

n d (p -1) 2 < 2 β+1 . (8)
The final modular reduction can be performed in many different ways (e.g., classical division, floating point multiplication by the inverse, Montgomery reduction, etc.), we just call the best one REDC here. At worst, it is equivalent to 1 machine division.

Fast Q-adic Transform

We represent modular polynomials of the form P = N i=0 a i X i by P = P i X d+1 i where the P i 's are polynomials of degree d stored in a single integer in the q-adic way.

Then a product P Q has the form

P Q = P i Q t-i X d+1 t ,
where each multiplication P i Q t-i is computed by Algorithm 1 on a single machine integer. The final reduction is performed by a tabulated REDQ and can also be delayed as long as conditions (3) are guaranteed.

Comparison

We use the following complexity model: multiplications and additions in the field are counted as atomic operations while the machine divisions are counted separately. For instance, we approximate REDC by one machine division and an axpy. We recall that REDQ k denotes a simultaneous reduction of k residues. In this complexity model a REDQ k thus requires 1 division and k/2 multiplications and additions. A REDQ with k residues is a Kronecker substitution with a polynomial of degree d = k -1. We also call k-FQT the use of a polynomial of degree d = k-1 for the q-adic substitution. Thus a multiplication P i Q t-i in a k-FQT requires the reduction of 2d+1 coefficients, i.e., a REDQ 2d+1 =REDQ 2k-1 .

Let P be a polynomial of degree N with indeterminate X. If we use a (d+1)-FQT, it will then become a polynomial of degree D q in the indeterminate Y = X d+1 , with

D q = N + 1 d + 1 -1.
Table 2 gives the respective complexities of both strategies. The values of n d and n q are given by Equations (3,8):

n d = 2 β+1 (p -1) 2 , n q = q (d + 1)(p -1) 2 with q = 2 m 2d+1 . Mul & Add Reductions Delayed (2N + 1) 2 (2N + 1) 2N +1 n d REDC d-FQT (2D q + 1) 2 (2D q + 1) 2Dq+1 nq

REDQ 2d+1

Table 2: Modular polynomial multiplication complexities.

Example 12. With p = 3, N = 500, n d is much larger than 2N + 1 and thus the classical delayed polynomial multiplication algorithm requires 10 6 multiplications and additions and 10 3 remainderings.

If we choose a double floating point representation and a 4-FQT (i.e 4 coefficients in a word, or a degree 3 substitution), the fully tabulated FQT boils down to 8.6 • 10 4 multiplications and additions and 5.7 • 10 3 divisions. On the one hand, the number of operations is therefore reduced by a factor close to 11. On the other hand the delayed code can compute every coefficient with a single reduction in this case, while the FQT code has to compute less coefficients, but breaks the pipeline.

Example 13. Even by switching to a larger mantissa, say e.g., 128 bits, so that the DQT multiplications are roughly 4 times as costly as double floating point operations, the FQT can still be useful.

Taking p = 1009 and choose d = 2, this still gives around 1.3 • 10 5 multiplications and additions over 128 bits and 7•10 3 divisions. The number of operations is still reduced by a factor of 7. This should therefore still be faster than the delayed multiplication over 32 bits.

On Figure 2, we compare our two implementations with that of NTL [START_REF] Shoup | NTL 5.4.1: A library for doing number theory[END_REF]. We see that the FQT is faster than NTL as long as the same algorithm is used. This shows that our strategy is very useful for small degrees and small primes; not only for the classical algorithm (left) but also for subquadratic ones (right): the use of FQT leads to a gain of an order of magnitude. Note however that in the special case p = 2, NTL offers a very optimized implementation which is still an order of magnitude faster than our general purpose implementation: specific binary routines, such as the ones proposed by [START_REF] Weimerskirch | Generic GF(2) arithmetic in software and its application to ECC[END_REF], enables to pack coefficients as bits of machine words. 6 Application 2: Small Finite Field Extensions

The isomorphism between finite fields of equal sizes gives a canonical representation: any finite field extension is viewed as the set of polynomials modulo a prime p and modulo an irreducible polynomial P of degree k. Clearly we can thus convert any finite field element to its q-adic expansion; perform the FQT between two elements and then reduce the polynomial thus obtained modulo P. Furthermore, it is possible to use floating point routines to perform exact linear algebra as demonstrated by [START_REF] Dumas | Dense linear algebra over wordsize prime fields: the FFLAS and FFPACK packages[END_REF].

We use the strategy of (Dumas et al., 2002, Algorithm 4.1): convert vectors over F p k to q-adic floating point; call a fast numerical linear algebra routine (BLAS); convert the floating point result back to the usual field representation. We improve all the conversion steps as follows:

1. replace the Horner evaluation of the polynomials, to form the q-adic expansion, by a single table look-up, recovering directly the floating point representation;

2. replace the radix conversion and the costly modular reductions of each polynomial coefficient, by a single REDQ operation;

3. replace the polynomial division by two table look-ups and a single field operation. This is presented in Algorithm 6. Line 1 is the table look-up of floating point values associated to elements of the field; line 2 is the numerical computation; line 3 is the first part of the REDQ reduction; lines 4 and 5 are a time-memory trade-off with two table accesses for the corrections of REDQ, combined with a conversion from polynomials to discrete logarithm representation; the last line 6 combines the latter two results, inside the field. A variant of REDQ is used in Algorithm 6, but u i still satisfies u i = 2k-2 j=i µ j q j-i mod p as shown in Theorem 5. Therefore the representations of µ i X j in the field can be precomputed and stored in two tables where the indexing will be made by (u 0 , . . . , u k-1 ) and (u k-1 , . . . , u 2k-2 ) and not by the µ i 's.

Algorithm 6 Fast Dot product over Galois fields via FQT and FQT inverse Input: A field F p k with elements represented as exponents of a generator of the field; Input: two vectors v 1 and v 2 of elements of F p k ; Input: a sufficiently large integer q. Output: R ∈ F p k , with R = v T 1 • v 2 . Tabulated q-adic conversion {Use conversion tables from exponent to floating point evaluation} 1: Set v 1 and v 2 to the floating point vectors of the evaluations at q of the elements of v 1 and v 2 .

The floating point computation

2: Compute r = v 1 T • v 2 ;
Computing a radix decomposition

3: r = REDQ COM P RESSION (r, p, q); Variant of REDQ CORRECTION {µ i is such that µ i = µ i mod p for r = 2k-2 i=0 µ i q i } 4: Set L = representation( k-2 i=0 µ i X i ). 5: Set H = representation(X k-1 × 2k-2 i=k-1 µ i X i-k+1 ).
Reduction in the field

6: Return R = H + L ∈ F p k ;
Thus, this algorithm approaches the performance of the prime field wrapping also for small extension fields. Indeed, suppose the internal representation of the extension field is already by discrete logarithms and uses conversion tables from polynomial to index representations (see e.g., [START_REF] Dumas | Efficient dot product over finite fields[END_REF] for details). Then we choose a time-memory trade-off for the REDQ operation of the same order of magnitude, that is to say p k . The overall memory required by these new tables only doubles and the REDQ requires only 2 accesses. Moreover, in the small extension, the polynomial multiplication must also be reduced by an irreducible polynomial, P. This reduction can be precomputed in the REDQ table look-up and is therefore almost free. Moreover, many things can be factorized if the field representation is by discrete logarithms. For instance, the elements are represented by their discrete logarithm with respect to a generator of the field, instead of polynomials. In this case there are already some table accesses for many arithmetic operations, see e.g., (Dumas, 2004, §2.4) for details.

Theorem 14. Algorithm 6 is correct.

PROOF. We have to prove that it is possible to compute L and H from the u i 's. We have µ 2k-2 = u 2k-2 and µ i = u i -qu i+1 mod p, for i = 0, . . . , 2k -3. Therefore a precomputed table of p k entries, indexed by (u 0 , . . . , u k-1 ), can provide the representation of

L = k-2 i=0 (u i -qu i+1 mod p)X i .
Another table with p k entries, indexed by (u k-1 , . . . , u 2k-2 ), can provide the representation of

H = u 2k-2 X 2k-2 + 2k-3 i=k-1 (u i -qu i+1 mod p)X i . Finally R = X k-1 × 2k-2 i=k-1 µ i X i-k+1 + k-2
i=0 µ i X i needs to be reduced modulo the irreducible polynomial used to build the field. But, if we are given the representations of H and L in the field, R is then equal to their sum inside the field, directly using the internal representations.

Table 3 recalls the respective complexities of the conversion phase in both algorithms. Here, q is a power of two and the REDQ division is computed via the floating point routines of Section 4.

Alg. 1 Alg. 6 Memory 3p k 4p k + 2 1+k⌈log 2 p⌉ Axpy 0 k Div 2k -1 0 Table 0 3 Red ≥ 5k 1
Table 3: Complexity of the back and forth conversion between extension field and floating point numbers

Figure 3 shows the speed of the conversion after the floating point operations. The log scales prove that for q ranging from 2 1 to 2 26 our new implementation is two to three times as fast as the previous one3 .

Furthermore, these improvements allow the extension field routines to reach the speed of 7800 millions of F 9 operations per second4 as shown on Figure 4 5 . The speed-up obtained with these new implementations in also shown on this Figure . It represents a reduction from the 15% overhead of the previous implementation to less than 4% now, when compared over F 11 . We now extend the results of [START_REF] Dumas | Compressed modular matrix multiplication[END_REF] with the REDQ algorithm.

The idea is to use Kronecker substitution to pack several matrix entries into a single machine word. We explore the possibilities of packing on the left or on the right only, together with packing on both matrices of a matrix multiplication.

Middle Product Algorithm

In this section, we show how a dot product of vectors of size d + 1 can be recovered from a polynomial multiplication performed by a single machine word multiplication. This extends to matrix multiplication by compressing both matrices first. We first illustrate the idea for 2 × 2 matrices and

d = 1. The product a b c d × e f g h = ae + bg af + bh ce + dg cf + dh is recovered from Qa + b Qc + d × e + Qg f + Qh = * + (ae + bg)Q + * Q 2 * + (af + bh)Q + * Q 2 * + (ce + dg)Q + * Q 2 * + (cf + dh)Q + * Q 2 ,
where the character * denotes other coefficients. In general, A is an m × k matrix to be multiplied by a k × n matrix B, the matrix A is first compressed into a m × ⌈k/(d + 1)⌉ CompressedRowMatrix, CA, and B is transformed into a ⌈k(d + 1)⌉×n CompressedColumnMatrix, CB. The compressed matrices are then multiplied and the result can be extracted from there. This is depicted on Fig. In terms of number of arithmetic operations, the matrix multiplication CA× CB can save a factor of d + 1 over the multiplication of A × B as shown on the 2 × 2 case above.

The computation has three stages: compression, multiplication and extraction of the result. The compression and extraction are less demanding in terms of asymptotic complexity, but can still be noticeable for moderate sizes. For this reason, compressed matrices are often reused and it is more informative to distinguish the three phases in an analysis. This is done in Section 7.6 (Table 5), where the actual matrix multiplication algorithm is also taken into account.

Partial compression. Note that the last column of CA and the last row of B might not have d + 1 elements if d + 1 does not divide k. Thus one has to artificially append some zeroes to the converted values. On CB this means just do nothing. On CA whose compression is reversed, this means multiplying by Q several times.

Available Mantissa and Upper Bound on Q for the Middle Product

If the product CA × CB is performed with floating point arithmetic we just need that the coefficient of degree d fits in the β bits of the mantissa. Writing CA × CB = c H Q d + c L , we see that this implies that c H , and only c H , must have entries that remain smaller than 2 β . It can then be recovered exactly by multiplication of CA × CB with the correctly precomputed and rounded inverse of Q d as shown e.g., in (Dumas, 2008, Lemma 2). With delayed reduction this means that

d i=0 k d + 1 (i + 1)(p -1) 2 Q d-i < 2 β .
On the other hand, delay reduction requires (cf. Eq (3))

k(p -1) 2 < Q. ( 9
)
Thus the recovery is possible if

Q d+1 < 2 β . ( 10
)
and a single reduction has to be made at the end of the dot product as follows:

Element& init( Element& rem, const double dp) const { double r = dp; // Multiply by the inverse of Q^d with correct rounding r *= _inverseQto_d; // Now we just need the part less than Q=2^t unsigned long rl( static_cast<unsigned long>(r) ); rl &= _QMINUSONE; // And we finally perform a single modular reduction rl %= _modulus; return rem = static_cast<Element>(rl); } Note that one can avoid the multiplication by the inverse of Q when Q is a power of 2, say 2 t : by adding Q 2d+1 to the final result one is guaranteed that the t(d + 1) high bits represent exactly the d + 1 high coefficients. On the one hand, the floating point multiplication is then replaced by an addition. On the other hand, this doubles the size of the dot product and thus reduces by a factor of d+1 √ 2 the largest possible dot product size k.

Middle Product Performance

On Figure 6 we compare our compression algorithm to the numerical double floating point matrix multiplication dgemm of GotoBlas by Goto and van de Geijn (2002) and to the fgemm modular matrix multiplication of the FFLAS-LinBox library by [START_REF] Dumas | Finite field linear algebra subroutines[END_REF]. For the latter we show timings using dgemm and also sgemm over single floating points. This figure shows that the compression (d + 1) is very effective for small primes: the gain over the double floating point routine is quite close to d.

Observe that the curve of fgemm with underlying arithmetic on single floats oscillates and drops sometimes. Indeed, the matrix begins to be too large and modular reductions are now required between the recursive matrix multiplication steps. Then the floating point BLAS routines are used only when the sub-matrices are small enough. One can see the subsequent increase in the number of classical arithmetic steps on the drops around 2048, 4096 and 8192. On Table 4, we show the compression factors modulo 3, with Q a power of 2 to speed up conversions. For a dimension n ≤ 256 the compression is at a factor of five and the time to perform a matrix multiplication is slightly more than a millisecond. Then from dimensions from 257 to 2048 one has a factor of 4 and the times are roughly 16 times the time of the four times smaller matrix. The next stage, from 2048 to 32768 is the one that shows on Figure 5.

Figure 6 shows the dramatic impact of the compression dropping from 4 to 3 between n = 2048 and n = 2049. It would be interesting to compare the multiplication of 3-compressed matrices of size 2049 with a decomposition of the same matrix into matrices of sizes 1024 and 1025, thus enabling 4-compression also for matrices larger than 2048, but with more modular reductions.

Right or Left Compressed Matrix Multiplication

Another way of performing compressed matrix multiplication is to multiply an uncompressed m × k matrix to the right by a row-compressed k × n/(d + 1) matrix. We illustrate the idea on 2 × 2 matrices:

a b c d × e + Qf g + Qh = (ae + bg) + Q(af + bh) (ce + dg) + Q(cf + dh)
The general case is depicted on Fig. 7, center. This is called Right Compressed 9) and ( 10).

The major difference with the Compressed Matrix Multiplication lies in the reductions. Indeed, now one needs to reduce simultaneously the d+1 coefficients of the polynomial in Q in order to get the results. This simultaneous reduction can be made by the REDQ algorithm.

When working over compressed matrices CA and CB, a first step is to uncompress CA, which has to be taken into account when comparing methods. Thus the whole right compressed matrix multiplication is the following algorithm A = Uncompress(CA); CC = A × CB; REDQ(CC)

We see on Figure 8 that the used of REDQ instead of the middle product algorithm has a high benefit. Indeed for small matrices, the conversion can represent 30% of the time and any improvement there has a high impact.

Full Compression

It is also possible to compress simultaneously both dimensions of the matrix product (see Fig. 7,right). This is achieved by using polynomial multiplication with two variables Q and Θ. Again, we start by an example in dimension 2:

a + Qc b + Qd × e + Θf g + Θh = (ae + bg) + Q(ce + dg) + Θ(af + bh) + QΘ(cf + dh) .
More generally, let d q be the degree in Q and d θ be the degree in Θ. Then, the dot product is:

a • b = [ dq i=0 a i0 Q i , . . . , dq i=0 a in Q i ] × [ d θ j=0 b 0j Θ j , . . . , d θ j=0 b nj Θ j ], = k l=0 ( dq i=0 a il )( d θ j=0 b lj )Q i Θ j = dq i=0 d θ j=0 ( k l=0 a il b lj )Q i Θ j .
In order to guarantee that all the coefficients can be recovered independently, Q must still satisfy Eq. ( 9) but then Θ must satisfy an additional constraint:

Q dq+1 ≤ Θ. ( 12
)
This imposes restrictions on d q and d θ :

Q (dq+1)(d θ +1) < 2 β . (13)

CMM Comparisons

In Table 5, we summarize the differences between the algorithms presented on Figures 5 and7. As usual, the exponent ω denotes the exponent in the complexity of matrix multiplication. Thus, ω = 3 for the classical matrix multiplication, while ω < 3 for faster matrix multiplications, as used in (Dumas et al., 2009, §3.2). For products of rectangular matrices, we use the classical technique of first decomposing the matrices into square blocks and then using fast matrix multiplication on those blocks.

Compression Factor

The costs in Table 5 are expressed in terms of a compression factor e, that we define as

e := β log 2 (Q) ,
where, as above, β is the size of the mantissa and Q is the integer chosen according to Eqs. ( 9) and ( 10), except for Full Compression where the more constrained Eq. ( 13) is used. Thus the degree of compression for the first three algorithms is just d = e-1, while it becomes only d = √ e -1 for the full compression algorithm (with equal degrees d q = d θ = d for both variables Q and Θ). In terms of asymptotic complexity, the cost in number of arithmetic operations is dominated by that of the product (column Operations in the table), while reductions and conversions are linear in the dimensions. This is well reflected in practice. For example, with algorithm Right Compression on matrices of sizes 10, 000×10, 000 it took 90.73 seconds to perform the matrix multiplication modulo 3 and 1.63 seconds to convert the resulting matrix. This is less than 2%. For 250 × 250 matrices it takes less than 0.00216 seconds to perform the multiplication and roughly 0.0016 seconds for the conversions. There, the conversions account for 43% of the time and it therefore of extremely high importance to optimize the conversions.

Algorithm

In the case of rectangular matrices, the second column of Table 5 shows that one should choose the algorithm depending on the largest dimension: CMM if the common dimension k is the largest, Right Compression if n if the largest and Left Compression if m dominates. The gain in terms of arithmetic operations is e ω-2 for the first three variants and e ω-1 2 for full compression. This is not only of theoretical interest but also of practical value, since the compressed matrices are then less rectangular. This enables more locality for the matrix computations and usually results in better performance. Thus, even if ω = 3, i.e., classical multiplication is used, these considerations point to a source of speed improvement.

The full compression algorithm seems to be the best candidate for locality and use of fast matrix multiplication; however the compression factor is an integer, depending on the flooring of either β log 2 (Q) or β log 2 (Q) . Thus there are matrix dimensions for which the compression factor of e.g., the right compression will be larger than the square of the compression factor of the full compression. There the right compression will have some advantage over the full compression.

If the matrices are square (m = n = k) or if ω = 3, the products all become the same, with similar constants implied in the O(), so that apart from locality considerations, the difference between them lies in the time spent in reductions and conversions. Since the REDQ e reduction is faster than e classical reductions [START_REF] Dumas | Q-adic transform revisited[END_REF], and since INIT e and EXTRACT e are roughly the same operations, the best algorithm would then be one of the Left, Right or Full compression. Further work would include implementing the Full compression and comparing the actual timings of conversion overhead with that of the Right algorithm and that of CMM.

Conclusion

We have proposed a new algorithm for simultaneous reduction of several residues stored in a single machine word. For this algorithm we also give a time-memory trade-off implementation enabling very fast running time if enough memory is available.

We have shown very effective applications of this trick for packing residues in large applications. This proves efficient for modular polynomial multiplication, extension fields conversion to floating point and linear algebra routines over small prime fields.

Further work is needed to compare of running times between different choices for q. Indeed our experiments were made with q a power of two and large table look-up. With q a multiple of p the table look-up is not needed but divisions by q i will be more expensive. A possibility would be taking q in the form q = p2 t , then only divisions by p or p i would be made.

It would also be interesting to see in practice how this trick extends to larger precision implementations: on the one hand the basic arithmetic slows down, but on the other hand the trick enables a more compact packing of elements (e.g., if an odd number of field elements can be stored inside two machine words, etc.).

Lemma 6 .

 6 For r ∈ N and a, b ∈ N * , Let k = ⌊r/ab⌋, so that kab ≤ r < (k + 1)ab. Then kb ≤ r/a < (k + 1)b and since kb is an integer it follows that kb ≤ ⌊r/a⌋ < (k + 1)b. Dividing by b yields ⌊⌊r/a⌋/b⌋ = k. The other side of the identity follows by exchanging a and b.

Figure 1 :

 1 Figure 1: REDQ 3 COMPRESSION with 2 axpy

Figure 2 :

 2 Figure 2: Number of classical and Karatsuba polynomial multiplications modulo 3 per second on a Xeon 3.6 GHz (logarithmic scale)

Figure 3 :

 3 Figure 3: Small extension field conversion speed on a Xeon 3.6GHz

  5

Figure 4 :Figure 5 :

 45 Figure 4: Speed of finite field Winograd matrix multiplication on a XEON, 3.6 GHz

Figure 6 :

 6 Figure 6: Compressed matrix multiplication compared with dgemm (the floating point double precision matrix multiplication of GotoBlas) and fgemm (the exact routine of FFLAS) with double or single precision.

Figure 7 :

 7 Figure 7: Left, Right and Full Compressions Matrix Multiplication. Left Compressed Matrix Multiplication is obtained by transposition. Here also Q and d must satisfy Eqs. (9) and (10).The major difference with the Compressed Matrix Multiplication lies in the reductions. Indeed, now one needs to reduce simultaneously the d+1 coefficients of the polynomial in Q in order to get the results. This simultaneous reduction can be made by the REDQ algorithm.When working over compressed matrices CA and CB, a first step is to uncompress CA, which has to be taken into account when comparing methods. Thus the whole right compressed matrix multiplication is the following algorithm A = Uncompress(CA); CC = A × CB; REDQ(CC) (11)

  matrix multiplication with Goto BLAS on a XEON, 3.6 GHz CMM Compressed double fgemm mod 3 Right Compressed double fgemm mod 3

Figure 8 :

 8 Figure 8: Right Compression and CMM.

Table 4 :

 4 Compression factors for different common matrix dimensions modulo 3, with 53 bits of mantissa and Q a power of 2.

	Compression	2 3..4 5..8	8	7	6	5	4	3
	Degree d	1	5	9	7	6	5	4	3	2
	Q-adic	2 3	2 4	2 5	2 6	2 7	2 8	2 10	2 13	2 17
	Dimensions	2 ≤ 4 ≤ 8 ≤ 16 ≤ 32 ≤ 64 ≤ 256 ≤ 2048 ≤ 32768

Table 5 :

 5 Number of arithmetic operations for the different algorithms 7.6.2 Analysis

	Operations	Reductions	Conversions
	CMM Right Comp. O mk n O mn k e e Left Comp. O nk m e Full Comp. O k mn e	ω-2 ω-2 ω-2 ω-1 2	m × n REDC e REDQ e m × n m e × n REDQ e m √ e × n √ e REDQ e	1 e mn INIT e e mn EXTRACT e 1 1 e mn EXTRACT e 1 e mn INIT e

On a Xeon 3.6GHz using doubles, addition, multiplication and axpy take roughly the same time, while division is 10 times slower and fmod (floating point remainder) again

2.5 times as long as division

How ties are broken is irrelevant here.

On a 32 bit Xeon.

On a XEON, 3.6 GHz, using Goto BLAS-1.09 dgemm as the numerical routine[START_REF] Goto | On reducing TLB misses in matrix multiplication[END_REF] and FFLAS fgemm for the fast prime field matrix multiplication[START_REF] Dumas | Dense linear algebra over wordsize prime fields: the FFLAS and FFPACK packages[END_REF].

The FFLAS routines are available within the LinBox 1.1.4 library (LinBox Group, 2007) and the FQT is in implemented in the givgfqext.h file of the Givaro 3.2.9 library[START_REF] Dumas | Givaro 3.2.9: C++ library for arithmetic and algebraic computations[END_REF].

* This research was partly supported by the French National Research Agency (ANR Safescale, ANR Gecko)