음이 아닌 실수 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 )

 

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

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

아래의 소스는 2.x 버전대의 파이썬에서 실행되도록 작성된 소스이다.

#!/usr/bin/env python
# -*- encoding: utf-8 -*-

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


import math

MAX_ITER = 20000
M_EPSILON = 1.0e-15

#
# Compute the n-th power of x to a given scale, x > 0.
#
def nPow(a, n):
    if n > 0:
        if n == 1:
            return a
        else:
            if a == 0.0 or a == 1.0:
                return a
            elif a == -1.0:
                if n % 2 == 1:
                    return -1.0
                else:
                    return 1.0
            elif a < 0.0:
                if n % 2 == 1:
                    return -nPow(-a, n)
                else:
                    return nPow(-a, n)
            else:
                y = 1.0
                for i in range(n):
                    y *= a
                return y
    elif n == 0:
        return 1.0
    else:      #  when n < 0
        if a == 0.0:
            raise Error,'Negative powering exception of zero.'
        else:
            if n == -1:
                return 1.0/a
            else:
                return 1.0/nPow(a, -n)

 

#
# Compute the n-th power of x to a given scale, x > 0.
#
def gPow(a, n):
    if n > 0:
        if n == 1:
            return a
        else:
            if a == 0.0 or a == 1.0:
                return a
            elif a == -1.0:
                if n % 2 == 1:
                    return -1.0
                else:
                    return 1.0
            elif a < 0.0:
                if n % 2 == 1:
                    return -gPow(-a, n)
                else:
                    return gPow(-a, n)
            else:
                y = 1.0
                r = a
                m = 8*4 - 1            #  8*sizeof(int) - 1;
                one = 1
                for i in range(m + 1):
                    if (n & one) == 0:
                        y *= 1.0
                    else:
                        y *= r
                    r = r*r
                    one <<= 1
                    if one > n:
                        break
                return y
    elif n == 0:
        return 1.0
    else:      #  when n < 0
        if a == 0.0:
            raise Error,'Negative powering exception of zero.'
        else:
            if n == -1:
                return 1.0/a
            else:
                return 1.0/gPow(a, -n)

 

#
# Compute the n-th power of x to a given scale, x > 0.
#
def mPow(a, n):
    if n > 0:
        if n == 1:
            return a
        else:
            if a == 0.0 or a == 1.0:
                return a
            elif a == -1.0:
                if n % 2 == 1:
                    return -1.0
                else:
                    return 1.0
            elif a < 0.0:
                if n % 2 == 1:
                    return -mPow(-a, n)
                else:
                    return mPow(-a, n)
            else:
                y = 1.0
                r = a
                m = n
                while m > 0:
                    if (m & 0x1) == 1:
                        y *= r
                    r = r*r
                    m >>= 1
                return y
    elif n == 0:
        return 1.0
    else:      #  when n < 0
        if a == 0.0:
            raise Error,'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):
    if a < 0.0:
        raise ValueError,'Cannot find the sqrt of a negative number.'
    elif a == 0.0 or 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
            counter += 1
            if counter > MAX_ITER:
                break
        if counter >= MAX_ITER:
            raise ValueError,'Inaccurate sqrt exception by too many iterations.'
        return x2


#
# Compute the cubic root of x to a given scale, x > 0.
#
def newtonCbrt(a):
    if a == 0.0 or a == 1.0 or a == -1.0:
        return a
    elif 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 abs(er) < abs(M_EPSILON*x1):
                break
            counter += 1
            if counter > MAX_ITER:
                break
        if counter >= MAX_ITER:
            raise Error,'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, a):
    if n == 0:
        return 1.0
    elif n == 1:
        return a
    elif n > 0:
        if a == 0.0 or a == 1.0:
            return a
        elif a == -1.0:
            if n % 2 == 1:
                return a
            else:
                raise ValueError,'Cannot find the even n-th root of a negative number.'
        elif a < 0.0:
            if n % 2 == 1:
                return -newtonNthRoot(n, -a)
            else:
                raise ValueError,'Cannot find the even n-th root of a negative number.'
        elif 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
                counter += 1
                if counter > MAX_ITER:
                    break
            if counter >= MAX_ITER:
                raise ValueError, 'Inaccurate n-th root exception by too many iterations.'
            return x2
    else:
        if a == 0.0:
            raise Error, 'Cannot find the negative n-th root of zero.'
        else:
            return 1.0/newtonNthRoot(-n, a)


