The Double Pendulum Equation

This program solves the double pendulum equation using Euler's method. We integrate the differential equation over a time of 20 seconds, with 50000 Euler steps per second. We print the final value of the coordinates (x2,y2) of the end of the second pendulum. Here are the results obtained with different mathematical libraries (ordered by increasing value of x2):

libraryvalue of x2value of y2comments
icc 19.0.4.243-1.3434690.625392with -O1
RedHat's libm-1.2481480.609395
glibc 2.13-1.2481480.609395idem with 2.14-2.16.0
icc 19.0.4.243-1.1826200.601759with -O2/-O3
glibc-2.29-0.6547060.631662idem with 2.30-2.31
icc 19.0.4.243-0.6546940.631660with -O0
musl-0.6545670.631743
OpenLibm 0.4.1-0.6545670.631743
glibc-2.17-0.6545620.631744idem with 2.18-2.28
AMD libm 3.5.0-0.6545620.631744
GNU MPFR-0.6545620.631744correctly rounded

The program using MPFR is available here: it simply replaces the sinf and cosf by variants calling the corresponding MPFR functions.