Ackermann Function 값 구하기

옛날도 그렇고 지금도 그렇지만 저는 삽질성 프로그래밍을 좀 많이 하는 편입니다. 흐윽. 언제나 제대로 된 프로그램을 짤 수 있을 지....

아무튼 그런 삽질의 결과로 탄생한 몇 가지 프로그램 중, 가장 최근에 만든 것이 Ackermann Function 값 구하는 프로그램입니다. Ackermann Function은 (정보과학 전공하셨다면 다 아시겠지만) 이렇게 정의됩니다.

A(m,n)
= n+1                 if m=0
= A(m-1, 1)           if m!=0 and n=0
= A(m-1, A(m, n-1))   otherwise

직접 계산해 보시면 알겠지만 A(0, n)부터 A(3, n)까지는 별 무리 없이 계산이 가능합니다. 하지만 A(4, n)은 n=2일 때 이미 10^20000에 가까이 가 있고, A(5, n)은 n=0일 때 빼고는 전혀 계산이 불가능합니다.

오랜만에 위키백과사전을 뒤적거리다가 Ackermann Function 페이지에 다다랐을 때 머릿속에 갑자기 A(4,3)을 구해 보면 어떨까! 하는 생각이 들었습니다. 물론 A(4,3)은 적어도 10^19000자리가 넘는 무지막지하게 큰 수입니다만, 이 우주 안에서 표현이 불가능한 숫자를 전부 구하려는 바보같은 짓은 절대 안 하죠. 대신 저는 이 숫자의 마지막 몇 자리를 구해 보려 했습니다.

A(4,3) = 2^(2^65536) - 3입니다. 마지막 N자리를 구한다 치면 (2^(2^65536) - 3) mod 10^N을 구하는 거니까 결국 정수론 문제로 변해 버리고 마는데, 기초적인 정수론(사실 전 정말 기초적인 것 밖에 모릅니다 :)의 정리를 사용하면 어렵지 않게 구할 수 있습니다.

다음은 (매우 더러운-_-) 파이썬 코드입니다:

## .tokinized Python: Ackermann Function
## June 26, 2004. Kang Seonghoon (Tokigun; 'zenith%ctokigun%cnet' % (64,46))
## distributable under LGPL.

import itertools, math

def ackermann(m, n):
	"""Calculate the Ackermann Function A(m,n)."""
	# error handling
	if m < 0: raise OverflowError, "first argument must be >= 0."
	# generalized format (0 <= m <= 3)
	elif m == 0: return n + 1
	elif m == 1: return n + 2
	elif m == 2: return 2 * n + 3
	elif m == 3: return 2**(n+3) - 3
	# as definition (m > 3)
	elif n == 0: return ackermann(m - 1, 1)
	else: return ackermann(m - 1, ackermann(m, n - 1))

def _primetest(x):
	"""Iterator for phi function."""
	yield 2 # two and -
	yield 3 # three is not 6n+1/6n-1 prime.
	p = 5
	while p <= x:
		yield p # 6n-1 prime
		yield p + 2 # 6n+1 prime
		p += 6

def _sqrt(x):
	"""Returns greater value than sqrt(x) where x is positive integer, for factorization."""
