Topological method for symmetric periodic orbits for maps with a reversing symmetry

We present a topological method of obtaining the existence of infinite number of symmetric periodic orbits for systems with reversing symmetry. The method is based on covering relations. We apply the method to a four-dimensional reversible map.


Introduction
The goal of this paper is to present a topological method, which allows to establish, in finite computation, the existence of an infinite number of symmetric periodic orbits for dynamical systems with a reversing symmetry.
The role and impact of time reversing symmetries in dynamical systems has been extensively covered in the literature, see [La,La1,D1,D2] and references given there. In this paper we will restrict ourselves to the discrete time case. Before we outline our results we need first to introduce a few definitions.
Definition 1 An invertible transformation M : Ω −→ Ω is called a symmetry of the discrete dynamical system induced by f : It is easy to see that a symmetry of a discrete dynamical system transforms trajectories into trajectories.
Alternatively, a dynamical system may admit transformations S : Ω → Ω that map trajectories into other trajectories reversing the time direction.
Definition 2 [La1] An invertible transformation S : Ω −→ Ω is called a reversing symmetry for the discrete dynamical system induced by f : Definition 3 Let ϕ : T × Ω → Ω (where T = R or T = Z) be a dynamical system. For each point x 0 ∈ Ω, we define its orbit Γ(x 0 ) by To describe the symmetry properties of orbits, following [GSS,La1] we introduce the notion of isotropy subgroup.
Definition 4 Consider a dynamical system on a phase space Ω with a symmetry (reversing symmetry) S.
If S(Γ(x 0 )) = Γ(x 0 ) for some x 0 ∈ Ω, then we say that the orbit Γ(x 0 ) is S-symmetric. If the symmetric orbit Γ(x 0 ) is periodic, then we will say that x 0 is a symmetric periodic point.
If S is clear from the context, then we will often drop it and speak of a symmetric orbit.
Let P : X → X be a map with a reversing symmetry S. Let Fix (S) = {x | S(x) = x} be the fixed point set for S.
A standard method for studying symmetric periodic orbits for systems with a reversing symmetry is the Fixed Set Iteration (FSI) method [LQ,La1] (also known as DeVogelaere method [DV]). This method is based on the intersections of the iterates of Fix (S) and Fix (P •S) (see [La1,Prop. 1.2.2] for more details). Below we give its simplified version.
Theorem 1 If y 1 ∈ Fix (S) and P k (y 1 ) ∈ Fix (S) for some k ∈ Z + , then y 1 is a symmetric periodic point for P and its principal period divides 2k.
Theorem 1 shows that to obtain a symmetric periodic orbit of period k, we have to examine the k/2-iterate of P , which is an apparent obstacle for obtaining infinite number of S-symmetric periodic orbits with unbounded periods in a finite computation.
One way to overcome this problem is the method proposed by Devaney in [D2]. To describe the Devaney's method, assume that P : R 2n → R 2n is a C 1 -map with an reversing symmetry S and Fix (S) is an n-dimensional manifold. Now, if we have p ∈ Fix (S), which is a fixed (periodic) point for P and the unstable manifold of p, W u (p), intersects transversally Fix (S), then by symmetry W s (p) = S (W u (p)), hence W s (p) intersects Fix (S) transversally at q. Now we apply P k to the neighborhood of q and it is easy to see, that for k large enough there exists a disk D k in Fix (S), such that P k (D k ) intersects Fix (S) transversally in the vicinity of p. Now by Theorem 1 one obtains periodic points of arbitrary high periods.
The method proposed in this paper bears some resemblance with the Devaney's method, as it also is based on some kind of transversality, which is preserved under the iteration of the map. Our method is purely topological and is based on the notion of covering relation (see [GiZ]) as a tool for the propagation of the topological transversality. Here we will informally outline our approach. We say that a cube N P -covers a cube M , denoted by N P =⇒ M , if the image of the cube N under the map P is stretched across the cube M in a topologically nontrivial manner (see Definition 7). These cubes together with the corresponding choices of coordinate systems will be referred to as h-sets (the letter h suggesting the hyperbolic-like directions). Now in each h-set N we can define horizontal and vertical disks (see Definitions 10 and 11). In Section 3 we prove (see Theorem 3) that if we have the chain of covering relations then for any horizontal disk H in N 0 and any vertical disk V in N k , there exists x ∈ H, such that P i (x) ∈ N i for i = 1, . . . , k and P k (x) ∈ V . Now if Fix (S) forms a horizontal disk in N and Fix (S) forms a vertical disk in M , then from Theorem 1 we obtain symmetric periodic point. Observe that if we can build a chain of covering relations (1) linking N and M of arbitrary length, then we will have symmetric periodic points of arbitrary high period. For example it was shown that this happens for suitable 2-dimensional Poincaré maps for the Michelson system arising from the Kuramoto-Sivashinsky PDE [W], the planar restricted three body problem modelling the motion of the Oterma comet in the Sun-Jupiter system [WZ2] or the Henon-Heiles Hamiltonian [AZ]. In the above mentioned applications we had one unstable and one stable directions. In this paper we apply our method to a four-dimensional reversible map to prove the existence of symbolic dynamics and an infinite number of symmetric periodic orbits. The main feature, which makes this example interesting, is the fact that both stable and unstable directions are two-dimensional. The proof is computer assisted, i.e., rigorous numerics is was used to verify assumptions of abstract theorems.
The proposed method was introduced in [W] in the planar case and direct coverings. In this case the proof of the transversality theorem (corresponding to Theorem 3) is very simple and is based on the connectivity argument. Unfortunately this proof cannot be generalized to higher dimension or to include coverings induced by inverse mappings.
The content of this paper can be described as follows. In Section 2 we define topological notions: h-sets, covering relations and backcoverings relations. In Section 3 we prove the main transversality theorem (Theorem 3). In Section 4 describe how it can be applied to reversible dynamical systems in general. In Section 5 we consider an application of our method to a four-dimensional map with an reversing symmetry. In Section 6 we describe how to verify the existence of covering relations by computer.

