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

 

Go 언어로 프로그래밍을 할 때 import 구문에 "math" 를 추가하고, 소스에서 math.Pow(밑수, 지수) 함수를 사용하여 지수 계산을 할 수 있다. 하지만 차후 필요한 데가 있을 것 같아서 이와 유사한 n 제곱 함수와 n 제곱근 함수를 구현해 보았다.

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

Go 언어로 작성된 소스를 실행시키고자 할 때, 옵션 run 을 사용하면 소스 파일을 바로 실행(실제로는 임시 폴더에 바이너리 파일을 생상하여 실행)할 수 있고, 옵션 build 를 사용하면 컴파일하여 실행가능 바이너리 파일을 만든다. (그런데 생성된 실행가능 바이너리 파일의 파일 크기가 꽤 크다. 간단히 Hello 한 즐을 출력하는 소스도 컴파일하면 생성된 실행가능 파일의 크기가 1메가를 넘는다.)

Go 언어에서 float64 라는 타입이 C, C++, Java, C# 에서의 double 타입이라고 하는 것과 동일하다.

Go 언어에서 val 에약어로 선언된 것은 상수이고 var 예약어로 선언된 것은 변수이다. (여기서 상수라 함은 저장된 값을 변경하지 못하는 변수라는 뜻으로 immutable 변수로 이해하면 됨.)

        val a = 1.5   // a 는 상수
        a = 2.7       //  에러

        var b = 1.5   // b 는 변수
        b = 2.7        //  b 에는 새로운 값 2,7 이 저장됨

        c := 2.1   // var c = 2.1 와 동일한 의미

Go 언어는 Java, C# 언어에서의 예외상황 처리를 지원하지 않는다. 그래서 예외상황을 처리하는 편법으로 모든 함수(거듭제곱을 구하는 함수와 거듭제곱근을 구하는 함수)는 리턴할 때 항상 계산된 값과 에러 두 개(소스에서 val 과 err 두 개)를 리턴하도록 하였다. 리턴된 에러가 nil 이면 에외상황이 없는 걸로 갖주된다. 이 때문에 함수의 시용이 번거러운 점은 있다. 예를 들어, 아래의 소스 중에 Main() 함수 내에서

       y, err = newtonCbrt(x)

로 된 것은 세제곱근 함수를 호출하여 err 가 nil 이면 리턴 받은 y 의 값이 유효하고 그렇지 않으면 y 의 값이 유효하지 않다는 뜻이다.

 

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

package main

import (
    "fmt"
    "math"
    "errors"
)

const (
    MAX_ITER = 20000
    M_EPSILON = 1.0e-15
)

/**
  * Compute the n-th root of x to a given scale, x > 0.
  */
func nPow(a float64, n int) (result float64, error error) {
    if (n > 0) {
        if (n == 1) {
            return a, nil 
        } else {
            if (a == 0.0 || a == 1.0) {
                return a, nil
            } else if (a == -1.0) {
                if (n % 2 == 1) {
                    return -1.0, nil
                } else {
                    return 1.0, nil
                }
            } else if (a < 0.0) {
                if (n % 2 == 1) {
                    val, err := nPow(-a, n)
                    return -val, err
                } else {
                    return nPow(-a, n)
                }
            } else {
                var y = 1.0
                var i int
                for i = 0; i < n; i++ {
                    y *= a
                }
                return y, nil
            }
        }
    } else if (n == 0) {
        return 1.0, nil
    } else {      //  when n < 0
        if (a == 0.0) {
            return -1.0, errors.New("Negative powering exception of zero.")
        } else {
            if (n == -1) {
                return 1.0/a, nil
            } else {
                val, err := nPow(a, -n)
                return 1.0/val, err
            }
        }
    }
    return -1.0, errors.New("No retured value.")
}

 

/**
  * Compute the n-th root of x to a given scale, x > 0.
  */
func gPow(a float64, n int) (result float64, error error) {
    if (n > 0) {
        if (n == 1) {
            return a, nil
        } else {
            if (a == 0.0 || a == 1.0) {
                return a, nil
            } else if (a == -1.0) {
                if (n % 2 == 1) {
                  return -1.0, nil
                } else {
                    return 1.0, nil
                }
            } else if (a < 0.0) {
                if (n % 2 == 1) {
                    val, err := gPow(-a, n)
                    return -val, err
                } else {
                    val, err := gPow(-a, n)
                    return val, err
                }
            } else {
                var y = 1.0
                var r = a
                var m = 8*4 -1      // 8*sizeof(int) - 1;
                var one = 1
                var i int
                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, nil
            }
        }
    } else if (n == 0) {
        return 1.0, nil
    } else {      //  when n < 0
        if (a == 0.0) {
            return -1.0, errors.New("Negative powering exception of zero.")
        } else {
            if (n == -1) {
                return 1.0/a, nil
            } else {
                val, err := gPow(a, -n)
                return 1.0/val, err
            }
        }
    }
    return -1.0, errors.New("No retured value.")
}

 

