multi precision arithmatic

머리가 안 따라주면 주말에도 놀지 말고 일해야 한다. 어렸을 적에 짰던 코드와 하등 다를 것이 없는 것을 짰다. 그런데 누쓰 알고리즘이 어떻게 되는지 잊어버렸다. 책을 찾아봐야 하는데, 그의 책은 옛날 옛날 먼 옛날에 고구마 궈 먹으려고 들판에서 태워버렸다. n-bit precision arithmatic 중 사칙연산과 논리연산만 있으면 되는데 토요일 하루 종일 안 돌아가는 머리로 코드를 짰다. 곱셈과 나눗셈 계산 하는 방법을 잊어먹은 것이다. 제길헐... 간신히 되긴 했는데 내가 왜 이 코드를 짜게 되었는지 그 '목적'을 잊어 버렸다. OTL

지금은 어떻게 연산 속도를 높일 수 있을까 궁리를 하다가... 속도가 별 문제되지 않는다는 점을 새삼스럽게 깨달았다. 난 대체 뭐하는 놈일까...

typedef unsigned int Unit;

class Integer
{
public:

#ifdef _WIN32
enum { // ms vc 6.0 has some bug on static const init in class definition
bits = 2048,
stsize = sizeof(Unit) * 8,
bytes = bits / 8,
MaxSize = bits / stsize,
};
#else // for linux
static const int bits = 2048;
static const int stsize = sizeof(Unit) * 8;
static const int bytes = bits / 8;
static const int MaxSize = bits / stsize;
#endif

Unit d[MaxSize];
...
};

int Integer::bitcount()
{
int kn = 0;

for (int i = MaxSize - 1; i >= 0; i--) {
if (d[i] != 0)
break;
}

for (int j = stsize - 1; j >= 0; j--) {
if ((d[i] & (1 << j)) != 0)
return i * stsize + j;
}

return kn;
}


Integer& Integer::operator+ (const Integer& x)
{
result = *this;
unsigned __int64 carry = 0;

for (int i = 0; i < MaxSize; i++) {
carry += result.d[i];
carry += x.d[i];
result.d[i] = (Unit)carry;
carry >>= (stsize);
}
return result;
}

Integer& Integer::operator- (const Integer &x)
{
Integer s = x;;
s.neg(); // 2's complementary
result = *this + s;
return result;
}

Integer& Integer::operator* (const Integer& x)
{
Integer N(*this); // multiplee
Integer M(x); // multiplier

result.clear();
if (M.is_zero()) {
return result;
}

int i, kn = M.bitcount();

M = M << (bits-kn-1);
for (i = 0; i < kn + 1; i++) {
result.shift_left(0);
if (M.shift_left(0))
result = result + N;
}
return result;
}

void Integer::div(Integer D, int modf)
{
// N = Q * D + R,
// where N is given numerator, Q = quotient, D = divisor, R = remainder

Integer N(*this);
Integer R;
Integer Q;
int i, carry, kn = bitcount();

if (N.is_zero() || D.is_zero()) { // D is zero: divide by zero
*this = R; // R or Q = 0
return;
}

if (N < D) {
* this = (modf) ? N : R; // N = remainder, R = 0
return;
}

N = N << (bits-kn-1); // left justify the numerator
for (i = 0; i < kn+1; i++) {
carry = N.shift_left(0);
R.shift_left(carry); // rotate left with carry
if (R >= D) {
Q.shift_left(1);
R = R - D;
}
else {
Q.shift_left(0);
}
}

* this = (modf) ? R : Q;
return;
}

Integer& Integer::operator/ (const Integer& x)
{
result = *this;
result.div(x, 0); // get quotient
return result;
}
Integer& Integer::operator% (const Integer& x)
{
result = *this;
result.div(x, 1); // get remainder
return result;
}

...

void main(void)
{
Integer x = "0x1234`5678`9abc`def0`1234`5678`9abc`def0";
Integer y = "0b00011001`11010101";
Integer z = "0x7E315816B6472299DF";
Integer p = "2327845370165345819103"; // same as z

cout << "x * y = " << (x * y).string(16) << "\n"; // represent value as hex string
cout << "x / y = " << (x / y).string(16) << "\n";
cout << "x << r = " << (x << p).string(2) << "\n"; // bin string
}