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

 

OCaml 언어에는 부동소수점수의 지수연산자 ** 가 있다. 하지만 차후 필요한 데가 있을 것 같아서 이와 유사한 n 제곱 함수와 n 제곱근 함수를 구현해 보았다.

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

 

1. OCaml 언어는 (절차적 프로그래밍도 지원하는) 함수형 언어이다.

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

2. OCaml 언어는 (C 언어 처럼 소스에 들여쓰기 규칙을 지키지 않는다.

3. OCaml 언어에서는 상수든 변수든 함수든 처음 선언할 때는 구문 선두에 let 예약어를 반드시 붙인다.

     let a = 2.7 in               // a (통용범위가 한정된) 변수
     let b = ref 5.3 in          // b 는 (통용범위가 한정된) 참조 변수 
     let g x = x*x  in           // g 는 (통용범위가 한정된) 함수

아래의 소스 첫 부분에

        let nMAX_ITER = 20000 ;;
        let nM_EPSILON = 1.0e-15 ll

라고 선언하였으니, nMAX_ITER 와 nM_EPSILON 는 (재정의 되지 않는 한) 상수로 사용된다. 상수나 변수명은 대문자로 시작되지 않아야 한다.

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

5. 예외상황 처리를 위해 예외 던지기 구문 raise ... 과 예외 받기 구문 try ... with ...을 이용하였다. (해당 부분 주석을 제거하면 확인 가능)

6. OCaml 언어에서 부동소수점수의 부호를 반대로 바꾸는 단항 연산자 기호는 - 가 아니고 -. 이다.

 

(*
 * Filename: testNthRoot.ml
 *
 *            Approximate square roots, cubic roots and n-th roots of a given number.
 *
 * Execute: ocaml testNthRoot.ml
 *  Or
 * Compile: ocamlc -o testNthRoot.exe testNthRoot.ml
 * Execute: testNthRoot
 *
 * Date: 2013. 2. 6.
 * Copyright (c) 2013  pkim __AT__ scripts.pe.kr
 *)

open Printf ;;

exception InvalidPower ;;
exception SqrtOfNegative ;;
exception EventhRootOfNegative ;;
exception TooManyInterations ;;

let nMAX_ITER = 20000 ;;
let nM_EPSILON = 1.0e-15 ;;


(*
   Compute the n-th root of x to a given scale, x > 0.
*)
let rec nPow a n =
    if n > 0 then begin
        if n = 1 then
            a
        else begin
            if a = 0.0 or a = 1.0 then
                a
            else if a = -1.0 then begin
                if n mod 2 = 1 then
                    -1.0
                else
                    1.0
            end
            else if a < 0.0 then
                if n mod 2 = 1 then
                    -. (nPow (-. a) n)
                else
                    nPow (-. a) n
            (*
            *)
            else begin
                let y = ref 1.0 in
                for i = 1 to n do
                    y := !y *. a
                done ;
                !y
            end
        end
    end
    else if n = 0 then begin
        1.0
    end
    else begin     (*  when n < 0  *)
        if a = 0.0 then
            raise InvalidPower
        else begin
            if n = -1 then
                1.0 /. a
            else
               1.0 /. (nPow a (-n))
        end
    end ;;
(*
    end
    5.0 ;;
*)

(*
  Compute the n-th root of x to a given scale, x > 0.
*)
let rec gPow a n =
    if n > 0 then
        if n = 1 then
            a
        else
            if a = 0.0 || a = 1.0 then
                a
            else if a = -1.0 then
                if n mod 2 = 1 then
                    -1.0
                else
                    1.0
            else if a < 0.0 then
                if n mod 2 = 1 then
                    -. (gPow (-. a) n)
                else
                    gPow (-.a) n
            else
                let y = ref 1.0 in
                let r = ref (0.0 +. a) in
                let m = 8*4 - 1 in           (*  8*sizeof(int) - 1; *)
                let one = ref 1 in
                for i = 0 to m do
                    if (n land !one) = 0 then
                        y := !y *. 1.0
                    else
                        y := !y *. !r ;
                    r := !r *. !r ;
                    one := (!one lsl 1) ;    (*  one = one << 1  *)
                done ;
                !y
    else if n = 0 then
        1.0
    else      (*  when n < 0  *)
        if a = 0.0 then
            raise InvalidPower
        else
            if n = -1 then
                1.0 /. a
            else
                1.0 /. (gPow a (-n))  ;;

 

(*
   Compute the n-th root of x to a given scale, x > 0.
*)
let rec mPow a n =
    if n > 0 then
        if n = 1 then
            a
        else
            if a = 0.0 || a = 1.0 then
                a
            else if a = -1.0 then
                if n mod 2 = 1 then
                    -1.0
                else
                    1.0
            else if a < 0.0 then
                if n mod 2 = 1 then
                    -. (mPow (-. a) n)
                else
                    mPow (-. a) n
            else
                let y = ref 1.0 in
                let r = ref (0.0 +. a) in
                let m = ref (0 + n)  in
                while !m > 0 do
                    if (!m land 0x1) = 1 then
                        y := !y *. !r ;
                    r := !r *. !r ;
                    m := (!m lsr 1) ;   (* m <- (m >>> 1) *)
                done ;
                !y
    else if n = 0 then
        1.0
    else      (*  when n < 0  *)
        if a = 0.0 then
            raise InvalidPower
        else
            if n = -1 then
                1.0 /. a
            else
                1.0 /. (mPow a (-n)) ;;


(*
   Compute the square root of x to a given scale, x > 0.
*)
let rec heronSqrt a =
    if a < 0.0 then
        raise SqrtOfNegative
    else if a = 0.0 || a = 1.0 then
        a
    else
        let x1 = ref (0.0 +. a) in
        let x2 = ref ((!x1 +. a /. !x1) /. 2.0) in
        let er = ref (!x1 -. !x2) in
        let counter = ref 0 in
        let not_stop = ref true in
        while (((!x1 +. !er) != !x1) && !not_stop) do
            x1 := !x2 ;
            x2 := (!x1 +. a /. !x1) /. 2.0 ;
            er := !x1 -. !x2 ;
            if (abs_float !er) < (abs_float (nM_EPSILON *. !x1)) then
                not_stop := false ;
            counter := !counter + 1 ;
            if !counter > nMAX_ITER then
                not_stop := false ;
        done ;
        if !counter >= nMAX_ITER then
            raise TooManyInterations ;
        !x2 ;;


(*
  Compute the cubic root of x to a given scale, x > 0.
*)
let rec newtonCbrt a =
    if a = 0.0 || a = 1.0 || a = -1.0 then
        a
    else if a < 0.0 then
        -. newtonCbrt (-. a)
    else
        let x1 = ref (0.0 +. a) in
        let x2 = ref ((2.0 *. !x1 +. a /. (!x1 *. !x1)) /. 3.0) in
        let er = ref (!x1 -. !x2) in
        let counter = ref 0 in
        let not_stop = ref true in
        while (!x1 +. !er <> !x1) && !not_stop do
            x1 := !x2 ;
            x2 := (2.0 *. !x1 +. a /. (!x1 *. !x1)) /. 3.0 ;
            er := !x1 -. !x2 ;
            if (abs_float !er) < (abs_float (nM_EPSILON *. !x1)) then
                not_stop := false ;
            counter := !counter + 1 ;
            if !counter > nMAX_ITER then
                not_stop := false
        done ;
        if !counter >= nMAX_ITER then
            raise TooManyInterations ;
        !x2 ;;


(*
  Compute the n-th root of x to a given scale, x > 0.
*)
let rec newtonNthRoot n a =
    if n = 0 then
        1.0
    else if n = 1 then
        a
    else if n > 0 then
        if a = 0.0 || a = 1.0 then
            a
        else if a = -1.0 then
            if n mod 2 = 1 then
                a
            else
                raise EventhRootOfNegative
        else if a < 0.0 then
            if n mod 2 = 1 then
                -. (newtonNthRoot n (-. a))
            else
                raise EventhRootOfNegative
        else if a < 1.0 then
            1.0 /. (newtonNthRoot n (1.0 /. a))
        else
            let x1 = ref (0.0 +. a) in
            let xn = ref (mPow !x1 (n - 1)) in
            let x2 = ref ((((float_of_int n) -. 1.0) *. !x1 +. a /. !xn)  /. (float_of_int n)) in
            let er = ref (!x1 -. !x2) in
            let counter = ref 0 in
            let not_stop = ref true in
            while !x1 +. !er <> !x1 && !not_stop do
                x1 := !x2 ;
                xn := mPow !x1 (n - 1) ;
                x2 := (((float_of_int n) -. 1.0) *. !x1 +. a /. !xn) /. (float_of_int n) ;
                er := !x1 -. !x2 ;
                if (abs_float !er) < (abs_float nM_EPSILON *. !x1) then
                    not_stop := false ;
                counter := !counter + 1 ;
                if !counter > nMAX_ITER then
                    not_stop := false ;
            done ;
            if !counter >= nMAX_ITER then
                raise TooManyInterations ;
            !x2
    else
        if a = 0.0 then
            raise EventhRootOfNegative
        else
            1.0 /. (newtonNthRoot (-n) a) ;;


 

(* ---- Begin here ---- *)
let x = ref 16.0 in
let u = ref (sqrt !x) in

printf "[ Testing (heronSqrt double) ]--------------------\n" ;
printf "x = %g\n" !x ;
printf "u = sqrt(%g) = %g\n" !x !u ;
let y = ref (heronSqrt !x) in
printf "y = (heronSqrt %g) = %g\n" !x !y ;
printf "y*y = %g\n" (!y *. !y) ;
printf "\n" ;

printf "[ Testing (newtonCbrt double) ]--------------------\n" ;
x := 729000000000.0 ;
printf "x = %g\n" !x ;
printf "exp(log(x)/3.0) = %g\n" (exp ((log !x) /. 3.0)) ;
let w = ref (newtonCbrt !x) in
printf "w = (newtonCbrt %g) = %g\n" !x !w ;
printf "w*w*w = %g\n"  (!w *. !w *. !w)  ;
printf "\n" ;

printf "[ Testing (newtonNthRoot int double) ]--------------------\n" ;
let z = ref (newtonNthRoot 3 !x) in
printf "x = %g\n" !x ;
printf "z = (newtonNthRoot 3 %g) = %g\n" !x !z ;
printf "z*z*z = %g\n" (!z *. !z *. !z) ;
printf "\n" ;

x := 12960000000000000000.0 ;
z := newtonNthRoot 4 !x ;
printf "x = %g\n" !x ;
printf "z = (newtonNthRoot 4 x) = (newtonNthRoot 4 %g) = %g\n" !x !z ;
printf "z*z*z*z = %g\n" (!z *. !z *. !z *. !z) ;
printf "\n" ;

x := 1.0 /. 12960000000000000000.0 ;
z := newtonNthRoot 4 !x ;
printf "x = %g\n" !x ;
printf "exp(log(x)/4.0) = %g\n" (exp ((log !x) /. 4.0)) ;
printf "z = (newtonNthRoot 4 x) = (newtonNthRoot 4 %g) = %g\n" !x !z ;
printf "z*z*z*z = %g\n" (!z *. !z *. !z *. !z) ;
printf "\n" ;


(*
try
    x := -4.0 ;
    printf "[ Test Exception (heronSqrt double) ]--------------------\n" ;
    printf "x = %g\n" !x ;
    printf "Calculating (heronSqrt %g)\n" !x ;
    y := heronSqrt !x ;
    printf "y = (heronSqrt %g) = %g\n" !x !y ;
    printf "y*y = %g\n" (!y *. !y) ;
    printf "\n" ;
with
    | SqrtOfNegative -> printf "Tried to get a square root of a negative number in calculating (heronSqrt %g)\n" !x ;
             printf "\n" ;
    | TooManyInterations -> printf "Too many iterations in calculating (heronSqrt %g)\n" !x ;
             printf "\n" ;

try
    x := -4.0 ;
    printf "[ Test Exception (newtonCbrt double) ]--------------------\n" ;
    printf "x = %g\n" !x ;
    printf "Calculating (newtonCbrt %g)\n" !x ;
    y := newtonCbrt !x ;
    printf "y = newtonCbrt(%g) = %g\n" !x !y ;
    printf "y*y*y = %g\n" (!y *. !y *. !y) ;
    printf "\n" ;
with
    | TooManyInterations -> printf "Too many iterations in calculating (heronSqrt %g)\n" !x ;
             printf "\n" ;
*)

printf "[ Test calculations by powering ]-----------------------------\n" ;
x := 200.0 ;
z := newtonNthRoot 10  !x ;
printf "x = %g\n" !x ;
printf "exp(log(x)/10.0) = %g\n" (exp ((log !x) /. 10.0)) ;
printf "z = (newtonNthRoot 10 x) = (newtonNthRoot 10 %g) = %g\n" !x !z ;
printf "z**10.0 = (mPow z 10) = %g\n" (mPow !z 10) ;
printf "z**10.0 = z ** 10.0 = %g\n" (!z ** (float_of_int 10)) ;
printf "\n" ;

x := 3001.0 ;
z := newtonNthRoot 99 !x ;
printf "x = %g\n" !x ;
printf "exp(log(x)/99.0) = %g\n" (exp ((log !x) /. 99.0)) ;
printf "z = (newtonNthRoot 99 x) = (newtonNthRoot 99 %g) = %g\n" !x !z ;
printf "z**99.0 = (mPow z 99) = %g\n" (mPow !z 99) ;
printf "z**99.0 = z ** 99.0 = %g\n" (!z ** (float_of_int 99)) ;
printf "\n" ;

x := 3001.0 ;
z := newtonNthRoot (-99) !x ;
printf "x = %g\n" !x ;
printf "exp(log(x)/-99.0) = %g\n" (exp ((log !x) /. (-99.0))) ;
printf "z = (newtonNthRoot (-99) x) = (newtonNthRoot (-99) %g) = %g\n" !x !z ;
printf "mPow z (-99) = %g\n" (mPow !z (-99)) ;
printf "1.0/(z**99.0) = %g\n" (1.0 /. (!z ** 99.0)) ;
printf "z**(-99.0) = %g\n" (!z ** (-99.0)) ;
printf "\n" ;


printf "2.1**2.1 = pow(2.1, 2.1) = %g\n" (2.1**2.1) ;
printf "2.1**(-2.1) = pow(2.1, -2.1) = %g\n" (2.1**(-2.1)) ;
printf "2.1**2.1 * 2.1**(-2.1) = pow(2.1, 2.1) * pow(2.1, -2.1) = %g\n" (2.1**2.1 *. 2.1**(-2.1)) ;
printf "2.1**2.1 = exp(2.1*log(2.1)) = %g\n" (exp (2.1 *. (log 2.1))) ;
printf "2.1**(-2.1) = exp(-2.1*log(2.1)) = %g\n" (exp (-2.1 *. (log 2.1))) ;
printf "2.1**2.1 * 2.1**(-2.1) = exp(2.1*log(2.1)) * exp(-2.1*log(2.1)) = %g\n"  (2.1**2.1 *. 2.1**(-2.1)) ;
printf "\n" ;


let k = ref 301 in
x := -1.029 ;
let t1 = ref (nPow !x !k) in
let t2 = ref (gPow !x !k) in
let t3 = ref (mPow !x !k) in
let tt = ref (!x ** (float_of_int !k)) in
printf "%g ** %d = %g\n" !x !k !tt ;
printf "t1 = (nPow %g %d) = %g\n" !x !k !t1 ;
printf "t2 = (gPow %g %d) = %g\n" !x !k !t2 ;
printf "t3 = (mPow %g %d) = %g\n" !x !k !t3 ;
printf "t1 / t2 = %g\n" (!t1 /. !t2) ;
printf "t1 - t2 = %g\n" (!t1 -. !t2) ;
printf "t1 == t2 ? " ;
if t1 = t2 then printf "yes\n"
           else printf "no\n" ;
printf "t1 / t3 = %g\n" (!t1 /. !t3) ;
printf "t1 - t3 = %g\n" (!t1 -. !t3) ;
printf "t1 == t3 ? "  ;
if t1 = t3 then printf "yes\n"
           else printf "no\n" ;
printf "t2 / t3 = %g\n" (!t2 /. !t3) ;
printf "t2 - t3 = %g\n" (!t2 -. !t3) ;
printf "t2 == t3 ? "  ;
if t2 = t3 then printf "yes\n"
           else printf "no\n" ;
printf "\n" ;


printf "Done.\n" ;

 

(*
Output:
[ Testing (heronSqrt double) ]--------------------
x = 16
u = sqrt(16) = 4
y = (heronSqrt 16) = 4
y*y = 16

[ Testing (newtonCbrt double) ]--------------------
x = 7.29e+011
exp(log(x)/3.0) = 9000
w = (newtonCbrt 7.29e+011) = 9000
w*w*w = 7.29e+011

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

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

x = 7.71605e-020
exp(log(x)/4.0) = 1.66667e-005
z = (newtonNthRoot 4 x) = (newtonNthRoot 4 7.71605e-020) = 1.66667e-005
z*z*z*z = 7.71605e-020

[ Test calculations by powering ]-----------------------------
x = 200
exp(log(x)/10.0) = 1.69865
z = (newtonNthRoot 10 x) = (newtonNthRoot 10 200) = 1.69865
z**10.0 = (mPow z 10) = 200
z**10.0 = z ** 10.0 = 200

x = 3001
exp(log(x)/99.0) = 1.08424
z = (newtonNthRoot 99 x) = (newtonNthRoot 99 3001) = 1.08424
z**99.0 = (mPow z 99) = 3001
z**99.0 = z ** 99.0 = 3001

x = 3001
exp(log(x)/-99.0) = 0.922308
z = (newtonNthRoot (-99) x) = (newtonNthRoot (-99) 3001) = 0.922308
mPow z (-99) = 3001
1.0/(z**99.0) = 3001
z**(-99.0) = 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

-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-011
t1 == t2 ? no
t1 / t3 = 1
t1 - t3 = 6.18456e-011
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 )

 

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

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

아래의 소스는 대부분 버전의 Lua 에서 실행되도록 작성된 소스이다.

예외상황 처리에 대해서는 (Java 언어에서 처럼) throw ... 구문으로 예외상황 던지기를 하고, try ... catch ... 구문으로 예와상황 받기를 한다.

<?php

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

 

$NO_EXCEPTION         =  (0);
$CANNOT_EVALUATE   =  (100);

$MAX_ITER = 20000;
$M_EPSILON = 1.0e-15;

 


/**
  * Compute the n-th root of x to a given scale, x > 0.
  */
function nPow($a, $n) {
    if ($n > 0) {
        if ($n == 1)
            return $a;
        else {
            if ($a == 0.0 || $a == 1.0) {
                return $a;
            }
            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 {
                $y = 1.0;
                for ($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.
  */
function gPow($a, $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 {

                $y = 1.0;
                $r = $a;
                $m = 8*4 - 1;            // 8*sizeof(int) - 1;  ignore sign bit, which is MSB
                $one = 1;
                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;
            }
        }
    }
    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.
  */
function mPow($a, $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 {

                $y = 1.0;
                $r = $a;
                $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 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.
  */
function heronSqrt($a) {
    global $MAX_ITER, $M_EPSILON;

    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 {
        $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++;
            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.
  */
function newtonCbrt($a) {
    global $MAX_ITER, $M_EPSILON;

    if ($a == 0.0 || $a == 1.0 || $a == -1.0) {
        return $a;
    }
    else if ($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++;
            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.
  */
function newtonNthRoot($n, $a) {
    global $MAX_ITER, $M_EPSILON;

    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 {
            $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++;
                if ($counter > $MAX_ITER)
                    break;
            }
            if ($counter >= $MAX_ITER) {
                throw new Exception("Inaccurate n-th 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);
        }
    }
}

 

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

printf("[ Testing newtonCbrt(double) ]--------------------\n");
$x = -216.0;
printf("x = %g\n", $x);
printf("-exp(log(-x)/3.0) = %g\n", -exp(log(-$x)/3.0));
$w = newtonCbrt($x);
printf("w = newtonCbrt(%g) = %g\n", $x, $w);
printf("w*w*w = %g\n", $w*$w*$w);
printf("\n");

$x = 729000000000.0;
printf("x = %g\n", $x);
printf("exp(log(x)/3.0) = %g\n", exp(log($x)/3.0));
$w = newtonCbrt($x);
printf("w = newtonCbrt(%g) = %g\n", $x, $w);
printf("w*w*w = %g\n", $w*$w*$w);
printf("\n");

printf("[ Testing newtonNthRoot(int, double) ]--------------------\n");
$z = newtonNthRoot(3, $x);
printf("x = %g\n", $x);
printf("z = newtonNthRoot(3, %g) = %g\n", $x, $z);
printf("z*z*z = %g\n", $z*$z*$z);
printf("\n");

$x = 12960000000000000000.0;
$z = newtonNthRoot(4, $x);
printf("x = %g\n", $x);
printf("z = newtonNthRoot(4, x) = newtonNthRoot(4, %g) = %g\n", $x, $z);
printf("z*z*z*z = %g\n", $z*$z*$z*$z);
printf("\n");

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


try {
    $x = -4.0;
    printf("[ Test Exception heronSqrt(double) ]--------------------\n");
    printf("x = %g\n", $x);
    printf("Calculating heronSqrt(%g)\n" , $x);
    $y = heronSqrt($x);
    printf("y = heronSqrt(%g) = %g\n", $x, $y);
    printf("y*y = %g\n", $y*$y);
    printf("\n");
}
catch (Exception $err) {
    printf("%s\nCaught some exception in calculating heronSqrt(%g)\n", $err->getMessage(), $x);
    printf("\n");
}


try {
    $x = -4.0;
    printf("[ Test Exception in newtonCbrt(double) ]--------------------\n");
    printf("x = %g\n", $x);
    printf("Calculating newtonCbrt(%g)\n", $x);
    $y = newtonCbrt($x);
    printf("y = newtonCbrt(%g) = %g\n", $x, $y);
    printf("y*y*y = %g\n", $y*$y*$y);
    printf("\n");
}
catch (Exception $err) {
    printf("%s\nCaught some exception in calculating newtonCbrt(%g)\n", $err->getMessage(), $x);
    printf("\n");
}


printf("[ Test calculations by powering ]-----------------------------\n");
$x = 200.0;
$z = newtonNthRoot(10, $x);
printf("x = %g\n", $x);
printf("exp(log(x)/10.0) = %g\n", exp(log($x)/10.0));
printf("z = newtonNthRoot(10, x) = newtonNthRoot(10, %g) = %g\n", $x,  $z);
printf("pow(z, 10) = %g\n", pow($z, 10));
printf("\n");

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

$x = 3001.0;
$z = newtonNthRoot(-99, $x);
printf("x = %g\n", $x);
printf("exp(log(x)/-99.0) = %g\n", exp(log($x)/-99.0));
printf("z = newtonNthRoot(-99, x) = newtonNthRoot(-99, %g) = %g\n", $x, $z);
printf("1.0/pow(z, 99) = %g\n", 1.0/pow($z, 99));
printf("\n");


printf("2.1**2.1 = pow(2.1, 2.1) = %g\n", pow(2.1, 2.1));
printf("2.1**(-2.1) = pow(2.1, -2.1) = %g\n", pow(2.1, -2.1));
printf("2.1**2.1 * 2.1**(-2.1) = pow(2.1, 2.1) * pow(2.1, -2.1) = %g\n", pow(2.1, 2.1)*pow(2.1, -2.1));
printf("2.1**2.1 = exp(2.1*log(2.1)) = %g\n", exp(2.1*log(2.1)));
printf("2.1**(-2.1) = exp(-2.1*log(2.1)) = %g\n", exp(-2.1*log(2.1)));
printf("2.1**2.1 * 2.1**(-2.1) = exp(2.1*log(2.1)) * exp(-2.1*log(2.1)) = %g\n", exp(2.1*log(2.1)) * exp(-2.1*log(2.1)));
printf("\n");

 

$k = 301;
$x = -1.029;
$t1 = nPow($x, $k);
$t2 = gPow($x, $k);
$t3 = mPow($x, $k);
printf("t1 = nPow(%g, %d) = %g\n", $x, $k, $t1);
printf("t2 = gPow(%g, %d) = %g\n", $x, $k, $t2);
printf("t3 = mPow(%g, %d) = %g\n", $x, $k, $t3);
printf("t1 / t2 = %g\n", $t1 / $t2);
printf("t1 - t2 = %g\n", $t1 - $t2);
printf("t1 == t2 ? %s\n",  $t1 == $t2 ? "yes" : "no");
printf("t1 / t3 = %g\n", $t1 / $t3);
printf("t1 - t3 = %g\n", $t1 - $t3);
printf("t1 == t3 ? %s\n",  $t1 == $t3 ? "yes" : "no");
printf("t2 / t3 = %g\n", $t2 / $t3);
printf("t2 - t3 = %g\n", $t2 - $t3);
printf("t2 == t3 ? %s\n",  $t2 == $t3 ? "yes" : "no");
printf("\n");


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
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-5
z = newtonNthRoot(4, x) = newtonNthRoot(4, 7.71605e-20) = 1.66667e-5
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

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 )

 

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

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

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

예외상황 처리룰 하기 위해, C 언어에서는 sejmp.h 를 인클루드하고, 예외상황을 던질 때는 longjmp() 함수를, 예외상황을 받을 때는 setjmp() 함수와 switch{ ... } 문을 사용한다.

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


#include <stdio.h>
#include <math.h>
#include <setjmp.h>

jmp_buf g_stkBuf;

#define NO_EXCEPTION          (0)
#define CANNOT_EVALUATE    (100)

#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)
{
    int i;
    double y;

    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 {
                y = 1.0;
                for (i = 0; i < n; i++) {
                    y *= a;
                }
                return y;
            }
        }
    }
    else if (n == 0) {
        return 1.0;
    }
    else {      //  when n < 0
        if (a == 0.0) {
            printf("Negative powering exception of zero.\n");
            longjmp( g_stkBuf, CANNOT_EVALUATE );  // setjmp로 설정한 g_stkBuf(스택정보)를 가지고 해당 시점으로 복귀
            return -1.0;
        }
        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)
{
    int i;
    double y;
    double r;
    int m;
    int one;

    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 {

                y = 1.0;
                r = a;
                m = 8*sizeof(int) - 1;    // 8*sizeof(int) - 1;  ignore sign bit, which is MSB
                one = 1;
                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;
            }
        }
    }
    else if (n == 0) {
        return 1.0;
    }
    else {      //  when n < 0
        if (a == 0.0) {
            printf("Negative powering exception of zero.\n");
            longjmp( g_stkBuf, CANNOT_EVALUATE );  // setjmp로 설정한 g_stkBuf(스택정보)를 가지고 해당 시점으로 복귀
            return -1.0;
        }
        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) {
            printf("Negative powering exception of zero.\n");
            longjmp( g_stkBuf, CANNOT_EVALUATE );  // setjmp로 설정한 g_stkBuf(스택정보)를 가지고 해당 시점으로 복귀
            return -1.0;
        }
        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)
{
    double x1, x2, er;
    int counter;

    if (a < 0.0) {
        printf("Cannot find the sqrt of a negative number.\n");
        longjmp( g_stkBuf, CANNOT_EVALUATE );  // setjmp로 설정한 g_stkBuf(스택정보)를 가지고 해당 시점으로 복귀
                                                 // errNo는 setjmp로 설정한 지점으로 복귀할 때 리턴 값
        return -1.0;
    }
    else if (a == 0.0 || a == 1.0) {
        return a;
    }
    else {
        x1 = a;
        x2 = (x1 + a/x1)/2.0;
        er = x1 - x2;
        counter = 0;
        while (x1 + er != x1) {
            x1 = x2;
            x2 = (x1 + a/x1)/2.0;
            er = x1 - x2;
            if (fabs(er) < fabs(M_EPSILON*x1))
                break;
            counter++;
            if (counter > MAX_ITER)
                break;
        }
        if (counter >= MAX_ITER) {
            printf("Inaccurate sqrt exception by too many iterations.\n");
            longjmp( g_stkBuf, CANNOT_EVALUATE );  // setjmp로 설정한 g_stkBuf(스택정보)를 가지고 해당 시점으로 복귀
            return -1.0;
        }
        return x2;
    }
}

/**
  * Compute the cubic root of x to a given scale, x > 0.
  */
double newtonCbrt(double a)
{
    double x1, x2, er;
    int counter;

    if (a == 0.0 || a == 1.0 || a == -1.0) {
        return a;
    }
    else if (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 (fabs(er) < fabs(M_EPSILON*x1))
                break;
            counter++;
            if (counter > MAX_ITER)
                break;
        }
        if (counter >= MAX_ITER) {
            printf("Inaccurate cbrt exception by too many iterations.\n");
            longjmp( g_stkBuf, CANNOT_EVALUATE );  // setjmp로 설정한 g_stkBuf(스택정보)를 가지고 해당 시점으로 복귀
            return -1.0;
        }
        return x2;
    }
}

/**
  * Compute the n-th root of x to a given scale, x > 0.
  */
double newtonNthRoot(int n, double a)
{
    double x1, x2, xn, er;
    int counter;

    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 {
                printf("Cannot find the even n-th root of a negative number.\n");
                longjmp( g_stkBuf, CANNOT_EVALUATE );  // setjmp로 설정한 g_stkBuf(스택정보)를 가지고 해당 시점으로 복귀
                return -1.0;
            }
        }
        else if (a < 0.0) {
            if (n % 2 == 1)
                return -newtonNthRoot(n, -a);
            else {
                printf("Cannot find the even n-th root of a negative number.\n");
                longjmp( g_stkBuf, CANNOT_EVALUATE );  // setjmp로 설정한 g_stkBuf(스택정보)를 가지고 해당 시점으로 복귀
                return -1.0;
            }
        }
        else if (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 (fabs(er) < fabs(M_EPSILON*x1))
                    break;
                counter++;
                if (counter > MAX_ITER)
                    break;
            }
            if (counter >= MAX_ITER) {
                printf("Inaccurate n-th exception by too many iterations.\n");
                longjmp( g_stkBuf, CANNOT_EVALUATE );  // setjmp로 설정한 g_stkBuf(스택정보)를 가지고 해당 시점으로 복귀
                return -1.0;
            }
            return x2;
        }
    }
    else {
        if (a == 0.0) {
            printf("Cannot find the negative n-th root of zero.\n");
            longjmp( g_stkBuf, CANNOT_EVALUATE );  // setjmp로 설정한 g_stkBuf(스택정보)를 가지고 해당 시점으로 복귀
            return -1.0;
        }
        else {
            return 1.0/newtonNthRoot(-n, a);
        }
    }
}

 

int main()
{
    double x = 16.0;
    double u = sqrt(x);

    int k;
    double t1, t2, t3;
    double y, z, w;

    int errNo = NO_EXCEPTION;

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

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

    x = 729000000000.0;
    printf("x = %g\n", x);
    printf("exp(log(x)/3.0) = %g\n", exp(log(x)/3.0));
    w = newtonCbrt(x);
    printf("w = newtonCbrt(%g) = %g\n", x, w);
    printf("w*w*w = %g\n", w*w*w);
    printf("\n");

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

    x = 12960000000000000000.0;
    z = newtonNthRoot(4, x);
    printf("x = %g\n", x);
    printf("z = newtonNthRoot(4, x) = newtonNthRoot(4, %g) = %g\n", x, z);
    printf("z*z*z*z = %g\n", z*z*z*z);
    printf("\n");

    x = 1.0/12960000000000000000.0;
    z = newtonNthRoot(4, x);
    printf("x = %g\n", x);
    printf("exp(log(x)/4.0) = %g\n", exp(log(x)/4.0));
    printf("z = newtonNthRoot(4, x) = newtonNthRoot(4, %g) = %g\n", x, z);
    printf("z*z*z*z = %g\n", z*z*z*z);
    printf("\n");


    errNo = setjmp( g_stkBuf );  // longjmp로 복귀할 지점을 설정, g_stkBuf에 스택정보 저장
    switch ( errNo )
    {
        case NO_EXCEPTION:
            x = -4.0;
            printf("[ Test Exception heronSqrt(double) ]--------------------\n");
            printf("x = %g\n", x);
            printf("Calculating heronSqrt(%g)\n" , x);
            y = heronSqrt(x);
            printf("y = heronSqrt(%g) = %g\n", x, y);
            printf("y*y = %g\n", y*y);
            printf("\n");
            break;
        case CANNOT_EVALUATE:
            printf("Caught some exception in calculating heronSqrt(%g)\n", x);
            printf("\n");
            errNo = NO_EXCEPTION;
            break;
    }


    errNo = setjmp( g_stkBuf );  // longjmp로 복귀할 지점을 설정, g_stkBuf에 스택정보 저장
    switch ( errNo )
    {
        case NO_EXCEPTION:
            x = -4.0;
            printf("[ Test Exception in newtonCbrt(double) ]--------------------\n");
            printf("x = %g\n", x);
            printf("Calculating newtonCbrt(%g)\n", x);
            y = newtonCbrt(x);
            printf("y = newtonCbrt(%g) = %g\n", x, y);
            printf("y*y*y = %g\n", y*y*y);
            printf("\n");
            break;
        case CANNOT_EVALUATE:
            printf("Caught some exception in calculating newtonCbrt(%g)\n", x);
            printf("\n");
            errNo = NO_EXCEPTION;
            break;
    }

    printf("[ Test calculations by powering ]-----------------------------\n");
    x = 200.0;
    z = newtonNthRoot(10, x);
    printf("x = %g\n", x);
    printf("exp(log(x)/10.0) = %g\n", exp(log(x)/10.0));
    printf("z = newtonNthRoot(10, x) = newtonNthRoot(10, %g) = %g\n", x,  z);
    printf("pow(z, 10) = %g\n", pow(z, 10));
    printf("\n");

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

    x = 3001.0;
    z = newtonNthRoot(-99, x);
    printf("x = %g\n", x);
    printf("exp(log(x)/-99.0) = %g\n", exp(log(x)/-99.0));
    printf("z = newtonNthRoot(-99, x) = newtonNthRoot(-99, %g) = %g\n", x, z);
    printf("1.0/pow(z, 99) = %g\n", 1.0/pow(z, 99));
    printf("\n");

    printf("2.1**2.1 = pow(2.1, 2.1) = %g\n", pow(2.1, 2.1));
    printf("2.1**(-2.1) = pow(2.1, -2.1) = %g\n", pow(2.1, -2.1));
    printf("2.1**2.1 * 2.1**(-2.1) = pow(2.1, 2.1) * pow(2.1, -2.1) = %g\n", pow(2.1, 2.1)*pow(2.1, -2.1));
    printf("2.1**2.1 = exp(2.1*log(2.1)) = %g\n", exp(2.1*log(2.1)));
    printf("2.1**(-2.1) = exp(-2.1*log(2.1)) = %g\n", exp(-2.1*log(2.1)));
    printf("2.1**2.1 * 2.1**(-2.1) = exp(2.1*log(2.1)) * exp(-2.1*log(2.1)) = %g\n", exp(2.1*log(2.1)) * exp(-2.1*log(2.1)));
    printf("\n");


    k = 301;
    x = -1.029;
    t1 = nPow(x, k);
    t2 = gPow(x, k);
    t3 = mPow(x, k);
    printf("t1 = nPow(%g, %d) = %g\n", x, k, t1);
    printf("t2 = gPow(%g, %d) = %g\n", x, k, t2);
    printf("t3 = mPow(%g, %d) = %g\n", x, k, t3);
    printf("t1 / t2 = %g\n", t1 / t2);
    printf("t1 - t2 = %g\n", t1 - t2);
    printf("t1 == t2 ? %s\n",  t1 == t2 ? "yes" : "no");
    printf("t1 / t3 = %g\n", t1 / t3);
    printf("t1 - t3 = %g\n", t1 - t3);
    printf("t1 == t3 ? %s\n",  t1 == t3 ? "yes" : "no");
    printf("t2 / t3 = %g\n", t2 / t3);
    printf("t2 - t3 = %g\n", t2 - t3);
    printf("t2 == t3 ? %s\n",  t2 == t3 ? "yes" : "no");
    printf("\n");

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

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

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

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

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

[ 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

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-011
t1 == t2 ? no
t1 / t3 = 1
t1 - t3 = 6.18456e-011
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 )

 

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

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

 

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

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

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

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

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

아래의 소스 첫 부분에

        val MAX_ITER = 20000
        val M_EPSILON = 1.0e-15
       

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

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

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

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

import java.io._

object testNthRoot {

    val MAX_ITER = 20000
    val M_EPSILON = 1.0e-15
   
   
    //
    // Compute the n-th root of x to a given scale, x > 0.
    //
    def nPow(a: Double, n: Int) : Double = {
        if (n > 0) {
            if (n == 1)
                return a
            else {
                if (a == 0.0 || a == 1.0)
                    return a
                else if (a == -1.0) {
                    if ((n % 2) == 1)
                        return -1.0
                    else
                        return 1.0
                }
                else if (a < 0.0) {
                    if ((n % 2) == 1)
                        return -nPow(-a, n)
                    else
                        return nPow(-a, n)
                }
                else {
                    var y = 1.0
                    for (i <- 1 to n) {
                        y = y * a
                    }
                    return y
                }
            }
        }
        else if (n == 0)
            return 1.0
        else {     //  when n < 0
            if (a == 0.0) {
                throw new RuntimeException("Negative powering exception of zero.")
            }
            else {
                if (n == -1)
                    return 1.0/a
                else
                   return 1.0/nPow(a, -n)
            }
        }
    }


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


    //
    // Compute the n-th root of x to a given scale, x > 0.
    //
    def mPow(a: Double, n: Int) : Double = {
        if (n > 0) {
            if (n == 1)
                return a
            else {
                if (a == 0.0 || a == 1.0)
                    return a
                else if (a == -1.0) {
                    if ((n % 2) == 1)
                        return -1.0
                    else
                        return 1.0
                }
                else if (a < 0.0) {
                    if ((n % 2) == 1)
                        return -mPow(-a, n)
                    else
                        return mPow(-a, n)
                }
                else {
                    var y = 1.0
                    var r = a
                    var m = n
                    while (m > 0) {
                        if ((m & 0x1) == 1) {
                            y = y * r
                        }
                        r = r*r
                        m = (m >> 1)
                    }
                    return y
                }
            }
        }
        else if (n == 0)
            return 1.0
        else {     //  when n < 0
            if (a == 0.0) {
                throw new RuntimeException("Negative powering exception of zero.")
            }
            else {
                if (n == -1)
                    return 1.0/a
                else
                    return 1.0/mPow(a, -n)
            }
        }
    }


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


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


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

 

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

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

            x = 729000000000.0
            printf("x = %g\n", x)
            printf("exp(log(x)/3.0) = %g\n", scala.math.exp(scala.math.log(x)/3.0))
            var w = newtonCbrt(x)
            printf("w = newtonCbrt(%g) = %g\n", x, w)
            printf("w*w*w = %g\n", w*w*w)
            printf("\n")

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

            x = 12960000000000000000.0
            z = newtonNthRoot(4, x)
            printf("x = %g\n", x)
            printf("z = newtonNthRoot(4, x) = newtonNthRoot(4, %g) =  %g\n", x, z)
            printf("z*z*z*z = %g\n", z*z*z*z)
            printf("\n")

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


            try {
                x = -4.0
                printf("[ Test Exception heronSqrt(double) ]--------------------\n")
                printf("x = %g\n", x)
                printf("Calculating heronSqrt(%g)\n", x)
                y = heronSqrt(x)
                printf("y = heronSqrt(%g) = %g\n", x, y)
                printf("y*y = %g\n", y*y)
                printf("\n")
            }
            catch {
                case ex : RuntimeException =>
                        printf("%s\nCaught some exception in calculating heronSqrt(%g)\n", ex.getMessage(),  x);
                        printf("\n")
            }

            try {
                x = -4.0
                printf("[ Test Exception newtonCbrt(double) ]--------------------\n")
                printf("x = %g\n", x)
                printf("Calculating newtonCbrt(%g)\n", x)
                y = newtonCbrt(x)
                printf("y = newtonCbrt(%g) = %g\n", x, y)
                printf("y*y*y = %g\n", y*y*y)
                printf("\n")
            }
            catch {
                case ex : RuntimeException =>
                        printf("%s\nCaught some exception in calculating newtonCbrtrt(%g)\n", ex.getMessage(),  x);
                        printf("\n")
            }


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

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

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

            printf("2.1**2.1 = pow(2.1, 2.1) = %g\n", scala.math.pow(2.1, 2.1))
            printf("2.1**(-2.1) = pow(2.1, -2.1) = %g\n", scala.math.pow(2.1, -2.1))
            printf("2.1**2.1 * 2.1**(-2.1) = pow(2.1, 2.1) * pow(2.1, -2.1) = %g\n",  scala.math.pow(2.1, 2.1) * scala.math.pow(2.1, -2.1))
            printf("2.1**2.1 = exp(2.1*log(2.1)) = %g\n", scala.math.exp(2.1*scala.math.log(2.1)))
            printf("2.1**(-2.1) = exp(-2.1*log(2.1)) = %g\n", scala.math.exp(-2.1*scala.math.log(2.1)))
            printf("2.1**2.1 * 2.1**(-2.1) = exp(2.1*log(2.1)) * exp(-2.1*log(2.1)) = %g\n", scala.math.pow(2.1, 2.1) * scala.math.pow(2.1, -2.1))
            printf("\n")


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

            printf("Done.\n")
    }

}


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

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

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

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

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

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

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

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

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

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

2.1**2.1 = pow(2.1, 2.1) = 4.74964
2.1**(-2.1) = pow(2.1, -2.1) = 0.210542
2.1**2.1 * 2.1**(-2.1) = pow(2.1, 2.1) * pow(2.1, -2.1) = 1.00000
2.1**2.1 = exp(2.1*log(2.1)) = 4.74964
2.1**(-2.1) = exp(-2.1*log(2.1)) = 0.210542
2.1**2.1 * 2.1**(-2.1) = exp(2.1*log(2.1)) * exp(-2.1*log(2.1)) = 1.00000

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

Done.
*/

 

 

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 )

 

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

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

 

1. F# 언어는 닷넷을 목표로 하고 닷넷 위에서 동작하는 (절차적 프로그래밍도 지원하는) 함수형 언어이다.

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

2. F# 언어는 (Python 언어 처럼 소스에 들여쓰기 규칙을 철저히 지킨다.

3. F# 언어에서는 상수든 변수든 함수든 처음 선언할 때는 구문 선두에 let 예약어를 반드시 붙인다.

     let a = 2.7                // a 상수
     let mutable b = 5.3    // b 는 변수 
     let g x = x*x             // g 는 함수

아래의 소스 첫 부분에

        let MAX_ITER = 20000
        let M_EPSILON = 1.0e-15

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

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

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

 

 (*
 * Filename: testNthRoot.fs
 *
 *            Approximate square roots, cubic roots and n-th roots of a given number.
 *
 * Compile: fsc --codepage:949 testNthRoot.fs
 * Execute: testNthRoot
 *
 * Date: 2013. 1. 7.
 * Copyright (c) 2013 PH Kim  (pkim __AT__ scripts.pe.kr)
 *)

# light

let MAX_ITER = 20000
let M_EPSILON = 1.0e-15


//
// Compute the n-th root of x to a given scale, x > 0.
//
let rec nPow(a: double, n: int) : double =
    if n > 0 then
        if n = 1 then
            a
        else
            if a = 0.0 || a = 1.0 then
                a
            elif a = -1.0 then
                if n % 2 = 1 then
                    -1.0
                else
                    1.0
            elif a < 0.0 then
                if n % 2 = 1 then
                    -nPow(-a, n)
                else
                    nPow(-a, n)
            else
                let mutable y = 1.0
                for i = 1 to n do
                    y <- y * a
                y
    elif n = 0 then
        1.0
    else      //  when n < 0
        if a = 0.0 then
            raise (new System.InvalidOperationException("Negative powering exception of zero."))
        else
            if n = -1 then
                1.0/a
            else
               1.0/nPow(a, -n)


//
// Compute the n-th root of x to a given scale, x > 0.
//
let rec gPow(a: double, n: int) : double =
    if n > 0 then
        if n = 1 then
            a
        else
            if a = 0.0 || a = 1.0 then
                a
            elif a = -1.0 then
                if n % 2 = 1 then
                    -1.0
                else
                    1.0
            elif a < 0.0 then
                if n % 2 = 1 then
                    -gPow(-a, n)
                else
                    gPow(-a, n)
            else
                let mutable y = 1.0
                let mutable r = a
                let mutable m = 8*4 - 1            //  8*sizeof(int) - 1;
                let mutable one = 1
                for i = 0 to m do
                    if (n &&& one) = 0 then
                        y <- y * 1.0
                    else
                        y <- y * r
                    r <- r*r
                    one <- (one <<< 1)
                y
    elif n = 0 then
        1.0
    else      //  when n < 0
        if a = 0.0 then
            raise (new System.InvalidOperationException("Negative powering exception of zero."))
        else
            if n = -1 then
                1.0/a
            else
                1.0/gPow(a, -n)

 

//
// Compute the n-th root of x to a given scale, x > 0.
//
let rec mPow(a: double, n: int) : double =
    if n > 0 then
        if n = 1 then
            a
        else
            if a = 0.0 || a = 1.0 then
                a
            elif a = -1.0 then
                if n % 2 = 1 then
                    -1.0
                else
                    1.0
            elif a < 0.0 then
                if n % 2 = 1 then
                    -mPow(-a, n)
                else
                    mPow(-a, n)
            else
                let mutable y = 1.0
                let mutable r = a
                let mutable m = n 
                while m > 0 do
                    if (m &&& 0x1) = 1 then
                        y <- y * r
                    r <- r*r
                    m <- (m >>> 1)
                y
    elif n = 0 then
        1.0
    else      //  when n < 0
        if a = 0.0 then
            raise (new System.InvalidOperationException("Negative powering exception of zero."))
        else
            if n = -1 then
                1.0/a
            else
                1.0/mPow(a, -n)

 

//
// Compute the square root of x to a given scale, x > 0.
//
let rec heronSqrt(a: double) : double =
    if a < 0.0 then
        raise (new System.InvalidOperationException("Cannot find the sqrt of a negative number."))
    elif a = 0.0 || a = 1.0 then
        a
    else
        let mutable x1 = a
        let mutable x2 = (x1 + a/x1)/2.0
        let mutable er = x1 - x2
        let mutable counter = 0
        let mutable not_stop = true
        while x1 + er <> x1 && not_stop do
            x1 <- x2
            x2 <- (x1 + a/x1)/2.0
            er <- x1 - x2
            if abs(er) < abs(M_EPSILON*x1) then
                not_stop <- false
            counter <- counter + 1
            if counter > MAX_ITER then
                not_stop <- false
        if counter >= MAX_ITER then
                raise (new System.InvalidOperationException("Inaccurate sqrt exception by too many iterations."))
        x2


//
// Compute the cubic root of x to a given scale, x > 0.
//
let rec newtonCbrt(a: double) : double =
    if a = 0.0 || a = 1.0 || a = -1.0 then
        a
    elif a < 0.0 then
        -newtonCbrt(-a)
    else
        let mutable x1 = a
        let mutable x2 = (2.0*x1 + a/(x1*x1))/3.0
        let mutable er = x1 - x2
        let mutable counter = 0
        let mutable not_stop = true
        while x1 + er <> x1 && not_stop do
            x1 <- x2
            x2 <- (2.0*x1 + a/(x1*x1))/3.0
            er <- x1 - x2
            if abs(er) < abs(M_EPSILON*x1) then
                not_stop <- false
            counter <- counter + 1
            if counter > MAX_ITER then
                not_stop <- false
        if counter >= MAX_ITER then
            raise (new System.InvalidOperationException("Inaccurate cbrt exception by too many iterations."))
        x2


//
// Compute the n-th root of x to a given scale, x > 0.
//
let rec newtonNthRoot(n: int, a: double) : double =
    if n = 0 then
        1.0
    elif n = 1 then
        a
    elif n > 0 then
        if a = 0.0 || a = 1.0 then
            a
        elif a = -1.0 then
            if n % 2 = 1 then
                a
            else
                raise (new System.InvalidOperationException("Cannot find the even n-th root of a negative number."))
        elif a < 0.0 then
            if n % 2 = 1 then
                -newtonNthRoot(n, -a)
            else
                raise (new System.InvalidOperationException("Cannot find the n-th root of a negative number."))
        elif a < 1.0 then
            1.0/newtonNthRoot(n, 1.0/a)
        else
            let mutable x1 = a
            let mutable xn = mPow(x1, n - 1)
            let mutable x2 = ((double n - 1.0)*x1 + a/xn)/(double n)
            let mutable er = x1 - x2
            let mutable counter = 0
            let mutable not_stop = true
            while x1 + er <> x1 && not_stop do
                x1 <- x2
                xn <- mPow(x1, n - 1)
                x2 <- ((double n - 1.0)*x1 + a/xn)/(double n)
                er <- x1 - x2
                if abs(er) < abs(M_EPSILON*x1) then
                    not_stop <- false
                counter <- counter + 1
                if counter > MAX_ITER then
                    not_stop <- false
            if counter >= MAX_ITER then
                raise (new System.InvalidOperationException("Inaccurate n-th root  exception by too many iterations."))
            x2
    else
        if a = 0.0 then
            raise (new System.InvalidOperationException("Cannot find the negative n-th root of zero."))
        else
            1.0/newtonNthRoot(-n, a)


 

 

let mutable x = 16.0
let mutable u = System.Math.Sqrt(x)

printfn "[ Testing heronSqrt(double) ]--------------------"
printfn "x = %g" x
printfn "u = sqrt(%g) = %g" x u
let mutable y = heronSqrt(x)
printfn "y = heronSqrt(%g) = %g" x y
printfn "y*y = %g" (y*y)
printfn ""

x <- 729000000000.0
printfn "x = %g" x
printfn "exp(log(x)/3.0) = %g" (System.Math.Exp (System.Math.Log x)/3.0)
let mutable w = (newtonCbrt x)
printfn "w = newtonCbrt(%g) = %g" x w
printfn "w*w*w = %g"  (w*w*w)
printfn ""

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

x <- 12960000000000000000.0
z <- newtonNthRoot(4, x)
printfn "x = %g" x
printfn "z = newtonNthRoot(4, x) = newtonNthRoot(4, %g) =  %g" x z
printfn "z*z*z*z = %g" (z*z*z*z)
printfn ""

x <- 1.0/12960000000000000000.0
z <- newtonNthRoot(4, x)
printfn "x = %g" x
printfn "exp(log(x)/4.0) = %g" (System.Math.Exp (System.Math.Log x)/4.0)
printfn "z = newtonNthRoot(4, x) = newtonNthRoot(4, %g) =  %g" x z
printfn "z*z*z*z = %g" (z*z*z*z)
printfn ""


try
    x <- -4.0
    printfn "[ Test Exception heronSqrt(double) ]--------------------"
    printfn "x = %g" x
    printfn "Calculating heronSqrt(%g)" x
    y <- heronSqrt(x)
    printfn "y = heronSqrt(%g) = %g" x y
    printfn "y*y = %g" (y*y)
    printfn ""
with
    | err -> printfn "%s\nCaught some exception in calculating heronSqrt(%g)" err.Message x
             printfn ""

try
    x <- -4.0
    printfn "[ Test Exception newtonCbrt(double) ]--------------------"
    printfn "x = %g" x
    printfn "Calculating newtonCbrt(%g)" x
    y <- newtonCbrt(x)
    printfn "y = newtonCbrt(%g) = %g" x y
    printfn "y*y*y = %g" (y*y*y)
    printfn ""
with
    | err -> printfn "%s\nCaught some exception in calculating newtonCbrtrt(%g)" err.Message x
             printfn ""


printfn "[ Test calculations by powering ]-----------------------------"
x <- 200.0
z <- newtonNthRoot(10,  x)
printfn "x = %g" x
printfn "exp(log(x)/10.0) = %g"  ((System.Math.Exp (System.Math.Log x))/10.0)
printfn "z = newtonNthRoot(10, x) = newtonNthRoot(10, %g) = %g" x z
printfn "z**10.0 = pow(z, 10.0) = %g" (z**10.0)
printfn ""

x <- 3001.0
z <- newtonNthRoot(99, x)
printfn "x = %g" x
printfn "exp(log(x)/99.0) = %g" (System.Math.Exp(System.Math.Log(x)/99.0))
printfn "z = newtonNthRoot(99, x) = newtonNthRoot(99, %g) = %g" x z
printfn "z**99.0 = pow(z, 99) = %g" (z**99.0)
printfn ""

x <- 3001.0
z <- newtonNthRoot(-99, x)
printfn "x = %g" x
printfn "exp(log(x)/-99.0) = %g" (System.Math.Exp(System.Math.Log(x)/(-99.0)))
printfn "z = newtonNthRoot(-99, x) = newtonNthRoot(-99, %g) = %g" x z
printfn "1.0/z**99.0 = 1.0/pow(z, 99) = %g" (1.0/z**99.0)
printfn ""

printfn "2.1**2.1 = pow(2.1, 2.1) = %g" (2.1**2.1)
printfn "2.1**(-2.1) = pow(2.1, -2.1) = %g" (2.1**(-2.1))
printfn "2.1**2.1 * 2.1**(-2.1) = pow(2.1, 2.1) * pow(2.1, -2.1) = %g" (2.1**2.1 * 2.1**(-2.1))
printfn "2.1**2.1 = exp(2.1*log(2.1)) = %g" (System.Math.Exp(2.1*System.Math.Log(2.1)))
printfn "2.1**(-2.1) = exp(-2.1*log(2.1)) = %g" (System.Math.Exp(-2.1*System.Math.Log(2.1)))
printfn "2.1**2.1 * 2.1**(-2.1) = exp(2.1*log(2.1)) * exp(-2.1*log(2.1)) = %g"  (2.1**2.1 * 2.1**(-2.1))
printfn ""


let mutable k = 301
x <- -1.029
let mutable t1 = nPow(x, k)
let mutable t2 = gPow(x, k)
let mutable t3 = mPow(x, k)
let mutable tt =  System.Math.Pow(x, double k)
printfn "System.Math.Pow(%g, %d) = %g" x k tt
printfn "t1 = nPow(%g, %d) = %g" x k t1
printfn "t2 = gPow(%g, %d) = %g" x k t2
printfn "t3 = mPow(%g, %d) = %g" x k t3
printfn "t1 / t2 = %g" (t1 / t2)
printfn "t1 - t2 = %g" (t1 - t2)
printf "t1 == t2 ? "
if t1 = t2 then printfn "yes"
           else printfn "no"
printfn "t1 / t3 = %g" (t1 / t3)
printfn "t1 - t3 = %g" (t1 - t3)
printf "t1 == t3 ? "
if t1 = t3 then printfn "yes"
           else printfn "no"
printfn "t2 / t3 = %g" (t2 / t3)
printfn "t2 - t3 = %g" (t2 - t3)
printf "t2 == t3 ? "
if t2 = t3 then printfn "yes"
           else printfn "no"
printfn ""


printfn "Done."

 

(*
Output:
[ Testing heronSqrt(double) ]--------------------
x = 16
u = sqrt(16) = 4
y = heronSqrt(16) = 4
y*y = 16

x = 7.29e+11
exp(log(x)/3.0) = 2.43e+11
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.92901e-20
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 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) = 20
z = newtonNthRoot(10, x) = newtonNthRoot(10, 200) = 1.69865
z**10.0 = pow(z, 10.0) = 200

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

x = 3001
exp(log(x)/-99.0) = 0.922308
z = newtonNthRoot(-99, x) = newtonNthRoot(-99, 3001) = 0.922308
1.0/z**99.0 = 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

System.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 )

 

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

 

 

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 )

 

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

 

 

 

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 + 1/(x**(n - 1))) / n    를 이용

