Cyrille Bertelle

Antoine Dutot

Michel Nabaa

Damien Olivier

Pascal Mallet
email: pascal.mallet@agglo-havraise.fr

MODELING OF THE VULNERABILITY RELATED TO THE DYNAMIC ROAD TRAFFIC

Keywords: Dynamic Vulnerability Map, Decision Support System, Risk, Dynamic Graph, Communities Detection, Self-Organization, GIS, Traffic Flow

The utilization of the road network by vehicles with different behaviors can generate a danger under normal and especially under evacuation situations. In Le Havre agglomeration (CODAH), there are 33 establishments classified SEVESO with high threshold. The modeling and assessment of the danger is useful when it intersects with the exposed stakes. The most important factor is people. In the literature, vulnerability maps are constructed to help decision makers assess the risk. These maps are based on several types of vulnerability: socio-demographic, biophysical and other different types of hazards. Nevertheless, such approaches remain static and do not take into account the population displacement in the estimation of the vulnerability. We propose a decision support system which consists in a dynamic vulnerability map based on the difficulty to evacuate different sectors in Le Havre agglomeration. This map is visualized using the Geographic Information System (GIS) of the CO-DAH and evolves according to the dynamic state of the road traffic through a detection of communities in a large graph. This detection is realized by an ant algorithm.

INTRODUCTION

Le Havre agglomeration is exposed to several types of natural and industrial hazards: 16 establishments are classified SEVESO 1 (EUROPEAN Community 1996) with very high threshold. The examination of impacted populations remains a difficult exercise. In this context, the Major Risk Management Direction team (DIRM) of Le Havre Agglomeration (CODAH) has developed a model which estimates the nocturnal and diurnal exposed population allocation PRET-RESSE [START_REF] Bourcier | Allocation spatio-temporelle de la population exposée aux risques majeurs. Contribution à l'expologie sur le bassin de risques majeurs de l'estuaire de la seine: modèle PRET-RESSE[END_REF]; the scale is the building. Although the model is able to locate the diurnal and nocturnal population, it remains static because it does not take into account the daily movement of people and the road network utilization. For a better people evacuation in a major risk case, we need to have detailed information about the state of road traffic network to determine how to allocate the vehicles on the road network and model the movement of these vehicles. In fact, the panic effect of some people and the redundant drivers behaviors can lead to accidents and traffic jam, this can be very grievous and spread quickly. In the literature, several models were developed to calculate a score for the vulnerability related to the road network utilization. Most of these models adopt a pessimistic approach to calculate this vulnerability: this case is met when a group of people in a hazardous area decide all to take the same route to evacuate this area, which unfortunately happens quite often in the real world evacuation situations. Although it helps decision-makers to estimate the risk by a census vulnerability map, this approach remains static and does not take into account the evolution of the road network traffic. In this paper, we propose a dynamic and pessimistic approach related to the road network utilization. To this end, we model the road network by a dynamic graph (the dynamics is due to the traffic evolution). A simple model based on traffic flow will also be proposed and the interaction between micro and macro traffic simulation will be discussed. Then, we apply a self-organization algorithm (belonging to the collective intelligence algorithms) to detect communities in the graph. These communities correspond to clusters in the graph which evolve along the time. Each community is the Community. a group of people in a neighborhood who share totally or partially, at a give time, the same route. The algorithm allows us to define the different vulnerable neighborhoods of the agglomeration in the case of a pessimistic evacuation due to a potential danger, while taking into account the evolution of the road network traffic due to the impact of the other flows coming from other neighborhoods. The result of this algorithm will be visualized into a GIS on a dynamic vulnerability map which categorizes various sectors depending on the difficulty of access to the road network. Finally, this vulnerability map will construct a tool to aid decision-makers in a better estimation of risk in the communes of the CODAH. This tool enriches the PRET-RESS static model developed at the CODAH, taking into account the mobility of the population.

STATE OF THE ART

Traditional methods evaluating the risk for population do not treat generally the behavioral of evacuee (e.g. initial response to an evacuation, travel speed, family interactions / group, and so on.); they describe prescriptive rules as the travel distance. These traditional methods are not very sensitive to human behavior for different emergency scenarios. The computerized models offer the potential to evaluate the evacuation of a neighborhood in emergency situations and overcome these limitations [START_REF] Castel | Developing a prototype agent-based pedestrian evacuation model to explore the evacuation of King's Cross St Pancras underground station. Center for advanced spatial analysis[END_REF]. Recently, some interesting applications have been developed by including the population dynamics, the models of urban growth patterns and land use. For computer modelers, this integration provides the ability to have computing entities as agents that are linked to real geographical locations. For GIS users, it provides the ability to model the emergence of phenomena by various interactions of agents in time and space by using a GIS [START_REF] Najlis | Repast for GIS[END_REF]. So, combining several layers as houses, road network, population... allows us to model different types of agents into a GIS environment. In [START_REF] Cutter | Revealing the vulnerability of people and places: a case study of Georgetown County, South Carolina[END_REF], the authors present a method to spatially estimate the vulnerability and treats the biophysical and social aspects (access to resources, people with evacuation special needs, people with reduced mobility ...). Several layers are created in the GIS (a layer by a danger), and all these layers are combined into one composed of intersecting polygons to build a generic vulnerability map.

To complete this, it was necessary to take into account the infrastructure and various possible routes of evacuation. So, a new map has been constructed and a new layer has been incorporated. This work has been applied to the George Town canton in which we find various natural and industrial risks, and where there are different types of people. In the neighborhood evacuation cases on a micro scale, a number of studies based on micro simulation have been developed. In their paper [START_REF] Church | Mapping evacuation risk on transposion networks using a spatial optimization model[END_REF], the authors presented a model to estimate the necessary time to evacuate a neighborhood according to the effective of the population, the number of vehicles, roads capacity and the number of vehicles per minute. The model is based on the optimization in order to find the critical area around a point at a potential danger in a pessimistic way. This model has been coupled with a GIS (ArcInfo) to visualize the results (identify evacuation plans) and construct an evacuation vulnerability map for the city (Santa Barbara). Cova and Church [START_REF] Church | Modelling community evacuation vulnerability using GIS[END_REF] opened the way on the study based on geographic information systems to evacuate people. Their study identified the communities that may face transport difficulties during an evacuation. Research has modelled the population by lane occupation during an evacuation emergency using the city of Santa Barbara. An optimization based model (graph partitioning problem) was realized to find the neighborhood that causes the highest vulnerability around each node in the graph and a vulnerability map around nodes in the city was constructed.

A constructive heuristic has been used to calculate the best cluster around each node. This heuristic was developed in C and the result was displayed on a map (with ArcInfo). Nevertheless, in this approach, we predefine the maximum number of nodes in a neighborhood, which may not always be realistic and does not take into account the traffic evolution during the calculation of critical neighborhoods. So, the vulnerable neighborhoods don't evolve according to traffic state. Since an individual panic under evacuation situation may cause a collective panic and the changes that may occur in the environment (buildings collapse, reverse route direction), the evolution of the system is unpredictable; so, in our work, we try to build a dynamic vulnerability map evolving with the traffic dynamics, in which the nodes number of in a critical neighborhood, is not predetermined and can change depending on actual traffic state.

DYNAMIC MODEL

Problem description

In this paper, the term "vulnerability" depends on the access to the road network. To address this vulnerability, we have to finely represent the population and the dynamic state of road traffic. In PRET-RESSE model developed within the major risks management team of CODAH, we have ventilated the day / night population at the building scale. The model was able to locate people during the day both in their workplace and their residence (the unemployed and retirees). It has been estimated that people will be in their residence during the night. PRET-RESS will be enriched by our model that will try to dynamically assess the vulnerability related to the road traffic evolution; so we are interested in the allocation of the vehicles on the road network, the regulation of the traffic flow and trying to dynamically estimate the vulnerability according to the traffic.

System architecture

Our system consists of two modules as shown in the figure 1.

Figure 1: System architecture The simulation module contains three components:

• The dynamic graph extracted from the road network layer and detailed in the following section,

• The flow management component consists of vehicles flow simulator applied on the graph,

• The communities detection component, detailed in the "Communities Detection" section. Its input is the extracted graph and the current flow. It returns the communities that are formed according to the current state of road traffic.

The visualization module consists of the road network layer integrated into the GIS. This module communicates with the simulation module: the graph is constructed from this module, which in turn get the simulation result and visualize it as a dynamic vulnerability map.

Environment modeling

The road network is integrated as a layer in the Geographic Information System (GIS). From this layer, we extract the data by using the open source java GIS toolkit Geotools. This toolkit provides several methods to manipulate geospatial data and implements Open Geospatial Consortium (OGC) specifications, so we can read and write to ESRI shapefile format. Once data road network are extracted, we use the GraphStream tool [START_REF] Dutot | Graph-Stream: A Tool for bridging the gap between Complex Systems and Dynamic Graphs[END_REF]) developed within LITIS laboratory of Le Havre to construct a graph corresponding to the GIS network layer. This tool is designed for modeling; processing and visualizing graphs.

Figure 2 represents the network layer into the GIS and the figure 3 is the graph extracted with geotools and visualized with GraphStream by using the same points coordinates.

The data extracted from network layer contains the roads circulation direction, roads id, roads type, their lengths and geometry.

The extracted multigraph G (t) = (V (t), E (t)) represents the road network at time t where V (t) is the set of nodes at t and E (t) the set of arcs at t. We deal with a multigraph because we have sometimes more than one oriented arc in the same direction between two adjacent nodes due to multiple routes between two points in Le Havre road network. GraphStream facilitates this task because it is adapted to model and visualize multigraphs. In the constructed multigraph:

• The nodes represent roads intersections,

• The arcs represent the roads taken by vehicles,

• The weight on each arc represents the needed time to cross this arc, depending on the current load of the traffic

• Dynamic aspect relates to the weights of the arcs, which can evolve in time, according to the evolution of the fluidity of circulation.

We have also constructed a Voronoi tessellation (Thiessen polygon) around nodes and projected the population in buildings on these nodes. The population in buildings is extracted from PRET-RESS model.

Vehicles flow

For a better people evacuation in a major risk situation, we need to know the state of the road network to determine how to allocate the vehicles on this network and to model the movement of these vehicles. Different types of models can be adopted:

• The microscopic model details the behavior of each individual vehicle by representing interactions with other vehicles and in general by using a spatialization. It is used on the scale of a sector or a neighborhood. It has the advantage to model vehicle behavior in an evacuation of a neighborhood, people panic, interactions between vehicles, accidents. . . So, if we have vehicles with different behaviors and interacting between them, these vehicles sometimes self organize and with a top down approach, we can examine the global behavior of the system and try to locally modify the environment, when necessary, to ameliorate the system,

• The macroscopic model is based on the analogy between vehicular traffic and the fluid flow within a canal. It allows us to visualize the flow on the roads rather than individual vehicles. It is used at many sectors or the entire city scale,

• The hybrid model allows coupling the two types of dynamics flow models within the same simulation. Several works have already borrowed this direction [START_REF] Hennecke | Macroscopic simulation of open systems and micro-macro link[END_REF][START_REF] Bourrel | Mixing micro and macro representation of trafic flow : a first theoretical step[END_REF][START_REF] Magne | Toward an hybrid macro and micro traffic flow simulation model[END_REF], however, this approach is relatively new and very few have adopted it to our knowledge [START_REF] Hman | Simulation hybride de flux de trafic basée sur les systèmes multi-agents[END_REF]. The use of a hybrid model is very important to us: changing the scale from micro to macro in a region where we haven't a crisis situation (everything is normal) allows to economize the computation and the In this paper, we used a simple model of macroscopic flow:

