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

Question

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.

Answer 1

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.

Answer 2

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.

Answer 3

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.

Answer 4

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).

Answer 5

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.

Answer 6

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"}]

Answer 7

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
quadruple*	real128	128 bit	1.94361921524798815846937773902061750

---

* 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).

Answer 8

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.

Answer 9

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!