if __name__ == "__main__":

    x = 16.0
    u = math.sqrt(x)

    print "[ Testing heronSqrt(double) ]--------------------"
    print "x = %g" % x
    print "u = sqrt(%g) = %g" % (x, u)
    y = heronSqrt(x)
    print "y = heronSqrt(%g) = %g" % (x, y)
    print "y*y = %g" % (y*y)
    print

    print "[ Testing newtonCbrt(double) ]--------------------"
    x = -216.0
    print "x = %g" % x
    print "-exp(log(-x)/3.0) = %g" % -math.exp(math.log(-x)/3.0)
    w = newtonCbrt(x)
    print "w = newtonCbrt(%g) = %g" % (x, w)
    print "w*w*w = %g" % (w*w*w)
    print

    x = 729000000000.0
    print "x = %g" % x
    print "exp(log(x)/3.0) = %g" % math.exp(math.log(x)/3.0)
    w = newtonCbrt(x)
    print "w = newtonCbrt(%g) = %g" % (x, w)
    print "w*w*w = %g" % (w*w*w)
    print

    print "[ Testing newtonNthRoot(int, double) ]--------------------"
    z = newtonNthRoot(3, x)
    print "x = %g" % x
    print "z = newtonNthRoot(3, %g) = %g" % (x, z)
    print "z*z*z = %g" % (z*z*z)
    print

    x = 12960000000000000000.0
    z = newtonNthRoot(4, x)
    print "x = %g" % x
    print "z = newtonNthRoot(4, x) = newtonNthRoot(4, %g) =  %g" % (x, z)
    print "z*z*z*z = %g" % (z*z*z*z)
    print

    x = 1.0/12960000000000000000.0
    z = newtonNthRoot(4, x)
    print "x = %g" % x
    print "exp(log(x)/4.0) = %g" % math.exp(math.log(x)/4.0)
    print "z = newtonNthRoot(4, x) = newtonNthRoot(4, %g) =  %g" % (x, z)
    print "z*z*z*z = %g" % (z*z*z*z)
    print


    try:
        x = -4.0
        print "[ Test Exception heronSqrt(double) ]--------------------"
        print "x = %g" % x
        print "Calculating heronSqrt(%g)" % x
        y = heronSqrt(x)
        print "y = heronSqrt(%g) = %g" % (x, y)
        print "y*y = %g" % (y*y)
        print
    except ValueError, ex:
        print "%s\nCaught some exception in calculating heronSqrt(%g)" % (ex, x)
        print


    try:
        x = -4.0
        print "[ Test Exception in newtonCbrt(double) ]--------------------"
        print "x = %g" % x
        print "Calculating newtonCbrt(%g)" % x
        y = newtonCbrt(x)
        print "y = newtonCbrt(%g) = %g" % (x, y)
        print "y*y*y = %g" % (y*y*y)
        print
    except ValueError, ex:
        print "%s\nCaught some exception in calculating newtonCbrt(%g)" % (ex, x)
        print


    print "[ Test calculations by powering ]-----------------------------"
    x = 200.0
    z = newtonNthRoot(10, x)
    print "x = %g" % x
    print "exp(log(x)/10.0) = %g" % math.exp(math.log(x)/10.0)
    print "z = newtonNthRoot(10, x) = newtonNthRoot(10, %g) = %g" % (x, z)
    print "pow(z, 10) = %g" % math.pow(z, 10)
    print

    x = 3001.0
    z = newtonNthRoot(99, x)
    print "x = %g" % x
    print "exp(log(x)/99.0) = %g" % math.exp(math.log(x)/99.0)
    print "z = newtonNthRoot(99, x) = newtonNthRoot(99, %g) = %g" % (x, z)
    print "pow(z, 99) = %g" % math.pow(z, 99)
    print

    x = 3001.0
    z = newtonNthRoot(-99, x)
    print "x = %g" % x
    print "exp(log(x)/-99.0) = %g" % math.exp(math.log(x)/-99.0)
    print "z = newtonNthRoot(-99, x) = newtonNthRoot(-99, %g) = %g" % (x,z)
    print "1.0/pow(z, 99) = %g" % (1.0/math.pow(z, 99))
    print

    print "2.1**2.1 = pow(2.1, 2.1) = %g" % math.pow(2.1, 2.1)
    print "2.1**(-2.1) = pow(2.1, -2.1) = %g" % math.pow(2.1, -2.1)
    print "2.1**2.1 * 2.1**(-2.1) = pow(2.1, 2.1) * pow(2.1, -2.1) = %g" % (math.pow(2.1, 2.1)*math.pow(2.1, -2.1))
    print "2.1**2.1 = exp(2.1*log(2.1)) = %g" % math.exp(2.1*math.log(2.1))
    print "2.1**(-2.1) = exp(-2.1*log(2.1)) = %g" % math.exp(-2.1*math.log(2.1))
    print "2.1**2.1 * 2.1**(-2.1) = exp(2.1*log(2.1)) * exp(-2.1*log(2.1)) = %g" % (math.exp(2.1*math.log(2.1)) * math.exp(-2.1*math.log(2.1)))
    print


    k = 301
    x = -1.029
    t1 = nPow(x, k)
    t2 = gPow(x, k)
    t3 = mPow(x, k)
    print "math.pow(%g, %d) = %g" % (x, k,  math.pow(x, k))
    print "t1 = nPow(%g, %d) = %g" % (x, k,  t1)
    print "t2 = gPow(%g, %g) = %g" % (x, k, t2)
    print "t3 = mPow(%g, %g) = %g" % (x, k, t3)
    print "t1 / t2 = %g" % (t1 / t2)
    print "t1 - t2 = %g" % (t1 - t2)
    print "t1 == t2 ?", 
    if t1 == t2: print "yes"
    else: print "no"
    print "t1 / t3 = %g" % (t1 / t3)
    print "t1 - t3 = %g" % (t1 - t3)
    print "t1 == t3 ?", 
    if t1 == t3: print "yes"
    else: print "no"
    print "t2 / t3 = %g" % (t2 / t3)
    print "t2 - t3 = %g" % (t2 - t3)
    print "t2 == t3 ?", 
    if t2 == t3: print "yes"
    else: print "no"
    print

    print "Done."

