H Lapierre

C Basile
email: cbasile@ujf-grenoble.fr

T Berly

E Canard

Potassic late orogenic Stephanian volcanism in the South West French Massif central (Decazeville, Figeac, Lacapelle-Marival basins): an example for mantle metasomatism along strike-slip faults?

Keywords: late orogenic magmatism, shoshonitic magma, metasomatism, Sillon Houiller fault, Variscan magmatisme tardi-orogénique, shoshonitique, métasomatisme, Sillon Houiller, Varisque

In the Southwestern part of the French Massif central (Decazeville basin, at the Sillon Houiller fault termination; Figeac and Lacapelle-Marival basins along the Argentat fault),

Stephanian volcanism exhibits shoshonitic affinities. Their chondrite-normalized Rare Earth Element (REE) patterns are enriched in light REE, but almost flat for heavy REE, with marked negative Eu anomalies. Primitive mantle-normalized element spectra show negative Nb, Ta, P, Sm, Ti, and positive Th, U, Pb anomalies, respectively. ε Nd values are negative and homogeneous (-6 to -4). This volcanism shares the same geochemical patterns as the late-orogenic Stephanian-Permian magmatism from the southern part of the Variscan belt (Pyrénées, Alps, Sardinia). We explain these geochemical characteristics as resulting from the partial melting of a metasomatised mantle. We propose a new mechanism to explain this melting process: horizontal displacement along the main late-orogenic strike-slip faults might bring into contact a hydrated lower crust with the lithospheric mantle. Mantle metasomatism within the strike-slip fault zone may then induce partial melting.

Résumé

Le volcanisme Stéphanien du Sud Ouest du Massif central français, observé dans les bassins de Decazeville à l'extrémité du Sillon Houiller, de Figeac et Lacapelle-Marival le long de la faille d'Argentat, présente des caractéristiques géochimiques similaires. Il s'agit d'un magmatisme d'affinité shoshonitique, présentant des spectres de terres rares normalisés aux chondrites enrichis en terre rares légères, mais presque plat pour les terres rares lourdes, avec une anomalie négative en Eu marquée. Les spectres élémentaires normalisés au manteau primitif montrent des anomalies négatives en Nb, Ta, P, Sm et Ti, et des anomalies positives en Th, U et Pb. Les valeurs d'ε Nd sont négatives et également très homogènes (-6 à -4). Ce volcanisme présente les mêmes caractéristiques que le volcanisme tardi-orogénique stéphano-permien du Sud de la chaîne varisque (Pyrénées, Alpes, Sardaigne). Ces caractéristiques sont interprétées comme résultant de la fusion partielle d'une source mantellique métasomatisée. Il est proposé un nouveau mécanisme provoquant cette fusion, par l'intermédiaire du déplacement horizontal des failles décrochantes tardi-orogéniques majeures, qui pourrait mettre en contact une croûte inférieure hydratée et le manteau lithosphérique. La métasomatose du manteau dans la zone de faille décrochante provoquerait sa fusion.

Introduction

Late orogenic times are key periods in mountain belt evolution, when the continental lithosphere acquires a new structure that will control further evolution, including its reworking in new orogens. The continental crust is particularly modified by late orogenic magmatism and tectonism. Bimodal magmatism is widespread, and results primarily from partial melting of the crust and underlying mantle [START_REF] Bonin | Do coeval mafic and felsic magmas in post-collisional to within-plate regimes necessarily imply two constrasting, mantle and crustal, sources? A review[END_REF]. The magmas evolving from mantle melts are commonly potassium-rich [High Potassium Calc Alcaline: HKCA, e.g. Liégeois et al., 1998], but the mechanisms associated with this enrichment are still a matter of discussion [START_REF] Gibson | Ultrapotassic magmas along the flank of the Oligomiocene Rio-Grande rift, USA -Monitors of the zone of lithospheric mantle extension and thinning beneath a continental rift[END_REF]Liégeois et al., 1998;[START_REF] Miller | Postcollisional potassic and ultrapotassic magmatism in SW Tibet: geochemical and Sr-Nd-Pb-O isotopic constraints for mantle source characteristics and petrogenesis[END_REF]Pecerillo and Martinotti, 2006]. This magmatism is coeval with a change in the tectonic regime, from regional compression to strike-slip, then extensional regime, that precedes the tectonic quiescence [e.g. Blès et al., 1989]. The strike-slip regime generates numerous vertical faults that cross-cut the orogenic structures.

At the end of the Variscan orogeny, Westphalian synkinematic plutonism accommodated extension parallel to the mountain belt [START_REF] Faure | Late orogenic carboniferous extensions in the Variscan French Massif central[END_REF]. Before the widespread post-orogenic Permian extension, the Stephanian period represents a transition with numerous strike-slip faults and associated intracontinental coal-bearing basins. In the French Massif central, the Sillon Houiller is one of these strike-slip faults (Fig. 1), coeval to the formation of numerous deep and narrow Stephanian basins first infilled with potassic volcanism [START_REF] Letourneur | Le grand Sillon houiller du Plateau Central français[END_REF].

Potassic volcanism is also known in other Stephanian basins, as for example along the Argentat fault in the western part of the Massif central (Fig. 1). These two faults meet together in the southwestern Massif central, where a comparison of these late orogenic potassic volcanism in coeval strike-slip and extensional basins can be done.

Geological setting

The present study focuses on the volcanism located in the Decazeville basin, at the southern termination of the Sillon Houiller, and in Figeac and Lacapelle-Marival basins along the Argentat fault.

The Sillon Houiller fault is a 270 km-long, nearly N-S to NNE-SSW striking left-lateral fault, which crosses the entire Massif central (Fig. 1). The cumulated horizontal displacement has been estimated as ranging from 70 to 100 kilometres [START_REF] Grolier | L'évolution tectonique du grand Sillon houiller du Massif central français[END_REF]. This fault has been considered as a transfer fault accommodating differences between NW-SE and NE-SW extension directions on its western and eastern sides [START_REF] Burg | Le Sillon houiller du Massif central français: faille de transfert pendant l'amincissement crustal de la chaîne varisque?[END_REF]. However, this interpretation is not supported by the age of extension, which is Namurian-Westphalian for the western NW-SE trend and Stephanian-Autunian for the NE-SW extension [START_REF] Faure | Late orogenic carboniferous extensions in the Variscan French Massif central[END_REF]. Nevertheless, due to its NNE-SSW trend and connection with coeval Stephanian basins, the Sillon Houiller fault can be considered either as a transfer [START_REF] Faure | Late orogenic carboniferous extensions in the Variscan French Massif central[END_REF] or transform fault [START_REF] Basile | A new interpretation of Stephanian deformation in the Decazeville basin (Massif central, France): consequences on late Variscan tectonism[END_REF] in Stephanian times. The formation of the Decazeville basin was tectonically controlled by its horsetail splay termination in a transtensional regime [START_REF] Basile | A new interpretation of Stephanian deformation in the Decazeville basin (Massif central, France): consequences on late Variscan tectonism[END_REF]. In the lowermost part of the basin infilling, a thick volcanic complex overlies coarse conglomerates (Fig. 1). The upper part of the 1 800 m-thick sedimentary cross-section is almost devoid of volcanism, expressed only by a few thin layers of tuffs intercalated in lacustrine to deltaic continental sediments.

The Argentat fault is a 160 km-long curved fault in the western part of the Massif central, running from the Bosmoreau Carboniferous basin to the north and the Decazeville basin to the south (Fig. 1). This fault strikes N-S in its northern part and NW-SE in its southern part. The Argentat fault first acted as a right-lateral ductile normal fault, then as a brittle strike-slip fault during the Late Carboniferous [START_REF] Feix | Arguments en faveur d'un épisode majeur en cisaillement dextre le long de la faille d'Argentat (Massif central)[END_REF][START_REF] Roig | -Modélisation 3D du paléochamp hydrothermal As + Au (330-300 Ma) le long de la faille d'Argentat (Massif central français)[END_REF]. Figeac, St Perdoux and Lacapelle-Marival represent the southernmost basins associated with this fault (Fig. 1). They are not directly bound by the Argentat fault, but appear as the Argentat basin located on a basement flexure [START_REF] Roig | Le bassin houiller d'Argentat (Massif central français), conséquence topographique d'un plissement de son substratum varisque[END_REF]. Due to of Mesozoic sedimentary cover, the connection between Figeac and Lacapelle-Marival basins is presently hidden (Fig. 1), but both basins probably belong to a single one whose remnants were drilled below the Mesozoic sediments [START_REF] Lefavrais- | Notice explicative[END_REF]. As for the Decazeville basin, the 300 to 400 m-thick volcanic complex is located at the base of the basin infilling, overlain by coarse detrital sediments. Numerous NW-SE-trending feeder dikes also cross-cut the basement close to the Lacapelle-Marival basin (Fig. 1), providing evidence for a NE-SW extension during volcanism and basin subsidence.

In these basins, sediments are middle Stephanian (Stephanian B) in age, with early

Stephanian ages (A) suspected in the St Perdoux basin, and the base of late Stephanian (C) observed in the Decazeville basin (Fig. 1) [START_REF] Vetter P | Géologie et paléontologie des bassins houillers de Decazeville, de Figeac et du détroit de Rodez[END_REF]. Permian sediments unconformably overlay Stephanian beds in the eastern part of the Decazeville basin. They were also found in the Figeac basin [START_REF] Lefavrais- | Notice explicative[END_REF]. The age of volcanism is also assummed to be Stephanian, but this age is controversial in the Decazeville basin. From palynology, the older sediments were dated as the base of middle Stephanian [START_REF] Vetter P | Géologie et paléontologie des bassins houillers de Decazeville, de Figeac et du détroit de Rodez[END_REF], whereas U-Pb dating on unlocated volcanic zircons provided an older Visean age [333 Ma, Bruguier et al., 1998].

Volcanism is expressed as basalts, intermediate or felsic rocks, mainly as lava flows, but also as ejecta and numerous dikes in the crystalline basement. [START_REF] Morre | Etude pétrographique sur les roches volcaniques carbonifères de la région de Figeac[END_REF] and [START_REF] Morre | -Etude pétrographique sur les roches volcaniques carbonifères de Lacapelle-Marival[END_REF] emphasized the potassium-rich character of these lava, especially in intermediate rocks in which K 2 O ranges from 6 to 11%.

-In the northern part of Decazeville basin, samples were collected from two areas (Fig. 1), the Bouran section along the D963 road (samples CLB95-1 to -8), and the Lot river banks on either side of Livinhac bridge (CLB95-9 to -17, and CLB96-1). The Bourran section consists of volcaniclastic rocks and volcanic flows conformably overlain by detrital sediments [Campagnac unit, Vetter, 1968]. The graded bedded volcaniclastic rocks rest unconformably on Variscan micaschists, and consist of breccias and lapilli-tuffs, interlayered with massive volcanic flows. The breccias are composed of meter-sized volcanic blocks caught in a tuffaceous matrix while the lapilli-tuffs include centimeter-sized angular volcanic fragments.

The lapilli-tuffs locally grade to millimeter-sized crystal tuffs. Near Livinhac, massive basalts are exposed on both sides of the Lot river. The northern outcrop is cross-cutted by a porphyric dike (CLB95-17).

-In the Figeac area, two sections were sampled, East (Célé) and North (Planioles) of the city, respectively (Fig. 1). Along the Planioles section, hyalopilitic basalts crop out (CLB96-3). They are successively overlain by vitric tuffs, very altered rhyolites (CLB96-2) and detrital sediments. In the Célé valley, two outcrops of massive flows of intermediate composition (CLB96-5) can be found.

