음이 아닌 실수 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.
*/
'프로그래밍 > Groovy' 카테고리의 다른 글
| IntelliJ IDEA 를 이용하여 작성한 Groovy 애플리케이션 (1) | 2014.01.25 |
|---|---|
| 이진 파일을 읽어서 16진수로 보여주는 HexView 소스 with Groovy (0) | 2013.08.05 |
| Groovy 언어로 역삼각함수, 역쌍곡선함수 값을 구하는 예제 (0) | 2013.01.02 |
| Tomcat 6 또는 Tomcat 7에서 Groovlet을 reload할 때 404 에러가 나는 경우 (0) | 2012.11.20 |
| Groovy 언어의 지수 연산자 ** 의 연산 진행 순서에 관한 문제 (0) | 2010.07.04 |