Topological tools: h-sets and covering relations
In this section we present main topological tools used in this paper. The crucial notion is that of covering relation [GiZ].

h-sets
Notation: For a given norm in R n by B n (c, r) we will denote an open ball of radius r centered at c ∈ R n . When the dimension n is obvious from the context we will drop the subscript n. Let S n (c, r) = ∂B n+1 (c, r), by the symbol S n we will denote S n (0, 1). We set For a given set Z, by int Z, Z, ∂Z we denote the interior, the closure and the boundary of Z, respectively. For the map h : [0, 1] × Z → R n we set h t = h(t, ·). By Id we denote the identity map. For a map f , by dom(f ) we will denote the domain of f . Let f : Ω ⊂ R n → R n be a continuous map, then we will say that X ⊂ dom (f −1 ) if the map f −1 : X → R n is well defined and continuous. For N ⊂ Ω, N -open and c ∈ R n by deg(f, N, c) we denote the local Brouwer degree. For the properties of this notion we refer the reader to [L] (see also Appendix in [GiZ]).
Definition 5 [GiZ, Definition 1] A h-set, N , is the object consisting of the following data Hence a h-set, N , is a product of two closed balls in some coordinate system. The numbers, u(N ) and s(N ), stand for the dimensions of nominally unstable and stable directions, respectively. The subscript c refers to the new coordinates given by homeomorphism c N . Observe that if u(N ) = 0, then N − = ∅ and if s(N ) = 0, then N + = ∅. In the sequel to make notation less cumbersome we will drop the bars in the symbol |N | and we will use N to denote both the h-sets and its support.
Definition 6 [GiZ,Definition 3] Let N be a h-set. We define a h-set N T as follows Observe that N T,+ = N − and N T,− = N + . This operation is useful in the context of inverse maps.
Moreover, we require that Intuitively, N f =⇒ M if f stretches N in the 'nominally unstable' direction, so that its projection onto 'unstable' direction in M covers in topologically nontrivial manner projection of M . In the 'nominally stable' direction N is contracted by f . As a result N is mapped across M in the unstable direction, without touching M + . It is also very helpful to note that the degree w in the covering relation depends only on A |∂Bu(0,1) .
Definition 8 [GiZ,Definition 7] Assume N, M are h-sets, such that u(N ) = u(M ) = u and s(N ) = s(M ) = s. Let g : Ω ⊂ R n → R n . Assume that g −1 : |M | → R n is well defined and continuous. We say that N g,w The following theorem was proved in [GiZ].
Theorem 2 [GiZ,Theorem 9] Assume N i , i = 0, . . . , k, N k = N 0 are h-sets and for each i = 1, . . . , k we have either Then there exists a point x ∈ int |N 0 |, such that Obviously we cannot make any claim about the uniqueness of x in Theorem 2.
Theorem 2 shows that both the direct and the inverse covering relation can be treated on the same footing. This justifies the following definition.
Definition 9 Assume N, M are h-sets and P is a continuous map. We say that N P,w

