Nicolas Bouleau

IRREGULAR AND SIMULATABLE FUNCTIONALS ON WIENER SPACE in Probability Numeriaques INRIA 1992

We discuss the difficulties of effective statistical simulation by studying the case of functionals of multivariate Brownian motion. This leads us to a proposal of definition of s-functionals which could be a step toward a statement like a church thesis for simulatable random variates. The neologism "simulatable" is taken here in the sense of "•effective for statistical simulation".

I. Different ways of simulation of Brownian functionals

In this part, we discuss the main methods for computing the expectation of Brownian functionals by simulation with a particular interest on a.s. and L p approximation.

A. Almost sure versus L n approximation.

Let us consider for the discussion the case of a stochastic differential equation of Ito's type

ds Jo Jo where x € IR", (Z? ( )t>o is a d-dimensional standard Brownian motion starting at 0, and where a : HT -+ lR n * d and~6 : IFT -* IR" are suitably regular. It is well known (see especially [16] [10] [27]) that if the Frobenius commutativity condition Vr € IR n rr'(x)[(T(x)tl]t> = <r'{x)[cr(x)v}u \/u,v € B. d is fulfilled then a version of the solution can he chosen such that the Ito mapping (ar,w)e lR n xC([0, l],IR d )-» A'*(w)6C([0,l],IR B )

be locally Lipschitz. This occurs in particular when d -1. In that case it is possible to construct almost sure and L p approximations X x ' n (u>) of X x (u>) by approximating the Brownian path u> and using the continuity of the Ito mapping (see e.g. [6] [28] [21]). Actually, with or without the Frobenius assumption, direct discretization schemes X x ' n can be constructed which converge in V and for some schemes almost surely, let us quote among a large literature on this subject [17] [7] [22] [19] [29] [11]. Now there are essentially two ways of using these approximations to compute numerically the expectation of the functional t'{J) = f{X x (u>)) where / is a regular functional : 1. By L p -approximation results, a functional F e is first chosen, defined on a finite (eventually high) dimensional space, such that \\F-F t \\ L i<e and then E[F S ] is computed by simulation that is to say

where the u; n 's are picked out of fi independently, -either by an ergodic transform using Birskhoff's theorem EF e = lim I £ F.(T») "Too yv ^rj (see for instance [START_REF] Arous | Flots et series de Taylor stochastiques[END_REF] [2]) -or using equidistributed sequences. 2. An approximate value of F(u) is computed by taking f{X x,n {u>)) for n sufficiently large where X x,n is an almost sure approximation of the solution of the S.D.E. This procedure is repeated either for independently picked w n 's n > 1 and then averaging or for a sequence of points r n (u) (ergodic theorem) etc. approximation in£" average of direct \i^) computations of F{<J) for a sequence of points ut's" picked in Q computation of lEF e by simulation approximate value of EF

We emphasize the following points : only the second method uses almost sure approximation results.The first one only needs L p approximation or even approximation in law, (see part B below). In the second step of the first method, because the functional F e is chosen to be regular and defined on a finite dimensional space, there exists in general stopping criterion for the computation of JEF t by simulation. On the contrary, in the second method such stopping criteria are not always available. Let us remark, at least, that if, in the second method, we stop at an N which is selected independently of iw, we obtain actually a copy of the first method.

Let us consider a very simple example : Let

F{u) = 1 + (Jo' M*)l<fc) a be defined on C([0,1],IR). More F is continuous and even Lipschitz with ratio 1 on C([0,1],1R) equipped with the uniform norm. Thus the first method applies without difficulty : If where the constant A" can be explicitely computed. Now, a problem with the second method is to write down a stopping criterion which garrantees that G n (u) is near u> using only informations concerning G n (u>). We shall come back to this key point in part II and part IV.

B. About the approximation in law. Especially in the case of Wiener functionals, the approximation in law was largely studied in the literature, it is based on an invariance principle such as the Doob-Donsker theorem or on some of it generalisations to diffusion processes (see [START_REF] Knuth | The complexity of non uniform random number generation, in Algorithms and complexity[END_REF], [START_REF] Milstein | Approximate integration of stochastic differential equations[END_REF], [START_REF] Paley | Fourier transform in the complex domain[END_REF]). For a bounded continuous functional F or Riemann-integrable with respect to Wiener measure, it allows to approximate TEF by the expectation of F(Z.) where {Z t )t<o is a process which approximates weakly the Brownian motion. Just as there are results on the speed of convergence in the central limit theorem (cf. [START_REF] Platen | An approximate method for a class of Ito processes with jump component[END_REF]), estimates of the speed of convergence in the approximation of W,F by JEF(Z.) can be obtained. Thus in the discussion of the preceding paragraph the approximation in law plays a similar role as the approximation in L p : It yields a functional F e (on a different probability space) such that IE/'j is near IEF and can be computed by simulation.

Let us remark nevertheless that the weak convergence can often be improved into stronger convergences. On one hand there is a well known general theorem of Skorohod which gives a framework where it becomes a convergence in probability, but on the other hand in the concrete cases of explicit settings this transformation can be done very naturally in general. So, in finite dimension for example, where the changes of spaces are easy, one has practically never to use weak convergence for computing an expectation.

II. Non-Riemann-integrability of multiple Wiener integrals.

When the Frobenius commutativity condition fails, the solution of SDE's with regular coefficients are not regular. To make the discussion simple we take the case of Jo where w = (u; 1 ,^2) € Q = C([0,1],IR 2 ). (We know, see especially [?], [START_REF] Azencott | Formule de Taylor stochastique et developpement asymptotique d'integrales de Feynmann[END_REF], [START_REF] Alaya | Contribution in this volume[END_REF], that multiple Wiener integrals are generic examples of diffusions at least for SDE's with analytical coefficients). In this case, F has no continuous version (see e.g. [START_REF] Talay | Resolution trajectorielle et analyse numerique des equations differentielles stochastiques[END_REF] and [START_REF] Newton | An asymptotic efficient difference formula for solving Ito stochastic differential equations[END_REF] for stronger results) and it is not difficult to prove that F is discontinuous at every point of C([0,1],IR JO Jo Let K(t) be the increasing process K(t) = /W) 2 + (wj)> Jo then /J WJCL; 2 -/J u^cL;] can be written (cf. [12]) £(A'(1)) where (B t ) t > 0 is a Brownian motion independent of ((u;/) 2 + (w 2 ) 2 ),>o-So thr probability to be evaluated is equal to

IK[/(A'(l))l|M|<r]
with f(t) = lP{B(t) > a) and is strictly positive since / > 0 and IP{|| u ||< r} > 0.

b) Let G be an open set. By density there is a Cameron-Martin function h € G. By Cameron-Martin theorem under P/, (which is equivalent to IP) u> -h is a Brownian motion to which part a) applies, which completes the proof. D Let us quote here the important result of Stroock and Varadhan [START_REF] Platen | Higher weak approximation of Ito diffusions by Markov chains, AdW der DDR IMATH[END_REF] that under suitable regularity assumptions on the coefficient of the SDE, the solution can be redefined to be approximately continuous at points u> belonging to C°°.

By the lemma if <p is a continuous non constant bounded function, the map u> -* y>(/o u>]<&*; 2 ) is not Riemann integrable. This is a fortiori still true when the Wiener space is equipped with weaker topologies such as the topology of pointwise convergence on dyadic points which we shall use later on.

Let us now emphasize that an almost sure approximation scheme is useless for direct simulation unless one has constructed a pointwise stopping criterion. Explicit examples are given in [START_REF] Arous | Flots et series de Taylor stochastiques[END_REF].

III. /^-approximation results.

The preceding discussion shows the importance of L p -approximation results and we quote here for completeness a result which can be used in applications, (see [START_REF] Kushner | Probabilistic methods for approximations in stochastic control and for elliptic equations[END_REF] [32] [START_REF] Sussmann | On the gap between deterministic and stochastic differential equations[END_REF] [24] [START_REF] Doss | Liens entre equations differentielles stochastiques et ordinaires[END_REF] for related results).

Consider the strong solution of the equation A-,=.r + r <T[\'"s)dB. + r b{x" S )ds

Y 0 = •*' Y t = Y tk +<r(Y tk ,t k )(B t -B tk ) + b(Y tk ,t k )(t-t k ) iott e [**,*Ar+i[. Then for p € [1, oof, there is a constant G = C{p, K(T)) such that || sup IK, -X t \ ||"< T a "l\ <€[0,T1
As pointed out by 0. Faure [START_REF] Doss | Liens entre equations differentielles stochastiques et ordinaires[END_REF] this implies that, taking a dyadic partition of [0, T], the corresponding Y^ approximation converges almost surely to X for the uniform norm on [0,T]. A similar result holds for the piecewise linear Euler Scheme Z defined by

Zt = Y ik p±^ + Y tk+l -l-£- for *€[**, W
which allows a implementation by recursive simulation of the Brownian path. By the preceding discussion, the practical use of this fact is, however, questionable in general unless for particular functional as we shall see below.

Remark. What makes problem in the use of almost sure approximation is not that it converges a priori only outside a negligible set (it. can be shown moreover, under suitable regularity assumptions on the coefficients of the EDS, that the Euler scheme converges also in the quasieverywhere sense, that is to say outside a set of zero capacity for the Dirichlet form associated with the Ornstein-Uhlenbeck process on Wiener space [START_REF] Bouleau | On effective computation of expectations in large or infinite dimension[END_REF]). Even if the approximation would converge everywhere the irregularity of the functional puts a stopping criterion problem.

IV. Simulatable functionals.

In this part, we leave the Wiener space for more general frameworks and we draw from the preceding discussion a definition of functionals able to be pointwise simulated. We introduce first the ideas in the case of discrete probabilities.

A. Simulatable functionals on {0,1} N .

Let (x n ) n gN the coordinate mappings of {0,1 } N into each factor. As usual, we introduce the a-fields An = <r{x m m <n) and the product probability IP = (}-6 0 + -SA on Aoo.

Definition 1 . A Random variable from {0,1} N into IR is said to be simple if there exists an {A n )-stopping times T P-almost surely finite such that F be Aj-fneasurable.

Often F has to be supposed moreover bounded (cf. [START_REF] Arous | Flots et series de Taylor stochastiques[END_REF] 1.4 a)). By the fact F = J2^=o ^•l{r=n}P s * the law of F is necessarily discrete :

A-l
It can be shown that an ,4 n -stopping times S and a (,4s)-measurable random variable G with law fi can be constructed such that IE.S be minimal. For x € [0,1], let {x} be the fractional part of x, let us put u(x) -££Lo jsr--Then the minimal value of ES (finite or infinite) is (cf. [START_REF] Ikeda | Stochastic differential equations and diffusion processes[END_REF]) We adopt now the hypothesis that it is possible to pick out of [0,1] a sequence of points according to Lebesgue measure and independently for reasoning on the simulation algorithms as does for example L. Devroye in his book [9].

Let (U n ) n >o, be the coordinate mappings of [0,1] N , F n = o(U m ,m < n), IP = dx® N . We are going to consider random variables F defined on [0,1] N and Ff-measurable for some (F n )stopping time T IP-almost surely finite. But which regularity has to be asked on F to get an interesting definition ?

At least, F has to be supposed IP-Riemann integrable (with [0,1] N equipped with the product topology) such that, outside a negligible set, the value of F at a point x can be approximated by an approximate knowledge of x. However, it is not difficult to define, quite explicitely a Borelian mapping ,\ from [0,1] into [0,1] N such that x De continuous outside a negligible set and that the image of the measure dx be cta® N : If

oo X k = E^T **€{0,1} Jfc=0 2*+i
is the canonical dyadic expansion of x € [0,1], define y = \(x) to be y = (y") n >o with Thus, if F is only supposed to be P-Riemann integrable, it is possible to come back to [0,1] by considering Fox-But, by the serious fact that the digits of y n depend on too faraway digits of x, this is useless practically for the simulation of functionals of random processes 1 .

So, we have to ask more regularity on F than the only Riemann-integrability. The following definition seems reasonable to be taken in consideration :

Definition 2 A random variable F from [0,1] N into
JR, d is said to be simple if there is an (Fn)stopping time T dx**-almost surely finite such that the sets {T = n} have negligible boundaries, F is Tr-measurable and, for every n, F is continuous in the inside of {T = n}.

As before, F has often to be supposed also bounded. It follows from this definition that F is continuous on a set of measure 1 and is therefore Riemann-integrable if it is bounded. T is an (,F n )-stopping time with geometric law, the sets {T -n} have negligible boundaries and X is simple in the sense of definition 2. More generally Definition 3 . Given a probability measure fi on (1R, N ,#(]R. N )) and a Polish space (E,S), we shall say that a measurable mapping G from (1R N ,51R N ) into (E,E) is simple if there is a stopping time S for the filtration (Q n ) of the coordinates (X n ) o/IR. N , (i-almost surely finite, such that the boundaries of the sets {S = n} be negligible and such that G be Q n -measurable and continuous on the inside of {T = n) for every n.

C. Regularity properties of simulatable functionals.

If the stopping time used in the definition is not only almost surely finite but everywhere finite, the laws of simple random variables have the property to depend continuously on the perturbations of the random number generator and analytically if the stopping time possesses an exponential moment. Let us take first the case of {0,1} N and suppose that the generator yields independent bits but with law pSi + (1 -p)6o for p € [0,1] : # p (i) is therefore the sum of a series of holomorphic functions in {z : \z -1| < o} uniformly convergent in every compact of this domain and the result follows.

Proposition 1 . 1) Let F be simple on {0,1 } N with respect to the stopping time T which is finite, then the distribution function H p (t) = JP P {F < t} under probability JP

P = {p8\ + (1 - p)^o)® N is continuous in p € [0,1]. 2) If
• These properties extend to ^V-measurable functionals an in particular to simple functionals on [0,1] N : let h be a bounded Borelian function on [0,1] such that /J h(x)dx = 0, and let us put C =|| h ||oo. Suppose that our generator yield independent points distributed according to the law (1 + Xh{x)).dx for A € [-£, £]. Then

Proposition 2 . 1) Let F be an T?-measurable function from [0,1] N into JR. where T is an {Fn)-stopping time which is almost surely finite for all the probability measures JP\

