Scala 언어로 평방근, 입방근, 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 )

 

Scala 언어에는 지수 계산을 위한 함수 scala.math.pow(밑수, 지수) 가 이미 구현되어 있다. 하지만 차후 필요한 데가 있을 것 같아서 이와 유사한 n 제곱 함수와 n 제곱근 함수를 구현해 보았다.

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

 

1. Scala 언어는 JVM(자바가상기계) 위에서 동작하는 (절차적 프로그래밍도 지원하는) 함수형 언어이다.

   하지만, 아래의 소스에서는 함수형 언어의 기능은 거의 사용하지 않앗고, 절차적 언어의 기능을 위주로 사용하였다. 하나의 함수를 구현하는데 함수형 기능과 절차적 기능을 섞어서 사용하는 것은 좋지 않은 습관이다.

2. Scala  언어는 타입에 동적인(dynamic) 언어이다. 변수 선언시에 타입을 지정하지 않아도 된다.

3. Scala  언어에서는 val 로 선언 된 것은 상수, var 로 선언된 것은 변수이다.

     val a = 2.7      // a 상수
     var b = 5.3     // b 는 변수 

아래의 소스 첫 부분에

        val MAX_ITER = 20000
        val M_EPSILON = 1.0e-15
       

라고 선언하였으니, MAX_ITER 와 M_EPSILON 는 상수로 선언되었다.

4. Scala 언어에는 return 이라는 예약어가 없다. 그냥 리턴될 값만 수식으로 표현해 놓으면, Scala 컴파일러(또는 Scala 인터프리터)가 리턴될 값을 알아서 인식하고 처리해 준다.

5. 예외상황 처리를 위해 예외 던지기 구문 throw ... 와 예외 받기 구문 try ... catch ...을 이용하였다.

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

import java.io._

object testNthRoot {

