As Xander Henderson and others have mentioned, this operation converges on the solution of the equation $\cos x = x$, and have explained why it converges, so I won't cover that fact in this answer. I'll merely mention that if you look at the graphs of $y = \cos x$ and $y = x$ it's pretty obvious that there's one and only one real solution to $\cos x = x$. However, I noticed that none of the existing answers actually give a value for this $x$, sometimes known as Dottie's Number, or mention a faster way to evaluate it, and I intend to remedy that situation.
It has also been mentioned that there's no closed-form solution to this equation using elementary functions. It converges reasonably quickly, as you've discovered, (and as I discovered several decades ago, when scientific calculators first appeared :) ), but we can easily use a little bit of calculus to find an algorithm that converges much more rapidly.
This is a classic application of the well-known Newton-Raphson method.
We want to solve $y = x - \cos x$ for $y = 0$
$$y = x - \cos x$$
Differentiating,
$$dy/dx = 1 + \sin x$$
Now $\frac{dy}{dx} \approx \Delta y / \Delta x$
So we can let
$$\Delta x = \frac{\Delta y}{1 + \sin x}$$
Now
$$\Delta y = -y = \cos x - x$$
Thus
$$\Delta x = \frac{\cos x - x}{1 + \sin x}$$
We could stop there, but if we continue we get a nice cancellation.
Let $x'$ be the new approximation for $x$. That is
$$\begin{align}
x' & = x + \Delta x\\
& = x + \frac{\cos x - x}{1 + \sin x}\\
& = \frac{x + x \sin x + \cos x - x}{1 + \sin x}\\
x' & = \frac{x \sin x + \cos x}{1 + \sin x}\\
\end{align}$$
Here's a Python program which uses that formula. We only need to do a few loops to reach the limit of precision of 64 bit floating-point numbers.
from math import sin, cos
x = 1
for i in range(4):
print(i, x)
s = sin(x)
x = (cos(x) + s * x) / (s + 1)
print(x, cos(x))
###output
0 1
1 0.7503638678402439
2 0.7391128909113617
3 0.7390851333852839
0.7390851332151607 0.7390851332151607
Using the arbitrary precision mpmath library, we can easily get a lot more digits with only a few more loops.
from mpmath import mp
def print_digits(s, colsize=5, rowsize=50, blocksize=250):
''' Print high precision decimal string in columns, rows, & blocks '''
# Split off integer part & decimal point
*parts, s = s.partition('.')
print(''.join(parts))
# Format digits after the decimal place
cols = []
for i in range(0, len(s), colsize):
j = i + colsize
# Columns are separated by a space,
# rows are separated by a single newline,
# blocks are separated by two newlines.
sep = (' ' if j % rowsize
else '\n' if j % blocksize else '\n\n')
cols.extend([s[i:j], sep])
print(''.join(cols))
# Set the precision to 2000 decimal places
mp.dps = 2000
sin, cos = mp.sin, mp.cos
x = 1
for i in range(11):
s = sin(x)
x = (cos(x) + s * x) / (s + 1)
print_digits(str(x))
###output
0.
73908 51332 15160 64165 53120 87673 87340 40134 11758 90075
74649 65680 63577 32846 54883 54759 45993 76106 93176 65318
49801 24664 39871 63027 71490 36913 08420 31578 04405 74620
77868 85249 03891 53928 94388 45095 23480 13356 31276 77223
15809 56353 77657 24512 04373 41993 64335 12538 40978 00343
40646 70047 94021 43478 08027 18018 83771 13613 82042 06631
63350 37277 99169 67312 23230 06138 86582 03621 77081 09978
97062 68424 05880 94898 68326 18606 00485 89895 85487 25736
76401 50752 27608 18039 14595 18101 62815 91200 96461 64606
75440 51326 41517 10644 66281 10936 08258 48783 71383 95555
61751 41494 71593 90062 77527 56325 86349 38869 73014 08366
51525 11520 42678 85153 02529 41718 03651 76420 17708 60718
99276 01609 87432 71545 52267 56579 82462 97611 77553 96166
99549 31115 85665 34834 95383 85231 59636 02527 49955 87252
50666 64013 13187 40139 25388 88055 20618 69859 21392 52528
54154 11079 10029 98282 92986 40521 69046 55473 66968 71438
73564 60065 21225 46891 49975 92096 99758 50136 42495 08565
04732 49725 84248 37155 48364 83437 27583 74675 25453 35800
66420 04788 39718 85848 90145 31155 06041 78123 37047 77395
34717 10345 11958 54600 72656 14647 21419 78753 73880 23680
29553 44127 94853 01620 77437 43315 90133 91933 23148 76628
28552 17782 70052 31111 78246 86229 57127 86199 58490 58929
78171 80601 56715 85092 53714 04181 46882 85824 54046 44526
55883 15798 59786 67282 99052 07226 86870 94531 30864 95350
44481 38762 32367 76569 23613 25971 52294 15582 29334 15223
69636 98322 65805 17766 85366 37759 37066 43679 29565 98287
11924 91109 47930 11267 60115 22614 29243 71121 71487 02935
43102 93038 78065 42303 10930 07600 02409 80335 56727 30891
51766 68247 56247 72917 20259 45634 73838 58189 99548 42071
81825 61281 93120 90775 73740 22398 85853 96422 00631 32270
25328 45511 60110 76298 67410 60230 99696 24600 46288 53977
84470 79428 66907 48442 61986 19229 57526 39675 15917 84325
95550 97753 94210 80967 31399 81533 82554 94805 27727 29806
59555 18624 34873 99989 44589 81750 86210 57831 26013 92783
08279 72906 91694 42203 86806 52977 23102 01474 61968 67501
16979 00223 36525 95179 73670 54931 31976 05150 85764 33674
63449 04471 87196 10601 81776 20236 18885 83971 51415 31833
42332 51271 72162 70117 51797 43026 02487 76822 00067 45587
38185 02965 94214 70461 70423 42065 97465 13966 68395 61557
29100 71089 57281 30968 30388 12301 43385 35545 47878 98962
We can simplify this formula. The simplified version doesn't converge as rapidly, but it's easier to compute, and it still converges faster than iterating $x = \cos x$
The sine of Dottie's number is $\approx 0.673612$. We can round that to $\frac{2}{3}$ and plug that into our formula. It will still converge even if we start with $x = 1$, but it's better to start with a closer approximation, say $x = \frac{3}{4}$. Note that $\left(\frac{2}{3}\right)^2 + \left(\frac{3}{4}\right)^2 = \frac{64 + 81}{144} = \frac{145}{144} \approx 1$.
$$x' = \frac{x \sin x + \cos x}{1 + \sin x}$$
Substituting $\sin x = \frac{2}{3}$
$$\begin{align}
x' & = \frac{\frac{2}{3} x + \cos x}{1 + \frac{2}{3}}\\
& = \frac{\frac{2}{3} x + \cos x}{\frac{5}{3}}\\
& = \frac{2 x + 3 \cos x}{5}\\
x' & = 0.4 x + 0.6 \cos x
\end{align}$$
In other words, we perform a weighted mean of $x$ and $\cos x$. This still converges to the correct value. At convergence,
$$\begin{align}
x & = \frac{2 x + 3 \cos x}{5}\\
5x & = 2x + 3 \cos x\\
x & = \cos x
\end{align}$$
And here's a short demo in Python.
from math import cos
x = 0.75
for i in range(8):
y = cos(x)
print(i, x, y)
x = 0.4 * x + 0.6 * y
###output
0 0.75 0.7316888688738209
1 0.7390133213242926 0.7391335046629345
2 0.7390854313274777 0.7390849324030849
3 0.739085131972842 0.7390851340520015
4 0.7390851332203376 0.7390851332116734
5 0.7390851332151391 0.7390851332151751
6 0.7390851332151607 0.7390851332151607
7 0.7390851332151607 0.7390851332151607