음이 아닌 실수 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.
--]]
'소개와 설치 > Lua' 카테고리의 다른 글
| 최신 (윈도우 용) Lua 5.1.4 다운로드 및 설치 그리고 간단한 테스트 (0) | 2011.08.26 |
|---|
