# C, Julia, Python, Maxima, Mathematica, ChatGPT and numerical errors

I am completely stunned how numerical errors can diverge for so innocent programs.

In Python 3.11.7 the program

import math
u = 5.0
for i in range(73):
u = 6*math.cos(u)
print(u)


displays $$-3.8591995867613633$$.

In Julia 1.9.2, the program

u = 5.0
for i=0:72
u = 6.0*cos(u)
end
print(u)


displays $$-3.770177449898713$$.

For Maxima 5.45.1, the variable un in the code,

u0:5.0;
un:u0;
for k:0 while k<73 do (
un:6.0*cos(un));
print(un)


displays the same value as the Python code.

In the C language, the program

#include <stdio.h>
#include <math.h>

int main() {
float u = 5.0;
for (int i=0;i<73;i++){
u = 6.0*cos(u);
}
printf("%f",u);
}


displays $$1.989200$$.

In Mathematica, the program

u=5.0
For[i=0,i<73,i++,u = 6.0*Cos[u]]
Print[u]


prints the same as the Python's or Maxima's code.

I am not saying Julia or C is false, because I can't say that the value given by Maxima, Mathematica or Python is correct. It seems different version of Python 3 could give different results (on the same architecture). For instance, see the output of ChatGPT here (that's unbelievable):

It should come from the different ways Python and Julia are representing float number in memory, but still that's very strange.

• In C, can you try using double?
– NNN
Commented Feb 15 at 13:41
• Generally, ChatGPT produces a conversation piece, a good-looking answer, there is no guarantee for internal consistency, here for instance that the code displayed is actually executed. Commented Feb 15 at 19:40
• I recommend you read up on chaos theory. This is quite similar to a logistic map and is therefore expected behavior. The problem is not computer math libraries but an intrinsically chaotic system with sensitive dependence on initial conditions. Commented Feb 16 at 1:05
• Out of curiosity, I've computed this using PARI/GP with large precision (first 1,000 digits, then 10,000 digits). AFAIK they use MPFR so I would trust their results. Thus, the first few digits of the (supposedly "correct") result should be: 2.7767584084827734649172502497247816. Commented Feb 16 at 2:09
• Nothing ChatGPT generates should be taken seriously at any time. Commented Feb 16 at 10:07

This function, as any similar hill, tent or hat functions, implements the expansion-and-folding scheme. This means that in the most benign interpretation you lose in each iteration one bit from the initial mantissa from the front. This is then filled from the back by the quasi-random noise of the truncation errors of the floating-point operations, to the size of 2.5 bit per iteration.

Thus one can expect that latest after 53 iterations the result is as good as random, the more so after 73 iterations. So to get a display precision of 30 bit, 9 decimal digits, after 73 iterations, one would need a working precision of $$73+30+2.5*73=286$$ bit in some multi-precision package.

The non-linear nature of the function might impose some attractor like a periodic orbit that reduces the randomness of the result. This is for the actual map of floating point numbers, so it may be different in period and distribution of the cycle for different compilers, math libraries and optimization options.

