상세 컨텐츠

본문 제목

[진짜 BigFloat 구현하기] 3. 곱셈 구현하기

프로그래밍/Multi precision

by ∫2tdt=t²+c 2013. 11. 26. 12:05

본문

앞서서 덧셈과 뺄셈을 구현해보았는데요, 이미 둘을 구현했다면 곱셈을 구현하는 것은 생각보다 간단합니다. 알고리즘에 대해서는 다음을 참조하시면 됩니다.


2012/10/18 - [프로그래밍/Multi precision] - BigFloat로 Pi를 구해보자-3. 곱셈 구현 + 10진수로 출력하기


제일 단순한 방법으로 곱셈을 구현하면 n길이인 정수에 대해 O(n^2) 시간에 계산을 완료할 수 있습니다. 더 빠른 알고리즘 구현은 일단 추후로 미루고ㅋ



당연하겠지만 지수부와 정수배열로 표현된 두 수의 곱은 정수배열끼리는 곱하고 지수부는 더하면 됩니다. 부호연산은 자명하니 설명이 필요없구요. 바로 코드를 보시죠.



1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
// xl_mul : a와 b의 곱을 구하여 dest에 저장한다
void xl_mul(PXLARGE dest, PXLARGE a, PXLARGE b)
{
    if(xl_isnan(a) || xl_isnan(b)) return xl_setnan(dest);
    XL_UINT32 aLen = _XL_GetLength(a);
    XL_UINT32 bLen = _XL_GetLength(b);
    xl_init(dest, 0);
//연산 결과는 a의 길이와 b의 길이의 합입니다.
//처음 크기가 1이었으므로 aLen + bLen - 1를 뒤쪽에 넣으면
//전체 길이는 aLen + bLen이 되겠죠.
    _XL_PushBackN(dest, aLen + bLen - 1);
     
// [0, aLen)과 [0, bLen)인 두 정수에 대해 루프를 돌립니다.
    for(XL_UINT32 i = 0; i < aLen; ++i) for(XL_UINT32 j = 0; j < bLen; ++j)
    {
        XL_UINT32 ad = *_XL_GetFrac(a, i);
        XL_UINT32 bd = *_XL_GetFrac(b, j);
// 둘 중에 하나라도 0이면 결과가 0이므로 더할 필요가 없음
        if(ad && bd)
        {
// 위에서 dest의 데이터가 일렬로 배치되도록 했기에,
// 두 덩어리로 나뉘는 경우를 고려할 필요가 없습니다.
            i_addu64(&dest->frac[dest->begin + i + j], dest->end - dest->begin - i - j, __emulu(ad, bd));
        }
    }
    _XL_Shrink(dest);
    dest->sign = a->sign ^ b->sign;
// 부호는 xor연산
    dest->exp += a->exp + b->exp;
// 지수부는 덧셈 연산
}


곱셈만 구현하고 끝낼수 없죠. 곱셈한 결과를 더해주는 FMA(Fused Multiply-Add)를 구현해봅시다. 이 연산이 은근히 유용합니다.



수식으로 정리하면 위 처럼 되겠죠? 실제 구현에서는 절대값을 더해주는 경우와 빼주는 경우로 나눠서 구현을 해야합니다.


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
inline void _Abs_MulAdd(PXLARGE dest, PXLARGE a, PXLARGE b)
{
    XL_UINT32 aLen = _XL_GetLength(a);
    XL_UINT32 bLen = _XL_GetLength(b);
    XL_INT32 tExp = a->exp + b->exp;
    XL_INT32 tBack = tExp + aLen + bLen + 1; // carry를 대비한 여분의 +1
    if(dest->exp > tExp) _XL_PushFrontN(dest, dest->exp - tExp);
    XL_INT32 dBack = dest->exp + _XL_GetLength(dest);
    if(dBack < tBack) _XL_PushBackN(dest, tBack - dBack);
 
    XL_UINT32 orgDestBegin = dest->begin;
    if(tExp > dest->exp) dest->begin += tExp - dest->exp;
    if(dest->begin > dest->size) dest->begin -= dest->size;
 
// dest 데이터가 일렬로 배치된 경우
    if(dest->begin < dest->end)
    {
        for(XL_UINT32 i = 0; i < aLen; ++i) for(XL_UINT32 j = 0; j < bLen; ++j)
        {
            XL_UINT32 ad = *_XL_GetFrac(a, i);
            XL_UINT32 bd = *_XL_GetFrac(b, j);
            if(ad && bd)
            {
                i_addu64(&dest->frac[dest->begin + i + j], dest->end - dest->begin - i - j, __emulu(ad, bd));
            }
        }
    }
// dest 데이터가 두 덩어리로 나눠진 경우
    else
    {
        for(XL_UINT32 i = 0; i < aLen; ++i) for(XL_UINT32 j = 0; j < bLen; ++j)
        {
            XL_UINT32 ad = *_XL_GetFrac(a, i);
            XL_UINT32 bd = *_XL_GetFrac(b, j);
            if(ad && bd)
            {
                if(dest->begin + i + j < dest->size)
                {
                    i_addu64_2(&dest->frac[dest->begin + i + j], dest->size - dest->begin - i - j, dest->frac, dest->end, __emulu(ad, bd));
                }
                else
                {
                    i_addu64(&dest->frac[dest->begin + i + j - dest->size], dest->end + dest->size - dest->begin - i - j, __emulu(ad, bd));
                }
            }
        }
    }
    dest->begin = orgDestBegin;
    _XL_Shrink(dest);
}
 