"""
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+11
exp(log(x)/3.0) = 9000
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, x) = newtonNthRoot(4, 1.296e+19) =  60000
z*z*z*z = 1.296e+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
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

math.pow(-1.029, 301) = -5457.93
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-11
t1 == t2 ? no
t1 / t3 = 1
t1 - t3 = 6.18456e-11
t1 == t3 ? no
t2 / t3 = 1
t2 - t3 = 0
t2 == t3 ? yes

Done.
"""

 

 

Posted by Scripter
,

음이 아닌 실수 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 )

 

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

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

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

소스 첫 부분에

        private const int MAX_ITER = 20000;
        private const double M_EPSILON = 1.0e-15;

라고 선언하였으니 변수 MAX_ITER 와 M_EPSILON 는 (타입을 갖는) 상수로 선언되었다. const 예약어가 붙으면 static 예약어는 동시에 붙일 수 없다. Java 언어로는 상수를 선언할 방법이 없지만 C# 언어로는 이와 같이 const 예약어(키워드)를 이용하여 상수를 선언할 수 있다.

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


using System;

namespace TestApproximate
{

    public class TestNthRootApp {

        private const int MAX_ITER = 20000;
        private const double M_EPSILON = 1.0e-15;

        /**
         * Compute the n-th root of x to a given scale, x > 0.
         */
        public static double nPow(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 -nPow(-a, n);
                        else
                            return nPow(-a, n);
                    }
                    else {
                        double y = 1.0;
                        for (int i = 0; i < n; i++) {
                            y *= a;
                        }
                        return y;
                    }
                }
            }
            else if (n == 0) {
                return 1.0;
            }
            else {      //  when n < 0
                if (a == 0.0)
                    throw new Exception("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.
         */
        public static double gPow(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 -gPow(-a, n);
                        else
                            return gPow(-a, n);
                    }
                    else {

                        double y = 1.0;
                        double r = a;
                        int m = 8*4 - 1;            ///  8*sizeof(int) - 1;
                        int one = 1;
                        for (int 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)
                    throw new Exception("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.
         */
        public static 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) == 1) {
                                y *= r;
                            }
                            r = r*r;
                            m >>= 1;
                        }
                        return y;
                    }
                }
            }
            else if (n == 0) {
                return 1.0;
            }
            else {      //  when n < 0
                if (a == 0.0)
                    throw new Exception("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.
         */
        public static double heronSqrt(double a) {
            if (a < 0.0) {
                throw new Exception("Cannot find the sqrt of a negative number.");
            }
            else if (a == 0.0 || a == 1.0) {
                return a;
            }
            else {
                double x1 = a;
                double x2 = (x1 + a/x1)/2.0;
                double er = x1 - x2;
                int counter = 0;
                while (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)
                    throw new Exception("Inaccurate sqrt exception by too many iterations.");
                return x2;
            }
        }

        /**
         * Compute the cubic root of x to a given scale, x > 0.
         */
        public static double newtonCbrt(double a) {
            if (a == 0.0 || a == 1.0 || a == -1.0) {
                return a;
            }
            else if (a < 0.0) {
                return -newtonCbrt(-a);
            }
            else {
                double x1 = a;
                double x2 = (2.0*x1 + a/(x1*x1))/3.0;
                double er = x1 - x2;
                int counter = 0;
                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))
                        break;
                    counter++;
                    if (counter > MAX_ITER)
                        break;
                }
                if (counter >= MAX_ITER)
                    throw new Exception("Inaccurate cbrt exception by too many iterations.");
                return x2;
            }
        }

        /**
         * Compute the n-th root of x to a given scale, x > 0.
         */
        public static double newtonNthRoot(int n, double a) {
            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
                        throw new Exception("Cannot find the even n-th root of a negative number.");
                }
                else if (a < 0.0) {
                    if (n % 2 == 1)
                        return -newtonNthRoot(n, -a);
                    else
                        throw new Exception("Cannot find the even n-th root of a negative number.");
                }
                else if (a < 1.0) {
                    return 1.0/newtonNthRoot(n, 1.0/a);
                }
                else {
                    double x1 = a;
                    double xn = mPow(x1, n - 1);
                    double x2 = ((n - 1)*x1 + a/xn)/n;
                    double er = x1 - x2;
                    int counter = 0;
                    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))
                            break;
                        counter++;
                        if (counter > MAX_ITER)
                            break;
                    }
                    if (counter >= MAX_ITER)
                        throw new Exception("Inaccurate n-th root exception by too many iterations.");
                    return x2;
                }
            }
            else {
                if (a == 0.0) {
                    throw new Exception("Cannot find the negative n-th root of zero.");
                }
                else {
                    return 1.0/newtonNthRoot(-n, a);
                }
            }
        }


        public static void Main(string[] args) {

            double x = 16.0;
            double u = Math.Sqrt(x);

            Console.WriteLine("[ Testing heronSqrt(double) ]--------------------");
            Console.WriteLine("x = " + x );
            Console.WriteLine("u = Sqrt(" + x + ") = " + u );
            double y = 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) );
            double w = 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) ]--------------------" );
            double z = 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 (Exception ex) {
                Console.WriteLine(ex.Message + "\n" + "Caught some exception in calculating heronSqrt(" + x + ")");
                Console.WriteLine();
            }


            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 (Exception ex) {
                Console.WriteLine(ex.Message + "\n" + "Caught some exception in calculating newtonCbrt(" + x + ")");
                Console.WriteLine();
            }


            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();


            int k = 301;
            x = -1.029;
            double t1 = nPow(x, k);
            double t2 = gPow(x, k);
            double t3 = 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.WriteLine("t1 == t2 ? " + ((t1 == t2) ? "yes" : "no") );
            Console.WriteLine("t1 / t3 = " + (t1 / t3) );
            Console.WriteLine("t1 - t3 = " + (t1 - t3) );
            Console.WriteLine("t1 == t3 ? " + ((t1 == t3) ? "yes" : "no") );
            Console.WriteLine("t2 / t3 = " + (t2 / t3) );
            Console.WriteLine("t2 - t3 = " + (t2 - t3) );
            Console.WriteLine("t2 == t3 ? " + ((t2 == t3) ? "yes" : "no") );
            Console.WriteLine();

            Console.WriteLine("Done.");
        }
    }
}