n = 2 인 경우에는 Newton-Raphson 의 방법이 Heron 의 방법과 동일하다.

(참조. http://en.wikipedia.org/wiki/Newton's_method )

 

Groovy 언어로는 Java 언어 처럼 java.lang 패키지의 지수 연산 함수 Math.pow(밑수, 지수) 를 사용하면 된다. 아니면 Python 언어 처럼 (밑수)**(지수) 의 형식의 구문을 사용해도 된다. 하지만 차후 필요한 데가 있을 것 같아서 이와 유사한 n 제곱 함수와 n 제곱근 함수를 구현해 보았다.

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

아래의 소스는 대부분 버전의 Groovy 에서 실행되도록 작성된 소스이다.

Groovy 언어는 내부적으로 수치를 다룰 때 int 타입과 BigInteger 타입 간에 그리고 double 타입과 BigDecimal 타입 간에 자동으로 변환시키기 때문에 수치 리터럴을 파싱하고 수치를 출력하는데 Java 와 다를 수 있다. 예를 들어, 아래의 소스 중에서

//---------------------------
/// x = 1.0/12960000000000000000.0            //  Groovy fails to parse this. Is It a bug of Groovy?
x = 1.0/(12960000000000000000.0).toDouble()     // Okay!! This should work in Groovy.
z = newtonNthRoot(4, x)

로 작성된 부분은 Groovy 가  수치 리터럴 12960000000000000000.0 은 제대로 파싱하면서 그 역수 1.0/12960000000000000000.0 은 0,0 으로 인식하는 문제가 잇어 BigDecimal 을 명시적으로 double 타입으로 타입변환하는 부분 (수치).toDouble() 을 추가하여

x = 1.0/(12960000000000000000.0).toDouble()

로 작성하게 되었다.

 

/*
 * Filename: testNthRoot.groovy
 *
 *            Approximate square roots, cubic roots and n-th roots of a given number.
 *
 * Execute: groovy testNthRoot.groovy
 *
 * Date: 2013. 1. 7.
 * Copyright (c) 2013 PH Kim  (pkim __AT__ scripts.pe.kr)
 */
 
MAX_ITER = 20000
M_EPSILON = 1.0e-15

//
// Compute the n-th root of x to a given scale, x > 0.
//
def 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 {
                def y = 1.0
                for (i in 0..<n) {
                    y *= a
                }
                return y
            }
        }
    }
    else if (n == 0)
        return 1.0
    else {      //  when n < 0
        if (a == 0.0)
            throw new RuntimeException("Negative powering exception of zero.");
        else {
            if (n == -1)
                return 1.0/a
            else
                return 1.0/nPow(a, -n)
        }
    }
}