⇐⇒ M
If one of the two following conditions is satisfied We would like to stress, that the relation N P,w ⇐⇒ M is not symmetric.

The topological transversality theorem
The goal of this section is to state and prove the main topological transversality theorem for chain of covering relations. For this end we need first to define the notions of vertical and horizontal disks in an h-set.
We would like to remark here that the horizontal disk in N can be at the same time also vertical in N . An example of such disk is shown on Fig. 1. In case homotopies used in the definitions of horizontal and vertical disks are different. The existence of such disks, which are both vertical and horizonal will play very important role in our method for detection of an infinite number symmetric periodic orbits for maps with reversal symmetry. Now we are ready to state and prove the main topological transversality theorem. A simplified version of this theorem was given in [W] for the case of one unstable direction and covering relations chain without backcoverings. The argument in [W], which was quite simple and was based on the connectivity only, cannot be carried over to larger number of unstable directions or to the situation when both covering and backcovering relations are present.
Proof: Without any loss generality we can assume that Then We define g i = f −1 i , for those i for which we have the back-covering relation Notice that from the definition of covering relation, it follows immediately that there are u ≥ 0, s ≥ 0, such that u(N i ) = u and s(N i ) = s, for all i = 0, . . . , k.
The idea of the proof is to rewrite our problem as a zero finding problem for a suitable map, then to compute its local Brouwer degree to infer the existence of a solution.
As a tool for keeping track of the occurrences of coverings and backcoverings, we define the map dir : In the case of a direct covering (i.e. dir (i) = 1), the homotopy h i satisfies In the case of a backcovering (i.e. dir (i) = 0), the homotopy h i satisfies Let h t and h z be the homotopies appearing in the definition of the horizontal and vertical disk for b 0 and b e , respectively. It is enough to prove that there exists t ∈ B u (0, 1), z ∈ B s (0, 1) and We will treat (27) as a multidimensional system of equations to be solved. To this end, let us define We define a map F = (F 1 , . . . , F k ) : Π → R (u+s)k as follows: for i = 2, . . . , k − 1 we set With this notation, solving the system (27) is equivalent to solving the equation F (x) = 0 in int Π.
We define a homotopy H = (H 1 , . . . , H k ) : [0, 1] × Π → R (u+s)k as follows. For i = 2, . . . , k − 1 we set Notice that H(0, x) = F (x). The assertion of the theorem is a consequence of the following two lemmas, which will be proved after we complete the current proof.
Proof of Lemma 4: From the homotopy property of the local Brouwer degree (see Appendix in [GiZ]) it is enough to prove that H(λ, x) = 0, for all x ∈ ∂(Π) and λ ∈ [0, 1].
Proof: From the homotopy property of the local degree (see Appendix in [GiZ]), it follows that it is enough to prove that Let us take x = (t, p 1 , q 1 , . . . , p k−1 , q k−1 , z) ∈ ∂Π. One of the following conditions holds true t ∈ S u , z ∈ S s , p i ∈ S u , for some i = 2, . . . , k − 1 q i ∈ S s , for some i = 2, . . . , k − 1.
Now we turn to the computation of the degree of C(1, ·). Observe that C(1, ·) has the following form: for i = 2, . . . , k − 1 for i = 1p and for i = kp From the product property of the degree (see Appendix in [GiZ]) it follows that In the formula above if dir −1 (j) = ∅ (for j = 0, 1), then the corresponding product is set to be equal to 1. Similarly if u = 0 or s = 0, then the corresponding product is also set equal to 1.
From Collorary 18 in [GiZ] it follows that This finishes the proof.

Reversing symmetry and covering relations.
In this section we apply the tools developed in previous sections to the study of symmetric periodic orbits for reversible maps.
Theorem 7 Let S be a reversing symmetry for the local dynamical system induced by the map P . Assume that Proof: From Theorem 3 it follows that there exists x ∈ Fix (S), such that The assertion now follows from the reversing symmetry of P .