/*
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
t1 == t3 ? no
t2 / t3 = 1
t2 - t3 = 0
t2 == t3 ? yes

Done.
*/

 

 

 

Posted by Scripter
,

음이 아닌 실수 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 )

 

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

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

아래의 소스는 대부분 버전의 JVM(자바가상기계) 위에서 컴파일 되고 실행되게 작성된 소스이다.

소스 첫 부분에

    final private static int MAX_ITER = 20000;
    final private static double M_EPSILON = 1.0e-15;

라고 선언하였으니 변수 MAX_ITER 와 M_EPSILON 는 상수는 아니지만 상수와 거의 같은 효과를 같는 클래스 소속 변수(스태틱 변수)이다. Java 언어에는 C 언어나 C++ 언어에서 말하는 상수 선언이 없다.

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


public class TestNthRootApp {

    final private static int MAX_ITER = 20000;
    final private static double M_EPSILON = 1.0e-15;

    /**
     * Compute the n-th root of x to a given scale, x > 0.
     */
    public static double nPow(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 -nPow(-a, n);
                    else
                        return nPow(-a, n);
                }
                else {
                    double y = 1.0;
                    for (int i = 0; i < n; i++) {
                        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.
     */
    public static double gPow(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 -gPow(-a, n);
                    else
                        return gPow(-a, n);
                }
                else {

                    double y = 1.0;
                    double r = a;
                    int m = 8*4 - 1;            ///  8*sizeof(int) - 1;
                    int one = 1;
                    for (int 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)
                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.
     */
    public static 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) == 1) {
                            y *= r;
                        }
                        r = r*r;
                        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.
     */
    public static double heronSqrt(double a) {
        if (a < 0.0) {
            throw new RuntimeException("Cannot find the sqrt of a negative number.");
        }
        else if (a == 0.0 || a == 1.0) {
            return a;
        }
        else {
            double x1 = a;
            double x2 = (x1 + a/x1)/2.0;
            double er = x1 - x2;
            int counter = 0;
            while (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)
                throw new RuntimeException("Inaccurate sqrt exception by too many iterations.");
            return x2;
        }
    }

    /**
     * Compute the cubic root of x to a given scale, x > 0.
     */
    public static double newtonCbrt(double a) {
        if (a == 0.0 || a == 1.0 || a == -1.0) {
            return a;
        }
        else if (a < 0.0) {
            return -newtonCbrt(-a);
        }
        else {
            double x1 = a;
            double x2 = (2.0*x1 + a/(x1*x1))/3.0;
            double er = x1 - x2;
            int counter = 0;
            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))
                    break;
                counter++;
                if (counter > MAX_ITER)
                    break;
            }
            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.
     */
    public static double newtonNthRoot(int n, double a) {
        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
                    throw new RuntimeException("Cannot find the even n-th root of a negative number.");
            }
            else if (a < 0.0) {
                if (n % 2 == 1)
                    return -newtonNthRoot(n, -a);
                else
                    throw new RuntimeException("Cannot find the even n-th root of a negative number.");
            }
            else if (a < 1.0) {
                return 1.0/newtonNthRoot(n, 1.0/a);
            }
            else {
                double x1 = a;
                double xn = mPow(x1, n - 1);
                double x2 = ((n - 1)*x1 + a/xn)/n;
                double er = x1 - x2;
                int counter = 0;
                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))
                        break;
                    counter++;
                    if (counter > MAX_ITER)
                        break;
                }
                if (counter >= MAX_ITER)
                    throw new RuntimeException("Inaccurate n-th root exception by too many iterations.");
                return x2;
            }
        }
        else {
            if (a == 0.0) {
                throw new RuntimeException("Cannot find the negative n-th root of zero.");
            }
            else {
                return 1.0/newtonNthRoot(-n, a);
            }
        }
    }


    public static void main(String[] args) {

        double x = 16.0;
        double u = Math.sqrt(x);

        System.out.println("[ Testing heronSqrt(double) ]--------------------");
        System.out.println("x = " + x );
        System.out.println("u = sqrt(" + x + ") = " + u );
        double y = heronSqrt(x);
        System.out.println("y = heronSqrt(" + x + ") = " + y );
        System.out.println("y*y = " + y*y );
        System.out.println();

        System.out.println("[ Testing newtonCbrt(double) ]--------------------" );
        x = -216.0;
        System.out.println("x = " + x );
        System.out.println("-exp(log(-x)/3.0) = " + -Math.exp(Math.log(-x)/3.0) );
        double w = newtonCbrt(x);
        System.out.println("w = newtonCbrt(" + x + ") = " + w );
        System.out.println("w*w*w = " + w*w*w );
        System.out.println();

        x = 729000000000.0;
        System.out.println("x = " + x );
        System.out.println("exp(log(x)/3.0) = " + Math.exp(Math.log(x)/3.0) );
        w = newtonCbrt(x);
        System.out.println("w = newtonCbrt(" + x + ") = " + w );
        System.out.println("w*w*w = " + w*w*w );
        System.out.println();

        System.out.println("[ Testing newtonNthRoot(int, double) ]--------------------" );
        double z = newtonNthRoot(3, x);
        System.out.println("x = " + x );
        System.out.println("z = newtonNthRoot(3, " + x + ") = " + z );
        System.out.println("z*z*z = " + z*z*z );
        System.out.println();

        x = 12960000000000000000.0;
        z = newtonNthRoot(4, x);
        System.out.println("x = " + x );
        System.out.println("z = newtonNthRoot(4, x) = newtonNthRoot(4, " + x + ") = " + z );
        System.out.println("z*z*z*z = " + z*z*z*z );
        System.out.println();

        x = 1.0/12960000000000000000.0;
        z = newtonNthRoot(4, x);
        System.out.println("x = " + x );
        System.out.println("exp(log(x)/4.0) = " + Math.exp(Math.log(x)/4.0) );
        System.out.println("z = newtonNthRoot(4, x) = newtonNthRoot(4, " + x + ") = " + z );
        System.out.println("z*z*z*z = " + z*z*z*z );
        System.out.println();


        try {
            x = -4.0;
            System.out.println("[ Test Exception heronSqrt(double) ]--------------------" );
            System.out.println("x = " + x );
            System.out.println("Calculating heronSqrt(" + x + ")" );
            y = heronSqrt(x);
            System.out.println("y = heronSqrt(" + x + ") = " + y );
            System.out.println("y*y = " + y*y );
            System.out.println();
        }
        catch (Exception ex) {
            System.out.println(ex.getMessage() + "\n" + "Caught some exception in calculating heronSqrt(" + x + ")" );
            System.out.println();
        }


        try {
            x = -4.0;
            System.out.println("[ Test Exception in newtonCbrt(double) ]--------------------" );
            System.out.println("x = " + x );
            System.out.println("Calculating newtonCbrt(" + x + ")" );
             = newtonCbrt(x);
            System.out.println("y = newtonCbrt(" + x + ") = " + y );
            System.out.println("y*y*y = " + y*y*y );
            System.out.println();
        }
        catch (Exception ex) {
            System.out.println(ex.getMessage() + "\n" + "Caught some exception in calculating newtonCbrt(" + x + ")");
            System.out.println();
        }


        System.out.println("[ Test calculations by powering ]-----------------------------" );
        x = 200.0;
        z = newtonNthRoot(10, x);
        System.out.println("x = " + x );
        System.out.println("exp(log(x)/10.0) = " + Math.exp(Math.log(x)/10.0) );
        System.out.println("z = newtonNthRoot(10, x) = newtonNthRoot(10, " + x + ") = " + z );
        System.out.println("pow(z, 10) = " + Math.pow(z, 10) );
        System.out.println();

        x = 3001.0;
        z = newtonNthRoot(99, x);
        System.out.println("x = " + x );
        System.out.println("exp(log(x)/99.0) = " + Math.exp(Math.log(x)/99.0) );
        System.out.println("z = newtonNthRoot(99, x) = newtonNthRoot(99, " + x + ") = " + z );
        System.out.println("pow(z, 99) = " + Math.pow(z, 99) );
        System.out.println();

        x = 3001.0;
        z = newtonNthRoot(-99, x);
        System.out.println("x = " + x );
        System.out.println("exp(log(x)/-99.0) = " + Math.exp(Math.log(x)/-99.0) );
        System.out.println("z = newtonNthRoot(-99, x) = newtonNthRoot(-99, " + x + ") = " + z );
        System.out.println("1.0/pow(z, 99) = " + 1.0/Math.pow(z, 99) );
        System.out.println();


        System.out.println("2.1**2.1 = pow(2.1, 2.1) = "  + Math.pow(2.1, 2.1) );
        System.out.println("2.1**(-2.1) = pow(2.1, -2.1) = "  + Math.pow(2.1, -2.1) );
        System.out.println("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) );
        System.out.println("2.1**2.1 = exp(2.1*log(2.1)) = "  + Math.exp(2.1*Math.log(2.1)) );
        System.out.println("2.1**(-2.1) = exp(-2.1*log(2.1)) = " + Math.exp(-2.1*Math.log(2.1)) );
        System.out.println("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)) );
        System.out.println();


        int k = 301;
        x = -1.029;
        double t1 = nPow(x, k);
        double t2 = gPow(x, k);
        double t3 = mPow(x, k);
        System.out.println("t1 = nPow(" + x + ", " + k + ") = " + t1 );
        System.out.println("t2 = gPow(" + x + ", " + k + ") = " + t2 );
        System.out.println("t3 = mPow(" + x + ", " + k + ") = " + t3 );
        System.out.println("t1 / t2 = " + (t1 / t2) );
        System.out.println("t1 - t2 = " + (t1 - t2) );
        System.out.println("t1 == t2 ? " + ((t1 == t2) ? "yes" : "no") );
        System.out.println("t1 / t3 = " + (t1 / t3) );
        System.out.println("t1 - t3 = " + (t1 - t3) );
        System.out.println("t1 == t3 ? " + ((t1 == t3) ? "yes" : "no") );
        System.out.println("t2 / t3 = " + (t2 / t3) );
        System.out.println("t2 - t3 = " + (t2 - t3) );
        System.out.println("t2 == t3 ? " + ((t2 == t3) ? "yes" : "no") );
        System.out.println();

        System.out.println("Done.");
    }
}

/*
[ Testing heronSqrt(double) ]--------------------
x = 16.0
u = sqrt(16.0) = 4.0
y = heronSqrt(16.0) = 4.0
y*y = 16.0

[ Testing newtonCbrt(double) ]--------------------
x = -216.0
-exp(log(-x)/3.0) = -6.000000000000001
w = newtonCbrt(-216.0) = -6.0
w*w*w = -216.0

x = 7.29E11
exp(log(x)/3.0) = 9000.000000000004
w = newtonCbrt(7.29E11) = 9000.0
w*w*w = 7.29E11

[ Testing newtonNthRoot(int, double) ]--------------------
x = 7.29E11
z = newtonNthRoot(3, 7.29E11) = 9000.0
z*z*z = 7.29E11

x = 1.296E19
z = newtonNthRoot(4, x) = newtonNthRoot(4, 1.296E19) = 60000.0
z*z*z*z = 1.296E19

x = 7.716049382716049E-20
exp(log(x)/4.0) = 1.666666666666666E-5
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
pow(z, 10) = 199.9999999999999

x = 3001.0
exp(log(x)/99.0) = 1.0842361893258805
z = newtonNthRoot(99, x) = newtonNthRoot(99, 3001.0) = 1.0842361893258805
pow(z, 99) = 3000.9999999999955

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

2.1**2.1 = pow(2.1, 2.1) = 4.749638091742242
2.1**(-2.1) = pow(2.1, -2.1) = 0.21054235726688475
2.1**2.1 * 2.1**(-2.1) = pow(2.1, 2.1) * pow(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)) = 1.0

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.
*/


 

 

Posted by Scripter
,

음이 아닌 실수 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 )

 

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

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

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

