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

 

Julia 언어에는 지수 연산자 ^ 를 (밑수)^(지수) 의 형식으로 언어 자체에서 지원하고 있다. 하지만 차후 필요한 데가 있을 것 같아서 이와 유사한 n 제곱 함수와 n 제곱근 함수를 구현해 보았다. (참고로 Julia 언어에는 pow() 함수가 정의되어 있지 않다.)

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

 

#!/usr/bin/env julia

# Filename: testNthRoot.jl
#
#            Approximate square roots, cubic roots and n-th roots of a given number.
#
# Execute: julia testNthRoot.jl
#
# Compile: ipy %IRONPYTHON_HOME%\tools\scripts\pyc.py /main:testNthRoot.py /target:exe
# Execute: testNthRoot
#
# Date: 2013. 3. 7.
# Copyright (c) 2013 PH Kim  (pkim __AT__ scripts.pe.kr)

MAX_ITER = 20000
M_EPSILON = 1.0e-15

#
# Compute the n-th root of x to a given scale, x > 0.
#
function nPow(a, n)
    if n > 0
        if n == 1
            return a
        else
            if a == 0.0 || a == 1.0
                return a
            elseif a == -1.0
                if n % 2 == 1
                    return -1.0
                else
                    return 1.0
                end
            elseif a < 0.0
                if n % 2 == 1
                    return -nPow(-a, n)
                else
                    return nPow(-a, n)
                end
            else
                y = 1.0
                for i = 1:n
                    y *= a
                end
                return y
            end
        end
    elseif n == 0
        return 1.0
    else      #  when n < 0
        if a == 0.0
            println("Error")   # throws "Negative powering exception of zero."
        else
            if n == -1
                return 1.0/a
            else
                return 1.0/nPow(a, -n)
            end
        end
    end
end

#
# Compute the n-th root of x to a given scale, x > 0.
#
function gPow(a, n)
    if n > 0
        if n == 1
            return a
        else
            if a == 0.0 || a == 1.0
                return a
            elseif a == -1.0
                if n % 2 == 1
                    return -1.0
                else
                    return 1.0
                end
            elseif a < 0.0
                if n % 2 == 1
                    return -gPow(-a, n)
                else
                    return gPow(-a, n)
                end
            else
                y = 1.0
                r = a
                m = 8*4 - 1            #  8*sizeof(int) - 1;
                one = 1
                for i = 0:m
                    if (n & one) == 0
                        y *= 1.0
                    else
                        y *= r
                    end
                    r = r*r
                    one <<= 1
                    if one > n
                        break
                    end
                end
                return y
            end
        end
    elseif n == 0
        return 1.0
    else      #  when n < 0
        if a == 0.0
            throw("Negative powering exception of zero.")
        else
            if n == -1
                return 1.0/a
            else
                return 1.0/gPow(a, -n)
            end
        end
    end
end

#
# Compute the n-th root of x to a given scale, x > 0.
#
function mPow(a, n)
    if n > 0
        if n == 1
            return a
        else
            if a == 0.0 || a == 1.0
                return a
            elseif a == -1.0
                if n % 2 == 1
                    return -1.0
                else
                    return 1.0
                end
            elseif a < 0.0
                if n % 2 == 1
                    return -mPow(-a, n)
                else
                    return mPow(-a, n)
                end
            else
                y = 1.0
                r = a
                m = n
                while m > 0
                    if (m & 0x1) == 1
                        y *= r
                    end
                    r = r*r
                    m >>= 1
                end
                return y
            end
        end
    elseif n == 0
        return 1.0
    else      #  when n < 0
        if a == 0.0
            throw("Negative powering exception of zero.")
        else
            if n == -1
                return 1.0/a
            else
                return 1.0/mPow(a, -n)
            end
        end
    end
end

#
# Compute the square root of x to a given scale, x > 0.
#
function heronSqrt(a)
    if a < 0.0
        throw("Cannot find the sqrt of a negative number.")
    elseif 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 abs(er) < abs(M_EPSILON*x1)
                break
            end
            counter += 1
            if counter > MAX_ITER
                break
            end
        end
        if counter >= MAX_ITER
            throw("Inaccurate sqrt exception by too many iterations.")
        end
        return x2
    end
end