Now we turn our attention to the action of symmetry on h-sets and covering relations.
Definition 12 Let N be an h-set in R n . Let L : R n → R n be a homeomorphism.
We define an h-set L * N as follows We define an h-set L T * N by Informally speaking, L * N is just a natural symmetric image of N and L T * N is the symmetric image of N , but we additionally switch the 'expanding' and 'contracting' directions. We have the following The reversing symmetry maps, in a natural way, horizontal disks to vertical disks and vice versa. Namely, we have the following obvious lemma.
Lemma 9 Let S : R n → R n a homeomorphism, N be an h-set and γ be a horizontal (vertical) disk in N .
Then S(γ) is a vertical (horizontal) disk in S T * N .
In the context of proving the existence of an infinite number of symmetric periodic orbits symmetric h-sets are of special importance.
Definition 13 Let S be a reversing symmetry. We say that an h-set N is It is easy to see, that if N is S T -symmetric h-set, then u(N ) = s(N ) and the dimension of the phase space n = u(N ) + s(N ) must be even.
The following theorem, which is an easy consequence of Theorem 7, illustrates our method of proving of the existence of an infinite number of symmetric periodic orbits.
Proof: From Lemma 9 it follows that in M 0 and M 2 there exist vertical disks contained in Fix (S). We can now apply Theorem 7 to the chain of covering relations M α0 P ⇐⇒ M α1 P ⇐⇒ · · · P ⇐⇒ M αn 5 One four-dimensional reversible example.
In this section we present the application of the method introduced throughout the paper to a four-dimensional reversible map. As a consequence we obtain the existence of chaotic dynamics and the existence of an infinite number of symmetric periodic orbits for a certain iteration of such map. The proof is computer assisted, i.e., rigorous numerics is used to verify assumptions of abstract theorems. The main feature, which makes this example interesting is the fact that both stable and unstable directions are two-dimensional and the map itself is not close to a product of two two-dimensional maps, with unstable dimension each.

An example of four-dimensional reversible map.
Let f : R n → R n , n > 0 be a continuous map, F : R 2n → R 2n be a map defined by Therefore F is a reversible homeomorphism of R 2n . In suitable coordinates we may rewrite F as follows and the reversing symmetry S will be given by S(x, y) = (−x, y). Let us fix n = 2 and let f : In the remainder of this section we will investigate the map F given by (52), where f is as above. The map S(x 1 , x 2 , y 1 , y 2 ) = (−x 1 , −x 2 , y 1 , y 2 ) is a reversing symmetry of F . It is easy to verify that each solution of F (x 1 , x 2 , y 1 , y 2 ) = (x 1 , x 2 , y 1 , y 2 ) (54) satisfies x 1 = x 2 = 0, hence all fixed points of F are symmetric. Solving of Eq. (54) leads to the following system of equations y 2 1 + (y 2 + 1) 2 = 9 (y 1 + 1) 2 − y 2 2 = 1 describing an intersection of a circle with a hyperbola in four points as was shown in Fig. 2. Fixed points P 1 , P 2 are hyperbolic with real eigenvalues, the fixed point P 3 is hyperbolic with four complex eigenvalues, the point P 4 possesses two complex eigenvalues on the unit circle and two real eigenvalues, i.e., it is of elliptic-hyperbolic type. The fixed points may be exactly computed (for example using Mathematica). However, for our further consideration it is sufficient to use two approximate fixed points, which we will still denote by P i , given by P 2 = (0, 0, 2.199939462565084, −3.0396731015162355).
We will show that in the vicinity of P 1 and P 2 the map F 7 has symbolic dynamics on two symbols and there exist an infinite number of symmetric periodic points of an arbitrary large period. Before we proceed with the statement of main results for F we need to discuss how we represent h-sets.

Representation of h-sets in R n .
To define a h-set N we need to specify a homeomorphism c N of R n and two numbers u(N ) and s(N ) (see Def. 5). Since we will use the computer in order to verify covering relations, the homeomorphism must be representable by the machine. The simplest case is to take an affine map. We will use the maximum norm on R n , i.e. x = max i |x i | and we treat vectors as columns.
We define a h-set N = h(x, u 1 , . . . , u k , s 1 , . . . , s n−k ) as follows Hence, the h-set N defined above is a parallelepiped centered at x. In the sequel we will work in R 2k with s(N ) = u(N ) = k. In this case we will use also the notation N = h(x, M ), where M is a linear isomorphism of R 2k . The first k columns of M correspond to unstable directions and the last k columns of M correspond to stable directions.

