C++/Algorithms

C++/알고리듬에 대해 올리는 글입니다.

Visual C++ 2010/2013의 linker bug

Visual c++ 2010, 2013에는 “unresolved external symbol” linking bug가 있다 -_-;; 그렇다, 2013에서도 사라지지 않은 문제이다.

다음과 같은 세 가지 케이스가 아니라면 & 어쩔 땐 되고 어쩔 땐 안 된다면 linking bug를 의심해보자.

1. Template class의 method definition을 .h가 아니라 .cpp에 넣어두었다. <= template class의 method definition은 .h에 있어야 한다. Visual c++에선 definition이 cpp에 있어도 가끔씩 링킹이 된다 -_-;

2. signature가 안 맞는다. 특히 function name의 오타자 / scope / “&” / const 위치 주의

 

아무리 해봐도 linking error가 풀리지 않는다면.. 내가 찾은 야매를 공개하자면 :

1. 해당 method가 정의된 .cpp 파일을 오른쪽 클릭한 뒤 “컴파일” 해 본다.

2. 솔루션 clear 한 다음 rebuild / build 해본다.

3. 이건 본인이 아주 즐겨쓰는 방법인데 ㅡㅡ; method가 정의된 cpp 파일에다 “using namespace std;”를 추가한다.

링킹 에러가 사라졌다가 다시 생긴다면 “using namespace std;”를 주석처리해본다.(…)

 

3번 방법은 특히 아주 짜세이다..

C++0x의 shared_ptr

C++은 garbage collection을 지원을 하지 않기 때문에 생성한 heap 객체를 언제 free해 줘야 할 지가 중요한 issue가 된다.(그런데 정말로 C++가 garbage collection을 지원하지 않나..?) 이는 ‘잘 쓰면’ 문제가 없지면 문제는 어떤 객체의 수명이 언제까지인지, 그리고 누가 들고 있는지를 고려하면서 프로그래밍을 하는 것은 때때론 지나치게 어렵다는 것이다..

JAVA에서 어떤 객체가 garbage collect될 조건은 바로 그 객체를 가리키는 reference가 하나도 없을떄이다. 그리고 gc는 주기적으로 이것을 체크한다. C++0x에도 이와 똑같은 방식으로 heap memory를 계속 해재해 줄 수 있도록 std::shared_ptr<T>를 지원한다.(원래는 boost library에 있었다 한다)

1. 사용법

#include<memory>
#include<iostream>
using namespace std;

class A{
public:
	int x;
};

int main(){
	shared_ptr<int> p1(new int(4));
	shared_ptr<int> p2(p1); // reference count를 1 늘린다.
	shared_ptr<int> p3(new int(4));
	shared_ptr<A> p2(new A());

	(*p2).x = 123;
	cout << "p1 == p2 : " << (p1 == p2) << endl; // result : 1
	cout << "p1 == p3 : " << (p1 == p3) << endl; // result : 0 ; == is pointer-equivalence
	cout << "*p1 == *p3 : " << (*p1 == *p3) << endl; // result : 1
	cout << "p2->x : " << p2->x << endl;
	return 0;
}

원래 pointer T*가  있을 자리에 shared_ptr<T>를 다 넣기만 하면 된다 ㅇㅇ