예외상황 처리 때문에, Visual C++ 로 컴파일 할 때는 /EHsc 옵션을 붙여 주어야 하고, g++ 로 컴파일할 때는 -fexceptions 옵션을 반드시 붙여 주어야 한다.

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

#include <stdio.h>
#include <iostream>
#include <cmath>
#include <cstring>
// #include <memory.h>
#include <exception>

using namespace std;

class CalcException: public exception
{
  virtual const char* what() const throw()
  {
    return "Not supported calculation exception";
  }
} cannotCalcEx;

class NotDefinedException: public exception
{
  virtual const char* what() const throw()
  {
    return "Not defined calculation exception";
  }
} notDefinedEx;

class MayInaccurateException: public exception
{
  virtual const char* what() const throw()
  {
    return "Inaccurate approximation by too many iterrations.";
  }
} inaccurateApproxEx;

 

#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)
{
    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 {
                double y = 1.0;
                for (int i = 0; i < n; i++) {
                    y *= a;
                }
                return y;
            }
        }
    }
    else if (n == 0) {
        return 1.0;
    }
    else {      //  when n < 0
        if (a == 0.0)
            // throw "Negative powering exception of zero.";
            throw cannotCalcEx;
        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)
{
    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 {


                double y = 1.0;
                double r = a;
                int m = 8*sizeof(int) - 1;   // ignore the most significant bit which means the +/- sign.
                int one = 1;
                for (int 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)
            // throw "Negative powering exception of zero.";
            throw cannotCalcEx;
        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)
            // throw "Negative powering exception of zero.";
            throw cannotCalcEx;
        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)
{
    if (a < 0.0) {
        // throw "Cannot find the sqrt of a negative number.";
        throw notDefinedEx;
    }
    else if (a == 0.0 || a == 1.0) {
        return a;
    }
    else {
        double x1 = a;
        double x2 = (x1 + a/x1)/2.0;
        double er = x1 - x2;
        int 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;
            counter++;
            if (counter > MAX_ITER)
                break;
        }
        if (counter >= MAX_ITER)
            // throw "Inaccurate sqrt exception by too many iterations.";
            throw inaccurateApproxEx;
        return x2;
    }
}

