ROOT  6.06/09
Reference Guide
sqrt.h
Go to the documentation of this file.
1 /*
2  * sqrt.h
3  * Implementations born on the Quake 3 fast inverse square root
4  * function.
5  * http://en.wikipedia.org/wiki/Fast_inverse_square_root
6  *
7  * Created on: Jun 24, 2012
8  * Author: Danilo Piparo, Thomas Hauth, Vincenzo Innocente
9  */
10 
11 /*
12  * VDT is free software: you can redistribute it and/or modify
13  * it under the terms of the GNU Lesser Public License as published by
14  * the Free Software Foundation, either version 3 of the License, or
15  * (at your option) any later version.
16  *
17  * This program is distributed in the hope that it will be useful,
18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20  * GNU Lesser Public License for more details.
21  *
22  * You should have received a copy of the GNU Lesser Public License
23  * along with this program. If not, see <http://www.gnu.org/licenses/>.
24  */
25 
26 #ifndef SQRT_H_
27 #define SQRT_H_
28 
29 #include "vdtcore_common.h"
30 
31 namespace vdt{
32 
33 //------------------------------------------------------------------------------
34 
35 
36 /// Sqrt implmentation from Quake3
37 inline double fast_isqrt_general(double x, const uint32_t ISQRT_ITERATIONS) {
38 
39  const double threehalfs = 1.5;
40  const double x2 = x * 0.5;
41  double y = x;
42  uint64_t i = details::dp2uint64(y);
43  // Evil!
44  i = 0x5fe6eb50c7aa19f9ULL - ( i >> 1 );
45  y = details::uint642dp(i);
46  for (uint32_t j=0;j<ISQRT_ITERATIONS;++j)
47  y *= threehalfs - ( x2 * y * y ) ;
48 
49  return y;
50 }
51 
52 //------------------------------------------------------------------------------
53 
54 /// Four iterations
55 inline double fast_isqrt(double x) {return fast_isqrt_general(x,4);}
56 
57 /// Three iterations
58 inline double fast_approx_isqrt(double x) {return fast_isqrt_general(x,3);}
59 
60 //------------------------------------------------------------------------------
61 
62 /// For comparisons
63 inline double isqrt (double x) {return 1./std::sqrt(x);}
64 
65 //------------------------------------------------------------------------------
66 
67 /// Sqrt implmentation from Quake3
68 inline float fast_isqrtf_general(float x, const uint32_t ISQRT_ITERATIONS) {
69 
70  const float threehalfs = 1.5f;
71  const float x2 = x * 0.5f;
72  float y = x;
73  uint32_t i = details::sp2uint32(y);
74  i = 0x5f3759df - ( i >> 1 );
75  y = details::uint322sp(i);
76  for (uint32_t j=0;j<ISQRT_ITERATIONS;++j)
77  y *= ( threehalfs - ( x2 * y * y ) );
78 
79  return y;
80 }
81 
82 //------------------------------------------------------------------------------
83 
84 /// Two iterations
85 inline float fast_isqrtf(float x) {return fast_isqrtf_general(x,2);}
86 
87 /// One (!) iterations
88 inline float fast_approx_isqrtf(float x) {return fast_isqrtf_general(x,1);}
89 
90 //------------------------------------------------------------------------------
91 
92 /// For comparisons
93 inline float isqrtf (float x) {return 1.f/std::sqrt(x);}
94 
95 //------------------------------------------------------------------------------
96 
97 // void isqrtv(const uint32_t size, double const * __restrict__ iarray, double* __restrict__ oarray);
98 // void fast_isqrtv(const uint32_t size, double const * __restrict__ iarray, double* __restrict__ oarray);
99 // void fast_approx_isqrtv(const uint32_t size, double const * __restrict__ iarray, double* __restrict__ oarray);
100 // void isqrtfv(const uint32_t size, float const * __restrict__ iarray, float* __restrict__ oarray);
101 // void fast_isqrtfv(const uint32_t size, float const * __restrict__ iarray, float* __restrict__ oarray);
102 // void fast_approx_isqrtfv(const uint32_t size, float const * __restrict__ iarray, float* __restrict__ oarray);
103 
104 } // end namespace vdt
105 
106 #endif /* SQRT_H_ */
double fast_isqrt_general(double x, const uint32_t ISQRT_ITERATIONS)
Sqrt implmentation from Quake3.
Definition: sqrt.h:37
float uint322sp(int x)
Converts an int to a float.
uint32_t sp2uint32(float x)
Converts a float to an int.
double sqrt(double)
static const double x2[5]
Double_t x[n]
Definition: legend1.C:17
double fast_approx_isqrt(double x)
Three iterations.
Definition: sqrt.h:58
float fast_isqrtf_general(float x, const uint32_t ISQRT_ITERATIONS)
Sqrt implmentation from Quake3.
Definition: sqrt.h:68
float fast_isqrtf(float x)
Two iterations.
Definition: sqrt.h:85
uint64_t dp2uint64(double x)
Converts a double to an unsigned long long.
double isqrt(double x)
For comparisons.
Definition: sqrt.h:63
double uint642dp(uint64_t ll)
Converts an unsigned long long to a double.
Double_t y[n]
Definition: legend1.C:17
float fast_approx_isqrtf(float x)
One (!) iterations.
Definition: sqrt.h:88
Definition: asin.h:32
float isqrtf(float x)
For comparisons.
Definition: sqrt.h:93
double fast_isqrt(double x)
Four iterations.
Definition: sqrt.h:55