diff --git a/src/decn/proto/CMakeLists.txt b/src/decn/proto/CMakeLists.txt index 4790df0..9fc63e1 100644 --- a/src/decn/proto/CMakeLists.txt +++ b/src/decn/proto/CMakeLists.txt @@ -7,3 +7,6 @@ target_link_libraries(ln mpfr) add_executable(exp exp.cpp) target_link_libraries(exp mpfr) + +add_executable(sqrt recip_sqrt.cpp) +target_link_libraries(sqrt mpfr) diff --git a/src/decn/proto/recip_sqrt.cpp b/src/decn/proto/recip_sqrt.cpp new file mode 100644 index 0000000..0a8d41b --- /dev/null +++ b/src/decn/proto/recip_sqrt.cpp @@ -0,0 +1,138 @@ +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU 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 General Public License +// along with this program. If not, see . + + +//prototype for calculating reciprocal square root 1/sqrt(x) +//based on 0x5F3759DF fast inverse square root: +// calculate using newton-raphson iterations, +// get initial estimate using linear approximation + +#include +#include +#include +#include +using namespace boost::multiprecision; + +// #define DEBUG + +using std::cout; +using std::endl; + + +static const unsigned int CALC_PRECISION = 18; + +int get_exp(mpfr_float& x){ + mpfr_float x_exp(log10(x), CALC_PRECISION); + + return x_exp.convert_to(); +} + +double get_significand(mpfr_float& x){ + int exp = get_exp(x); + double signif = x.convert_to(); + signif = signif / pow(10.0, exp); + + return signif; +} + + +int main(void){ + cout << std::scientific << std::setprecision(CALC_PRECISION); + + //loop through values to test +#ifdef DEBUG + mpfr_float x(0.1, CALC_PRECISION); + { +#else + for (mpfr_float x(1e-99, CALC_PRECISION); x < 1e99; x *= 1.03){ +#endif + //build mpf values + mpfr_float x_orig_2(x/2, CALC_PRECISION); //copy of original x / 2 + mpfr_float est(0, CALC_PRECISION); //current estimate for 1/sqrt(x) + mpfr_float accum(0, CALC_PRECISION); //working register + mpfr_float_1000 actual; + actual = 1 / sqrt(x); + + //get decimal floating point representation + int x_exp = get_exp(x); + double x_signif = get_significand(x); + //normalize: make sure significand is between 1 and 10 + while (x_signif < 1.0){ + x_signif *= 10; + x_exp--; + } +#ifdef DEBUG + printf("x = %f*10^%d\n", x_signif, x_exp); +#endif + //get initial estimate for 1/sqrt(x) == 10^(-0.5 * log(x)): + // approximate significand == 10^(-0.5 * log(x_signif)) + // with linear approximation: -0.18 * x_signif + 2.5 + // new exponent part is (10^(-0.5 * log(10^x_exp))) + // == 10^(-0.5 * x^exp) + double est_signif; + if (x_exp & 0x1){ //odd + //increment x_exp and + //approximate estimate significand as (-0.056*x_signif + 0.79) * 10^0.5 + // == -0.18 * x_signif + 2.5 + x_exp++; + est_signif = -0.18 * x_signif + 2.5; + } else { //even + //keep x_exp as is and approximate estimate significand as + // -0.056*x_signif + 0.79 + est_signif = -0.056 * x_signif + 0.79; + } + int est_exp = -x_exp / 2; + est = est_signif * pow(10.0, est_exp); +#ifdef DEBUG + printf("x: %f (%d)\n", x_signif, x_exp); + printf("est: %f (%d)\n", est_signif, est_exp); + cout << "initial estimate for " << x << ": " << est << endl; +#endif + + //newton-raphson iterations + for (int i = 0; i < 6; i++){ + accum = est * est; + accum *= x_orig_2; //accum = x/2 * est * est + accum = -accum; + accum += 1.5; //accum = 3/2 - x/2 * est * est + accum *= est; //accum = 0.5 * est * (3 - x * est * est) + est = accum; +#ifdef DEBUG + printf("%2d: ", i); + cout << est << endl; +#endif + } + + + //calculate relative error + mpfr_float_1000 calc, diff; + calc = accum; + diff = abs((actual - calc)/actual); + +#ifdef DEBUG_ALL + if (1){ +#else + if (diff > 5e-17){ +#endif + cout << x << ": "; + cout << std::setprecision(18) << accum << ", "; +#if 1 + cout << actual << ", "; +#endif + cout << diff << endl; + } + } + + return 0; +} + diff --git a/src/decn/proto/recip_sqrt.py b/src/decn/proto/recip_sqrt.py new file mode 100644 index 0000000..4be28b7 --- /dev/null +++ b/src/decn/proto/recip_sqrt.py @@ -0,0 +1,24 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Wed Dec 11 00:56:01 2019 + +@author: jeffrey +""" + +import numpy as np +import matplotlib.pyplot as plt + +xvals = np.arange(1, 10, 0.01) +yvals = 10 ** (-0.5 * np.log10(xvals)) +#linear fit for significand of 1/sqrt(x) +coef = np.polyfit(xvals, yvals, 1) +poly1d_fn = np.poly1d(coef) + +plt.plot(xvals, yvals, 'b') +plt.plot(xvals, poly1d_fn(xvals), '--k') +plt.xlabel('significand') +plt.ylabel('1/sqrt(significand)') +plt.show() +print("linear approx coeff:", coef) +print(" * 3.16:", coef * 3.16227766) \ No newline at end of file