#
# Compute the cubic root of x to a given scale, x > 0.
#
function newtonCbrt(a)
    if a == 0.0 || a == 1.0 || a == -1.0
        return a
    elseif a < 0.0
        return -newtonCbrt(-a)
    else
        x1 = a
        x2 = (2.0*x1 + a/(x1*x1))/3.0
        # println("    x2 = $x2")
        er = x1 - x2
        counter = 0
        while x1 + er != x1
            x1 = x2
            x2 = (2.0*x1 + a/(x1*x1))/3.0
            # println("        x2 = $x2")
            er = x1 - x2
            # println("        abs(er) = $(abs(er))")
            # println("        abs(M_EPSILON*x1) = $(abs(M_EPSILON*x1))")
            if abs(er) < abs(M_EPSILON*x1)
                break
            end
            counter += 1
            if counter > MAX_ITER
                break
            end
        end
        if counter >= MAX_ITER
            throw("Inaccurate cbrt exception by too many iterations.")
        end
        return x2
    end
end

#
# Compute the n-th root of x to a given scale, x > 0.
#
function newtonNthRoot(n, a)
    if n == 0
        return 1.0
    elseif n == 1
        return a
    elseif n > 0
        if a == 0.0 || a == 1.0
            return a
        elseif a == -1.0
            if n % 2 == 1
                return a
            else
                throw("Cannot find the even n-th root of a negative number.")
            end
        elseif a < 0.0
            if n % 2 == 1
                return -newtonNthRoot(n, -a)
            else
                throw("Cannot find the even n-th root of a negative number.")
            end
        elseif 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 abs(er) < abs(M_EPSILON*x1)
                    break
                end
                counter += 1
                if counter > MAX_ITER
                    break
                end
            end
            if counter >= MAX_ITER
                throw("Inaccurate n-th root exception by too many iterations.")
            end
            return x2
        end
    else
        if a == 0.0
            throw("Cannot find the negative n-th root of zero.")
        else
            return 1.0/newtonNthRoot(-n, a)
        end
    end
end

x = 16.0
u = sqrt(x)

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

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

x = 729000000000.0
@printf("x = %f\n", x)
@printf("exp(log(x)/3.0) = %f\n", exp(log(x)/3.0))
w = newtonCbrt(x)
@printf("w = newtonCbrt(%f) = %f\n", x, w)
@printf("w*w*w = %f\n", w*w*w)
println()

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

x = 12960000000000000000.0
z = newtonNthRoot(4, x)
println("x = $x")
println("z = newtonNthRoot(4, x) = newtonNthRoot(4, $x) = $z")
@printf("z*z*z*z = %f\n", z*z*z*z)
println()

x = 1.0/12960000000000000000.0
z = newtonNthRoot(4, x)
println("x = $x")
@printf("exp(log(x)/4.0) = %f\n", exp(log(x)/4.0) )
println("z = newtonNthRoot(4, x) = newtonNthRoot(4, $x) = $z")
println("z*z*z*z = $(z*z*z*z)" )
println()

try
        x = -4.0
        println("[ Test Exception heronSqrt(double) ]--------------------")
        println("x = $x")
        println("Calculating heronSqrt($x)")
        y = heronSqrt(x)
        println("y = heronSqrt($x) = $y")
        println("y*y = $(y*y)")
        println()
catch ex
        println("$ex\nCaught some exception in calculating heronSqrt($x)")
        println()
end

try
        x = -4.0
        println("[ Test Exception in newtonCbrt(double) ]--------------------")
        println("x = $x")
        println("Calculating newtonCbrt($x)")
        y = newtonCbrt(x)
        println("y = newtonCbrt($x) = $y")
        println("y*y*y = $(y*y*y)")
        println()
catch ex
        println("$ex\nCaught some exception in calculating newtonCbrt$x)")
        println()
end

println("[ Test calculations by powering ]-----------------------------")
x = 200.0
z = newtonNthRoot(10, x)
println("x = $x")
println("exp(log(x)/10.0) = $(exp(log(x)/10.0))")
println("z = newtonNthRoot(10, x) = newtonNthRoot(10, $x) = $z")
println("z^10 = $(z^10)")
println()

x = 3001.0
z = newtonNthRoot(99, x)
println("x = $x")
println("exp(log(x)/99.0) = $(exp(log(x)/99.0))")
println("z = newtonNthRoot(99, x) = newtonNthRoot(99, $x) = $z")
println("z^99 = $(z^99)")
println()

x = 3001.0
z = newtonNthRoot(-99, x)
println("x = $x")
println("exp(log(x)/-99.0) = $(exp(log(x)/-99.0))")
println("z = newtonNthRoot(-99, x) = newtonNthRoot(-99, $x) = $z")
println("1.0/z^99 = $(1.0/(z^99))")
println("z^(-99) = $(z^(-99))")
println()