-Near Lacapelle-Marival, two volcanic levels are interbedded with Carboniferous shales, sandstones and conglomerates. Sampling concerns massive flows of mafic to intermediate composition in the northern part of the basin (CLB96-13 to -17) and small outcrops of intermediate to acidic lavas to the south (CLB96-6 to -12, Fig. 1).

Petrography and mineralogy

The studied volcanic rocks are frequently severely weathered. Plagioclase is often replaced by albite or sericite, or even quartz and epidote. Olivine is transformed to iddingsite, serpentine or smectites. Orthopyroxene and clinopyroxene remain more commonly fresh.

Biotite and hornblende are often chloritized. When altered, Ti-rich oxides are filled with titanite. Mafic glass has crystallized into chlorite ± smectites while vesicles are filled with calcite or smectites. The groundmass of intermediate lavas and rhyolites consists of fine grained quartz associated with alkali feldspar (albite or adularia). As no metamorphism has been observed in the sediments overlying the lavas, the occurence of chlorite probably indicates hydrothermal alteration coeval with lava eruption.

In both Decazeville and Lacapelle-Marival basins, mafic rocks are commonly glomeroporphyritic with clots of large zoned plagioclase, orthopyroxene and clinopyroxene. These phenocrysts are associated with either olivine or amphibole and biotite ± quartz. The groundmass includes plagioclase microliths and oxides. Quartz xenocrysts rimmed by orthopyroxene may occur. In Lacapelle-Marival area, some basalts differ by their trachytoid texture that consists of plagioclase laths (30 %) and less abundant larger phenocrysts (~ 10 %). The intermediate lavas exhibit seriate textures with millimeter-sized alkali feldspar phenocrysts. Their groundmass includes sandstone inclusions and quartz-filled flattenned vesicles. The rhyolites display quartz and feldspar phenocrysts (20 %) and biotite needles (~ 10 %), that includes apatite and zircon. The Figeac rhyolites differ from those from Lacapelle-Marival by the presence of green hornblende.

When preserved (CLB95-5, CLB95-16, CLB96-1, CLB96-3), plagioclase composition of mafic volcanic rocks varies from Ca-rich cores (An 98-60 ) to Na-rich rims (An 50-40 ). Preserved alkali feldspar (CLB96-11) gives sanidine compositions (Ab 57-66 ). Clinopyroxene phenocrysts (CLB96-1, CLB95-16, CLB95-5) show homogeneous augite compositions [Wo 39-44 En 44-52 Fs 12-21 ; [START_REF] Morimoto | -Nomenclature of pyroxenes[END_REF], while preserved orthopyroxene (CLB95-5, CLB95-16) has bronzite compositions (Wo 2-2.5 En 77-78 Fs 19-21 ).

Geochemistry

Analytical techniques

Thirty one samples have been analysed for whole rock major-and trace-element concentrations. Among them, twenty one have been analysed for Nd, and nine for Sr isotopic chemistry (tabl. I). Major and compatible trace-elements have been measured by ICP-AES at the Centre de Recherches Pétrographiques et Géochimiques de Nancy (samples CLB95) or by X-ray fluorescence (XRF) at the Earth Sciences Department of the Australian National University at Canberra (Australia) (samples CLB96). The trace-element concentrations have been determined by inductively coupled plasma-mass spectrometry (ICP-MS, VG-PQ 2+ ) at the Laboratoire de Géodynamique des Chaînes Alpines de l'Université Joseph Fourier (Grenoble), after acid dissolution, using the procedures described by [START_REF] Barrat | Determination of rare earth elements in sixteen silicate reference samples by ICP-MS after Tm addition and ion exchange separation. -Geost[END_REF].

Standards used for the analyses were JB2, WSE, BIR-1, JR1, and BHVO.

Sr and Nd isotopic ratios were measured at the Laboratoire de Géochimie isotopique de l'Université Paul Sabatier (Toulouse), on a Finnigan MAT261 multicollector mass spectrometer using the analytical procedures of [START_REF] Lapierre | -Is the lower Duarte Igneous Complex (Hispaniola) a remnant of the Caribbean plume-generated oceanic plateau?[END_REF]. Results on La Jolla Nd standards yielded 143 Nd/ 144 Nd = 0.511850 ± 8 (mean on 39 runs). Results on NBS 987 Sr standard yielded 87 Sr/ 86 Sr = 0.710250 ± 5 (mean on 200 runs). The isotopic data have been corrected for "in situ decay" for an age of 300 Ma.

Alteration and elemental mobility

Most of the analysed rocks have a loss on ignition lower than 6% (LOI, tabl. I), due to the presence of calcite-or smectites-filled vesicules. Moreover, the altered nature of the studied Carboniferous volcanic rocks means that before any petrological inferences can be drawn from their chemistry, possible effects of subsolidus transformations must be considered.

As Zirconium (Zr) is considered to be immobile during low-grade alteration of igneous rocks of mafic to felsic composition by hydrothermal fluids [START_REF] Gibson | The trace-element composition of the lavas and dikes from a 3-km vertical section through the lava pile of eastern Iceland[END_REF], it has been used ot test the mobility of minor (TiO 2 ), and incompatible trace elements (Rb, Sr, Pb, La, Nb, Y, U and Th) in figure 2. Rb, Sr, Pb and La do not display any correlation with Zr. This implies that Large Ion Lithophile Elements (LILE) have been extensively mobilised. For example, some samples are enriched in Pb (others in Ba), probably in relation with hydrothermal circulations [START_REF] Marignac | Ore deposits of the French Massif central: insight into the metallogenesis of the Variscan collision belt[END_REF]. Thus, variations in LILE are not expected to reflect the composition of the sources (Sr isotopic ratios) or the magmatic processes. These variations will not be discussed further on. Similarly, the behaviour of K and Ca also exclude their use for sample classification.

In contrast, other elements such as U, Th, Nb or TiO 2 correlate rather well with Zr, suggesting relative immobility. However, some samples (CLB95-17, 95-15, 95-1, 96-5) differ from the correlation trend on the Nb-Zr plot. The same samples, together with CLB95-6, 95-12, 96-2, 96-8 show a negative correlation on the TiO 2 -Zr diagram, while all other samples exhibit a positive correlation. CLB95-17 is the most mafic sample (SiO 2 51.9%, MgO 7,7%, tabl. I); CLB95-1, 95-6, 95-12, 95-15, 96-2, 96-5, 96-8 represent the most SiO 2 -rich and MgO-poor samples (SiO 2 > 65% and MgO < 1.3%, tabl. I).

In the Y-Zr diagram, each sampling site corresponds to distinct Y-Zr correlation trends. This suggests that the Carboniferous volcanic rocks are not related by simple fractionation processes, but rather reflect variations of partial melting ratios, and/or complex fractionation processes dominated by clinopyroxene or amphibole accumulation and/or removal.

Volcanic rocks terminology

Because of the alteration experienced by the Carboniferous volcanic rocks, it was difficult to name the different facies on the basis of mineralogy only, especially for the mafic rocks.

The various types of volcanic rocks have been discriminated using Ti and immobile trace elements Zr, Nb and Y (figure 3). The most mafic rock (CLB95-17) plots in the andesite/basalt field. The samples that exhibit a positive correlation in the TiO 2 -Zr diagram (figure 2) have a constant Zr/Ti ratio and plot near the fence separating the field or rhyodacites from that of andesites, with the exception of sample 95-13 which plots in the alkali basalt field. Among the most SiO2-rich samples that show a negative correlation in the TiO 2 -Zr diagram (figure 2), CLB95-1 plots in the trachy-andesite field, CLB95-6, 95-12, 96-2 and 96-8 in the rhyodacite field, CLB96-5 in the rhyolite field and CLB95-15 at the fence between trachyte and comendite.

Because K is highly mobile during alteration and low grade metamorphism processes, we cannot rely on the K contents of the Carboniferous volcanic rocks (half of the samples have K contents from 6 to 12%) to infer calc-alkaline or shoshonitic affinities. [START_REF] Müller | Geochemical discrimination between shoshonitic and potassic volcanic rocks in different tectonic settings: a pilot study[END_REF] proposed a Th/Yb versus Ta/Yb correlation diagram to discriminate arc-tholeiites, calcalkaline and shoshonitic suites. As shown in figure 4, most of the Carboniferous volcanic rocks fall in the shoshonitic field, with the exception of 95-17 andesite/basalt which fall in the calc-alkaline field. In this plot the most SiO2-rich samples exhibit a more shoshonitic affinity.

Trace-element geochemistry

Figure 5 illustrates the chondrite-normalized [START_REF] Sun | Chemical and isotopic systematics of oceanic basalts: implications for mantle composition and processes[END_REF]] Rare Earth Elements (REE) patterns of the Carboniferous volcanic rocks. All the rocks are enriched in Light (L) REE relative to the heavy (H) REE [2.8 ≤ (La/Yb) n ≤ 31.9]. However, HREE patterns are either flat or gently sloping, suggesting that LREE enrichment reflects a metasomatised source [START_REF] Keppler | Constraints from partitioning experiments on the composition of subduction-zone fluids[END_REF][START_REF] Tatsumi | Trace element transport during dehydration processes in the subducted oceanic crust: 2. Origin of chemical and physical characteristics in arc magmatism[END_REF][START_REF] Beccaluva | Coexisting anorogenic and subduction-related metasomatism in mantle xenoliths from the Betic Cordillera (southern Spain)[END_REF]. All samples also show negative Eu anomalies (0.38 < Eu/Eu* < 0.84). These anomalies suggest an early fractionation of plagioclase during magma differentiation. The most SiO2-rich rocks (e.g. samples CLB95-1 and -6 in Bourran) are commonly less enriched than mafic rocks (e.g. samples CLB95-4 and -8 in Bourran), especially for HREE (Figure 5). This suggests again that these lavas cannot result from fractional crystallisation of basalts.

There are no significative differences between the sampling sites, with the exception of negative Ce anomalies (0.67 ≤ Ce/Ce* ≤ 0.38) for four rhyodacites from the southern part of Lacapelle-Marival basin (CLB96-6, CLB96-10 to -12). Two rhyodacites from Decazeville Figure 6 shows the primitive mantle-normalized spiderdiagrams [START_REF] Sun | Chemical and isotopic systematics of oceanic basalts: implications for mantle composition and processes[END_REF]]. Whatever their location or mafic character, Carboniferous volcanic rocks exhibit similar trace element patterns, characterized by negative Nb, Ta, P, Sm, Ti and positive Th, U and Pb anomalies, respectively. As shown previously (Fig. 2), Sr vary from 10 to 330 ppm.

The most mafic rocks exhibit a homogeneous composition for all trace elements, with the exception of CLB95-17 that differs from the other mafic rocks by a Zr and Hf depletion. CLB95-17 as other samples (CLB95-1, 95-13, 95-15, 96-5) that define a negative trend in TiO2/Zr correlation (Fig. 2) are all characterized by a negative Y anomaly.

Nd and Sr isotopic geochemistry

The Carboniferous volcanic rocks span a rather small range of ε Nd values, from -5.8 to -3.9, with the relative exception of CLB96-5 rhyolite (ε Nd = -2.2). All ε Nd are negative, and do not present significative variations with rock type or rock chemistry (tabl. I). These very homogeneous values suggest that these lavas are cogenetic, and indicate a crustal signature in the melting source. The 87 Sr/ 86 Sr ratio is less homogeneous and ranges between 0.71128 and 0.70415; ε Sr ranges from 0 to 325 (tabl. I). As exposed above, Rb and Sr mobility during hydrothermalism and alteration precludes any use of these elements as markers of the magmatic sources.

Discussion

Geochemistry