= [(1 + A.h(x))<ix]® N for A € [-a, a] (a < ^) then the distribution function H\(t) = JPA{-F < *} w continuous in A on [-a,a].
2) Suppose moreover 1E 0 [(1 + ac) T ] < +oo for some a > 0 fAen //\(0 IS analytical in {A:|A|<a}.

The proof is similar to preceding one.

D. Simulatable functionals on a general probability space

We are going to propose a definition of simulatable functionals in general.

Let us consider a Polish space W equipped with its Borelian <7-field W and with a probability m. For example (W, W,m) could be the Wiener space or a canonical space of some stochastic process.

As the preceding discussion shows, the topology on the space W is in general too strong to be used in simulation. Indeed to be effectively used the topology has to be such that a fundamental system of neighbourhoods of a point be composed of sets described by a finite number of rational numbers i.e. by an integer. On the Wiener space, for instance, it could be the simple convergence on rational numbers (cf. part V below).

Then we define a presentation of (W, W, m) as a copy of its measurable structure equipped with a weaker topology : Definition 4 . Let (W,W,m) be a polish space, a presentation of(W,W,m) in a mapping <pofW mR N such that -tp is an isomorphism of measurable spaces between (W, W,m) and (lR N ,5(IR N ),/x) [i.e. fi = <p+m and 3N,m(N) = 0, SN^niN 1 ) = 0 such that <p : W\N -• U^N 1 be one to one and bimeasurable],

• y> is continuous.

Then we need to be able to simulate progressively the coordinate process X = (A" n ) n > 0 of (IR. N , 5(Irt N ), fi) so that if we are given a simple functional G on (IR N , 5(1R N ), fi) we can simulate G too.

Definition 5 . Let X = {X n ) n >o the coordinate o/IR N and (G n ) = a(X m m < n), we shall say that ij>: [0,1] N -• IR N is a progressive simulation of X = {X n ) n > 0 under ft if -rj> is measurable from ([0,1] N ,5[0,1] N ) into (IR N ,S(|IR N |
) and ^(dx N ) = fi -There exists an increasing sequence of (T n )-stopping times Tk,k > 0, dx®**-almost surely finite such that V& the map (X 0 ,..., A'fc)oi/' from [0,1] N into JR. k+l be F^-measurable.

At last, we shall say that a functional F defined on (W,W,m) with values in the Polish space (E,S) is an s-functional if a presentation <p of (W, W,m), a simple functional G from (1R 1N The preceding definitions will be enlightened by their application to the case of Wiener space.

V. The case of Wiener space

Let (W,W,m) be the space of standard Brownian motion vanishing at 0, where W = Co([0,1],IR ), W is the Borelian <x-field on W equipped with the uniform convergence, m is the Wiener measure.

We shall look at several presentations of the Wiener space.

A. Dyadic presentation

Let, as before, X n be the coordinate mappings of 1R N let us put

X 0 o<p(io) = tt'(l)
X 2 n +k o<p(iv) = t»(2ftf) fc = 0,l,...,2 n -l;n>0.

This defines a continuous mapping v? from W into 1R N (equipped with the product topology). Under the probability measure v?*m the coordinates [X n ) constitute a Gaussian process and the law of the vector (Xo,...,Xin.x) is a permutation of the one of (tw(i),... ,iw(£),... ,w(l)). Clearly <p is an isomorphism of measurable spaces.

It is easy here to explicit a pmynssivt simulation ip of the process X = (X n ) on the space ([0,1] N , J5[0, l] N ,ota N ). This consists of simulating AV* AV+i,..., AV+i-i given Xo,..., AV-i, and that can be done by the cclcbraicd recursive definition of the Brownian motion of Paul Levy. We dont write the details.

Thus we have the diagram Proof. It is easily seen that we have to proved the properties of definition 3 for F e and T e when W is equipped with the topology of simple convergence on dyadic points and with the <r-fields

H 2 » = *(w(±)k = l,...,2 n ). a)
Clearly T e is an (?i2 n )-s topping time m-almost surely finite, and F e is Tir,-measurable. Let us put

k k -1 b) Let us put A 2 n(W) = sup | w (--) _ w (_-)|.
F n = E[F|W 2 «] = £ g 2 g (% -u>i. 0% -"fc) of course F n -• F almost surely and in L p p € [1, oof. Hi) E[(F-F n )W = ~E 1 ^~ Jfc=0 fc + 1 * + 6 2»
Then putting V^ for this last expression, let us define

S t = inf{n > 0 : V 2 < e 2 }.
Then by the same argument as in the preceding example Fs, is an s-functional. There obviously holds E[(F -Fs t ) 2 \H2s t ] < e 2 hence taking the expectation \\F-F s ,\\<e.

Then the simulation of h(F e ) to compute the expectation of h(F), for h Lipschitz and bounded say, is a manner to give an effective sense to the second method exposed in the part I. Remark. If, instead of f 0 l u\dui 2 3 we would have taken the Levy area itself L(t) = \ fo(uldu 2 -wjdu}]) it would be easy to find explicitely a continuous random variable with the same law as L(t) on the Wiener space of three dimensional Brownian motion. This is due to the independence property used in the proof of the lemma of part II. Thus the random variable L(t) can be reduced in the sense of [START_REF] Arous | Flots et series de Taylor stochastiques[END_REF] section 1.3.

B.Presentation with general partitions

Instead of choosing dyadic partitions, it is possible to choose all the rational numbers. Let r : IN -» Q D [0,1] be a numbering of the rational numbers E € [0,1] and let us put

X n c*p(io) = w(r(n))
with this presentation it is possible to define s-functionals which depend on the points of a partition of [0,1] whose "deepth" is random as precedingly but now whose "thiness" can vary from place to place in [0,1] in function of the values at points already picked. This family includes the discretization in space for diffusion studied by 0. Faure in this volume. Of course a lot of other presentations of Wiener space can be constructed giving each time its own family of s-functionals. I am gratefull to W. Kendall for comments and very usefull discussion on the concept of effectiveness in simulation. In the oral lecture, s-functionals were called graphic functionals. It is indeed clear now, that if one consider a graphic simulation of, e.g. a Brownian path as on the picture below, this simulation, that is to say the program which yields the picture, defines a map from the Wiener space into the Polish space of piecewise affine lines of the plane which is an s-functional.

  linear dyadic approximation of the order n of u>, and if we take F n (u>) = F(G n (u>)), we have |EF -EF n \ <|| sup Ms) -G»(*)| || 2 < JL

Lemma 1 .

 1 For every open set G C C([0,1], IR 2 ) and a € 1R JP[{F > a} D G) > 0 and P[{F < a} n G] > 0. Proof a) Let us first show that Va € IR, Vr > 0 IP{ f u\du] > a, II w ||< r} > 0. Jo For this, by the relation u*{\)u 2 (\) = f »\<L>] + f JM, Jo Jo it is sufficient to prove Va € IR, Vr > 0 JP{[ 1 u\dw] -C u^du] > a, || u ||< r} > 0.

  a : IR" x IK + -IK" X/ b: lR n x 1R+ -* IR n and Vx,y 6 IR" and s,t < T \<T(X,S) -<r(y,0| + \b(.v,s) -b(y,t)\ < K(T){\x -y\ + \t-s\ a ) for an a > 0, where |.| is an Euclidean norm. Let 0 < t\ < ... < t p < ik+\ < . • • < T be a partition of [0, T] and r = sup \t k +i -t k \, the k Euler scheme Y t is defined by induction on k by

fc=l

  In general however, the laws of the simple random variables to be simulated, are not known. It is precisely for this reason that the simulation is usefull. Examples are numerous : -solving the Dirichlet problem in an open set of JR d for large d by spatial discretization and simulation of a symmetric random walk -pricing of European or American options by discretization of the Black-Scholes model, cf. [8] -modelling of queuing systems, etc. B. Simulatable functionals on [0,1] N .

1

  The mapping \ and its reciprocal are the transformations by fusion and splitting of Paley and Wiener [?] and used in the probabilistic interpretation of quantum mechanics of Wiener and Siegel [?] [?]. The physical meaning of the probabilistic hidden variable of this model was criticized for the same reason as here (see [?] p 149).

  Example. Let X be a random variable with values in JR d with density / continuous (or Riemann integrable) on IR d . If A' is simulated by the rejection method with a random variable Y with continuous density g such that / < kg, and if Y is simulated by Y = (p(x) x € [0,1] with <p continuous, this can be writen T = mf{2n+l:kU 2n +igotp(U 2n )<fo<p(U2n)} and .V = ^(UT-I).