Important h-sets and covering relations between them
As was mentioned before P 1 and P 2 are good numerical approximations to two hyperbolic fixed points with two-dimensional stable and unstable manifolds. We choose u i j , s i j ∈ R 4 , i, j = 1, 2 to be a good numerical approximations of unstable and stable eigenvectors of DF (P i ). Put We define two S T symmetric h-sets centered at P 1 and P 2 by The important remark is that the sets |N 1 | and |N 2 | are disjoint. The projection of |N 1 | and |N 2 | onto (y 1 , y 2 ) coordinates is presented in Fig. 3. The following lemma was proved with a computer assistance.

Lemma 11
The following covering relations hold: The details of the proof will be presented in Section 6. Remark 12 From above lemma and Theorem 2 it follows immediately, that there exists fixed points in N 1 and N 2 . However, we cannot claim that they are unique in N i nor S-symmetric.
Now we will construct a dynamical link between sets N 1 and N 2 . More precisely, we construct a chain of covering relations connecting N 1 and N 2 . For this purpose we look for a point Q 1 in the neighborhood of P 1 on the unstable manifold of P 1 and such that F k (Q 1 ) (here k = 10) is close to P 2 . Then we define additional points Q i by taking some forward iterates of Q 1 .

Lemma 13
The following covering relations hold The details of the proof will be presented in Section 6.
Let us comment briefly about the spatial relations between N i and H j . The set H 1 is close to N 1 , this is the reason for using the same stable and unstable directions in the definitions of these sets. Sets H 2 and H 3 are close to N 2 and as in the previous case we used the same stable and unstable directions for them.
The numerical evidence of the existence of covering relation H 1 Fig. 4.
projected onto unstable directions of (H 2 ) c is outside the unit ball, (right) the set F 4 c ((H 1 ) c ) projected onto stable directions of (H 2 ) c is inside the unit ball.

Chaotic dynamics of F .
Theorem 14 The discrete dynamical system induced by the map F 7 is semiconjugated with the full shift on two symbols, i.e. for an arbitrary (i j ) j∈Z ∈ {1, 2} Z there exists a point x 0 ∈ |N i0 | such that Moreover, if the sequence (i j ) j∈Z is periodic, then the point may be chosen as a periodic point for F 7 with the same principal period.
Proof: From Lemma 13 we obtain that Since the h-sets N 1 and N 2 are symmetric (by their definition), the reversing symmetry property of F implies From Lemma 11 we get Let (i j ) j∈Z ∈ {1, 2} Z be a periodic sequence od symbols, i.e. i j+k = i j for j ∈ Z and a certain k > 0. Let be a periodic sequence of covering relations, where by N ij ⇐⇒ N ij+1 we mean a corresponding sequence (61-64). Now, Theorem 2 implies that there exists a k-periodic point x 0 ∈ |N i0 | for F 7 , such that assertion (60) is satisfied. Let (i j ) j∈Z ∈ {1, 2} Z be a nonperiodic sequence of symbols. We define the periodic sequences Since N i0 is a compact set, we can find a condensation point x 0 ∈ |N i0 | of {x k } k>0 . Obviously, x 0 satisfies assertion (60).

