C 언어로 평방근, 입방근, n제곱근 구하는 함수를 구현하고 테스트하기

음이 아닌 실수 A 의 평방근 sqrt(A) 를 구하는 Heron 의 방법:

        반복함수  g(x) = (x + A/x) / 2   를 이용

 

실수 A 의 n제곱근 root(n, A) 를 구하는 Newton-Raphson 의 방법

        반복함수  g(x) = ((n-1)*x + A/(x**(n - 1))) / n    를 이용

n = 2 인 경우에는 Newton-Raphson 의 방법이 Heron 의 방법과 동일하다.

(참조. http://en.wikipedia.org/wiki/Newton's_method )

 

(math.h 를 인클루드하여 사용하는) C 언어의 표준 수학 라이브러리에 지수 계산 함수 pow() 가 이미 구현되어 있다. 하지만 차후 필요한 데가 있을 것 같아서 이와 유사한 n 제곱 함수와 n 제곱근 함수를 구현해 보았다.

지수가 정수인 거듭제곱을 계산하는  함수도 nPow(), gPow, mPow() 세 개 구현해 놓았는데, 이들 세 함수는 절차적 언어의 성능상 재귀호출이 아니고 단순 반복 기법을 사용하는 함수이다. 이 세 함수 중 mPow() 의 성능이 가장 우수하다. 큰 지수의 경우 for 반복문의 반복회수를 따져 보면 성능 비교를 할 수 있을 것이다. (성능 비교를 위해 세 가지를 모두 소스에 남겨 두었다.) mPow() 함수는 n 제곱근을 구하는 재귀함수 newtonNthRoot(int, double) 의 구현에 사용되기도 한다. if ... else ... 구문이 많아 소스가 복잡하게 보일지 모르겠으나 이는 밑수나 지수가 음수이거나 0인 경우의 처리를 위함이다. 구현된 모든 함수의 구현에는 예외상황(예를 들어, 음수의 짝수 제곱근 같은 예외상황) 처리 과정이 있다.

아래의 소스는 Visual C++ 외에 gcc 로 컴파일해도 수정 없이 잘 되리라고 본다.

예외상황 처리룰 하기 위해, C 언어에서는 sejmp.h 를 인클루드하고, 예외상황을 던질 때는 longjmp() 함수를, 예외상황을 받을 때는 setjmp() 함수와 switch{ ... } 문을 사용한다.

// Filename: testNthRoot.c
//
//            Approximate square roots, cubic roots and n-th roots of a given number.
//
// Compile: cl testNthRoot.c
// Execute: testNthRoot
//
// Date: 2013. 1. 6.
// Copyright (c) 2013 PH Kim  (pkim __AT__ scripts.pe.kr)

#include <stdio.h>
#include <math.h>
#include <setjmp.h>

jmp_buf g_stkBuf;

#define NO_EXCEPTION          (0)
#define CANNOT_EVALUATE    (100)

#define MAX_ITER 20000
#define M_EPSILON 1.0e-15

 

/**
  * Compute the n-th root of x to a given scale, x > 0.
  */
double nPow(double a, int n)
{
    int i;
    double y;

    if (n > 0) {
        if (n == 1)
            return a;
        else {
            if (a == 0.0 || a == 1.0) {
                return a;
            }
            else if (a == -1.0) {
                if (n % 2 == 1)
                    return -1.0;
                else
                    return 1.0;
            }
            else if (a < 0.0) {
                if (n % 2 == 1)
                    return -nPow(-a, n);
                else
                    return nPow(-a, n);
            }
            else {
                y = 1.0;
                for (i = 0; i < n; i++) {
                    y *= a;
                }
                return y;
            }
        }
    }
    else if (n == 0) {
        return 1.0;
    }
    else {      //  when n < 0
        if (a == 0.0) {
            printf("Negative powering exception of zero.\n");
            longjmp( g_stkBuf, CANNOT_EVALUATE );  // setjmp로 설정한 g_stkBuf(스택정보)를 가지고 해당 시점으로 복귀
            return -1.0;
        }
        else {
            if (n == -1)
                return 1.0/a;
            else
                return 1.0/nPow(a, -n);
        }
    }
}

 

/**
  * Compute the n-th root of x to a given scale, x > 0.
  */
double gPow(double a, int n)
{
    int i;
    double y;
    double r;
    int m;
    int one;

    if (n > 0) {
        if (n == 1)
            return a;
        else {
            if (a == 0.0 || a == 1.0) {
                return a;
            }
            else if (a == -1.0) {
                if (n % 2 == 1)
                    return -1.0;
                else
                    return 1.0;
            }
            else if (a < 0.0) {
                if (n % 2 == 1)
                    return -gPow(-a, n);
                else
                    return gPow(-a, n);
            }
            else {

                y = 1.0;
                r = a;
                m = 8*sizeof(int) - 1;    // 8*sizeof(int) - 1;  ignore sign bit, which is MSB
                one = 1;
                for (i = 0; i < m; i++) {
                    if ((n & one) == 0) {
                        y *= 1.0;
                    }
                    else {
                        y *= r;
                    }
                    r = r*r;
                    one <<= 1;
                    if (one > n)
                        break;
                }
                return y;
            }
        }
    }
    else if (n == 0) {
        return 1.0;
    }
    else {      //  when n < 0
        if (a == 0.0) {
            printf("Negative powering exception of zero.\n");
            longjmp( g_stkBuf, CANNOT_EVALUATE );  // setjmp로 설정한 g_stkBuf(스택정보)를 가지고 해당 시점으로 복귀
            return -1.0;
        }
        else {
            if (n == -1)
                return 1.0/a;
            else
                return 1.0/gPow(a, -n);
        }
    }
}

 

/**
  * Compute the n-th root of x to a given scale, x > 0.
  */
double mPow(double a, int n)
{
    if (n > 0) {
        if (n == 1)
            return a;
        else {
            if (a == 0.0 || a == 1.0) {
                return a;
            }
            else if (a == -1.0) {
                if (n % 2 == 1)
                    return -1.0;
                else
                    return 1.0;
            }
            else if (a < 0.0) {
                if (n % 2 == 1)
                    return -mPow(-a, n);
                else
                    return mPow(-a, n);
            }
            else {

                double y = 1.0;
                double r = a;
                int m = n;
                while (m > 0) {
                    if ((m & 0x1) != 0) {
                        y *= r;
                    }
                    r = r*r;
                    m >>= 1;
                }
                return y;
            }
         }
    }
    else if (n == 0) {
        return 1.0;
    }
    else {      //  when n < 0
        if (a == 0.0) {
            printf("Negative powering exception of zero.\n");
            longjmp( g_stkBuf, CANNOT_EVALUATE );  // setjmp로 설정한 g_stkBuf(스택정보)를 가지고 해당 시점으로 복귀
            return -1.0;
        }
        else {
            if (n == -1)
                return 1.0/a;
            else
                return 1.0/mPow(a, -n);
        }
    }
}

 

/**
  * Compute the square root of x to a given scale, x > 0.
  */
double heronSqrt(double a)
{
    double x1, x2, er;
    int counter;

    if (a < 0.0) {
        printf("Cannot find the sqrt of a negative number.\n");
        longjmp( g_stkBuf, CANNOT_EVALUATE );  // setjmp로 설정한 g_stkBuf(스택정보)를 가지고 해당 시점으로 복귀
                                                 // errNo는 setjmp로 설정한 지점으로 복귀할 때 리턴 값
        return -1.0;
    }
    else if (a == 0.0 || a == 1.0) {
        return a;
    }
    else {
        x1 = a;
        x2 = (x1 + a/x1)/2.0;
        er = x1 - x2;
        counter = 0;
        while (x1 + er != x1) {
            x1 = x2;
            x2 = (x1 + a/x1)/2.0;
            er = x1 - x2;
            if (fabs(er) < fabs(M_EPSILON*x1))
                break;
            counter++;
            if (counter > MAX_ITER)
                break;
        }
        if (counter >= MAX_ITER) {
            printf("Inaccurate sqrt exception by too many iterations.\n");
            longjmp( g_stkBuf, CANNOT_EVALUATE );  // setjmp로 설정한 g_stkBuf(스택정보)를 가지고 해당 시점으로 복귀
            return -1.0;
        }
        return x2;
    }
}

/**
  * Compute the cubic root of x to a given scale, x > 0.
  */
double newtonCbrt(double a)
{
    double x1, x2, er;
    int counter;

    if (a == 0.0 || a == 1.0 || a == -1.0) {
        return a;
    }
    else if (a < 0.0) {
        return -newtonCbrt(-a);
    }
    else {
        x1 = a;
        x2 = (2.0*x1 + a/(x1*x1))/3.0;
        er = x1 - x2;
        counter = 0;
        while (x1 + er != x1) {
            x1 = x2;
            x2 = (2.0*x1 + a/(x1*x1))/3.0;
            er = x1 - x2;
            if (fabs(er) < fabs(M_EPSILON*x1))
                break;
            counter++;
            if (counter > MAX_ITER)
                break;
        }
        if (counter >= MAX_ITER) {
            printf("Inaccurate cbrt exception by too many iterations.\n");
            longjmp( g_stkBuf, CANNOT_EVALUATE );  // setjmp로 설정한 g_stkBuf(스택정보)를 가지고 해당 시점으로 복귀
            return -1.0;
        }
        return x2;
    }
}

/**
  * Compute the n-th root of x to a given scale, x > 0.
  */
double newtonNthRoot(int n, double a)
{
    double x1, x2, xn, er;
    int counter;

    if (n == 0) {
        return 1.0;
    }
    else if (n == 1) {
        return a;
    }
    else if (n > 0) {
        if (a == 0.0 || a == 1.0) {
            return a;
        }
        else if (a == -1.0) {
            if (n % 2 == 1)
                return a;
            else {
                printf("Cannot find the even n-th root of a negative number.\n");
                longjmp( g_stkBuf, CANNOT_EVALUATE );  // setjmp로 설정한 g_stkBuf(스택정보)를 가지고 해당 시점으로 복귀
                return -1.0;
            }
        }
        else if (a < 0.0) {
            if (n % 2 == 1)
                return -newtonNthRoot(n, -a);
            else {
                printf("Cannot find the even n-th root of a negative number.\n");
                longjmp( g_stkBuf, CANNOT_EVALUATE );  // setjmp로 설정한 g_stkBuf(스택정보)를 가지고 해당 시점으로 복귀
                return -1.0;
            }
        }
        else if (a < 1.0) {
            return 1.0/newtonNthRoot(n, 1.0/a);
        }
        else {
            x1 = a;
            xn = mPow(x1, n - 1);
            x2 = ((n - 1)*x1 + a/xn)/n;
            er = x1 - x2;
            counter = 0;
            while (x1 + er != x1) {
                x1 = x2;
                xn = mPow(x1, n - 1);
                x2 = ((n - 1)*x1 + a/xn)/n;
                er = x1 - x2;
                if (fabs(er) < fabs(M_EPSILON*x1))
                    break;
                counter++;
                if (counter > MAX_ITER)
                    break;
            }
            if (counter >= MAX_ITER) {
                printf("Inaccurate n-th exception by too many iterations.\n");
                longjmp( g_stkBuf, CANNOT_EVALUATE );  // setjmp로 설정한 g_stkBuf(스택정보)를 가지고 해당 시점으로 복귀
                return -1.0;
            }
            return x2;
        }
    }
    else {
        if (a == 0.0) {
            printf("Cannot find the negative n-th root of zero.\n");
            longjmp( g_stkBuf, CANNOT_EVALUATE );  // setjmp로 설정한 g_stkBuf(스택정보)를 가지고 해당 시점으로 복귀
            return -1.0;
        }
        else {
            return 1.0/newtonNthRoot(-n, a);
        }
    }
}

 

int main()
{
    double x = 16.0;
    double u = sqrt(x);

    int k;
    double t1, t2, t3;
    double y, z, w;

    int errNo = NO_EXCEPTION;

    printf("[ Testing heronSqrt(double) ]--------------------\n");
    printf("x = %g\n", x);
    printf("u = sqrt(%g) = %g\n", x, u);
    y = heronSqrt(x);
    printf("y = heronSqrt(%g) = %g\n", x, y);
    printf("y*y = %g\n", y*y);
    printf("\n");

    printf("[ Testing newtonCbrt(double) ]--------------------\n");
    x = -216.0;
    printf("x = %g\n", x);
    printf("-exp(log(-x)/3.0) = %g\n", -exp(log(-x)/3.0));
    w = newtonCbrt(x);
    printf("w = newtonCbrt(%g) = %g\n", x, w);
    printf("w*w*w = %g\n", w*w*w);
    printf("\n");

    x = 729000000000.0;
    printf("x = %g\n", x);
    printf("exp(log(x)/3.0) = %g\n", exp(log(x)/3.0));
    w = newtonCbrt(x);
    printf("w = newtonCbrt(%g) = %g\n", x, w);
    printf("w*w*w = %g\n", w*w*w);
    printf("\n");

    printf("[ Testing newtonNthRoot(int, double) ]--------------------\n");
    z = newtonNthRoot(3, x);
    printf("x = %g\n", x);
    printf("z = newtonNthRoot(3, %g) = %g\n", x, z);
    printf("z*z*z = %g\n", z*z*z);
    printf("\n");

    x = 12960000000000000000.0;
    z = newtonNthRoot(4, x);
    printf("x = %g\n", x);
    printf("z = newtonNthRoot(4, x) = newtonNthRoot(4, %g) = %g\n", x, z);
    printf("z*z*z*z = %g\n", z*z*z*z);
    printf("\n");

    x = 1.0/12960000000000000000.0;
    z = newtonNthRoot(4, x);
    printf("x = %g\n", x);
    printf("exp(log(x)/4.0) = %g\n", exp(log(x)/4.0));
    printf("z = newtonNthRoot(4, x) = newtonNthRoot(4, %g) = %g\n", x, z);
    printf("z*z*z*z = %g\n", z*z*z*z);
    printf("\n");

    errNo = setjmp( g_stkBuf );  // longjmp로 복귀할 지점을 설정, g_stkBuf에 스택정보 저장
    switch ( errNo )
    {
        case NO_EXCEPTION:
            x = -4.0;
            printf("[ Test Exception heronSqrt(double) ]--------------------\n");
            printf("x = %g\n", x);
            printf("Calculating heronSqrt(%g)\n" , x);
            y = heronSqrt(x);
            printf("y = heronSqrt(%g) = %g\n", x, y);
            printf("y*y = %g\n", y*y);
            printf("\n");
            break;
        case CANNOT_EVALUATE:
            printf("Caught some exception in calculating heronSqrt(%g)\n", x);
            printf("\n");
            errNo = NO_EXCEPTION;
            break;
    }

    errNo = setjmp( g_stkBuf );  // longjmp로 복귀할 지점을 설정, g_stkBuf에 스택정보 저장
    switch ( errNo )
    {
        case NO_EXCEPTION:
            x = -4.0;
            printf("[ Test Exception in newtonCbrt(double) ]--------------------\n");
            printf("x = %g\n", x);
            printf("Calculating newtonCbrt(%g)\n", x);
            y = newtonCbrt(x);
            printf("y = newtonCbrt(%g) = %g\n", x, y);
            printf("y*y*y = %g\n", y*y*y);
            printf("\n");
            break;
        case CANNOT_EVALUATE:
            printf("Caught some exception in calculating newtonCbrt(%g)\n", x);
            printf("\n");
            errNo = NO_EXCEPTION;
            break;
    }

    printf("[ Test calculations by powering ]-----------------------------\n");
    x = 200.0;
    z = newtonNthRoot(10, x);
    printf("x = %g\n", x);
    printf("exp(log(x)/10.0) = %g\n", exp(log(x)/10.0));
    printf("z = newtonNthRoot(10, x) = newtonNthRoot(10, %g) = %g\n", x,  z);
    printf("pow(z, 10) = %g\n", pow(z, 10));
    printf("\n");

    x = 3001.0;
    z = newtonNthRoot(99, x);
    printf("x = %g\n", x);
    printf("exp(log(x)/99.0) = %g\n", exp(log(x)/99.0));
    printf("z = newtonNthRoot(99, x) = newtonNthRoot(99, %g) = %g\n", x, z);
    printf("pow(z, 99) = %g\n", pow(z, 99));
    printf("\n");

    x = 3001.0;
    z = newtonNthRoot(-99, x);
    printf("x = %g\n", x);
    printf("exp(log(x)/-99.0) = %g\n", exp(log(x)/-99.0));
    printf("z = newtonNthRoot(-99, x) = newtonNthRoot(-99, %g) = %g\n", x, z);
    printf("1.0/pow(z, 99) = %g\n", 1.0/pow(z, 99));
    printf("\n");

    printf("2.1**2.1 = pow(2.1, 2.1) = %g\n", pow(2.1, 2.1));
    printf("2.1**(-2.1) = pow(2.1, -2.1) = %g\n", pow(2.1, -2.1));
    printf("2.1**2.1 * 2.1**(-2.1) = pow(2.1, 2.1) * pow(2.1, -2.1) = %g\n", pow(2.1, 2.1)*pow(2.1, -2.1));
    printf("2.1**2.1 = exp(2.1*log(2.1)) = %g\n", exp(2.1*log(2.1)));
    printf("2.1**(-2.1) = exp(-2.1*log(2.1)) = %g\n", exp(-2.1*log(2.1)));
    printf("2.1**2.1 * 2.1**(-2.1) = exp(2.1*log(2.1)) * exp(-2.1*log(2.1)) = %g\n", exp(2.1*log(2.1)) * exp(-2.1*log(2.1)));
    printf("\n");

    k = 301;
    x = -1.029;
    t1 = nPow(x, k);
    t2 = gPow(x, k);
    t3 = mPow(x, k);
    printf("t1 = nPow(%g, %d) = %g\n", x, k, t1);
    printf("t2 = gPow(%g, %d) = %g\n", x, k, t2);
    printf("t3 = mPow(%g, %d) = %g\n", x, k, t3);
    printf("t1 / t2 = %g\n", t1 / t2);
    printf("t1 - t2 = %g\n", t1 - t2);
    printf("t1 == t2 ? %s\n",  t1 == t2 ? "yes" : "no");
    printf("t1 / t3 = %g\n", t1 / t3);
    printf("t1 - t3 = %g\n", t1 - t3);
    printf("t1 == t3 ? %s\n",  t1 == t3 ? "yes" : "no");
    printf("t2 / t3 = %g\n", t2 / t3);
    printf("t2 - t3 = %g\n", t2 - t3);
    printf("t2 == t3 ? %s\n",  t2 == t3 ? "yes" : "no");
    printf("\n");

    printf("Done.\n");

}

/*
Output:
[ Testing heronSqrt(double) ]--------------------
x = 16
u = sqrt(16) = 4
y = heronSqrt(16) = 4
y*y = 16

[ Testing newtonCbrt(double) ]--------------------
x = -216
-exp(log(-x)/3.0) = -6
w = newtonCbrt(-216) = -6
w*w*w = -216

x = 7.29e+011
exp(log(x)/3.0) = 9000
w = newtonCbrt(7.29e+011) = 9000
w*w*w = 7.29e+011

[ Testing newtonNthRoot(int, double) ]--------------------
x = 7.29e+011
z = newtonNthRoot(3, 7.29e+011) = 9000
z*z*z = 7.29e+011

x = 1.296e+019
z = newtonNthRoot(4, x) = newtonNthRoot(4, 1.296e+019) = 60000
z*z*z*z = 1.296e+019

x = 7.71605e-020
exp(log(x)/4.0) = 1.66667e-005
z = newtonNthRoot(4, x) = newtonNthRoot(4, 7.71605e-020) = 1.66667e-005
z*z*z*z = 7.71605e-020

[ Test Exception heronSqrt(double) ]--------------------
x = -4
Calculating heronSqrt(-4)
Cannot find the sqrt of a negative number.
Caught some exception in calculating heronSqrt(-4)

[ Test Exception in newtonCbrt(double) ]--------------------
x = -4
Calculating newtonCbrt(-4)
y = newtonCbrt(-4) = -1.5874
y*y*y = -4

[ Test calculations by powering ]-----------------------------
x = 200
exp(log(x)/10.0) = 1.69865
z = newtonNthRoot(10, x) = newtonNthRoot(10, 200) = 1.69865
pow(z, 10) = 200

x = 3001
exp(log(x)/99.0) = 1.08424
z = newtonNthRoot(99, x) = newtonNthRoot(99, 3001) = 1.08424
pow(z, 99) = 3001

x = 3001
exp(log(x)/-99.0) = 0.922308
z = newtonNthRoot(-99, x) = newtonNthRoot(-99, 3001) = 0.922308
1.0/pow(z, 99) = 3001

2.1**2.1 = pow(2.1, 2.1) = 4.74964
2.1**(-2.1) = pow(2.1, -2.1) = 0.210542
2.1**2.1 * 2.1**(-2.1) = pow(2.1, 2.1) * pow(2.1, -2.1) = 1
2.1**2.1 = exp(2.1*log(2.1)) = 4.74964
2.1**(-2.1) = exp(-2.1*log(2.1)) = 0.210542
2.1**2.1 * 2.1**(-2.1) = exp(2.1*log(2.1)) * exp(-2.1*log(2.1)) = 1

t1 = nPow(-1.029, 301) = -5457.93
t2 = gPow(-1.029, 301) = -5457.93
t3 = mPow(-1.029, 301) = -5457.93
t1 / t2 = 1
t1 - t2 = 6.18456e-011
t1 == t2 ? no
t1 / t3 = 1
t1 - t3 = 6.18456e-011
t1 == t3 ? no
t2 / t3 = 1
t2 - t3 = 0
t2 == t3 ? yes

Done.
*/