cbrt.h

Go to the documentation of this file.

/*******************************************************************************
 * This file is part of SWIFT.
 * Copyright (c) 2016 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
 *                    Matthieu Schaller (schaller@strw.leidenuniv.nl)
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published
 * by the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 ******************************************************************************/
#ifndef SWIFT_CBRT_H
#define SWIFT_CBRT_H
 
/* Config parameters. */
#include <config.h>
 
/* Some standard headers. */
#include <math.h>
 
/* Local headers. */
#include "inline.h"
 
__attribute__((always_inline)) INLINE static float icbrtf(float x_in) {
 
  union {
    float as_float;
    unsigned int as_uint;
    int as_int;
  } cast;
 
  /* Extract the exponent. */
  cast.as_float = x_in;
  const int exponent = ((cast.as_int & 0x7f800000) >> 23) - 127;
 
  /* Clear the exponent and sign to get the mantissa. */
  cast.as_uint = (cast.as_uint & ~0xff800000) | 0x3f800000;
  const float x_norm = cast.as_float;
 
  /* Multiply by sqrt(1/2) and subtract one, should then be in the
     range [sqrt(1/2) - 1, sqrt(2) - 1). */
  const float x = x_norm * (float)M_SQRT1_2 - 1.0f;
 
  /* Compute the polynomial interpolant. */
  float res =
      9.99976591940035e-01f +
      x * (-3.32901212909283e-01f +
           x * (2.24361110929912e-01f +
                x * (-1.88913279594895e-01f + x * 1.28384036492344e-01f)));
 
  /* Compute the new exponent and the correction factor. */
  int exponent_new = exponent;
  if (exponent_new < 0) exponent_new -= 2;
  exponent_new = -exponent_new / 3;
  const int exponent_rem = exponent + 3 * exponent_new;
  cast.as_uint = (exponent_new + 127) << 23;
  static const float scale[3] = {8.90898718140339e-01f, 7.07106781186548e-01f,
                                 5.61231024154687e-01f};
  const float exponent_scale = cast.as_float * scale[exponent_rem];
 
  /* Scale the result and set the correct sign. */
  res = copysignf(res * exponent_scale, x_in);
 
  /* One step of Newton iteration to refine the result. */
  res *= (1.0f / 3.0f) * (4.0f - x_in * res * res * res);
 
  /* We're done. */
  return res;
}
 
#endif /* SWIFT_CBRT_H */