OCaml 언어로 평방근, 입방근, n제곱근 구하는 함수를 구현하고 테스트하기

음이 아닌 실수 A 의 평방근 sqrt(A) 를 구하는 Heron 의 방법:

        반복함수  g(x) = (x + A/x) / 2   를 이용

 

실수 A 의 n제곱근 root(n, A) 를 구하는 Newton-Raphson 의 방법

        반복함수  g(x) = ((n-1)*x + A/(x**(n - 1))) / n    를 이용

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

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

 

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