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

 

Visual BASIC 언어에는 System 모듈에 지수 계산 함수 Math.Pow(double, double) 가 이미 구현되어 있다. 하지만 차후 필요한 데가 있을 것 같아서 이와 유사한 n 제곱 함수와 n 제곱근 함수를 구현해 보았다.

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

아래의 소스는 대부분 버전의 비쥬얼 스튜디오의 Visual BASIC 컴파일러로 컴파일 되고 실행되게 작성된 소스이다.

소스 첫 부분에

    Private Const MAX_ITER As Integer = 20000
    Private Const M_EPSILON As Double = 1.0e-15

라고 선언하였으니 변수 MAX_ITER 와 M_EPSILON 는 (타입을 갖는) 상수로 선언되었다. Java 언어로는 상수를 선언할 방법이 없지만 Visual BASIC 언어로는 이와 같이 Const 예약어(키워드)를 이용하여 상수를 선언할 수 있다.

예외상황 처리를 위해 예외 던지기 구문 Throw ... 과 예외 받기 구문 Try ... Catch ... End Try 구문을 사용히였다. 그리고 여러 줄을 주석문으로 처리하기 위해 (Visual Basic 구문은 아니고) Visual Basic 컴파일러에게 지시하는 구문이지만

#If False Then
..............................................................
..... 이곳에 여러 줄 주석을 달 수 있음 .......
..............................................................
#End If

를 사용하였다.

 

REM ============================================================================
REM Filename: TestNthRootApp.vb
REM
REM            Approximate square roots, cubic roots and n-th roots of a given number.
REM
REM Compile: vbc TestNthRootApp.vb
REM Execute: TestNthRootApp
REM
REM Date: 2013. 1. 9.
REM Copyright (c) 2013 PH Kim  (pkim __AT__ scripts.pe.kr)
REM ============================================================================

