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

/* ----------------------------- acbrt ------------------------------ */

/* This is a novel and fast routine for the approximate cube root of
an IEEE float (single precision). It is derived from a similar program
for the approximate square root.
   The relative error ranges from 0 to +0.00103.
   Caution:
Result for -0 is NaN.
Result for 0 is 1.24e-13.
For denorms it is either within tolerance or gives a result < 2.1e-13.
Gives the correct result (inf) for x = inf.
Gives the correct result (NaN) for x = NaN. */

float acbrt(float x) {
   float y;
   int i = *(int *)&x;          // View x as an int.
   i = i/4 + i/16;              // Approximate divide by 3.
   i = i + i/16;
   i = i + i/256;
   i = 0x2a5137a0 + i;          // Initial guess.
   y = *(float *)&i;            // View i as float.
   y = 0.33333333f*(2.0f*y + x/(y*y));  // Newton step.
   return y;
}

/* ----------------------------- acbrt1 ----------------------------- */

/* This is acbrt with an additional step of the Newton iteration, for
increased accuracy.
   The relative error ranges from 0 to +0.00000116. */

float acbrt1(float x) {
   float y;
   int i = *(int *)&x;          // View x as an int.
   i = i/4 + i/16;              // Approximate divide by 3.
   i = i + i/16;
   i = i + i/256;
   i = 0x2a5137a0 + i;          // Initial guess.
   y = *(float *)&i;            // View i as float.
   y = 0.33333333f*(2.0f*y + x/(y*y));  // Newton step.
   y = 0.33333333f*(2.0f*y + x/(y*y));  // Newton step again.
   return y;
}

/* ----------------------------- acbrt2 ----------------------------- */

/* This is a very approximate but very fast version of acbrt. It is just
two integer instructions (shift right and divide), plus instructions
to load the constant.
   The constant 0x2a51067f balances the relative error at +-0.0316. */

float acbrt2(float x) {
   int i = *(int *)&x;          // View x as an int.
   i = 0x2a51067f + i/3;        // Initial guess.
   x = *(float *)&i;            // View i as float.
   return x;
}

/* ---------------------------- acbrt2a ----------------------------- */

/* This is a very approximate but very fast version of acbrt. It is just
eight integer instructions (shift rights and adds), plus instructions
to load the constant.
   1/3 is approximated as 1/4 + 1/16 + 1/64 + 1/256 + ... + 1/65536.
   The constant 0x2a511cd0 balances the relative error at +-0.0321. */

float acbrt2a(float x) {
   int i = *(int *)&x;          // View x as an int.
   i = i/4 + i/16;              // Approximate divide by 3.
   i = i + i/16;
   i = i + i/256;
   i = 0x2a511cd0 + i;          // Initial guess.
   x = *(float *)&i;            // View i as float.
   return x;
}

/* ---------------------------- acbrt2b ----------------------------- */

/* This is a very approximate but very fast version of acbrt. It is just
six integer instructions (shift rights and adds), plus instructions
to load the constant.
   1/3 is approximated as 1/4 + 1/16 + 1/64 + 1/256.
   The constant 0x2a6497f8 balances the relative error at +-0.151. */

float acbrt2b(float x) {
   int i = *(int *)&x;          // View x as an int.
   i = i/4 + i/16;              // Approximate divide by 3.
   i = i + i/16;
   i = 0x2a6497f8 + i;          // Initial guess.
   x = *(float *)&i;            // View i as float.
   return x;
}

/* ------------------------------ main ------------------------------ */

int main(int argc, char *argv[]) {
   float x, yt, ya;

   if (argc != 2) {
      printf("Need exactly one argument, a floating-point number.\n");
      return 1;
   }

   x = atof(argv[1]);
   yt = cbrtf(x);               // From the math library.
   ya = acbrt(x);               // Change to acbrt1 or acbrt2 etc. to
                                // test the other routines.
   printf("x = %g %08x\nTrue cbrt = %g, approx = %g, rel error = %.7f\n",
      x, *(int *)&x, yt, ya, (ya - yt)/yt);
   return 0;
}