/**
  * Compute the n-th root of x to a given scale, x > 0.
  */
func mPow(a float64, n int) (result float64, error error) {
    if (n > 0) {
        if (n == 1) {
            return a, nil
        } else {
            if (a == 0.0 || a == 1.0) {
                return a, nil
            } else if (a == -1.0) {
                if (n % 2 == 1) {
                    return -1.0, nil
                } else {
                    return 1.0, nil
                }
            } else if (a < 0.0) {
                if (n % 2 == 1) {
                    val, err := mPow(-a, n)
                    return -val, err
                } else {
                    val, err := mPow(-a, n)
                    return val, err
                }
            } else {
                var y = 1.0
                var r = a
                var m = n
                for ;m > 0; {
                    if ((m & 0x1) != 0) {
                        y *= r
                    }
                    r = r*r
                    m >>= 1
                }
                return y, nil
            }
        }
    } else if (n == 0) {
        return 1.0, nil
    } else {      //  when n < 0
        if (a == 0.0) {
      return -1.0, errors.New("Negative powering exception of zero.")
        } else {
            if (n == -1) {
                return 1.0/a, nil
            } else {
                val, err := mPow(a, -n)
                return 1.0/val, err
            }
        }
     }
    return -1.0, errors.New("No retured value.")
}

 

/**
  * Compute the square root of x to a given scale, x > 0.
  */
func heronSqrt(a float64) (result float64, error error) {
    if (a < 0.0) {
        return -1.0, errors.New("Cannot find the sqrt of a negative number.")
    } else if (a == 0.0 || a == 1.0) {
        return a, nil
    } else {
        var x1 = a
        var x2 = (x1 + a/x1)/2.0
        var er = x1 - x2
        var counter int = 0
        for ;(x1 + er != x1); {
            x1 = x2
            x2 = (x1 + a/x1)/2.0
            er = x1 - x2
            if (math.Abs(er) < math.Abs(M_EPSILON*x1)) {
                break
            }
            counter++
            if (counter > MAX_ITER) {
                break
            }
        }
        if (counter >= MAX_ITER) {
            return -1.0, errors.New("Inaccurate sqrt exception by too many iterations.")
        }
        return x2, nil
    }
    return -1.0, errors.New("No retured value.")
 }

/**
  * Compute the cubic root of x to a given scale, x > 0.
  */
func newtonCbrt(a float64) (result float64, error error) {
    if (a == 0.0 || a == 1.0 || a == -1.0) {
        return a, nil
    } else if (a < 0.0) {
        val, err := newtonCbrt(-a)
        return -val, err
    } else {
        var x1 = a
        var x2 = (2.0*x1 + a/(x1*x1))/3.0
        var er = x1 - x2;
        var counter int = 0
        for ;(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)) {
                break
            }
            counter++
            if (counter > MAX_ITER) {
                break
            }
        }
        if (counter >= MAX_ITER) {
            return -1.0, errors.New("Inaccurate cbrt exception by too many iterations.")
        }
        return x2, nil
    }
    return -1.0, errors.New("No retured value.")
}

/**
  * Compute the n-th root of x to a given scale, x > 0.
  */