inline void _Abs_MulSub(PXLARGE dest, PXLARGE a, PXLARGE b)
{
    XL_UINT32 aLen = _XL_GetLength(a);
    XL_UINT32 bLen = _XL_GetLength(b);
    XL_INT32 tExp = a->exp + b->exp;
    XL_INT32 tBack = tExp + aLen + bLen;
    if(dest->exp > tExp) _XL_PushFrontN(dest, dest->exp - tExp);
    XL_INT32 dBack = dest->exp + _XL_GetLength(dest);
    if(dBack < tBack) _XL_PushBackN(dest, tBack - dBack);
 
    XL_UINT32 orgDestBegin = dest->begin;
    if(tExp > dest->exp) dest->begin += tExp - dest->exp;
    if(dest->begin > dest->size) dest->begin -= dest->size;
 
// 부호가 뒤집혔는 지 여부를 나타내는 변수
    bool invert = false;
 
// dest 데이터가 일렬로 배치된 경우
    if(dest->begin < dest->end)
    {
        for(XL_UINT32 i = 0; i < aLen; ++i) for(XL_UINT32 j = 0; j < bLen; ++j)
        {
            XL_UINT32 ad = *_XL_GetFrac(a, i);
            XL_UINT32 bd = *_XL_GetFrac(b, j);
            if(ad && bd)
            {
// 부호가 뒤집히지 않았으면 i_subu64, 뒤집혔으면 i_addu64 호출
                if((invert ? i_addu64 : i_subu64)(&dest->frac[dest->begin + i + j], dest->end - dest->begin - i - j, __emulu(ad, bd)))
                {
                    // 오버플로우가 발생한 경우 부호가 뒤집힌것
                    invert = true;
                    dest->sign = !dest->sign;
                    if(orgDestBegin < dest->end)
                    {
                        i_cmplt(&dest->frac[orgDestBegin], dest->end - orgDestBegin, 0);
                    }
                    else
                    {
                        int carry = i_cmplt(&dest->frac[orgDestBegin], dest->size - orgDestBegin, 0);
                        if(dest->end) i_cmplt(dest->frac, dest->end, carry);
                    }
                     
                }
            }
        }
    }
// dest 데이터가 두 덩어리로 나눠진 경우
    else
    {
        for(XL_UINT32 i = 0; i < aLen; ++i) for(XL_UINT32 j = 0; j < bLen; ++j)
        {
            XL_UINT32 ad = *_XL_GetFrac(a, i);
            XL_UINT32 bd = *_XL_GetFrac(b, j);
            if(ad && bd)
            {
                if(dest->begin + i + j < dest->size)
                {
// 부호가 뒤집히지 않았으면 i_subu64_2, 뒤집혔으면 i_addu64_2 호출
                    if((invert ? i_addu64_2 : i_subu64_2)(&dest->frac[dest->begin + i + j], dest->size - dest->begin - i - j, dest->frac, dest->end, __emulu(ad, bd)))
                    {
                        // 부호 뒤집기
                        invert = true;
                        dest->sign = !dest->sign;
                        int carry = i_cmplt(&dest->frac[orgDestBegin], dest->size - orgDestBegin, 0);
                        if(dest->end) i_cmplt(dest->frac, dest->end, carry);
                    }
                }
                else
                {
// 부호가 뒤집히지 않았으면 i_subu64, 뒤집혔으면 i_addu64 호출
                    if((invert ? i_addu64 : i_subu64)(&dest->frac[dest->begin + i + j - dest->size], dest->end + dest->size - dest->begin - i - j, __emulu(ad, bd)))
                    {
                        // 부호 뒤집기
                        invert = true;
                        dest->sign = !dest->sign;
                        int carry = i_cmplt(&dest->frac[orgDestBegin], dest->size - orgDestBegin, 0);
                        if(dest->end) i_cmplt(dest->frac, dest->end, carry);
                    }
                }
            }
        }
    }
    dest->begin = orgDestBegin;
    _XL_Shrink(dest);
}


진짜 FMA를 하는 함수는 위 함수를 이용해서 다음과 같이 심플하게 짤수 있겠죠.


1
2
3
4
5
6
7
8
9
10
11
void xl_muladd(PXLARGE dest, PXLARGE a, PXLARGE b)
{
    if(xl_isnan(dest) || xl_isnan(a) || xl_isnan(b)) return xl_setnan(dest);
    return (dest->sign == (a->sign ^ b->sign) ? _Abs_MulAdd : _Abs_MulSub)(dest, a, b);
}
 
void xl_mulsub(PXLARGE dest, PXLARGE a, PXLARGE b)
{
    if(xl_isnan(dest) || xl_isnan(a) || xl_isnan(b)) return xl_setnan(dest);
    return (dest->sign != (a->sign ^ b->sign) ? _Abs_MulAdd : _Abs_MulSub)(dest, a, b);
}


이제 곱셈이 완료되었으니, 10진수로 출력하는 함수를 만들수 있습니다. 그건 다음으로 패스!




관련글 더보기