• To illustrate the point, with Intel Fortran and compiling with /fp:strict, using double precision (REAL*8) yields -3.8591995867613633, using quadruple precision (REAL*16) yields 1.9436192152480203, and using a multi-precision library (Brent's MP configured for 350-bit precision) yields 2.7767584084827734. Commented Feb 15 at 22:38
• Is the "expansion-and-folding scheme" a term? I take it that you mean the scaling of $u$ by a factor larger than one (I assume that's "expansion") and the subsequent compression by the $\cos$ to $[-1,1]$ ("folding")? Commented Feb 16 at 15:01
• @njuffa Assuming their appearance in the source code, REAL*8 and REAL*16 to set a precision sound so old compared to use a selected_real_kind, or to pick one from the intrinsic fortran_iso_env module today. Old in a tune similar to depreciated (example) for new code to be written and a detail to consider a modernization. Commented Feb 16 at 15:09
• @njuffa The gfortran compiler of course accepts a REAL*8 declaration. Reading it is a bit like going back to the 1990s, when textbooks presented both FORTRAN77 and Fortran90/95 side by side; a bit different to today (e.g., Curcic's Modern Fortran with a free excerpt). Good to see a comparison with Intel's ifort (and now ifx). Commented Feb 16 at 19:45

I thought it would be interesting to see how the number of bits of precision affects the error. I wrote the following using arb, a library for interval arithmetic (and in particular a library that tracks precision loss due to machine arithmetic). Instead of working with numbers, arb works with intervals, and when a machine instruction causes a loss of precision the intervals get larger.

// arblib in cpp
#include <iostream>
#include "arb.h"

int main() {
for (slong prec : { 32, 64, 128, 256, 512 } ){
arb_t x;
arb_init(x);
arb_set_si(x, 5);      // x == 5

arb_printn(x, 50, 0);  // print x to 50 decimal digits if available
std::cout << "\n";

for (int i = 0; i < 73; i++) {
arb_cos(x, x, prec);        // set x = cos(x), computed with prec precision
arb_mul_ui(x, x, 6, prec);  // x = 6 * x
std::cout << i << "  ";
arb_printn(x, 50, 0);
std::cout << "\n";
}
arb_clear(x);
}
}


With 32 bits of precision, the first lines of output look like

5.0000000000000000000000000000000000000000000000000
0  [1.70197311 +/- 3.48e-9]
1  [-0.78480545 +/- 8.22e-9]
2  [4.2451546 +/- 6.14e-8]
3  [-2.702513 +/- 4.55e-7]
4  [-5.43086 +/- 2.41e-6]


By the 10th iteration, it's completely hopeless.

9  [5.7 +/- 0.0451]
10  [5e+0 +/- 0.191]
11  [2e+0 +/- 0.809]
12  [+/- 7.01]


With 64 bits of precision, it works a bit better.

5.0000000000000000000000000000000000000000000000000
0  [1.701973112779357587 +/- 5.61e-19]
1  [-0.78480545237593126 +/- 2.53e-18]
2  [4.2451546016343455 +/- 2.28e-17]
3  [-2.7025128979996596 +/- 9.55e-17]
4  [-5.430859490891791 +/- 5.27e-16]
...
22  [5.0 +/- 0.0819]
23  [2e+0 +/- 0.325]
24  [+/- 3.62]


Even though this is a different library, the loss of precision due to arithmetic precision will always be present.

With 128 bits of precision, we have

5.0000000000000000000000000000000000000000000000000
0  [1.70197311277935758679983502908134385000 +/- 9.16e-39]
1  [-0.7848054523759312613033562894196438204 +/- 9.13e-38]
2  [4.245154601634345486781290325089041750 +/- 7.45e-37]
3  [-2.70251289799965963668557149619904016 +/- 4.82e-36]
...
46  [5.66 +/- 8.81e-3]
47  [4.9 +/- 0.0749]
48  [+/- 1.14]


Much further! With 256 bits of precision, we have

5.0000000000000000000000000000000000000000000000000
0  [1.7019731127793575867998350290813438500065355535133 +/- 4.34e-51]
1  [-0.78480545237593126130335628941964382041108157876186 +/- 3.49e-51]
2  [4.2451546016343454867812903250890417494462853687055 +/- 1.75e-50]
3  [-2.7025128979996596366855714961990401635527246365811 +/- 1.16e-50]
...
70  [5.758713845316529420111 +/- 2.17e-22]
71  [5.193532385598491745861 +/- 8.95e-22]
72  [2.77675840848277346492 +/- 6.34e-21]


That is, the first several digits on an exact calculator should read 2.776757... But even using 256 bits only gives 21 decimal digits after all the iterations.

With 512 bits of precision (which we now expect to be far more than enough), we confirm

5.0000000000000000000000000000000000000000000000000
0  [1.7019731127793575867998350290813438500065355535133 +/- 4.34e-51]
1  [-0.78480545237593126130335628941964382041108157876186 +/- 3.49e-51]
2  [4.2451546016343454867812903250890417494462853687055 +/- 1.75e-50]
...
70  [5.7587138453165294201111165200927171983191371489036 +/- 1.65e-50]
71  [5.1935323855984917458612965663072073056805703605339 +/- 3.50e-51]
72  [2.7767584084827734649172502497247816301002663316858 +/- 1.81e-50]


The e-50 error terms here are because I was printing the first 50 digits, and are not an indicator of the actual precision bound in the end.

Other answers have explained the results from Python, Julia, Maxima, C, and Mathematica. I'll explain the result you got from ChatGPT.

Generally, the way that ChatGPT works is that you ask a question, and then it gives you something resembling an answer to a question resembling the question you asked. Sometimes, you get lucky and the thing you get back really is the answer to the question you asked. Just as often, you get back something which looks like an answer, but is actually completely made up and wrong.

It looks like you're using ChatGPT 3.5. Well, ChatGPT 3.5 is completely unable to run Python or Julia code, so every single statement that it made about running the code is false. It's also unable to do the calculation in the program, so if it tries to print the answer, it will always print a guess instead of actually calculating the answer.

For the first response, ChatGPT almost certainly generated a supposed answer at random and then falsely wrote that its randomly-generated answer was the result of running the Python code. If you ask ChatGPT the same question several times in different chat sessions, you'll probably get a different random result every time.

For the second response, ChatGPT almost certainly copied the answer from its previous response with a random tiny change, and falsely stated that it got the answer by running that Julia code. It was merely guessing that Python and Julia would both return almost exactly the same answer.

For the third response, you must have told ChatGPT what answer you were expecting. ChatGPT almost certainly copied the number from your question into its own answer and falsely stated that it re-ran the Python code and got that answer.

As you can see, ChatGPT is not reliable at all. You should assume that all of its answers are always randomly generated and false unless you have a good reason to believe otherwise.

Perhaps the reason you are stunned is that you don't understand that what you are asking to calculate is impossible to ever calculate exactly. So if you want the answer to some precision, you need to know if the calculator you are using has that precision.

The cos() function basically takes a number of radians as input, and walks around a 2D unit circle that many radians, and tells you where on the 1D x-axis you ended up after the walk.

The problem with this is... the unit circle has 2*pi radians, and pi is an irrational transcendental number. Thus the cos result is irrational unless you have moved in units of pi, which you don't. So in a practical sense: the position you end up as has to be estimated.

Then you take that estimated distance, and use it to walk around the circle again, and again your position has to be estimated, so now your position on the circle is wrong by two estimates. And it accumulates with each iteration.

Then you must remember that you are using 32/64bit computer floating point. So you can guestimate that you are going to have 1 bit of error with each iteration, so that by 73, you surely have nothing but error, and your position on the unit circle is going to be anywhere from -1 to 1, seemingly randomly, depending only on the method of the estimation (giving your result as anywhere from -6 to 6).

• Great answer. I'd put even more emphasis on the meaning of "wrong by two estimates. And it accumulates with each iteration". As you say in the next paragraph, there is about one bit of error with each iteration, which means the error accumulates in an exponential sense, not in an additive sense. In other words, the error on an expression like 6cos(6cos(...6cos(5.0))) should be expected to be much much worse than the error on an expression like 6cos(1) + 6cos(2) + ... + 6cos(73).
– Stef
Commented Feb 17 at 12:26

Unfortunately this comes down to the math library. IEEE754 does not mandate exactly how the cosine should be computed (technically one doesn't even need to include a function to compute the cosine to be standard-compliant). In particular, implementations can differ in their last digits and do not need to be exactly rounded.

The function you are computing is very sensitive to rounding errors in the first steps (the word "chaotic map" comes to mind but I'm not an expert in the area). So that's why in the end the final results are very different.

That said, I am a little surprised that Python and Julia do not use the exact same math library for this computation.

EDIT: Matlab gives the same result as Julia. Could this be due to both of them using the vector math library from MKL? I'd be curious to see if the result is the same on a non-Intel processor.

• Julia uses a pure Julia implementation of cos that is based on fdlibm Commented Feb 15 at 15:51
• Thanks, very interesting! This would maybe deserve an answer of its own. For the curious, I have located the source here. Commented Feb 15 at 16:18
• we use our own implimentation of all the libm type things (because that way you can make sure that different OSes get the same answer and windows libm is often pretty slow) Commented Feb 15 at 19:04
• Thanks! When you say 'we', you mean that you are a Julia developer? Or in a different context? Commented Feb 16 at 10:20
• I am a julia developer. Commented Feb 16 at 16:15

We could use exact arithmetic to compute the correct value. A range of exact numbers $$e$$ correspond to the same $$a=5.0$$ representation in float32 IEEE754, plot over this range:

You can see that after 12 iterations, function value becomes sensitive to the choice of $$e$$ in the valid range which is approximately $$(4.9999990, 5.000001)$$. This choice is left up to your implementation , so you'll get different answers for different implementations

ClearAll["Global*"];
wp = 100;
iters = 20;
mantissa = 24;
targetValue = 5;
SF = StringForm;
f[x_?NumericQ, iters_] := Nest[6   Cos[#] &, x, iters];
interval =
Interval@SetPrecision[targetValue, Floor[mantissa*Log[10, 2]]];
minmax = SetPrecision[interval[[1]], wp];
iterList = {5, 12, 13, 14};
Plot @@ {f[x, #] & /@ iterList, {x}~Join~minmax,
WorkingPrecision -> wp, PlotLabel -> TraditionalForm[x == 6 Cos[x]],
PlotLegends -> (SF[" iterations", #] & /@ iterList)}


Traditional approach to "bound the damage" is to do interval arithmetic:

intervals = NestList[6 Cos[#] &, Interval[5.], 73];
points = Reverse[MinMax @@ #] & /@ intervals;
ListPlot[points\[Transpose], Filling -> {1 -> {2}},
PlotLegends -> {"upper", "lower"},
PlotLabel -> "Iterating 6cos(x) in 64-bit precision"]
ListPlot[Subtract @@@ points, PlotLabel -> "Range of correct answers",
AxesLabel -> {"iteration #", "range"}]

• Interestingly, the error on iterated 6cos(x) blows up between iterations 8 and 13; whereas the error on iterated cos(x) (without the factor 6) decreases with each iteration. I think this is because the sequence given by u(n) = 6 cos(u(n-1)) converges towards a region where cos has a high slope, whereas the sequence given by v(n) = cos(v(n-1)) converges towards a region where cos has a small slope.
– Stef
Commented Feb 17 at 12:24
• @Stef to useful theorems here are Banach's fixed point theorem and conditions for stable fixed point Commented Feb 17 at 21:08
• @Stef : the cosine iteration has a named limit, the "dottie number". // It is highly unlikely that the factor 6 sequence converges. Due to the finite number of floating point numbers in the interval the numerical sequence will be periodic, but likely with a very large period. It will go through all kinds of points, with high and low slopes. Commented Feb 18 at 10:16
• Interesting... does anyone know the reason why cos seems to have 2.5 to 3 bits error in most implementations? Is this a matter of trading off more terms for performence, or simply the order of magnitude different in the terms? Commented Feb 19 at 5:32

Intrigued by the answer by Lutz Lehmann, I gave Fortran with the gfortran compiler (version 13.2.0) a spin because for one it wasn't considered by the OP. For two, since 2008 the intrinsic iso_fortran_env module offers an explicit adjustment of the working precision (wp below) of floating numbers (reference):

program test
use iso_fortran_env, only: wp => real64
implicit none

real (kind=wp) :: u
integer :: i
u = 5.0_wp

do i = 1, 73
u = 6*cos(u)
print *, u
end do
end program test


with the following results

precision keyword storage result
single real32 32 bit -2.75651002
double real64 64 bit -3.8591995867613633

* In 2017, Steve Lionel informed

there’s one compiler where REAL128 is stored in 128 bits but is really the 80-bit “extended precision” real used in the old x87 floating point stack registers

to recommend an approach which then defines quadruple precision by

integer, parameter :: wp = selected_real_kind(33, 4931)


instead (reference, and a related discussion). In a test run, to drop the line iso_fortran_env in favor of this other approach however did not yield a different numeric result (and still does not suffice here).

• All more precision gives you is a different and longer random number, unfortunately. But I think between your Fortran and my Cobol programs, we've proved empirically what better mathematicians seem to be claiming above - that the results might as well be random. More precision doesn't seem to narrow in on an increasingly accurate result. You just get a different random number between -6 and +6. Commented Feb 18 at 15:49
• @CalebFuller From the answers especially by Lutz Lehmann, Yaroslav Bulatov, and davidlowryduda I retain passing replacing 32bit by 64bit, or 128bit precision did not improve the result per se about the 73rd iteration. Accuracy was lost at some point earlier -- when depending on the precision provided. Which let me wonder if instead of a numeric solution a symbolic one were (again?) a better one. But I have no working experience with Maple, SymPy etc. to substantiate this hunch with evidence. Commented Feb 18 at 16:05
• I actually tried it on an 8-bit computer running BASIC (emulation obviously). The first 10 results were about the same, the next 5 diverged slightly but less than 0.1, then result 16 onward it was suddenly a random grab-bag of numbers with no relation to the same Python code. Commented Feb 18 at 16:49

bc can calculate this to an arbitrary number of decimal places, and is installed on most Linux systems:

scale = 50+73
u = 5
for (i = 0; i < 73; i++) {
u = 6 * c(u)
}
u
quit


By saving this in a file a.bc and running bc -ql a.bc, we see:

\$ bc -ql a.bc
2.776758408482773464917250249724781630100266331685781980683517358426\
674759157161878144773316800820496120395277514821538061684


of which the first 50 digits should be accurate, assuming a digit of information loss per iteration.

• Nice. Wish I'd thought of that! Commented Feb 18 at 20:49

I ran it in Python 3.12.1 and the result is -3.7701774498987155. Interestingly, it is now identical to Julia until the very final digits. Seems there was an update to the maths libraries?

Seeing as Fortran has already been tried above, I thought I'd roll out Grandma Cobol to see what she made of it. Cobol has two interesting features. First, numbers are encoded as BCDs, not floating point binary. Second, you can specify EXACTLY how many significant digits come after the decimal. The implementation I used (https://onecompiler.com/cobol/) can take up to 38 digits so up to 37 after the decimal place.

The following code was duly fed into the punch card reader: 😜

identification division.
program-id. numerr.
data division.
working-storage section.
01 u pic s9v9(13) value 5.0 usage is computational-3.
01 i pic 99 value 0.
procedure division.
perform varying i from 0 by 1 until i > 72
compute u = 6 * function cos(u)
end-perform.
display u.
stop run.


The results are interesting, and basically confirms what some posters have said - that this algorithm might as well be a random number generator!

Merely by changing the number in brackets after the s9v9 on the 5th line, you adjust the decimal places "u" has, from 1 to 37. Every single time this number is changed, the answer changes completely. The only rule is that it is always a number in the range of -6 to +6, but within a single additional decimal place, the answer can jump wildly around. It is not moving towards a more accurate or refined answer, it is just random, wild jumps. All more precision is giving you is a longer random number!