CC BY
36Russian Journal of Nonlinear Dynamics, 2022, vol. 18, no. 1, pp. 83-102. Full-texts are available at http://nd.ics.org.ru DOI: 10.20537/nd220106
MATHEMATICAL PROBLEMS OF NONLINEARITY
MSC 2010: 68Q05, 70F45, 37D50, 68Q80, 37E10, 33E05
The Elliptic Integral Machine: A Collision-based Model of Computation
I. A. C. Melnik
In this work we will show how any elliptic integral can be computed by analyzing the asymptotic behavior of idealized mechanical models. Specifically, our results reveal how a set of circular billiard systems computes the canonical set of three elliptic integrals defined by Legendre. We will treat these Newtonian systems as a particular application of the billiard-ball model, a ballistic computer idealized by Eduard Fredkin and Tommaso Toffoli. Initially, we showed how to define the initial conditions in order to encode the computation of a set of integral functions. We then combined our first conclusions with results established in the 18th and 19th centuries mostly by Euler, Lagrange, Legendre and Gauss in developing the theory of integral functions. In this way, we derived collision-based methods to compute elementary functions, integrals functions and mathematical constants. In particular, from the Legendre identity for elliptic integrals, we were able to define a new collision-based method to compute the number n, while an identity demonstrated by Gauss revealed a new method to compute the arithmetic-geometric mean. In order to explore the computational potential of the model, we admitted a hypothetical device that measures the total number of collisions between the balls and the boundary. There is even the possibility that the methods we are about to describe could one day be experimentally applied using optical phenomena, as recent studies indicate that it is possible to implement collision-based computation with solitons.
Keywords: collision-based computing, billiard, physical models of computation, elliptic integral, arithmetic geometric mean
1. Introduction
Computation is a natural process and computers are machines designed to explore physical phenomena in order to accomplish desired computational tasks. An example is the Antikythera mechanism, which was designed by the ancient Greek civilization with the aim of predicting astronomical events, and that is considered the first known analog computer. Much later, in the
Received July 07, 2021 Accepted January 28, 2022
Igor A. C. Melnik igorcancela@gmail.com
Universidade Federal do Pampa
Av. Pedro Anunciagäo, 111, Cagapava do Sul, 96570000, Brazil
early 1800s, Charles Babbage designed a purely mechanical device designed to tabulate logarithms and trigonometric functions. His design originated the concept of a digital programmable computer. Electrical components can also be combined to form a purely electrical computer, a textbook example is the RLC circuit which shares its mathematical description with the massspring system. Some devices explore, simultaneously, the continuum change in electrical and mechanical phenomena. These machines are classified as electromechanical analog computers and examples are spacecraft navigation instruments, torpedo fire-control and analog watches.
In the 1930s, the German engineer Konrad Zuse created the first digital computer, Z1, by combining electric and mechanical components. Its successors, the Z2 and Z3, were improvements made upon the Z1 design by replacing mechanical components by relays. The analog computers are not as versatile as the digital ones. The number of analog components needed varies according to the complexity of the problem that will be computed. However, as Zuse pointed out in [1], noncontinuous processes are not well reproducible with analog instruments, while continuous processes are not perfectly reproducible with digital ones. Despite that, the use of digital-based computers progressed largely, while its computational power increased exponentially in accordance with Moore's law. Hybrid computers combine the qualities of both types of devices, thus, most of today's PCs are digital-based, but work with hybrid components1 and are built with an analog-to-digital conversion system.
Physical phenomena that can be explored to perform computation depend on environmental conditions. Recently, the concept of an almost mechanical-driven rover was proposed to explore Venus [2]. It would be controlled by a wind-powered computer which is designed to support the extreme conditions found in that planet. The project considers that a small set of high temperature electronics will be needed to carry out a specialized series of actions.
The fact that a computer is a scenario where physical phenomena are explored in order to perform computational tasks is not necessarily clear from the perspective of a user of a mathematical theory of computation. A circuit designer, for example, can systematically think in terms of Boolean logic through logical primitives, with any concerns about the physical realiz-ability of the model he designs. Given the physical implementation of a digital logic family, they are treated as axioms by the user of the theory, who does not need to be aware of its physical meaning [3]. In the case of digital computers, the computational model may be specified by the concepts established by the Turing Machine, where an algorithmic procedure must be executed by an automaton in a finite number of steps. Thus, digital computers can be adapted for a general purpose use, but deal with continuous functions only via approximations [4].
Today, the AND and the OR logical gates are macroscopic dissipative devices made of no fewer than 109 to 1010 atoms [5]. Besides attempts to build smaller components and improve hardware designs, new computational paradigms are being considered. Researchers seek to develop novel methods to explore the computational resources offered by nature, and find a vast field of natural phenomena that can be explored in order to perform computational tasks. For a complete summary of results in this field, the reader is referred to the monumental book of Gregorz Rozenberg [6]. Also, some general considerations in order to find a physical system that is suitable to perform a specific computational task through experimental methods were established in [4].
Of particular importance to the present work is the mechanical computer described in the seminal work of Edward Fredkin and Tommaso Toffoli [3]. This conservative model is based on elastic collisions of identical balls and is named billiard ball computer. The authors demonstrated
1 An example is a piezoelectric sensor.
that it can be used to simulate any turing machine, revealing the computational potential of Newton's laws of motion. In Fredkin own words, Quite literally, the functional behavior of a general-purpose digital computer can be reproduced by a perfect gas placed in a suitably shaped container and given appropriate initial conditions [3]. In this billiard ball model (hereafter BBM), the computation is encoded in the initial conditions and is performed by the dynamics of the collisions [7]. A computer must perform three primitive functions: remember, communicate and compute. Fredkin states that, in his model, the particles remember by existing, communicate by moving, and compute by interacting. The Turing completeness of the BBM was explored to define collision-based methods to perform elementary arithmetic operations. In [8] it is described how additions, subtraction and multiplication methods can be embedded in such a system, while in [5] a collision-based division algorithm is presented. The computational power of Newton's laws of motion was discussed in [4, 9, 10] in the context of the Newtonian machine these authors designed. Also, the extensive book of Andrew Adamatzky [11] presents a diversity of models that derived from the BBM. A detailed description of these results is far from the scope of this
In this paper we will describe a mechanical machine designed to compute any elliptic integral through collisional processes. We established theorems and corollaries that specify how the initial conditions of the balls enclosed by a circular boundary encode the computation of the canonical set of the three elliptic integrals defined by Legendre. From an initial set of results, the theory of the elliptic functions was evoked in order to expand the computational capabilities of the model. In this way we built collision-based methods to compute elementary functions and mathematical constants.
Cellular automaton (hereafter CA) is a discrete model of computation first proposed by von Neumann and Stanislaw Ulam. It is specified by simple and local rules that act in a discrete time and space. Any model defined in the BBM can be represented by a CA, where the states represent the ball positions; the updating rules of each cell's neighborhood should be defined in order to encode the collision process and thus support the dynamics of the balls [5]. Details regarding the implementation of ball movements and collisions trough a CA were established by Margolus in [18]. Despite their simplicity, CA can produce complex behavior, and for this reason they can be exploited to model systems in fields as diverse as biological and economic sciences. Several examples of applications of these systems are given in [13]. CA models constitute an alternative to the traditional approach of using partial differential equation to model dynamical systems. The expectation is that, eventually, a physicist will be able to formulate CA rules in a manner analogous to how they derive differential equations [14]. For this, we should be able to establish the initial conditions and the local transition rules of a CA from a given prescribed global situation. This problem is traditionally termed in the literature the inverse problem [13]. However, there has been no work to date that analytically derives a subset of the CA rules in order to display a desired global behavior; moreover, the inverse problem is extremely difficult to solve for general cases [13]. For that reason, CA models are still being used for just a few computing tasks [16]. Furthermore, a classification of CA models according to its computational capabilities, which are defined by their dynamics, was proposed in [15].
In [17] Toffoli argued that the development of mathematical physics was influenced by the nature of the available computational tools at that time, which are based on continuous variables. Following his point of view, the discrete features of a CA are evidence that time and space are themselves discrete. The provocative title of one of his articles, Cellular automata as an alternative to (rather than an approximation of) differential equations in modeling physics, make his view clear. Similar ideas were also proposed by Fredkin [19] and Konrad Zuse [1], and a part
of these discussions is synthesized by the Fredkin finite nature hypothesis in digital philosophy. Fredkin also proposed that cellular automata can model not only observed dynamical systems, but also the laws of physics themselves [19].
Due to its characteristics, CA models can lead to the development of practical computer architectures. These models were implemented in hardware, and as an example we have the cellular automata machine family designed mainly by Toffoli and Margolus at the Massachusetts Institute of Technology [7]. The CAM-6 was even commercially produced. A recent overview of the existing physical implementations of CA is given by [16]. There is also a possibility to implement some classes of CA using photons [5]. In these physical systems, optical soliton collisions occur in photorefractive crystals and fibers [5, 11]. These machines operate in a high degree of parallelism, and thus can simulate models with a performance several orders of magnitude greater than a general-purpose computer programmed to do the same task [20]. A very clear introduction to these systems was given in [7], which Toffoli puts as an introduction and orchestration manual for composers of cellular-automaton universes. He concludes that the invention of useful models that requires CA machines will be a stimulus for their evolution, which we believe is the main contribution of the present work.
2. Collision-based method to compute elliptic integrals
Legendre showed that integrals whose integrand can be expressed as a rational function R ^t, y/P(t)j, where P(t) is a polynomial of degree three or four in t, can be reduced using a canonical set of three special functions. They are elliptic integrals of the first, second and third kind, and in the Legendre canonical forms they are defined, in that order, by
mk)= 1 de, (2.1)
J J\- k2 sin2 0
0
<P
E(</>, k) = J Vl - k2 sin2 e do (2.2) 0
and
11(0, n, k) = [-1 . = d0, (2.3)
0 (1 - n sin2 0)V1 - k2 sin2 0
where k ^ 1, (f) e [0, and n G M. In this work we will present a collision-based method to compute the functions Fk), E(0, k) and n(0, n, k) with a set of conservative dynamical systems. The procedure we elaborated is based on a set of circular billiard tables, each one with a set of point-like balls on it. In this scenario, the computation of the elliptic integrals of the first, second and third type are encoded by the initial conditions of the balls and are performed by elastic collisions. Our first six results condensed the information necessary to build collision-based methods to compute any elliptic integral function, elementary functions and some mathematical constants. They define a logic family that are applied using results of the theory of the elliptic functions to further explore the capabilities of the model.
The method we designed can be viewed as a particular case in which we can compute with Newton's laws of motion. From this perspective, our conclusions are merely manifestations of the
fact that the laws of physics were computationally universal [19]. Other examples of computations performed from Newton's laws of motion were given in [4, 9, 10, 12]. In particular, the method we are about to describe shares some characteristics with the method recently idealized by Galperin [12]. His model received much attention recently, even from nonspecialized media2, and was later pushed further, while an explicit solution for the balls positions and velocities as a function of the number of collisions was obtained in [29], where even some considerations about the experimental realizations of the method are presented.
From another point of view, in this work we created the rules and established the initial conditions of a CA family that computes a set of integral functions defined by the theory of the elliptic functions. We then solved an example of what is traditionally called the inverse problem [13]. Our method establishes a family of CA models that can be used to compute (or as an alternative to [17]) a set of integrals of mathematical physics. In this way, we linked some results given mostly by Gauss, Euler and Lagrange, while they developed the theory of the elliptic functions with the BBM idealized by Fredkin, and thus to the CA model of computation proposed by von Neumann.
2.1. The method
Here we assume that the reader is familiar with the basic definitions of billiards. For a formal mathematical treatment of these systems the reader is referred to [26]. Also, the method we are about to describe assumes that the balls do not interact with each other, and we can avoid this assumption by considering, alternatively, an equivalent set of independent billiard systems, each one having one ball on it.
We will derive analytical expressions for several integral functions in terms of the mean collision rate per particle by analyzing the asymptotic behavior of circular billiard systems.
2.1.1. Main results
The next results show how the initial positions and velocities of a set of point-like balls encode the computation of the integral function F($, k), a class of its integrals, and of the function n(0, n, k).
• Computing F (§, k)
Theorem 1. Consider a circular billiard table of radius R composed of N noninteracting billiard balls. Lets assume that, initially,
• the balls are uniformly distributed along a circle of radius r = kR, with k ^ 1, and that
• the velocity vector of each ball has isotropic orientation and modulus v.
Let /3C fc) be the total number of collisions measured after a time t = m
2R
v
Then, as N —> oo and m oo,
2See, for example, https://wordplay.blogs.nytimes.com/2014/03/10/pi/.
Fig. 1. Panel (a) represents the initial state of the system established by Theorem 1 and by Corollary 1, where the initial position and velocity of a generic ball i is evidenced as a function of 61i and 62i. Panel (b) shows the set x ©2 for the conditions established in Theorem 1 and in Corollary 2
Proof.
The initial conditions of the system defined above are represented in Fig. 1a. It is trivial to show that in a circular billiard table each incident angle, zi, is constant throughout the whole trajectory. Consequently, the collision rate of each ball is constant, and is given by fi = j-, where Li = 2R cos zi is the distance between two successive collisions. Considering a system of N noninteracting balls, it is also trivial to show that, as m ^ <x>, the mean collision rate converges to
ac t n 7 \ N N N
i=1 i=1 ' i=1 '
Dividing both sides of Eq. (2.5) by N, we conclude that the mean collision rate per particle is given by
= (26) Nm N ^ cos z/ v ' 7
On the left-hand side of the equality we have the time average of the collision rate per particle, while on the right-hand side we have the average of a propriety associated to each particle. The second will be written as a space average when a given set of initial conditions is defined. The equivalence between these two parameters is due to the fact that the circular billiard table defines a measure-preserving dynamical system [26]. In particular, we admit the initial conditions established in the theorem. The angles d1i and d2i in Fig. 1a define the initial position and velocity vectors of a generic ball i. The ball lies at a distance r = kR from the center of the circle. From the law of sines we have
sin^ = sin^ = fcsin/^ = fcsin(|6>2i -6>H|), (2.7)
with p = -\e2l — ell\\. Thus,
(2.8)
Nm N
Now define 01 and 02 as the set of all possibles values for eii and d2i, respectively. As specified by the initial conditions, 01 x 62 is the set of all points that lie inside the 2* side square represented in Fig. 1b. Thus, as N ^-<x>, the Riemann sum on the right-hand side of Eq. (2.8) converges to the integral l[ given by
2n 2n
i f 1 = de.de,. (2.9)
\*k2) I I \Ji-k2sm2(e.2~el)
As a consequence of the *-periodicity of sin2(e2 — d1), the integrand of l[ is periodic along some of the triangles evidenced in Fig. 1b. By a careful analysis we can conclude that can also be written as
ii={t~2 )[ [ , 1 = de.de 2, (2.io)
V27F ] I ej2-, \/l "fc2 sin2(02
which, in its one-dimensional form in the variable e = e2 — ei, is related
to F{ f, k) by
0
This concludes the proof. □
• Computing F(0, k)
Corollary 1. Consider the system established in Theorem 1 and assume that, initially,
• the balls are uniformly distributed along a circle of radius r = kR, with k ^ 1,
• the velocity vector of each ball has modulus v, and that
• 01 x 02 e R2 is the set of all points 0 ^ eii ^ 2* and 0 ^ e2i ^ 2* that, given a real number (f> £ [0, , satisfies one of the following conditions:
- * + 0 < e2l — eH;
- 0 < e2i — eii < *;
- —* + 0 < e2i — eii < 0;
- —2* + 0 < e2l — eu < —*•
Let, /3C fc) be defined by Theorem 1 and pc((f), k) be the total number of collisions measured after a time t = m^. Then, as N —> oo and m —> oo,
*> = (m) [r (f ■k) - i1 -1) ^ k\
(2.12)
Fig. 2. Geometrical representation of the set ©t x ©2 f°r the conditions established by Corollaries 1 and 3 Proof.
From the same reasons given in the proof of Theorem 1, we have
im i
In this case B1 x B2 depends on 0 as it is defined as the set of all points that lie simultaneously inside the square of side 2n and outside the dashed region represented in Fig. 2. Thus, as N — to, the Riemann sum on the right-hand side of the last equation converges to a bidimensional integral defined in this region. As a consequence of the n-periodicity of the sin2 d function, IFF can be written in a one-dimensional form in the variable d = d2 — d1 as
TF 1 r 1
¿2=7-IT / / .=<№, (2-14)
(*■-<£) / y/l - k2 sin2 e
which is related to F(0, k) by
¿2F = ¿f (0, k) = (2F k) - F(0, kj) . (2.15)
Now we will show that the initial conditions encode the computation of some integrals of the integral function F(0, k).
k
• Computing / F (§, k) kl+ldk _o_
Corollary 2. Consider a circular billiard table of radius R composed of N noninteracting billiard balls. Let us assume that, initially,
• the balls are radially distributed inside a disc of radius rd = kR, with k ^ 1,
• the number of balls per unit area is a(r) a rl, with l e N and r ^ rd, and that
• the velocity vector of each ball has the same modulus v and has isotropic orientation.
Let (3d k) be the total number of collisions that have occurred after a time t = mThen, as N ^^ and m ^ x>,
k
J F (I} ki+1 dk = 0
Proof.
The angles eii and e2i shown in Fig. 3 define the initial position and velocity vectors of a generic ball i. The ball lies at a distance ri = k^^R ^ rd from the center of the circle, and the incident angle is represented by zi.
*
2(l + 2)Nm
kl+2p d
*
(2.16)
2
Fig. 3. Representation of the initial state of the system established by Corollaries 2 and 3. The initial position and velocity of a generic ball i are evidenced as functions of eii and e2i
From the same reasons given in the proof of Theorem 1, we have
(2.17)
13d (§, k) _ 1 A 1
Nm N
¿=i \J l-^sin2«^-^)'
Considering the radial symmetry established by the corollary and the result established by Theorem 1, it is natural that, as N ^-<x>, the Riemann sum in Eq. (2.17) converges to the integral IF defined by
rd rd
Sl[(k)a(r)dA f F (§, k) a(r) dA % = -^ (|) ^-'
fa(r) dA V 7 fa(r) dA 00
where dA is the area element in polar coordinates (dA = 2*rdr).
Assuming the power-law density profile, it is trivial to conclude that
T F T3
(2.19)
□
• Computing J F(0, k)kl+1 dk _o_
Corollary 3. Consider the system defined in Corollary 2 and assume that, initially,
• the balls are radially distributed inside a disc of radius rd = kR, with k ^ 1,
• the number of balls per unit area is a(r) a rl, with l £ N and r ^ rd,
• the velocity vector of each ball has modulus v, and that
• 01 x ©2 £ R2 is the set of all points 0 ^ O1i ^ 2n and 0 ^ O2i ^ 2n that, given a real number (p £ [0, , satisfies one of the following conditions:
- n + 0 < 02i - du;
- 0 < K - On < n;
- -n + 0 < O2i - Ou < 0;
- -2n + 0 < O2i - Ou < -n.
Let /3d k) be defined by Corollary 2 and fid(4>, k) be the total number of collisions measured after a time t = m ^. Then, as N —> oo and m —> oo,
F (0, k) kl+1 dk =
n
(l + 2) Nm
kl+2
P
dfn A h 0 \ ad
(2.20)
Proof.
From the same reasons given in the proof of Corollary 2, we have
Pd(<t>, k)
Nm
1 N -Y
N
1
i ^i-fc2sin2 {e2l-eu)
(2.21)
In this case 01 x02 depends on 0 as it is defined as the set of all points that lie simultaneously inside the square of side 2n and outside the dashed region represented in Fig. 2.
Considering the radial symmetry established by the corollary and the result established by Corollary 1, it is natural that, as N — to, the Riemann sum on the right-hand side of Eq. (2.21) converges to the integral TF defined by
tf _ 0
/ TF(0, k)a(r) dA
f a(r) dA
(n -
f [2F(f, k) — F((f), k)] a(r) dA
f a(r) dA 0
(2.22)
i
d
0
k
0
r
r
1
4
r
and assuming the power-law density profile, we conclude that
Id
tF_,1 + 2\ / R\'-~
d0
□
Now we will show how the initial conditions encode the computation of the function n(0, n, k).
• Computing II n, k)
Corollary 4. Consider the system established in Theorem 1 and assume that, initially,
• the balls are uniformly distributed along a circle of radius r = kR, with k ^ 1, and that
• the velocity vector of each ball has isotropic orientation and modulus vi(n) = v0[1 — —n sin2(e2i — eii)]"1, with positive v0 e R and negative n e R.
Let, ¡3C n, k) be the total number of collisions measured after a time t = m^. Then, as N ^^ and m ^ x>,
• Computing n(0, n, k)
Corollary 5. Consider the system established in Theorem 1 and assume that, initially,
• the balls are uniformly distributed along a circle of radius r = kR, with k ^ 1,
• the velocity vector of each ball has modulus vi(n) = v0 [1 — n sin2(e2i — eH)] 1, with positive v0 e R, negative n e R, and that
• 01 x 02 e R2 is the set of all points 0 ^ eii ^ 2* and 0 ^ e2i ^ 2* that, given a real number <j) £ [0, satisfies one of the following conditions:
- * + 0 < e2i — eii;
- 0 < e2i — eii < *;
- —* + 0 < e2i — eu < 0;
- —2* + 0 < e2l — eu < —*.
Let, /3C n, k) be defined by Corollary 4 and pc((f), n, k) be the total number of collisions measured after a time t = m^. Then, as N —> oo and m —> oo,
w, n, k)=y (i'k) - -1)?^ n>k)
(2.25)
The proofs of Corollaries 4 and 5 are analogous to the proofs of Theorem 1 and Corollary 1, and for that reason they will be omitted.
Figure 4 shows a geometrical representation of the trajectories in the systems defined in Theorem 1 and in Corollary 1. The upper panels show possible trajectories for a set of balls under the conditions established by Theorem 1 — panel (a) and by Corollary 1 — panel (b). The bottom panels show the same, but considering a family of balls uniformly distributed along a circle.
Fig. 4. Some of the possibles trajectories until the first collision for the billiard system defined in the text. The upper panels show also some orientations of the velocity vectors at a given initial position for 0 = 0, in panel (a), and for 0 > 0, in panel (b). The bottom panels show the trajectories assuming a set of point-like balls uniformly distributed along a circle of radius r = 0.4R. In panel (c) 0 = 0, while in panel (d) 0 > 0
2.1.2. Further results
In this section we will combine the previous methods with results established in the theory of the elliptic functions. In this way we will build collision-based methods to compute elementary functions, integral functions and mathematical constants, exploring further the computational capabilities of the model. The equations we evoked can be found in [21] or in [22], where an enormous set of differential and integral relations between the functions F(0, k), E(0, k) and n(0, n, k) is listed. From a hardware designer perspective, our first six results define a logic family, while the theory of the elliptic functions is a mathematical model of computation. First we will reveal a collision-based method to compute the number n. After that, we will show how the elliptical integral of the second kind and some of its integrals can be computed by a similar mechanism. From these results we will show that the Legendre identity reveals a second collision-based method to compute n. Along the discussion we will show three collision-based methods to compute some elementary functions. Finally, we will show how the machine computes the arithmetic geometric mean from which the Gauss constant will be computed. Some results presented next are written in terms of the complementary modulus k' = \/l — k2.
• Computing 7r — method one
Corollary 6. Consider the results and the definitions established by Corollary 2 considering rd = k'R and in the particular case where l = -1. Then, as N — to and m — to,
* = {-k Ml-1)- <*■*>
The proof follows by combining the fact that
J F k') dk =
n
T
(2.27)
with the results given by Corollary 2 in the case where k = 1 and l = —1.
This result resembles the collision-based method recently developed by Galperin in [12], who showed that n can be computed from a billiard system of just two particles. In both methods it is admitted that we can count the total number of collision and, through this measurement, they show two examples of how the number n can be computed from Newton's laws of motion.
Similarly, by replacing k by k in Eq. (2.27), we can compute Catalan's constant G* once
F
k ) dk = 2G.
(2.28)
The next results show how the billiard computes the logarithmic function, given cos 1 x is precomputed.
Computing In y — method one
Corollary 7. Consider the results and the definitions established by Theorem 1 and by
2
VI+1
Corollary 1. Choose y G R and compute <fi = cos 1 ' ——
Then, as N —> oo, m
oo
and k 1-
Iny = F((p, k) = (f, k) - 1 - I) k)
(2.29)
The proof follows by combining the fact that. F(<f>, 1) = In j with the result estab-
lished in Corollary 1.
Now we will show how to combine the previous results in order to define a method to
compute the integral functions E , k) and E(<fi, k). In order to compute E(<fi, k) we will assume that sin 0 is precomputed.
• Computing E , k) and E(<p, k)
Corollary 8. Consider the results and the definitions established by Theorem 1 and Corollary 2. Let 0 ^ </> ^ f and 0 ^ k ^ 1. Then, as N —> oo and m —> oo,
t n \ { k2 od tn
(2.30)
Corollary 9. Consider the results and the definitions established by Theorem 1 and Corollaries 1, 2 and 3. Let 0 ^ </> ^ § and 0 ^ k ^ 1. Then, as N —> oo and m —> oo,
+
+(k')2 p
where f(<f>, k) = (l — \/l — k'2 sin2 <f>j cot q').
c 1f' k) ~ I1" i' k)
+ f (0,k), (2.31)
1
0
1
0
The proof of both corollaries follows by combining the fact that k
J F(0, k)kdk = E(0, k) - (k')2F(0, k) - f (0, k) (2.32)
0
with the results given in Theorem 1 and Corollaries 1, 2 and 3 with l = 0. • Computing sin</>
Corollary 10. Consider the results and the definitions established by Corollaries 2 and 3
7T 2 '
and set (f) such that 0 ^ </> ^ Then, as N —> oo and m —> oo
^ = IT^y (2-33)
where
»M*^) K-0 - H)
(2.34)
The proof of this corollary follows by combining the fact that E(0, 1) = sin 0 with the results given by Corollary 9 with k = 1. Alternatively, the proof follows by combining the fact that
1
[ F(<f>, k)kdk = 1 ~cos0 (2.35)
sin 0
0
with the result established in Corollary 3 with k = 1 and l = 0. • Computing 7r — method two
Corollary 11. Consider the results and the definitions established by Theorem 1 and Corollary 2 with 0 ^ k ^ 1. Then, as N — to and m — to,
_4 N2m2_
fc2/3d(§, fc)/?c(§, k') + {k')2!3d{l, k>)pc(^, k)'
The proof of this corollary follows by combining the results of Theorem 1 and Corollary 8 with the Legendre relation for complete elliptic integrals, which is given by
= l (,37)
This equation was also utilized by Salamin [24], while he established his famous formula to compute n.
Now we will establish methods to compute some integrals and derivatives of the functions F(0, k) and E(0, k).
k k
• Computing f E k) kdk and f E(<f>, k)kdk _o_o_
k
The following results shows how the machine computes the integral functions f E k) kdk
02
k
and f E(0, k)k dk. 0
Corollary 12. Consider the results and the definitions established by Theorem 1 and Corol-
lary 2. Let 0 ^ (j) ^ j and 0 ^ k ^ 1. Then, as N —> to and m —> to, k
0
(2.38)
Corollary 13. Consider the results and the definitions established by Theorem 1 and Corollaries 1, 2 and 3. Lei 0 ^ (f) ^ f and 0 ^ k ^ 1. Then, as N —> oo and m —> oo,
E(6, k)kdk =
/ic (f, fc) - (l - fc)]}
+(k')2 k
+
k2f (6, k) . + C (2.39)
3
The proof of both corollaries follows by combining the fact that
E(0, k)kdk = - [(1 + k2)E(0, k) - {k'fF{0, k) - /(</>, k)]
(2.40)
with the results established in Corollaries 1 and 9.
In [21, 22] a set of equations is listed which can be combined with the results presented here in order to design the collision-based method of computation for a set of integrals of the functions F(0, k), E(0, k) and n(0, n, k). We explore here only some of these possibilities.
Next, we will show a method to compute some derivatives. In particular, we will show how the derivatives of E k) and E(0, k) with respect to the modulus k can be computed from the collisional process.
• Computing derivatives of E (5, k) and E(0, k)
Corollary 14. Consider the results and the definitions established by Theorem 1 and Corollary 2. Set k such that 0 ^ k ^ 1. Then, as N — to and m — to,
f * \Jßd&k)
\2 Nm) \ 2
- ß
c'\k
(2.41)
Corollary 15. Consider the results and the definitions established by Theorem 1 and Corollaries 1, 2, 3 and 9. Set k such that 0 ^ k ^ 1. Then, as N — to and m — to,
v =
+ (2.42)
k
0
k
0
The proof of both corollaries follows by combining the results established in Corollaries 1 and 9 with the fact that
d 1
—£(</>, k) = -kDM, k) = ^(£(</>, k) - F(</>, *)), (2.43)
where D(0, k) is another frequently encountered elliptic integral defined by
D(<f>, k) = B1Jf ° de = k) - E(<f>, k)). (2.44)
J VI - k2 sin2 9 k 0
In a similar way we can also derive collision-based methods to compute the second-order derivative of E(0, k), and the first and second derivatives of F(0, k).
2.1.3. The arithmetic geometric mean
The arithmetic-geometric mean (AGM) between two positive real numbers, a0 and b0, is the common limit of the interdependent sequences {ai} and {bi}, defined by
ai+i = + h) (2-45)
and
h+i = (2-46)
with i £ N. This limit is usually denoted by M(a0, b0). Also, a related and useful auxiliary sequence is defined by
c2 = |a2 - b21. (2.47)
Gauss established a correlation between elliptic integrals and the AGM by showing that the symmetric form of the complete elliptic integral of the first kind, I(a0,b0), is invariant under
a specific change of variables. Given A;2 = 1 —= 1 — (A/)2, he showed that.
a0
F (I, k) = O0/(O0, W = (I) = (I) jjJ^j, (2.48)
and that
/ ro \
n , \ / _ , „in
2' J Ci r \ 2
\ i=0 /
Thus, Gauss revealed that the theory of elliptic functions can be founded on the function M(a0, b0) [27]. The next corollary shows how a set of circular billiards computes the arithmetic-geometric mean between any two positive numbers.
• Computing the arithmetic-geometric mean
Corollary 16. Consider the system, the definitions, and the initial conditions established in Theorem 1. Let M(a0, b0) be the arithmetic-geometric mean of two given positive real numbers a0
2 b2
and b0, with a0 ^ b0. Also define k = 1--Then, as N —> oo and m —> oo,
(Nm \
The proof of this corollary follows by combining the result established by Theorem 1 with the integral-form expression for M(a0, b0) derived by Gauss and with the symmetrical form of the elliptic integral of the first kind. If a,0 = \[2 and b0 = 1, we conclude that, the Gauss constant can be computed by
c_ 1 _Pc(l,k) M (y/2, 1) y/2Nm '
where k = N —> oo and m —> oo.
The fast convergence of the AGM iteration, together with its integral form derived by Gauss, makes it useful to write efficient computational algorithms. The remarkable papers cited in [24] and in [25] established fast computations of elementary functions and of the number n using results from elliptic integral theory. A collection of similar results is given in [23], where AGM-ba-sed methods to compute elementary functions, mathematical constants and relevant functions of mathematical physics are presented. Some of these results can be combined with the results presented here in order to expand the computational capabilities of the model we are proposing. As an example we will present another algorithm to compute logarithmic functions.
The following result is due to Salamin [24], who established an asymptotic formula to compute lny from the AGM recurrence.
• Computing ln y — method two
Set y > 4 and define
* = + (2.52)
y
Then, if k — 1,
fc)-lny, (2.53)
which may then be combined with Theorem 1. In this paper we will not try to assess more results that can be associated to our method. However, the AGM iteration lies at the heart of theta function theory (see [23, 27]) and a detailed analysis may lead to further applications of the model.
3. Discussion
In our method we admitted that the number of particles, N, goes to infinity. From a physical perspective, this hypothesis is satisfied by considering a continuous mass distribution, and we may think in terms of an expanding gas. If we imagine that the balls represent photons, we may hypothesize that the initial conditions can be realized from a cylindrical source of radiation that generates a cylindrical wavefront. Also, note that we admitted a boundary that acts as a fixed mirror, and it makes the model physically unrealistic. The momentum is not being conserved as the system behavior requires a circular boundary infinitely massive. This imposes an external architecture onto the computing substrate, which is not consistent with the architectureless concept in collision-based computing [30]. For a rich discussion about the possibility to build a collision-based computer based on optical phenomena the reader is referred to the work [5] and references therein, where the concept of a soliton machine is presented.
Suppose that a dynamical variable is described in terms of a set of elliptic integrals that are computable from the model presented above. A textbook example is the exact solution to
the plane pendulum, where the period is written in terms of the function F , k), with k being defined from the initial conditions of the mass [28]. Given this, our results reveal an alternative description for that variable, a description which is based on the collisional process occurring in a set of idealized dynamical systems.
From the framework of Automata Theory, this work defines a collection of block cellular automata that, through the concept of the Margolus neighborhood [7], can be used to compute some key functions of the theory of elliptic functions. The states of the cells in each CA model represent idealized particles, and its rules are designed in order to encode the particle movement inside a circular boundary, with which the balls collide elastically. The theorems and their corollaries specify how the initial conditions of the cells must be defined in order to achieve a set of desired computational tasks, like computing elementary and integral functions and mathematical constants. Results from elliptical integral theory were evoked in constructing these results, and in this context it served as a mathematical model of computation. The continuum character of these functions was reconciled with the discrete nature of CA considering three limiting cases. The number of balls of the system goes to infinity, the dimension of the cells (the balls radius) goes to zero, while the execution time goes to infinity. From this, the model we outlined here provides an alternative CA description for any dynamical variables which can be described in terms of any combination of the functions that can be computed from the previous results. This contrasts with most of the recent studies found in the literature which in general model specific dynamical systems with CA models [13].
While the physical implementation of the dynamical systems considered here is a technical challenge, the implementation of these methods from dynamical billiards simulations is trivial, and this in fact guided the author during this study.
The present work could also be presented within the framework of measure theory, on which the dynamical billiard description is based [26]. In our analysis, we replaced the time average by a space average, which could be analyzed within the framework of ergodic theory. However, we preferred a presentation more directed towards practical applications.
We hope that the present work may serve as a small contribution to the development of a formal method for creating CA models that perform a desired computation, a subject that is extremely difficult and that has no ideal mathematical treatment yet [13, 33].
4. Conclusions and some final remarks
We have found it surprising that our analysis, which is built upon the simplest billiard table, revealed such a large number of supposedly new results. However, we have not been able to find any studies that linked a circular billiard table with the theory of elliptic functions. The circular billiard is considered to be a trivial subject of study from a more specialized mathematical perspective [31]. Despite that, we have showed that these systems have computational capabilities which, to the best of our knowledge, remain unexplored.
This fact is in line with the view of the authors of a similar study [4], which claimed that a neutral intellectual space in the interface between physical theories and computability theory is needed. Also, we inevitably remember some words written by Konrad Zuse in the preface of his essay Calculating Space [1]. There, he highlights that the close interplay between mathematicians, physicists and data processing specialists may have a favorable effect on the development of models in these fields. At that time, he had already identified that these areas were splitting in a crescent number of subfields, each one of them demanding years of specialized
study. Perhaps this effect explains why the results presented here have remained unexplored for such a long time.
A natural question that occurs to us is whether there will be a redefining in our understanding of the physical world if one day we are able to formulate CA that computes not only elliptic integrals, but all special functions of mathematical physics and even that serves as an alternative to differential equations. Some of these ideas have been extensively discussed for years [7, 14, 17] and we will finish our discussion reproducing the words of one of the main actors in this field.
We hypothesize that there will be found a single cellular automaton rule that models all of microscopic physics, and models it exactly. We call this field DM, for digital mechanics. (Eduard Fredkin, in [33].)
The author would like to thank the anonymous referees who provided valuable suggestions.
References
[1] Zuse, K., Rechnender Raum, Braunschweig: Vieweg, 1969.
[2] Sauder, J., Hilgemann, E., Johnson, M., Parness, A., Bienstock, B., Hall, J., Kawata, J., and Stack, K., Automaton Rover for Extreme Environments, Technical Report, California: California Institute of Technology, 2017.
[3] Fredkin, E. and Toffoli, T., Conservative Logic, Int. J. Theor. Phys., 1982, vol. 21, nos. 3-4, pp. 219253.
[4] Beggs, E. J. and Tucker, J. V., Embedding Infinitely Parallel Computation in Newtonian Kinematics, Appl. Math. Comput, 2006, vol. 178, no. 1, pp. 25-43.
[5] Jakubowski, M. J., Computing with Solitons in Bulk Media, PhD Thesis, Princeton, N.J., Princeton Univ., 1999, 158 pp.
[6] Rozenberg, G., Bäck, Th., and Kok, J.N., Handbook of Natural Computing, Berlin: Springer, 2012.
[7] Toffoli, T. and Margolus, N., Cellular Automata Machines: A New Environment for Modeling, Cambridge: MIT Press, 1987.
[8] Richard, K. S. and Steiglitz, K., Programmable Parallel Arithmetic in Cellular Automata Using a Particle Model, Complex Syst., 1994, vol. 8, no. 5, pp. 311-323.
[9] Beggs, E. J. and Tucker, J. V., Experimental Computation of Real Numbers by Newtonian Machines, Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci, 2007, vol. 463, no. 2082, pp. 1541-1561.
[10] Beggs, E.J. and Tucker, J. V., Can Newtonian Systems, Bounded in Space, Time, Mass and Energy Compute All Functions?, Theoret. Comput. Sci., 2007, vol. 371, nos. 1-2, pp. 4-19.
[11] Collision-Based Computing, A.Adamatzky (Ed.), London: Springer, 2002.
[12] Galperin, G. A., Playing Pool with n (the Number n from a Billiard Point of View), Regul. Chaotic Dyn, 2003, vol. 8, no. 4, pp. 375-394.
[13] Ganguly, N., Sikdar, B.K., Deutsch, A., Canright, G., and Chaudhuri, P.A., Survey on Cellular Automata, Technical Report No. 9, Centre for High Performance Computing, Dresden University of Technology, 2003, 30 pp.
[14] Smith, M. A., Cellular Automata Methods in Mathematical Physics, PhD Thesis, Cambridge, Mass., MIT, 1994, 243 pp.
[15] Wolfram, S., Universality and Complexity in Cellular Automata, Phys. D, 1984, vol. 10, nos. 1-2, pp. 1-35.
[16] Dascalu, M., Cellular Automata Hardware Implementations: An Overview, Rom. J. Inf. Sci. Tech-nol., 2016, vol. 19, no. 4, pp. 360-368.
[17] Toffoli, T., Cellular Automata As an Alternative to (Rather Than an Approximation of) Differential Equations in Modeling Physics, Phys. D, 1984, vol. 10, nos. 1-2, pp. 117-127.
[18] Margolus, N., Physics-Like Models of Computation, Phys. D, 1984, vol. 10, nos. 1-2, pp. 81-95.
[19] Fredkin, E., An Introduction to Digital Philosophy, Internat. J. Theoret. Phys, 2003, vol. 42, no. 2, pp.189-247.
[20] Toffoli, T., CAM: A High-Performance Cellular-Automaton, Phys. D, 1984, vol. 10, nos. 1-2, pp. 195-204.
[21] Byrd, P. F. and Friedman, M. D., Handbook of Elliptic Integrals for Engineers and Scientists, 2nd ed., rev., Grundlehren Math. Wiss., vol. 67, New York: Springer, 1971.
[22] Hancock, H., Elliptic Integrals, New York: Dover, 1958.
[23] Borwein, J. M. and Borwein, P. B., Pi and the AGM: A Study in Analytic Number Theory and Computational Complexity, New York: Wiley, 1987.
[24] Salamin, E., Computation of n Using Arithmetic-Geometric Mean, Math. Comp., 1976, vol. 30, no. 135, pp. 565-570.
[25] Brent, J. M., Fast Multiple-Precision Evaluation of Elementary Functions, J. Assoc. Comput. Mach., 1976, vol. 23, no. 2, pp. 242-251.
[26] Chernov, N. and Markarian, R., Chaotic Billiards, Math. Surv. Monogr., vol. 127, Providence, R.I.: AMS, 2006.
[27] Bellman, R., A Brief Introduction to Theta Functions, New York: Holt, Rinehart & Winston, 1961.
[28] Thornton, S.T. and Marion, J.B., Classical Dynamics of Particles and Systems, 5th ed., Pacific Grove, Calif.: Brooks/Cole, 2004.
[29] Aretxabaleta, X. M., Gonchenko, M., Harshman, N.L., Jackson, S.G., Olshanii, M., and As-trakharchik, G.E., The Dynamics of Digits: Calculating Pi with Galperin's Billiards, Mathematics, 2020, vol. 8, no. 4, Art. 509, 32 pp.
[30] Zhang, L., Computing Naturally in the Billiard Ball Model, in UC'09: Proc. of the 8th Internat. Conf. on Unconventional Computation (Ponta Delgada, Portugal, Sept 2009), Berlin: Springer,
2009, pp. 277-286.
[31] Berger, B., Geometry Revealed: A Jacob's Ladder to Modern Higher Geometry, Berlin: Springer,
2010.
[32] Chang, S.-J. and Friedberg, R., Elliptical Billiards and Poncelet's Theorem, J. Math. Phys., 1988, vol. 29, no. 7, pp. 1537-1550.
[33] Fredkin, E., Digital Mechanics: An Informational Process Based on Reversible Universal Cellular Automata, Phys. D, 1990, vol. 45, nos. 1-3, pp. 254-270.