println("2.1^2.1 = $(2.1^ 2.1)")
println("2.1^2.1 * 2.1^(-2.1) = $(2.1^2.1 * 2.1^(-2.1))")
println("2.1^2.1 = exp(2.1*log(2.1)) = $(exp(2.1*log(2.1)))")
println("2.1^(-2.1) = exp(-2.1*log(2.1)) = $(exp(-2.1*log(2.1)))")
println("2.1^2.1 * 2.1^(-2.1) = exp(2.1*log(2.1)) * exp(-2.1*log(2.1)) = $(exp(2.1*log(2.1))) * $(exp(-2.1*log(2.1))) = $( exp(2.1*log(2.1)) * exp(-2.1*log(2.1)) )")
println()

k = 301
x = -1.029
t1 = nPow(x, k)
t2 = gPow(x, k)
t3 = mPow(x, k)
println("$x^$k = $(x^k)")
println("t1 = nPow($x, $k) = $t1")
println("t2 = gPow($x, $k) = $t2")
println("t3 = mPow($x, $k) = $t3")
println("t1 / t2 = $(t1 / t2)")
println("t1 - t2 = $(t1 - t2)")
print("t1 == t2 ? ")
if t1 == t2
    print("yes")
else
    print("no")
end
println()
println("t1 / t3 = $(t1 / t3)")
println("t1 - t3 = $(t1 - t3)")
print("t1 == t3 ? ")
if t1 == t3
    print("yes")
else
    print("no")
end
println()
println("t2 / t3 = $(t2 / t3)")
println("t2 - t3 = $(t2 - t3)")
print("t2 == t3 ? ")
if t2 == t3
    print("yes")
else
    print("no")
end
println()
println()
println("Done.")

"""
[ Testing heronSqrt(double) ]--------------------
x = 16.000000
u = sqrt(16.000000) = 4.000000
y = heronSqrt(16.000000) = 4.000000
y*y = 16.000000

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

x = 729000000000.000000
exp(log(x)/3.0) = 9000.000000
w = newtonCbrt(729000000000.000000) = 9000.000000
w*w*w = 729000000000.000000

[ Testing newtonNthRoot(int, double) ]--------------------
x = 729000000000.000000
z = newtonNthRoot(3, 729000000000.000000) = 9000.000000
z*z*z = 729000000000.000000

x = 1.296e19
z = newtonNthRoot(4, x) = newtonNthRoot(4, 1.296e19) = 60000.0
z*z*z*z = 12960000000000000000.000000

x = 7.716049382716049e-20
exp(log(x)/4.0) = 0.000017
z = newtonNthRoot(4, x) = newtonNthRoot(4, 7.716049382716049e-20) = 1.6666666666
666667e-5
z*z*z*z = 7.716049382716051e-20

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

[ Test Exception in newtonCbrt(double) ]--------------------
x = -4.0
Calculating newtonCbrt(-4.0)
y = newtonCbrt(-4.0) = -1.5874010519681994
y*y*y = -3.999999999999999

[ Test calculations by powering ]-----------------------------
x = 200.0
exp(log(x)/10.0) = 1.6986464646342472
z = newtonNthRoot(10, x) = newtonNthRoot(10, 200.0) = 1.6986464646342472
z^10 = 199.9999999999999

x = 3001.0
exp(log(x)/99.0) = 1.0842361893258805
z = newtonNthRoot(99, x) = newtonNthRoot(99, 3001.0) = 1.0842361893258805
z^99 = 3000.999999999995

x = 3001.0
exp(log(x)/-99.0) = 0.9223082662659932
z = newtonNthRoot(-99, x) = newtonNthRoot(-99, 3001.0) = 0.9223082662659932
1.0/z^99 = 3001.000000000005
z^(-99) = 3001.0000000000045

2.1^2.1 = 4.749638091742242
2.1^2.1 * 2.1^(-2.1) = 0.9999999999999999
2.1^2.1 = exp(2.1*log(2.1)) = 4.749638091742242
2.1^(-2.1) = exp(-2.1*log(2.1)) = 0.21054235726688478
2.1^2.1 * 2.1^(-2.1) = exp(2.1*log(2.1)) * exp(-2.1*log(2.1)) = 4.74963809174224
2 * 0.21054235726688478 = 1.0

-1.029^301 = -5457.928015771622
t1 = nPow(-1.029, 301) = -5457.92801577163
t2 = gPow(-1.029, 301) = -5457.928015771692
t3 = mPow(-1.029, 301) = -5457.928015771692
t1 / t2 = 0.9999999999999887
t1 - t2 = 6.184563972055912e-11
t1 == t2 ? no
t1 / t3 = 0.9999999999999887
t1 - t3 = 6.184563972055912e-11
t1 == t3 ? no
t2 / t3 = 1.0
t2 - t3 = 0.0
t2 == t3 ? yes

Done.
"""