Module ApproximateModule

    Private Const MAX_ITER As Integer = 20000
    Private Const M_EPSILON As Double = 1.0e-15

    Private Const CRLF As String = Chr(13) & Chr(10)

    '
    ' Compute the n-th root of x to a given scale, x > 0.
    '
    Public Function nPow(a As Double, n As Integer) As Double
        Dim y As Double

        If (n > 0) Then
            If (n = 1) Then
                return a
            Else
                If (a = 0.0 Or a = 1.0) Then
                    return a
                ElseIf (a = -1.0) Then
                    If (n Mod 2 = 1) Then
                        return -1.0
                    Else
                        return 1.0
                    End If
                ElseIf (a < 0.0) Then
                    If (n Mod 2 = 1)Then
                        return -nPow(-a, n)
                    Else
                        return nPow(-a, n)
                    End If
                Else
                    y = 1.0
                    For i As Integer = 1 to n
                        y = y * a
                    Next
                    return y
                End If
            End If
        ElseIf (n = 0) Then
            return 1.0
        Else       '  when n < 0
            If (a = 0.0) Then
                Throw New Exception("Negative powering exception of zero.")
            Else
                If (n = -1) Then
                    return 1.0/a
                Else
                    return 1.0/nPow(a, -n)
                End If
            End If
        End If
    End Function

 

    '
    ' Compute the n-th root of x to a given scale, x > 0.
    '
    Public Function gPow(a As Double, n As Integer) As Double
        Dim  y As Double
        Dim r As Double
        Dim one As Integer
        Dim m As Integer

        If (n > 0) Then
            If (n = 1) Then
                return a
            Else
                If (a = 0.0 Or a = 1.0) Then
                    return a
                ElseIf (a = -1.0) Then
                    If (n Mod 2 = 1) Then
                        return -1.0
                    Else
                        return 1.0
                    End If
                ElseIf (a < 0.0) Then
                    If (n Mod 2 = 1) Then
                        return -gPow(-a, n)
                    Else
                        return gPow(-a, n)
                    End If
                Else
                    y = 1.0
                    r = a
                    m = 8*4 - 1           '  8*sizeof(int) - 1;
                    one = 1
                    For i As Integer = 1 To m
                        If ((n And one) = 0) Then
                            y = y * 1.0
                        Else
                            y = y * r
                        End If
                        r = r*r
                        one = one << 1
                        if (one > n) Then
                            Exit For
                        End if
                    Next
                    return y
                End If
            End If
        ElseIf (n = 0) Then
            return 1.0
        Else       '  when n < 0
            If (a = 0.0) Then
                Throw New Exception("Negative powering exception of zero.")
            Else
                if (n = -1) Then
                    return 1.0/a
                Else
                    return 1.0/gPow(a, -n)
                End if
            End If
        End If
    End Function

    '
    ' Compute the n-th root of x to a given scale, x > 0.
    '
    Public Function mPow(a As Double, n As Integer) As Double
        Dim  y As Double
        Dim r As Double
        Dim m As Integer

        If (n > 0) Then
            If (n = 1) Then
                return a
            Else
                If (a = 0.0 Or a = 1.0)  Then
                    return a
                ElseIf (a = -1.0)  Then
                    If (n Mod 2 = 1) Then
                        return -1.0
                    Else
                        return 1.0
                    End If
                ElseIf (a < 0.0) Then
                    If (n Mod 2 = 1) Then
                        return -mPow(-a, n)
                    Else
                        return mPow(-a, n)
                    End If
                Else
                    y = 1.0
                    r = a
                    m = n
                    Do While (m > 0)
                        If ((m And 1) = 1) Then
                            y = y * r
                        End If
                        r = r*r
                        m = m >> 1
                    Loop
                    return y
                End If
            End If
        ElseIf (n = 0) Then
            return 1.0
        Else       '  when n < 0
            If (a = 0.0) Then
                Throw New Exception("Negative powering exception of zero.")
            Else
                If (n = -1) Then
                    return 1.0/a
                Else
                    return 1.0/mPow(a, -n)
                End If
            End If
        End If
    End Function

 

    '
    ' Compute the square root of x to a given scale, x > 0.
    '
    Public Function heronSqrt(a As Double) As Double
        Dim x1 As Double
        Dim x2 As Double
        Dim er As Double
        Dim counter As Integer

        If (a < 0.0) Then
            Throw New Exception("Cannot find the Sqrt of a negative number.")
        ElseIf (a = 0.0 Or a = 1.0) Then
            return a
        Else
            x1 = a
            x2 = (x1 + a/x1)/2.0
            er = x1 - x2
            counter = 0
            Do While (x1 + er <> x1)
                x1 = x2
                x2 = (x1 + a/x1)/2.0
                er = x1 - x2
                If (Math.Abs(er) < Math.Abs(M_EPSILON*x1)) Then
                    Exit Do
                End If
                counter = counter + 1
                If (counter > MAX_ITER) Then
                    Exit Do
                End If
            Loop
            If (counter > MAX_ITER) Then
                Throw New Exception("Inaccurate sqrt exception by too many iterations.")
            End If
            return x2
        End If
    End Function

    '
    ' Compute the cubic root of x to a given scale, x > 0.
    '
    Public Function newtonCbrt(a As Double) As Double
        Dim x1 As Double
        Dim x2 As Double
        Dim er As Double
        Dim counter As Integer

        If (a = 0.0 Or a = 1.0 Or a = -1.0) Then
            return a
        ElseIf (a < 0.0) Then
            return -newtonCbrt(-a)
        Else
            x1 = a
            x2 = (2.0*x1 + a/(x1*x1))/3.0
            er = x1 - x2
            counter = 0
            Do While (x1 + er <> x1)
                x1 = x2
                x2 = (2.0*x1 + a/(x1*x1))/3.0
                er = x1 - x2
                If (Math.Abs(er) < Math.Abs(M_EPSILON*x1)) Then
                    Exit Do
                End if
                counter = counter + 1
                if (counter > MAX_ITER) Then
                    Exit Do
                End if
            Loop
            If (counter >= MAX_ITER) Then
                Throw New Exception("Inaccurate cbrt exception by too many iterations.")
            End If
            return x2
        End If
    End Function

    '
    ' Compute the n-th root of x to a given scale, x > 0.
    '
    Public Function newtonNthRoot(n As integer, a As Double) As Double
        Dim x1 As Double
        Dim x2 As Double
        Dim xn As Double
        Dim er As Double
        Dim counter As Integer

        If (n = 0) Then
            return 1.0
        ElseIf (n = 1) Then
            return a
        ElseIf (n > 0) Then
            If (a = 0.0 Or a = 1.0) Then
                return a
            ElseIf (a = -1.0) Then
                If (n Mod 2 = 1) Then
                    return a
                Else
                    Throw New Exception("Cannot find the even n-th root of a negative number.")
                End If
            ElseIf (a < 0.0) Then
                if (n Mod 2 = 1) Then
                    return -newtonNthRoot(n, -a)
                Else
                    Throw New Exception("Cannot find the even n-th root of a negative number.")
                End If
            ElseIf (a < 1.0) Then
                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
                Do While (x1 + er <> x1)
                    x1 = x2
                    xn = mPow(x1, n - 1)
                    x2 = ((n - 1)*x1 + a/xn)/n
                    er = x1 - x2
                    If (Math.Abs(er) < Math.Abs(M_EPSILON*x1)) Then
                        Exit Do
                    End if
                    counter = counter + 1
                    if (counter > MAX_ITER) Then
                        Exit Do
                    End if
                Loop
                if (counter >= MAX_ITER) Then
                    Throw New Exception("Inaccurate n-th root exception by too many iterations.")
                End If
                return x2
            End If
        Else
            If (a = 0.0) Then
                Throw New Exception("Cannot find the negative n-th root of zero.")
            Else
                return 1.0/newtonNthRoot(-n, a)
            End If
        End If
    End Function

    Sub Main()

        Dim x As Double = 16.0
        Dim u As Double = Math.Sqrt(x)

        Console.WriteLine("[ Testing heronSqrt(double) ]--------------------")
        Console.WriteLine("x = " & x )
        Console.WriteLine("u = Sqrt(" & x & ") = " & u )
        Dim y As Double = heronSqrt(x)
        Console.WriteLine("y = heronSqrt(" & x & ") = " & y )
        Console.WriteLine("y*y = " & y*y )
        Console.WriteLine()

        Console.WriteLine("[ Testing newtonCbrt(double) ]--------------------" )
        x = -216.0
        Console.WriteLine("x = " & x )
        Console.WriteLine("-Exp(Log(-x)/3.0) = " & -Math.Exp(Math.Log(-x)/3.0) )
        Dim w As Double = newtonCbrt(x)
        Console.WriteLine("w = newtonCbrt(" & x & ") = " & w )
        Console.WriteLine("w*w*w = " & w*w*w )
        Console.WriteLine()

        x = 729000000000.0
        Console.WriteLine("x = " & x )
        Console.WriteLine("Exp(Log(x)/3.0) = " & Math.Exp(Math.Log(x)/3.0) )
        w = newtonCbrt(x)
        Console.WriteLine("w = newtonCbrt(" & x & ") = " & w )
        Console.WriteLine("w*w*w = " & w*w*w )
        Console.WriteLine()

        Console.WriteLine("[ Testing newtonNthRoot(int, double) ]--------------------" )
        Dim z As Double = newtonNthRoot(3, x)
        Console.WriteLine("x = " & x )
        Console.WriteLine("z = newtonNthRoot(3, " & x & ") = " & z )
        Console.WriteLine("z*z*z = " & z*z*z )
        Console.WriteLine()

        x = 12960000000000000000.0
        z = newtonNthRoot(4, x)
        Console.WriteLine("x = " & x )
        Console.WriteLine("z = newtonNthRoot(4, x) = newtonNthRoot(4, " & x & ") = " & z )
        Console.WriteLine("z*z*z*z = " & z*z*z*z )
        Console.WriteLine()

        x = 1.0/12960000000000000000.0
        z = newtonNthRoot(4, x)
        Console.WriteLine("x = " & x )
        Console.WriteLine("Exp(Log(x)/4.0) = " & Math.Exp(Math.Log(x)/4.0) )
        Console.WriteLine("z = newtonNthRoot(4, x) = newtonNthRoot(4, " & x & ") = " & z )
        Console.WriteLine("z*z*z*z = " & z*z*z*z )
        Console.WriteLine()

        Try
            x = -4.0
            Console.WriteLine("[ Test Exception heronSqrt(double) ]--------------------" )
            Console.WriteLine("x = " & x )
            Console.WriteLine("Calculating heronSqrt(" & x & ")" )
            y = heronSqrt(x)
            Console.WriteLine("y = heronSqrt(" & x & ") = " & y )
            Console.WriteLine("y*y = " & y*y )
            Console.WriteLine()
        Catch ex As Exception
            Console.WriteLine(ex.Message & CRLF & "Caught some exception in calculating heronSqrt(" & x & ")" )
            Console.WriteLine()
        End Try

        Try
            x = -4.0
            Console.WriteLine("[ Test Exception in newtonCbrt(double) ]--------------------" )
            Console.WriteLine("x = " & x )
            Console.WriteLine("Calculating newtonCbrt(" & x & ")" )
            y = newtonCbrt(x)
            Console.WriteLine("y = newtonCbrt(" & x & ") = " & y )
            Console.WriteLine("y*y*y = " & y*y*y )
            Console.WriteLine()
        Catch ex As Exception
            Console.WriteLine(ex.Message & CRLF & "Caught some exception in calculating newtonCbrt(" & x & ")" )
            Console.WriteLine()
        End Try

        Console.WriteLine("[ Test calculations by powering ]-----------------------------" )
        x = 200.0
        z = newtonNthRoot(10, x)
        Console.WriteLine("x = " & x )
        Console.WriteLine("Exp(Log(x)/10.0) = " & Math.Exp(Math.Log(x)/10.0) )
        Console.WriteLine("z = newtonNthRoot(10, x) = newtonNthRoot(10, " & x & ") = " & z )
        Console.WriteLine("Pow(z, 10) = " & Math.Pow(z, 10) )
        Console.WriteLine()

        x = 3001.0
        z = newtonNthRoot(99, x)
        Console.WriteLine("x = " & x )
        Console.WriteLine("Exp(Log(x)/99.0) = " & Math.Exp(Math.Log(x)/99.0) )
        Console.WriteLine("z = newtonNthRoot(99, x) = newtonNthRoot(99, " & x & ") = " & z )
        Console.WriteLine("Pow(z, 99) = " & Math.Pow(z, 99) )
        Console.WriteLine()

        x = 3001.0
        z = newtonNthRoot(-99, x)
        Console.WriteLine("x = " & x )
        Console.WriteLine("Exp(Log(x)/-99.0) = " & Math.Exp(Math.Log(x)/-99.0) )
        Console.WriteLine("z = newtonNthRoot(-99, x) = newtonNthRoot(-99, " & x & ") = " & z )
        Console.WriteLine("1.0/Pow(z, 99) = " & 1.0/Math.Pow(z, 99) )
        Console.WriteLine()

        Console.WriteLine("2.1**2.1 = Pow(2.1, 2.1) = "  & Math.Pow(2.1, 2.1) )
        Console.WriteLine("2.1**(-2.1) = Pow(2.1, -2.1) = "  & Math.Pow(2.1, -2.1) )
        Console.WriteLine("2.1**2.1 * 2.1**(-2.1) = Pow(2.1, 2.1) * Pow(2.1, -2.1) = " & Math.Pow(2.1, 2.1)*Math.Pow(2.1, -2.1) )
        Console.WriteLine("2.1**2.1 = Exp(2.1*Log(2.1)) = "  & Math.Exp(2.1*Math.Log(2.1)) )
        Console.WriteLine("2.1**(-2.1) = Exp(-2.1*Log(2.1)) = " & Math.Exp(-2.1*Math.Log(2.1)) )
        Console.WriteLine("2.1**2.1 * 2.1**(-2.1) = Exp(2.1*Log(2.1)) * Exp(-2.1*Log(2.1)) = "  & Math.Exp(2.1*Math.Log(2.1)) * Math.Exp(-2.1*Math.Log(2.1)) )
        Console.WriteLine()

        Dim k As Integer = 301
        x = -1.029
        Dim t1 As Double = nPow(x, k)
        Dim t2 As Double = gPow(x, k)
        Dim t3 As Double = mPow(x, k)
        Console.WriteLine("t1 = nPow(" & x & ", " & k & ") = " & t1 )
        Console.WriteLine("t2 = gPow(" & x & ", " & k & ") = " & t2 )
        Console.WriteLine("t3 = mPow(" & x & ", " & k & ") = " & t3 )
        Console.WriteLine("t1 / t2 = " & (t1 / t2) )
        Console.WriteLine("t1 - t2 = " & (t1 - t2) )
        Console.Write("t1 == t2 ? ")
        If (t1 = t2) Then Console.WriteLine("yes") Else Console.WriteLine("no")
        Console.WriteLine("t1 / t3 = " & (t1 / t3) )
        Console.WriteLine("t1 - t3 = " & (t1 - t3) )
        If (t1 = t3) Then Console.WriteLine("yes") Else Console.WriteLine("no")
        Console.WriteLine("t2 / t3 = " & (t2 / t3) )
        Console.WriteLine("t2 - t3 = " & (t2 - t3) )
        If (t2 = t3) Then Console.WriteLine("yes") Else Console.WriteLine("no")
        Console.WriteLine()

        Console.WriteLine("Done.")
    End Sub

