Deutsche Version

The Spiral of Theodorus

The Spiral of Theodorus starts with an isosceles right triangle with both legs of length 1. More right triangles are added, one leg the hypotenuse of the previous triangle, the other, outside leg, always of length 1.

The Pythagorean Theorem tells us that the hypotenuses of these triangles have lengths etc. .

This Geogebra Applet dealing with the Spiral of Theodorus and its interpolating curve uses methods developed here.

The Spiral of Theodorus is more than a nice application of the Pythagorean Theorem. If we continue to plot the spiral there are more questions which need more sophisticated mathematics to be answered.

  1. The distances between the windings look very regular and seem to approach π. So do we have approximately an Archimedean Spiral, which can be described in polar coordinates by the equation r(φ) = a · φ + b ?
  2. If yes, which are the values of the parameters a and b?
  3. The first revolution is completed with 17 triangles, the second with 54 etc.. Is there any formula to describe this sequence 17, 54, 110, 186, 281.....?

If you add the displayed nth tringle, the increase of the angle is = .
The corresponding increase of the radius is.

For large n is about and about half the size, namely , which can be shown easily with the Taylor series expansion.

For each full revolution, the increase of the angle is 2π, thus the increase of the radius is about half this size, namely π.

This argumentation only shows point 1. To make this more precise, we we use the Euler–Maclaurin formula.

The Euler–Maclaurin formula provides a way to evaluate finite sums using integrals.



Here Br are the Bernoulli numbers and R(k,a,b) a remainder term, which depends on a, b and the parameter k >=1, which can be choosen freely.

We have

R(k,a,b) =                                                                                 (2)

Bk(x) is the k-th Bernoulli polynomial.
We have for example B1(x) = x - 1/2, B2(x) = x2 - x + 1/6 und B3(x) = x3 - 3/2 x2 + 1/2 x.

We use (1) to determine the total angle φ(n-1) of the first n - 1 triangles of the square root spiral. We have

φ(n -1) = , with f(x) = , n>=1.

In (1) and (2) we therefore select a=1 and b=n.
The hypotenuse associated with the total angle φ(n-1) has the length r = .

Using an estimation of D.H. Lehmer (American Mathematical Monthly, volume 47, pages 533–538 (1940)) which shows that |Bk(x)| in [0,1] essentially is limited by , we can conclude that with f(x) = the limit of R(k,1,n) in (2) exists for n → ∞ and moreover we have
|R(k,1,∞)| < || .

This means that there is a constant C (depending on k) such that φ(n -1) can be approximated with an arbitrary degree of precision by


for sufficiently large n.

We now make a Taylor expansion of (3) at n = ∞ to the order g and replace by r. It can be shown that the terms up to order g are independent of k as long as kg.

With Mathematica the computation can be performed with the following program:

g=6;(*Example for order 6*)

Expanding to the order g we then get, with a universal constant K

φ(r)= +.... +O(1/r2g+1)       (4)

Here r = is the length of the hypotenuse of the n th triangle and φ(r) the corresponding angle sum of the first n triangles..
With an expansion to order g the highest exponent of r is 2g - 1. In (4) for example we see the terms for g = 6.

We show below that K = - 2.157782996659446220929...
E. Hlawka (Gleichverteilung und Quadratwurzelschnecke, Monatsh. Math., 89 (1980) 19-44) describes this constant K as "Schneckenkonstante".

In the zeroth-order approximation we get φ(r) ≈ K + 2r or r(φ) ≈ φ/2 - K/2.

This allows us to answer the second and third question above. For r → ∞ we have an Archimedean Spiral with a = 1/2 and b = - K/2.

In the first-order approximation we get φ(r) ≈

Even for small r this approximation is very good.
The blue curve gives ( r | K+2r+1/(6r) ) in polar coordinates.

To compute K, we observe that φ(r) is the sum of the angles of the first r2 -1 triangles,
φ(r) =.

This leads to


If we take into account 1/r199 (order g=100) for the highest power of r, r = 1000 suffices to compute K with an accuracy of more than 500 decimal places:


Alternatively, K can be calculated elegantly with an infinite series using the Riemann zetafunction
(D. Brink, The spiral of Theodorus and sums of zeta-values at the half-integers, The American Mathematical Monthly, Vol. 119, No. 9 (November 2012), pp. 779-786).

φ(r) = +..........

is also a good starting point for a differentiable continuation for all real r>0 (see D. Gronau: The spiral of Theodorus. Amer. Math. Monthly 111, 2004, 230-237).

While the Gronaus representation for φ(r) converges very slowly, the representation (4) expanded to a certain order g gives an arbitrarily exact approximation for all rRg, where Rg depends on the order g.

Particulary interesting is the interval 0< r<1, but the representation (4) is inappropriate here.
But we have

and with a sufficiently large n we can use approximation (4) for 0< r< 1.

In the Geogebra Applet we use terms up to the order g=7 and n=10 . This guarantees an error < 10-10 for φ(r) for all r>0.

The following functions are very good approximations and lower and upper bounds for φ(r) with an error of order O(1/r3):

φL(r) =

φU(r) =

To see this, we compare the Taylor series of φL(r) and φU(r), expanded at r = ∞, with equation (4), which includes the term 1/(120 r3):

φL(r) =

φU(r) =

Surprisingly we even have

φL(r) < φ(r) < φU(r)        for all r≥0,

demonstrated by the following graphs:


                φ(r) - φL(r)                                                         φU(r) - φ(r)

φL and φU have simple inverse functions.

With a = (φ - K)/2 we have:

rL(φ) := φU-1(φ) = , and

rH(φ) := φL-1(φ) = and hence

rL(φ) < r(φ) < rH(φ),      φ > Sqrt[6]/3 + K

Sequence A072895 from the On-Line Encyclopedia of Integer Sequences is defined as the least number of triangles to complete k revolutions. Its elements are

17, 54, 110, 186, 281, 396, 532, 686, 861, 1055, 1269, 1503, 1757, 2030, 2323, 2636, 2968, 3320, 3692, 4084, 4495, 4927, 5377, 5848, 6338.....

We have

rL(2kπ) < r(2kπ) < rH(2kπ) and because n = r2 -1

Floor[ rL(2kπ)2 ] ≤ A072895(k) ≤ Floor[ rH(2kπ)2 ]

For all k≤1 000 000 000 Floor[ rL(2kπ)2 ] = Floor[ rH(2kπ)2 ] holds, so we have

A072895(k) = Floor[ rL(2kπ)2 ] = Floor[ (kπ - K/2)2 - 1/6 ] , k≤1 000 000 000

and also

A072895(k) = Floor[ rH(2kπ)2 ] = Floor[ (kπ - K/2)2 - 1/6 + 1/(12kπ - 6K)2 ], k≤1 000 000 000


The following Mathematica Program computes A072895(k) :

a[k_]:=Module[{ a= -K/2 + kπ, b}, b = a2 - 1/6; If [Floor[b]==Floor[b+1/(12 a)2], Floor[b], Undefined]]

Under the assumption of an equidistribution mod 1, the probability that Floor[b] != Floor[b+1/(12 a)2] is 1/(12kπ - 6K)2 for each k.

The probability that a[k] is defined for all k > 1 000 000 000 is

Product [1 - 1/(12kπ - 6K)2,{k,1000000001,∞}] ≈0.99999999999929638

With high probability, therefore, a[k] = A072895(k) for all k.


< Home >

© 2018  Herbert Kociemba