//
// Compute the n-th root of x to a given scale, x > 0.
//
def gPow(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 {
                def y = 1.0
                def r = a
                def m = 8*4 - 1           //  8*sizeof(int) - 1;
                def one = 1
                for (i in 0..m) {
                    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.
//
def 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 {
                def y = 1.0
                def r = a
                def 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.
//
def 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 {
        def x1 = a
        def x2 = (x1 + a/x1)/2.0
        def er = x1 - x2
        def 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 += 1
            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.
//
def newtonCbrt(double a) {
    if (a == 0.0 || a == 1.0 || a == -1.0)
        return a
    else if (a < 0.0) {
        return -newtonCbrt(-a)
    }
    else {
        def x1 = a
        def x2 = (2.0*x1 + a/(x1*x1))/3.0
        def er = x1 - x2
        def 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 += 1
            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.
//
def 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 {
            def x1 = a
            def xn = mPow(x1, n - 1)
            def x2 = ((n - 1)*x1 + a/xn)/n
            def er = x1 - x2
            def 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 += 1
                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)
    }
}

 

x = 16.0
u = Math.sqrt(x)

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

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

x = 729000000000.0
printf("x = %g\n", x)
printf("exp(log(x)/3.0) = %g\n", Math.exp(Math.log(x)/3.0))
def w = newtonCbrt(x)
printf("w = newtonCbrt(%g) = %g\n", x, w)
printf("w*w*w = %g\n", w*w*w)
print "\n"

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

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

//---------------------------
/// x = 1.0/12960000000000000000.0     // Groovy fails to parse this. Is It a bug of Groovy?
x = 1.0/(12960000000000000000.0).toDouble()   // Okay!! This should work in Groovy.
z = newtonNthRoot(4, x)
printf("x = %g\n", x)
printf("exp(log(x)/4.0) = %g\n", Math.exp(Math.log(x)/4.0))
printf("z = newtonNthRoot(4, x) = newtonNthRoot(4, %g) =  %g\n", x, z)
printf("z*z*z*z = %g\n", z*z*z*z)
print "\n"


try {
    x = -4.0
    print "[ Test Exception heronSqrt(double) ]--------------------\n"
    printf("x = %g\n", x)
    printf("Calculating heronSqrt(%g)\n", x)
    y = heronSqrt(x)
    printf("y = heronSqrt(%g) = %g\n", x, y)
    printf("y*y = %g\n", y*y)
    print "\n"
}
catch (RuntimeException ex) {
    printf("%s\nCaught some exception in calculating heronSqrt(%g)\n", ex.getMessage(), x)
    print "\n"
}


try {
    x = -4.0
    print "[ Test Exception in newtonCbrt(double) ]--------------------\n"
    printf("x = %g\n", x)
    printf("Calculating newtonCbrt(%g)\n", x)
    y = newtonCbrt(x)
    printf("y = newtonCbrt(%g) = %g\n", x, y)
    printf("y*y*y = %g\n", y*y*y)
    print "\n"
}
catch (RuntimeException ex) {
    printf("%s\nCaught some exception in calculating newtonCbrt(%g)\n", ex.getMessage(), x)
    print "\n"
}


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

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

x = 3001.0
z = newtonNthRoot(-99, x)
printf("x = %g\n", x)
printf("exp(log(x)/-99.0) = %g\n", Math.exp(Math.log(x)/-99.0))
printf("z = newtonNthRoot(-99, x) = newtonNthRoot(-99, %g) = %g\n", x, z)
printf("1.0/z**99 = 1.0/pow(z, 99) = %g\n", 1.0/z**99)
print "\n"

printf("2.1**2.1 = pow(2.1, 2.1) = %g\n", 2.1**2.1)
printf("2.1**(-2.1) = pow(2.1, -2.1) = %g\n", 2.1**(-2.1))
printf("2.1**2.1 * 2.1**(-2.1) = pow(2.1, 2.1) * pow(2.1, -2.1) = %g\n", ((2.1** 2.1)*2.1**( -2.1)))
printf("2.1**2.1 = exp(2.1*log(2.1)) = %g\n", Math.exp(2.1*Math.log(2.1)))
printf("2.1**(-2.1) = exp(-2.1*log(2.1)) = %g\n", Math.exp(-2.1*Math.log(2.1)))
printf("2.1**2.1 * 2.1**(-2.1) = exp(2.1*log(2.1)) * exp(-2.1*log(2.1)) = %g\n", (Math.exp(2.1*Math.log(2.1)) * Math.exp(-2.1*Math.log(2.1))))
print "\n"


k = 301
x = -1.029
t1 = nPow(x, k)
t2 = gPow(x, k)
t3 = mPow(x, k)
printf("%g**%d = %g\n", x, k,  x**k )   // pow(x, k)
printf("t1 = nPow(%g, %d) = %g\n", x, k,  t1)
printf("t2 = gPow(%g, %d) = %g\n", x, k,  t2)
printf("t3 = mPow(%g, %d) = %g\n", x, k,  t3)
printf("t1 / t2 = %g\n", t1 / t2)
printf("t1 - t2 = %g\n", t1 - t2)
print "t1 == t2 ? "
if (t1 == t2)
    print "yes\n"
else
   print "no\n"
printf("t1 / t3 = %g\n", t1 / t3)
printf("t1 - t3 = %g\n", t1 - t3)
print "t1 == t3 ? "
if (t1 == t2)
    print "yes\n"
else
   print "no\n"
printf("t1 / t3 = %g\n", t1 / t3)
printf("t1 - t3 = %g\n", t1 - t3)
print "t1 == t3 ? "
if (t1 == t3)
    print "yes\n"
else
   print "no\n"
printf("t2 / t3 = %g\n", t2 / t3)
printf("t2 - t3 = %g\n", t2 - t3)
print "t2 == t3 ? "
if (t2 == t2)
    print "yes\n"
else
   print "no\n"
print "\n"


print "Done.\n"

 

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

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

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

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

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

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

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

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

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

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

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

2.1**2.1 = pow(2.1, 2.1) = 4.74964
2.1**(-2.1) = pow(2.1, -2.1) = 0.210542
2.1**2.1 * 2.1**(-2.1) = pow(2.1, 2.1) * pow(2.1, -2.1) = 1.00000
2.1**2.1 = exp(2.1*log(2.1)) = 4.74964
2.1**(-2.1) = exp(-2.1*log(2.1)) = 0.210542
2.1**2.1 * 2.1**(-2.1) = exp(2.1*log(2.1)) * exp(-2.1*log(2.1)) = 1.00000

-1.02900**301 = -5457.93
t1 = nPow(-1.02900, 301) = -5457.93
t2 = gPow(-1.02900, 301) = -5457.93
t3 = mPow(-1.02900, 301) = -5457.93
t1 / t2 = 1.00000
t1 - t2 = 6.18456e-11
t1 == t2 ? no
t1 / t3 = 1.00000
t1 - t3 = 6.18456e-11
t1 == t3 ? no
t1 / t3 = 1.00000
t1 - t3 = 6.18456e-11
t1 == t3 ? no
t2 / t3 = 1.00000
t2 - t3 = 0.00000
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 )

 

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

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

아래의 소스는 대부분 버전의 Lua 에서 실행되도록 작성된 소스이다.

#!/usr/bin/env lua

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

local MAX_ITER = 20000
local M_EPSILON = 1.0e-15

local BASE36 = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ"

function itoa(num, radix) 
    local isNegative = false 
    if num < 0 then 
        isNegative = true 
        num = -num 
    end

    local arr = {} 
    local q, r = num, 0 
    while q >= radix do 
        r = math.fmod(q, radix) 
        q = math.floor(q / radix) 
        table.insert(arr, string.sub(BASE36, r+1, r+1)) 
    end

    table.insert(arr, string.sub(BASE36, q+1, q+1)) 
    if isNegative then 
        table.insert(arr, "-") 
    end 
    local ret = "" 
    for i = #arr, 1, -1 do 
        ret = ret .. arr[i] 
    end 
    return ret
end


--
-- Compute the n-th root of x to a given scale, x > 0.
--
function nPow(a, n)
    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 % 2 == 1 then
                    return -1.0
                else
                    return 1.0
                end
            elseif a < 0.0 then
                if n % 2 == 1 then
                    return -nPow(-a, n)
                else
                    return nPow(-a, n)
                end
            else
                local y = 1.0
                for i = 1, n do
                    y = y * a
                end
                return y
            end
        end
    elseif n == 0 then
        return 1.0
    else      --  when n < 0
        if a == 0.0 then
            error('Negative powering exception of zero.')
        else
            if n == -1 then
                return 1.0/a
            else
                return 1.0/nPow(a, -n)
            end
        end
    end
end


--
-- Compute the n-th root of x to a given scale, x > 0.
--
function gPow(a, n)
    if n > 0 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 % 2 == 1 then
                    return -1.0
                else
                    return 1.0
                end
            elseif a < 0.0 then
                if n % 2 == 1 then
                    return -gPow(-a, n)
                else
                    return gPow(-a, n)
                end
            else
                local y = 1.0
                local r = a
                local m = 8*4 - 1           --  8*sizeof(int) - 1;
                local one = 1
                local sn = itoa(n, 2)
                for i = 0, m-1 do
                    if string.sub(sn, string.len(sn) - i, string.len(sn) - i) == '0' then          
                        y = y * 1.0
                    else
                        y = y * r
                    end
                    r = r*r
                    one = one * 2
                    if one > n then
                        break
                    end
                end
                return y
            end
        end
    elseif n == 0 then
        return 1.0
    else      --  when n < 0
        if a == 0.0 then
            error('Negative powering exception of zero.')
        else
            if n == -1 then
                return 1.0/a
            else
                return 1.0/gPow(a, -n)
            end
        end
    end
end


--
-- Compute the n-th root of x to a given scale, x > 0.
--
function mPow(a, n)
    if n > 0 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 % 2 == 1 then
                    return -1.0
                else
                    return 1.0
                end
            elseif a < 0.0 then
                if n % 2 == 1 then
                    return -mPow(-a, n)
                else
                    return mPow(-a, n)
                end
            else
                y = 1.0
                r = a
                m = n
                while m > 0 do
                    if (m % 2) == 1 then
                        y = y * r
                    end
                    r = r*r
                    m = math.floor(m / 2)
                end
                return y
            end
        end
    elseif n == 0 then
        return 1.0
    else      --  when n < 0
        if a == 0.0 then
            error('Negative powering exception of zero.')
        else
            if n == -1 then
                return 1.0/a
            else
                return 1.0/mPow(a, -n)
            end
        end
    end
end


--
-- Compute the square root of x to a given scale, x > 0.
--
function heronSqrt(a)
    if a < 0.0 then
        error('Cannot find the sqrt of a negative number.')
    elseif a == 0.0 or a == 1.0 then
        return a
    else
        local x1 = a
        local x2 = (x1 + a/x1)/2.0
        local er = x1 - x2
        local counter = 0
        while x1 + er ~= x1 do
            x1 = x2
            x2 = (x1 + a/x1)/2.0
            er = x1 - x2
            if math.abs(er) < math.abs(M_EPSILON*x1) then
                break
            end
            counter = counter + 1
            if counter > MAX_ITER then
                break
            end
        end
        if counter >= MAX_ITER then
            error('Inaccurate sqrt exception by too many iterations.')
        end
        return x2
    end
end


--
-- Compute the cubic root of x to a given scale, x > 0.
--
function newtonCbrt(a)
    if a == 0.0 or a == 1.0 or a == -1.0 then
        return a
    elseif a < 0.0 then
        return -newtonCbrt(-a)
    else
        local x1 = a
        local x2 = (2.0*x1 + a/(x1*x1))/3.0
        local er = x1 - x2
        local counter = 0
        while x1 + er ~= x1 do
            x1 = x2
            x2 = (2.0*x1 + a/(x1*x1))/3.0
            er = x1 - x2
            if math.abs(er) < math.abs(M_EPSILON*x1) then
                break
            end
            counter = counter + 1
            if counter > MAX_ITER then
                break
            end
        end
        if counter >= MAX_ITER then
            error('Inaccurate cbrt exception by too many iterations.')
        end
        return x2
    end
end


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

 

local x = 16.0
local u = math.sqrt(x)


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

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

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

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

x = 12960000000000000000.0
local z = newtonNthRoot(4, x)
print(string.format("x = %g", x))
print(string.format("z = newtonNthRoot(4, x) = newtonNthRoot(4, %g) =  %g", x, z))
print(string.format("z*z*z*z = %g",  z*z*z*z))
print()

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


x = -4.0
print("[ Test Exception heronSqrt(double) ]--------------------")
print(string.format("x = %g", x))
print(string.format("Calculating heronSqrt(%g)", x))
local y  -- = heronSqrt(x)
status, value = pcall(heronSqrt, x)
if status then
    y = value
    print(string.format("y = heronSqrt(%g) = %g", x, y))
    print(string.format("y*y = %g", y*y))
    print()
else
    print(string.format("%s\nCaught some exception in calculating heronSqrt(%g)", value, x))
    print()
end

x = -4.0
print("[ Test Exception in newtonCbrt(double) ]--------------------")
print(string.format("x = %g", x))
print(string.format("Calculating newtonCbrt(%g)", x))
status, value = pcall(newtonCbrt, x)
if status then
    y = value
    print(string.format("y = newtonCbrt(%g) = %g", x, y))
    print(string.format("y*y*y = %g", y*y*y))
    print()
else
    print(string.format("%s\nCaught some exception in calculating newtonCbrtrt(%g)", value, x))
    print()
end


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

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

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


print(string.format("pow(2.1, 2.1) = %g", math.pow(2.1, 2.1)))
print(string.format("pow(2.1, -2.1) = %g", math.pow(2.1, -2.1)))
print(string.format("pow(2.1, 2.1) * pow(2.1, -2.1) = %g", math.pow(2.1, 2.1)*math.pow(2.1,  -2.1)))
print(string.format("2.1**2.1 = exp(2.1*log(2.1)) = %g", math.exp(2.1*math.log(2.1))))
print(string.format("2.1**(-2.1) = exp(-2.1*log(2.1)) = %g", math.exp(-2.1*math.log(2.1))))
print(string.format("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(string.format("math.pow(%g, %d) = %g", x, k,  math.pow(x, k)))   -- pow(x, k)
print(string.format("t1 = nPow(%g, %d) = %g", x, k,  t1))
print(string.format("t2 = gPow(%g, %d) = %g",x, k,  t2))
print(string.format("t3 = mPow(%g, %d) = %g", x, k,  t3))
print(string.format("t1 / t2 = %g", t1 / t2))
print(string.format("t1 - t2 = %g", t1 - t2))
io.write("t1 == t2 ? ")
if t1 == t2 then
    print "yes"
else
   print "no"
end
print(string.format("t1 / t3 = %g", t1 / t3))
print(string.format("t1 - t3 = %g", t1 - t3))
io.write("t1 == t3 ? ")
if t1 == t3 then
    print "yes"
else
   print "no"
end
print(string.format("t2 / t3 = %g", t2 / t3))
print(string.format("t2 - t3 = %g", t2 - t3))
io.write("t2 == t3 ? ")
if t2 == t3 then
    print "yes"
else
   print "no"
end
print()

print "Done."


--[[
# System.out.println("sizeof(double) = " + sizeof(double) );;
# System.out.println("sizeof(int) = " + sizeof(int) );
## print "sizeof(double) = %s = %d" % ("%X" % -1.0, len("%X" % -1.0))
## print "sizeof(int) = %s = %d" % ("%X" % -1, len("%X" % -1))
# print "sizeof(double)*4 = " + String.format("%X", Double.doubleToRawLongBits(-1.0)).length()*4 );;
# print sizeof(int)*4 = " + String.format("%X", -1).length()*4 );
# print sys.getsizeof(int)   # See: http://stackoverflow.com/questions/449560/how-do-i-determine-the-size-of-an-object-in-python
# print sys.getsizeof(-1.0)   # See: http://stackoverflow.com/questions/449560/how-do-i-determine-the-size-of-an-object-in-python
# print sys.getsizeof(-1.0)*4   # See: http://stackoverflow.com/questions/449560/how-do-i-determine-the-size-of-an-object-in-python
# print sys.getsizeof(-1.0)*4   # See: http://stackoverflow.com/questions/449560/how-do-i-determine-the-size-of-an-object-in-python
## print getsizeof(-1.0)*4   # See: http://stackoverflow.com/questions/449560/how-do-i-determine-the-size-of-an-object-in-python
--]]


--[[
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+011
exp(log(x)/3.0) = 9000
w = newtonCbrt(7.29e+011) = 9000
w*w*w = 7.29e+011

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

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

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

[ Test Exception heronSqrt(double) ]--------------------
x = -4
Calculating heronSqrt(-4)
testNthRoot.lua:228: 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

pow(2.1, 2.1) = 4.74964
pow(2.1, -2.1) = 0.210542
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-011
t1 == t2 ? no
t1 / t3 = 1
t1 - t3 = 6.18456e-011
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 )

 

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

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

아래의 소스는 대부분 버전의 Ruby 에서 실행되도록 작성된 소스이다.

참고로, Ruby 언어는 (다른 언어와 달리) 대분자만으로 쓰여진 것은 구문상 변수가 아니라 상수를 의미한다. 잠간 irb(인터랙티브 Ruby)를 실행하여 확인해 보자.

irb(main):001:0> ABC = 5
=> 5
irb(main):002:0> ABC = 3
(irb):2: warning: already initialized constant ABC
=> 3

그러므로 아래의 소스 첫 부분에

MAX_ITER = 20000
M_EPSILON = 1.0e-15

로 선언한 것은 (두 전역변수가 아니라) 두 상수를 선언한 것이며, 이는 소스 전체에서 통용된다.

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

# Filename: testNthRoot.rb
#
#            Approximate square roots, cubic roots and n-th roots of a given number.
#
# Execute: ruby testNthRoot.rb
#
# Date: 2013. 1. 6.
# Copyright (c) 2013 PH Kim  (pkim __AT__ scripts.pe.kr)

MAX_ITER = 20000
M_EPSILON = 1.0e-15

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


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

 

#
# Compute the n-th root of x to a given scale, x > 0.
#
def mPow(a, n)
    if n > 0
        if n == 1
            return a
        else
            if a == 0.0 or a == 1.0
                return a
            elsif a == -1.0
                if n % 2 == 1
                    return -1.0
                else
                    return 1.0
                end
            elsif a < 0.0
                if n % 2 == 1
                    return -mPow(-a, n)
                else
                    return mPow(-a, n)
                end
            else
                y = 1.0
                r = a
                m = n
                while m > 0 do
                    if (m & 0x1) == 1
                        y *= r
                    end
                    r = r*r
                    m >>= 1
                end
                return y
            end
        end
    elsif n == 0
        return 1.0
    else      #  when n < 0
        if a == 0.0
            raise RuntimeError,'Negative powering exception of zero.'
        else
            if n == -1
                return 1.0/a
            else
                return 1.0/mPow(a, -n)
            end
        end
    end
end


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


#
# 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
    elsif 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 do
            x1 = x2
            x2 = (2.0*x1 + a/(x1*x1))/3.0
            er = x1 - x2
            if (er).abs < (M_EPSILON*x1).abs
                break
            end
            counter += 1
            if counter > MAX_ITER
                break
            end
        end
        if counter >= MAX_ITER
            raise Error,'Inaccurate sqrt exception by too many iterations.'
        end
        return x2
    end
end


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


x = 16.0
u = Math.sqrt(x)

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

print "[ Testing newtonCbrt(double) ]--------------------\n"
x = -216.0
print "x = %g\n" % x
print "-exp(log(-x)/3.0) = %g\n" % -Math::exp(Math::log(-x)/3.0)
w = newtonCbrt(x)
print "w = newtonCbrt(%g) = %g\n" % [x, w]
print "w*w*w = %g\n" % (w*w*w)
print "\n"

x = 729000000000.0
print "x = %g\n" % x
print "exp(log(x)/3.0) = %g\n" % Math::exp(Math::log(x)/3.0)
w = newtonCbrt(x)
print "w = newtonCbrt(%g) = %g\n" % [x, w]
print "w*w*w = %g\n" % (w*w*w)
print "\n"

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

x = 12960000000000000000.0
z = newtonNthRoot(4, x)
print "x = %g\n" % x
print "z = newtonNthRoot(4, x) = newtonNthRoot(4, %g) =  %g\n" % [x, z]
print "z*z*z*z = %g\n" % (z*z*z*z)
print "\n"

x = 1.0/12960000000000000000.0
z = newtonNthRoot(4, x)
print "x = %g\n" % x
print "exp(log(x)/4.0) = %g\n" % Math::exp(Math::log(x)/4.0)
print "z = newtonNthRoot(4, x) = newtonNthRoot(4, %g) =  %g\n" % [x, z]
print "z*z*z*z = %g\n" % (z*z*z*z)
print "\n"


begin
    x = -4.0
    print "[ Test Exception heronSqrt(double) ]--------------------\n"
    print "x = %g\n" % x
    print "Calculating heronSqrt(%g)\n" % x
    y = heronSqrt(x)
    print "y = heronSqrt(%g) = %g\n" % [x, y]
    print "y*y = %g\n" % (y*y)
    print "\n"
rescue RuntimeError => ex
    print "%s\nCaught some exception in calculating heronSqrt(%g)\n" % [ex, x]
    print "\n"
end


begin
    x = -4.0
    print "[ Test Exception in newtonCbrt(double) ]--------------------\n"
    print "x = %g\n" % x
    print "Calculating newtonCbrt(%g)\n" % x
    y = newtonCbrt(x)
    print "y = newtonCbrt(%g) = %g\n" % [x, y]
    print "y*y*y = %g\n" % (y*y*y)
    print "\n"
rescue RuntimeError => ex
    print "%s\nCaught some exception in calculating newtonCbrtrt(%g)\n" % [ex, x]
    print "\n"
end


print "[ Test calculations by powering ]-----------------------------\n"
x = 200.0
z = newtonNthRoot(10, x)
print "x = %g\n" % x
print "exp(log(x)/10.0) = %g\n" % Math::exp(Math::log(x)/10.0)
print "z = newtonNthRoot(10, x) = newtonNthRoot(10, %g) = %g\n" % [x, z]
print "pow(z, 10) = %g\n" % z**10
print "\n"

x = 3001.0
z = newtonNthRoot(99, x)
print "x = %g\n" % x
print "exp(log(x)/99.0) = %g\n" % Math::exp(Math::log(x)/99.0)
print "z = newtonNthRoot(99, x) = newtonNthRoot(99, %g) = %g\n" % [x, z]
print "pow(z, 99) = %g\n" % z**99
print "\n"

x = 3001.0
z = newtonNthRoot(-99, x)
print "x = %g\n" % x
print "exp(log(x)/-99.0) = %g\n" % Math::exp(Math::log(x)/-99.0)
print "z = newtonNthRoot(-99, x) = newtonNthRoot(-99, %g) = %g\n" % [x, z]
print "1.0/pow(z, 99) = %g\n" % (1.0/z**99)
print "\n"

print "2.1**2.1 = pow(2.1, 2.1) = %g\n" % 2.1**2.1
print "2.1**(-2.1) = pow(2.1, -2.1) = %g\n" % 2.1**(-2.1)
print "2.1**2.1 * 2.1**(-2.1) = pow(2.1, 2.1) * pow(2.1, -2.1) = %g\n" % ((2.1** 2.1)*2.1**( -2.1))
print "2.1**2.1 = exp(2.1*log(2.1)) = %g\n" % Math::exp(2.1*Math::log(2.1))
print "2.1**(-2.1) = exp(-2.1*log(2.1)) = %g\n" % 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\n" % (Math::exp(2.1*Math::log(2.1)) * Math::exp(-2.1*Math::log(2.1)))
print "\n"


k = 301
x = -1.029
t1 = nPow(x, k)
t2 = gPow(x, k)
t3 = mPow(x, k)
print "math.pow(%g, %d) = %g\n" % [x, k,  x**k ]   # pow(x, k)
print "t1 = nPow(%g, %d) = %g\n" % [x, k,  t1]
print "t2 = gPow(%g, %d) = %g\n" % [x, k,  t2]
print "t3 = mPow(%g, %d) = %g\n" % [x, k,  t3]
print "t1 / t2 = %g\n" % (t1 / t2)
print "t1 - t2 = %g\n" % (t1 - t2)
print "t1 == t2 ? "
if t1 == t2
    print "yes\n"
else
   print "no\n"
end
print "t1 / t3 = %g\n" % (t1 / t3)
print "t1 - t3 = %g\n" % (t1 - t3)
print "t1 == t3 ? "
if t1 == t3
    print "yes\n"
else
   print "no\n"
end
print "t2 / t3 = %g\n" % (t2 / t3)
print "t2 - t3 = %g\n" % (t2 - t3)
print "t2 == t3 ? "
if t2 == t3
    print "yes\n"
else
   print "no\n"
end
print "\n"

print "Done.\n"


=begin
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.
=end

 

 

 

Posted by Scripter

댓글을 달아 주세요