2

 2 One has // p (/) = £Ei[M n l {F < t} l {t=n} ]. n=0 Let us put q = p -4, the expectation El [M n l {F < 0 l {T =n}] = E 11(1 + 2q(2x k -l))l {F <t}l{T=n} LJt=0 is bounded in absolute value by (1 + 2|^|)" +1 IP^[{T = n}].

  ,5(]R ),y>*m) into (E,£) and a progressive simulation xp were found such that F = Gap. {W,W,m) (IR N ,8(IR N ),V>.m) _£ .

(

  W,W,m) (IR N ,5(IR N ),^m) _JL " (E,S)([0,l] N ,S([0,l] N ),dx N )We shall see now that there are simple functional G which define interesting s-functionalsF = Goip. Example 1.For e > 0 let us define 7; = inf{2": sup |«,(1) -to(^i)| < e} and let us denote by Z,2"(u>, 0 the piecewise linear dyadic approximation of order n of w € W. So, w -• L^n(w,t) is an mapping with values in CQ([0,1], 1R. ) equipped with its Borelian sets. This is a Polish space. Let us put Fg = LT, • Then F e is an s-functional : it factorizes in F e = G € cnp and we have Lemma 2 . G s is a simple functional (definition 3).

2 .

 2 Let us remark that the set A € = {w : A Tt (w) < e} is open for the topology of simple convergence on dyadic points. Indeed if wo 6 A e and if the T e (w 0 ) points Nsfr), k = l,...,Tg(w 0 )) are sufficiently closed respectively to the points k ( w o(^r>-r),*s l,...,T e (u;o)) then T e (w) = T e (tw 0 ) and hence Lj t {io) is closed to L,T t (w 0 ). It follows that F e is continuous on the inside of {T e = 2 n } and it. can be proved that the boundary of {T e = 2 n } is negligible, d In the same way the following functional F} = L T} where T e l = inf{2" : 5>(£) -"(^l > \) F 2 = L n where T 2 = inf {2" : £ |w( A) -»(^r )| a 6 [1 -e, 1 + «]} are s-functionals. Example Let us take again the example of part II F(w)= fu\dw]. Jo Let us notice the following facts which are simple to prove. t) E[/>i<KK = (a,,6i) f a» ( = (02,63)] = ^(fc -6 X ) ») E[(J»£ -^(6 2 -6,)) 2 K = (o,,6,),w* = (02,62)] = ^((«2-«.) 2 + (62-6i) 2 ) + ii T li


Nicolas Bouleau

On effective computation of expectations in large or infinite dimension

Keywords: Monte Carlo, quasi-Monte Carlo, Riemann integrable, discrepancy, effectivity, infinite dimension, Bernoulli shift, pointers

Resume

Cette etude est une analyse des difficultes inherentes a l'integration par des methodes de Monte Carlo ou de quasi-Monte Carlo. Contrairement a ce qui est ecrit parfois ces methodes ne sont praticables que dans certains cas. Pour le probleme important du calcul de l'esperance de fonctionnelles de processus stochastiques, on presente les avantages d'une methode de simulation en dimension grande ou infinie fondee sur une implementation de l'operateur de shift na.r des pointeurs.

On one hand simulation methods are often presented for solving stochastic problems, as the only ones which work when no direct computation is possible, but on the other hand they are known to be slow. Some reservations have to be brought to these two statements, first because these methods dont work in any cases but need precise assumptions which are often not satisfied, second because several improvements allow in some cases significant speeding up. The first part of this work is an attempt to make this analysis more accurate to clarify what succeeds and what fails when computing the expectation of a random variable by simulation, i.e. by Monte Carlo or Quasi-Monte Carlo sampling and averaging. For quasi-Monte Carlo methods this depends highly on the dimension of the space on which the random variable is defined. That leads to the concept of simulation morphisms. This is presented in paragraphs 1 and 2. But for both Monte Carlo and quasi-Monte Carlo cases a more fundamental difficulty comes from the fact that some random variables (often even very simple ones) are not Riemann integrable in the setting where they are naturally defined (see examples 2,3 and 4 of §3). Finding an other random variable with the same law and Riemann integrable, is in general a difficult problem which remains open in much cases. As exposed in §4 it is nevertheless possible to compute the expectation of non Riemann integrable random variables without proceeding to such a reduction to the Riemann integrable case, if the random variable can be shown to be effectively in L 1 , that is, can be effectively approximated for the jL 1 -norm by bounded continuous functions. In this part many examples are given and some open questions are set.

The second part is devoted to the infinite dimensional case which is shown to be very usual when computing expectations of functionals of stochastic processes such as Markov chains or diffusions. This is just for introducing the third part in which the shift method is exposed. Using Birkhoff 's pointwise ergodic theorem instead of the law of large numbers is quite a natural idea, but what is more interesting is that, despite the fact that there is no standard speed of convergence for this theorem, this method runs surprisingly fast in usual cases, so much so that it is the only available one in some situations. This is due to the efficient use of pointers that this method allows. As the law of large numbers the shift method has its Monte Carlo and its quasi-Monte Carlo version, the latter being in connexion with the notion of rapidely normal numbers or sequences which replace here low discrepancy sequences. The shift method was exposed in a course of the author at Universite Paris 6 in 1987-88, and in two lectures in European universities in 1988. At present time, this method asks several questions to fundamental research especially on inequalities for the speed of convergence and whether random sequences are the fastest or not in infinite dimension, a result of Flajolet, Kirschenhofer and Tichy [START_REF] Flajolet | Deviations from uniformity in random strings[END_REF] showing that, in a certain sense, they are not far to be the best ones.

The word effective in the title and along the paper is only a preformalized and naive version of recursive predicate theory.

I Difficulties of simulation 1.1 Fast integration in small dimension

Our study being mainly devoted to other cases (integration in large or infinite dimension) we recall this case just for completeness.

Here X is a real random variable defined on the space ([0,1]*, #([0,1]*), dx) and s is a small natural integer. Several papers tends to made the word "small" more precise (cf for instance Warnock [START_REF] Warnock | Application of number theory to numerical analysis[END_REF] or Sarkar and Prasad [START_REF] Sarkar | A comparative study of pseudo and quasi random sequences for the solution of integral equations[END_REF]), (roughly, s < 15 say). More over X is supposed to be Riemann integrable i.e. can be inserted between two continuous functions whose integrals do not differ from much than e, for every e > 0.

In this case which is quite usual for applications for engineers there is better to do than to pick randomly independently points U n in [0, l] s and to compute EX by the law of large numbers:

1 a EX = lim -V X(U n ) dx -a.s.
n=l that is to say, in practice, by using pseudo-random sequences. Indeed, if f = (£ n ) n eN* is a uniformly distributed sequence on [0, l] s , i.e. if the measure fx^ = jj ^2n=i ^U converges narrowly to the Lebesgue measure dx, it is well known, and easy to check, that X being Riemann integrable, it holds :

/ X dfijy -• / X dx.
But for some sequences this convergence is faster than with independent uniform random samples. The discrepancy being defined by

£>*({) = sup |// n ([0,Z!] x ... x [0,x s ]) -x 1 ...x s \ y
for independent samples it cannot be asymptotically better than J °s2 °6n , by a form of the iterated logarithm law, while several sequences are known for which it is of the order c ' ogn ' . then£>;(0 < iffi;! 1 ^ as soon as n > r s _ x .

cf Hammersley [START_REF] Hammersley | MonteCarlo methods for solving multivariate problems[END_REF] c) Irrational translations of the torus. Let a t G (0,1), i = l,...,s, be such that 1, ax,..., a s be independent over Q, then the sequence £ n = ({nai},...,{nor,})

(where {x} = x mod 1 ) satisfies 1 Ve>0, 3c(f,e), V7V J\T({) < <#, e) TV 1 -'

cf Schmidt [START_REF] Schmidt | Simultaneous approximation to agebraic numbers by rationals[END_REF] With such sequences, moreover, the stopping criterions are deterministic when using the discrepancy with the Koksma-Hlawka inequality

\mx-Jxd» N \<v(x)D* N (()

cf [START_REF] Koksma | Een algemeene stelling uit de theorie der gelikmatige Verdeeling modulo 1[END_REF], [START_REF] Hlawka | Funktionen von beschrankter variation in der theorie Gleichverteilung[END_REF] where V(x) is the variation of X, in the sense of Hardy-Krauze (what is modulo some technical points related to the mass on the boundary of the cube, the total mass of the measure v if it exists such that ^([0, xi] x ... x [0,:r s ]) = X(x) -IEX where x = (aj 1? ..., x s ) £ [0, l] s .

Other sequences are known (cf for instance [START_REF] Loo | Application of number theory to numerical analysis[END_REF]) satisfying the estimate

ND* N (t) < C s (log N) s + 0((log N)*-1 ) N>2 (1)
some with constants C s as small as possible such as Faure sequences [START_REF] Faure | Discrepance de suites associees a un systeme de numeration (en dimension s[END_REF] themselves improved by Niederreiter [START_REF] Niederreiter | Low-discrepancy and low-dispersion sequences[END_REF]. But the practical difficulty comes from the fact that in (1) the righthand side increases rapidely with the dimension s (despite the fact that C s itself could be decreasing). That is so true that for s = 20 several millions of iterations are necessary to reach similar results as with random sampling (see [START_REF] Warnock | Application of number theory to numerical analysis[END_REF], [START_REF] Sarkar | A comparative study of pseudo and quasi random sequences for the solution of integral equations[END_REF]).

Remark. It is nevertheless interesting to note that equidistributed sequences allow the use of rejection method. More precisely, if Y is a real random variable which can be simulated by rejection method using a dominating random variable simulated on [0, l] p and if the random variable X whose expectation is to be computed is given by '

X = F(Y,U u ...,U q )
where the U^s are independent uniform variables independent with Y and F is supposed to be Riemann integrable, then 1EX can be computed by an equidistributed sequence on [0, l] p+q , see for instance [START_REF] Bouleau | Probabilites de l'ingenieur, variables aleatoires et simulation[END_REF].

1.2 Changes of dimension.

Let (fii, Ai, Pi) and (f^, >^2 5 P2) be two Polish spaces equipped with their Borel (j-field and probability measures Pi, P2, then if $ is a mapping from Hi into f^ such that i) Pi {a: : $ is not continuous at x} = 0 ii) $*Pi = P 2 then if / is Riemann integrable from Q2 into IR, / o $ is Riemann integrable over Hi and if £ = (£ n ) is Pi-uniformly distributed

*> 1 A r r \ r dJ p T7£/°*«»)-//°*<flPl= // iV n=l J J hence $(£ n ) is P2-uniformly distributed
x = £ £ «» € {0,1} with V7V £ a n ? 0.

n=l l n=N

Let us define 1 _ n _ X(x) = sup -V a k then X is not Riemann integrable. Indeed, let be This theorem is effective, its proof gives h explicitely when fi is known ,see [START_REF] Bouleau | Probabilites de l'ingenieur, variables aleatoires et simulation[END_REF]. A more general theorem exists for stochastic processes ,see [START_REF] Skorohod | Studies in the theory of random processes[END_REF]. c) These facts lead us to introduce the following notion:

\\y -\x n \\\h < E 3-i=n+l l

Hence to compute WY with the precision e > 0, it is sufficient to chose n large enough so that Y£f+i ^ f anc ^ ^nen ^° compute by simulation E|Jf n | with the accuracy |. rj On the classical probability spaces used in simulation which are compact, lEF can be computed as soon as a sequence of continuous functions f n is effectively known such that ||y -/ n ||^i -> 0.

In practice on [0, l] s , this property of being effectively in L 1 splits into two properties :

i) An effective bound decreasing to zero is known showing that lim f \Y\ dJP = 0 fcToo./|r|>fc ii) For each k a sequence f n ^ of polynomials with rational coefficients is known such that the limit

\im\\f n , k -(YAk)V(-k)\\ L i=0
is proved with effective bound. Remark : If the belonging of Y to L 1 is known only by Fatou's lemma, the variable Y does not need to be effectively in L 1 . In general the almost sure limits of non uniformly integrable positive supermartingales are not effectively in L 1 unless they can be obtained also by an other way. Even the dominated convergence theorem does not yield the f n ,k's-

II Integration of functionals of stochastics processes

In situations related to stochastic processes, one has often to do with integration in infinite dimension. When computing (expectations of functions of) entrance times or entrance points of a stochastic process in a set or in a event, or more general functionals of the sample paths, one is dealing with random variables which can be bounded and Riemann integrable but defined on ([0,1]™, £([0,1] N ), dx®™).

Let us consider as a generic example the case of Markov chains (to which reduces also the case of diffusions by discretization of SDE's) and let us look how the classical simulation method works.

Let X n be a Markov chain defined by

X n+1 = F(X n , n, U n+1 ) ,X 0 = x, ( 2
)
where be uniformly distributed over [0, l] q then (UJ) is uniformly distributed over [0,1]^ and in the narrow sense. With Halton sequences, we get infinite dimensional Halton sequences which where studied by Sobol' [START_REF] Sobol | Onthe convergence of infinite dimensional cubature and simulation of Markov chains[END_REF], For such sequences (UJ) the average 1 N iV n=l converges to Eg for smooth g i.e. for "almost cylindrical" functions g. Nevertheless these sequences are unusable in practice because, even for cylindrical functions, the convergence is two slow (or more exactly two late) in high dimension as we have already seen in part I.

F : IR d x IN x [0, l] k -• IR d is
Ill The Bernoulli shift method.

The infinite or high dimensional case turns out to have its specific features which allow to avoid picking completely different sample paths and doing the average on them.

III.l An efficient method.

Although the pointwise ergodic theorem has no typical speed of convergence (such as the iterated logarithm law) the method, we shall explain now, seems to be the most efficient one at present for computing high dimensional functionals of stochastic processes: The implementation of this algorithm takes great advantage of the notion of pointers (or equivalent device of piles management).

Let us consider the preceding Markov chain

X n +i = F(X ni n,U n ii) ,X 0 = x
where the (U n ) are i.i.d. random variables with law dx on [0,1]. As we saw the process (X n ) n >i can be represented on ([0,1]^, dx® 1 *). For the notations to be simple, it is actually convenient here to represent it on ([0,1]^ , dx® 1 * ). The f/ n 's being the coordinate mappings, n > 1.

Let us take, for instance, an entrance time T, T = inf{n > 1 : X n E A} and a functional G(X T ,T).

The expectation JEG(X t ,T) can be computed by Birkhoff's theorem EG(X r ,r)=liml i £ I G(X j ,T)or a.s. Suppose the sequence to = (U\, U2,...) has been picked out, and the pseudo-random numbers Ui, U 2 ,... have been put in pointers as indicated below : The following sequence 0{UJ) = (C/ 2 , C/3,...) is partially already chosen, either it is long enough or it must be lengthened in what case we have the new scheme : In these pictures the arrows ~» represent computations to do, that is, mainly, to compute

X n +i = F{X n ,n,U n+ i)
but in this computation of F(x, n, y) only the variable x (and n) is new and partial computations depending only on y can be stored in the box of U n +i in such a way that only easy computations remain to be done. This storage has of course to be made during the lengthening of the sequence when it occurs as explained above.

III.2 Speed of convergence

Let us quote the following two results on speed of convergence for the Bernoulli shift ([0,1] N , dx®**), cf Krengel [START_REF] Krengel | Ergodic theorems de Gruyter[END_REF].

Proposition 3 For every sequence (ct n ) n£ f^, a n > 0, a n -> 0, there exists a continuous function f from [0, l]

of o(i).

Nevertheless, a speed of convergence can be obtained if the function whose expectation is to be computed fulfills additional assumptions. Some of these results will be published elsewhere.

III.3 Connexion with normal numbers

Let us replace [0, l] w by {0, l} 1 ^ for simplicity of the discussion.

The best points for integration by the shift method are fast normal numbers. Let us recall that a real number x in [0,1], or a sequence of binary digits, is called normal if the sequence xo9 n is uniformly distributed over {0,1}^ equipped with the measure (|<S 0 + \&\)®^. Explicit normal numbers are known, the Champernowne sequence for instance. There are other known normal sequences which are generalisations of Champernowne's one (see for instance Rauzy [15] [16] or Dumont [START_REF] Dumont | Une modification multiplicative des nombres g-normaux[END_REF]). Some of these constructions make use of (.s, ra)-de Bruijn sequences which are periodic mappings u from IN into {0,... ,m -1} of period m s such that the m s strings of length s (u(n),..., u(n + 5 -1)) n = 0,1,..., m s -1 are all different. Algorithms are known for generating infinite de Bruijn sequences (see Ivanyi [START_REF] Ivanyi | Construction of infinite de Bruijn arrays[END_REF] and its references), that is to say infinite sequences whose initial sections of length m s , considered as loops, are (s,m)-de Bruijn sequences for every s.

The efficiency of these sequences is difficult to evaluate, and is largely an open subject. One of the main difficulty is to summarize the performance of a sequence by a single scale while there is naturally an infinity of discrepancies one on each finite product space.

On the space {0,1}^ Flajolet et al. [START_REF] Flajolet | Deviations from uniformity in random strings[END_REF] give an asymptotic result which shows that in a certain sense, random sequences are not far to be the best ones. Ending remarks. a) What was explained for a Markov chain is of course valid for a large class of stochastic processes or stochastic fields which, in usual cases, can be defined in a natural way on ([0, l]^, B([0, l] 1 ^), dx®^) after discretization if necessary. b) It can happen that the process, whose functional is studied, be itself mixing or simply ergodic. In this case the shift method can be applied to the process itself instead of to the underlying coordinates. This works as soon as a sample path of arbitrary length of the process can be simulated. c) It is well known that pointers are particularly convenient to represent tree structures, for instance binary trees. Then the space which takes the place of the product [0,1]^ is the space [0,1] B or a space of the form A B where A is some probability space and B = {/, r}^' (finite sequences of /'s and r's). Such a space with a suitable extension of Birkhoff's theorem is an interesting device for dealing with branching processes which are of great practical interest. Nevertheless the "linear" above explained shift method works actually for branching processes as well, and after our experiments this "linear" implementation runs even faster. d) The shift method can be slightly improved by using the two-sided Bernoulli shift, which allows to never throw a box away. The idea is to simulate two independent copies of the process, one with the classical right to left shift the other with a left to right shift beginning with the box left by the first one. But the benefit is negligible in high dimension.

n

  Examples a) Halton sequences Let r be an integer > 2. If the r-adic expansion of n is n = a 0 + air + ... + a m r m a,-G {0,1,..., r -1} then $ r (n) is defined by belongs to [0,1], and the Halton sequence is the following £ n = ($ ri (n),...,$ r ,(n)) where ri,..., r s are the s first prime numbers. Then it holds [1] Dl(£) < -TT -, as soon as n > r s n t\ lo & r i b) Hammersley sequences. The second time a simulation is run, it is generally possible to improve it by using finite sequences, that is the case for Halton sequences. If (£m)> 1 < m < n is defined by U = (f,$ ri (m),...,$r,-i(™))

