diff options
Diffstat (limited to 'test/monniaux/math/test_gsl_pow.c')
-rw-r--r-- | test/monniaux/math/test_gsl_pow.c | 44 |
1 files changed, 44 insertions, 0 deletions
diff --git a/test/monniaux/math/test_gsl_pow.c b/test/monniaux/math/test_gsl_pow.c new file mode 100644 index 00000000..54e53889 --- /dev/null +++ b/test/monniaux/math/test_gsl_pow.c @@ -0,0 +1,44 @@ +#include <math.h> +#include <stdio.h> + +double gsl_pow_uint(double x, unsigned int n); + +double gsl_pow_int(double x, int n) +{ + unsigned int un; + + if(n < 0) { + x = 1.0/x; + un = -n; + } else { + un = n; + } + + return gsl_pow_uint(x, un); +} + +double gsl_pow_uint(double x, unsigned int n) +{ + double value = 1.0; + + /* repeated squaring method + * returns 0.0^0 = 1.0, so continuous in x + */ + do { + if(n & 1) value *= x; /* for n odd */ + n >>= 1; + x *= x; + } while (n); + + return value; +} + +int main() { + double y, y_expected; + for (int n = -9; n < 10; n++) { + y = gsl_pow_int (-3.14, n); + y_expected = pow (-3.14, n); + printf("%d %g %g\n", n, y, y_expected); + } + return 0; +} |