Whatever their location or petrology, all the analysed samples display high (La/Yb) n ratios but almost flat HREE patterns, marked negative Eu anomalies, homogeneous negative ε Nd values (between -5.8 and -3.9) and in primitive mantle-normalized plots, negative Nb, Ta, P, The occurence of mafic rocks implies partial melting of the mantle. As the most SiO2-rich rocks display REE spectrums that are less enriched in LREE than mafic rocks, they can not be linked together only by crystal fractionation processes. However, as they present similar characteristics, all these lavas may have been formed from similar sources and similar mechanisms. All these rocks belong to a shoshonitic suite, characterized by the absence of Ti and Fe enrichment during crystal fractionation. This is a classical feature of hydrous melts [START_REF] Gill | Orogenic andesites and plate tectonics[END_REF]. Similarly, as Nb and Ta are insoluble in water, their negative anomalies in primitive mantle-normalized plots can also result from a mantle source metasomatised by water.

As shown by the occurence of sedimentary xenoliths in these lavas, the positive U, Th and Pb anomalies may be explained by contamination by continental crust, while the negative Ce anomaly may indicate contamination by oceanic rocks (sediments or altered oceanic crust) [START_REF] Bohrson | Genesis of silicic peralkaline rocks in an ocean island setting by crustal melting and open-system processes: Socorro Island, Mexico[END_REF]. These contaminations may occur during magma ascent, but more probably reflect the direct contribution of hydrous crustal materials to mantle melting, as suggested by the negative ε Nd values, the homogeneous shape of REE patterns, and the Ti/Y versus Zr/Y and Nb/La versus Th/La diagrams (Fig. 7). All samples display a positive linear trend, except one basalt (CLB95-17) in both diagrams and the most SiO 2 -rich rocks in Ti/Y versus Zr/Y diagram. This suggests that the contamination by hydrous crustal material is strictly linked to melting of the source, i.e. that the mantle protolith was metasomatised by hydrous crustal materials, that participated in its melting. Among the analysed rocks, the Decazeville acidic lavas have the highest Ti/Y, Zr/Y, Nb/La and Th/La ratios, suggesting that the proportion of crustal material was higher in the source of these volcanic rocks than for

Structural control on magmatism

Despite differences in structural setting, there are no marked differences in the petrology or geochemistry of rocks from Decazeville, Figeac and Lacapelle-Marival basins. Thus, the faults and the tectonic regime control the formation of these Stephanian basins, where the volcanism took place, but they do not influence its petrology or geochemistry. The NW-SE striking dikes found in the basement of Lacapelle-Marival basin, and the occurence of volcanism as the first infilling of all the basins along the Sillon Houiller fault indicate a strikeslip tectonic regime at the French Massif central scale, with a NE-SW stretching direction and a NW-SE shortening direction. This regime is associated to the left-lateral slip on the Sillon Houiller fault, and to extension or right-lateral transtension along the Argentat fault. This seems to be the main tectonic regime at the time of the infilling of the Stephanian basins along these faults [START_REF] Basile | A new interpretation of Stephanian deformation in the Decazeville basin (Massif central, France): consequences on late Variscan tectonism[END_REF]. As the same volcanism systematically occurred in all the basins along the Sillon Houiller fault [START_REF] Letourneur | Le grand Sillon houiller du Plateau Central français[END_REF], it can be proposed that the Lacapelle-Marival volcanism, located only 25 km West of the Sillon Houiller fault, came from the same source as the Decazeville volcanism. The local variations in the tectonic regime only controlled the volume and surface covered by volcanic lava flows.

Late-orogenic Variscan magmatism

Stephanian volcanism from the French Massif central belongs to the extended K-rich igneous province emplaced during Stephanian to Early Permian times in the southern part of the European Variscan belt [Pyrénées: [START_REF] Innocent | Sr-Nd isotope and trace-element geochemistry of late Variscan volcanism in the Pyrénées: magmatism in post-orogenic extension?[END_REF]Western Alps: Cannic et al., 2002;Central Alps: Schaltegger et al., 1991;Southern Alps: Rottura et al., 1998[START_REF] Cortesogno | The Variscan post-collisional volcanism in Late Carboniferous-Permian sequences of Ligurian Alps, Southern Alps and Sardinia (Italy): a synthesis[END_REF]Corsica and Sardinia: Tommasini et al., 1995]. In the Pyrénées [START_REF] Bixel | -Magmatisme, tectonique et sédimentation dans les fossés stéphano-permiens des Pyrénées occidentales[END_REF] and Western Alps [START_REF] Capuzzo | High-precision dating and origin of synsedimentary volcanism in the Late Carboniferous Salvan-Dorenaz basin (Aiguilles Rouges Massif, Western Alps). -Schweiz[END_REF], this volcanism is located in basins clearly controlled by strike-slip faults. Together with the K-rich character, these igneous provinces share similar geochemical patterns with the Stephanian volcanism for the SW Massif central: in chondrite-normalized diagrams, LREE-enriched patterns, flat HREE spectra and negative Eu anomalies; in primitive mantle normalized diagrams, negative Nb, Ta and Ti anomalies, and positive Th and U anomalies; and negative ε Nd values (-3 to -7), with the exception of the most recent (Permian) lava that have slightly positive ε Nd values [START_REF] Innocent | Sr-Nd isotope and trace-element geochemistry of late Variscan volcanism in the Pyrénées: magmatism in post-orogenic extension?[END_REF][START_REF] Rottura | Origin and significance of the Permian high-K calc-alkaline magmatism in the central-eastern Southern Alps, Italy[END_REF][START_REF] Cannic | Late orogenic evolution of the Variscan lithosphere: Nd isotopic constraints from the western Alps. -Schweiz[END_REF]]. It is noteworthy that these geochemical characteristics were homegeneous during a large time span and in various parts of the Variscan belt.

The origin of potassic late-orogenic magmatism is still a matter of debate [START_REF] Gibson | Ultrapotassic magmas along the flank of the Oligomiocene Rio-Grande rift, USA -Monitors of the zone of lithospheric mantle extension and thinning beneath a continental rift[END_REF]Liégeois et al., 1998;[START_REF] Miller | Postcollisional potassic and ultrapotassic magmatism in SW Tibet: geochemical and Sr-Nd-Pb-O isotopic constraints for mantle source characteristics and petrogenesis[END_REF]Pecerillo and Martinotti, 2006].