/**
  * Compute the cubic root of x to a given scale, x > 0.
  */
double newtonCbrt(double a)
{
    if (a == 0.0 || a == 1.0 || a == -1.0) {
        return a;
    }
    else if (a < 0.0) {
        return -newtonCbrt(-a);
    }
    else {
        double x1 = a;
        double x2 = (2.0*x1 + a/(x1*x1))/3.0;
        double er = x1 - x2;
        int counter = 0;
        while (x1 + er != x1) {
            x1 = x2;
            x2 = (2.0*x1 + a/(x1*x1))/3.0;
            er = x1 - x2;
            if (abs(er) < abs(M_EPSILON*x1))
                break;
            counter++;
            if (counter > MAX_ITER)
                break;
        }
        if (counter >= MAX_ITER)
            // throw "Inaccurate cbrt exception by too many iterations.";
            throw inaccurateApproxEx;
        return x2;
    }
}

/**
  * Compute the n-th root of x to a given scale, x > 0.
  */
double newtonNthRoot(int n, double a)
{
    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
                // throw "Cannot find the even n-th root of a negative number.";
                throw notDefinedEx;
        }
        else if (a < 0.0) {
            if (n % 2 == 1)
                return -newtonNthRoot(n, -a);
            else
                // throw "Cannot find the even n-th root of a negative number.";
                throw notDefinedEx;
        }
        else if (a < 1.0) {
            return 1.0/newtonNthRoot(n, 1.0/a);
        }
        else {
            double x1 = a;
            double xn = mPow(x1, n - 1);
            double x2 = ((n - 1)*x1 + a/xn)/n;
            double er = x1 - x2;
            int 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;
                counter++;
                if (counter > MAX_ITER)
                    break;
            }
            if (counter >= MAX_ITER)
                // throw "Inaccurate n-th root exception by too many iterations.";
                throw inaccurateApproxEx;
            return x2;
        }
    }
    else {
        if (a == 0.0) {
            // throw "Cannot find the negative n-th root of zero.";
            throw cannotCalcEx;
        }
        else {
            return 1.0/newtonNthRoot(-n, a);
        }
    }
}

 

