/* orbit_double.c with checkpoints, for the divergence figure: the same loop body,
   run in stretches between checkpoints twenty to a decade from 1 to 10^8 steps,
   printing the state and the digest at each one. Build it the four ways the
   table does; the last line of each equals that build's orbit_double.c result. */
#include <stdio.h>
#include <stdint.h>
#include <string.h>
#include <math.h>
int main(void) {
    const double GM = 398600441800000.0;
    double x = 6778137.0, y = 0.0, vx = 0.0, vy = 7668.558;
    double rmin = 1e300, rmax = 0.0;
    uint64_t digest = 0;
    long i = 0;
    for (int k = 0; k <= 160; k++) {
        long check = (long)llround(pow(10.0, k / 20.0));
        if (check <= i) continue;
        for (; i < check; i++) {
            double r2 = x * x + y * y;
            double r = sqrt(r2);
            double g = GM / r2;
            double ax = -g * x / r;
            double ay = -g * y / r;
            vx += ax; vy += ay;
            x += vx; y += vy;
            if (r < rmin) rmin = r;
            if (r > rmax) rmax = r;
            uint64_t bx, by, bvx, bvy;
            memcpy(&bx, &x, 8); memcpy(&by, &y, 8); memcpy(&bvx, &vx, 8); memcpy(&bvy, &vy, 8);
            digest = (digest * 1099511628211ULL) ^ bx ^ (by << 21) ^ (bvx << 42) ^ (bvy << 7);
        }
        printf("%ld %.17g %.17g %.17g %.17g %llu\n", i, x, y, vx, vy, (unsigned long long)digest);
    }
    return 0;
}