diff options
author | Erik Schnetter <schnetter@gmail.com> | 2012-11-30 15:50:03 -0500 |
---|---|---|
committer | Erik Schnetter <schnetter@gmail.com> | 2012-11-30 15:50:03 -0500 |
commit | d2614759a1d542c41af59b04d8711246d2a1e876 (patch) | |
tree | 2689209034d9ccc3d29208d50b82b4f89116aa1b /mathfuncs_rcp.h | |
download | vecmathlib-d2614759a1d542c41af59b04d8711246d2a1e876.zip vecmathlib-d2614759a1d542c41af59b04d8711246d2a1e876.tar.gz |
Import initial version
Diffstat (limited to 'mathfuncs_rcp.h')
-rw-r--r-- | mathfuncs_rcp.h | 48 |
1 files changed, 48 insertions, 0 deletions
diff --git a/mathfuncs_rcp.h b/mathfuncs_rcp.h new file mode 100644 index 0000000..1704d8f --- /dev/null +++ b/mathfuncs_rcp.h @@ -0,0 +1,48 @@ +// -*-C++-*- + +#ifndef MATHFUNCS_RCP_H +#define MATHFUNCS_RCP_H + +#include "mathfuncs_base.h" + +#include <cassert> +#include <cmath> + + + +namespace vecmathlib { + + template<typename realvec_t> + realvec_t mathfuncs<realvec_t>::vml_rcp(realvec_t x) + { + // Handle negative values + realvec_t x0 = x; + x = fabs(x); + + // Initial guess + assert(all(x > RV(0.0))); + intvec_t ilogb_x = ilogb(x); + // For stability, choose a starting value that is below the result + realvec_t r = scalbn(RV(0.5), -ilogb_x); + + // Iterate + int const nmax = 7; + for (int n=1; n<nmax; ++n) { + // Step + assert(all(x > RV(0.0))); + // Newton method: + // Solve f(r) = 0 for f(r) = x - 1/r + // r <- r - f(r) / f'(r) + // r <- 2 r - r^2 x + r *= RV(2.0) - r * x; + } + + // Handle negative values + r = copysign(r, x0); + + return r; + } + +}; // namespace vecmathlib + +#endif // #ifndef MATHFUNCS_RCP_H |