1 n 1 ATheorem 1

 111 = {x : lim -£ a k = -and a\ = 0}, n k=i 2 clearly for Lebesgue measure IP (A) = | and if x G A then X(rc) < 1. But in every neighbourhood of x there is a y with X(y) = 1. rj Example 3. Here is a more probabilistic example. The X n 's being i.i.d. with law f on [-1,1], let S n = X\ -f • • • + X n be the symmetrical random walk. Let us put X = ((sup , , f n , J A 2 ) V 0 VVn>2|V2nloglogny / then X is bounded but not Riemann integrable. I~l b) To see where the difficulty lies, let us recall what can be called the fundamental theorem of simulation For every probability measure \i on IR d with compact support, there is a Riemann integrable random variable H with value in IR defined on ([0,1], dx) with law ft.

  given and where the C/»'s are i.i.d. with law ebon [0,1]*. If the following quantities have to be computed o E(T), JE(Xr), JE[G{X T ,T)] for T stopping time such as T = inf{n : (X n ,n) eA}, o 1EYlk=i C(Xk,Xk+i) for C : IR 2 -> IR representing costs, o or more generally E[H(X 9 )] where H : IR™ -» IR, by (2) X n can actually be defined on the probability space ([0, l] k )^ with Ui being the coordinates, and the problem is to find good double sequences m = (tioi,...,w,i,...) e ([o,i]*) N Uj = (tic,.,...,^,..-!) e ([0,1]^ in such a way that N v 1 E^(^.( u i)) i=i converges rapidely to IE[i7(X # )]. a) Let us take k = 1 for simplicity. Then let us remark at first that if we take Uij = Vi + (j mod 1 where (Vf)*€W i s i-i-d. on [0,1] and (Cj)j>i is a (low-discrepancy) uniformly distributed sequence on [0,1], one gets a sequence (iij) which is not uniformly distributed over ([0,1]) N for dx®^. b) On the contrary if (u{j) is chosen in such a way that for every q E IN u) = (u oji ... i u qj )

1

 1 G(XT,T) £ L 1 , where 0 is the shift operator satisfying \/UJ e [0,1] N * U n o 0(u>) = ^n+1 (w) Vn > The use of pointers can be explained as follows. For applying formula (3) we have to compute G(XT, T) on successive points u>, 0(LO), 0 2 (LO), ... of [0, l] 1 ^ , these points are sequences of points in [0,1] « = ([/,(u;),t/ 2 (u,),...,t4(u,),...) *(w) = (U 2 (u),U 3 (u),...,U M (u),...) Then if the stopping time T is finite , what we suppose here, only a finite length of each of these sequences have to be computed : the sequence LO being computed until T(u) and this is told by the test jtest^ Xk(u) still outside A =» k < T(LJ) X^LO) 6 A for the first time =>• k = T(u)

Figure 1 X

 1 Figure 1X \/\y*X 2 \/\^+ ---v/\/+ X k