• A flow of cars moving from one arc to another adjacent one,

• The arcs are limited in capacity,

• The flow can be broken and two or more flows can gather on a node,

• Traffic jams may appear in certain places of the road network; those places will be more vulnerable than others.

We have adopted a macroscopic model in which flows circulate normally (no accidents) because the goal now is to establish a dynamic pessimistic vulnerability map which is not always the case in the real world (e.g. 90% of people takes an exit route and the rest takes another route for example). Hence, it is important to have in the near future a micro approach with a change of scale (from micro to macro and vice versa during the simulation) to simulate scenarios of danger in real time (accidents, behavior of drivers, vehicles interactions ...), a study on which we are working actually.

COMMUNITIES DETECTION

Our aim is to identify communities in graphs, i.e. dense areas strongly linked to each other and more weakly linked to the outside world. If the concept of communities in a graph is difficult to define formally, it can be seen as a set of nodes whose internal connections density is higher than the outside density without defining formal threshold [START_REF] Pons | Détection de structures de communautés dans les grands réseaux d'interactions[END_REF]. Thus the goal is to find a partition of nodes in communities according to a certain predefined criteria without fixing the number of such communities or the number of nodes in a community. Interesting works were developed in the literature on the detection of structure in large communities in graphs [START_REF] Clauset | Finding community structure in very large networks[END_REF], Newman 2004a;b, Pons and Latapy 2006).

In our problem, we look for a self-organization in networks with an algorithm close to the detection of communities in large graphs and belonging to collective intelligence algorithms. Organizations connect elements, events or individuals by interrelations so that they become components of a whole. They assume the solidarity and robustness of these links, and ensure that the system will eventually be long lasting despite random perturbations. The organizations, therefore : transform, produce, tie and maintain. Time is present generating the dynamic and we will try to fight against these organizations, in the case of risk management, to avoid bottlenecks which do not facilitate the evacuations. Thus we will considered the graph at time t, G(t) = (V (t), E(t)) where the edges are weighted, this weight being noted |e| for the edge e and represents the needed time The algorithm used can be referred to as an ant algorithm [START_REF] Dorigo | Ant Colony Optimization[END_REF]. Our ant algorithm use several colonies of ants, each of a distinct color. Ants travel inside the graph an lay down pheromones, information that can be detected by other ants. Pheromones are also colored. Ants tend to be repulsed by pheromones of other colors. Furthermore, ants tend to favor edges with important weights. The colored dynamic graph precedently mentionned is defined such that:

• V (t) is the set of vertices at time t. Each vertex v is characterized by:

a color c ∈ C(t),

• E(t) is the set of edges at time t. Each edge e is characterized by:

a weight |e| ∈ N + that corresponds to interaction importance between the elements at each end of edge e.

a quantity of pheromones of each color.

• C(t) is a set of colors representing the ant colonies at time t.

The algorithm principle is to color the graph using pheromones. Each colony will collaborate to colonize zones, whereas colonies compete to maintain their own colored zone (see figure 4). Solutions will therefore emerge and be maintained by the ant behavior (see figure 5). The solutions will be the color of each vertex in the graph. Indeed, colored pheromones are deposited by ants on edges. The color of a vertex is obtained from the color having the largest proportion of pheromones on all incident edges (see algorithm 1).

We have an interaction between each two local adjacent nodes according to the attraction force that exists between them. This force depends in our case on the report of the number of vehicles on the arc between 2 nodes neighbors / vehicles capacity of the arc. This report was chosen because, in every community, we will have a large number of vehicles which all decide to exit through a single road in the case of a potential danger; this responds well to one of the purposes listed in beginning to have a pessimistic approach in the calculation of vulnerability. The algorithm has the advantage of not allowing the breaking of a link between 2 adjacent nodes to maintain the structure of the road network. When the traffic evolves, the algorithm detects that and communities can change or disappear as a result of local forces that change between the nodes locally.

Dynamic vulnerability map

At each simulation time step, the flow on the arcs changes following traffic conditions and the attraction may change also. Once communities are formed on the graph, the result will be transmitted to the road network layer into the GIS to Choose an edge to cross in a weighted random fashion, using edges weight (if available) Lay down a small amount of pheromone of the ant color on this edge n ← vertex at the other end of the chosen edge be visualized.

CONCLUSION

In this paper, we have proposed a decision support system to assess the danger depending on the road network usage by the vehicles population using the network. This tool enables decision makers to visualize, on a geographic information system, a dynamic vulnerability map related to the difficulty of evacuating the various streets in the metropolitan area of Le Havre agglomeration. We simulated road network traffic by using a simple model of vehicles flow. A communities detection algorithm in the large graphs was adopted. It enabled us to form communities in a graph thanks to local forces propagation rule between adjacent nodes. The communities evolve according to the current state of the road network traffic. The result of the evolution of communities is visualized by using a GIS. The adopted approach allowed us to estimate the risk due to the use of the road network by vehicles and categorize Le Havre agglomeration areas by their vulnerability. We will complete our work by using real traffic data retrieved from a displacements survey with Le Havre population which will help us to better locate people during the day and therefore having a more realistic vulnerability dynamic map.

Figure 2 :

 2 Figure 2: The initial network layer

Figure 3 :

 3 Figure 3: The extracted graph

Figure 5 :

 5 Figure 5: Example of a dynamic evolution and communities detection.

Figure 4 :

 4 Figure 4: Communities detection in a graph


Claire David
email: claire.david@sorbonne-universite.fr

Michel L Lapidus
email: lapidus@math.ucr.edu

New Insights for Fractal Zeta Functions Polyhedral Neighborhoods vs Tubular Neighborhoods

Keywords:

, we introduce new fractal zeta functions, associated with polyhedral neighborhoods, better suited to fractals when exact expressions for the volume of tubular neighborhoods cannot be computed. Accordingly, in the model -and very significant -case of the Weierstrass Curve Γ W , we give exact expressions for the volume of polyhedral neighborhoods for the sequence of prefractal graphs which converge to Γ W -the so-called Weierstrass Iterated Fractal Drums (in short, Weierstrass IFDs). Those IFDs are associated with a suitable (and geometrically meaningful) sequence of small parameters tending to zero, also known as the cohomology infinitesimals, due to their connections with fractal cohomology. We also introduce the associated local and global polyhedral fractal zeta functions. The local fractal zeta functions consist in the sequence of zeta functions associated with the sequence of polyhedral neighborhoods, and satisfy a recurrence relation, which enables us to prove that the poles of the limit fractal zeta function -the global zeta function, associated with the limit fractal object -are exactly the same as the Complex Dimensions of the Weierstrass function itself. This result makes the connection with fractal cohomology, where, for any nonnegative integer m, the m t h cohomology group is comprised of continuous functions which possess a generalized Taylor expansion, with fractional derivatives of orders the underlying -and actual -Complex Dimensions. By using the aforementioned exact expressions of the polyhedral neighborhoods, we also revisit the computation of the box-counting (or Minkowski) dimension of the Weierstrass Curve, in a fully rigorous manner and therefore prove part of Mandelbrot's conjecture concerning the fractal dimension of the Weierstrass Curve.

Introduction

Up to now, tubular neighborhhoods were a compulsory and unavoidable step when it came to computing fractal zeta functions (see [START_REF] Lapidus | Fractal Geometry, Complex Dimensions and Zeta Functions: Geometry and Spectra of Fractal Strings[END_REF], [START_REF]Fractal Zeta Functions and Fractal Drums: Higher-Dimensional Theory of Complex Dimensions[END_REF], [START_REF]Fractal tube formulas for compact sets and relative fractal drums: Oscillations, complex dimensions and fractality[END_REF]). However, one cannot always obtain the exact expression of the associated volumes. For instance, this occurs when the fractals involved are defined by means of nonlinear and noncontractive iterated function systems (i.f.s.). A model -and very significant -object in this case is the (nowhere differentiable) Weierstrass Curve Γ W , where, in addition, nonlinearity makes the geometry extremely complicated.

A careful look at the sequence of prefractal, polygonal graph approximations which converge to Γ W shows that not only polyhedral neighborhoods appear as perfectly suited to this fractal curve, but, even more important, we can obtain the exact expressions of the associated volumes. This comes from the fact that we know the coordinates of the vertex points which constitute the prefractal approximations. More precisely, those coordinates involve values of the Weierstrass function W.

A very interesting -and rather a priori unexpected -feature then arises. Indeed, the Weierstrass function W itself possess a Complex Dimensions series expansion. Recall that the theory of Complex Dimensions, developed for many years now by Michel L. Lapidus and his collaborators in [START_REF] Lapidus | Fractal drum, inverse spectral problems for elliptic operators and a partial resolution of the Weyl-Berry conjecture[END_REF], [START_REF]Vibrations of fractal drums, the Riemann hypothesis, waves in fractal media and the Weyl-Berry conjecture[END_REF], [START_REF] Lapidus | The Riemann zeta-function and the one-dimensional Weyl-Berry conjecture for fractal drums[END_REF], [START_REF] Lapidus | The Riemann hypothesis and inverse spectral problems for fractal strings[END_REF], [START_REF] Lapidus | A tube formula for the Koch snowflake curve, with applications to complex dimensions[END_REF], [START_REF] Lapidus | Search of the Riemann Zeros: Strings, Fractal Membranes and Noncommutative Spacetimes[END_REF], [START_REF] Lapidus | Fractal Geometry, Complex Dimensions and Zeta Functions: Geometry and Spectra of Fractal Strings[END_REF], [START_REF] Lapidus | Distance and tube zeta functions of fractals and arbitrary compact sets[END_REF], [START_REF]Fractal Zeta Functions and Fractal Drums: Higher-Dimensional Theory of Complex Dimensions[END_REF], [START_REF]An overview of complex fractal dimensions: From fractal strings to fractal drums, and back[END_REF], [START_REF] Herichi | Quantized Number Theory, Fractal Strings and the Riemann Hypothesis: From Spectral Operators to Phase Transitions and Universality[END_REF], [START_REF]From Complex Fractal Dimensions and Quantized Number Theory To Fractal Cohomology: A Tale of Oscillations, Unreality and Fractality[END_REF], makes the connection between the geometry of an object and its differentiability properties, by means of geometric (or fractal) zeta functions, which stand for the trace of a differential operator at a complex order s. The poles of those fractal zeta functions are called the Fractal Complex Dimensions. The existence of nonreal Complex Dimensions is a characteristic of fractality and gives rise (via explicit formulas) to the oscillations that are intrinsic to fractal geometries.

By considering the fractal zeta functions associated with the sequence of polygonal neighborhoods, we came across the recurrence relation between consecutive fractal zeta functions. To our knowledge, this is the first time that such a result is obtained. It is all the more important that it enables us to prove the existence of the limit fractal zeta function -the global zeta function, associated with the limit fractal object, the Weierstrass Curve here, or rather, the Weierstrass iterated fractal drum (IFD). Moreover, this global zeta function has the same (possible) poles as the Complex Dimensions involved in the aforementioned series expansion of the Weierstrass function.

Our main results in the present setting can be found in the following places: i. In Theorem 20, on page 11, where we obtain recurrence relations satisfied by the coefficients of the exact fractal (or Complex Dimensions) expansion for the complexified and the ordinary Weierstrass functions.

ii. In Definition 23, on page 12, where we introduce the sequence of polyhedral neighborhoods.

iii. In Theorem 24, on page 12, where we show that the polyhedral and tubular neighborhoods are nested.

iv. In Theorem 28, on page 16, where we obtain an exact fractal power series expansion for the volume of suitable m t h polyhedral neighborhoods of the Weierstrass Curve.