func newtonNthRoot(n int, a float64) (result float64, error error) {
    if (n == 0) {
        return 1.0, nil
    } else if (n == 1) {
        return a, nil
    } else if (n > 0) {
        if (a == 0.0 || a == 1.0) {
      return a, nil
        } else if (a == -1.0) {
            if (n % 2 == 1) {
                return a, nil
            } else {
                return -1.0, errors.New("Cannot find the even n-th root of a negative number.")
            }
        }else if (a < 0.0) {
            if (n % 2 == 1) {
                val, err := newtonNthRoot(n, -a)
                return -val, err
            } else {
                return -1.0, errors.New("Cannot find the even n-th root of a negative number.")
            }
        } else if (a < 1.0) {
            val, err := newtonNthRoot(n, 1.0/a)
            return 1.0/val, err
        } else {
            var x1 = a
            var xn, err = mPow(x1, n - 1)
            if err != nil {
                return -1.0, err
            }
            var x2 = (float64(n - 1) *x1 + a/xn)/float64(n)
            var er = x1 - x2
            var counter int= 0
            for ;(x1 + er != x1); {
                x1 = x2;
                xn, err = mPow(x1, n - 1)
                if err != nil {
                    return -1.0, err
                }
                x2 = (float64(n - 1)*x1 + a/xn)/float64(n)
                er = x1 - x2
                if (math.Abs(er) < math.Abs(M_EPSILON*x1)) {
                    break
                }
                counter++
                if (counter > MAX_ITER) {
                    break
                }
            }
            if (counter >= MAX_ITER) {
                return -1.0, errors.New("Inaccurate n-th root exception by too many iterations.")
            }
            return x2, nil
        }
    } else {
        if (a == 0.0) {
            return -1.0, errors.New("Cannot find the negative n-th root of zero.")
        } else {
            val, err := newtonNthRoot(-n, a)
            return 1.0/val, err
        }
    }
    return -1.0, errors.New("No retured value.")
}

 