int main()
{
    std:cout.precision(16);

    double x = 16.0;
    double u = sqrt(x);


    cout << "[ Testing heronSqrt(double) ]--------------------" << endl;
    cout << "x = " << x << endl;
    cout << "u = sqrt(" << x << ") = " << u << endl;
    double y = heronSqrt(x);
    cout << "y = heronSqrt(" << x << ") = " << y << endl;
    cout << "y*y = " << y*y << endl;
    cout << endl;

    cout << "[ Testing newtonCbrt(double) ]--------------------" << endl;
    x = -216.0;
    cout << "x = " << x << endl;
    cout << "-exp(log(-x)/3.0) = " << -exp(log(-x)/3.0) << endl;
    double w = newtonCbrt(x);
    cout << "w = newtonCbrt(" << x << ") = " << w << endl;
    cout << "w*w*w = " << w*w*w << endl;
    cout << endl;

    x = 729000000000.0;
    cout << "x = " << x << endl;
    cout << "exp(log(x)/3.0) = " << exp(log(x)/3.0) << endl;
    w = newtonCbrt(x);
    cout << "w = newtonCbrt(" << x << ") = " << w << endl;
    cout << "w*w*w = " << w*w*w << endl;
    cout << endl;

    cout << "[ Testing newtonNthRoot(int, double) ]--------------------" << endl;
    double z = newtonNthRoot(3, x);
    cout << "x = " << x << endl;
    cout << "z = newtonNthRoot(3, " << x << ") = " << z << endl;
    cout << "z*z*z = " << z*z*z << endl;
    cout << endl;

    x = 12960000000000000000.0;
    z = newtonNthRoot(4, x);
    cout << "x = " << x << endl;
    cout << "z = newtonNthRoot(4, x) = newtonNthRoot(4, " << x << ") = " << z << endl;
    cout << "z*z*z*z = " << z*z*z*z << endl;
    cout << endl;

    x = 1.0/12960000000000000000.0;
    z = newtonNthRoot(4, x);
    cout << "x = " << x << endl;
    cout << "exp(log(x)/4.0) = " << exp(log(x)/4.0) << endl;
    cout << "z = newtonNthRoot(4, x) = newtonNthRoot(4, " << x << ") = " << z << endl;
    cout << "z*z*z*z = " << z*z*z*z << endl;
    cout << endl;


    try {
        x = -4.0;
        cout << "[ Test Exception heronSqrt(double) ]--------------------" << endl;
        cout << "x = " << x << endl;
        cout << "Calculating heronSqrt(" << x << ")" << endl;
        y = heronSqrt(x);
        cout << "y = heronSqrt(" << x << ") = " << y << endl;
        cout << "y*y = " << y*y << endl;
        cout << endl;
    }
    // catch( char *str )    {
    catch( exception& ex )    {
        cout << ex.what() << endl << "Caught some exception in calculating heronSqrt(" << x << ")" << endl;
        cout << endl;
    }

    try {
        x = -4.0;
        cout << "[ Test Exception in newtonCbrt(double) ]--------------------" << endl;
        cout << "x = " << x << endl;
        cout << "Calculating newtonCbrt(" << x << ")" << endl;
        y = newtonCbrt(x);
        cout << "y = newtonCbrt(" << x << ") = " << y << endl;
        cout << "y*y*y = " << y*y*y << endl;
        cout << endl;
    }
    // catch( char *str )    {
    catch( exception& ex )    {
        cout << ex.what() << endl << "Caught some exception in calculating newtonCbrtrt(" << x <<  ")" << endl;
        cout << endl;
    }

    cout << "[ Test calculations by powering ]-----------------------------" << endl;
    x = 200.0;
    z = newtonNthRoot(10, x);
    cout << "x = " << x << endl;
    cout << "exp(log(x)/10.0) = " << exp(log(x)/10.0) << endl;
    cout << "z = newtonNthRoot(10, x) = newtonNthRoot(10, " << x << ") = " << z << endl;
    cout << "pow(z, 10) = " << pow(z, 10) << endl;
    cout << endl;

    x = 3001.0;
    z = newtonNthRoot(99, x);
    cout << "x = " << x << endl;
    cout << "exp(log(x)/99.0) = " << exp(log(x)/99.0) << endl;
    cout << "z = newtonNthRoot(99, x) = newtonNthRoot(99, " << x << ") = " << z << endl;
    cout << "pow(z, 99) = " << pow(z, 99) << endl;
    cout << endl;

    x = 3001.0;
    z = newtonNthRoot(-99, x);
    cout << "x = " << x << endl;
    cout << "exp(log(x)/-99.0) = " << exp(log(x)/-99.0) << endl;
    cout << "z = newtonNthRoot(-99, x) = newtonNthRoot(-99, " << x << ") = " << z << endl;
    cout << "1.0/pow(z, 99) = " << 1.0/pow(z, 99) << endl;
    cout << endl;

    cout << "2.1**2.1 = pow(2.1, 2.1) = "  << pow(2.1, 2.1) << endl;
    cout << "2.1**(-2.1) = pow(2.1, -2.1) = "  << pow(2.1, -2.1) << endl;
    cout << "2.1**2.1 * 2.1**(-2.1) = pow(2.1, 2.1) * pow(2.1, -2.1) = "  << pow(2.1, 2.1)*pow(2.1, -2.1) << endl;
    cout << "2.1**2.1 = exp(2.1*log(2.1)) = "  << exp(2.1*log(2.1)) << endl;
    cout << "2.1**(-2.1) = exp(-2.1*log(2.1)) = "  << exp(-2.1*log(2.1)) << endl;
    cout << "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)) << endl;
    cout << endl;


    int k = 301;
    x = -1.029;
    double t1 = nPow(x, k);
    double t2 = gPow(x, k);
    double t3 = mPow(x, k);
    cout << "t1 = nPow(" << x << ", " << k << ") = " << t1 << endl;
    cout << "t2 = gPow(" << x << ", " << k << ") = " << t2 << endl;
    cout << "t3 = mPow(" << x << ", " << k << ") = " << t3 << endl;
    cout << "t1 / t2 = " << (t1 / t2) << endl;
    cout << "t1 - t2 = " << (t1 - t2) << endl;
    cout << "t1 == t2 ? " << ((t1 == t2) ? "yes" : "no") << endl;
    cout << "t1 / t3 = " << (t1 / t3) << endl;
    cout << "t1 - t3 = " << (t1 - t3) << endl;
    cout << "t1 == t3 ? " << ((t1 == t3) ? "yes" : "no") << endl;
    cout << "t2 / t3 = " << (t2 / t3) << endl;
    cout << "t2 - t3 = " << (t2 - t3) << endl;
    cout << "t2 == t3 ? " << ((t2 == t3) ? "yes" : "no") << endl;
    cout << endl;

    cout << "Done." << endl;
}

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

