n = 2 인 경우에는 Newton-Raphson 의 방법이 Heron 의 방법과 동일하다.
지수가 정수인 거듭제곱을 계산하는 함수도 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.
*/
?>
