BigFloat로 Pi를 구해보자-3. 곱셈 구현 + 10진수로 출력하기

저번 글에서는 간단하게 덧셈과 곱셈을 구현해보았는데요, 이번에는 곱셈을 구현해보고 합니다. 이번에도 곱셈을 처음배우는 아이에게 곱셈을 가르쳐주듯, 컴퓨터에게 차근차근 가르쳐봅시다.

a = a_0 . a_1 a_2

b = b_0 . b_1 b_2 라고 했을때 이 두 수의 곱셈은 다음과 같이 계산됩니다.

예전에 곱셈공부할때가 새록새록 생각나지요??

c = ab라고 하면,

c_0 = a_0 b_0

c_1 = a_0 b_1 + a_1 b_0

c_2 = a_0 b_2 + a_1 b_1 + a_2 b_0

c_3 = a_1 b_2 + a_2 b_1

c_4 = a_2 b_2

가 됩니다. (편의상 각 자리에서 일어나는 받아올림은 생략했습니다... 만 물론 당연히 이것도 계산해주어야겠지요.)

n자리 짜리 두 수를 곱하는데는 총 n*n번의 곱셈이 필요하다는 것과, 곱셈 결과는 자릿수가 늘어난다는것 (거의 2배로 늘어나지요...) 을 알수 있었습니다. 그럼 알고리즘을 어떻게 세울지 구상을 해봅시다.

1. 곱셈결과를 저장할 c을 만들고 깨끗하게 0으로 세팅합니다.

2. 모든 i,j쌍에 대해 a_i 와 b_j를 가져와서 곱하고, c_(i+j)자리에 더해줍니다.

3. c를 출력하면 끝

자 이제 코드를 작성해볼까요?

/* RRU32는 부호없는 32비트 정수, RRU64는 부호없는 64비트 정수입니다. UINT32, UINT64로 했더니 다른 플랫폼에서는 막 충돌하더군요. 흔한 이름으 폐해...ㅠㅠ */ inline int i_addu32(RRU32* dest, RRU32 len, RRU32 src, RRU32 pos) { /* 이 함수는 dest배열에서 pos위치에 src를 더해주는 역할을 합니다. */ int carry = 0; RRU32 i = pos; dest[i] += src; if(dest[i--] < src) carry = 1; for(; i != (RRU32)(-1); --i) { if(carry) { dest[i]++; carry = dest[i] == 0; } } return carry; } inline int i_addu64(RRU32* dest, RRU32 len, RRU64 src, RRU32 pos) { /* 이 함수는 dest배열에서 pos위치 src를 더해주는 역할을 합니다. 그런데 src가 64비트라서, 상위32비트, 하위32비트 둘로 나누어서 더해줍니다. */ int carry = 0; RRU32 i = pos; if(len < pos) return 0; if(len == pos) { return i_addu32(dest, len, (RRU32)(src >> 32), pos - 1); } dest[i] += (RRU32)(src & 0xFFFFFFFF); if(dest[i--] < (RRU32)(src & 0xFFFFFFFF)) carry = 1; if(i == (RRU32)(-1)) return carry; dest[i] += (RRU32)(src >> 32); if(carry) { dest[i]++; carry = dest[i--] <= (RRU32)(src >> 32); } else { carry = dest[i--] < (RRU32)(src >> 32); } for(; i != (RRU32)(-1); --i) { if(carry) { dest[i]++; carry = dest[i] == 0; } } return carry; }

이번에도 도우미함수부터 정의했습니다~~

구현을 보면 알수있지만, 역시나 계산결과가 자리수를 넘어가면... 과감하게 무시할수 밖에 없습니다. (지못미 ㅠㅠ)

이제 진짜 곱셈함수를 구현해볼까요