func main() {
    x := 16.0
    u := math.Sqrt(x)

    fmt.Printf("[ Testing heronSqrt(double) ]--------------------\n")
    fmt.Printf("x = %g\n", x)
    fmt.Printf("u = sqrt(%g) = %g\n", x, u)

    var y, err = heronSqrt(x)
    if (err == nil) {
        fmt.Printf("y = heronSqrt(%g) = %g\n", x,  y)
        fmt.Printf("y*y = %g\n", y*y)
        fmt.Printf("\n");
    }

    fmt.Printf("[ Testing newtonCbrt(double) ]--------------------\n")
    x = -216.0
    fmt.Printf("x = %g\n", x)
    fmt.Printf("-exp(log(-x)/3.0) = %g\n", -math.Exp(math.Log(-x)/3.0))
    var w float64
    w, err = newtonCbrt(x)
    if (err == nil) {
        fmt.Printf("w = newtonCbrt(%g) = %g\n", x, w)
        fmt.Printf("w*w*w = %g\n", w*w*w);
    }
    fmt.Printf("\n");

    x = 729000000000.0
    fmt.Printf("x = %g\n", x)
    fmt.Printf("exp(log(x)/3.0) = %g\n", math.Exp(math.Log(x)/3.0))
    w, err = newtonCbrt(x)
    if (err == nil) {
        fmt.Printf("w = newtonCbrt(%g) = %g\n", x, w)
        fmt.Printf("w*w*w = %g\n", w*w*w);
    }
    fmt.Printf("\n");

    fmt.Printf("[ Testing newtonNthRoot(int, double) ]--------------------\n");
    var z float64
    z, err = newtonNthRoot(3, x)
    fmt.Printf("x = %g\n", x)
    if (err == nil) {
        fmt.Printf("z = newtonNthRoot(3, %g) = %g\n", x, z)
        fmt.Printf("z*z*z = %g\n", z*z*z);
    }
    fmt.Printf("\n");

    x = 12960000000000000000.0
    z, err = newtonNthRoot(4, x)
    fmt.Printf("x = %g\n", x)
    if (err == nil) {
        fmt.Printf("z = newtonNthRoot(4, %g) = %g\n", x, z)
        fmt.Printf("z*z*z*z = %g\n", z*z*z*z);
    }
    fmt.Printf("\n");

    x = 1.0/12960000000000000000.0
    z, err = newtonNthRoot(4, x)
    fmt.Printf("x = %g\n", x)
    if (err == nil) {
        fmt.Printf("exp(log(x)/4.0) = %g\n", math.Exp(math.Log(x)/4.0))
        fmt.Printf("z = newtonNthRoot(4, x) = newtonNthRoot(4, %g) = %g\n", x, z)
        fmt.Printf("z*z*z*z = %g\n", z*z*z*z);
    }
    fmt.Printf("\n");

    x = -4.0
    fmt.Printf("[ Test Exception heronSqrt(double) ]--------------------\n")
    fmt.Printf("x = %g\n", x)
    fmt.Printf("Calculating heronSqrt(%g)\n", x)
    y, err = heronSqrt(x)
    if (err == nil) {
        fmt.Printf("y = heronSqrt(%g) = %g\n", x, z)
        fmt.Printf("y*y = %g\n", y*y);
    } else {
        fmt.Printf("%s\nCaught some exception in calculating heronSqrt(%g)\n", err, x)
    }
    fmt.Printf("\n");

    x = -4.0
    fmt.Printf("[ Test Exception in newtonCbrt(double) ]--------------------\n")
    fmt.Printf("x = %g\n", x)
    fmt.Printf("Calculating newtonCbrt(%g)\n", x)
    y, err = newtonCbrt(x)
    if (err == nil) {
        fmt.Printf("y = newtonCbrt(%g) = %g\n", x, z)
        fmt.Printf("y*y*y = %g\n", y*y*y);
    } else {
        fmt.Printf("%s\nCaught some exception in calculating newtonCbrt(%g)\n", err, x)
    }
    fmt.Printf("\n");

    fmt.Printf("[ Test calculations by powering ]-----------------------------\n")
    x = 200.0
    z, err = newtonNthRoot(10, x)
    fmt.Printf("x = %g\n", x)
    if (err == nil) {
        fmt.Printf("exp(log(x)/10.0) = %g\n", math.Exp(math.Log(x)/10.0))
        fmt.Printf("z = newtonNthRoot(10, x) = newtonNthRoot(10, %g) = %g\n", x, z)
        fmt.Printf("pow(z, 10) = %g\n", math.Pow(z, 10));
    }
    fmt.Printf("\n");

    x = 3001.0
    z, err = newtonNthRoot(99, x)
    fmt.Printf("x = %g\n", x)
    if (err == nil) {
        fmt.Printf("exp(log(x)/99.0) = %g\n", math.Exp(math.Log(x)/99.0))
        fmt.Printf("z = newtonNthRoot(99, x) = newtonNthRoot(99, %g) = %g\n", x, z)
        fmt.Printf("pow(z, 99) = %g\n", math.Pow(z, 99));
    }
    fmt.Printf("\n");

    x = 3001.0
    z, err = newtonNthRoot(-99, x)
    fmt.Printf("x = %g\n", x)
    if (err == nil) {
        fmt.Printf("exp(log(x)/-99.0) = %g\n", math.Exp(math.Log(x)/-99.0))
        fmt.Printf("z = newtonNthRoot(-99, x) = newtonNthRoot(-99, %g) = %g\n", x, z)
        fmt.Printf("1.0/pow(z, 99) = = %g\n", 1.0/math.Pow(z, 99));
    }
    fmt.Printf("\n");

    fmt.Printf("2.1**2.1 = pow(2.1, 2.1) = %g\n", math.Pow(2.1, 2.1))
    fmt.Printf("2.1**(-2.1) = pow(2.1, -2.1) = %g\n" , math.Pow(2.1, -2.1))
    fmt.Printf("2.1**2.1 * 2.1**(-2.1) = pow(2.1, 2.1) * pow(2.1, -2.1) = %g\n", math.Pow(2.1, 2.1)*math.Pow(2.1, -2.1))
    fmt.Printf("2.1**2.1 = exp(2.1*log(2.1)) = %g\n", math.Exp(2.1*math.Log(2.1)))
    fmt.Printf("2.1**(-2.1) = exp(-2.1*log(2.1)) = %g\n", math.Exp(-2.1*math.Log(2.1)))
    fmt.Printf("2.1**2.1 * 2.1**(-2.1) = exp(2.1*log(2.1)) * exp(-2.1*log(2.1)) = %g\n", math.Exp(2.1*math.Log(2.1)) * math.Exp(-2.1*math.Log(2.1)))
    fmt.Printf("\n");

    fmt.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.000000000000002
w = newtonCbrt(-216) = -6
w*w*w = -216

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

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

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

x = 7.716049382716049e-20
exp(log(x)/4.0) = 1.6666666666666654e-05
z = newtonNthRoot(4, x) = newtonNthRoot(4, 7.716049382716049e-20) = 1.6666666666
666667e-05
z*z*z*z = 7.716049382716051e-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.6666666666666667e-05
y*y*y = -3.999999999999999

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

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

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

2.1**2.1 = pow(2.1, 2.1) = 4.749638091742242
2.1**(-2.1) = pow(2.1, -2.1) = 0.21054235726688478
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.749638091742241
2.1**(-2.1) = exp(-2.1*log(2.1)) = 0.2105423572668848
2.1**2.1 * 2.1**(-2.1) = exp(2.1*log(2.1)) * exp(-2.1*log(2.1)) = 1

Done.
*/