##	return sqrt(x) # it can't handle long integer greater than 10^300
	y = str(x)
	if len(y) % 2 == 0:
		return 4 * long(y[:len(y)//2])
	else:
		return long(y[:len(y)//2 + 1])

def phi(x):
	"""Calculate Euler's totient function phi(x)."""
	# error handling
	if x <= 0: raise OverflowError, "phi(x) is defined for only positive integer x."
	
	totient = 1
	# for p which can be prime...
	for p in _primetest(_sqrt(x)):
		if x == 1: break
		# squeezing factors from x
		for n in itertools.count(0):
			if x % p == 0: x /= p
			else: break
		if n: totient *= (p-1) * p**(n-1)
	if x > 1: totient *= (x-1) # last factor
	return totient

def ackermann43(modulo, exponent10 = False):
	"""Calculate A(4,3) mod modulo using binary exponentiation. Thanks sanxiyn!"""
	terms = 7

	table = []
	if exponent10:
		# if modulo = 10^n, use phi^k(10^n) = 4^k * 10^(n-k).
		for x in xrange(0, terms + 1):
			table.insert(0, 4**x * modulo)
			modulo /= 10
	else:
		# else, use phi() function directly.
		for x in xrange(0, terms + 1):
			table.insert(0, modulo)
			modulo = phi(modulo)

	# calculate 2^(2^65536) mod modulo
	value = 1
	for x in xrange(0, terms):
		nvalue = value or table[x] # use p'=phi^k(10^n) if p=0.
		value = 1
		xvalue = 2 # 2^(2^k) mod modulo
		
		# binary exponentiation :)
		while nvalue:
			if nvalue % 2 == 1:
				value = (value * xvalue) % table[x+1]
			xvalue = (xvalue * xvalue) % table[x+1]
			nvalue = nvalue // 2

	# return (caculated value - 3).
	return (value - 3) % table[terms]

그리고 이걸로 ackermann43(10**2000, True)을 돌려 보면 (한참 후에) 다음과 같은 결과를 볼 수 있습니다. (첫째 argument는 A(4,3) mod m에서 m이고, 둘째 argument는 만약 m이 10의 거듭제곱 꼴일 경우 phi 계산을 생략하기 위해서 사용합니다.)

9192132667165851739267747048363889971566474914789067281426098995458694086764808807544415691412982362
9110560621500811980107332716373418540279342485096947298001269869994808709402939338274893473418136111
3646401608308716210604409754283027068790799721793634881311838403823265803485772391173918165833850930
7638773547842138398060195398672570286590870276078354922574731350253294242293585733901610724033405676
3233529824919500040840321441839262330478250767026489225446213936088996969372446928099642097643019211
0057975587360653163884791784727199840841778905347406677492320860307729730947221864787735071490594650
6001702669020411954599896851969532951089871664497072033893127451412555703792429465436153394490326612
5069117392015084367100890381248414333948851056327215992347214521240320673051342772087091544388305781
3979169328016186478173447140453174024691778124452031714615477067879000271075122364637024069641146831
8560596942867019060923922316502312982917708537460914981991319872431911819505667712358719823507473930
3207928664210514616925528169933437379078816780166602241399010598575782940409162406071882282742866488
4861265957725351502055606977843955657920709710656110708833854369329199976598310055799973996293977636
4453388722785985282815452157895371901370150488103364176421206789956233833382973419460918643809670075
3749403489246900618452027009588493693610707006376835247092315503180198793807275909812992686824617045
4054263099718005810182760277154054388182976416247682177558427667501992257454342390708885312186636048
3277277637118951766333897207632034303964686516279422202507409904914983002184459254460065477443775334
2503765846404126179638487650727694969306582301859277658341138177266102786709078054836212701794623555
7527699032203945992283624109828456932287631569981908293938478077944285174453688567242382204873780761
8793483615188054047198307281493106948514510174580999746807069993441788170742347723593009002851818204
9563524291540826682468708727402008345317538245636184782113173338559542427968873213010766659621748733

제 프로그램에서는 phi 함수와 binary exponentiation을 사용합니다만 상당히 비효율적으로 짜인 것 같습니다. 흑흑. 그 증거로 위의 2000자리 출력 결과는 제 컴퓨터에서 7분 걸렸습니다.

뭐 오랜만에 이렇게 삽질 해 보니까 재밌긴 한데, 혹시 이런 거 직접 해 보신 분 계신지요? :) 혹시나 저기서 더 고칠 데가 있을까요? (근데 이런 짓을 왜 하나 몰라... -_-;;)

음... 나중에 시간 나면 A(4,4)나 A(5,1) 같은 것도 한 번 해 봐야 겠습니다. :) (그런데 시험 기간이라서 상당히 압박스럽군요. 하루 투자해서 만든 건데... 쿨럭.)

- 토끼군

참고: http://tokigun.net/python.ackermann.php에 제가 이 프로그램에서 사용한 (안 좋은) 알고리즘에 대해 설명되어 있습니다.