(*
 * ====================================================
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 *
 * Developed at SunPro, a Sun Microsystems, Inc. business.
 * Permission to use, copy, modify, and distribute this
 * software is freely granted, provided that this notice
 * is preserved.
 * ====================================================
 *
 * Optimized by Bruce D. Evans.
 *
 * Ported to pascal for single datatype by Jarl <slacky> Holta
 *)

const
  B1: UInt32 = 709958130; // (127-127/3-0.03306235651)*2^23
  B2: UInt32 = 696219795; // (127-127/3-54/3-0.03306235651)*2^23

  P0: Single =  1.87595182427177;  // Polynomial coefficients (single precision)
  P1: Single = -1.88497979543377;
  P2: Single =  1.62142972010535;
  P3: Single = -0.75839793477877;
  P4: Single =  0.14599619288661;

function cbrt(x: Single): Single; inline;
var
  hx, sign, approx: UInt32;
  t, r, s, w: Single;
  u: record
    case Boolean of
      True: (value: Single);
      False: (bits: UInt32);
  end;
begin
  u.value := x;
  hx := u.bits;

  sign := hx and UInt32($80000000);
  hx   := hx and UInt32($7FFFFFFF);

  if hx >= $7F800000 then Exit(x + x);
  if hx = 0 then          Exit(x);

  if hx < $00800000 then
  begin
    u.value := x * Single(1.6777216e+7);
    approx := u.bits;
    approx := approx div 3 + B2;
    u.bits := approx;
    t := u.value * Single(5.9604645e-8);
  end
  else
  begin
    // Compute the initial approximation for normal numbers
    approx := hx div 3 + B1;
    u.bits := sign or approx; // Combine sign and approximation
    t := u.value;
  end;

  r := (t * t) * (t / x);
  t := t * ((P0 + r * (P1 + r * P2)) + ((r * r) * r) * (P3 + r * P4));

  u.value := t;
  u.bits := (u.bits + UInt32($00400000)) and UInt32($FF800000);
  t := u.value;

  s := t * t;
  r := x / s;
  w := t + t;
  r := (r - t) / (w + r);
  t := t + t * r;

  Result := t;
end;