STL container도 다 쓸 수 있는데(set, map, vector..) 이유는 comparator가 잘 정의가 되어 있기 때문이다. == 연산자는 ‘주소가 같은가?’를 비교하고 <, > 연산자는 주소값을 비교한다.(shared_ptr<T>간의 비교는 T의 comparator<를 부르지 않는다!!) (http://en.cppreference.com/w/cpp/memory/shared_ptr/operator_cmp)

위 코드에서 memory leak은 없다. 이유는 shared_ptr이 reference count를 체크하면서 자기가 빠졌을 때 count가 0인 애들에 대해 저절로 free를 다 해주기 때문이다.

2. 주의할 점

다음과 같은 코드는 망한다

#include<memory>
#include<iostream>
using namespace std;

int main(){
	int p = new int(4);
	{
		shared_ptr<int> p1(p);
		cout << *p1 << endl;
	}
	*p = 123;
	cout << *p << endl;
	return 0;
}

아래도 망한다. 이 때는 weak_ptr을 써야 한다고 하는데 어떻게 하는지는 찾아봐야 할 듯 하다

#include<memory>
using namespace std;

class Faulty{
public:
	shared_ptr<Faulty> toSharedPtr()
	{ return shared_ptr<Faulty>(this); }
};

3. 그 외

Thread-safe하다고 한다. 어떤 개체에 대해 여러 개의 shared_ptr<T>를 실시간으로 여러 개 생성하고 파괴해도 안전하다. 단, 인자를 &로 전달했을 경우엔 nonsafe할 수 있다고 한다.

shared_ptr<T> 간의 복사는 약간의 copy load가 있다고 한다.

shared_ptr<T>간의 static_cast, dynamic_cast, const_cast는 static_pointer_cast(), dynamic_pointer_cast(), const_pointer_cast() 함수를 불러서 해결할 수 있다.(http://en.cppreference.com/w/cpp/memory/shared_ptr/pointer_cast)

GDB 사용법/팁 모음

1. 기본적인 사용법

터미널에서 시작할 때

$ gdb 프로그램이름

한 뒤, 적당히 필요한 옵션들을 준다.

ex) b f 라고 한다면 함수 f를 시작할 때 breakpoint가 걸린다

그 다음 r 인자1 인자2 … 와 같이 실행하면 인자1 인자2 ..를 준 상태에서 프로그램 실행이 시작됨.

 

2. gdb로 STL 컨테이너들 보기

http://www.yolinux.com/TUTORIALS/GDB-Commands.html#STLDEREF

개인적으로 Wndows에서 mingw에 있던 mingw32-gdb.exe를 쓰는지라, .gdbinit 파일을 어디 두어야 하는지를 계속 찾다가 그냥 실행하는 폴더에 두니까 잘 되더라 –;

예를 들면, 내가 gdb로 소스들을 디버깅할 위치인 c:\users\aqjune\Desktop\workspace에다가 .gdbinit를 두고 gdb를 실행하니까 pvector 명령어가 잘 먹더라 이말씀

 

3. assert가 떴을 때 멈추게 하기

gdb에서는 자동으로 멈추게 하더라.

이게 assert가 떴을 때 gdb에서 break이 걸리지 않고 그냥 프로그램이 끝나는 경우가 있는데, 컴파일할 때 -g 옵션을 안 주면 안 멈추더라 이말씀

만약 -g 옵션을 안 주고 컴파일했을 때도 assert에서 프로그램이 멈추게 하고 싶으면, “b abort”라고 입력하면 된다.

 

4. Break 걸렸을 때, call stack에서 왔다갔다 하기

b를 걸거나 assert에서 걸리거나 해서 멈추었을 경우, info stack 이라고 치면 call stack이 쭉 나온다.

현재 관찰하고자 하는 call frame에 가고 싶을 때는, frame 14와 같이 입력하면 call stack에서 14번째 함수로 이동함.

여기서 info local 이라고 하면 지역변수 값이 쭉 나오고, p 변수 하면 변수 값이 나온다.

STL set/priority_queue의 속도 비교

STL을 조금 사용한 분이라면, STL의 set의 속도가 꽤나 느리다는 것을 알고 있을 것이다.

난 이 점을 고려해서 set 대신 priority_queue를 쓸 수 있는 상황이라면 priority_queue를 항상 써왔었는데, 아는 후배로부터 priority_queue가 오히려 set보다 느리다는 이야기를 듣게 되어 한번 실험을 해 보게 되었다.

실험 내용은 단순하다, set으로, priority_queue로 한번씩 정렬을 해 본 다음 거기에 걸리는 시간을 비교하는 것이다. 코드 내용은 아래와 같다.

#include<stdio.h>
#include<queue>
#include<set>
#include<stdlib.h>
#include<time.h>
using namespace std;
#define datatype pair<double, double>

set<datatype> s;
priority_queue<datatype> p;
vector<datatype> data;
datatype tmpdata;

void data_gen(){
	double dd = 1.0;
	int i;
	datatype d;
	for(i = 0; i < 1000000; i++){
		dd = (rand() % 1000000) * (rand() % 1000000);
		if(dd > 10.0 || dd < -10.0)
			dd /= 100;

		d.first=  dd;
		dd = 1.5 + dd * dd - 1.0 / dd;
		if(dd > 10.0 || dd < -10.0)
			dd /= 100;

		d.second = dd;

		data.push_back(d);
	}
}
void sort_pq(){
	printf("sort_q beg\n");
	int i;
	for(i = 0; i < data.size(); i++)
		p.push(data[i]);
	while(!p.empty())
	{
		tmpdata  = p.top();
		p.pop();
	}
}
void sort_set(){
	printf("sort_set beg\n");
	int i;
	for(i = 0; i < data.size(); i++)
		s.insert(data[i]);
	set<datatype>::const_iterator iend = s.end();
	for(set<datatype>::const_iterator itr = s.begin(); itr != iend; itr++){
		tmpdata  = *itr;
	}
	s.clear();
}
void run(void (*func)()){
	clock_t startt = clock();
	func();
	clock_t endt = clock();
	printf("time diff : %lf\n", (endt- startt) / (double)CLOCKS_PER_SEC);
}
int main(){
	data_gen();
	printf("datagen end..\n");

	run(sort_set);
	run(sort_pq);
	return 0;
}

정렬을 할 데이터 타입은 pair<double, double>로 했는데, 이유는 다음과 같다.

(1) set이나 priority_queue에서 pair간의 비교를 하려면 오버로드된 operator를 불러야 하고, double간의 대소 비교 또한 코스트가 크다. 이는 지나치게 비교를 많이 하는 data structure에 강한 패널티를 줄 것이다.
이것은 set에 약간 더 큰 패널티를 줄 것으로 예상하는데, 이유는 priority_queue가 Binary heap으로 되어 있을 경우 높이가 정확히 log_2 N 밖에 되지 않지만 set은 RBTree로 이루어 졌을 경우 height가 더 높을 것이기 때문이다.
(2) integer과는 달리 pair<double, double>은 만약 value copy가 일어날 경우 16바이트나 복사를 해야 하는데, 이는 적어도 여러 instruction을 거쳐야 이루어진다. 이는 값 복사가 너무 많이 일어나는 data structure에 강한 패널티를 줄 것이다.
이것은 priority_queue에 좀 더 큰 패널티를 줄 것으로 예측하는데, RBtree는 노드 삽입을 할 때만 value copy가 일어나고 order를 맞출 때는 전혀 value copy를 쓰지 않지만 Binary heap은 (priority_queue가 배열로 구현되어 있다고 가정할 시) heap 내 order를 맞출 때 여러 번의 value copy가 일어날 수 있기 때문이다.

 

set으로 sorting을 할 때는 iterator를처음부터 끝까지 돌면서 global 변수에 값복사를 한 번 하는 행동을 반복 하였다.

priority_queue로 sorting을 할 때는 top(), pop()을 번갈아가면서 부르면서 아까 말한 global 변수에 값복사를 한 번 하는 행동을 반복하였다.
둘 다 sorting이 끝났을 때 텅 비어있도록 하기 위해서 마지막에 set.clear()를 호출해 주었다.

g++, g++ -O2, VC++2010 Debug, VC++2010 Release에서 각각 돌려보았는데, 결과는 상당히 재미있었다 -_-;

 

1. G++

datagen end..
sort_set beg
time diff : 2.788000
sort_q beg
time diff : 3.670000

보다시피 priority_queue가 더 느렸다.

2. G++ -O2

datagen end..
sort_set beg
time diff : 1.628000
sort_q beg
time diff : 1.283000

둘 다 빨라졌으나 set이 약간 더 느렸다.

3. MSVC++2010 Debug

datagen end..
sort_set beg
time diff : 26.297000
sort_q beg
^C

priority_queue를 10여분동안 켜두었으나 무려 끝나지가 않았다..

4. MSVC++ 2010 Release

datagen end..
sort_set beg
time diff : 1.721000
sort_q beg
time diff : 0.585000

priority_queue가 3배 가량 빨랐다.

 

요약을 하자면, 컴파일러가 optimize를 하지 않았을 땐 set이 더 빠르나 optimize를 했을 땐 priority_queue가 더 빠른 속도를 보였다.

하긴 priority_queue는 set보다 할 수 있는게 더 적은데 얘가 더 느리면 안되지 -_-;;

여하튼 독특한 경험이었다.

Chu-Liu/Edmonds’ algorithm in C++ / POJ3164

Directed graph에서 MST를 구하는 알고리즘인 Chu-Liu / Edmonds’ algorithm입니다. POJ 3164(http://poj.org/problem?id=3164)를 푸는데 필요한 핵심 알고리즘이기도 합니다.

http://en.wikipedia.org/wiki/Edmonds’_algorithm 의 설명을 참조하여 만들었습니다.

알고리즘의 procedure은 간단히 설명하자면 다음과 같습니다.

(1) 주어진 graph G(V, E)에서, 어떤 노드 v로 들어오는 edge 중 가장 cost가 작은 edge를 (\pi(v), v)로 하고, 이 edge들로 이루어진 새로운 그래프를 G'(V, E’)라 합니다.

(2) 만약 E’에 cycle이 없다면 만들어진 graph는 SRT(shortest route tree)를 나타내며 이는 G의 MST가 됩니다.

(3) 만약 cycle이 존재한다면, 이 cycle을 하나의 node로 한 새로운 graph G”에 대해서 (1)을 다시 수행합니다. 단, 이 과정에서 edge들의 cost가 바뀌어야 되는 경우가 존재하는데, 이곳에 적기엔 상당히 복잡하고 위키피디아 문서에 자세히 설명이 되어 있습니다.

// Time complexity O(N^2)
// Assumes node 0 is the root node.
// i_map[N][N] : input map,
// SRT_pi[v] : output. (MST에서, 노드 v의 parent.)
// node 0이 MST의 root라고 가정함.
// Accepted : POJ 3164
#define MAX_NODE 100
#define NO_EDGE numeric_limits<double>::infinity()

double edmonds_map[MAX_NODE * 3][MAX_NODE * 3]; // input.
int org_N; // input에서 주어진 노드 수

bool node_isused[MAX_NODE * 3];
int cur_N; // node 수가 늘어나면서 된 노드 수
int SRT_pi[MAX_NODE * 3]; // SRT_pi[v] : pi[v] ; v의 parent. 답이 저장되는 곳.

bool dfs_visited[MAX_NODE * 3];
bool dfs_total_visited[MAX_NODE * 3];
int dfs_history[MAX_NODE * 3];
int dfs_history_size;

bool dfs_pi_findCycle(int node){
	int i;

	if(SRT_pi[node] != -1){
		i = SRT_pi[node];
		assert(edmonds_map[SRT_pi[node]][node] != NO_EDGE);
		assert(node_isused[SRT_pi[node]]);

		if(dfs_visited[i]){
			dfs_history[dfs_history_size++] = i;
			// cycle found!!
			return true;
		}
		else if(!dfs_total_visited[i]){
			dfs_total_visited[i] = true;
			dfs_visited[i] = true;
			dfs_history[dfs_history_size++] = i;
			if(dfs_pi_findCycle(i))
				return true;
			dfs_visited[i] = false;
			dfs_history_size--;
		}
	}
	return false;
}

void calc_SRT_pi()
{
	int i;
	SRT_pi[0] = -1;
	for(i = 1; i < cur_N; i++){
		if(!node_isused[i])
			continue;
		SRT_pi[i] = -1;

		int j;
		for(j = 0; j < cur_N; j++){
			if(!node_isused[j])
				continue;
			if(edmonds_map[j][i] == NO_EDGE)
				continue;
			if(SRT_pi[i] == -1)
				SRT_pi[i] = j;
			else if(edmonds_map[SRT_pi[i]][i] > edmonds_map[j][i])
				SRT_pi[i] = j;
		}
	}
}

void edmonds_algorithm(){
	int i;
	assert(cur_N >= org_N);
	assert(cur_N < org_N * 3);

	calc_SRT_pi();

	bool cycle_found = false;
	for(i = 0; i < cur_N; i++)
		dfs_total_visited[i] = dfs_visited[i] = false;

	for(i = 0; i < cur_N; i++){
		if(!node_isused[i])
			continue;
		else if(dfs_total_visited[i])
			continue;

		dfs_total_visited[i] = dfs_visited[i] = true;
		dfs_history[0] = i;
		dfs_history_size = 1;
		if(dfs_pi_findCycle(i)){
			cycle_found = true;
			break;
		}
		dfs_visited[i] = false;
	}

	if(cycle_found){
		int cycle_start = -1;
		for(i = 0; i < dfs_history_size - 1; i++){
			if(dfs_history[i] == dfs_history[dfs_history_size - 1]){
				cycle_start = i;
				break;
			}
		}
		assert(cycle_start != -1);

		for(i = cycle_start; i < dfs_history_size - 1; i++){
			assert(node_isused[dfs_history[i]]);
			assert(dfs_history[i] != 0);// not root.
			node_isused[dfs_history[i]] = false;
		}
		int new_node = cur_N;

		int j;
		for(j = 0; j < cur_N; j++){
			if(!node_isused[j])
				continue;

			int k;
			for(k = cycle_start; k < dfs_history_size - 1; k++){
				// j -> dfs_history[k] :
				int v = dfs_history[k];
				if(edmonds_map[j][v] != NO_EDGE){
					double cost= edmonds_map[j][v] - edmonds_map[SRT_pi[v]][v];
					if(edmonds_map[j][new_node] == NO_EDGE)
						edmonds_map[j][new_node] = cost;
					else if(edmonds_map[j][new_node] > cost)
						edmonds_map[j][new_node] = cost;
				}
				if(edmonds_map[v][j] != NO_EDGE){
					if(edmonds_map[new_node][j] == NO_EDGE)
						edmonds_map[new_node][j] = edmonds_map[v][j];
					else if(edmonds_map[new_node][j] > edmonds_map[v][j])
						edmonds_map[new_node][j] = edmonds_map[v][j];
				}
			}
		}

		int cycle_nodes[MAX_NODE * 3];
		int SRT_pi_copied[MAX_NODE * 3];
		int cycle_nodes_size =0 ;
		for(j = cycle_start; j < dfs_history_size - 1; j++)
			cycle_nodes[cycle_nodes_size++] = dfs_history[j];
		for(j = 0; j < cur_N; j++)
			SRT_pi_copied[j] = SRT_pi[j];

		node_isused[new_node] = true;
		cur_N++;

		edmonds_algorithm();
		// 이제 SRT_pi[new_node], SRT_pi[v] = new_node인 v에 대해서 처리.

		cur_N--;
		node_isused[new_node] = false;

		// SRT_pi[new_node] -> new_node.
		// SRT_pi[new_node]는 cycle 밖에 있음.
		assert(node_isused[SRT_pi[new_node]]);
		for(j = 0; j < cycle_nodes_size; j++){
			if(edmonds_map[SRT_pi[new_node]][cycle_nodes[j]] - edmonds_map[SRT_pi_copied[cycle_nodes[j]]][cycle_nodes[j]]
					== edmonds_map[SRT_pi[new_node]][new_node]){
				SRT_pi[cycle_nodes[j]] = SRT_pi[new_node];
				break;
			}
		}
		int vv = -1 ;// SRT_pi[vv] == new_node인 vv ; new_node -> vv를 cycle에 있는 node -> vv로 바꾸어야 함.
		for(i = 0; i < cur_N; i++){
			if(!node_isused[i])
				continue;
			if(SRT_pi[i] == new_node){
				vv = i;
				for(j = 0; j < cycle_nodes_size; j++){
					if(edmonds_map[new_node][vv] == edmonds_map[cycle_nodes[j]][vv]){
						SRT_pi[vv] = cycle_nodes[j];
						break;
					}
				}
				assert(SRT_pi[vv] != new_node);
			}
		}
		for(j = 0; j < cycle_nodes_size; j++){
			assert(node_isused[cycle_nodes[j]] == false);
			node_isused[cycle_nodes[j]] = true;
		}
		node_isused[cur_N] = false;
	}
}

void init_edmonds_algorithm(){
	cur_N = org_N = N;
	int i, j;
	for(i = 0; i < N; i++)
		for(j = 0; j < N; j++)
			edmonds_map[i][j] = i_map[i][j];
	for(i = 0; i < 3 * N; i++)
		for(j = 0; j< 3 * N; j++)
			if(!(i < N && j < N))
				edmonds_map[i][j] = NO_EDGE;
	for(i = 0; i < N; i++)
		node_isused[i] = true;
	for(i = N; i < N * 3; i++)
		node_isused[i] = false;
}

http://poj.org/problem?id=3164

Trie 헤더 파일

연구참여 하면서 VIsual c++로 간단히 작성했던 Trie. 묵혀두기 아까워서 공유 -_-aㅋㅋㅋ

string을 key로 하는 dictionary type을 만들고 싶을 때, STL의 map보다 월등히 높은 성능(특히 메모리)을 나타냅니다.

iterator, get/put/erase 모두 지원합니다 :)

Memory leak이 없는 것은 확인 하였고, Class를 value로 하는 타입에 대해서는 프로그램이 제대로 동작하지 않을 수도 있습니다 -,.-

Key로 한국어도 가능합니다.

 

Trie.h :

/*
author : akdal
https://agidari.wordpress.com
*/
#ifndef TRIE_H

#define TRIE_H

#include<vector>
#include<assert.h>
#include<algorithm>
#include<iostream>
using namespace std;

template<class T>
class Trie{
public:
        //static int NODE_ALLOCCNT;
        //static int NODE_DEALLOCCNT;
protected:
        class Node{
        public:
                T value;
                bool hasValue;
                vector<pair<char, Node *> > *children;
                bool isLeaf;
                bool isConstant;

                Node(T &value){
                        this->hasValue = true;
                        this->value = value;
                        this->isLeaf = true;
                        this->chlidren = NULL;
                        this->isConstant = false;
                        //Trie<T>::NODE_ALLOCCNT++;
                }
                Node(){
                        this->hasValue = false;
                        this->isLeaf = true;
                        this->children = NULL;
                        this->isConstant = false;
                        //Trie<T>::NODE_ALLOCCNT++;
                }
                Node(const Node &othernode){
                        this->value = othernode.value;
                        this->hasValue = othernode.hasValue;
                        this->isLeaf = othernode.isLeaf;
                        this->isConstant = othernode.isConstant;
                        this->children = NULL;
                        if(othernode.children != NULL){
                                this->children = new vector<pair<char, Node *> >();
                                for(vector<pair<char, Node *> >::iterator itr = othernode.children->begin(); itr != othernode.children->end(); itr++){
                                        pair<char, Node *> &p = *itr;
                                        pair<char, Node *> newp;
                                        newp.first = p.first;
                                        newp.second = new Node(*p.second);
                                        this->children->push_back(newp);
                                }
                        }
                        //Trie<T>::NODE_ALLOCCNT++;
                }
                ~Node(){
                        int i;
                        if(!isLeaf){
                                for(i = 0; i < children->size(); i++)
                                        delete (*children)[i].second;
                                delete children;
                        }
                        //Trie<T>::NODE_DEALLOCCNT++;
                }
                void makeConstant(){
                        this->isConstant = true;
                        if(isLeaf)
                                return;
                        int i, j;
                        for(i = 0; i < this->children->size(); i++){
                                for(j = i + 1 j < this->children->size(); j++){
                                        if(children->at(i).first > children->at(j).first)
                                                swap(children->at(i), children->at(j));
                                }
                        }
                        for(i = 0; i < this->children->size(); i++)
                                this->children->at(i).second.makeConstant();
                }
                Node *getNode(char c){
                        if(this->isLeaf)
                                return NULL;

                        if(this->isConstant){
                                int pos;
                                int sz = children->size();
                                int st= 0, ed = sz - 1;
                                while(true){
                                        if(st == ed){
                                                pair<char, Node *> &p = children->at(st);
                                                if(p.first == c)
                                                        return p.second;
                                                else
                                                        return NULL;
                                        }else if(st == ed - 1){
                                                pair<char, Node *> &p = children->at(st);
                                                if(p.first == c)
                                                        return p.second;
                                                else
                                                        st = ed;
                                        }else{
                                                int piv = (st + ed) / 2;
                                                pair<char, Node *> &p = children->at(piv);
                                                if(p.first == c)
                                                        return p.second;
                                                else if(p.first < c){
                                                        ed= piv - 1;
                                                }else{
                                                        st = piv + 1;
                                                }
                                        }
                                }
                                return NULL;
                        }else{
                                vector<pair<char, Node*> >::iterator iend = children->end();
                                for(vector<pair<char, Node *> >::iterator itr = children->begin(); itr != iend; itr++){
                                        pair<char, Node *> &p = *itr;
                                        if(p.first == c)
                                                return p.second;
                                }
                                return NULL;
                        }
                }
                void deleteNode(char c){
                        assert(!isConstant);
                        if(this->isLeaf){
                                assert(false);
                                return;
                        }
                        vector<pair<char, Node *> >::iterator iend =children->end();
                        for(vector<pair<char, Node *> >::iterator itr = children->begin(); itr != iend; itr++){
                                pair<char, Node *> &p = *itr;
                                if(p.first == c){
                                        pair<char, Node *> p = *itr;
                                        children->erase(itr);
                                        assert(p.second != NULL);
                                        delete p.second;
                                        break;
                                }
                        }
                        if(this->children->size() == 0)
                                this->isLeaf = true;
                }
                Node *childAt(int i)
                { return this->children->at(i).second; }
                char edgeAt(int i)
                { return this->children->at(i).first; }
                int size()
                { return this->children->size(); }
                Node *addNode(char c){
                        assert(!isConstant);
                        if(this->isLeaf){
                                this->isLeaf = false;
                                this->children = new vector<pair<char, Node *> >();
                        }
                        Node *n = new Node();
                        this->children->push_back(pair<char, Node *>(c, n));
                        return n;
                }
                void print(ostream &ostr, int space){
                        int i;
                        if(this->hasValue){
                                for(i = 0; i < space; i++)
                                        ostr << (" ");
                                ostr << "VALUE='" << this->value << "'" << endl;
                        }
                        if(!this->isLeaf)
                                for(i = 0; i < this->children->size(); i++){
                                        int j;
                                        for(j = 0; j < space; j++)
                                                ostr << (" ");
                                        ostr << "edge '" << this->children->at(i).first << "'" << endl;
                                        this->children->at(i).second->print(ostr, space + 1);
                                }
                }
        };

public:
        class iterator{
        private:
                char *trace;
                int *stack;
                Node **nodeStack;
                //int stackSize;
                int nodeStackSize;
                int maxStackSize;
                bool end;
                Node *root;

        protected:

                void _nextItr(){
                        if(end)
                                return;
                        if(this->nodeStackSize == 0){
                                this->nodeStackSize = 1;
                                this->nodeStack[0] = root;
                                this->trace[0] = 0;
                                return;
                        }
                        if(!this->nodeStack[this->nodeStackSize - 1]->isLeaf){
                                // ?????!
                                this->nodeStack[this->nodeStackSize] = this->nodeStack[this->nodeStackSize - 1]->childAt(0);
                                this->trace[this->nodeStackSize - 1] = this->nodeStack[this->nodeStackSize - 1]->edgeAt(0);
                                this->trace[this->nodeStackSize] = 0;
                                this->stack[this->nodeStackSize -1] = 0;
                                this->nodeStackSize++;
                        }else{
                                // ?????!
                                while(true){
                                        this->nodeStackSize--;
                                        if(this->nodeStackSize == 0){
                                                end = true;
                                                break;
                                        }
                                        // this->stack[this->nodeStackSize - 1]?? ????. this->nodeStackSize[this->nodeStackSize - 1]->children size??
                                        // ????? ???? fold.
                                        this->stack[this->nodeStackSize - 1]++;
                                        int idx = this->stack[this->nodeStackSize - 1];
                                        if(idx < this->nodeStack[this->nodeStackSize - 1]->children->size()){
                                                this->nodeStack[this->nodeStackSize] = this->nodeStack[this->nodeStackSize - 1]
                                                                ->childAt(idx);
                                                this->trace[this->nodeStackSize - 1] = this->nodeStack[this->nodeStackSize - 1]->edgeAt(idx);
                                                this->trace[this->nodeStackSize] = 0;
                                                this->nodeStackSize++;
                                                break;
                                        }
                                }
                        }
                }
        public:
                iterator(int maxStackSize, Node *root){
                        this->maxStackSize = maxStackSize;
                        this->nodeStackSize = 0;
                        //this->stackSize = 0;

                        this->trace = new char[this->maxStackSize + 1];
                        this->stack = new int[this->maxStackSize];
                        this->nodeStack = new Trie<T>::Node *[this->maxStackSize + 1];
                        this->end = false;
                        this->root = root;
                }
                iterator(const iterator &other){
                        this->maxStackSize = other.maxStackSize;
                        this->nodeStackSize = other.nodeStackSize;
                        this->root = other.root;
                        this->end = other.end;

                        this->trace = new char[this->maxStackSize + 1];
                        memcpy(this->trace, other.trace, this->maxStackSize + 1);;
                        this->stack = new int[this->maxStackSize];
                        memcpy(this->stack, other.stack, this->maxStackSize * sizeof(int));
                        this->nodeStack = new Node*[this->maxStackSize + 1];
                        memcpy(this->nodeStack, other.nodeStack, (this->maxStackSize + 1) * sizeof(int));
                }

                iterator& operator=(const iterator &other){
                        if(this->nodeStack)
                                delete this->nodeStack;
                        if(this->stack)
                                delete this->stack;
                        if(this->trace)
                                delete this->trace;

                        this->maxStackSize = other.maxStackSize;
                        this->nodeStackSize = other.nodeStackSize;
                        this->root = other.root;
                        this->end = other.end;

                        this->trace = new char[this->maxStackSize + 1];
                        memcpy(this->trace, other.trace, this->maxStackSize + 1);;
                        this->stack = new int[this->maxStackSize];
                        memcpy(this->stack, other.stack, this->maxStackSize * sizeof(int));
                        this->nodeStack = new Node*[this->maxStackSize + 1];
                        memcpy(this->nodeStack, other.nodeStack, (this->maxStackSize + 1) * sizeof(int));

                        return *this;
                }
                iterator(iterator &&other){
                        this->trace = other.trace;
                        this->stack = other.stack;
                        this->nodeStack = other.nodeStack;
                        other.trace = NULL;
                        other.stack = NULL;
                        other.nodeStack = NULL;

                        this->maxStackSize = other.maxStackSize;
                        this->nodeStackSize = other.nodeStackSize;
                        this->root = other.root;
                        this->end = other.end;
                }
                ~iterator(){
                        if(trace)
                                delete trace;
                        if(stack)
                                delete stack;
                        if(nodeStack)
                                delete nodeStack;
                }

                bool isEnd()
                { return this->end; }
                T *next(){
                        if(end)
                                return NULL;

                        do{
                                _nextItr();
                        }while(!end && !this->nodeStack[this->nodeStackSize - 1]->hasValue);

                        T *answer;
                        if(this->nodeStackSize == 0)
                                answer = NULL;
                        else
                                answer=   &this->nodeStack[this->nodeStackSize - 1]->value;

                        return answer;
                }
                const char *getCurrentTrace()
                { return this->trace; }
        };

private:
        Node *rootNode;
        int maxDepth;
        int elemSize;
        bool doConstant;
public:
        Trie(){
                this->rootNode = new Node();
                this->maxDepth = 1;
                this->doConstant = false;
                this->elemSize= 0;
        }
        Trie(const Trie &orgtrie){
                this->rootNode = new Node(*(orgtrie.rootNode));
                this->maxDepth = orgtrie.maxDepth;
                this->elemSize = orgtrie.elemSize;
                this->doConstant = false;
        }
        void makeConstant()
        { this->doConstant = true; }
        void put(const char *letters, const T& value);
        void put(const string &str, const T& value);
        void erase(const string &str){
                assert(!doConstant);
                this->_erase(str.c_str(), this->rootNode);
                this->elemSize--;
        }
        T *get(const char *letters);
        T *get(const string &str);
        typename Trie<T>::iterator itr();
        void clear();
        int size() const
        { return this->elemSize; }

        void print(ostream &st)
        {this->rootNode->print(st, 0); }

protected:
        void _put(const char *letters, Node *focus, const T& t);
        // returns : node???? *chrs?? ????? ???? ???? node ????? ?????? ???°? ???°??
        bool _erase(const char *chrs, Node *node){
                if(*chrs == 0){
                        node->hasValue = false;
                        return node->isLeaf;
                }
                Trie<T>::Node *n = node->getNode(*chrs);
                assert(n != NULL);
                if(this->_erase(chrs + 1, n) == true){
                        node->deleteNode(*chrs);
                }
                if(node->isLeaf && !node->hasValue)
                        return true;
                return false;
        }
        T *_get(const char *letters, Node *focus);
        T *_get(const string &str, int idx, int len, Node *focus);
};

#include"Trie.h"
#include<assert.h>

template<class T>
void Trie<T>::_put(const char *letters, Node *focus, const T& t){
        if(letters[0] == 0){
                focus->hasValue = true;
                focus->value = t;
        }else{
                Trie<T>::Node *n = focus->getNode(letters[0]);
                if(n == NULL)
                        this->_put(letters + 1, focus->addNode(letters[0]), t);
                else
                        this->_put(letters + 1, n, t);
        }
}
template<class T>
void Trie<T>::put(const char *letters, const T &t)
{
        assert(!this->doConstant);
        int len = strlen(letters) + 1;
        this->elemSize++;
        if(len > this->maxDepth)
                this->maxDepth = len;
        this->_put(letters, this->rootNode, t);
}
template<class T>
void Trie<T>::put(const string &str, const T& value){
        assert(!this->doConstant);
        this->put(str.c_str(), value);
}

template<class T>
T* Trie<T>::_get(const char *letters, Node *node){
        if(letters[0] == 0)
                return node->hasValue ? &(node->value) : NULL;
        Trie<T>::Node *n = node->getNode(letters[0]);
        if(n == NULL)
                return NULL;
        return this->_get(letters + 1, n);
}
template<class T>
T* Trie<T>::_get(const string &str, int idx, int len, Node *node){
        if(idx == len)
                return node->hasValue ? &(node->value) : NULL;
        Trie<T>::Node *n = node->getNode(str.at(idx));
        if(n == NULL)
                return NULL;
        return this->_get(str, idx + 1, len, n);
}

template<class T>
T *Trie<T>::get(const char *letters)
{ return this->_get(letters, this->rootNode); }
template<class T>
T *Trie<T>::get(const string &letters)
{ return this->_get(letters, 0, letters.length(), this->rootNode); }
/*
template<class T>
void Trie<T>::erase(const string &str){
        this->_erase(s.c_str(), this->rootNode);
}
*/
template<class T>
typename Trie<T>::iterator Trie<T>::itr(){
        return Trie<T>::iterator(this->maxDepth, tihis->rootNode);
}

template<class T>
void Trie<T>::clear(){
        assert(!this->doConstant);
        delete this->rootNode;
        this->maxDepth = 0;
        this->rootNode = new Node();
        this->elemSize = 0;
}

#endif

Trie_test.cpp :

/*
 * Author : akdal
 * https://agidari.wordpress.com
 * */
#include"Trie.h"
#include<iostream>
using namespace std;

int main(){
	Trie<int> t;
	Trie<bool> tt;
	t.put("aaaaa", 5);
	t.put("a", 1);
	t.put("ab", 2);
	t.put("bc", 2);
	t.put("b", 1);
	t.put("테스트", 4); // Korean is also possible.
	t.put("테스블", 5);

	int *p_aaa = t.get("aaa");
	int *p_a = t.get("a");
	int *p_dd = t.get("dd");
	int *p_bc = t.get("bc");
	int *p_ab = t.get("ab");
	int *p_b = t.get("b");

	t.print(cout);
	Trie<int>::iterator itr = t.itr();
	while(!itr.isEnd()){
		int *val = itr.next();
		if(val == NULL)
			break;
		cout << itr.getCurrentTrace() << " : " << *val << endl;
	}
	t.erase("bc");
	t.print(cout);
	cout<<"--"<<endl;
	t.erase("a");
	t.erase("aaaaa");
	t.print(cout);
	cout<<"--"<<endl;
	return 0;
}

Explicit method와 Implicit method

numerical analysis의 한 영역인 differential equation solving에 대한 explicit method와 implicit method 두 방식에 대한 간단한 글입니다.

시간이 난다면, 실험과 함께 포스팅을 이어가도록 하겠습니다.

특정한 differential equation을 풀기 위한 방법은 수없이 많다.(자기가 지어내면 그것이 또 다른 한 방법 ^^) 그 많은 방법을 분류할 수 있는 기준 중에 강력한 하나가 바로 explicit 형식이냐, 아니면 implicit형식이냐 이다.

풀고자 하는 1계 differential equation을

f(t, y(t)) = dy/dt (t)

이라고 할 때, explicit method는 dt시간 후의 y값을 구하기 위해

y(t + dt) = y(t) + dy/dt(t) dt = y(t) + f(t, y(t))dt

에서 y(t + dt)를 구하고, implicit method는

y(t + dt) = y(t) + dy/dt(t + dt)dt = y(t) + f(t + dt, y(t + dt))dt

에서 y(t + dt)를 구한다.

보다시피, implicit method의 경우 좌변의 argument로 아직 값이 ‘결정되지 않은’ y(t + dt)꼴이 들어가기 때문에, 이를 newton method 혹은 더욱 간단한 근사를 사용해 나타내기도 한다.

여기에서, implicit method의 경우에는 근사를 explicit method보다 ‘어쩔 수 없이’ 한 번 더 사용할 수밖에 없구나! 라는 것을 알게 된다. 그렇다면 그 결과도 implicit method가 explicit method보다 더 못하리라는 예측을 할 수 있을것이다.

하지만 실상은 정 반대이다.

많은 예제에서 explicit method를 통해 구한 해는 발산한다. 가장 간단한 예로는 용수철상수가 매우 큰 용수철의 시뮬레이션이 있다. dt를 어느 정도 작게 하지 않는 한, |y(t + dt)| > |y(t)|가 되버린다. 이 외에도 전혀 ‘발산하리라고 생각지 않았던’ 부분에서, explicit method는 시간이 가면서 아주 빠르게 발산하는 것을 볼 수 있다. 값이 발산해 버리는 것은 시뮬레이터의 입장으로써는 바람직하지 않을 것이다.

하지만 implicit method에서는 반대로 수렴한다. 그 이유는 다음과 같다 : 이론적으로, implicit method는 시간을 거꾸로 돌리는 explicit method와 동일하기 때문이다. (물론 y(t + dt)를 구하면서 근사의 방법을 쓰기에 time-reversed explicit method와 완전히 동일하지는 않겠지만) 시간이 0으로 갈수록 발산하기에, 시간이 infinite로 갈수록 수렴하게 된다.(실제 implicit method로 시뮬레이션을 step by step으로 보면 아마 납득이 갈텐데.. 참고로 정말 신기하다)

implicit method로 구한 답은 (정답에 가까워지든 멀어지든) 어느 방향으로든 수렴하기 때문에 시뮬레이터로써는 안정성을 갖출 수 있다.

여기까지가 어느 정도 정리한 글이다.

그럼 여기서 질문 :

‘답에 멀어지든 답으로 수렴하든, explicit method로 풀었을 때 어떻게든 수렴하는 미분방정식은 없을까?’

‘그래서, implicit method로 풀었을 경우 발산해버리는 미분방정식은 없을까?’