    val MAX_ITER = 20000
    val M_EPSILON = 1.0e-15
   
   
    //
    // Compute the n-th root of x to a given scale, x > 0.
    //
    def nPow(a: Double, n: Int) : Double = {
        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 {
                    var y = 1.0
                    for (i <- 1 to n) {
                        y = y * a
                    }
                    return y
                }
            }
        }
        else if (n == 0)
            return 1.0
        else {     //  when n < 0
            if (a == 0.0) {
                throw new RuntimeException("Negative powering exception of zero.")
            }
            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.
    //
    def gPow(a: Double, n: Int) : Double = {
        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 {
                    var y = 1.0
                    var r = a
                    var m = 8*4 - 1            //  8*sizeof(int) - 1;
                    var one = 1
                    for (i <- 0 to m) {
                        if ((n & one) == 0)
                            y = y * 1.0
                        else {
                            y = y * r
                        }
                        r = r*r
                        one = (one << 1)
                        if (one > n) {
                            return y
                        }
                    }
                    return y
                }
            }
        }
        else if (n == 0)
            return 1.0
        else {     //  when n < 0
            if (a == 0.0) {
                throw new RuntimeException("Negative powering exception of zero.")
            }
            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.
    //
    def mPow(a: Double, n: Int) : Double = {
        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 {
                    var y = 1.0
                    var r = a
                    var m = n
                    while (m > 0) {
                        if ((m & 0x1) == 1) {
                            y = y * r
                        }
                        r = r*r
                        m = (m >> 1)
                    }
                    return y
                }
            }
        }
        else if (n == 0)
            return 1.0
        else {     //  when n < 0
            if (a == 0.0) {
                throw new RuntimeException("Negative powering exception of zero.")
            }
            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.
    //
    def heronSqrt(a: Double) : Double = {
        if (a < 0.0) {
            throw new RuntimeException("Cannot find the sqrt of a negative number.")
        }
        else if (a == 0.0 || a == 1.0)
            a
        else {
            var x1 = a
            var x2 = (x1 + a/x1)/2.0
            var er = x1 - x2
            var counter = 0
            var not_stop = true
            while (x1 + er != x1 && not_stop) {
                x1 = x2
                x2 = (x1 + a/x1)/2.0
                er = x1 - x2
                if (scala.math.abs(er) < scala.math.abs(M_EPSILON*x1)) {
                    not_stop = false
               }
                counter = counter + 1
                if (counter > MAX_ITER) {
                    // break
                    not_stop = false
                }
            }
            if (counter >= MAX_ITER) {
                    throw new RuntimeException("Inaccurate sqrt exception by too many iterations.")
            }
            x2
        }
    }

    //
    // Compute the cubic root of x to a given scale, x > 0.
    //
    def newtonCbrt(a: Double) : Double = {
        if (a == 0.0 || a == 1.0 || a == -1.0 )
            a
        else if (a < 0.0) {
            return -newtonCbrt(-a)
        }
        else {
            var x1 = a
            var x2 = (2.0*x1 + a/(x1*x1))/3.0
            var er = x1 - x2
            var counter = 0
            var not_stop = true
            while (x1 + er != x1 && not_stop) {
                x1 = x2
                x2 = (2.0*x1 + a/(x1*x1))/3.0
                er = x1 - x2
                if (scala.math.abs(er) < scala.math.abs(M_EPSILON*x1)) {
                    not_stop = false
                }
                counter = counter + 1
                if (counter > MAX_ITER) {
                    not_stop = false
                }
            }
            if (counter >= MAX_ITER) {
                throw new RuntimeException("Inaccurate cbrt exception by too many iterations.")
            }
            return x2
        }
    }

    //
    // Compute the n-th root of x to a given scale, x > 0.
    //
    def newtonNthRoot(n: Int, a: Double) : Double = {
        if (n == 0)
            1.0
        else if (n == 1)
            a
        else if (n > 0) {
            if (a == 0.0 || a == 1.0)
                a
            else if (a == -1.0) {
                if ((n % 2) == 1)
                    a
                else {
                    throw new RuntimeException("Cannot find the even n-th root of a negative number.")
                }
            }
            else if (a < 0.0) {
                if ((n % 2) == 1)
                    -newtonNthRoot(n, -a)
                else {
                    throw new RuntimeException("Cannot find the even n-th root of a negative number.")
                }
            }
            else if (a < 1.0)
                1.0/newtonNthRoot(n, 1.0/a)
            else {
                var x1 = a
                var xn = mPow(x1, n - 1)
                var x2 = ((n - 1)*x1 + a/xn)/n
                var er = x1 - x2
                var counter = 0
                var not_stop = true
                while (x1 + er != x1 && not_stop) {
                    x1 = x2
                    xn = mPow(x1, n - 1)
                    x2 = ((n - 1)*x1 + a/xn)/n
                    er = x1 - x2
                    if (scala.math.abs(er) < scala.math.abs(M_EPSILON*x1)) {
                        not_stop = false
                    }
                    counter = counter + 1
                    if (counter > MAX_ITER) {
                        not_stop = false
                    }
                }
                if (counter >= MAX_ITER) {
                    throw new RuntimeException("Inaccurate n-th root exception by too many iterations.")
                }
                x2
            }
        }
        else {
            if (a == 0.0) {
                throw new RuntimeException("Cannot find the negative n-th root of zero.")
            }
            else
                1.0/newtonNthRoot(-n, a)
        }
    }

 

    def main(args: Array[String]) {
            var x : Double = 16
            var u = scala.math.sqrt(x)

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

            x = 729000000000.0
            printf("x = %g\n", x)
            printf("exp(log(x)/3.0) = %g\n", scala.math.exp(scala.math.log(x)/3.0))
            var 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")
            var 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", scala.math.exp(scala.math.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")

            try {
                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")
            }
            catch {
                case ex : RuntimeException =>
                        printf("%s\nCaught some exception in calculating heronSqrt(%g)\n", ex.getMessage(),  x);
                        printf("\n")
            }

            try {
                x = -4.0
                printf("[ Test Exception 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")
            }
            catch {
                case ex : RuntimeException =>
                        printf("%s\nCaught some exception in calculating newtonCbrtrt(%g)\n", ex.getMessage(),  x);
                        printf("\n")
            }

            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",  scala.math.exp(scala.math.log(x)/10.0))
            printf("z = newtonNthRoot(10, x) = newtonNthRoot(10, %g) = %g\n", x, z)
            printf("pow(z, 10.0) = %g\n", scala.math.pow(z, 10.0))
            printf("\n")

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

            x = 3001.0
            z = newtonNthRoot(-99, x)
            printf("x = %g\n", x)
            printf("exp(log(x)/-99.0) = %g\n", scala.math.exp(scala.math.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/scala.math.pow(z, 99.0))
            printf("\n")

            printf("2.1**2.1 = pow(2.1, 2.1) = %g\n", scala.math.pow(2.1, 2.1))
            printf("2.1**(-2.1) = pow(2.1, -2.1) = %g\n", scala.math.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",  scala.math.pow(2.1, 2.1) * scala.math.pow(2.1, -2.1))
            printf("2.1**2.1 = exp(2.1*log(2.1)) = %g\n", scala.math.exp(2.1*scala.math.log(2.1)))
            printf("2.1**(-2.1) = exp(-2.1*log(2.1)) = %g\n", scala.math.exp(-2.1*scala.math.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", scala.math.pow(2.1, 2.1) * scala.math.pow(2.1, -2.1))
            printf("\n")

            var k : Int = 301
            x = -1.029
            var t1 = nPow(x, k)
            var t2 = gPow(x, k)
            var t3 = mPow(x, k)
            var tt =  scala.math.pow(x, k)
            printf("scala.math.sqrt(%g, %d) = %g\n", x, k, tt)
            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 ? " )
            if (t1 == t2) printf("yes\n") else printf("no\n")
            printf("t1 / t3 = %g\n", t1 / t3)
            printf("t1 - t3 = %g\n", t1 - t3)
            printf("t1 == t3 ? " )
            if (t1 == t3) printf("yes\n") else printf("no\n")
            printf("t2 / t3 = %g\n", t2 / t3)
            printf("t2 - t3 = %g\n", t2 - t3)
            printf("t2 == t3 ? " )
            if (t2 == t3) printf("yes\n") else printf("no\n")
            printf("\n")

            printf("Done.\n")
    }

}

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

x = 7.29000e+11
exp(log(x)/3.0) = 9000.00
w = newtonCbrt(7.29000e+11) = 9000.00
w*w*w = 7.29000e+11

[ Testing newtonNthRoot(int, double) ]--------------------
x = 7.29000e+11
z = newtonNthRoot(3, 7.29000e+11) = 9000.00
z*z*z = 7.29000e+11

x = 1.29600e+19
z = newtonNthRoot(4, x) = newtonNthRoot(4, 1.29600e+19) =  60000.0
z*z*z*z = 1.29600e+19

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

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

[ Test Exception newtonCbrt(double) ]--------------------
x = -4.00000
Calculating newtonCbrt(-4.00000)
y = newtonCbrt(-4.00000) = -1.58740
y*y*y = -4.00000

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

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

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

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.00000
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.00000

scala.math.sqrt(-1.02900, 301) = -5457.93
t1 = nPow(-1.02900, 301) = -5457.93
t2 = gPow(-1.02900, 301) = -5457.93
t3 = mPow(-1.02900, 301) = -5457.93
t1 / t2 = 1.00000
t1 - t2 = 6.18456e-11
t1 == t2 ? no
t1 / t3 = 1.00000
t1 - t3 = 6.18456e-11
t1 == t3 ? no
t2 / t3 = 1.00000
t2 - t3 = 0.00000
t2 == t3 ? yes

Done.
*/