Symmetric periodic points for F .
Theorem 15 There exists an infinite number of symmetric periodic points for F with an arbitrary large principal periods.
The proof of Theorem 15 is a direct consequence of the following lemma.
Then there exists a symmetric periodic point x 0 of F , such that Proof: Recall, that the sets N 1 , N 2 are symmetric and defined by vectors u j i , s j i , i, j = 1, 2 -see (57) and (59). For i = 1, 2 we define a map b i : We define the homotopy h i by It is easy to see, that conditions (11-13) from Definition 10 are satisfied. Namely, we have This proves that b i is a horizontal disk in N i . Let us remind the reader that S(P i ) = P i and S(u i j ) = s i j . Hence we obtain S(b i (p, q)) = S(k i (pu i 1 + qu i 2 + ps i 1 + qs i 2 ) + P i ) = = k i pS(u i 1 ) + qS(u i 2 ) + pS(s i 1 ) + qS(s i 2 ) + S(P i ) = b i (p, q).
where k 1 = 0.012, k 2 = 0.31 are the coefficients used in (58) to define N i . This proves that b i are contained in Fix (S). We will show that b i is also vertical disk in N i . Namely, if follows from Lemma 9 that S(b i ) = b i is a vertical disk in S T * N i = N i .
Since V 0 , V k ∈ {N 1 , N 2 } we conclude that there exists a horizontal disk contained in Fix (S)∩|V 0 | and there exists vertical disk contained in Fix (S)∩|V k |. Now, the assertion is a direct consequence of Theorem 3.
Proof of Theorem 15: The assertion follows from the fact that the sets |N 1 | and |N 2 | are disjoint and we can construct an arbitrary number of sequences satisfying assumptions of Lemma 16, for example for k > 0.
6 How to verify covering relations with a computer assistance.
In this section we discuss some numerical aspects of the verification of covering relations. Let N, M be a h-sets in R n such that u(N ) = u(M ) = u and s(N ) = s(M ) = s and let f : |N | → R n be a continuous. In order to prove that the covering relation N f =⇒ M holds, it is necessary to find the homotopy h : [0, 1] × N c → R u × R s and a map A : R u → R u satisfying conditions (2-5).
Since in our example f c = c M • f • c −1 N is a diffeomorphism we can try to find a homotopy between f c and its derivative computed in the center of set and projected onto unstable directions, i.e., we define where π u : R n → R u is a projection onto first u variables. We require A to be an isomorphism, then we have deg(A, B u (0, 1), 0) = sgn (det A) = ±1. Now we define the homotopy between f c and (p, q) → (A(p), 0) by h(t, p, q) = (1 − t)f c + t(A(p), 0), for (p, q) ∈ B u (0, 1) × B s (0, 1).
Below we describe precise algorithms.
Definition 14 Let U ⊂ R n be a bounded set. We say that G ⊂ 2 R n is a grid of U if 1. G is a finite set 2. U ⊂ G∈G G

each G ∈ G can be represent in a computer
Definition 15 Let U ⊂ R be a bounded set. By (U ) I we denote the interval enclosure of the set, i.e., the set (U ) I is the smallest representable interval containing U or [−∞, ∞] if there is not a representable interval containing U . Let U ⊂ R n be a bounded set. By (U ) I we denote (π 1 (U )) I × · · · × (π n (U )) I where π i is a projection onto i-th variable.
In the algorithms presented bellow all computations are performed in interval arithmetic [Mo].
First we discuss how we check condition (3).
Lemma 17 Assume N , M be a h-sets and f : |N | → R n be such that f c is smooth. Let G 1 be a grid of ∂B u (0, 1) and let G 2 be a grid of B s (0, 1). If Algorithm 1 is called with arguments (G 1 , G 2 ) and returns True then the homotopy defined in (67) satisfies condition (3).
Proof: Let (p, q) ∈ N − c . Since G 1 is a grid of ∂B u (0, 1) and G 2 is a grid of B s (0, 1) then G 1 × G 2 := {G 1 × G 2 | G 1 ∈ G 1 , G 2 ∈ G 2 } is a grid of N − c . Therefore (p, q) ∈ G 1 × G 2 for some G 1 ∈ G 1 , G 2 ∈ G 2 . Since the Algorithm stops and returns True the condition π u ((f c (G 1 × G 2 ) ∪ LX) I ) ⊂ R u \ B u (0, 1) is satisfied, which implies that for t ∈ [0, 1], h(t, p, q) / ∈ M c . Now we discuss how we verify condition (4). The main point of our approach is that it is enough to compute f c (∂N c ).
Lemma 18 Assume N , M be a h-sets, f : R n → R n be such that f c is a diffeomorphism. Let G be a grid of ∂N c . If Algorithm 2 is called with argument G and returns True then the homotopy defined by (67) satisfies condition (4).

Technical data.
The grids used in the numerical proof of Lemma 11 and Lemma 13 always consist of "boxes", i.e., products of representable intervals. The total number of boxes used in the proof is approximately 2.2 · 10 8 . The numerical proof of Lemma 11 and Lemma 13 took approximately 36 minutes on 2.4GHz processor under PLD Linux Distribution. The C++ sources with a short description how to run the program are available at [W1].