End Module

 

#If False Then
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 = 729000000000
Exp(Log(x)/3.0) = 9000
w = newtonCbrt(729000000000) = 9000
w*w*w = 729000000000

[ Testing newtonNthRoot(int, double) ]--------------------
x = 729000000000
z = newtonNthRoot(3, 729000000000) = 9000
z*z*z = 729000000000

x = 1.296E+19
z = newtonNthRoot(4, x) = newtonNthRoot(4, 1.296E+19) = 60000
z*z*z*z = 1.296E+19

x = 7.71604938271605E-20
Exp(Log(x)/4.0) = 1.66666666666667E-05
z = newtonNthRoot(4, x) = newtonNthRoot(4, 7.71604938271605E-20) = 1.66666666666
667E-05
z*z*z*z = 7.71604938271605E-20

[ 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.5874010519682
y*y*y = -4

[ Test calculations by powering ]-----------------------------
x = 200
Exp(Log(x)/10.0) = 1.69864646463425
z = newtonNthRoot(10, x) = newtonNthRoot(10, 200) = 1.69864646463425
Pow(z, 10) = 200

x = 3001
Exp(Log(x)/99.0) = 1.08423618932588
z = newtonNthRoot(99, x) = newtonNthRoot(99, 3001) = 1.08423618932588
Pow(z, 99) = 3001

x = 3001
Exp(Log(x)/-99.0) = 0.922308266265993
z = newtonNthRoot(-99, x) = newtonNthRoot(-99, 3001) = 0.922308266265993
1.0/Pow(z, 99) = 3001.00000000001

2.1**2.1 = Pow(2.1, 2.1) = 4.74963809174224
2.1**(-2.1) = Pow(2.1, -2.1) = 0.210542357266885
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.74963809174224
2.1**(-2.1) = Exp(-2.1*Log(2.1)) = 0.210542357266885
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.92801577163
t2 = gPow(-1.029, 301) = -5457.92801577169
t3 = mPow(-1.029, 301) = -5457.92801577169
t1 / t2 = 0.999999999999989
t1 - t2 = 6.18456397205591E-11
t1 == t2 ? no
t1 / t3 = 0.999999999999989
t1 - t3 = 6.18456397205591E-11
no
t2 / t3 = 1
t2 - t3 = 0
yes

Done.
#End If