Preperation for Class exam (C) -> n-th root with self-written functions (no math.h or something)

Question

the objective is easy:

find the geometrical middle of X1, X2, ... Xn (in my case n = 3)

BUT I have to write own functions and aren't able to use pow(), exp(), log2(), etc.

so I tried to calculate it with paper first, before I start coding. I used 125 as result of (a * b * c) because I know the 3rd root of 125 is 5

I thought to use "125 e(1/n)" But I'm really stuck at calculating this exp becouse I simply dont know how... Google isn't really helpfull though..

This is just a task to learn for the exam...

Answer 1

There is no general algorithm to find the n-th root of a number. But as we have a computer at hands we can use numerical approximations.

We have two generic methods to find the root of a monotonic function

• dichotomy: if you have one value below and one beyond, take the average and proceed. This is a very robust method but its convergence is slow
• Newton: You use the derivative formula to find a better value from an initial guess. This one is very fast, provided you start close enough to the root, but can be poor if you start with a bad guess. BTW we all know that the derivate for x -> x^n in x -> n * x^(n - 1)...

So a rule of thumb is to start with a dichotomy to find an acceptable guess and then go with Newton to find a very precise approximation

As we get a number of value we want to take the geometric mean, we know that the result is greater than the minimal value and less than the maximal one, so we have what we need to initialize a dichotomy.

A possible code could be:

#include <stdio.h>
#include <stdlib.h>
#include <float.h>

// a trivial implementation for x -> x^n
double ipow(double x, int n) {
    double val = 1;
    int i;
    for (i = 0; i < n; i++) val *= x;
    return val;
}

int main() {
    //double *arr;
    int i, n;
    printf("Enter the number of values: ");
    for (;;) {
        if (1 == scanf("%d", &n) && n > 1) break;
        printf("Invalid input, try again\n");
        int c;
        while ((c = fgetc(stdin)) != EOF && c != '\n');
        if (c == EOF) return 1;
    }
    double min = DBL_MAX, max = -DBL_MAX, a = 1;
    //arr = malloc(n * sizeof(*arr));
    printf("Enter %d values: ", n);
    for (i = 0; i < n; i++) {
        double val;
        for (;;) {
            if (1 == scanf("%lg", &val) && (val > 0)) break;
            printf("Invalid input, try again\n");
            int c;
            while ((c = fgetc(stdin)) != EOF && c != '\n');
            if (c == EOF) return 1;
        }
        a *= val;
        if (min > val) min = val;
        if (max < val) max = val;
    }
    // we want x, ipow(x, n) "close" to a, we know min <= x <= max
    // first a dichotomy
    double x, eps = 1 / 10.;
    int nd = 0, nn = 0; // will trace the number of dich and Newton iterations
    for (;;) {
        ++nd;
        x = (min + max) / 2;
        if (max - min < x * eps) {
            break;
        }
        if (ipow(x, n) < a) min = x;
        else max = x;
    }
    eps = 1e-10;
    // let's go with Newton
    for (;;) {
        ++nn;
        double x1 = x;
        x = x1 + (a - ipow(x1, n)) / n / ipow(x1, n - 1);
        if (x1 - x > -eps * x && x1 - x < eps * x) break;
    }
    printf("sqrt%d(%g) = %g (%d dichotomy, %d Newton)\n", n, a, x, nd, nn);
    return 0;
}

It is probably not the most efficient way, but it is robust because it is able to compute the geometric mean of 1e-300 and 1e300 while 1e300 * 1e300 overflows to inf... (but as fortunately inf is greater that any float value, the dichotomy part works fine)

Answer 2

Just like "abelenky" wrote in his (deleted) answer, you may also use a binary search if you want to calculate the inverse function of a strictly increasing function and you know that the result is positive:

Searching the n-th root of a number y is the inverse function of the function y=f(x, n); for x>0, this function is strictly increasing.

Because writing the C code is for your homework practice, I will not show working C code here but only explain you the principle:

Before doing the binary search, you have to perform two checks:

• If the number is zero, the solution (the root) is also zero.
• If the number is negative, you may calculate -root(-y, n) if n is odd; if n is even, there is no solution.
  (-root(-y, n) means: You calculate the n-th root of the absolute value of y and change the sign of the result.)

Now search the starting number g:

double g = 1;
while(f(g) > y) { g /= 2; }
while(f(g) < y) { g *= 2; }

In the case of the n-th root, f(x) is pow(x, n), which you can easily replace by a self-written function (using a for loop) if n is known to be a positive integer.

Then perform the binary search:

double x;
int i;
for(i=0; i<60; i++)
{
    if(f(x+g) <= y) { x+=g; }
    g /= 2;
}

Due to the precision of double variables, increasing the number of loop runs will not change the result. If you are using more precise data types (such as long double), you need more loop runs to get the best result possible.

Answer 3

You can use Newton method to compute nth root of x.

Newton method is

y = y - f(y) / f'(y)

with

f(y) = y^n - x

This gives the following iteration:

y = (n - 1) * y / n + x / n * y^(1-n)

For an initial y greater than x this sequence is convergent (see nth-root-iteration )

In plain C this gives:

#include <assert.h>
#include <stdio.h>

double nthpower(int n, double x)
{
  if (n < 0)
    return nthpower(-n, 1 / x);

  double y = 1;

  for (int i = 0; i < n; i++)
  {
    y = y * x;
  }

  return y;
}

int close_to_zero(double x)
{
  const double eps = 1e-10;

  return (-eps < x) && (x < eps);
}

double nthroot(int n, double x)
{
  assert(x >= 0);
  assert(n >= 0);

  switch (n)
  {
    case 0:
      return 1;

    case 1:
      return x;

    default:
      double yp, y = x;

      do
      {
        yp = y;
        y = (n - 1) * y / n + x / n * nthpower(1 - n, y);
      } while (!close_to_zero(yp - y));

      return y;
  }
}

double geometric_mean(double* x, int n)
{
  double p = 1;

  for (int i = 0; i < n; i++)
  {
    p *= x[i];
  }

  return nthroot(n, p);
}

int main(void)
{
  double x[6] = {2, 3, 4, 5, 6, 7};

  printf("GM %f\n", geometric_mean(x, 0));
  printf("GM %f\n", geometric_mean(x, 1));
  printf("GM %f\n", geometric_mean(x, 6));

  return 0;
}

which prints:

GM 1.000000

GM 2.000000

GM 4.140681

There is room for improvement, by example, nth power can be computed more effectively by computing x.x, then x^2.x^2, then x^4.x^4... But I think that the main idea (using the Newton method) is correctly illustrated.