/* */ void RR_mul(PRR dest, PRR a, PRR b) { memset(dest->frac, 0, SIZE_TO_BYTE(dest->alloced)); //1. 먼저 dest를 0으로 세팅. for(RRU32 j = b->size - 1; j != (RRU32)(-1); --j) { for(RRU32 i = a->size - 1; i != (RRU32)(-1); --i) { //2. 모든 i, j에 대해서 if(a->frac[i] && b->frac[j]) { //곱셈을 실시하고, i+j 자리에 더한다. i_addu64(dest->frac, dest->alloced, __emulu(a->frac[i], b->frac[j]), i + j); } } } dest->sign = a->sign ^ b->sign; //소소한 부호결정 코드. xor연산자를 이용해 부호가 다른경우만 음수가 되도록ㅋ } /* 이 함수는 RealReal포맷과 32비트 정수를 곱해줍니다. 당연히 위 함수보다는 빠르겠지요. */ int RR_mulu32(PRR dest, PRR a, RRU32 b) { int carry = 0; memset(dest->frac, 0, SIZE_TO_BYTE(dest->alloced)); for(RRU32 i = a->size - 1; i < (RRU32)(-1); --i) { if(a->frac[i]) { carry += i_addu64(dest->frac, dest->alloced, __emulu(a->frac[i], b), i); } } dest->sign = a->sign; return carry; }

중간에 이상한 함수가 하나 불쑥 들어가 있습니다. __emulu 함수인데, 이 함수는 어셈블리어 명령 하나와 대응됩니다. 이런함수를 intrinsic함수라고 합니다. 원래 32비트 정수 둘을 곱하면 64비트의 결과가 나오게 되지요. 하지만 c/c++언어의 기본 곱셈연산자는 상위 32비트는 버리고, 하위 32비트만 결과로 남깁니다. __emulu함수는 두 32비트 정수를 곱해서 64비트 결과를 내주는 함수이지요. 얘는 부호없는 버젼이고, 부호있는 정수일 경우에는 __emul을 사용하면 됩니다.

이런 intrinsic함수를 사용하려면 <intrin.h> 헤더를 include해주고, #pragma intrinsic(__emulu)와 같은 전처리표시자를 넣어주어야합니다.

32비트 x86의 어셈블리 연산 mul은 두 32비트 정수를 입력받아, 결과값을 edx레지스터에 상위32비트, eax레지스터에 하위32비트로 나누어 구해줍니다. 만약 intrinsic함수를 사용할수 없는 경우에는 __emulu함수 대신 (unsigned long long)으로 캐스팅해주고 곱셈을 해도 됩니다. 요즘 컴파일러는 충분히 똑똑해서 저 연산 역시 mul 명령어 하나로 번역해주니깐요.

또 보면 for문에서 0 에서 size-1까지 가는게 아니라, size-1에서 0까지 거슬러 올라오고 있는데, (순서는 결과에 영향을 미치지 않습니다만) 이는 우리가 제일 윗자리를 0번 인덱스에, 아랫자리일수록 인덱스 번호가 커지도록 정의했기때문입니다. 아랫자리부터 곱해서 더해나가야지 혹시나 일어날 받아올림을 최소화할수 있기때문이지요.

곱셈을 구현했으니, 이제 10진수로 출력을 할수 있어요~~!

/* 이 함수는 rr을 precision의 정밀도로 출력해줍니다. 만약 precision이 0이라면 디폴트값. */void RR_print(PRR rr, FILE* f, int precision) { if(precision == 0) { precision = rr->size * 9; //32비트 정수는 대략 10진수 유효자리 9개를 가집니다. } if(rr->sign) fprintf(f, "-"); fprintf(f, "%u.", rr->frac[0]);//정수부분 출력 PRR x = RR_create(rr->size); RR_move(x, rr); // x = rr PRR y = RR_create(rr->size); /* 여기서부터 조금 계산이 들어갑니다. 소수부분에 1000000000를 곱하면, 소수점아래 9자리까지가 정수부분으로 올라오게 됩니다. 그러면 그 9자리를 포착해서 출력해주는것이지요. */ for(;precision >= 9; precision -= 9) { x->frac[0] = 0; //정수부분을 날려버립니다. RR_mulu32(y, x, 1000000000); // y = x * rr fprintf(f, "%09u", y->frac[0]); PRR t = x; x = y; y = t; //x,y를 서로 바꾸면서 계속 반복. } x->frac[0] = 0; switch(precision) { case 1: RR_mulu32(y, x, 10); fprintf(f, "%01u", y->frac[0]); break; case 2: RR_mulu32(y, x, 100); fprintf(f, "%02u", y->frac[0]); break; case 3: RR_mulu32(y, x, 1000); fprintf(f, "%03u", y->frac[0]); break; case 4: RR_mulu32(y, x, 10000); fprintf(f, "%04u", y->frac[0]); break; case 5: RR_mulu32(y, x, 100000); fprintf(f, "%05u", y->frac[0]); break; case 6: RR_mulu32(y, x, 1000000); fprintf(f, "%06u", y->frac[0]); break; case 7: RR_mulu32(y, x, 10000000); fprintf(f, "%07u", y->frac[0]); break; case 8: RR_mulu32(y, x, 100000000); fprintf(f, "%08u", y->frac[0]); break; } /* 고집이 세서 스위치문으로 구현했습니다. 좀더 우아한 구현은 과제로...ㅋ */ RR_delete(x); RR_delete(y); }

보면 알수있듯이 2^32진수로 저장되어있는 수치를 10진수로 바꾸는 일은 금방 끝나는 일이 아닙니다... (정밀도가 높아질수록 곱하는 횟수도 늘어나야하지요... ㅠㅠ) 곱셈대신 비트연산자와 덧셈을 이용하는 방법도 있습니다. 이런건 나중에 천천히 알아보도록 하지요.

전체코드보기

//// RealReal.cpp 파일 #include "stdafx.h" #include <intrin.h> #include "RealReal.h" #include "RealRealUtils.h" #pragma intrinsic(__emulu) PRR RR_create(RRU32 size) { PRR rr = (PRR)malloc(sizeof(RR)); rr->alloced = MAX(size, 4); rr->frac = (RRU32*)malloc(SIZE_TO_BYTE(rr->alloced)); memset(rr->frac, 0, SIZE_TO_BYTE(rr->alloced)); rr->size = rr->alloced; rr->sign = 0; return rr; } void RR_delete(PRR rr) { free(rr->frac); free(rr); } void RR_move(PRR dest, PRR src) { dest->sign = src->sign; if(dest->size > src->size) { memcpy(dest->frac, src->frac, SIZE_TO_BYTE(src->size)); RRU32 l = dest->size - src->size; memset(dest->frac + dest->size, 0, SIZE_TO_BYTE(l)); } else { memcpy(dest->frac, src->frac, SIZE_TO_BYTE(dest->size)); } } int RR_init(PRR rr, double value) { if(value > 0xFFFFFFFF) return -1; if(value < 0) { rr->sign = 1; value = -value; } else { rr->sign = 0; } rr->frac[0] = (RRU32)value; value -= rr->frac[0]; value *= 4294967296.; rr->frac[1] = (RRU32)value; value -= rr->frac[1]; value *= 4294967296.; rr->frac[2] = (RRU32)value; memset(rr->frac + 3, 0, SIZE_TO_BYTE(rr->size - 3)); return 0; } void RR_initu(PRR rr, RRU32 value, RRU32 sign) { rr->sign = !!sign; rr->frac[0] = value; memset(rr->frac + 1, 0, SIZE_TO_BYTE(rr->size - 1)); } void RR_printHex(PRR rr, FILE* f) { if(rr->sign) fprintf(f, "-"); fprintf(f, "%x.", rr->frac[0]); for(RRU32 i = 1; i < rr->size; ++i) { fprintf(f, "%08x", rr->frac[i]); } } inline void i_RR_add(PRR dest, PRR src) { i_adds(dest->frac, src->frac, MIN(dest->size, src->size)); } inline void i_RR_sub(PRR dest, PRR src) { if(i_subs(dest->frac, src->frac, MIN(dest->size, src->size))) { dest->sign = ! dest->sign; i_cmplt(dest->frac, dest->size); } } void RR_add(PRR dest, PRR src) { if(dest->sign == src->sign) { i_RR_add(dest, src); } else { i_RR_sub(dest, src); } } void RR_sub(PRR dest, PRR src) { if(dest->sign == src->sign) { i_RR_sub(dest, src); } else { i_RR_add(dest, src); } } void RR_neg(PRR dest) { dest->sign = !dest->sign; } void RR_mul(PRR dest, PRR a, PRR b) { memset(dest->frac, 0, SIZE_TO_BYTE(dest->alloced)); for(RRU32 j = b->size - 1; j != (RRU32)(-1); --j) { for(RRU32 i = a->size - 1; i != (RRU32)(-1); --i) { if(a->frac[i] && b->frac[j]) { i_addu64(dest->frac, dest->alloced, __emulu(a->frac[i], b->frac[j]), i + j); } } } dest->sign = a->sign ^ b->sign; } int RR_mulu32(PRR dest, PRR a, RRU32 b) { int carry = 0; memset(dest->frac, 0, SIZE_TO_BYTE(dest->alloced)); for(RRU32 i = a->size - 1; i < (RRU32)(-1); --i) { if(a->frac[i]) { carry += i_addu64(dest->frac, dest->alloced, __emulu(a->frac[i], b), i); } } dest->sign = a->sign; return carry; } void RR_print(PRR rr, FILE* f, int precision) { if(precision == 0) { precision = rr->size * 9; } if(rr->sign) fprintf(f, "-"); fprintf(f, "%u.", rr->frac[0]); PRR x = RR_create(rr->size); RR_move(x, rr); PRR y = RR_create(rr->size); for(;precision >= 9; precision -= 9) { x->frac[0] = 0; RR_mulu32(y, x, 1000000000); fprintf(f, "%09u", y->frac[0]); PRR t = x; x = y; y = t; } x->frac[0] = 0; switch(precision) { case 1: RR_mulu32(y, x, 10); fprintf(f, "%01u", y->frac[0]); break; case 2: RR_mulu32(y, x, 100); fprintf(f, "%02u", y->frac[0]); break; case 3: RR_mulu32(y, x, 1000); fprintf(f, "%03u", y->frac[0]); break; case 4: RR_mulu32(y, x, 10000); fprintf(f, "%04u", y->frac[0]); break; case 5: RR_mulu32(y, x, 100000); fprintf(f, "%05u", y->frac[0]); break; case 6: RR_mulu32(y, x, 1000000); fprintf(f, "%06u", y->frac[0]); break; case 7: RR_mulu32(y, x, 10000000); fprintf(f, "%07u", y->frac[0]); break; case 8: RR_mulu32(y, x, 100000000); fprintf(f, "%08u", y->frac[0]); break; } RR_delete(x); RR_delete(y); }

//// RealReal.h 파일 #pragma once typedef unsigned char RRU8; typedef unsigned short RRU16; typedef unsigned long RRU32; typedef unsigned long long RRU64; typedef struct tagRealReal { RRU32* frac; // fraction part RRU32 alloced; // allocated RRU32 counts of fraction part RRU32 size; // RRU32 counts RRU32 sign; } RR; typedef RR* PRR; #define SIZE_TO_BYTE(size) ((size)*sizeof(RRU32)) #define MAX(a, b) (((a) > (b)) ? (a) : (b)) #define MIN(a, b) (((a) < (b)) ? (a) : (b)) PRR RR_create(RRU32 size); void RR_delete(PRR rr); void RR_move(PRR dest, PRR src); // dest = src int RR_init(PRR rr, double value); // rr = value // If value is too large to put rr, returns -1, otherwise returns 0 void RR_initu(PRR rr, RRU32 value, RRU32 sign); void RR_printHex(PRR rr, FILE* f); void RR_print(PRR rr, FILE* f, int precision); void RR_add(PRR dest, PRR src); // dest += src void RR_sub(PRR dest, PRR src); // dest -= src void RR_neg(PRR dest); // dest = -dest void RR_mul(PRR dest, PRR a, PRR b); // dest = a*b int RR_mulu32(PRR dest, PRR a, RRU32 b); // dest = a*b

//// RealRealUtils.h 파일 #pragma once inline int i_adds(RRU32* dest, const RRU32* src, RRU32 len) { int carry = 0; for(RRU32 i = len - 1; i != (RRU32)(-1); --i) { dest[i] += src[i]; if(carry) { dest[i]++; carry = dest[i] <= src[i]; } else { carry = dest[i] < src[i]; } } return carry; } inline int i_addu32(RRU32* dest, RRU32 len, RRU32 src, RRU32 pos) { int carry = 0; RRU32 i = pos; dest[i] += src; if(dest[i--] < src) carry = 1; for(; i != (RRU32)(-1); --i) { if(carry) { dest[i]++; carry = dest[i] == 0; } } return carry; } inline int i_addu64(RRU32* dest, RRU32 len, RRU64 src, RRU32 pos) { int carry = 0; RRU32 i = pos; if(len < pos) return 0; if(len == pos) { return i_addu32(dest, len, (RRU32)(src >> 32), pos - 1); } dest[i] += (RRU32)(src & 0xFFFFFFFF); if(dest[i--] < (RRU32)(src & 0xFFFFFFFF)) carry = 1; if(i == (RRU32)(-1)) return carry; dest[i] += (RRU32)(src >> 32); if(carry) { dest[i]++; carry = dest[i--] <= (RRU32)(src >> 32); } else { carry = dest[i--] < (RRU32)(src >> 32); } for(; i != (RRU32)(-1); --i) { if(carry) { dest[i]++; carry = dest[i] == 0; } } return carry; } inline int i_subs(RRU32* dest, const RRU32* src, RRU32 len) { int carry = 0; for(RRU32 i = len - 1; i != (RRU32)(-1); --i) { if(carry) { carry = dest[i] <= src[i]; dest[i]--; } else { carry = dest[i] < src[i]; } dest[i] -= src[i]; } return carry; } inline void i_cmplt(RRU32* dest, RRU32 len) { RRU32 i; for(i = len - 1; i != (RRU32)(-1); --i) { if(dest[i] > 0) { dest[i]--; break; } } for(; i != (RRU32)(-1); --i) { dest[i] = ~dest[i]; } }