sqrt.h

Go to the documentation of this file.

/*
 * sqrt.h
 * Implementations born on the Quake 3 fast inverse square root
 * function.
 * http://en.wikipedia.org/wiki/Fast_inverse_square_root
 *
 *  Created on: Jun 24, 2012
 *      Author: Danilo Piparo, Thomas Hauth, Vincenzo Innocente
 */

/*
 * VDT is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser 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 Lesser Public License for more details.
 *
 * You should have received a copy of the GNU Lesser Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

#ifndef SQRT_H_
#define SQRT_H_

#include "vdtcore_common.h"

namespace vdt{

//------------------------------------------------------------------------------


/// Sqrt implmentation from Quake3
inline double fast_isqrt_general(double x, const uint32_t ISQRT_ITERATIONS) {

  const double threehalfs = 1.5;
  const double x2 = x * 0.5;
  double y  = x;
  uint64_t i  = details::dp2uint64(y);
  // Evil!
  i  = 0x5fe6eb50c7aa19f9ULL  - ( i >> 1 );
  y  = details::uint642dp(i);
  for (uint32_t j=0;j<ISQRT_ITERATIONS;++j)
      y *= threehalfs - ( x2 * y * y ) ;

  return y;
}

//------------------------------------------------------------------------------

/// Four iterations
inline double fast_isqrt(double x) {return fast_isqrt_general(x,4);}

/// Three iterations
inline double fast_approx_isqrt(double x) {return fast_isqrt_general(x,3);}

//------------------------------------------------------------------------------

/// For comparisons
inline double isqrt (double x) {return 1./std::sqrt(x);}

//------------------------------------------------------------------------------

/// Sqrt implmentation from Quake3
inline float fast_isqrtf_general(float x, const uint32_t ISQRT_ITERATIONS) {

   const float threehalfs = 1.5f;
   const float x2 = x * 0.5f;
   float y  = x;
   uint32_t i  = details::sp2uint32(y);
   i  = 0x5f3759df - ( i >> 1 );
   y  = details::uint322sp(i);
   for (uint32_t j=0;j<ISQRT_ITERATIONS;++j)
      y  *= ( threehalfs - ( x2 * y * y ) );

   return y;
}

//------------------------------------------------------------------------------

/// Two iterations
inline float fast_isqrtf(float x) {return fast_isqrtf_general(x,2);}

/// One (!) iterations
inline float fast_approx_isqrtf(float x) {return fast_isqrtf_general(x,1);}

//------------------------------------------------------------------------------

/// For comparisons
inline float isqrtf (float x) {return 1.f/std::sqrt(x);}

//------------------------------------------------------------------------------

// void isqrtv(const uint32_t size, double const * __restrict__ iarray, double* __restrict__ oarray);
// void fast_isqrtv(const uint32_t size, double const * __restrict__ iarray, double* __restrict__ oarray);
// void fast_approx_isqrtv(const uint32_t size, double const * __restrict__ iarray, double* __restrict__ oarray);
// void isqrtfv(const uint32_t size, float const * __restrict__ iarray, float* __restrict__ oarray);
// void fast_isqrtfv(const uint32_t size, float const * __restrict__ iarray, float* __restrict__ oarray);
// void fast_approx_isqrtfv(const uint32_t size, float const * __restrict__ iarray, float* __restrict__ oarray);

} // end namespace vdt

#endif /* SQRT_H_ */