v. In Theorem 32, on page 20, where we introduce the local and global effective polyhedral zeta function and show that the global zeta function is well defined, meromorphic in all of C, and is given by an explicit fractal power series.

vi. In Theorem 33, on page 24, where we obtain the exact Complex Dimensions of the Weierstrass Curve (or of the Weierstrass IFD), precisely of the same form as the possible Complex Dimensions occurring in the fractal cohomology theory developed in [START_REF]Weierstrass fractal drums -II -Towards a fractal cohomology[END_REF].

vii. In Theorem 37, on page 26, where we introduce a new and completely rigorous proof of the computation of the box-counting dimension (or, equivalently, of the Minkowski dimension) of the Weierstrass Curve Γ W , by simply using the covers of Γ W by polyhedral neighborhoods -therefore fully establishing part of Mandelbrot's conjecture [START_REF] Mandelbrot | The Fractal Geometry of Nature[END_REF] about the fractal dimensions of Γ W .

Geometric Framework

Henceforth, we place ourselves in the Euclidean plane, equipped with a direct orthonormal frame. The usual Cartesian coordinates are denoted by (x, y). The horizontal and vertical axes are respectively referred to as (x ′ x) and (y ′ y). We also let N = {0, 1, 2, . . .} and

N ⋆ = {1, 2, . . .}. Given a, b with a, b ∈ R ∪ {-∞,
< λ < 1 , N b ∈ N ⋆ and λ N b > 1 . (♣)
Note that this implies that N b > 1 (i.e., N b ≥ 2).

As explained in [START_REF]On fractal properties of Weierstrass-type functions[END_REF], we deliberately made the choice to introduce the notation N b which replaces the initial number b, in so far as, in Hardy's paper [START_REF] Hardy | Weierstrass's Non-Differentiable Function[END_REF] (in contrast to Weierstrass's original article [START_REF] Weierstrass | Über continuirliche Funktionen eines reellen Arguments, die für keinen Werth des letzteren einen bestimmten Differential quotienten besitzen[END_REF]), b is any positive real number satisfying λ b > 1, whereas we deal here with the specific case of a natural integer, which accounts for the natural notation N b . Definition 2 (Weierstrass Function, Weierstrass Curve). We consider the Weierstrass function W (also called, in short, the W-function) defined, for any real number x, by

W(x) = ∞ n=0 λ n cos 2 π N n b x . ( 1
)
We call the associated graph the Weierstrass Curve, and denote it by Γ W .

Due to the one-periodicity of the W-function, since N b is an integer, from now on, and without loss of generality, we restrict our study to the interval [0, 1[= [0, 1). We also identify the points (0, W(0)) and (1, W(1)) = (1, W(0)).

Definition 3 (Weierstrass Complexified Function).

We introduce the Weierstrass Complexified function W comp , defined, for any real number x, by

W comp (x) = ∞ n=0 λ n e 2 i π N n b x .
Notation 4 (Logarithm). Given y > 0, ln y denotes the natural logarithm of y, while, given a > 1, ln a y = ln y ln a denotes the logarithm of y in base a; so that, in particular, ln = ln e .

Notation 5 (Minkowski Dimension and Hölder Exponent).

For the parameters λ and N b satisfying condition (♣) (see Notation 1, on page 3), we denote by

D W = 2 + ln λ ln N b = 2 -ln N b 1 λ ∈ ]1, 2[
the box-counting dimension (or Minkowski dimension) of the Weierstrass Curve Γ W , which happens to be equal to its Hausdorff dimension [START_REF] Kaplan | The Lyapunov dimension of a nowhere differentiable attracting torus[END_REF], [START_REF] Bara Ńsky | On the dimension of the graph of the classical Weierstrass function[END_REF], [START_REF] Shen | Hausdorff dimension of the graphs of the classical Weierstrass functions[END_REF], [START_REF] Keller | A simpler proof for the dimension of the graph of the classical Weierstrass function[END_REF]. We point out that the results in our present paper -and its long version [START_REF]Polyhedral neighborhoods vs tubular neighborhoods: New insights for fractal zeta functions and complex dimensions[END_REF] -also provide a direct geometric proof of the fact that D W , the Minkowski dimension (or box-counting dimension) of Γ W , exists and takes the above values, and that our results in our previous paper [START_REF]Weierstrass fractal drums -I -A glimpse of complex dimensions[END_REF] show that W is Hölder continuous with optimal Hölder exponent 2

-D W = - ln λ ln N b = ln N b 1 λ .
Note that the latter result was also obtained by G. H. Hardy in [START_REF] Hardy | Weierstrass's Non-Differentiable Function[END_REF], although new key information is obtained in [START_REF]Weierstrass fractal drums -I -A glimpse of complex dimensions[END_REF] where we establish a new discrete analog of the reverse Hölder estimates (of order D W ) which is used in an essential way in the present work; see Proposition 16, on page 8 and Theorem 17, on page 9.

Proposition 6 (Nonlinear and Noncontractive Iterated Function System (IFS)). Following our

previous work [START_REF] David | Bypassing dynamical systems: A simple way to get the box-counting dimension of the graph of the Weierstrass function[END_REF], we approximate the restriction Γ W to [0, 1[×R, of the Weierstrass Curve, by a sequence of graphs, built via an iterative process. For this purpose, we use the nonlinear (and noncontractive) iterated function system (IFS) of the family of C

∞ maps from R 2 to R 2 denoted by

T W = T 0 , ⋯, T N b -1 ,
where, for any integer i belonging to {0, ⋯, N b -1} and any point

(x, y) of R 2 , T i (x, y) = x + i N b , λ y + cos 2 π x + i N b .

Theorem 7 (Attractor of the IFS [2], [3]

). The Weierstrass Curve Γ W is the attractor of the IFS T W , and hence, is the unique nonempty compact subset

K of R 2 satisfying K = N b -1 ⋃ i =0 T i (K);
in particular, we have that

Γ W = N b -1 ⋃ i =0 T i (Γ W ).
Notation 8 (Fixed Points). For any integer i belonging to {0, ⋯, N b -1}, we denote by

P i = (x i , y i ) = i N b -1 , 1 1 -λ cos 2 π i N b -1
the unique fixed point of the map T i ; see [START_REF]On fractal properties of Weierstrass-type functions[END_REF].

Definition 9 (Sets of Vertices, Prefractals). We denote by V 0 the ordered set (according to increasing abscissae) of the points

P 0 , ⋯, P N b -1 .
The set of points V 0 -where, for any integer i in {0, ⋯, N b -2}, the point P i is linked to the point P i +1 -constitutes an oriented finite graph, ordered according to increasing abscissae, which we will denote by Γ W 0 . Then, V 0 is called the set of vertices of the graph Γ W 0 .

For any positive integer m, i.e., for m ∈ N

⋆ , we set V m = N b -1 ⋃ i =0 T i (V m-1 ).
The set of points V m , where two consecutive points are linked by an edge, is an oriented finite graph, ordered according to increasing abscissa, called the m t h order W-prefractal. Then, V m is called the set of vertices of the prefractal Γ W m ; see Figure 2, on page 7. We call Weierstrass Iterated Fractal Drums (IFD) the sequence of prefractal graphs which converge to the Weierstrass Curve. Proposition 12. [START_REF] David | Bypassing dynamical systems: A simple way to get the box-counting dimension of the graph of the Weierstrass function[END_REF] For any m ∈ N, the following statements hold:

i . V m ⊂ V m+1 . i i . #V m = (N b -1) N m b + 1, where #V m
(N b -1) k (N b -1) N m b , W (N b -1) k (N b -1) N m b , 1 ≤ k ≤ N m b -1 .
Hence, the total number of junction points is N m b -1. For instance, in the case N b = 3, the polygons are all triangles; see Figure 1.

In the sequel, we will denote by P 0 the initial polygon, whose vertices are the fixed points of the maps T i , 0 ≤ i ≤ N b -1, introduced in Notation 8 and Definition 9, on page 5, i.e., P 0 , ⋯, P N b -1 ; see, again, Figure 1, on page 6. ). For any m ∈ N, we introduce the following polygonal sets

P m = P m,k , 0 ≤ k ≤ N m b -1 and Q m = Q m,k , 0 ≤ k ≤ N m b -2 .
M j ,m = j (N b -1) N m b , W j (N b -1) N m b .
We also introduce, for any integer j in 1, ⋯,

(N b -1) N m b -1 :
i. the elementary horizontal lengths:

L m = j (N b -1) N m b ;
ii. the elementary heights:

h j -1, j ,m = W j (N b -1) N m b -W j -1 (N b -1) N m b ;
iii. the geometric angles:

θ j -1, j ,m = (y ′ y), M j -1,m M j ,m ,

which yield the following value of the geometric angle between consecutive edges,

namely, M j -1,m M j ,m , M j ,m M j +1,m , with arctan = tan -1 : θ j -1, j ,m + θ j , j +1,m = arctan L m |h j -1, j ,m | + arctan L m |h j , j +1,m |
.

Proposition 15 (Scaling Properties of the Weierstrass Function, and Consequences [5]

). Since,

for any real number x, W(x) = ∞ n=0 λ n cos 2 π N n b x , one also has W(N b x) = ∞ n=0 λ n cos 2 π N n+1 b x = 1 λ ∞ n=1 λ n cos 2 π N n b x = 1 λ (W(x) -cos (2 π x)) ,
which yields, for any strictly positive integer m and any j in {0, . . . , #V m },

W j (N b -1) N m b = λ W j (N b -1) N m-1 b + cos 2 π j (N b -1) N m b .
By induction, one then obtains that

W j (N b -1) N m b = λ m W j (N b -1) + m-1 k=0 λ k cos 2 π N k b j (N b -1) N m b .

Proposition 16 (Lower Bound and Upper Bound for the Elementary Heights [5]

). For any strictly positive integer m and any j in 0, ⋯, (N b -1) N m b , we have the following estimates:

C i n f L 2-D W m ≤ |W ((j + 1) L m ) -W ( j L m )| ≤ C sup L 2-D W m , m ∈ N, 0 ≤ j ≤ (N b -1) N m b
where the finite and positive constants C i n f and C sup are given by

C i n f = (N b -1) 2-D W min 0≤ j ≤N b -1, W j +1 N b -1 ≠W j N b -1 W j + 1 N b -1 -W j N b -1 and C sup = (N b -1) 2-D W max 0≤ j ≤N b -1 W j + 1 N b -1 -W j N b -1 + 2 π (N b -1) (λ N b -1)
.

One should note, in addition, that these constants C i n f and C sup depend on the initial polygon P 0 .

Theorem 17 (Sharp Local Discrete Reverse Hölder Properties of the Weierstrass Function).

For any natural integer m, let us consider a pair of real numbers

(x, x ′ ) such that x = (N b -1) k + j (N b -1) N m b = ((N b -1) k + j ) L m , x ′ = (N b -1) k + j + (N b -1) N m b = ((N b -1) k + j + ) L m , where 0 ≤ k ≤ N b -1 m -1, and i. if the integer N b is odd, 0 ≤ j < N b -1 2 and 0 < j + ≤ N b -1 2 or N b -1 2 ≤ j < N b -1 and N b -1 2 < j + ≤ N b -1 ; ii. if the integer N b is even, 0 ≤ j < N b 2 and 0 < j + ≤ N b 2 or N b 2 + 1 ≤ j < N b -1 and N b 2 + 1 < j + ≤ N b -1 ⋅
This means that the points (x, W(x)) and x ′ , W(x ′ ) are vertices of the polygon P m,k (see Property 12, on page 5 above), both located on the left-side of the polygon, or both located on the right-side; see Figure 3, on page 9.

Then, one has the following (discrete, local) reverse-Hölder inequality, with sharp Hölder exponent -

ln λ ln N b = 2 -D W , C i n f |x ′ -x| 2-D W ≤ W(x ′ ) -W(x) ⋅ Left -

Sketch of the proof.

We simply use Proposition 16, on page 8.

Cohomology Infinitesimal -Complex Dimensions Series Expansion of the

m m = (ε m ) m = 1 N b -1 1 N m b .
Observe that, clearly, ε m itself -and not just ε In the sequel, it is also useful to keep in mind that the sequence of positive numbers (ε m )

∞ m=0 itself satisfies ε m ∼ 1 N b , as m → ∞ ; i.e., ε m → 1 N b
, as m → ∞. In particular, ε m / → 0, as m → ∞, but, instead, ε m tends to a strictly positive and finite limit.

We also introduce, given any m ∈ N, the m t h intrinsic cohomology infinitesimal, denoted by ε

m > 0, such that ε m = 1 N m b , where ε = 1 N b .
We call ε the intrinsic scale, or intrinsic subdivision scale.

Note that

ε m m = ε m N b -1
and that the m t h intrinsic cohomology infinitesimal ε m is asymptotic equal (when m tends to ∞)

to the m t h cohomology infinitesimal ε m m .
Remark 19. We note that the choice of the m t h intrinsic cohomology infinitesimal, instead of the m t h cohomology infinitesimal, as is done in [START_REF]Weierstrass fractal drums -I -A glimpse of complex dimensions[END_REF] and [START_REF]Weierstrass fractal drums -II -Towards a fractal cohomology[END_REF], will significantly help the computation of the polyhedral prefractal volumes, the polyhedral fractal zeta functions and hence also, of the Complex Dimensions, without any loss of information.

+ i k k, j ,m p, with 0 ≤ k ≤ m, W comp j ε m m = W comp j ε m N b -1 = m k=0 c k, j ,m ε k (2-D W ) ε i k, j ,m p , ( 2
)
where, for 0 denoting the oscillatory period of the Weierstrass Curve, as introduced in [START_REF]Weierstrass fractal drums -I -A glimpse of complex dimensions[END_REF] and where:

≤ k ≤ m, ε k is the k t h intrinsic
i . k, j ,m ∈ Z is arbitrary. i i . c m, j ,m = W comp j N b -1 and, for 0 ≤ k ≤ m -1, c k, j ,m ∈ C is given by c k, j ,m = exp 2 i π N b -1 j ε m-k . ( 3
)
For any m ∈ N, the complex numbers c 0, j ,m+1 , ⋯, c m+1, j ,m+1

and the integers 0, j ,+1 , ⋯, m+1, j ,m+1 respectively satisfy the following recurrence relations:

c m+1, j ,m+1 = W j N b -1 = c m, j ,m (4)
and

∀ k ∈ {1, . . . , m} ∶ c k, j ,m+1 = c k, j ,m , k, j ,m+1 = k, j ,m , (5)
In addition, since relation ( 2) is valid for any m ∈ N ⋆ , we note that the associated Complex Dimensions (i.e., in fact, the Complex Dimensions associated with the Weierstrass function) are

D W -k (2 -D W ) + i k, j ,m p where 0 ≤ j ≤ #V m -1 and 0 ≤ k ≤ m.
This immediately ensures, for the Weierstrass function (i.e., the real part of the Weierstrass complexified function W comp ), that, for any strictly positive integer m and any j in {0, ⋯, #V m },

W j ε m m = 1 2 m k=0 ε k (2-D W ) c k, j ,m ε i k, j ,m p k + c k, j ,m ε -i k, j ,m p k . ( 6
)
Remark 21. The recurrence relations (4) and (5) stated in part i i . of Theorem 20 just above are new and will play a key role in order to establish the recurrence relations satisfied by the local fractal zeta functions; see relation [START_REF] Keller | A simpler proof for the dimension of the graph of the classical Weierstrass function[END_REF] in Theorem 32, on page 20 below. contains the set of vertices V m+k . In particular, the density property 10 (Proposition 10, on page 5), also ensures that the Weierstrass Curve is contained in Γ W m :

D Γ W m , ε m m = M = (x, y) ∈ R 2 , d M , Γ W m ≤ ε
∀ m ≥ k 0 ∶ Γ W ⊂ D Γ W m .

Notation 26 (Lebesgue Measure on R

2 ). In the sequel, we denote by µ L the Lebesgue measure on R 2 .

Definition 27 (The Weierstrass Complex Curve). We place ourselves in the Riemann sphere (or complex projective line) P

Γ W,comp = x, W comp (x) , x ∈ [0, 1] = Γ W + i Γ W,I ,
where Γ W,I is the graph (in R

2 ) associated with the imaginary part W I of the Weierstrass complefixied function W comp and defined, for any real number x, by

W I (x) = ∞ n=0 λ n sin 2 π N n b x . ( 7
V m (ε m m ) = µ L D Γ W m = ε m N m b -1 j =0 N b q=0 α q (N b ) W (N b -1) j + q (N b -1) N m b = ε m N m b -1 j =0 N b q=0 α q (N b ) m k=0 ε k (2-D W ) Re c m,k,(N b -1) j +q ε i k k,(N b -1) j +q,m p , ( 8
)
where

ε m = (N b -1) ε m m = 1 N m b
is the m t h intrinsic cohomology infinitesimal introduced in Definition 18, on page 10, while the real numbers α q (N b ), with 0 ≤ q ≤ N b , are coefficients which depend on N b and where the complex numbers c m,k,(N b -1) j +q ∈ C have been introduced in part i i . of Theorem 20, on page 11, while the integers m,k,(N b -1) j +q ∈ Z are arbitrary (see also Theorem 20, on page 11), with a. When the integer N b is odd:

α 0 (N b ) = α N b -1 (N b ) = -α N b -1 2 (N b ) = -α N b -1 2 +N b -1 (N b ) = N b -2 2 (N b -1) and for 1 ≤ q ≤ N b -2, α N b -1 2 +q (N b ) = -α q (N b ) = -α N m b ,(N b -1) N m b -(N b -1)+q (N b ) = 1 N b -1 , along with α N m b ,(N b -1) N m b -(N b -1) (N b ) = α N m b ,1 (N b ) = N b -2 2 (N b -1) . b. When the integer N b is even: α 0 (N b ) = α N b -1 (N b ) = -α N b 2 (N b ) = -α N b 2 +N b -1 (N b ) = N b -2 2 (N b -1) and for 1 ≤ q ≤ N b -2, α N b 2 +q (N b ) = -α q (N b ) = -α N m b ,(N b -1) N m b -(N b -1)+q (N b ) = - 1 N b -1 , along with α N m b ,(N b -1) N m b -(N b -1) (N b ) = α N m b ,1 (N b ) = N b -2 2 (N b -1)
.

Given the form of the expression in relation (8) just above, it is natural to introduce, for any m ∈ N

⋆ sufficiently large, the associated m t h complex natural polyhedral vol-

ume V m,comp (ε m m ), such that V m,comp (ε m m ) = ε m N m b -1 j =0 N b q=0 α q (N b ) m k=0 c m,k,(N b -1) j +q ε k (2-D W )+i k k,(N b -1) j +q,m p . ( 9
(C) = C ∪ ∞.
Observe that much as for its Euclidean counterpart obtained in our earlier work, [START_REF]Weierstrass fractal drums -I -A glimpse of complex dimensions[END_REF], the m t h polyhedral fractal formula in Theorem 28, on page 16, is expressed at the m t h cohomology infinitesimal ε m m , for all m ∈ N sufficiently large (instead of being expressed at any real number > 0 sufficiently small).

Notation 29 (Minimal and Maximal Values of the Weierstrass Function

W on [0, 1]).
We set

m W = min t ∈ [0,1] W(t) = - 1 1 -λ , M W = max t ∈ [0,1] W(t) = 1 1 -λ .
M (N b -1) j ,m , . . . , M (N b -1) j +N b -1,m
is obtained by substracting from the area of the trapezoid

H (N b -1) j ,m M (N b -1) j ,m M (N b -1) j +N b -1,m H (N b -1) j +N b -1,m
the respective areas of the following N b -1 trapezoids

H (N b -1) j ,m M (N b -1) j ,m M (N b -1) j +1,m H (N b -1) j +1,m ⋮ H (N b -1) j +N b -2,m M (N b -1) j +N b -2,m M (N b -1) j +N b -1,m H (N b -1) j +N b -1,m ; i.e., µ L P m, j +1 = (N b -1) ε m m 2 W (N b -1) j (N b -1) N m b + W (N b -1) j + N b -1 (N b -1) N m b - ε m m 2 N b -2 q=0 W (N b -1) j + q (N b -1) N m b + W (N b -1) j + q + 1 (N b -1) N m b = (N b -1) ε m m 2 W (N b -1) j (N b -1) N m b + W (N b -1) j + N b -1 (N b -1) N m b - ε m m 2 N b -2 q=1 2 W (N b -1) j + q (N b -1) N m b - ε m m 2 W (N b -1) j (N b -1) N m b .
The two-dimensional Lebesgue measure (i.e., area) of the non necessarily convex N bgon Q m, j +1 is obtained in a similar manner.

Lemma 30 (Natural Polyhedral Volume Extension Formula).

Given m ≥ 1 sufficiently large, we set

V m,comp (ε m m ) = ε m V par t i al ,m,comp ε m , ( 10
)
where

V par t i al ,m,comp ε m = N m b -1 j =0 N b q=0 m k=0 α q (N b ) c m,k,(N b -1) j +q ε k (2-D W )+i k k,(N b -1) j +q,m p
and where

ε m = (N b -1) ε m m = 1 N m b is

the m t h intrinsic cohomology infinitesimal introduced in

Definition 18, on page 10, while the real numbers α q (N b ), with 0 ≤ q ≤ N b , have been introduced in Proposition 28, on page 16, and where the complex numbers c m,k,(N b -1) j +q ∈ C have been introduced in part i i . of Theorem 20, on page 11, while the integers m,k,(N b -1) j +q ∈ Z are arbitrary (see Theorem 20,on page 11). We have used the notation V par t i al ,m,comp ε m , since, as mentioned previously, the sum

N m b -1 j =0 N b q=0 m k=0 α q (N b ) c m,k,(N b -1) j +q ε k (2-D W )+i k k,(N b -1) j +q,m p
can be considered as a function of the m t h intrinsic cohomology infinitesimal ε m .

We then introduce, for all sufficiently large m ∈ N ⋆ , Ṽpar t i al ,m,comp as the continuous function defined for all t ∈ [0, ε] by substituting t for ε in the right-hand side of the expression for V par t i al ,m,comp ε m , in relation [START_REF] Hardy | Weierstrass's Non-Differentiable Function[END_REF].

As is explained in [START_REF]Weierstrass fractal drums -I -A glimpse of complex dimensions[END_REF] Then, for all m ∈ N sufficiently large, ζe m admits a (necessarily unique) meromorphic extension to all of C, given, for all s ∈ C, by the following explicit expression:

ζe m (s) = ε m N m b -1 j =0 N b q=0 α q (N b ) m k=0 c k,(N b -1) j +q,m ε s-2+k (2-D W )+i k,(N b -1) j +q,m p s -2 + k (2 -D W ) + i k,(N b -1) j +q,m p , ( 12
)
where ε is the intrinsic scale introduced in Definition 18, on page 10), and where, for 0 ≤ j ≤ N m b -1 and 0 ≤ q ≤ N b , the coefficients α q (N b ) ∈ R have been introduced in Theorem 28, on page 16, while the coefficients c k,(N b -1) j +q,m ∈ C, along with the arbitrary integers k,(N b -1) j +q,m ∈ Z, have been introduced in Theorem 20, on page 11.

More specifically, still for all m ∈ N

⋆ sufficiently large, the function ζe m is well defined and meromorphic in all of C. Furthermore, its (necessarily unique) meromorphic extension (still denoted by ζe m ) is given, for all s ∈ C by the expressions given in the last two equalities of relation [START_REF] Kaplan | The Lyapunov dimension of a nowhere differentiable attracting torus[END_REF] just above.

Moreover, the associated sequence ζe m m ∈ N -initially given (for Re(s) > D W ) by the truncated Mellin transform in relation [START_REF] Herichi | Quantized Number Theory, Fractal Strings and the Riemann Hypothesis: From Spectral Operators to Phase Transitions and Universality[END_REF], on page 20 -satisfies the following recurrence relation, for all values of the positive integer m sufficiently large, and for all s ∈ C,

ζe m+1 (s) = ε ζe m (s) + ε m+1 N m b -1 j =0 N b q=0 α q (N b ) c m+1,(N b -1) j +q,m+1 ε s-2+(m+1) (2-D W )+i m+1,(N b -1) j +q,m+1 p s + 2 (m + 1 -1) -(m + 1) D W + i m+1,(N b -1) j +q,m+1 p +ε m+1 N m+1 b -1 j =N m b N b q=0 α q (N b ) m+1 k=0 c k,(N b -1) j +q,m+1 ε s-2+k (2-D W )+i k,(N b -1) j +q,m+1 p s + 2 (k -1) -k D W + i k,(N b -1) j +q,m+1 p . ( 13
)
This ensures the existence of the limit fractal zeta function ζe W , i.e., the fractal zeta function associated with the Weierstrass Curve Γ W (or, rather, with the Weierstrass IFD), called the global polyhedral effective zeta function and given by ζe

W = lim m→∞ ζe m ,
where the convergence is locally uniform on C, along with the existence of an integer m 0 ∈ N such that, for all m ≥ m 0 , the set of poles of ζe W consists of simple poles and contains the poles of the m t h fractal effective polyhedral zeta function ζe m . More specifically, ζe W is meromorphic in all of C and its meromorphic extension (still denoted ζe

W ) is given, for all s ∈ C, by ζe W (s) = ∞ m=m 0 ε m-m 0 +1 ζe m (s) + ∞ m=m 0 ε m+1 N m b -1 j =0 N b q=0 α q (N b ) c m+1,(N b -1) j +q,m+1 ε s-2+(m+1) (2-D W )+i m+1,(N b -1) j +q,m+1 p s + 2 m -(m + 1) D W + i m+1,(N b -1) j +q,m+1 p + ∞ m=m 0 ε m+1 N m+1 b -1 j =N m b N b q=0 α q (N b ) m+1 k=0 c k,(N b -1) j +q,m+1 ε s-2+k (2-D W )+i k,(N b -1) j +q,m+1 p s + 2 (k -1) -k D W + i k,(N b -1) j +q,m+1 p . ( 14
)
As was mentioned in the introduction, we note that our result is stronger than the one previously obtained in [START_REF]Weierstrass fractal drums -I -A glimpse of complex dimensions[END_REF], where, in particular, the values of the possible Complex Dimensions of the Weierstrass IFD included -2, 0 and 1 -2 k, with k ∈ N arbitrary. As we can see in relation [START_REF] Lapidus | Fractal drum, inverse spectral problems for elliptic operators and a partial resolution of the Weyl-Berry conjecture[END_REF] just above, the poles of the limit effective fractal zeta function ζe W are exactly the same as the Complex Dimensions of the Weierstrass function itself; see Theorem 20, on page 11. Note that, in [START_REF]Weierstrass fractal drums -I -A glimpse of complex dimensions[END_REF], the Complex Dimensions are defined in terms of the volume of the tubular (rather than polyhedral) neighborhood of the Weierstrass IFD.

Sketch of the proof.

i . We first give the explicit expression for the m t h local effective polyhedral zeta function ζe m .

We restrict ourselves to sufficienly large values of m ∈ N, i.e., m ≥ m 0 , for some suitable integer m 0 ∈ N.

We then have that, for Re (s)

> D W , ζe m (s) = ε 0 t s-3 Ṽm,comp (t) d t (for Re (s) > D W ) = ε m N m b -1 j =0 N b q=0 α q (N b ) m k=0 c k,(N b -1) j +q,m ε s-2+k (2-D W )+i k,(N b -1) j +q,m p s + 2 (k -1) -k D W + i k,(N b -1) j +q,m p .
We then have that, for any m ∈ N, We set U = {s ∈ C , 1 < Re(s) < 3}. We note that, for all m ≥ m 0 , we have that

ζe m+1 (s) = ε m+1 N m+1 b -1 j =0 N b q=0 α q (N b ) m+1 k=0 c k,(N b -1) j +q,m+1 ε s-2+k (2-D W )+i k,(N b -1) j +q,m+1 p s + 2 (k -1) -k D W + i k,(N b -1) j +q,m+1 p = ε ζe m (s) + ε m+1 N m b -1 j =0 N b q=0 α q (N b ) c m+1,(N b -1) j +q,m+1 ε s-2+(m+1) (2-D W )+i m+1,(N b -1) j +q,m+1 p s + 2 m -(m + 1) D W + i m+1,(N b -1) j +q,m+1 p +ε m+1 N m+1 b -1 j =N m b N b q=0 α q (N b ) m+1 k=0 c k,(N b -1) j +q,m+1 ε s-2+k (2-D W )+i k,(N b -1) j +q,m+1 p s + 2 (k -1) -k D W + i k
P ζe m 0 ⊂ P ζe m ⊂ U . The series ∞ m=m 0 ε m-m 0 +1 ζe m (s) + ∞ m=m 0 ε m+1 N m b -1 j =0 N b q=0 α q (N b ) c m+1,(N b -1) j +q,m+1 ε s-2+(m+1) (2-D W )+i m+1,(N b -1) j +q,m+1 p s + 2 m -(m + 1) D W + i m+1,(N b -1) j +q,m+1 p + ∞ m=m 0 ε m+1 N m+1 b -1 j =N m b N b q=0 α q (N b ) m+1 k=0 c k,(N b -1) j +q,m+1 ε s-2+k (2-D W )+i k,(N b -1) j +q,m+1 p s + 2 (k -1) -k D W + i k,(N b -1) j +q,m+1 p (16)
is (locally) normally convergent, and, hence, uniformly convergent on U. This ensures the existence of the limit effective global fractal zeta function, i.e., the fractal zeta function associated with the Weierstrass Curve Γ W (or with the Weierstrass IFD), given by ζe

W (s) = lim m→∞ ζe m (s) = ∞ m=m 0 ε m-m 0 +1 ζe m (s) + ∞ m=m 0 ε m+1 N m b -1 j =0 N b q=0 α q (N b ) c m+1,(N b -1) j +q,m+1 ε s-2+(m+1) (2-D W )+i m+1,(N b -1) j +q,m+1 p s + 2 m -(m + 1) D W + i m+1,(N b -1) j +q,m+1 p + ∞ m=m 0 ε m+1 N m+1 b -1 j =N m b N b q=0 α q (N b ) m+1 k=0 c k,(N b -1) j +q,m+1 ε s-2+k (2-D W )+i k,(N b -1) j +q,m+1 p s + 2 (k -1) -k D W + i k,(N b -1) j +q,m+1 p .
Here, and in the remainder of this proof, a (complex-valued) meromorphic function f is viewed as a continuous function with values in P 1 (C), equipped with the chordal metric, and such that, for any pole ω of f , f (ω) takes the value ∞ (for instance, as in [START_REF] Lapidus | Fractal Geometry, Complex Dimensions and Zeta Functions: Geometry and Spectra of Fractal Strings[END_REF], Section 3. 4 and Appendix C ).

More precisely, if P 1 (C) = C ∪ ∞ denotes the Riemann sphere (or complex projective line), we can show that, for the chordal metric, defined, for all (z 1 , z 2 ) ∈ P 1 (C) The sum of this series, i.e., the (uniform) limit fractal zeta function ζe W , is holomorphic on U. We can then deduce that, for all m ≥ m 0 , the zeta function ζe m is meromorphic on C.

2 by ∥z 1 , z 2 ∥ = |z 1 -z 2 | 1 + |z 2 1 | 1 + |z 2 2 | , if z 1 ≠ ∞ and z 2 ≠ ∞ and ∥z 1 , ∞∥ = 1 1 + |z 2 1 | , if z 1 ≠ ∞ ,
Moreover, the counterpart of the results obtained in [START_REF]Weierstrass fractal drums -I -A glimpse of complex dimensions[END_REF] for the sequence of tube zeta functions associated with the Weierstrass IFD, which admit a meromorphic continuation to all of C, obviously hold for the sequence of polyhedral tube zeta functions: hence, ζe m is meromorphic on C, with only simple poles, as specified in Theorem 33, on page 24.

In the next new result, which is now an immediate corollary of Theorem 32, on page 20, we define implicitly the Complex Dimensions of the Weierstrass Curve Γ W (or, rather, of the

D W -k (2 -D W ) + i p , ( 17
)
where the integers k ∈ N and ∈ Z are arbitrary (compare the corresponding result in Theorem 20, on page 11).

We note that they only partly coincide with the possible poles of the tube fractal zeta function ζe,tube For the exceptional cases, we refer to [START_REF]Weierstrass fractal drums -I -A glimpse of complex dimensions[END_REF] for a closely related discussion.

Remark 34. Note that, in Theorem 33 just above, the fact that -in light of relation (32) in Theorem 12, on page 20, and of the expression of the coefficients α q and c k, j m given, respectively, in relation [START_REF]Polyhedral neighborhoods vs tubular neighborhoods: New insights for fractal zeta functions and complex dimensions[END_REF] in Theorem 28, on page 16 and in part i i . of Theorem 20, on page 11 -the possible Complex Dimensions are all actual Complex Dimensions of the Weierstrass IFD follows from the fact that the coefficients of the fractal power series for ζe W in relation ( 14) are nonzero, which implies that the residues of ζe W at each possible pole is nonzero.

In the theory of Complex Dimensions (see, e.g., [START_REF] Lapidus | Fractal Geometry, Complex Dimensions and Zeta Functions: Geometry and Spectra of Fractal Strings[END_REF], [START_REF]Fractal Zeta Functions and Fractal Drums: Higher-Dimensional Theory of Complex Dimensions[END_REF], [START_REF]An overview of complex fractal dimensions: From fractal strings to fractal drums, and back[END_REF]), a geometric object is said to be fractal if it admits at least one nonreal Complex Dimension (defined as a pole of the associated geometric or fractal zeta function.) [In the present context, the fractal zeta function is the global polyhedral effective zeta function in the case of the

= D W -k (2 -D W ) + i p ,
where k ∈ N, 0 ≤ k ≤ m and ∈ Z is arbitrary.)

Revisiting the Computation of the Minkowski Dimension

Along the lines of our polyhedral fractal zeta functions, we hereafter revisit the computation of the box-counting (or box) dimension D W of the Weierstrass Curve Γ W . Contrary to (classical) covers of the fractal under study by means of balls (in dimension 2, disks or squares; see Definition 36, on page 25 below, we simply use our polyhedral neighborhood (see Definition 23, on page 12), which is undoubtedly the most natural cover of Γ W .

We therefore announce a fully rigorous proof of the fact that the Minkowski dimension -which, as is well known, coincides with the box (or box-counting) dimension D W of the Weierstrass Curve Γ W -is given by the expected formula, originally conjectured by Benoît Mandelbrot in [START_REF] Mandelbrot | The Fractal Geometry of Nature[END_REF] (in Chapter XI, on top of page 390). Note that B. Mandelbrot only mentioned he was talking about the fractal dimension and did not specify wether it was the Hausdorff or the Minkowski dimension of the curve. In our present context (see Theorem 37, on page 26 below), we have that D W = 2 - Note that earlier proofs of this fact were either not fully rigorous or complete. In addition, they rely on unexplicit estimates, whereas we work with explicit estimates that had not been obtained before, in relation with he local Hölder and reverse Hölder continuity of the Weierstrass function W.

Definition 36 (Box-Counting Dimension).

As can be found, for instance, in [START_REF] Falconer | The Geometry of Fractal Sets[END_REF], we recall that (when it exists) the box-counting dimension (or box dimension, in short), of Γ W , is given by

D W = -lim δ→0 + ln N δ (Γ W ) ln δ , (⋄)
where N δ (Γ W ) stands for any of the following quantities: i. the smallest number of sets (here, subsets of R

2 ) of diameter at most δ that cover Γ W on [0, 1[ ; ii. the smallest number of closed balls (disks, here) of radius δ that cover Γ W on [0, 1[ ; iii. the smallest number of cubes (squares, here) of side δ that cover Γ W on [0, 1[; iv. the number of δ-mesh cubes (squares, here) that intersect Γ W on [0, 1[; v. the largest number of disjoint balls (disks, here) of radius δ with centers in Γ W on [0, 1[. Furthermore, for the Weierstrass Curve Γ W -as, more generally, for any bounded subset of Euclidean space -the box-counting dimension coincides with the Minkowski dimension. Note that we thus dispose, with the sequence of polyhedral domains D Γ W m m≥m 0 , of a nonusual (but admissible, in the sense of Definition 36, on page 25 above) sequence of covers.

However, D Γ W m m≥m 0 is the most natural and optimal sequence of covers of the Weier- strass Curve Γ W , insofar that for any m ≥ m 0 , each domain D Γ W m contains the Weierstrass Curve Γ W , while, at the same time, the sequence D Γ W m m≥m 0 converges to Γ W . This, in particular, means that when m → ∞, the two-dimensional Lebesgue measure (i.e., area) of each polygon belonging to D Γ W m tends to 0. , where the diameter of a polygon is to be understood in the sense of the largest distance between any pair of vertices of the considered polygon.

In our context, the cohomology infinitesimal ε

Notation 1 (

 1 ∞}, we write, e.g., [a, b[ = [a, b) and similarly for other intervals. Weierstrass Parameters). In the sequel, λ and N b are two real numbers such that 0

Proposition 10 (

 10 Density of the Set V dense in the Weierstrass Curve Γ W . Definition 11 (Adjacent Vertices, Edge Relation). For any m ∈ N, the prefractal graph Γ W m is equipped with an edge relation ∼ m , as follows: two vertices X and Y of Γ W m (i.e. two points belonging to V m ) are said to be adjacent (i.e., neighboring or junction points) if and only if the line segment [X , Y ] is an edge of Γ W m ; we then write X ∼ m Y . Note that this edge relation depends on m, which means that points adjacent in V m might not remain adjacent in V m+1 .

  denotes the number of elements in the finite set V m . i i i . The prefractal graph Γ W m has exactly (N b -1) N m b edges. iv. The consecutive vertices of the prefractal graph Γ W m are the vertices of N m b simple nonregular polygons P m,k with N b sides (or N b -gons). For any strictly positive integer m, the junction point between two consecutive polygons is the point

1 yFigure 1 . 1 , 1 , P 1 , 2 , in the case when λ = 1 2 and N b = 3 . (See also Figure 2 , on page 7

 1111122327 Figure 1. The initial polygon P 0 , and the polygons P 1,0 , P 1,1 , P 1,2 , in the case when λ = 1 2 and N b = 3. (See also Figure 2, on page 7.)

1 yFigure 2 . 1 , Γ W 2 , Γ W 3 , Γ W 4 , Γ W 5 , in the case when λ = 1 2 and N b = 3 .

 121234523 Figure 2. The prefractal graphs Γ W 0 , Γ W 1 , Γ W 2 , Γ W 3 , Γ W 4 , Γ W 5 , in the case when λ = 1 2 and N b = 3. For example, Γ W 1 is on the right side of the top row, while Γ W 4 is on the left side of the bottom row.

side vertices Right -side vertices Bottom vertex 1 x -1 1 yFigure 3 .

 13 Figure 3. The left-side and right-side vertices.

  m m -depends on m. In addition, since N b > 1, ε m m satisfies the following asymptotic behavior, ε m m → 0 , as m → ∞, which, naturally, results in the fact that the larger m, the smaller ε m m . It is for this reason that we call ε m m -or rather, the infinitesimal sequence ε m m ∞ m=0 of positive numbers tending to zero as m → ∞, with ε m m = (ε m ) m , for each m ∈ N -an infinitesimal. Note that this m t h cohomology infinitesimal is the one naturally associated to the scaling relation of Proposition 15, on page 8.

Theorem 20 (

 20 Complex Dimensions Series Expansion of the Weierstrass Complexified Function W comp [6], [7], [8] and of the Weierstrass function W). For any sufficiently large positive integer m and any j in {0, ⋯, #V m -1}, we have the following exact expansion, indexed by the Complex Codimensions k (D W -2)

  cohomology infinitesimal, introduced in Definition 18, on page 10, with p = 2 π ln N b

  Neighborhood [5]). Given m ∈ N sufficiently large (so that ε m m be a sufficiently small positive number), we define the (m, ε m m )-neighborhood of the m t h prefractal approximation Γ W m as follows:

Definition 23 (

 23 Sequence of Domains Delimited by the Weierstrass IFD -Polyhedral Neighborhood of the Weierstrass Curve [4]). We introduce the sequence of domains delimited by the Weierstrass IFD, or polygonal neighborhood of the Weierstrass Curve, as the sequence D Γ W m m ∈ N of open, connected polygonal sets (P m ∪ Q m ) m ∈ N , where, for each m ∈ N, P m and Q m respectively denote the polygonal sets introduced in Definition 13, on page 6. Given ∈ N, we call D Γ W m the m t h polyhedral neighborhood (of the Weierstrass Curve Γ W ); see Figure 4, on page 12. Theorem 24 (The Nested Neighborhoods). i . Given m ∈ N sufficiently large, there exists k 1 ∈ N such that, for all k ≥ k 1 , the polygonal neighborhood D Γ W m introduced in Definition 23, on page 12 contains, but for a finite number of wedges, the (m + k, ε m+k m+k )-neighborhood; see Figures 4-7, on pages 12-14. i i . Given m ∈ N sufficiently large, there exists k 2 ∈ N such that, for all k ≥ k 2 , the (m, ε m m )neighborhood contains the polygonal neighborhood D Γ W m+k introduced in Definition 23, on page 12; see Figure 8, on page 14 and Figure 9, on page 15. i i i . Given m ∈ N sufficiently large, there exists k 3 ∈ N such that, for all k ≥ k 3 , the polygonal neighborhood D Γ W m introduced in Definition 23, on page 12 contains the polygonal neighborhood D Γ W m+k .

Figure 4 .

 4 Figure 4. The polygonal neigborhood D Γ W 2 , in the case when λ = 1 2 and N b = 3.

Figure 5 .

 5 Figure 5. The exterior boundary of the polygonal neigborhood D Γ W 2 (in red), and the tubular neighborhood D Γ W 7 , ε 7 7 , in the case when λ = 1 2 and N b = 3.

Figure 6 . 2 and

 62 Figure 6. The polygonal neigborhood D Γ W 3 , in the case when λ = 1 2 and N b = 3.

Figure 7 .

 7 Figure 7. The exterior boundary of the polygonal neigborhood D Γ W 3 (in red), and the tubular neighborhood D Γ W 7 , ε 7 7 , in the case when λ = 1 2 and N b = 3.

Figure 8 . 2 and N b = 3 .Figure 9 . 3 and 2 and N b = 3 .

 8239323 Figure 8. The polygonal neigborhood D Γ W 5 , in the case when λ = 1 2 and N b = 3.

1

  (C) = C ∪ ∞. We define the Weierstrass Complex Curve as the graph in P 1 (C) = C ∪ ∞ associated with the Weierstrass complefixied function W comp ; i.e., the set denoted by Γ W,comp such that

Figure 10 .

 10 Figure 10. The orthogonal projections H j ,m of the vertices M j ,m , for 0 ≤ j ≤ #V m -1, onto the horizontal line y = m W , where m W is introduced in Notation 29, on page 17.

  ,(N b -1) j +q,m+1 p . (15) i i . For any integer m ≥ m 0 , we denote by P ζe m ⊂ C the set of poles of the zeta function ζe m .

  we have, thanks to the local uniform convergence of the series on C, lim m→∞ ∥ ζe m , ζe W ∥ = 0 . Indeed, for any η > 0 and any compact set K ⊂ C, we can choose m 0 ∈ N ⋆ such that, for all m ≥ m 0 and s ∈ K, we have that | ζe m (s) -ζe W (s)| ≤ η , and hence, still for all m ≥ m 0 and s ∈ K, ∥ ζe m (s), ζe W (s)∥ ≤ | ζe m (s) -ζe W (s)| ≤ η .

m.

  In particular, we then deduce that the following possible Complex Dimensions previously obtained in[START_REF]Weierstrass fractal drums -I -A glimpse of complex dimensions[END_REF], i.e., 1 -2 k + i p, with k ∈ N and ∈ Z, along with -2 and 0, are not actual Complex Dimensions, in the sense introduced in the present paper.Moreover, the Complex Dimensions associated with DW , i.e., D W -k (2 -D W ) + i p ,where k ∈ N and ∈ Z are arbitrary, are actual Complex Dimensions of Γ W . Hence, all of the possible Complex Dimensions of Γ W , as given by relation[START_REF]An overview of complex fractal dimensions: From fractal strings to fractal drums, and back[END_REF] above, are exact and they are all simple.

Theorem 37 .

 37 The box-counting dimension -or, equivalently, the Minkowski dimension -of the Weierstrass Curve Γ W exists and is equal to DW = 2 + ln λ ln N b .Sketch of the proof.We simply apply the result given in Corollary 25, on page 14. We then have the existence of aninteger m 0 ∈ N such that ∀ m ≥ m 0 ∶ Γ W ⊂ D Γ W m .In other words, for all integers m ≥ m 0 , the polyhedral neighborhood D Γ W m covers the Weierstrass Curve Γ W .

mm

  plays the role of the elementary diameter δ in Definition 36, on page 25. Moreover, the area of each polygon belonging to D Γ W m also corresponds to the area of the reunion of (N b -1) h m ε m m elementary and smaller polygons, each of diameter ε m m

Weier- strass Function Definition 18 (m t h Cohomology Infinitesimal [5], [6] and m t h Intrinsic Cohomology Infinites- imal

	number ε	t h cohomology infinitesimal (of Γ W ) the m > 0 which also corresponds to the elementary horizontal length introduced in part i . m
	of Definition 14, on page 7; i.e., ε

). From now on, given any m ∈ N, we will call m

  )

	Theorem 28 (Exact Expression for the Volume of the m Natural Polyhedral Volume)). Given m ∈ N ⋆ sufficiently large, the volume (or two-dimensional t h t h Polyhedral Neighborhood (or m
	Lebesgue measure) V m (ε	m m ) of the m	t h polygonal neighborhood D Γ W m , or m	t h natural polyhe-
	dral volume, is given by

  )

	By considering, in the same manner, given m ∈ N ⋆ sufficiently large, the volume (or two-
	dimensional Lebesgue measure) V I,m (ε ciated with the m t h prefractal approximation to the graph Γ W,I of the imaginary part W I of the m m ) of the m t h polygonal neighborhood D Γ W,I m , asso-
	Weierstrass complefixied function W comp (see Definition 27, on page 15 above), we note that the complex volume V m,comp (ε m m ) can be envisioned as the volume associated with the Weierstrass
	complex curve Γ W,comp (see again Definition 27, on page 15) in the Riemann sphere (or complex projective line) P 1

Local and Global Polyhedral Effective Zeta Functions: Fractal Power Expansions and

  Ṽpar t i al ,m,comp associated with the m t h natural polyhedral complex volume V m,comp introduced in Theorem 28, on page 16. Alternatively, Ṽm,comp will be called the m

				(in the case of the ordinary Euclidean volume), one can think of
	ε	m Ṽpar t i al ,m,comp t	m as being the effective polyhedral volume of the m	t h prefractal approxi-
	mation to the Weierstrass Curve.
	Notation 31 (Natural Polyhedral Complex Volume Extension).
		For the sake of simplicity, given m ∈ N sufficiently large, we will from now on call
	the m	t h natural polyhedral complex volume extension, the volume extension func-
	tion Ṽm,comp = ε
					t h
	effective polyhedral complex volume.
					t h natural volume extension,
	the volume extension function tively, Ṽ t ube m will be called the m Ṽ t ube m t h effective tubular volume. associated with the the m	t h tubular volume V m t ube	. Alterna-
	Theorem 32 (ε	t	s-3 ε
				0

m

In the same way, as is done in

[START_REF]Weierstrass fractal drums -I -A glimpse of complex dimensions[END_REF]

, and given m ∈ N, we call the m Recurrence Relations). Given m ∈ N sufficiently large, we introduce the m t h local polyhedral effective zeta function ζe m , such that, for all s ∈ C with Re (s) > D W , ζe m (s) = m Ṽm,comp (t) d t , (11) where Ṽm = ε m Ṽpar t i al ,m,comp is the m t h effective polyhedral complex volume, introduced in Notation 31, on page 20 above, with ε denoting the m t h intrinsic cohomology infinitesimal introduced in Definition 18, on page 10.

  associated Weierstrass IFD) as the poles of the polyhedral global effective zeta function ζe W . The phrase "possible Complex Dimension" then refers to ω ∈ C which is either an actual ComplexDimension or appears as an ω-exponent in the fractal power series expansion[START_REF] Keller | A simpler proof for the dimension of the graph of the classical Weierstrass function[END_REF], on page 21, but is such that the corresponding residue vanishes: res ζe

	Theorem 33 (Complex Dimensions of the Weierstrass Curve). The possible Complex Dimen-
	sions of the Weierstrass Curve (or of the Weierstrass IFD) are given by Theorem 32, on page 20. Fur-
	thermore, if they are actual poles, they are simple poles. They are given as follows:

W , ω = 0.

Fractality of the The Weierstrass Curve and of the Prefractal Approximations).

  Weierstrass Curve Γ W . (resp., m t h local polyhedral effective zeta function of its m t h prefractal approximation, with m ≥ m 0 ).] In addition, given d ∈ R, it is fractal in dimension d if it has at least one nonreal Complex Dimension with real part d . (Note that nonreal Complex Dimensions come in complex conjugate pairs.) In particular, it is principally fractal if it is fractal in dimension d (here, D W ), the abscissa of convergence of the associated fractal zeta function, which is the largest possible value of d . We can now state the following corollary of Theorems 32, on page 20 and 33, on page 24. The same is true for the m t h prefractal approximation Γ W m to the Weierstrass Curve Γ W , with m ∈ N ⋆ sufficiently large, except for the fact that Γ W m is only fractal in finitely many dimensions -namely, d = D W -k (2 -D W ), with k ∈ N such that 0 ≤ k ≤ m. (Indeed, the Complex Dimensions of Γ W m are given by d

	Corollary 35 (

The Weierstrass Curve Γ W (or the Weierstrass IFD) is fractal, and even principally fractal, as well as fractal in infinitely many Complex Dimensions d -namely, d = D W -k (2 -D W ), with k ∈ N arbitrary).

Funding. The research of M. L. L. was supported by the Burton Jones Endowed Chair in Pure Mathematics, as well as by grants from the U. S. National Science Foundation.


Pascal Giorgi

Thomas Izard

Arnaud Tisserand

Comparison of Modular Arithmetic Algorithms on GPUs

Keywords: Multiple precision arithmetic, mathematical library, CUDA

We present below our first implementation results on a modular arithmetic library for cryptography on GPUs. Our library, in C++ for CUDA, provides modular arithmetic, finite field arithmetic and some ECC support. Several algorithms and memory coding styles have been compared: local, shared and register. For moderate sizes, we report up to 2.6 speedup compared to state-of-the-art library.

Introduction

Modular operations on large integers are used in many applications such as cryptography, coding and computer algebra. Efficient algorithms and implementations are required for a ± b mod p, a × b mod p where a, b and p are multiple precision integers and p is prime. Those operations are required in finite field arithmetic over F p and in elliptic curve cryptography (ECC) where sizes are about 200-600 bits. Graphic processor units (GPUs) are used in high-performance computing systems thanks to their massively multithreaded architectures. But due to their specific architecture and programming style, porting libraries to GPUs is not simple even using high-level tools such as CUDA [8].

This paper presents our first implementation results on modular arithmetic for large integers, arithmetic over F p and ECC scalar multiplication application on GPUs. This work is a part of a software library called PACE [START_REF] Imbert | A library for prototyping the computer arithmetic level in elliptic curve cryptography[END_REF]. This library is aimed at providing a very large set of mathematical objects, functions and algorithms to facilitate the writing of arithmetic applications. Below, we only deal with the use of GPUs as accelerators for parallel computations with different sets of data. In asymmetric cryptography applications, this kind of parallelism level is required for servers on which cyphering and/or digital signatures are computed for parallel sessions. This paper presents sequential arithmetic operations (one thread per operation) but for massively parallel computations on independent data sets. We compare our code with mpF q library [START_REF] Gaudry | The mpFq library and implementing curve-based key exchanges[END_REF] and the work from [START_REF] Szerwinski | Exploiting the power of gpus for asymmetric cryptography[END_REF].

Large Integers on GPU

We consider modular arithmetic for arbitrary values of the modulo (i.e. p) and size N of numbers in the range 160-384 bits. The operands a and b are integers in the range [0, p -1] (reduced values). Standard modular arithmetic algorithms and some implementation guidelines on standard processors may be found in good books such as [START_REF] Hankerson | Guide to Elliptic Curve Cryptography[END_REF]. But this is not the case for GPU implementation. The design and optimization of arithmetic operations is done using a complex trade-off between: number representation (width, radix), the algorithm(s) and some architecture constraints (type and number of functional units and memory latencies). Our goal is to define an efficient layer which can be used on a single thread, since our parallelism is only on the data. We will see that the memory mapping of the integers into various memories is a key element for this kind of GPU implementation.

Integers Representation and Memory

Large integers are usually stored into an array of words such as depicted below.

a = a n-1 a n-2 . . . a 0 = n=-1 i=0 a i β i
The word size depends on the functional units characteristics. This size leads to various values for the radix β. The choice of β and the word's datatype is clearly related to GPUs capabilities. According to NVIDIA CUDA programming guide [8], the possible native datatype for a word can be either 32-bit integer, 32-bit or 64-bit floating-point (FP) numbers. Considering current lower1 throughputs of the 64-bit FP numbers compare to 32-bit datatype on newest GPU (GT200 core), one can avoid for the time being the use of 64-bit FP number as word.

Considering 32-bit arithmetic units on GPU, one can issue one of the following operations in 4 clock cycles on a GPU multiprocessor: one exact multiply-and-add (MAD) operation on 10-bit operands with float, one exact multiplication on 16-bit operands with int and one exact addition with carry on 16-bit operands with int.

Looking at the school book multiplication It seems clear from this comparison that 32-bit integers with β = 2 16 is a better choice than 32-bit FP numbers. Moreover, a smaller word number saves memory. Furthermore, bit manipulations remain easier with integer representation.

The main difficulty with GPU is to fully benefit from high memory bandwidth within the global memory (the RAM of the GPU). In particular, one needs to design code such that many concurrent memory accesses from many concurrent threads of a multiprocessor can be coalesced into a single transaction. The basic idea is that the i-th thread of a multiprocessor needs to read/write the i-th data in a particular segment of memory. See [8,Chapter 5] for further details.

In our computational model, integers are not gathered between threads and using linear array to store integer's words leads to non-coalesced memory patterns. To fulfill this requirement, we provide functionality to load (resp. store) integer from (resp. to) global memory. These functions change the words order to ensure coalesced access. We define the following rule to reorder words of a same variable across all threads: i-th thread a = a 0 . . . a n-2 a n-1 (i + 1)-th thread a = ā0 . . . ān-2 ān-1

Global memory reordering a 0 ā0 . . . a n-2 ān-2 a n-1 ān-1 To ensure memory alignment, two consecutive words of one integer of one thread are distant of i word addresses, where i is the minimal power of 16 greater or equal to the total number of threads launched on the GPU.

Global memory accesses are quite costly and must be avoided. Therefore, one would prefer to store integers in a non-coalesced form in another memory region (local memory or shared memory). Fetching all data in local or shared memory at the beginning of a GPU kernel function will ensure that all reads from global memory are coalesced.

Large Integer in GPUs Registers

For efficiency reason it would be interesting to map large integers directly into GPU registers. This should save at least the latency of memory load/store instructions. However, as all hardware registers, GPU's registers are not indexable. This means that no arrays can be mapped into registers. The only way to map large integers into registers is to design specific structures with one variable per integer word. This approach is feasible but not used for generic and portable code.

Fortunately, C++ and template metaprogramming [START_REF] Alexandrescu | Modern C++ design: generic programming and design patterns applied[END_REF] can help. In particular, one can design a generic structure which defines fixed number of variables with compile time indexable access. t e m p l a t e <u i n t N> s t r u c t I n t e g e r R e g : I n t e g e r R e g <N-1> { u i n t word ; I n t e g e r R e g ( ) : I n t e g e r R e g <N-1>() { word = 0 ; } } ; t e m p l a t e <> s t r u c t I n t e g e r R e g <1> { u i n t word ; I n t e g e r R e g ( ) { word = 0 ; } } ; # d e f i n e AT( x , i ) x . I n t e g e r R e g <i +1 >:: word IntegerReg structure defined in Figure 2 uses C++ recursive inheritance to define a proper number of variables to store an N -word integer. One can access a particular variable from this structure by using the index given within the template parameter. The macro AT defines the proper accessor to the (i+1)-th word of the integer x. The only drawback of this structure is that all indices need to be known at compile time, meaning no runtime loop can be executed. Here again C++ templates come to rescue to define compile time loop as illustrated in Figure 3.

Compile time loop consists of completely unrolling loops by duplicating code with proper register values. This technique can be very efficient if the compiler is able to discover registers reusability. In case of GPU, nvcc compiler (the CUDA compiler)

Modular Arithmetic on Large Integers

We investigate GPU implementation of the 3 basic modular arithmetic operations on large integers: addition, subtraction and multiplication. We do not investigate division since it can be avoided in most of our target applications. Memory accesses within GPU can be very costly depending where data are localized. For instance one access to global memory costs around 400-600 cycles while one register or one shared memory access cost around 10 cycles. Minimizing data access and intermediate variables is thus critical.

For the modular addition/subtraction the basic algorithm can suffice since it requires only an integer addition (potentially with a carry) and a conditional subtraction. Moreover, this operation can be done inplace (a += b mod p) with only one extra register for the carry propagation. We use classical integer algorithm with carry propagation along limbs as described at [4, pp.30].

In the case of register version of integers, we have β = 2 32 and no instruction is available to get the carry of addition or subtraction of words. Thus, we need to calculate explicitly the carry using the following trick: the carry of a + b is equal to the result of the test a + b < a. A similar trick is available for subtraction.

Separated operation and reduction become too costly for multiplication. The product a × b is 2n words long. Then the reduction to a single n words number is close to a division by p. In that case, the operation and the reduction are interleaved using a bit-serial scan of one operand (usually, the multiplier a). At each step, the partial product a i,j × b (where a i,j is the jth bit of the word a i ) a is accumulated and reduced modulo p. We use the well known Montgomery's algorithm to avoid the use of multi-precision division [START_REF] Montgomery | Modular multiplication without trial division[END_REF]. In order to reduce memory usage we use an interleaved Montgomery's method which consists of doing multiplication and reduction at the same time. Our implementation is based on the Finely Integrated Operand Scanning (FIOS) method described in [START_REF] Kaya Koc | Analyzing and comparing montgomery multiplication algorithms[END_REF] which required only 3 extra words. Using β = 2 16 with 32-bit integers, we do not suffer from carry propagation along word additions.

We compare our arithmetic layer performance when data are localized in the different memory regions: local memory, shared memory and registers. In Table 1, we report computation times of modular operations for N in 160-384 bits. Operands are chosen randomly with full limb occupancy and the GPU kernel function consists in 1024 threads within 64 blocks. Each thread loads its two operands from the global memory to the chosen memory space (local, shared, register) and then compute a loop of 10 000 operations (using result in the loop as a new operand for the next operation). We perform our GPU computation on a Geforce 9800GX2 card. One can see from this table that GPU data location in memory is important, and register usage can improve speed. One can also remark that multiplication for N > 256 registers version becomes less efficient than shared version. This comes from the CUDA compilation chain which is not able to handle code with large register usage. Note that each thread can use at most 128 registers of 32-bit integers.

In order to demonstrate benefits of GPU for multiprecision modular arithmetic, we compare our layer with the mpF q library [START_REF] Gaudry | The mpFq library and implementing curve-based key exchanges[END_REF] which is currently the best library to perform modular arithmetic on CPU with modulo of moderate size (i.e. less than 600 bits). In Table 2 we report the comparison of our best GPU modular arithmetic implementation (register version) with mpF q library. Benchmarks are identical as in previous section, except that mpF q handles sequentially the 1024 threads on a Core(TM)2 Duo E8400-3GHz processor. For N less than 384 bits GPU beats CPU calculation. However, speedup factors are moderate since CPU processor has higher frequency than GPU multiprocessor. Moreover, mpF q uses well tuned code with SIMD SSE-2 instructions, which allows to be more efficient implementation of large integers than our GPU implementation.

ECC Application

In order to evaluate our library in a realistic application, we use it for scalar multiplications. This is the main operation required in ECC. For instance, servers performing many digital signatures have to compute parallel and independent scalar multiplications. Below, we consider an elliptic curve (EC) E(F p ) defined over the prime finite field F p (a is the parameter of the curve). See [START_REF] Hankerson | Guide to Elliptic Curve Cryptography[END_REF] for background, details and notations. In order to avoid modular inversions, we use the Jacobian coordinates system (P = (X, Y, Z) where X, Y and Z are in F p ). Montgomery representation is used for fast multiplications. All threads share p and a parameters but points P and Q are different. All operations presented in this section are used in the same way. Data are sent to the GPU global memory in coalesced form. The GPU kernel loads the selected memory (local, shared or registers) from the global memory, converts necessary values to the Montgomery form, performs the curve operation, and converts back to standard form. Then results are sent from the GPU to the CPU.

Point Addition and Doubling

Two basic operations are used for points P and Q on E(F p ): point addition P + Q (for P = Q) and point doubling 2P (specific addition for P = Q). Those point operations are defined using several additions, multiplications over F p . Point addition and doubling require respectively 2 and 3 extra integers (and 6 for P and Q). These intermediate values are stored in same memory than the other data.

Timings are reported in Table 3 for the three considered memory schemes: local, shared and registers. The computation was launched on 1024 threads (64 blocks of 16 threads) and 100 operations for each thread. Due to the limited number of registers for one thread, the CUDA compiler was unable to build the code for N greater than 192 bits (ror is reported for compilation aborted due to "run out of register" compiler error). For both operation, local memory is the worst due to its high latency. Except for point addition with N = 160, the shared memory implementation is the best. One drawback of shared memory is the fact it limits the number of threads per block to 16. In the shared memory version, the huge timing gap between N = 256 and N = 384 is due to the card occupancy which is twice for N = 256 (or less) than for N = 384. For point doubling, occupancy is the same for all N and the factor 2 is due to the computations. We notice that register version is slower than shared version due to the fact that the compiler was unable to use correctly all the registers it has and put data in local memory.

In Table 4, we compare our best implementation with the PACE library coupled to mpF q running on the CPU for the same data. The shared memory version is at least twice faster than the CPU version for point doubling, but is no more than twice for addition. The difference comes from the number of Montgomery multiplications required in each operation: 9 for point doubling and 16 for point addition and the fact that the speedup factor for one Montgomery multiplication is small.

Scalar Multiplication

Scalar multiplication Q = [k]P = P + • • • + P (with k additions of point P )
where k is an N -bit integer and P a point of E(F p ), is the main operation in ECC protocols (see [START_REF] Hankerson | Guide to Elliptic Curve Cryptography[END_REF] for details and algorithms). We use a right-to-left double-and-add algorithm since it does not require pre-computations and additional storage. It uses point addition and point doubling for each step of a loop over the bits of k. Furthermore, it allows some threads to follow the same execution pattern as the base point is doubled at each step of the loop. Table 5 reports times and operation throughput (in [k]P/s) for our best implementation (i.e., the shared version) and mpF q +PACE implementations. Our tests use the same N for k and F p and the Hamming weight of k equal to N/2. For N = 384, our scalar multiplication on GPU is slower than the CPU version (computed with mpF q +PACE). We noticed that the speedup factor decreases for larger numbers. This is mainly due to the small factor between the mpF q +PACE version and our shared memory version of the Montgomery multiplication (1.6 for N = 160 and 1.3 for N = 256). In [START_REF] Szerwinski | Exploiting the power of gpus for asymmetric cryptography[END_REF], a seminal GPU implementations of [k]P is provided with N = 224. Their throughput is about 1412.6 [k]P/s using mixed affine-Jacobian coordinates and left-toright double-and-add algorithm. Our implementation is a little better than this result.

w-NAF Implementation of Scalar Multiplication

Some signed-digit representations are used for recoding k in the [k]P . This is motivated by the fact that point subtraction on an EC is just as efficient as addition. Among these representations, w-digit windows non-adjacent forms (w-NAF) are frequently used in ECC [START_REF] Hankerson | Guide to Elliptic Curve Cryptography[END_REF]. Using w-NAF, k = l-1 i=0 k i 2 i where non-zero digits k i are odd, |k i | < 2 w-1 and at most one digit of any w-digit window is non-zero. Thus w-NAF recoding decreases the number of point additions/subtractions (≈ 1/(w + 1) w.r.t. 0.5 for binary). It is also used as a simple countermeasure against side channel attacks [START_REF] Byrne | Comparison of simple power analysis attack resistant algorithms for an ECC[END_REF].

Several versions have been implemented: w ∈ {2, 3} as well as basic and optimized storage. The basic storage uses a complete word for each w-NAF digit. The optimized storage uses the minimal number of bits for a w-NAF digit (i.e., 2 bits for w = 2 and 3 bits for w = 3) through a dedicated multiple precision storage. For w = 3, 3P is precomputed and stored. Addition of -P and -3P are obtained by subtracting respectively P and 3P . Recoding from integer k into w-NAF is performed on the GPU.

Due to nvcc compiler limitations, it was not possible to compile more than N = 256 bits with w = 2 and N = 160 for w = 3. Using a w = 3 leads to 10% speed improvement compared to w = 2. The number of operations is reduced, but the additional internal value 3P puts a too high pressure on scheduling. Surprisingly, basic and optimized storage of the w-NAF recoding gives very close results.

Conclusion and Future Prospects

In this work we report our first implementation results on modular arithmetic for large integers on GPUs which achieves a speedup of 2.6 compared to state-of-the-art library. We show that porting modular arithmetic algorithms on GPUs is not direct. Our long term goal is to design a high-performance arithmetic library for cryptography.

Figure 1 .

 1 Figure 1. Global memory organization for coalesced integer access

Figure 2 .

 2 Figure 2. A register-compliant structure for large integer with indexable access.

Table 1 .

 1 Computation times in ns for modular addition and multiplication.

			addition			multiplication
	N	local	shared register	local shared register
	160	9	2.3	0.7	88	40	22
	192	11	2.3	0.7	125	51	33
	224	23	5.0	1.1	172	107	55
	256	26	3.1	1.5	214	80	81
	384	38	7.4	3.9	673	221	261

Table 2 .

 2 Time comparison in ns of software library mpFq and our best GPU implementation.

			addition mod p		multiplication mod p
	N	our impl. mpfq speedup	our impl. mpfq speedup
	160	0.7	15	×21	22	64	×2.9
	192	0.7	16	×22	33	70	×2.1
	224	1.1	21	×19	55	105	×1.9
	256	1.5	21	×14	81	109	×1.3
	384	3.9	30	×7	261	210	×0.8

Table 3 .

 3 Computation timings in µs for point addition and doubling.

			point addition			point doubling
	N	local	shared register	local shared register
	160	2.57	0.78	0.70	1.64	0.50	0.54
	192	3.51	1.01	1.13	2.30	0.58	0.70
	224	4.41	1.95	ror	2.73	1.01	ror
	256	5.89	1.56	ror	3.71	1.09	ror
	384	13.9	7.50	ror	13.3	2.42	ror

Table 4 .

 4 Time [µs] for our GPU implementation and (PACE+mpFq) on CPU.

			point addition			point doubling
	N	our impl. mpfq+pace speedup	our impl. mpfq+pace speedup
	160	0.78	1.52	1.9	0.50	1.99	4.0
	192	1.01	1.91	1.9	0.58	1.99	3.4
	224	1.95	2.65	1.3	1.01	2.69	2.6
	256	1.56	2.65	1.7	1.09	2.65	2.4
	384	7.50	5.11	0.7	2.42	5.01	2.0

Table 5 .

 5 Scalar multiplication result for our GPU implementation and (PACE+mpFq) on CPU.

		computation time in µs	operation throughput in [k]P/s
	N	our impl. mpfq+pace	our impl. mpfq+pace	speedup
	160	179	464	5586	2155	2.6
	192	304	550	3289	1818	1.8
	224	507	878	1972	1138	1.7
	256	617	1003	1620	997	1.6
	384	4609	2941	216	340	0.6

Throughput of 64-bit float is barely equal to 1/8 of 32-bit float throughput.

School book multiplication requires N 2 multiplications and (N -1) 2 additions or N 2 MAD

http://www.open64.net/

Acknowledgments

The authors are grateful to the Nvidia donation program for its support with GPUs cards.