// Implement a floating-point cube root by iterating Newton's
// algorithm to a fixed point or periodic orbit.

// This is intended to be the simplest practically useful example of
// floating-point exact equality, for
// https://news.ycombinator.com/item?id=35883963.  On x87-equipped
// machines (not amd64) this code might hang unless compiled with
// -ffloat-store or -mfpmath=sse or equivalent.

// See also https://csclub.uwaterloo.ca/~pbarfuss/qbrt.pdf, "Computing
// a Real Cube Root", by Kahan, which gives not just the formula below
// but also more precise ways of computing the cube root.

// Note that some starting points, e.g., 29, will send this iteration
// into a periodic loop rather than finding a fixed point.  We do a
// marginal amount of extra computation to detect this.
#include <math.h>
#include <stdlib.h>
#include <stdio.h>


static inline float
m3form(float x, float y)
{
  float k = x*x*x;
  return x - (k - y)*x/(2*k + y);
}

float
cuberoot(float y)
{
  for (float hare, hop = y, tortoise = y;;) {
    for (size_t i = 0; i < 8; i++) {
      hare = hop;
      hop = m3form(hare, y);
      if (isnan(hop) || hop == hare || hop == tortoise) return hare;
    }

    tortoise = m3form(tortoise, y);
  }
}

#ifdef __GNUC__
#define weak __attribute__((weak))
#else
#define weak
#endif

weak int
main(int argc, char **argv)
{
  if (argc != 2) {
    fprintf(stderr, "Usage: %s 3.53\n(to compute ∛3.53)\n", argv[0]);
    return 1;
  }

  float y = atof(argv[1]);
  float x = cuberoot(y);
  printf("∛%g ≈ %g (err %g)\n", y, x, y-x*x*x);
  return 0;
}