test/monniaux/math/test_gsl_pow.c


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
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;
}