Figure 2 X \^v^X 2

 22 Figure 2

  Definition 2 A non Riemann integrable random variable X is said to be reduced if a Riemann integrable random variable Y can be exhibited with the same law as X.Most of almost sure convergence theorems of probability theory give examples of random variables whose reduction is known. The reduction is in general a less difficult operation than to find an analytical expression of the law even if analytical is understood in a large sense (allowing the use of non holomorphic usual functions). Example 4. Let (U n ) be the coordinate mappings of ([0,1] N , dx® 1 *). Let us consider the Markov chain {

	X = But the use of non bounded random variable is impossible in simulation if the position ~\ogU. of the poles are unknown. Even when this position is known, the result will depend on the speed with which the uniformly distributed sequence nears the poles. Example 5. To compute f 0 4| = 2 by simulation, let us use a uniformly distributed 1 limsup -J2-77= = +°°-n iV n =l VSn b) Nevertheless, it is often possible to compute the expectation of non bounded or even of non Riemann integrable functions (without to proceed to a reduction to the Riemann integrable case). Example 6. Let U n be the coordinates of ([0,l] w , dx® 1 *). Let us put » 977. _ i sequence £ = (£ 1 N and Y = limsup \X n \

x n = ;rb((™ + 1K-I + l{t/"<x n _ 1 })

[START_REF] Halton | On the efficiency of certain quasi-random sequences of points in evaluating multidimensional integrals[END_REF]

a* e (0,

[START_REF] Halton | On the efficiency of certain quasi-random sequences of points in evaluating multidimensional integrals[END_REF]

which is a manner of writing Polya's urn.

If T n -cr(Ui,i < n), (x n ) is an T nmartingale in [0,1] and if x^ = limsup n :c n , then x^ is not Riemann integrable. Nevertheless if x 0 = |,

the law of x^ is the Lebesgue measure on [0,1]. Hence x^ is reduced. 1.4 Computation of the expectations of non bounded or non Riemann integrable random variables: random variables effectively in L 1 . a) Non bounded random variables are very commonly used in simulation. For instance the log function on [0,1] is used for simulating the exponential law by n ) over (0,1) and let us consider if j_ If the sequence £ n is changed on a nul set A C IN (i.e. such that jj ]C n =i l^( n ) ~~¥ 0) the new sequence is still uniformly distributed'. Let us take £ p 2 = 4-then n

It can be proved that, Y is not bounded and that (7Al)V (-1) is not Riemann integrable. But the martingale X n converges almost surely and in L 2 and OO 1

  They show that the speed of convergence can be arbitrarily slow or arbitrarily near

	a,
	Proposition 4 For every sequence (c n ), c n > 0, increasing to infinity, with c\ > 2,
	there exists a Borel set A with IP (A) = | on which
	i n-l	-i
	1=0	<-n	Vn

1

^ into IR such that IjESb 1 / ° t -E/| +oo a.s.


member IEEE, H. Farhat member IEEE Y Nasser

J.-F Hélard Senior

Positioning Information Based Technique in Cooperative MIMO-OFDM Systems

Keywords: Positioning information, OFDM, MIMO, cooperative network, iterative receiver

Future Communication networks are tending towards a diverse wireless networking world where the positioning information (PI) could be helpful in different techniques like the dynamic resource allocation. On the other hand, the PI could be widely used for cooperative techniques in the relay and/or routing selection process. In this paper, we propose to use the PI in the selection of the relays and then to apply an efficient double layer distributed space time block code (DLSTBC) scheme between the different relays. Using the amplify and forward (AF) technique, we show that the proposed code is very efficient whatever the transmitted power is. Moreover, we show that the relay selection process based on PI yields very powerful results when compared to the random relay selection (RS) process.

INTRODUCTION

Future communication networks are tending towards a diverse wireless networking world, where scenarios define that the user will be able to attain any service, at any time on effectively any network that is optimized for the application at hand. In actual research works, one of the key parameters in the optimization strategies consists in considering the positioning information (PI) together with a database, such as used in finger printing techniques, to reduce the feedback overhead. Indeed, the actual transmission techniques consume a valuable bandwidth in the feedback from the mobile terminal (MT) to the base station (BS). Moreover, the localization information could be used efficiently in the dynamic resource allocation, in synchronization and for interference cancellation purposes.

In cooperative networks, the PI could be widely used in the relays and routing selection. Also, PI and tracking mechanisms allow for predictions of future geographic node positions, hence help in deciding whether to use an advanced communication technique. Moreover, it could be widely used in some critical MT situations like tunnels, canyon regions as well as indoor areas.

Recently, the cooperative techniques have taken a lot of attention due to their efficiency in different transmission scenarios. The first research results on the relaying networks were obtained in the seventies in [START_REF] Van Der Meulen | Three-terminal communication channels[END_REF][2] [START_REF] Cover | Capacity theorems for the relay channel[END_REF] and their interest was re-launched at the end of nineties [START_REF] Sendonaris | Increasing uplink capacity via user cooperation diversity[END_REF][5] [START_REF] Laneman | An efficient protocol for realizing cooperative diversity in wireless networks[END_REF] with the increase of services demands. It has then launched and activated a large amount of work in this domain. The cooperative techniques have widely used in different domains and their combination with other efficient techniques could be of great interest in the future wireless networks. Among these techniques, the combination of multiple input multiple output (MIMO) scheme with cooperative techniques, called hereafter COOPMIMO, has obtained a lot of attention. Contrarily to the classical techniques, the COOPMIMO system is based on the use of relays between the BS and the MT in such a way that the MIMO scheme is applied in a distributed manner among the different relays. In the literature, some works have been achieved in this area [START_REF] Ibrahim | Cooperative communications with relay-selection: when to cooperate and whom to cooperate with?[END_REF][10] [START_REF] Barbarossa | Cooperative wireless networks based on distributed space-time coding[END_REF]. In [START_REF] Ibrahim | Cooperative communications with relay-selection: when to cooperate and whom to cooperate with?[END_REF] and [START_REF] Liu | Cooperative Communications and Networking[END_REF], it has been proposed a general description of the COOPMIMO techniques with decode-and-forward [START_REF] Ng | Asynchronous space-time cooperative communications in sensor and robotic networks[END_REF] (DF), amplify-andforward [START_REF] Nosratinia | Cooperative communication in wireless networks[END_REF] (AF) algorithms and some analytical computations were given. However, the results with both algorithms were achieved without channel coding. The authors of [START_REF] Barbarossa | Cooperative wireless networks based on distributed space-time coding[END_REF] show how cooperation among nodes of a wireless network can be useful to reduce the overall radiated power necessary to guarantee reliable links among the network nodes. However, to the best knowledge of the authors, no contribution has been presented on COOPMIMO using PI. The reader may also refer to [START_REF] Anghel | Multi-user spacetime coding in cooperative networks[END_REF] [START_REF] Barbarossa | Distributed space-time coding strategies for wideband multihop networks: Regenerative vs. non-regenerative relays[END_REF] for other contributions.

The work presented in this paper has been carried out within the framework of the European FP7-IST-WHERE project. The contribution of this work is multifold. First, we propose PI based technique to select the relays among the competitive relays existing in the network. Cooperating nodes can be specifically assigned fixed, or portable relay nodes, or other user terminals that temporarily form an adhoc point-to-point link. Then, we propose to apply a distributed space-time (ST) encoding scheme between the selected relays to ensure the full COOPMIMO orthogonal frequency division multiplexing (OFDM) system and we propose a generalized transmission model from the BS to the MT using the AF technique. Therefore, we analyze and compare some of the most promising MIMO systems [14][15] found in the literature with a double layer ST (DLST) scheme proposed in our previous work [START_REF] Nasser | System level evaluation of innovative coded MIMO-OFDM systems for broadcasting digital TV[END_REF]. The DLST is selected due to its efficiency face to power imbalance between the transmitting antennas. It will be applied in a distributed way between different relays. This paper is structured as follows. Section 2 describes the PI based relay selection technique. In section 3 we present the COOPMIMO system model using the AF technique. In section 4, we present the receiver system model using iterative detector. The system performance is evaluated in section 4 where conclusions are drawn in section 5.

PI BASED RELAY SELECTION TECHNIQUE

The main problem in cooperative networks is to select the best relay node in such a way that the cooperative technique is efficient in complex situations like tunnels and canyons. Indeed, in real networks, when the MT is moving in urban or suburban areas through deep shadowing regions for a long time, the received signal power is really attenuated and then the connection would be lost. Therefore, the general layout of the movement areas as well as the PI could be of great help in ensuring the continuity of the signal power level and hence establishing the connection. In this work, we consider a network of N fixed relays uniformly distributed in a cell of diameter D meters. We assume that the transmission from the BS to the MT is achieved through a two-hop relay link by using R=2 or R=4 relays among the N competitive relays. We assume that each relay has one transmit/receive antenna where the destination has M=2 receiving antennas. Moreover, we assume that the MT is moving through a trajectory of K meters along the x-axis containing some deep shadowing regions like tunnels and that the PI of the N relays as well as of the MT are available at the BS (see Figure 1). Using the PI, the BS will select (and update each m meters) the R necessary relays to perform the COOPMIMO technique by computing the respective distances between the N relays and the MT each m meters. The R selected relays are chosen according to:

ݎ ൌ i R min ൫݀ ோ ,ெ் ൯ ൌ i R min ൬ට݀ ௌ,ோ ଶ ݀ ௌ,ெ் ଶ െ 2. ݀ ௌ,ெ் . ݀ ௌ,ோ cos ሺߠሻ൰ (1)
where ݀ ௌ,ெ் , ݀ ௌ,ோ and ݀ ோ ,ெ் are respectively the distances BS-MT, BS-relays and relays-MT.

Figure 1-Real transmission scenario scheme Once the choice of the cooperating relays is done, different techniques could then be applied in order to enhance the overall system performance in terms of bit rate, transmitted power, and connectivity. The problem turns out now to find an efficient COOPMIMO technique between the selected relays at each update distance. In the following, we will see how the number of competitive relays in the network as well as the cooperation technique will influence the results.

COOPMIMO SYSTEM MODEL

Transmission model

At the relay nodes, we assume that the AF technique is applied. In this technique, the relay amplifies and then retransmits to the destination the data it receives from the source. We assume that the transmission between different nodes is achieved by using the coded OFDM waveform with N s orthogonal sub-carriers. Figure 2 depicts the twohop relaying scheme used in our study. At the source node, the information bits b k are first channel encoded using a convolutional encoder of rate R c , randomly interleaved, and fed to a quadrature amplitude modulation (QAM) module. In this work, the transmission is modeled in frequency domain where ܺሾ݊ሿ is the complex transmitted symbol on the n th sub-carrier obtained at the output of the QAM mapper. We assume in this paper that the source symbols are normalized i.e. Εሼ|ܺሾ݊ሿ| ଶ ሽ ൌ 1 where E denotes the expectation operator.

The cooperation is done in two phases as follows (see Figure 2). Once the relays are selected based on the PI technique, the symbols are broadcasted to the MT and the selected relays. In phase 2, the selected relays retransmit the data to the MT. The update of the selection process will be done every m meters. That is, in phase 1 and at the receiving nodes (relays and MT), the frequency-domain received complex symbols at each relay node r are given by

ܻ ൌ ඥܲ ܪ ௦, ܺ ܸ (2)
and the received complex symbols, in phase 1, at each receiving antenna of the MT are given by:

ܷ ௦,ௗ ൌ ඥܲ ௗ ܪ ௦,ௗ ܺ ܸ ௦,ௗ (3)
where ܲ ௗ is the received power by the destination (MT), ܪ ௦,ௗ and ܪ ௦, are respectively the frequency channel coefficients between the source(BS) and the MT, and between the BS and the relay. ܸ ௦,ௗ and ܸ denote the additive white Gaussian noise (AWGN) at the MT and relay node respectively with variance ߪ ଶ . Taking into account the transmission path loss, the received powers depending on the position, expressed in dB, are given by:

ܲ ൌ 10݈01݃ሺܲ ௦ ሻ-32.4 10ߙlog ሺ݀ ௌ, ሻ ܲ ௗ ൌ 10݈01݃ሺܲ ௦ ሻ െ ൫െ32.4 10ߙlog ሺ݀ ௌ,ெ் ሻ൯ (4)
where ܲ ௦ is the transmitted power by the BS and ߙ is the propagation constant which depends on the transmission environment.

Figure 2-Two-hop COOPMIMO network

In phase 2, the R selected relays will retransmit the received data from the BS to the MT. Using the AF technique, we propose in this work to apply a distributed ST encoding scheme between the signals transmitted by the different relays. In our study, we limit the number of selected relays at each instant t to R=2 or 4 relays.

Relay processing

Using the AF technique, the selected relays normalize first the powers of the received data before applying the ST encoding scheme. The normalized received data obtained by (2) becomes: ൌ ܪ

ܻ ൌ 1 ඥܲ|ܪ | ଶ ߪ ଶ ൫√ܲܪ ܺ ܸ ൯ (

ܺ ܪ ܸ

Once the power normalized, each relay uses a ST encoder which takes Q data complex symbols and transforms them to a (R,T) output matrix according to the ST block coding (STBC) scheme. The STBC coding rate is then defined by L=Q/T. In this process, each group ܡ ܚ ൌ ሾܻ ଵ, , … ܻ ொ, ሿ (where ܻ , verifies (5)) of Q complex symbols at each relay is fed to the ST encoder. Let W r =[w i,t,r ] where w i,t,r (i=1,…,R; t=1,…,T) be the matrix corresponding to the output of the ST encoder at each relay r. In order to achieve the distributed MIMO transmission scheme, each relay r selects the r th vector of the encoding matrix W r and transmits it to the destination. The set of the R vectors transmitted by all the relays will be modeled by the matrix W having the same dimensions of W r , i.e. (R,T). In this model, the COOPMIMO system could be seen as a conventional MIMO scheme with the encoding matrix W. Moreover, in order to have a fair analysis between the different ST coding schemes, the output power of the ST encoder is normalized by R.

Let us now describe the transmission model between the relays and the destination independently of the ST scheme.

Let x=[X 1 ,…,X Q ] of Q complex symbols be the vector of the symbols transmitted by the source used to generate the vector ܡ ܚ ൌ ൣܻ ଵ, , … ܻ , , … , ܻ ொ, ൧ at the input of the ST encoder at the different relays. The relationship between x and ܡ ܚ can be deduced from (2) by:

ܡ ܚ ൌ ۶۾ ܚ ܠ ܞ ܚ r=1,..., R (6)
where ۶ ܚ is a (Q,Q) diagonal matrix composed of the channel coefficients between the source and each relay r and ۾ is the (Q,Q) matrix of the powers transmitted by the source (assumed identical for all symbols). By introducing an equivalent receive matrix at the second phase ܇ ܌ א ԧ ெൈ் whose elements are the complex received symbols from the relays, we can write the received signal on the receiving antennas at the destination as: ܇ ܌ ൌ ۶ ܌,ܚ ۾ ܌,ܚ ܅ ܄ ܌ [START_REF] Ng | Asynchronous space-time cooperative communications in sensor and robotic networks[END_REF] where H r,d is the (M,R) channel matrix whose components are the coefficients H j,r between the relays and the destination, W is a (R,T) complex matrix containing the transmitted symbols w i,t,r . ۾ ܌,ܚ is the (R,R) matrix of the powers received by the MT from the different relays and ܄ ܌ is a (M,T) complex matrix corresponding to the AWGN.

We separate the real and imaginary parts of the entries ܡ ܚ in (6), of the outputs W r of the ST encoder (of each relay r) as well as those of the channel matrix ۶ ܌,ܚ and the received signal ܇ ܌ . Let ܻ ,,Ը and ܻ ,,Ա be the real and imaginary parts of ܻ , . The main parameters of the ST code are given by its dispersion matrices U q and V q (q=1,…,Q), corresponding (not equal) respectively to the real and imaginary parts of W r . With these notations, W r is given by:

܅ ܚ ൌ ܻ ,,Ը ܃ ܙ ܻ ,,Ա ܄ ܙ ିۿ ܙୀ (8)
We separate the real and imaginary parts of , ܡ ܚ , ܇ ܌ and ܅ and stack them row-wise in vectors of dimensions (2Q,1), (2Q,1), (2MT,1) and (2RT,1) respectively. We obtain:

ܠ ൌ ൣܺ ଵ,Ը , ܺ ଵ,Ա , … , ܺ ொ,Ը , ܺ ொ,Ա ൧ ܂ ܡ ܚ ൌ ൣܻ ଵ,,Ը , ܻ ଵ,,Ա , … , ܻ ொ,,Ը , ܻ ொ,,Ա ൧ ܂ ܡ ܌ ൌ ൣܻ ଵ,Ը , ܻ ଵ,Ա , … , ܻ ்,Ը , ܻ ்,Ա , . . ܻ ெ்,Ը , ܻ ெ்,Ա ൧ ܂ ܟ ൌ ൣݓ ሺ1,1ሻ,Ը , ݓ ሺ1,1ሻ,Ա , … , ݓ ሺ2ܴ,ܶሻ,Ը , ݓ ሺ2ܴܶ,1ሻ,Ա ൧ ܂ (9)
where T holds for matrix transpose. Since, we use linear ST coding, vector ܟ can be written as:

ܟ ൌ ۴. ܡ (10)
where ܡ ൌ ሾܡ , … , ܡ ܚ , … , ܡ ܀ ሿ ܂ and F has the dimensions (2RT, 2RQ) and is obtained through the dispersion matrices of the real and imaginary parts of s.

The problem now is to achieve the relationship between the transmitted vector by the source, i.e. x, and the received vector by the destination ܡ ܌ . As we change the formulation of x, ܡ ܚ , ܇ ܌ and ܟ in ( 9), it can be shown that vectors x and y d are related through the matrix G and ۶ ܙ܍ of dimensions (2MT, 2RT) and (2RT, 2Q) such that:

ܡ ܌ ൌ ۵۰۴۶ ܙ܍ ܠ ۵۴۶ ܚܗܖ ܄ ܚ ܄ ܌ = ۵ ܙ܍ ܠ ܄ (11)
where matrix G is composed of blocks G j,r (j=1,…,M; r=1,…,R) each having (2T,2T) elements given by:

) 2 , 2 ( , , , , , , , , , , , , , , , , , , , , , , ,

H H H H H H H H H H H H

                          - - - = ℜ ℑ ℑ ℜ ℜ ℑ ℑ ℜ ℜ ℑ ℑ ℜ L O O M O L M L M L L r j, G (12)
where ܪ ,,Ը (respectively ܪ ,,Ա ) is the frequency channel coefficient between the r th relay and the j th receiving antenna. These channel coefficients are assumed to be invariant during T OFDM symbols. ۰ is a (2RT, 2RT) diagonal matrix corresponding to the powers transmitted by the relays. ۵ ܙ܍ is the equivalent channel matrix, of dimensions (2MT, 2Q), between the transmitted data vector by the source x and the received vector by the destination through the relays. ۶ ܙ܍ and ۶ ܚܗܖ of dimensions (2RQ, 2RQ) can be written in the same way of G. They are easily obtained through (5).

RECEIVER MODEL

Receiver processing: ST decoder

The detection problem at the second phase is to find the transmitted data x given the vector y d and to combine the estimated data at the second phase with that estimated data at the phase 1 in such a way to maximise the combining ratio. In our work, we combine the estimated values ܺ ,ଵ and ܺ ,ଶ as follows:

ܺ ൌ ߩ ଵ ܺ ,ଵ ߩ ଶ ܺ ,ଶ ߩ ଵ ߩ ଶ (13)
where ߩ ଵ and ߩ ଶ are the respective estimated SNR at each phase and ܺ ,ଵ (resp. ܺ ,ଶ ) is the estimated data at the first phase (resp. second phase). The question is now how to estimate the data at the second phase.

In order to estimate the received symbols in phase 2 and in the case of orthogonal STBC (OSTBC), the optimal receiver is made of a concatenation of ST decoder and channel decoder modules. In non-orthogonal STBC (NO-STBC) schemes, a sub-optimal solution composed of an iterative detector is proposed here as shown in Figure 3. At the first iteration, the demapper takes the estimated symbols ŝ , the knowledge of the channel G eq and of the noise variance, and computes the log likelihood ratio (LLR) values of each of the coded bits transmitted per channel use.

The estimated symbols ŝ are obtained via minimum mean square error (MMSE) filtering according to:

ܺ ,ଶ ሺଵሻ ൌ ܘ ܂ ൫۵ ܙ܍ . ۵ ܙ܍ ܂ ો ۷൯ ି ܡ ܌ (14)
where ܘ ܂ of dimension (2MT, 1) is the p th column of G eq (1≤ p≤ 2Q). ܺ ,ଶ ሺଵሻ is the estimation of the real part (p odd) or imaginary part (p even) of X q (1≤ q≤ Q) at the second phase. This estimation will be combined with the estimation of the transmitted symbol at the first phase according to [START_REF] Barbarossa | Distributed space-time coding strategies for wideband multihop networks: Regenerative vs. non-regenerative relays[END_REF].

Once the estimation of the different symbols X q is achieved by the soft mapper at the first iteration, we use this estimation for the next iterations process.

From the second iteration, we perform parallel interference cancellation (PIC) operation followed by a simple inverse filtering:

ܡ ො ܘ ൌ ܡ ܌ െ ۵ ܘ,ܙ܍ ܆ ෩ ሺℓିଵሻ ܆ ෩ ሺℓሻ ൌ 1 ܘ ܂ ۾ ܘ ܂ ܡ ො ܘ (15)
where eq,p G of dimension (2MT, 2Q-1) is the matrix eq G with its p th column removed, ܆ ෩ ,ଶ ሺℓିଵሻ of dimension (2Q-1, 1)

is the vector ܺ ෨ estimated by the soft mapper at the previous iteration with its p th entry removed.

Considered ST coding schemes

The aim of this section is to judiciously build the ST code so that the resulting MIMO scheme behaves efficiently in different transmission scenarios. In [START_REF] Nasser | System level evaluation of innovative coded MIMO-OFDM systems for broadcasting digital TV[END_REF], we have proposed a DLST code for broadcasting systems which is robust in unequal received power scenarios. We then exploit this conclusion to study the DLST in the COOPMIMO systems. In the sequel, we will first consider the well-known orthogonal Alamouti ST coding scheme [START_REF] Alamouti | A simple transmit diversity technique for wireless communications[END_REF] for its robustness and its simplicity. For NO schemes and R=2, we consider the full rate and the fully diverse Golden code [START_REF] Belfiore | The golden code: a 2 × 2 full-rate space-time code with non vanishing determinants[END_REF]. For R= 4 relays, we implement the DLST proposed in [START_REF] Nasser | System level evaluation of innovative coded MIMO-OFDM systems for broadcasting digital TV[END_REF]. This DLST is robust against unbalanced received power and is constructed by using the Alamouti code at the first layer and the Golden code at the second layer. It is given by: [START_REF] Nasser | System level evaluation of innovative coded MIMO-OFDM systems for broadcasting digital TV[END_REF] where

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )           
ߠ ൌ ଵା√ହ ଶ , ߠ ҧ ൌ 1 െ ߠ, ߙ ൌ 1 ݆ሺ1 െ ߠሻ, ߙ ഥ ൌ 1 ݆ሺ1 െ ߠ ҧ ሻ.

SYSTEM PERFORMANCE

In this section, we validate the proposed technique by simulation results. The simulation parameters considered in this paper are summarized in Table 1. In the simulations obtained hereafter, we consider that the links between the different nodes are achieved through a Rayleigh i.i.d frequency channel model. The cooperation is achieved with AF technique and various transmitted powers. The results are obtained with a spectral efficiency η= 4 b/s/Hz. These values are obtained as given in Table 2. We assume that the MT is traveling through a trajectory of 1 Km having a tunnel of 200 m. The propagation loss in the tunnel is characterized by two propagation regions as de described in [START_REF] Zhang | Novel model for propagation loss prediction in tunnels[END_REF]. The first region of 50 meters has a power loss of 15 dB. The second one is generally modeled as a wave guide having a power degradation of roughly 0.06 dB/meter. During its movement, we assume that the MT updates the selected relays according to (1) each 10 m. 5 give the measured bit error rate (BER) at the destination with respect to the number of competitive relays existing in the network using respectively two values of transmitted powers P s = 0dB and P s = 2dB. Figure 4 shows that the PI based technique gives better performance than the random relay selection (RS) based technique whatever the number of competitive relays is. Moreover, the last conclusion could be verified indirectly as well. Indeed, when the number of competitive relays increases, the MT could choose at each update time the closest R relay nodes and hence allows an increased received power. The system converges when the number of available relay nodes in the network is greater than 100. Eventually, we show that the proposed code DLST is very efficient when compared with Alamouti and Golden codes since it is designed to be robust face to unequal received powers [START_REF] Nasser | System level evaluation of innovative coded MIMO-OFDM systems for broadcasting digital TV[END_REF]. Figure 5 gives the same conclusions as in Figure 4 and shows the efficiency of our techniques. 6. CONCLUSION In this paper, we have proposed a generalized framework of the COOPMIMO system based on PI using the AF technique. We have shown that the PI information is powerful in the relay selection process. We have also proposed a very efficient double layer space time code to apply in a distributed way among the relays. This performance is verified in different scenarios whatever the transmitted power is. The DLST code is then very promising candidate for cooperation in OFDM systems. Moreover, its BER value decreases when it is used with PI based technique.

PIC detector

Figure 3 -

 3 Figure 3-Iterative receiver structure

Figure 4 -Figure 5 -

 45 Figure 4-BER measured at the destination, transmitted power P s = 0dB, η= 4 b/s/Hz

Source Relay 1 Destination H r H r,d H s,d STBC OFDM AF STBC OFDM Relay 2

  [START_REF] Gupta | The capacity of Wireless Networks[END_REF]

	y
	θ
	x
	MT at
	instant t 1
	MT at
	instant t 2
	relay

Table 1 -

 1 Simulations Parameters

	Simulation parameter	Value
	Bandwidth (f s =1/T s )	20 MHz
	Rate R c of conv. code	1/2, 2/3
	Polynomial code generator	(133,171) o
	Channel estimation	perfect
	Constellation	QPSK, 16-QAM, 64-QAM

Table 2 -

 2 Different MIMO schemes and efficiencies

	Spectral Efficiency	ST scheme	ST rate L	Constellation R
		Alamouti 1	64-QAM	2/3
	η=4 b/s/Hz	Golden	2	16-QAM	1/2
		DLST	2	16-QAM	1/2
	Figure 4 and Figure

ACKNOWLEDGMENTS

This work has been performed in the framework of the ICT project ICT-217033 WHERE, which is partly funded by the European Union.


Shub Blum

Jérôme Durand-Lose
email: durand-lose@univ-orleans.fr

Abstract geometrical computation: beyond the

Keywords: Abstract geometrical computation, Accumulations, Analog computation, BSS model, Signal machine

  L'archive ouverte pluridisciplinaire

Introduction

There is no agreed continuous/analog/R counterpart of the Church-Turing thesis. Relating the numerous models is crucial to understand the differences between their computing capabilities. For example, Bournez et al. [START_REF] Bournez | An analog characterization of elementarily computable functions over the real numbers[END_REF][START_REF] Bournez | The general purpose analog computer and computable analysis are two equivalent paradigms of analog computation[END_REF] have related Moore's recursion theory on R [START_REF] Moore | Recursion theory on the reals and continuous-time computation[END_REF], computable analysis [START_REF] Weihrauch | Introduction to computable analysis[END_REF] and the general purpose analog computer. The aim of the present paper is to relate two analog models of computation. One, abstract geometrical computation (AGC) deals with regular and automatic drawing on the Euclidean plane, while the second, the Blum, Shub and Smale model [START_REF] Blum | Complexity and real computation[END_REF] relies on algebraic computations over R. Bournez [Bou99] has already provided some relations between linear BSS and Piecewise constant derivative systems which also generates Euclidean drawings. In [START_REF] Durand-Lose | Abstract geometrical computation and the linear Blum, Shub and Smale model[END_REF], AGC without accumulation was proved equivalent to the linear BSS model (i.e. BSS restricted to multiplication only by constants) with an unbounded number of variables. In the present paper, the full BSS is related to AGC with accumulation.

Abstract geometrical computation (ACG) arises from the common use in cellular automata (CA) literature of Euclidean settings to explain an observed dynamics or to design a CA for a particular purpose. While CA operate in discrete time over discrete space, Euclidean geometry deals with both continuous time and space. This switch of context is justified by the scaling invariance of CA and relates to our preference and ability for thinking in continuous rather than discrete terms (for example just think how recent and complex is discrete geometry compared to the Euclidean one). Abstract geometrical computation works in a continuous setting: discrete signals/particles are dimensionless points; the local function of CA -computing the next state of a cell according to the states of neighbouring cells-is replaced by collision rules: which signals emerge from a collision of signals. Signals and rules define signal machines (SM).

This recent model, restricted to rational numbers, is able to carry out any (discrete) Turing-computation [START_REF] Durand-Lose | Abstract geometrical computation: Turing-computing ability and undecidability[END_REF]; even with the additional restriction of reversibility and conservativeness [START_REF] Durand-Lose | Reversible conservative rational abstract geometrical computation is turing-universal[END_REF]. With continuous time, comes Zeno effect (infinitely many discrete steps during a finite duration): not only are accumulations possible, but they can be used to decide recursively enumerable problems by using the black-hole scheme [START_REF] Durand-Lose | Abstract geometrical computation for black hole computation (extended abstract)[END_REF][START_REF] Durand-Lose | Abstract geometrical computation 1: embedding black hole computations with rational numbers[END_REF]. Although accumulations can easily be generated, they can hardly be foreseen [START_REF] Durand-Lose | Forcasting black holes in abstract geometrical computation is highly unpredictable[END_REF].

In the Blum, Shub and Smale model (BSS), machines compute over any ring. Roughly speaking, polynomial functions can be performed on variables as well as tests (according to some order) for branching. The dimension of the input is not bounded, a shift operator is provided in order to access any variable (finitely many variables are considered since only finitely many are accessed in finite time).

In [START_REF] Durand-Lose | Abstract geometrical computation and the linear Blum, Shub and Smale model[END_REF], AGC (without accumulation) and linear (multiplying two variables is forbidden) BSS over R with an unbounded number of variables are proved equivalent. This is done through linear real number unlimited register machines (the arguments of [START_REF] Novak | The real number model in numerical analysis[END_REF] for the equivalence of URM and BSS translate to the linear case).

With a reasonable handling of accumulations, restrictions on multiplication can be lifted thus achieving the full BSS computing capability. They are handled in the following way: there is no second (or higher) order accumulation; only finitely many signals leave any accumulation; and one signal appears where an accumulation takes place. Multiplication is embedded in the same context and the same encoding of real numbers as in [START_REF] Durand-Lose | Abstract geometrical computation and the linear Blum, Shub and Smale model[END_REF], so that addition, multiplication by constants and test do not have to be implemented again. Only the basis is recalled: a real number is encoded as the distance between two signals.

Multiplication of two real numbers, x and y, is done by producing the binary extension of y and according to y n adding or not x2 n . With integers, n is increasing from 0, with real numbers, n goes down to -∞ which explains the use of an accumulation. Each iteration divides x by 2, computes the next bit of y and updates the partial product accordingly. All the values are geometrically decreasing to zero and so are the time and space used by an iteration, so that there is an accumulation. The signal left by the accumulation is located exactly at the value of the product xy.

To prove that AGC with accumulation is, as expected, strictly more powerful than the BSS model, it is explained how to implement the square root. Each iteration only uses addition and multiplication by constants.

Since the reader might be more familiar with BSS than with ACG, more care and illustrations are given to ACG. Section 2 provides all the needed definitions. Section 3 recalls basic encoding and geometric constructions for BSS simulation. Section 4 provides the multiplication between variables. Section 5 explains how to build a signal machine able to compute the square root with an accumulation. Conclusion, remarks and perspectives are gathered in Section 6.

Definitions

Abstract geometrical computation. In this model, dimensionless objects are moving on the real axis. When a collision occurs they are replaced according to rules. This is defined by the following machines:

Definition 1 A signal machine with accumulation is defined by (M, S, R, µ a ) where M (meta-signals) is a finite set, S (speeds) a mapping from M to R, R (collision rules) a partial mapping from the subsets of M of cardinality at least two into subsets of M (speeds must differ in both domain and range) and µ a is a meta-signal, the one that comes out of any accumulation.

Each instance of a meta-signal is a signal. The mapping S assigns speeds to meta-signals. They correspond to the inverse slopes of the segments in spacetime diagrams. The collision rules, denoted ρ -→ρ + , define what emerge (ρ + ) from the collision of two or more signals (ρ -). Since R is a mapping, signal machines are deterministic. The extended value set, V , is the union of M and R plus two symbols: one for void, ⊘, and one for accumulation i. A configuration, c, is a total mapping from R to V such that the set

{ x ∈ R | c(x) = ⊘ } is finite.
A signal corresponding to a meta-signal µ at a position x, i.e. c(x) = µ, is moving uniformly with constant speed S(µ). A signal must start (resp. end) in the initial (resp. final) configuration or in a collision. A µ a signal may also start as the result of an accumulation. This corresponds to condition 2 in Def. 2. At a ρ -→ρ + collision signals corresponding to the meta-signals in ρ -(resp. ρ + ) must end (resp. start) and no other signal should be present (condition 3). Condition 4 deals with accumulations, the first line implies that sufficiently close to the accumulation, outside of the light cone, there is nothing but a µ a signal leaving the accumulation. The second line expresses that there is indeed an accumulation (this is not formalized since it would be too ponderous).

Let S min and S max be the minimal and maximal speeds. The causal past, or light-cone, arriving at position x and time t, J -(x, t), is defined by all the positions that might influence the information at (x, t) through signals, formally:

J -(x, t) = { (x ′ , t ′ ) | x -S max (t-t ′ ) ≤ x ′ ≤ x -S min (t-t ′ ) } .
Definition 2 The space-time diagram issued from an initial configuration c 0 and lasting for T , is a mapping c from [0, T ] to configurations (i.e. a mapping from R ×

[0, T ] to V ) such that, ∀(x, t) ∈ R × [0, T ] : 1. ∀t∈[0, T ], { x ∈ R | c t (x) = ⊘ } is finite, 2. if c t (x)=µ then ∃t i , t f ∈[0, T ] with t i <t<t f or 0=t i =t<t f or t i <t=t f =T s.t.: -∀t ′ ∈ (t i , t f ), c t ′ (x + S(µ)(t ′ -t)) = µ , -t i =0 or ( c ti (x i ) = ρ -→ρ + and µ ∈ ρ + ) or ( c ti (x i
)=i and µ=µ a ) where

x i =x + S(µ)(t i -t) , -t f =T or ( c t f (x f ) = ρ -→ρ + and µ ∈ ρ -) where x f =x + S(µ)(t f -t) ; 3. if c t (x)=ρ -→ρ + then ∃ε, 0<ε, ∀t ′ ∈[t-ε, t+ε] ∩ [0, T ], ∀x ′ ∈[x -ε, x + ε], -(x ′ , t ′ ) = (x, t) ⇒ c t ′ (x ′ ) ∈ ρ -∪ρ + ∪ {⊘}, -∀µ∈M , c t ′ (x ′ )=µ ⇔ or µ ∈ ρ -and t ′ < t and x ′ = x + S(µ)(t ′ -t) , µ ∈ ρ + and t < t ′ and x ′ = x + S(µ)(t ′ -t) . 4. if c t (x) = i then -∃ε > 0, ∀t ′ ∈[t-ε, t+ε] ∩ [0, T ], ∀x ′ ∈[x -ε, x + ε], (x ′ , t ′ ) / ∈ J -(x, t) ⇒ or c t ′ (x) = ⊘ and x ′ = x + S(µ a )(t ′ -t) c t ′ (x) = µ a and x ′ = x + S(µ a )(t ′ -t) , -∀ε > 0, there are infinitely many signals in J -(x, t).
On space-time diagrams, the traces of signals are line segments whose directions are defined by (S(.), 1) (1 is the temporal coordinate). Collisions correspond to the extremities of these segments. This definition can easily be extended to T = ∞. In the space-time diagrams, time increases upwards. To simplify, the same name is used for a signal throughout the computation, in fact, there is a different meta-signal for each speed. As a computing device, the input is the initial configuration and the output is the final configuration.

Blum, Shub and Smale model. (The reader is expected to be more familiar with the BSS model than with AGC, so this part is not very detailed.) A BSS machine is described by a (directed) graph. There are five kind of nodes: input, output (start and end of computation), set a variable to the value computed by a polynomial over neighboring variables, branch according to the value of a polynomial, and shift. There is exactly one input node. Outgoing edges indicate the next instruction. Each node has exactly one outgoing edge, except for output (no outgoing edge) and branch (two, the next instruction is chosen by a test of the form 0 ≤ P (x 0 , x 1 , . . . x n )) nodes. The shift operator allows to access any variable by shifting the context (this is useful since there is no indirect addressing).

Real number encoding.

A Real number is encoded as the distance from a ba signal to its pairing val signal. Since signal machines are scaleless, two sca signals whose distance amounts for a scale are provided as depicted on Fig. 1. All numbers use the same scale. For the value 0, the superposition of ba and val is encoded as a single signal nul. This value is never considered in the rest of the paper; the reader is invited to check that it can be easily covered. Geometric constructions. The construction of the multiplication, like the addition and the multiplication by a constant, relies on some basic geometric constructions. In each picture, the slopes of the line segments are indicated. It can easily be checked that it is scale invariant and provides the desired effect. Signals with equal speeds result in parallel segments, the key to many geometric properties. Figure 2(a) shows how a distance can be halved. This is done by starting two signals. The first one is going three times slower than the second one. The time the first one crosses half the way, the second one goes one full way and half way back, so that they meet exactly in the middle. Doubling is done the other way round.

Algorithm

The multiplication of two real numbers x and y uses the (possibly) infinite binary extension of y. This is done in two steps: normalization and then an infinite loop.

The first step starts by considering signs and special cases (i.e. multiplication by zero). Then, while y is greater than 2, x is multiplied by 2 and y divided by 2 (so that the product remains constant).

The second step carries out the multiplication. The binary extension y = y 0 .y 1 y 2 y 3 . . . (the initialization ensures that 0 < y < 2) is generated bit by bit. The underlying formula is:

xy = 0≤i y i x 2 i .
This is computed iteratively with the updating on Table 1. The two last cases correspond to y n = 1 and y n = 0 respectively. It is started with: p 0 = 0 (product), b 0 = 1 (bit test), x 0 = x and 0 < y 0 = y < 2b 0 = 2. The following invariants are satisfied:

-b n = 2 -n , -0 ≤ y n < 2b n , -x n = x2 -n

Initialisation

With our encoding, detecting whether x or y is zero is trivial. Detecting the signs and computing the sign of the product is also very easy. The following only deals with multiplication of positive values, the other cases are generated using absolute values and symmetry if the product is negative.

The above algorithm should be started with 0 < y 0 < 2b 0 = 2. So that there is a loop that multiplies x by 2 and divides y by 2 until y is small enough. This is illustrated on Fig. 3 (the first two space-time diagrams go one above the other).

The algorithm is sequential: signals for base (ba ), bit testers (b ), x (x ) and y (y ) are fix unless set on movement by some bouncing initialization signals (ini). All distances are from ba and are measured according to the "official" scale (not represented). The b signal stays at position 2 for testing the end of the loop (i.e. y < 2). The signal ini comes from the left. If it meets y before b, this means that the loop is finished (Fig. 3(c)). Otherwise it has to half y, to double x and to test again (figures 3(a) and 3(b)). At the end of initialization (which is always achieved in finite time and finitely many collisions), the b at position 2 is set at position 1 (halved as usual). The signal p amounting for p is generated. Since p is 0 at start, it is set on ba (this corresponds to a different meta-signal). And finally, ini turns to mul that handles the second step of the multiplication. Everything is depicted on Fig. 3(c).

Main loop

The multiplication follows the (possibly) infinite loop defined in Table 1. Basically, the things to do are additions, subtractions, divisions by 2 and tests. These are easy implementable inside signal machines. But, since the loop is infinite, a correct accumulation have to be generated. By correct, it is understood as at the right location and there is indeed an accumulation. The second point is not to be underestimated since, for example, if at each iteration some constant distance would have to be crossed, then there would be an infinite duration process but no accumulation.

The algorithm is driven by a mul signal that bounces between p and other signals. First of all, mul has to test to know which case to consider. This is done easily: going away from ba, if b is encountered before y this means that b n < y n otherwise y n < b n (when they are met simultaneously, i.e. they are equal, this leads to the end of the loop at the next iteration).

The simpler case is when y n < b n . There is nothing to do but to halve b n and x n , that is halve the distance from b and x to ba. Signals p and y are left untouched. This is done as depicted on the lower part of Fig. 4(b).

The other case b n < y n is depicted on Fig. 4(a) and on the upper part of Fig. 4(b). The addition of x n to p n is done by moving p to the location of x, meanwhile x is moved on the right by half the distance it has from p. The signal b is also moved on the right, at a distance from the new p that is half the previous distance. The signal y is moved to the right, at a distance from the new p that is equal to its original distance from b. To ensure accumulation, all lengths are geometrically scaled down and all the computations take place by p and are shifted according to the moves of p. This can be seen on Fig. 4. Each iteration leads to halving the distance between ba and both x and b. Since the distance from ba to y is at most twice the distance to b, this ensures that all distances are bounded by some Ms 2 n at the nth iteration. Clearly the duration of an iteration is bounded proportionally to the maximum distance between the signals. Thus it is also bounded by some Mt 2 n at the nth iteration. There is no time lag between two consecutive iterations, so that there is indeed an accumulation.

Square root

In this section, it is briefly argued why it is possible to compute the square root with an accumulation and how. Let a be any positive real number. The construc-tion follows the digit by digit approximation algorithm constructing a sequence b n such that:

b 2 n ≤ a < b n + 1 2 n 2 .
At each stage it should be tested whether b

n + 1 2 n+1 2 ≤ a. If it is the case then b n+1 = b n + 1 2 n+1 otherwise b n+1 = b n . Using the sequences: d n = a -b 2 n , e n = bn
2 n and f n = 1 4 n+1 , the test corresponds to computing the sign of a

-b n + 1 2 n+1 2 = a -b 2 n -bn 2 n -1 4 n+1 = d n -e n -f n .
The updating of all these sequences is shown on Table 2 (plus another helpful sequence g n = 1 2 n+1 ). The only operations used are additions, multiplications by constants and tests. Again the size of the computing part is decreasing geometrically and computation can be done by the signal encoding the value of b n and shifted with it. The algorithm starts by a pretreatment that finds the initial value for n. This value might be negative (e.g. n = -11 for a = 2 20 + 1).

Conclusion

In the present paper, AGC with accumulation is proved to be strictly more powerful than BSS. This is not very surprising because it is already known to decide in finite time any recursively enumerable problem (in the classical discrete setting). Like square root, a lot of functions should be computable with an accumulation, it would be interesting to identify them. Considering √ 2, an accumulation point of a rational signal machine can be irrational.

If the computation is stopped before the accumulation happens, then an approximation is generated. Computable analysis relies on the idea of an infinite approximating sequence both for representing real numbers and for computing (type-2 Turing machine needs an infinite number of iterations to compute a function on real numbers). The next step would be to relate these two models. One problem would be to miniaturize and to ensure the generation of an accumulation. Another one is that computable analysis only provides continuous functions, while in ACG, there is, for example, the sign function which is clearly not continuous. On the other side, Moore's recursion theory allows non continuous functions (even the characteristic function of rational numbers).

There might be many accumulations to simulate BSS, but none is of order two. Another issue is to consider nth order accumulation and connect with infinite Turing machines and ordinals [START_REF] David | A survey of infinite time turing machines[END_REF].

Fig. 1 .

 1 Fig. 1. Encoding: scale and positions of val for values -π, -1.5, 0, √ 2 and e.

  Figures 2(b) and 2(c) show two ways to halve a distance while shifting it. They also work with two signals emitted that change speed or direction at some point and meet at the desired location. Generating simple shifts is done by modifying only one slope ( 2 3 by 1 2 for Fig. 2(b) and 5 3 by 1 for Fig. 2(b)).

Fig. 2 .

 2 Fig. 2. Basic geometric constructions.

  , and xy = p n + xnyn bn . The last invariant is trivially true for n = 0 and preserved by the loop. Since x n yn bn < 2x n and x n = x2 -n , from the last invariant comes that lim n→∞ p n = xy.

Fig. 3 .

 3 Fig. 3. Preparing the data for multiplication.

Fig. 4 .

 4 Fig. 4. Multiplication.

Table 1 .

 1 Potentially infinite loop to compute the product.

		pn+1 xn+1 yn+1 bn+1
	if yn = 0	stop
	else if bn < yn pn + xn xn/2 yn -bn bn/2 else pn xn/2 yn bn/2

Table 2 .

 2 Infinite loop to compute the square root. < dnenfn bn + gn dnenfn en/2 + fn fn/4 gn/2

		bn+1	dn+1	en+1	fn+1 gn+1
	if dn -en -fn = 0 else if 0 else	bn	dn	stop	en/2	fn/4 gn/2