x = 729000000000
exp(log(x)/3.0) = 9000.000000000004
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+019
z = newtonNthRoot(4, x) = newtonNthRoot(4, 1.296e+019) = 60000
z*z*z*z = 1.296e+019

x = 7.716049382716049e-020
exp(log(x)/4.0) = 1.666666666666666e-005
z = newtonNthRoot(4, x) = newtonNthRoot(4, 7.716049382716049e-020) = 1.666666666
666667e-005
z*z*z*z = 7.716049382716052e-020

[ Test Exception heronSqrt(double) ]--------------------
x = -4
Calculating heronSqrt(-4)
Not defined calculation exception
Caught some exception in calculating heronSqrt(-4)

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

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

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

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

2.1**2.1 = pow(2.1, 2.1) = 4.749638091742242
2.1**(-2.1) = pow(2.1, -2.1) = 0.2105423572668848
2.1**2.1 * 2.1**(-2.1) = pow(2.1, 2.1) * pow(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.2105423572668848
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.928015771692
t3 = mPow(-1.029, 301) = -5457.928015771692
t1 / t2 = 0.9999999999999887
t1 - t2 = 6.184563972055912e-011
t1 == t2 ? no
t1 / t3 = 0.9999999999999887
t1 - t3 = 6.184563972055912e-011
t1 == t3 ? no
t2 / t3 = 1
t2 - t3 = 0
t2 == t3 ? yes

Done.
*/

 

 

Posted by Scripter
,