-As this magmatism indicates a hydrous melt, some authors [e.g. [START_REF] Stampfli | The Intra-Alpine terrain: a paleo-tethyan remnant in the Alpine Variscides[END_REF] proposed that it was generated in an active subduction or back-arc setting. However no subduction was active in this part of the Variscan belt at that time, the studied Stephanian volcanism can therefore not be directly linked to subduction processes.

-More frequently, the geochemical patterns of potassic late-orogenic magmatism have been explained by contamination, either by crustal contamination during magma ascent, or more probably by mantle metasomatism before melting, as indicated by the systematic and homogeneous character of contamination. Up to now, this metasomatism has been considered to be inherited from a suprasubduction position, either in the same orogenic cycle, or during previous subductions [e.g. [START_REF] Miller | Postcollisional potassic and ultrapotassic magmatism in SW Tibet: geochemical and Sr-Nd-Pb-O isotopic constraints for mantle source characteristics and petrogenesis[END_REF].

The southern part of the Variscan orogeny results from continental collision following a northward subduction, whose remnants are found in the French Massif central [START_REF] Matte | La chaîne varisque parmi les chaînes paléozoïques péri atlantiques, modèle d'évolution et position des grands blocs continentaux au Permo-Carbonifère[END_REF][START_REF] Matte | -The Variscan collage and orogeny (480-290 Ma) and the tectonic definition of the Armorica microplate: a review[END_REF]. However, the southernmost part of the French Massif central as well as the Pyrénées belong to the subducting plate, and its underlying mantle was not in a suprasubduction position, and could not be metasomatised by a Variscan subduction. However, this does not exclude older (Precambrian) metasomatism processes.

If the mantle was previously metasomatised, the main mechanism proposed to explain its melting in a post-orogenic setting is lithospheric thinning, inducing mantle decompression [e.g. [START_REF] Gibson | Ultrapotassic magmas along the flank of the Oligomiocene Rio-Grande rift, USA -Monitors of the zone of lithospheric mantle extension and thinning beneath a continental rift[END_REF][START_REF] Innocent | Sr-Nd isotope and trace-element geochemistry of late Variscan volcanism in the Pyrénées: magmatism in post-orogenic extension?[END_REF][START_REF] Rottura | Origin and significance of the Permian high-K calc-alkaline magmatism in the central-eastern Southern Alps, Italy[END_REF]]. However, the Stephanian volcanism occurred in a strike-slip tectonic regime in the French Massif central, and before the main Permian extension that ended the variscan orogenic cycle [Ménard and Molnar, 1988].

-We propose an alternative mechanism, controlled by lithospheric-scale strike-slip faulting, to explain both mantle metasomatism and mantle melting generating the potassic late-orogenic volcanism. At the end of the orogeny, the various parts of the mountain belt differ by their crustal and lithospheric thicknesses, thermal gradients, rock nature and geological history. As observed in the Ivrea body, hydrous metasediments can occur in the lower crust at the end of the Variscan orogeny [START_REF] Harlov | High-grade fluid metasomatism on both a local and a regional scale: the Seward Peninsula, Alaska, and the Val Strona di Omagna, Ivrea-Verbano zone, norther, Italy. Part I: Petrography and silicate mineral chemistry[END_REF].

As the late-orogenic strike-slip faults cross the orogenic zonation with some tens or hundreds kilometres of horizontal slipping, these bring into contact different parts of the mountain belt. For example, this horizontal displacement can bring into contact, and at the same depth, mantle materials on one side and lower crust, including still hydrated metasediments, on the other. The crustal fluids could then locally hydrate the adjacent mantle, and this metasomatism may induce mantle melting. The resulting magmas have geochemical patterns that indicate contamination by continental crust.

In this hypothetical mechanism, mantle metasomatism is not inherited from a previous suprasubduction zone, but is acquired by horizontal fluid transfer in or across lithosphericscale strike-slip faults. Consequently, mantle melting can occur without vertical displacements and adiabatic decompression, since metasomatism controlled by strike-slip faults drastically decreases its solidus temperature. Because fluid transfers across strike-slip faults and induced mantle metasomatism are only local and can occur only once, this mechanism can explain why this late-orogenic potassic magmatism represents small volumes of lava, and why it happens only at the beginning of late-orogenic strike-slip faulting. [START_REF] Winchester | Geochemical discrimination of different magma series and their differentiation products using immobile elements[END_REF]. [START_REF] Winchester | Geochemical discrimination of different magma series and their differentiation products using immobile elements[END_REF]. [START_REF] Müller | Geochemical discrimination between shoshonitic and potassic volcanic rocks in different tectonic settings: a pilot study[END_REF]. [START_REF] Müller | Geochemical discrimination between shoshonitic and potassic volcanic rocks in different tectonic settings: a pilot study[END_REF]. [START_REF] Sun | Chemical and isotopic systematics of oceanic basalts: implications for mantle composition and processes[END_REF].

Figure Captions

Sample sets as a function of location (cf. Fig. 1).

Figure 5: Spectres de terres rares normalisés aux chondrites [START_REF] Sun | Chemical and isotopic systematics of oceanic basalts: implications for mantle composition and processes[END_REF]. Les échantillons sont regroupés en fonction de leur localisation (cf. Fig. 1).

Figure 6: Primitive mantle-normalized spidergrams [START_REF] Sun | Chemical and isotopic systematics of oceanic basalts: implications for mantle composition and processes[END_REF]. Sample sets as a function of location (cf. Fig. 1).

Figure 6: Spidergrams normalisés au manteau primitif [START_REF] Sun | Chemical and isotopic systematics of oceanic basalts: implications for mantle composition and processes[END_REF]. Les échantillons sont regroupés en fonction de leur localisation (cf. Fig. 1).

Table Captions

Table I: Major-and trace-element analyses, 87 Sr/ 86 Sr and 143 Nd/ 144 Nd isotope ratios for studied samples. Major oxides are in wt%, Trace-elements are in ppm. Fe 2 O 3 : total iron as ferric oxide; LOI: loss on ignition; bdl: below detection limit; nd: not determined. Rock types are given from figure 3. Location on figure 1 (LM for Lacapelle-Marival).

Table I: Analyses chimiques en éléments majeurs et traces, rapports isotopiques en

(

  CLB95-8, CLB95-9) differ from the other rocks by flatter patterns [2.8 < (La/Yb) n < 4.9]. The trachyte from the Livinhac site (CLB95-15) differs from the other felsic rocks by a very steep REE pattern due to an important HREE depletion [(La/Yb) n = 31.9].

  Sm, Ti anomalies and positive Th, U and Pb anomalies. The most SiO 2 -rich and MgO-poor rocks and one basalt (CLB95-17) can be distinguished from other samples by a negative correlation in Ti/Zr diagram (figs. 2 & 3). These samples present additional geochemical characteristics, such as Y depletion and enhanced negative Ti anomalies, and for some samples negative Ce anomalies.

Figeac

  Figeac and Lacapelle-Marival magmas.

Figure 1 :

 1 Figure 1: Schematic geologic map of the area between Lacapelle-Marival and Decazeville

Figure 1 :

 1 Figure 1: Carte géologique schématique de la région comprise entre les bassins de

Figure 2 :

 2 Figure 2: Variations of minor-(TiO 2 ) and trace-elements (Th, U, Nb, La, Pb, Rb, Sr, Y)

Figure 2 :

 2 Figure 2: Variations des teneurs en éléments mineur (TiO 2 ) et traces (Th, U, Nb, La, Pb,

Fig

  Fig.1).

Figure 3 :

 3 Figure 3: Zr/Ti versus Nb/Y discrimination diagram[START_REF] Winchester | Geochemical discrimination of different magma series and their differentiation products using immobile elements[END_REF].

Figure 3 :

 3 Figure 3: Diagramme discriminant Zr/Ti par rapport à Nb/Y[Winchester and Floyd,

Figure 4 :

 4 Figure 4: Th/Yb versus Ta/Yb discrimination diagram[START_REF] Müller | Geochemical discrimination between shoshonitic and potassic volcanic rocks in different tectonic settings: a pilot study[END_REF].

Figure 4 :

 4 Figure4: Diagramme discriminant Th/Yb par rapport à Ta/Yb[START_REF] Müller | Geochemical discrimination between shoshonitic and potassic volcanic rocks in different tectonic settings: a pilot study[END_REF].

Figure 5 :

 5 Figure 5: Chondrite-normalized rare earth elements patterns[START_REF] Sun | Chemical and isotopic systematics of oceanic basalts: implications for mantle composition and processes[END_REF]].

Figure 7 :

 7 Figure 7: Variations of Ti/Y versus Zr/Y (left) and Nb/La versus Th/La (right). The

Figure 7 :

 7 Figure 7: Variations des rapports Ti/Y en fonction de Zr/Y (à gauche) et Nb/La en fonction

  Figeac

  Figeac

  Th U Nb Ta La Ce Pb Pr Sr P Nd Sm Zr Hf Eu Ti Gd Dy Y Yb Th U Nb Ta La Ce Pb Pr Sr P Nd Sm Zr Hf Eu Ti Gd Dy Y Yb Th U Nb Ta La Ce Pb Pr Sr P Nd Sm Zr Hf Eu Ti Gd Dy Y Yb Th U Nb Ta La Ce Pb Pr Sr P Nd Sm Zr Hf Eu Ti Gd Dy Y Yb Th U Nb Ta La Ce Pb Pr Sr P Nd Sm Zr Hf Eu Ti Gd Dy Y Yb

	Figure 5
				CLB96-13
	1000	Lacapelle-Marival North	CLB96-13 CLB96-14 CLB96-16
				CLB96-14 CLB96-17
				CLB96-16
	Rock / Chondrite	10 100		CLB96-17
		1
	La Ce P r Nd Sm Eu Gd Tb Dy Ho Er Tm Yb Lu 1000 CLB96-6 to -12 1000	CLB95-13 CLB95-15
			Lacapelle-Marival South	CLB96-6 CLB96-10			Livinhac	CLB95-9 CLB95-15 CLB95-9 to -12 CLB95-17
				CLB96-11			CLB95-10 to -14	CLB95-14
	Rock / Chondrite	10 100		CLB96-7 to -9 CLB96-12	10 100		CLB95-16 CLB96-1 CLB95-17	CLB95-16 CLB96-1
		1		Lacapelle-Marival North Lacapelle-Marival South 1	Bourran Livinhac Figeac
		La Ce P r Nd Sm Eu Gd Tb Dy Ho Er Tm Yb Lu CLB96-5 La Ce P r Nd Sm Eu Gd Tb Dy Ho Er Tm Yb Lu CLB95-1
	1000			1000 CLB96-2	CLB95-2
			Figeac	CLB96-2 CLB96-3	CLB96-3	Bourran	CLB95-1 CLB95-2 CLB95-4 to -8
				CLB96-5			CLB95-4
	Rock / Chondrite	10 100	Th/Yb		10 100		CLB95-5 CLB95-6 CLB95-7 CLB95-8
		1			1
		La Ce P r Nd Sm Eu Gd Tb Dy Ho Er Tm Yb Lu	La Ce P r Nd Sm Eu Gd Tb Dy Ho Er Tm Yb Lu
				1			1	10
				Ta/Yb

Acknowledgements

D. Arnaud, F. Senebier and F. Coeur are thanked for sample preparation, F. Keller and L.

Briqueu for the chemical analyses, F. Chauvet for helpfull discussions. We are particularly grateful to A. Demant and M. Faure for their reviews and comments, and to J. Allibon for his english tips and corrections. This work has been funded by the LGCA.


Factorization Algorithm

A MULTI-LEVEL BLOCKING DISTINCT DEGREE

Keywords: Complexité amortie, corps fini, exposant de Mersenne, factorisation en degrés distincts, factorisation de polynôme, trinôme irréductible, trinôme primitif

Un algorithme multi-étage pour la factorisation en degrés distincts de polynômes Résumé : Nous proposons un nouvel algorithme pour la factorisation en degrés distincts d'un polynôme P (x) sur GF(2), via une stratégie multi-étage. Le niveau supérieur remplace des calculs de pgcd par des multiplications, comme suggéré par Pollard (1975), von zur Gathen et Shoup (1992), et d'autres auteurs. L'originalité de notre approche tient dans un niveau inférieur, qui remplace des multiplications par des carrés, ce qui accélère le calcul de certains polynômes intervalles sur GF(2)[x]/P (x), quand P (x) est creux.

Comme application nous exhibons un algorithme rapide cherchant tous les trinômes irréductibles x r + x s + 1 de degré r sur GF(2), tout en produisant un certificat qui peut être vérifié plus rapidement. Les algorithmes naïfs coûtent O(r 2 ) par trinôme, soit O(r 3 ) pour tous les trinômes de degré r. Sous une hypothèse naturelle sur la distribution des facteurs de trinômes, le nouvel algorithme a une complexité O(r 2 (log r) 3/2 (log log r) 1/2 ) pour tester tous les trinômes de degré r. Notre implantation est 560 fois plus rapide que l'algorithme naïf dans le cas r = 24036583 (exposant de Mersenne).

Avec notre programme, nous avons trouvé deux nouveaux trinômes primitifs de degré 24036583 sur GF(2), le précédent record étant de degré 6972593.

Introduction

The problem of factoring a univariate polynomial P (x) over a finite field F often arises in computational algebra [START_REF] Cantor | A new algorithm for factoring polynomials over finite fields[END_REF][START_REF] Zur Gathen | Polynomial factorization over F 2[END_REF][START_REF] Gathen | Computing Frobenius maps and factoring polynomials[END_REF]. An important case is when F has small characteristic and P (x) has high degree but is sparse, that is P (x) has only a small number of nonzero terms.

To simplify the exposition we restrict attention to the case where F = GF(2) and P (x) is a trinomial P (x) = x r + x s + 1, r > s > 0, although the ideas apply more generally and should be useful for factoring sparse polynomials over fields of small characteristic.

Our aim is to give an algorithm with good amortized complexity, that is, one that works well on average. Since we are restricting attention to trinomials, we average over all trinomials of fixed degree r.

Our motivation is to speed up previous algorithms for searching for irreducible trinomials of high degree [START_REF] Brent | A fast algorithm for testing reducibility of trinomials mod 2 and some new primitive trinomials of degree 3021377[END_REF][START_REF] Brent | A primitive trinomial of degree 6972593[END_REF][START_REF] Kumada | New primitive t-nomials (t = 3, 5) over GF(2) whose degree is a Mersenne exponent[END_REF]. For given degree r, we want to find all irreducible trinomials x r + x s + 1.

In our examples the degree r is a Mersenne exponent, i.e., 2 r -1 is a Mersenne prime. In this case an irreducible trinomial of degree r is necessarily primitive. In general, without the restriction to Mersenne exponents, we would need the prime factorisation of 2 r -1 in order to test primitivity (see e.g., [START_REF] Zur Gathen | Modern Computer Algebra[END_REF]).

We are only interested in Mersenne exponents r = ±1 mod 8, because in other cases Swan's theorem [START_REF] Pellet | Sur la décomposition d'une fonction entière en facteurs irréductibles suivant un module premier p[END_REF][START_REF] Stickelberger | Über eine neue Eigenschaft der Diskriminanten algebraischer Zahlkörper[END_REF][START_REF] Swan | Factorization of polynomials over finite fields[END_REF] rules out irreducible trinomials of degree r (except for s = 2 or r -2, but these cases are usually easy to handle: for example if r = 13466917 or 20996011 we have r = 1 mod 3, so x r + x 2 + 1 is divisible by x 2 + x + 1).

Mersenne exponents can be found on the GIMPS website [START_REF] Woltman | Technopôle de Nancy-Brabois -Campus scientifique 615, rue du Jardin Botanique -BP 101 -54602 Villers-lès-Nancy Cedex[END_REF]. At the time of writing, the five largest known Mersenne exponents r satisfying the condition r = ±1 mod 8 are r = 6972593, 24036583, 25964951, 30402457 and 32582657. In the smallest case r = 6972593, a primitive trinomial was found by Brent, Larvala and Zimmermann [START_REF] Brent | A primitive trinomial of degree 6972593[END_REF] using an efficient implementation of the naive algorithm. However, it was not feasible to consider the larger Mersenne exponents r using the same algorithm, since the time complexity of this algorithm is roughly of order r 3 , and the next case r = 24036583 would take about 41 times longer than r = 6972593. With the new "fast" algorithm described in this paper we have been able to find two primitive trinomials of degree r = 24036583 in less time than the naive algorithm took for r = 6972593. The speedup over the naive algorithm for r = 24036583 is about a factor of 560.

If x r +x s +1 is reducible then we want to provide an easily-checked certificate of reducibility. The certificate can simply be an encoding of an irreducible factor f of x r + x s + 1. We choose the factor f of smallest degree d > 0. In case there are several factors of equal smallest degree d, we give the one that is least in lexicographic order, e.g.,

x 3 + x + 1 is preferred to x 3 + x 2 + 1.
1.1. Distinct degree factorization. Our basic algorithm performs distinct degree factorization [START_REF] Ph | The complete analysis of a polynomial factorization algorithm over finite fields[END_REF][START_REF] Zur Gathen | Modern Computer Algebra[END_REF][START_REF] Zur Gathen | Polynomial factorization over F 2[END_REF]. That is, if P (x) has several factors of the same degree d, the algorithm will produce the product of these factors. The Cantor-Zassenhaus algorithm is used to split this product into distinct factors. This is cheap because the product usually consists of just one irreducible factor or is a product of irreducible factors of small (equal) degree.

In the complexity analysis we only consider the time required to find one nontrivial factor (it will be a factor of smallest degree) or output "irreducible", since that is what is required in the search for irreducible trinomials.

1.2. Factorization over GF [START_REF] Bosma | Handbook of Magma Functions[END_REF]. It is well-known that x 2 d + x is the product of all irreducible polynomials of degree dividing d. For example,

x 2 3 + x = x(x + 1)(x 3 + x + 1)(x 3 + x 2 + 1).
Thus, a simple algorithm to find a factor of smallest degree of P (x) is to compute GCD(x 2 d + x, P (x)) for d = 1, 2, . . . The first time that the GCD is nontrivial, it contains a factor of minimal degree d. If the GCD has degree > d, it must be a product of factors of degree d. If no factor has been found for d ≤ r/2, where r = deg(P (x)), then P (x) must be irreducible. Some simplifications are possible when P (x) = x r + x s + 1 is a trinomial over GF(2) with r or s odd (otherwise P (x) is trivially reducible):

(1) We can skip the case d = 1 because a trinomial can not have a factor of degree 1.

(2) Since x r P (1/x) = x r + x r-s + 1, we only need consider s ≤ r/2.

(3) We can assume that P (x) is square-free. (4) By applying Swan's theorem, we can often show that the trinomial under consideration has an odd number of irreducible factors; in this case we only need check d ≤ r/3 before claiming that P (x) is irreducible.

Complexity of the algorithm

Note that x 2 d should not be computed explicitly; it is much better to compute x 2 d mod P (x) by repeated squaring. The complexity of squaring modulo a trinomial of degree r is only S(r) = O(r) bit-operations. INRIA 2.1. Complexity of polynomial multiplication and squaring. As well as performing GCD computations we need to perform multiplications in GF(2)[x]/P (x), and an important special case is squaring a polynomial modulo P (x), so we first consider the bit-complexity of these operations.

Multiplication of polynomials of degree r over GF(2) can be performed in time M (r) = O(r log r log log r). We have implemented an algorithm of Schönhage [START_REF] Schönhage | Schnelle Multiplikation von Polynomen über Körpern der Charakteristik 2[END_REF] that achieves this bound. The algorithm uses a radix-3 FFT and is different from the better-known Schönhage-Strassen algorithm [START_REF] Schönhage | Schnelle Multiplikation groβer Zahlen[END_REF]. We remark that the log log r term in the time-bound for the Schönhage-Strassen algorithm has been reduced by Fürer [START_REF] Fürer | Faster integer multiplication[END_REF], but it is not clear if a similar idea can be used to improve Schönhage's algorithm [START_REF] Schönhage | Schnelle Multiplikation von Polynomen über Körpern der Charakteristik 2[END_REF]. In any event the log log r term comes from the number of levels of recursion and is a small constant for the values of r that we are considering.

In practice, Schönhage's algorithm is not the fastest unless r is quite large. We have also implemented classical, Karatsuba and Toom-Cook algorithms that have M (r) = O(r α ), 1 < α ≤ 2, since these algorithms are easier to implement and are faster for small r. Our implementations of the Toom-Cook algorithms TC3 and TC4 are based on recent ideas of Bodrato [START_REF] Bodrato | Towards Optimal Toom-Cook Multiplication for Univariate and Multivariate Polynomials in Characteristic 2 and 0[END_REF].

For brevity we assume that r is large and Schönhage's algorithm is used. On a 64-bit machine the crossover versus TC4 occurs near degree r = 108000.

In the complexity estimates we assume that M (r) is a sufficiently smooth and wellbehaved function.

By Squaring we mean squaring a polynomial of degree < r and reduction mod P (x). Squaring in GF(2)[x]/P (x) can be performed in time S(r) = Θ(r) ≪ M (r) (assuming, as usual, that P (x) is a trinomial). Our algorithm takes advantage of the fact that squaring is much faster than multiplication.

Where possible we use the memory-efficient squaring algorithm of Brent, Larvala and Zimmermann [START_REF] Brent | A fast algorithm for testing reducibility of trinomials mod 2 and some new primitive trinomials of degree 3021377[END_REF], which in our implementation is about 2.2 times faster than the naive squaring algorithm.

Complexity of GCD. For GCDs we use a sub-quadratic algorithm that runs in time

G(r) = Θ(M (r) log r). More precisely, G(2r) = 2G(r) + Θ(M (r)), so for α > 1, M (r) = Θ(r α ) ⇒ G(r) = Θ(M (r)),
and

M (r) = Θ(r log r log log r) ⇒ G(r) = Θ(M (r) log r).
In practice, for r ≈ 2.4 × 10 7 and our implementation on a 2.2 Ghz Opteron, S(r) ≈ 0.005 second, M (r) ≈ 2 seconds, G(r) ≈ 80 seconds, so M (r)/S(r) ≈ 400, and G(r)/M (r) ≈ 40.

Avoiding GCD computations.

In the context of integer factorization, Pollard [START_REF] Pollard | A Monte Carlo method for factorization[END_REF] suggested a blocking strategy to avoid most GCD computations and thus reduce the amortized cost; von zur Gathen and Shoup [START_REF] Gathen | Computing Frobenius maps and factoring polynomials[END_REF] applied the same idea to polynomial factorization.

The idea of blocking is to choose a parameter ℓ > 0 and, instead of computing GCD(x

2 d + x, P (x)) for d ∈ [d ′ , d ′ + ℓ), compute GCD(p ℓ (x 2 d ′ , x), P (x))
, where the interval polynomial p ℓ (X, x) is defined by

p ℓ (X, x) = ℓ-1 j=0 X 2 j + x .
In this way we replace ℓ GCDs by one GCD and ℓ -1 multiplications mod P (x).

The drawback of blocking is that we may have to backtrack if P (x) has more than one factor with degree in the interval [d ′ , d ′ + ℓ), since the algorithm produces the product of these factors. Thus ℓ should not be too large. The optimal strategy depends on the expected size distribution of factors and the ratio of times for GCDs and multiplications.

2.4. Multi-level blocking. Our (apparently new) idea is to use a finer level of blocking to replace most multiplications by squarings, which speeds up the computation in GF(2)[x]/P (x) of the above interval polynomials. The idea is to split the interval

[d ′ , d ′ + ℓ) into k ≥ 2 smaller intervals of length m over which (1) p m (X, x) = m-1 j=0 X 2 j + x = m j=0
x m-j s j,m (X), where

(2) s j,m (X) = 0≤k<2 m , w(k)=j X k ,
and w(k) denotes the Hamming weight of k, that is the number of nonzero bits in the binary representation of k.

For example, for m = 3, we have:

p m (X, x) = x 3 + x 2 (X 4 + X 2 + X) + x(X 6 + X 5 + X 3 ) + X 7 , where s 0,3 (X) = 1, s 1,3 (X) = X 4 + X 2 + X, s 2,3 (X) = X 6 + X 5 + X 3 , and s 3,3 (X) = X 7 . Note that s j,m (X 2 ) = s j,m (X) 2 in GF(2)[x]/P (x).
Thus, p m (x 2 d , x) can be computed with cost m 2 S(r) if we already know s j,m (x 2 d-m ) for 0 < j ≤ m. (The constant polynomial s 0,m (X) = 1 is computed only once.) Continuing the example with m = 3, and assuming that we know s 1,3 (x 2 d-3 ), s 2,3 (x 2 d-3 ), and s 3,3 (x 2 d-3 ), squaring each of these m = 3 times gives s 1,3 (x 2 d ), s 2,3 (x 2 d ), and s 3,3 (x 2 d ), from which we can easily get p 3 (x 2 d , x) using the sum in Eq. [START_REF] Bodrato | Towards Optimal Toom-Cook Multiplication for Univariate and Multivariate Polynomials in Characteristic 2 and 0[END_REF].

In this way we replace m -1 multiplications and m squarings -if we used the product in Eq. (1) -by m 2 squarings. Each s j,m , 0 < j ≤ m, requires m squarings to be shifted from argument x Von zur Gathen and Gerhard [11, p. 1685] suggested using the same idea with m = 2 (thus reducing the number of multiplications by a factor of two), but did not consider choosing an optimal m > 2.

At first sight initialization of the polynomials s j,m (X) for X = x might appear to be expensive, since the definition (2) involves O(2 m ) terms. However, the polynomials s j,m (X) satisfy a "Pascal triangle" recurrence relation

s j,m (X) = s j,m-1 (X 2 ) + Xs j-1,m-1 (X 2 ) with boundary conditions s j,m (X) = 0 if j > m ≥ 0, 1 if m ≥ j = 0.
Using this recurrence, it is easy to compute s j,m (x) mod P (x) for 0 ≤ j ≤ m in time O(m 2 r). Thus, the initialization is cheap.

To summarise, we use two levels of blocking:

(1) The outer level replaces most GCDs by multiplications.

(2) The inner level replaces most multiplications by squarings.

(3) The parameter m ≈ M (r)/S(r) is used for the inner level of blocking. (4) A different parameter ℓ = km is used for the outer level of blocking. For example, suppose S = 1/400, M = 1, G = 40 (where we have normalised so M = 1). We could choose ℓ = 80 and m = 20. With no blocking, the cost for an interval of length 80 is 80G + 80S = 3200.2; with 1-level blocking the cost is G + 79M + 80S = 119.2; with 2-level blocking the cost is G + 3M + 1600S = 47.0.

2.5.

Sieving out small factors. We define a small factor to be one with degree d < 1 2 log 2 r, so 2 d < √ r. The constant 1 2 in the definition is arbitrary and could be replaced by any fixed constant in (0, 1). A large factor is a factor that is not small.

It would be inefficient to find small factors in the same way as large factors. Instead, let D = 2 d -1, r ′ = r mod D, s ′ = s mod D. Then P (x) = x r + x s + 1 = x r ′ + x s ′ + 1 mod (x D -1), so we only need compute GCD(x r ′ + x s ′ + 1, x D -1). Because r ′ , s ′ < D < √ r, the cost of finding small factors is negligible (both theoretically and in practice), so can be neglected.

2.6.

Outer level blocking strategy. The blocksize in the outer level of blocking is ℓ = km. We take a linearly increasing sequence of block sizes

k = k 0 j for j = 1, 2, 3, . . . ,
where the first interval starts at about log r (since small factors will have been found by sieving).

The choice k = k 0 j leads to a quadratic polynomial for the interval bounds; other possibilities are discussed by von zur Gathen and Gerhard [START_REF] Zur Gathen | Polynomial factorization over F 2[END_REF].

In principle, using the data that we have obtained on the distribution of degrees of smallest factors of trinomials (see 3), and assuming that this distribution is not very sensitive to the degree r, we could obtain a strategy that is close to optimal. However, the choice k 0 j with suitable k 0 is easy to implement and not too far from optimal. The number of GCD and sqr/mul operations is usually within a factor of 1.5 of the minimum possible in our experiments.

Distribution of degrees of factors

In order to predict the expected behaviour of our algorithm, we need to know the expected distribution of degrees of smallest irreducible factors. From Swan's theorem [START_REF] Swan | Factorization of polynomials over finite fields[END_REF], we know that there are significant differences between the distribution of factors of trinomials and of all polynomials of the same degree. Our complexity estimates are based on the heuristic assumption that this difference is not too large, in a sense made precise by Hypothesis 3.1.

Hypothesis 3.1. Over all trinomials x r + x s + 1 of degree r over GF(2), the probability π d that a trinomial has no nontrivial factor of degree ≤ d, 1 < d ≤ r, is at most c/d, where c is a constant. Hypothesis 3.1 implies that there are at most c irreducible trinomials of degree r. This is probably false, as there may well be a sequence of exceptional r for which the number of irreducible trinomials is unbounded. Thus, we may need to replace the constant c in Hypothesis 3.1 by a slowly-growing function c(r). Nevertheless, in order to give realistic complexity estimates that are in agreement with experiments, we assume below that Hypothesis 3.1 is correct. Under this assumption we use an amortized model to obtain the total complexity over all trinomials of degree r.

From Hypothesis 3.1, the probability that a trinomial does not have a small factor (as defined in 2.5) is O(1/ log r).

Table 1 gives the observed values of dπ d for r = 3021377, r = 6972593, and r = 24036583. The maximum values for each r are given in bold. The table shows that the values of dπ d are remarkably stable for small d, and bounded by 4 for large d (this is because there are four irreducible trinomials of degree 3021377 and also four of degree 24036583, when we count both trinomials x r + x s + 1 and their reciprocals x r + x r-s + 1).

3.1. Consequences of the hypothesis. Define p k = π d-1 -π d to be the probability that the smallest nontrivial factor f of a randomly chosen trinomial has degree d = deg(f ). In order to estimate the running time of our algorithm, we use the following Lemma, which gives the expectation E β of d β . Lemma 3.2. If β > 0 is constant and Hypothesis 3.1 holds, then Proof. We use summation by parts. Note that a trinomial has no factor of degree 1, so p 1 = 0 and π 0 = π 1 = 1. Thus

E β := r d=1 d β p d =    O(1) if β < 1, O(log r) if β = 1, O(r β-1 ) if β > 1.

INRIA

E β = r d=1 d β p d = r d=1 d β (π d-1 -π d ) = r-1 d=1 (d + 1) β -d β π d + π 0 -r β π r ≤ 1 + c r-1 d=1 (d + 1) β -d β d
(by Hypothesis 3.1)

≤ 1 + O r-1 d=1 d β-2
and the result follows.

The following Lemma gives a stronger result in the case β < 1.

Lemma 3.3. If 0 < β < 1, 0 < D ≤ r,

RR n 6331

Proof. The proof is similar to that of Lemma 3.2. We end with the upper bound r-1 d=D

(d + 1) β -d β d + D β π D-1 .
From Hypothesis 3.1, π D-1 = O(1/D), and the sum over d is O(D β-1 ), so the result follows.

Expected cost of sqr/mul and GCD

Recall that the inner level of blocking replaces m multiplications by m 2 squarings and one multiplication, where the choice m ≈ M (r)/S(r) makes the total cost of squarings about equal to the cost of multiplications.

For a smallest factor of degree d, the number of squarings is m(d

+ O( √ d))
, where the O( √ d) term follows from our choice of outer-level blocksizes (see 2.6). Averaging over all trinomials of degree r, the expected number of squarings is

O   m d≤r/2 (d + O( √ d))p d   ,
and from Lemma 3.2 this is O(m log r). Thus, the expected cost of sqr/mul operations per trinomial is

O S(r) log r M (r)/S(r) = O log r M (r)S(r) = O r(log r) 3/2 (log log r) 1/2 . ( 3
)
If we used only a single level of blocking, then the cost of multiplications would dominate that of squarings, with an expected cost per trinomial of O (log rM (r)) = O r(log r) 2 log log r .

(3) is correct as r → ∞. However, in practice, at least for r < 6.4 × 10 7 , our implementation of Schönhage's FFT-based polynomial multiplication algorithm [START_REF] Schönhage | Schnelle Multiplikation von Polynomen über Körpern der Charakteristik 2[END_REF] calls a different multiplication routine (usually TC4) to perform smaller multiplications, rather than recursively calling itself. TC4 has exponent α ′ = ln(7)/ ln(4) ≈ 1.4, so the effective exponent for FFT multiplication is α = (1 + α ′ )/2 ≈ 1.2 > 1. In this case, the expected cost of sqr/mul operations per trinomial is [START_REF] Brent | A fast algorithm for testing reducibility of trinomials mod 2 and some new primitive trinomials of degree 3021377[END_REF] O log r M (r)S(r) = O(r (1+α)/2 log r) = O(r 1. Thus the expected cost of GCDs per trinomial is (5) is asymptotically less than the expected cost (3) of sqr/mul operations. However, if M (r) = O(r α ) with α > 1, then the expected cost of GCDs is O(r α / log r), which is asymptotically greater than the expected cost (4) of sqr/mul operations. Note the expected cost of GCDs does not depend on whether we use one or two levels of blocking. For r ≈ 2.4 × 10 7 , GCDs take about 65% of the time versus 35% for sqr/mul.

4.2.

Comparison with previous algorithms. For simplicity we use the O notation which ignores log factors. For example, M (r) = O(r).

The "naive" algorithm, as implemented by Brent, Larvala and Zimmermann [START_REF] Brent | A fast algorithm for testing reducibility of trinomials mod 2 and some new primitive trinomials of degree 3021377[END_REF] and earlier authors, takes an expected time O(r 2 ) per trinomial, or O(r 3 ) to cover all trinomials of degree r.

The single-level blocking strategy and the new algorithm both take expected time O(r) per trinomial, or O(r 2 ) to cover all trinomials of degree r.

In practice, the new algorithm is faster over the naive algorithm by a factor of about 160 for r = 6972593, and by a factor of about 560 for r = 24036583. For r = 24036583, where sqr/mul operations take 35% of the total time in the new algorithm, and the corresponding speedup is about 10, this gives a global speedup of more than 4 over the single-blocking strategy.

4.3. Some details of our implementation. We first implemented the 2-level blocking strategy in NTL [START_REF] Shoup | NTL: A library for doing number theory[END_REF]. To get full efficiency, we rewrote all critical routines and tuned them efficiently on the target processors. Our squaring routine implements the algorithm described in [START_REF] Brent | A fast algorithm for testing reducibility of trinomials mod 2 and some new primitive trinomials of degree 3021377[END_REF], which is more than twice as fast as the corresponding optimized NTL routine for trinomials. Our multiplication routine implements Toom-Cook 3-way, 4way, and Schönhage's algorithm [START_REF] Schönhage | Schnelle Multiplikation von Polynomen über Körpern der Charakteristik 2[END_REF]. We also improved the basecase multiplication code; more details concerning efficient multiplication in GF(2)[x] will be published in [START_REF] Brent | Efficient Multiplication in GF (2)[END_REF]. Finally, we implemented a subquadratic GCD routine, since NTL only provides a classical GCD for binary polynomials. 4.4. Primitive trinomials. The largest published primitive trinomial is

x 6972593 + x 3037958 + 1,
found by Brent, Larvala and Zimmermann [START_REF] Brent | A fast algorithm for testing reducibility of trinomials mod 2 and some new primitive trinomials of degree 3021377[END_REF] in 2002 using a naive (but efficiently implemented) algorithm.

In March-April 2007, we tested our new program by verifying the published results on primitive trinomials for Mersenne exponents r ≤ 6972593, and in the process produced certificates of reducibility (lists of smallest factors for each reducible trinomial). These are available from the first author's website [START_REF] Brent | Search for primitive trinomials (mod 2)[END_REF].

In April-August 2007, we ran our new algorithm to search for primitive trinomials of degree r = 24036583. This is the next Mersenne exponent, apart from two that are trivial to exclude by Swan's theorem. It would take about 41 times as long as for r = 6972593 by the naive algorithm, but our new program is 560 times faster than the naive algorithm. Each trinomial takes on average about 16 seconds on a 2.2 Ghz Opteron.

The complete computation was performed in four months, using about 24 Opteron and Core 2 processors located at ANU and INRIA.

We found two new primitive trinomials of (equal) record degree:

(6)

x 24036583 + x 8412642 + 1 and (7)

x 24036583 + x 8785528 + 1.

4.5. Verification. Allan Steel [START_REF] Steel | [END_REF] kindly verified irreducibility of ( 6)-( 7) using Magma [START_REF] Bosma | Handbook of Magma Functions[END_REF].

Each verification took about 67 hours on an 2.4 GHz Core 2 processor. Independent verifications using our irred V3.15 program [START_REF] Brent | A fast algorithm for testing reducibility of trinomials mod 2 and some new primitive trinomials of degree 3021377[END_REF][START_REF] Brent | A primitive trinomial of degree 6972593[END_REF] took about 35 hours on a 2.2 Ghz Opteron.

The difference in speed is mainly due to the fast squaring algorithm implemented in irred.

Primitivity of ( 6)-( 7) follows from irreducibility provided that the degree 24036583 is a Mersenne exponent. We have not verified this, but rely on computations performed by the GIMPS project [START_REF] Woltman | Technopôle de Nancy-Brabois -Campus scientifique 615, rue du Jardin Botanique -BP 101 -54602 Villers-lès-Nancy Cedex[END_REF].

Reducibility of the remaining trinomials of degree 24036583 can be verified using the certificate (or extended log, a list of smallest irreducible factors) available from our website [START_REF] Brent | Search for primitive trinomials (mod 2)[END_REF]. The verification takes less than 10 hours using Magma on a 2.66 Ghz Core 2 processor.

Conclusion

The new double-blocking strategy, combined with fast multiplication and GCD algorithms, has allowed us to find new primitive trinomials of record degree.

The same ideas should work over finite fields GF(p) for small prime p > 2, and for factoring sparse polynomials P (x) that are not necessarily trinomials: all we need is that the time for p-th powers (mod P (x)) is much less than the time for multiplication (mod P (x)).

  2 d-m to argument x 2 d . The summation in Eq. (1) costs only O(mr), which INRIA is negligible. Choosing m ≈ M (r)/S(r) (about 20 if M (r)/S(r) ≈ 400), the speedup over single-level blocking is about m/2 ≈ 10 (not counting the cost of GCDs).

  and Hypothesis 3.1 holds, then r d=D d β p d = O D β-1 .

( 5 )

 5 O(G(r)/ log r) = O(M (r)) = O(r log r log log r).

Table 1 .

 1 dπ d for various degrees r.

	d	r = 3021377 r = 6972593 r = 24036583
	2	1.333	1.333	1.333
	3	1.429	1.429	1.429
	4	1.524	1.524	1.524
	5	1.536	1.536	1.536
	6	1.598	1.598	1.598
	7	1.600	1.600	1.600
	8	1.667	1.667	1.667
	9	1.642	1.642	1.642
	10	1.652	1.652	1.652
	100	1.763	1.771	1.770
	1000	1.783	1.756	1.786
	10000	1.946	1.873	1.786
	100000	1.986	1.606	1.880
	279383	1.480	2.084	1.813
	1000000	1.324	1.147	1.831
	10000000	-	-	1.664
	r -1	4.000	2.000	4.000

  ••• log r) 4.1. Expected cost of GCDs. Suppose that P (x) has a smallest factor of degree d.

	and by Lemma 3.3 this is
	1 + O	1 log r √	.
					The
	number of GCDs required to find the factor, using our (quadratic polynomial) blocking strategy, is at least 1, and O( √ d) if d is large. By Hypothesis 3.1, the expected number of
	GCDs for a trinomial with no small factor is 			
	1 + O		d 1/2 p d	 ,
		log 2 r<2d≤r
					INRIA

RR n 6331

Acknowledgements. We thank Allan Steel for verifying irreducibility of the trinomials ( 6)- (7), and Marco Bodrato, Pierrick Gaudry and Emmanuel Thomé for their assistance in implementing fast algorithms for multiplication of polynomials over GF [2]. ANU and IN-RIA provided computing facilities. The first author's research was supported by MASCOS and the Australian Research Council.


François Cuvelier
email: cuvelier@math.univ-paris13.fr

Exact integration for products of power of barycentric coordinates over d-simplexes in R n

Exact integral computation over a d-simplex in R n for products of powers of its barycentric coordinates is done in [9] by using mathematical induction and coordinate mappings. In this note we give a new proof using Laplace transformations without mathematical induction.

Local shape functions of a large variety of nite element on a d-simplex K Ă R n can be expressed in function of the barycentric coordinates tλ 0 , . . . , λ d u of K and their derivatives (see [START_REF] Arnold | Geometric decompositions and local bases for spaces of nite element dierential forms[END_REF] for examples).

In [START_REF] Vermolen | On an integration rule for products of barycentric coordinates over simplexes in rn[END_REF], the authors give a proof of the magic formula: let ν ν ν " pν 0 , . . . , ν d q P N d`1 , then ż K d ź i"0 λ νi i pq q qqdq q q " d!|K|

d ś i"0 ν i ! pd `d ř l"0 ν i q! ( 1
)
where |K| is the volume of K. In their proof, mathematical induction and coordinate mappings are mainly used. In this note we give a new proof of this formula using Laplace transformations without mathematical induction. Firstly we recall denitions of a d-simplex in R n and of its barycentric coordinates. Therafter we introduce Laplace transforms to compute the volume of the unit d-simplex K Ă R d and the magic formula (1) over K. In the last section, we propose to compute the gradients of the barycentric coordinates by solving linear systems. We also present the mapping of an integral over a d-simplex in R n to the reference unit d-simplex, allowing to proove [START_REF] Arnold | Geometric decompositions and local bases for spaces of nite element dierential forms[END_REF]. [START_REF] Arnold | Geometric decompositions and local bases for spaces of nite element dierential forms[END_REF]

Notations and denitions

Let n P N ˚be the space dimension and d P v0, nw. We recall the denition of a d-simplex in R n as well as its barycentric coordinates. Denition 1 (d-simplex) A d-simplex K Ă R n is the convex hull of pd `1q points q q q q q q q q q 0 , . . . ,q q q q q q q q q d of R n which form the vertices of K.

K " # q q q P R n | q q q " d ÿ i"0 θ i q q q q q q q q q i , with @i P v0, dw, θ i ě 0, and

d ÿ i"0 θ i " 1 + . (2)
For example, a 2-simplex is a triangle and a 3-simplex is a tetrahedron. It will be always assumed that a d-simplex is not degenerated, i.e., the set of vectors tq q q q q q q q q i ´q q q q q q q q q 0 u d i"1 is linearly independent.

Denition 2 (Barycentric coordinates) Let K Ă R n be a non-degenerate d-simplex and tq q q q q q q q q i u d i"0 its vertices. The parametrization of K with a convex combination of the vertices reads as follows K " # q q q P R n | q q q " d ÿ i"0 λ i pq q qqq q q q q q q q q i , with @i P v0, dw, λ i pq q qq ě 0, and

d ÿ i"0 λ i pq q qq " 1 + . (3)
The coecients λ 0 pq q qq, . . . , λ d pq q qq are called the barycentric coordinates on K of q q q.

As immediat property, the barycentric coordinates on K satisfy λ i pq q q q q q q q q j q " δ i,j , @pi, jq P v0, dw.

2 Some results on the unit d-simplex

The unit d-simplex Kd Ă R d is dened by the d `1 vertices q q q q q q q q q 0 , q q q q q q q q q 1 , ¨¨¨, q q q q q q q q q d ( " 0 0 0, ê ê ê1 , ¨¨¨, ê ê êd ( where ê ê ê1 , . . . , ê ê êd ( is the standard basis of R d . We have Kd "

# q q q P R d | q q q " d ÿ i"0
λi pq q qqq q q q q q q q q i , with λi pq q qq ě 0, and

d ÿ i"0
λi pq q qq " 1

+ . (5)
As immediat property, the barycentric coordinates p λi q d i"0 on Kd satisfy λi pq q q q q q q q q j q " δ i,j , @pi, jq P v0, dw.

and are explicitly given with q q q " px 1 , ¨¨¨, x d q t P Kd by λ0 pq q qq " 1 ´d ÿ

i"1

x i and @i P v1, dw, λi pq q qq " x i .

Indeed, as q q q q q q q q q 0 " 0 0 0, we have q q q " d ÿ i"0 λi pq q qqq q q q q q q q q i " d ÿ i"1 q q q q q q q q q i λi pq q qq From q q q q q q q q q i " ê ê êi , @i P v1, dw, we obtain d ÿ i"1 q q q q q q q q q i λi pq q qq " ¨q q q q q q q q q 1 ¨¨¨q q q q q q q q q d '¨λ 1 pq q qq . . . λd pq q qq ‹ '" I d ¨λ 1 pq q qq . . . λd pq q qq ‹ '" ¨λ 1 pq q qq . . . λd pq q qq ‹ ' and thus q q q " ¨x1 . . .

x d ‹ '" ¨λ 1 pq q qq . . . λd pq q qq ‹ '.

From (5), we have

d ÿ i"0
λi pq q qq " 1 and thus λ0 pq q qq " 1 ´d ÿ

i"1

λi pq q qq " 1 ´d ÿ

i"1

x i .

unit d-simplex volume

There are several ways to compute the volume | K| of the d-simplex K Ă R d which is given by the following integral:

| K| " ż K 1dq q q " ż 1 0 ż 1´x1 0 ż 1´x1´x2 0 . . . ż 1´px1`...`x d´1 q 0 1dx d . . . dx 3 dx 2 dx 1 .
An elegant way to perform this integration is explained in [START_REF] Jaynes | Probability Theory The Logic Of Science Edition 1[END_REF], section 18.10, and uses a Laplace transform. To use this method, we note that

K " R d `X t1 ´px 1 `. . . `xd q ě 0u. (8)
So we also have

| K| " ż R d `Xt1´px1`...`x d qě0u 1dx d . . . dx 1 .
By using a dirac function and extending the integration domain to R d`1 `, we also have

| K| " ż R d`1 `δpx 1 `. . . `xd `xd`1 ´1qdx d`1 dx d . . . dx 1
To use the Laplace transform theory, we dene the function f by

f ptq " ż R d`1 `δpx 1 `. . . `xd `xd`1 ´tqdx d`1 dx d . . . dx 1 so that | K| " f p1q. The Laplace transform of f is given by Lpf qpsq " ż 8 0 f ptqe ´st dt " ż R d`1 `ˆż 8 0 δpx 1 `. . . `xd `xd`1 ´tqe ´st dt ˙dx d`1 dx d . . . dx 1 " ż R d`1 `expp´s d`1 ÿ i"1 x i qdx d`1 dx d . . . dx 1 " d`1 ź i"1 ż 8 0 expp´sx i qdx i " 1 s d`1 .
By using the inverse Laplace transform table (see [START_REF] Schi | The Laplace Transform[END_REF] for example), we have

L -1 ps Þ Ñ d! s d`1 qptq " t d .
As f " L -1 ˝Lpf q and by linearity of the inverse Laplace transform we obtain

f ptq " t d d! .
So the volume of the unit d-simplex is

| K| " 1 d! (9)
2.2

Magic formula

Let ν ν ν " pν 0 , . . . , ν d q P N d`1 . The magic formula is given by

ż K d ź i"0 λνi i pq q qqdq q q " ś d i"0 ν i ! pd `řd i"0 ν i q! (10)
This formula is often used in P 1 -Lagrange nite element methods because P 1 -Lagrange basis functions on a d-simplex are the associated barycentic coordinates. For example, one can refer to [START_REF] Quarteroni | Numerical models for dierential problems[END_REF] (section 8.2.1, page 179, formula p8.3q), [START_REF] Vermolen | On an integration rule for products of barycentric coordinates over simplexes in rn[END_REF], [START_REF] Frey | Chapter 7. the nite element approximation[END_REF] section 7.3.3 page 126, [START_REF] Eisenberg | On nite element integration in natural co-ordinates[END_REF] for d P v1, 3w, [START_REF] Brenner | The Mathematical Theory of Finite Element Methods[END_REF] as exercise for d " 2 and d " 3. In this section, we propose a proof of this formula using Laplace transform theory. Let Îpνq denote the integral of (10). The barycentic coordinates λi are given in (7) and so with q q q " px 1 , . . . , x d q and using ( 8) we obtain

Îpνq " ż K p1 ´d ÿ i"1 x i q ν0 d ź i"1 x νi i dx d . . . dx 1 " ż R d `Xt1´px1`...`x d qě0u p1 ´d ÿ i"1 x i q ν0 d ź i"1
x νi i dx d . . . dx 1

From section 2.1, by using a dirac function and by extending the integration domain to R d`1 `we obtain with

ν d`1 " ν 0 Îpνq " ż R d`1 `δpx 1 `. . . `xd `xd`1 ´1qx ν0 d`1 d ź i"1 x νi i dx d`1 dx d . . . dx 1 " ż R d`1 `δpx 1 `. . . `xd `xd`1 ´1q d`1 ź i"1 x νi i dx d`1 dx d . . . dx 1
To use the Laplace transform theory, we dene the function f ν ν ν by

f ν ν ν ptq " ż R d`1 `δpx 1 `. . . `xd `xd`1 ´tq d`1 ź i"1
x νi i dx d`1 dx d . . . dx 1 so that Îpνq " f ν ν ν p1q. The Laplace transform of f ν ν ν is given by Lpf ν ν ν qpsq "

ż 8 0 f ν ν ν ptqe ´st dt " ż R d`1 `ˆż 8 0 δpx 1 `. . . `xd `xd`1 ´tqe ´st dt ˙d`1 ź i"1 x νi i dx d`1 dx d . . . dx 1 " ż R d`1 `expp´s d`1 ÿ i"1 x i q d`1 ź i"1 x νi i dx d`1 dx d . . . dx 1 " d`1 ź i"1 ż 8 0 x νi i expp´sx i qdx i " d`1 ź i"1 Lpt Þ Ñ t νi qpsq
In a classical Laplace transform table (see [START_REF] Schi | The Laplace Transform[END_REF] for example), we have

Lpt Þ Ñ t k k! qpsq " 1 s k`1
and by linearity of the Laplace transform

Lpt Þ Ñ t k qpsq " k! s k`1 .

So we obtain

Lpf ν ν ν qpsq " d`1 ź i"1 ν i ! s νi`1 " ś d`1 i"1 ν i ! s d`1`ř d`1 i"1 νi
By using the inverse Laplace transform table, we have

L -1 ps Þ Ñ 1 s k qptq " t k´1 k ´1 .
With the linearity of the inverse Laplace transform we obtain

f ν ν ν ptq " L -1 pLpf ν ν ν qpsqqptq " ś d`1 i"1 ν i ! pd `řd`1 i"1 ν i q! t d`ř d`1 i"1 νi .
As Îpνq " f ν ν ν p1q and ν d`1 " ν 0 , the equation ( 10) is proved.

3 Some results on a d-simplex in R n

3.1

Gradients of Barycentric coordinates on a d-simplex

Lemma 3 Let K Ă R n be a non-degenerate d-simplex and and tq q q q q q q q q i u d i"0 its vertices. The barycentric coordinates pλ i pq q qqq d i"0 are solution of the linear system

¨1 1 ¨¨¨1 0 . . . A t K A K 0 ‹ ‹ ‹ ' ¨λ0 pq q qq λ 1 pq q qq . . . λ d pq q qq ‹ ‹ ‹ ‹ ' "
¨1 A K pq q q ´q q q q q q q q q 0 q ‹ ‹ '

where A K P M n,d pRq is dened by A K " ¨q q q q q q q q q 1 ´q q q q q q q q q 0 ¨¨¨q q q q q q q q q d ´q q q q q q q q q 0 ' (12)

The barycentric coordinates are multivariate polynomials of rst degree and their gradients are given by

`∇ ∇ ∇ λ 1 pq q qq ¨¨¨∇ ∇ ∇ λ d pq q qq ˘" A K pA t K A K q -1 (13)
and

∇ ∇ ∇ λ 0 pq q qq " ´d ÿ i"1 ∇ ∇ ∇ λ i pq q qq. (14)
Proof: As ř d i"0 λ i pq q qq " 1, we have q q q " d ÿ i"0 λ i pq q qqq q q q q q q q q i ùñ q q q ´q q q q q q q q q 0 " d ÿ i"1 pq q q q q q q q q i ´q q q q q q q q q 0 qλ i pq q qq " A K ¨λ1 pq q qq . . .

λ d pq q qq ‹ '
Due to linear independence of tq q q q q q q q q i ´q q q q q q q q q 0 u d i"1 ,

H K def " A t K A K P M d,d pRq (15)
is a regular matrix and the barycentric coordinates are solution of the linear system

A t K A K
¨λ1 pq q qq . . . λ d pq q qq ‹ '" A t K pq q q ´q q q q q q q q q 0 q and d ÿ i"0 λ i pq q qq " 1.

In matrix form these equations can be written as (11) and we deduce that the barycentric coordinates λ i are multivariate polynomials of rst degree. So their gradients are constants on K.

The ane map/transformation F K from the unit d-simplex K Ă R d to K Ă R n is given by q q q " A K q q q `q q q q q q q q q 0 " F K pq q qq. (16) So we have A t K pq q q ´q q q q q q q q q 0 q " A t K A K q q q " H K q q q and thus F -1 K : K Ă R n ÝÑ K Ă R d is dened by q q q " H -1 K A t K pq q q ´q q q q q q q q q 0 q " F -1 K pq q qq. (17

)
So we have λ i pq q qq " p λi ˝F-1 K qpq q qq and λi pq q qq " pλ i ˝FK qpq q qq (18)

One can remark that if d " n then A K is a regular square matrix and H -1

K A t K " A -1 K .
Now, we may compute partial derivative of λ i and @i P v0, dw, @j P v1, nw, we obtain with q q q " px 1 , . . . , xd q and q q q " px 1 , . . . , x n q Bλ i Bx j pq q qq "

d ÿ l"1 B λi Bx j pF -1 K pq q qqq BF -1 K,l Bx j pq q qq From (17), denoting B K " H -1 K A t K P M d,m pRq gives BF -1 K,l
Bxj pq q qq " pB K q l,j . The barycentric coordinates are polynomials of rst degree, so their gradients are constants and we obtain

∇ ∇ ∇ λ i " B t K ∇ ∇ ∇ λi (in fact B K is the Jacobian matrix of F -1 K ).
The matrix H K is regular and symmetric, so B t K " A K H -1 K and we obtain

∇ ∇ ∇ λ i " A K H -1 K ∇ ∇ ∇ λi . (19)
From ( 7), we deduced that

´∇ ∇ ∇ λ1 ¨¨¨∇ ∇ ∇ λd ¯" I d (20)
and thus

`∇ ∇ ∇ λ 1 pq q qq ¨¨¨∇ ∇ ∇ λ d pq q qq ˘" A K H -1 K ´∇ ∇ ∇ λ1 ¨¨¨∇ ∇ ∇ λd " A K H -1 K .
As ř d i"0 λ i pq q qq " 1, we immediately have

∇ ∇ ∇ λ 0 pq q qq " ´d ÿ i"1 ∇ ∇ ∇ λ i pq q qq.
From (13) and ( 14), we immediatly have: Remark 4 The gradients of the barycentric coordinates are linear combinations of tq q q q q q q q q 1 ´q q q q q q q q q 0 , . . . ,q q q q q q q q q d ´q q q q q q q q q 0 u.

3.2

Integration over a d-simplex

If K is a non-degenerated d-simplex in R d , from (16) we have J F K pq q qq " A K .
Then A K is a regular square matrix and we have the classical formula:

ż K f pq q qqdq q q " |detpA K q| ż K f ˝FK pq q qqdq q q (21)
The following theorem extend this result to d-simplex in R n , with 1 ď d ď n.

Theorem 5 Let K Ă R n be a non-degenerated d-simplex and f : K ÝÑ R.

ż K f pq q qqdq q q " ˇˇdetpA t K A K q ˇˇ1{2 ż K f ˝FK pq q qqdq q q (22)
where K is the unit d-simplex in R n , A K P M d,n pRq is dened by A K " ´q q q q q q q q q 1 ´q q q q q q q q q 0 q q q q q q q q q 2 ´q q q q q q q q q 0 ¨¨¨q q q q q q q q q d ´q q q q q q q q q 0 ¯(23)

and F K : K ÝÑ K is given by F K pq q qq " A K q q q `q q q q q q q q q 0 (24)

Proof: The set v v v 1 , . . . , v v v d (
is linearly independent so we can extend it to a basis v v v 1 , . . . , v v v n ( . We denote by A P M n,n pRq the matrix such that the i-th column is the vector v v v i for all i P v1, nw. So we have

A " ´AK v v v d`1 ¨¨¨v v v n ¯(25)
By the QR-factorization theorem apply to the matrix A P M n pRq, there is an orthogonal matrix Q P M n pRq and a regular upper triangular matrix R P M n pRq such that A " QR So we have

Q t A " R
and we dene the matrix Ā P M n,d pRq to be the rst d columns of R:

Ā " Q t A K .
We can also note that Ā " ´q q q q q q q q q 1 ´q q q q q q q q q 0 q q q q q q q q q 2 ´q q q q q q q q q 0 ¨¨¨q q q q q q q q q d ´q q q q q q q q q 0 ¯" ´q q q q q q q q q 1 q q q q q q q q q 2 ¨¨¨q q q q q q q q q d ¯.

Let F : R n ÝÑ R n be the bijective function dened by Fpx x xq " Q t px x x ´q q q q q q q q q 0 q " x x x (26) and q q q q q q q q q i " Fpq q q q q q q q q i q " Q t pq q q q q q q q q i ´q q q q q q q q q 0 q, @i P v0, dw.

By construction q q q q q q q q q 0 " 0 0 0 and, @i P v1, dw, q q q q q q q q q i is the i-th column of the upper triangular matrix R. So we have @i P v0, dw, q q q q q q q q q i P Vectpe e e 1 , . . . , e e e d q

where e e e 1 , . . . , e e e n ( is the standard basis of R n . The set q q q q q q q q q 0 , . . . , q q q q q q q q q d ( are the vertices of the d-simplex K " FpKq and we deduce K Ă Vectpe e e 1 , . . . , e e e d q.

(27)

By change of variables, we obtain ż K f pq q qqdq q q " ż K f ˝F -1 pq q qq| detpJ F -1 pq q qqq|dq q q where J F -1 is the Jacobian matrix of F-1 . From (26), we have J F -1 pq q qq " Q and as Q is an orthogonal matrix, detpJ F -1 pq q qqq " 1. So we obtain ż K f pq q qqdq q q " ż K f ˝F -1 pq q qqdq q q.

(28)

Let P P M d,n pRq dened by P " `Id O d,n´d ȃnd @i P v0, dw, q q q q q q q q q i " Pq q q q q q q q q i P R d .

From (27), we deduce @i P v0, dw, q q q q q q q q q i " P t q q q q q q q q q i " ˆq q q q q q q q q i 0 0 0 ˙.

Let ḡ " f ˝F -1 and K be the d-simplex in R d with vertices q q q q q q q q q i , i P v0, dw. We denote by P : K Ă R d ÝÑ K Ă R n the application dened by Pp q q qq " P t q q q. We denote by ḡ : K ÝÑ R the application dened by ḡp q q qq " ḡ ˝Pp q q qq " ḡˆq q q 0 0 0 ˙.

3.3

Volume of a d-simplex

The volume/measure of the d-simplex K Ă R n is given by |K| " ż K 1dq q q (33)

Using formula (22) with f " 1 gives

|K| " |detpA t K A K q| 1{2 ż K 1dq q q " |detpA t K A K q| 1{2 | K|.
From ( 9), we nally obtain

|K| " |detpA t K A K q| 1{2 d! . (34)
In [START_REF] Hanson | Geometry for n-Dimensional Graphics[END_REF] this formula is proved with geometrical arguments. We can also remark that if d " n then A K is a square matrix and we obtain the classical formula

|K| " |detpA K q| d! . ( 35
)

Magic formula

In this section an exact computation of the integral over a d-simplex K Ă R n for products of power of its barycentric coordinates given by ( 1) is proved by using previous results obtained by Laplace transforms.

Using formula (22) with f pq q qq " ś d i"0 λ νi i pq q qq gives ż K d ź i"0 λ νi i pq q qqdq q q " |detpA t K A K q| 1{2 ż K d ź i"0 pλ i ˝FK pq q qqq νi dq q q

From (18) and (34), we obtain

ż K d ź i"0
λ νi i pq q qqdq q q " d!|K| ż K d ź i"0 λνi i pq q qqdq q q

Using formula (10) gives [START_REF] Arnold | Geometric decompositions and local bases for spaces of nite element dierential forms[END_REF]

This work was supported in part by the ANR project DEDALES under grant ANR-14-CE23-0005.

So we obtain ż K ḡpq q qqdq q q " ż K ḡp q q qqd q q q (29)

Let Ā P M d pRq be the matrix dened by Ā " ´q q q q q q q q q 1 ´q q q q q q q q q 0 q q q q q q q q q 2 ´q q q q q q q q q 0 ¨¨¨q q q q q q q q q d ´q q q q q q q q q 0 ¯.

We can remark that Ā " P Ā and Ā " P t Ā.

Let F : K ÝÑ K the bijective function dened by Fpq q qq " Āq q q `q q q q q q q q q 0

We can now apply the classical change of variables ż K ḡp q q qqd q q q " ż K ḡ ˝F pq q qq| detpJ F pq q qqq|dq q q "| detp Āq| ż K ḡ ˝F pq q qqdq q q To resume from ( 22) and (29), we have ż K f pq q qqdq q q " | detp Āq| ż K ḡ ˝F pq q qqdq q q (31)

We can note that ḡ ˝F " f ˝F -1 ˝P ˝F Let F K " F-1 ˝P ˝F , we have as expected F K pq q qq " F-1 ˝P ˝F pq q qq " F-1 pP t p Āq q qqq " F-1 p Āq q qq " Q Āq q q `q q q q q q q q q 0 " A K q q q `q q q q q q q q q 0 . and we obtain ż K f pq q qqdq q q " | detp Āq| ż K f ˝FK pq q qqdq q q (32)

To obtain formula (22), it remains to prove that | detp Āq| " | detpA t K A K q| 1{2 . We have

As Ā is a square matrix, we have detp Āt Āq " detp Āq 2 and thus