High-Precision Computation

TL;DR: Get the template at the end...

Definition

High-precision computation (Arbitrary-Precision Arithmetic), also called big-number (bignum) computation, uses algorithmic structures to support operations between larger integers (numbers whose sizes exceed the integer types built into the language).

Introduction

High-precision problems contain many small details, and their implementations also require considerable care.

So today we will implement a simple calculator together.

Storage

In ordinary implementations, high-precision numbers are represented by strings, with each character representing one decimal digit. Thus, high-precision numerical computation can actually be regarded as a special kind of string processing.

When a string is read, the most significant digit is at the beginning of the string (the position with the smaller index). By convention, however, the smallest index stores the least significant digit, that is, the string is stored in reverse. This is because the length of a number may change, while we want positions with the same place value to remain aligned (for example, all ones digits at index [0] and all tens digits at index [1]...); at the same time, addition, subtraction, and multiplication generally start from the ones digit (recall the column arithmetic taught in elementary school). These are all good reasons for storing the string in reverse.

We will follow this convention throughout. Define a constant LEN = 1004 as the maximum length accommodated by the program.

It is then easy to write code to read a high-precision number:

void clear(int a[]) {
  for (int i = 0; i < LEN; ++i) a[i] = 0;
}

void read(int a[]) {
  static char s[LEN + 1];
  scanf("%s", s);

  clear(a);

  int len = strlen(s);
  // As described above, reverse the string
  for (int i = 0; i < len; ++i) a[len - i - 1] = s[i] - '0';
  // s[i] - '0' is the digit represented by s[i]
  // Some readers may find it more intuitive to understand this as ord(s[i]) - ord('0')
}

Output also follows the reverse order of storage. Since we do not want to output leading zeroes, we search downward from the highest position for the first nonzero position and start outputting there. The termination condition i >= 1 rather than i >= 0 is because we still want to output a single character 0 when the entire number equals $0$.

void print(int a[]) {
  int i;
  for (i = LEN - 1; i >= 1; --i)
    if (a[i] != 0) break;
  for (; i >= 0; --i) putchar(a[i] + '0');
  putchar('\n');
}

Putting these together gives a complete copycat program.

#include <cstdio>
#include <cstring>

constexpr int LEN = 1004;

int a[LEN];

void clear(int a[]) {
  for (int i = 0; i < LEN; ++i) a[i] = 0;
}

void read(int a[]) {
  static char s[LEN + 1];
  scanf("%s", s);

  clear(a);

  int len = strlen(s);
  for (int i = 0; i < len; ++i) a[len - i - 1] = s[i] - '0';
}

void print(int a[]) {
  int i;
  for (i = LEN - 1; i >= 1; --i)
    if (a[i] != 0) break;
  for (; i >= 0; --i) putchar(a[i] + '0');
  putchar('\n');
}

int main() {
  read(a);
  print(a);

  return 0;
}

Four Arithmetic Operations

The four arithmetic operations also differ in difficulty. The simplest are high-precision addition and subtraction, followed by high-precision—single-precision (ordinary int) multiplication and high-precision—high-precision multiplication, and finally high-precision—high-precision division.

We will implement all the required functions in this order.

Addition

High-precision addition is really just column addition.

That is, starting from the least significant digit, add the digits at corresponding positions of the two addends and determine whether the result reaches or exceeds $10$. If it does, handle the carry: add $1$ to the result at the next higher position and subtract $10$ from the result at the current position.

void add(int a[], int b[], int c[]) {
  clear(c);

  // In a big-number implementation, LEN is generally set larger than the possible input
  // The last few iterations are then omitted, which avoids handling many boundary cases
  // Because the actual input does not exceed 1000 digits, iterating to LEN - 1 = 1003 here is sufficient
  for (int i = 0; i < LEN - 1; ++i) {
    // Add the digits at the corresponding positions
    c[i] += a[i] + b[i];
    if (c[i] >= 10) {
      // Carry
      c[i + 1] += 1;
      c[i] -= 10;
    }
  }
}

Combining this with the previous section gives an addition calculator.

#include <cstdio>
#include <cstring>

constexpr int LEN = 1004;

int a[LEN], b[LEN], c[LEN];

void clear(int a[]) {
  for (int i = 0; i < LEN; ++i) a[i] = 0;
}

void read(int a[]) {
  static char s[LEN + 1];
  scanf("%s", s);

  clear(a);

  int len = strlen(s);
  for (int i = 0; i < len; ++i) a[len - i - 1] = s[i] - '0';
}

void print(int a[]) {
  int i;
  for (i = LEN - 1; i >= 1; --i)
    if (a[i] != 0) break;
  for (; i >= 0; --i) putchar(a[i] + '0');
  putchar('\n');
}

void add(int a[], int b[], int c[]) {
  clear(c);

  for (int i = 0; i < LEN - 1; ++i) {
    c[i] += a[i] + b[i];
    if (c[i] >= 10) {
      c[i + 1] += 1;
      c[i] -= 10;
    }
  }
}

int main() {
  read(a);
  read(b);

  add(a, b, c);
  print(c);

  return 0;
}

Subtraction

High-precision subtraction is likewise column subtraction.

Subtract digit by digit starting from the ones digit; when the result is negative, borrow $1$ from the next higher position. The overall idea is exactly the same as for addition.

void sub(int a[], int b[], int c[]) {
  clear(c);

  for (int i = 0; i < LEN - 1; ++i) {
    // Subtract digit by digit
    c[i] += a[i] - b[i];
    if (c[i] < 0) {
      // Borrow
      c[i + 1] -= 1;
      c[i] += 10;
    }
  }
}

Replacing add() in the previous program with sub() gives a subtraction calculator.

#include <cstdio>
#include <cstring>

constexpr int LEN = 1004;

int a[LEN], b[LEN], c[LEN];

void clear(int a[]) {
  for (int i = 0; i < LEN; ++i) a[i] = 0;
}

void read(int a[]) {
  static char s[LEN + 1];
  scanf("%s", s);

  clear(a);

  int len = strlen(s);
  for (int i = 0; i < len; ++i) a[len - i - 1] = s[i] - '0';
}

void print(int a[]) {
  int i;
  for (i = LEN - 1; i >= 1; --i)
    if (a[i] != 0) break;
  for (; i >= 0; --i) putchar(a[i] + '0');
  putchar('\n');
}

void sub(int a[], int b[], int c[]) {
  clear(c);

  for (int i = 0; i < LEN - 1; ++i) {
    c[i] += a[i] - b[i];
    if (c[i] < 0) {
      c[i + 1] -= 1;
      c[i] += 10;
    }
  }
}

int main() {
  read(a);
  read(b);

  sub(a, b, c);
  print(c);

  return 0;
}

Try it: input 1 2—the output is /9999999. Hey, why did this OI Wiki give me fake code...

In fact, the code above can only handle the case where the minuend $a$ is greater than or equal to the subtrahend $b$. To handle the case where the minuend is smaller than the subtrahend, namely $aa$, call the sub function in the code above as sub(b,a,c). To obtain the value of $a-b$, simply add a minus sign before the result.

Multiplication

High-Precision—Single-Precision

High-precision multiplication is also column... wait, wait!

First consider a simple case: one of the factors is of the ordinary int type. Is there a simple way to handle it?

An intuitive idea is to multiply each digit of $a$ directly by $b$. Numerically, this method is correct, but it does not conform to decimal notation, so it must be rearranged into the normal form.

The rearrangement also processes carries upward digit by digit starting from the ones digit. However, the carry here may be very large, even much greater than $9$, because after multiplication each position may reach the order of $9b$. Therefore, the carry cannot simply be processed with a $-10$ operation; instead, it must be calculated using the quotient and remainder after division by $10$. See the comments in the code, or refer to the figure below, which shows the process of multiplying the high-precision number $1337$ by the single-precision number $42$.

For the same reason, this method requires special attention to the range of the multiplier $b$. If it is of the same order as $10^9$ (or the upper bound of the corresponding integer type), high-precision—single-precision multiplication should be used with caution.

void mul_short(int a[], int b, int c[]) {
  clear(c);

  for (int i = 0; i < LEN - 1; ++i) {
    // Directly multiply the digit at position i of a by the multiplier and add it to the result
    c[i] += a[i] * b;

    if (c[i] >= 10) {
      // Handle the carry
      // c[i] / 10, the quotient of the division, becomes the increment to the carry
      c[i + 1] += c[i] / 10;
      // c[i] % 10, the remainder of the division, becomes the value left at the current position
      c[i] %= 10;
    }
  }
}

High-Precision—High-Precision

If both factors are high-precision numbers, column multiplication can show its strength again.

Recall each step of column multiplication: it actually calculates the sum of several $a \times b_i \times 10^i$. For example, calculating $1337 \times 42$ means calculating $1337 \times 2 \times 10^0 + 1337 \times 4 \times 10^1$.

Thus, $b$ can be decomposed into all its digits, each of which is a single-precision number. Multiply each of them by $a$, shift the results left to their respective positions, and add them to obtain the answer. Finally, carries must also be handled in the same way as in the example above.

Note that this process is not exactly the same as column multiplication. In our algorithm, carries are not handled during each multiplication step; all results are kept in their corresponding positions and carries are handled together at the end. This does not affect the result.

void mul(int a[], int b[], int c[]) {
  clear(c);

  for (int i = 0; i < LEN - 1; ++i) {
    // Directly calculate the i-th position of the result from low to high, handling the carry at the same time
    // The i-th iteration adds to c[i] the sum of the products of all a[p] and b[q] satisfying p + q = i
    // This has the same effect as performing and summing the calculation in the figure directly, but is a shorter implementation
    for (int j = 0; j <= i; ++j) c[i] += a[j] * b[i - j];

    if (c[i] >= 10) {
      c[i + 1] += c[i] / 10;
      c[i] %= 10;
    }
  }
}

Division

One implementation of high-precision division is column long division.

Column long division can actually be viewed as a process of repeated subtraction. For example, the calculation of the tens digit of the quotient in the figure above can be understood as follows: starting with $45$, subtract $12$ three times until it becomes less than $12$ and cannot be subtracted again, so this digit is $3$.

To reduce redundant operations, we first obtain the lengths $l_a$ of the dividend and $l_b$ of the divisor, then calculate the quotient from high to low starting at index $l_a - l_b$. This is the same as aligning the highest digit of the first product with the highest digit of the dividend in hand calculation.

The reference program implements a function greater_eq() to determine whether the divisor can be subtracted again while keeping the dividend nonnegative when last_dg is the lowest position of the dividend. For each digit of the quotient, it repeatedly calls greater_eq() and, when the condition holds, uses high-precision subtraction to subtract the divisor from the remainder, thereby simulating long division.

// Whether dividend a, with last_dg as its lowest position, can subtract divisor b again while remaining nonnegative
// len is the length of divisor b, avoiding repeated calculation
bool greater_eq(int a[], int b[], int last_dg, int len) {
  // The remaining part of the dividend may be longer than the divisor; in this case it is at most one digit longer, so this check is sufficient
  if (a[last_dg + len] != 0) return true;
  // Compare digit by digit from high to low
  for (int i = len - 1; i >= 0; --i) {
    if (a[last_dg + i] > b[i]) return true;
    if (a[last_dg + i] < b[i]) return false;
  }
  // The equal case is also valid
  return true;
}

void div(int a[], int b[], int c[], int d[]) {
  clear(c);
  clear(d);

  int la, lb;
  for (la = LEN - 1; la > 0; --la)
    if (a[la - 1] != 0) break;
  for (lb = LEN - 1; lb > 0; --lb)
    if (b[lb - 1] != 0) break;
  if (lb == 0) {  // The divisor cannot be zero
    puts("> <");
    return;
  }

  // c is the quotient
  // d is the remaining part of the dividend and naturally becomes the remainder when the algorithm ends
  for (int i = 0; i < la; ++i) d[i] = a[i];
  for (int i = la - lb; i >= 0; --i) {
    // Calculate position i of the quotient
    while (greater_eq(d, b, i, lb)) {
      // Subtract if possible
      // This part is a big-number subtraction
      for (int j = 0; j < lb; ++j) {
        d[i + j] -= b[j];
        if (d[i + j] < 0) {
          d[i + j + 1] -= 1;
          d[i + j] += 10;
        }
      }
      // Increase this position of the quotient by 1
      c[i] += 1;
      // Return to the beginning of the loop and check again
    }
  }
}

The Introductory Section Is Complete!

Combining the implementations of the four arithmetic operations introduced above completes the calculator program mentioned at the beginning.

#include <cstdio>
#include <cstring>

constexpr int LEN = 1004;

int a[LEN], b[LEN], c[LEN], d[LEN];

void clear(int a[]) {
  for (int i = 0; i < LEN; ++i) a[i] = 0;
}

void read(int a[]) {
  static char s[LEN + 1];
  scanf("%s", s);

  clear(a);

  int len = strlen(s);
  for (int i = 0; i < len; ++i) a[len - i - 1] = s[i] - '0';
}

void print(int a[]) {
  int i;
  for (i = LEN - 1; i >= 1; --i)
    if (a[i] != 0) break;
  for (; i >= 0; --i) putchar(a[i] + '0');
  putchar('\n');
}

void add(int a[], int b[], int c[]) {
  clear(c);

  for (int i = 0; i < LEN - 1; ++i) {
    c[i] += a[i] + b[i];
    if (c[i] >= 10) {
      c[i + 1] += 1;
      c[i] -= 10;
    }
  }
}

void sub(int a[], int b[], int c[]) {
  clear(c);

  for (int i = 0; i < LEN - 1; ++i) {
    c[i] += a[i] - b[i];
    if (c[i] < 0) {
      c[i + 1] -= 1;
      c[i] += 10;
    }
  }
}

void mul(int a[], int b[], int c[]) {
  clear(c);

  for (int i = 0; i < LEN - 1; ++i) {
    for (int j = 0; j <= i; ++j) c[i] += a[j] * b[i - j];

    if (c[i] >= 10) {
      c[i + 1] += c[i] / 10;
      c[i] %= 10;
    }
  }
}

bool greater_eq(int a[], int b[], int last_dg, int len) {
  if (a[last_dg + len] != 0) return true;
  for (int i = len - 1; i >= 0; --i) {
    if (a[last_dg + i] > b[i]) return true;
    if (a[last_dg + i] < b[i]) return false;
  }
  return true;
}

void div(int a[], int b[], int c[], int d[]) {
  clear(c);
  clear(d);

  int la, lb;
  for (la = LEN - 1; la > 0; --la)
    if (a[la - 1] != 0) break;
  for (lb = LEN - 1; lb > 0; --lb)
    if (b[lb - 1] != 0) break;
  if (lb == 0) {
    puts("> <");
    return;
  }

  for (int i = 0; i < la; ++i) d[i] = a[i];
  for (int i = la - lb; i >= 0; --i) {
    while (greater_eq(d, b, i, lb)) {
      for (int j = 0; j < lb; ++j) {
        d[i + j] -= b[j];
        if (d[i + j] < 0) {
          d[i + j + 1] -= 1;
          d[i + j] += 10;
        }
      }
      c[i] += 1;
    }
  }
}

int main() {
  read(a);

  char op[4];
  scanf("%s", op);

  read(b);

  switch (op[0]) {
    case '+':
      add(a, b, c);
      print(c);
      break;
    case '-':
      sub(a, b, c);
      print(c);
      break;
    case '*':
      mul(a, b, c);
      print(c);
      break;
    case '/':
      div(a, b, c, d);
      print(c);
      print(d);
      break;
    default:
      puts("> <");
  }

  return 0;
}

Packed-Digit High Precision

Introduction

In ordinary high-precision addition, subtraction, and multiplication, we split the numbers involved in the operation into individual digits and perform the operation on them.

For example, when calculating $8192\times 42$, the high-precision-by-high-precision method actually calculates $(8000+100+90+2)\times(40+2)$.

When there are many digits, there are also many digits after splitting, and the efficiency of high-precision operations decreases.

Can this be optimized?

Notice that the way in which a number is split does not affect the final result, so we can combine several digits.

Process

Using the example above again, if we split the number into groups of two digits, we can write it as $(8100+92)\times 42$.

This splitting does not affect the final result, but because fewer numbers are produced, the calculation becomes more efficient.

From the perspective of positional numeral systems, we perform the calculation in a larger base (splitting a number every two digits above can be regarded as calculating in base $100$), thereby reducing the number of digits involved in the operation and improving efficiency.

This is the idea of packed-digit high precision.

Below is addition code for packed-digit high precision, which further illustrates its implementation:

// The arrays a, b, and c here are numbers in base p
// The number must be converted to decimal when the answer is finally output
void add(int a[], int b[], int c[]) {
  clear(c);

  for (int i = 0; i < LEN - 1; ++i) {
    c[i] += a[i] + b[i];
    if (c[i] >= p) {  // In ordinary big-number arithmetic, p=10
      c[i + 1] += 1;
      c[i] -= p;
    }
  }
}

Efficient Long Division with Packed-Digit High Precision

When using packed-digit high precision, if we still use the method introduced above for trial division, the number of trials will be large and the constant factor will be very high. For example, in base ten thousand, each digit needs an average of 5000 trial quotients, and this huge constant is unacceptable. Therefore, we need a more efficient trial-quotient method.

We can use double as an intermediate. Suppose the dividend has four digits, $a_4,a_3,a_2,a_1$, and the divisor has three digits, $b_3,b_2,b_1$. Then we only need to try one quotient digit: in base $base$, estimate the quotient using $\dfrac{a_4 base + a_3}{b_3 + b_2 base^{-1} + (b_1+1)base^{-2}}$. For multiple digits, use the one-digit form in a loop. Since three digits of the divisor are used for the estimate, the estimated quotient q' and actual quotient q satisfy $q-1 \le q' \le q$, so in the worst case each digit needs only two trial quotients. At the same time, $base^3$ must be within the effective precision of double, namely $base^3 < 2^{53}$. Therefore, when using this method, it is recommended not to exceed base 32768; otherwise, insufficient precision can easily cause errors and incorrect results.

In addition, since the estimated quotient is always less than or equal to the actual quotient, there is room for further optimization. In most cases, each digit is estimated only once. When estimating the next digit, the resulting quotient may be greater than or equal to base because of the error in the preceding digit, but this does not matter: a unified carry operation at the end is sufficient. For example, suppose base is 10 and we calculate $395081/9876$. The trial-quotient calculation is as follows:

1. First, trial division gives $3950/988=3$, so $395081-(9876 \times 3 \times 10^1) = 98801$. An error occurs at this step, but it can be ignored; continue to the next step.
2. Continue trial division on remainder 98801 to obtain $9880/988=10$, so $98801-(9876 \times 10 \times 10^0) = 41$, which is the final remainder.
3. Add the results of the trial-quotient process and handle the carry: $3 \times 10^1 + 10 \times 10^0 = 40$ is the exact quotient.

Although the method looks simple, its concrete implementation is easy to get wrong. The following is a carefully verified implementation for reference; details requiring attention are included in the comments.

// Complete template and implementation https://baobaobear.github.io/post/20210228-bigint1/
// Subtract b multiplied by mul and shifted left by offset, for division
BigIntSimple &sub_mul(const BigIntSimple &b, int mul, int offset) {
  if (mul == 0) return *this;
  int borrow = 0;
  // Unlike subtraction, borrow may be very large, so the subtraction form cannot be used
  for (size_t i = 0; i < b.v.size(); ++i) {
    borrow += v[i + offset] - b.v[i] * mul - BIGINT_BASE + 1;
    v[i + offset] = borrow % BIGINT_BASE + BIGINT_BASE - 1;
    borrow /= BIGINT_BASE;
  }
  // Continue processing if there is still a borrow
  for (size_t i = b.v.size(); borrow; ++i) {
    borrow += v[i + offset] - BIGINT_BASE + 1;
    v[i + offset] = borrow % BIGINT_BASE + BIGINT_BASE - 1;
    borrow /= BIGINT_BASE;
  }
  return *this;
}

BigIntSimple div_mod(const BigIntSimple &b, BigIntSimple &r) const {
  BigIntSimple d;
  r = *this;
  if (absless(b)) return d;
  d.v.resize(v.size() - b.v.size() + 1);
  // Precompute the reciprocal of the highest three digits of the divisor plus 1; if the highest three digits are a3, a2, a1
  // then db is the reciprocal of a3+a2/base+(a1+1)/base^2, and multiplication is finally used to estimate each quotient digit
  // This method can be used within the int32 range when BIGINT_BASE<=32768
  // Even with int64, it can be used only when BIGINT_BASE<=131072 (limited by the precision of double)
  // This guarantees that the estimated result q' and actual result q satisfy q'<=q<=q'+1
  // Thus, on average, only one trial division is needed for each digit; the carries can be handled uniformly afterward
  // If a larger base is to be used, another trial-quotient method is needed
  double t = (b.get((unsigned)b.v.size() - 2) +
              (b.get((unsigned)b.v.size() - 3) + 1.0) / BIGINT_BASE);
  double db = 1.0 / (b.v.back() + t / BIGINT_BASE);
  for (size_t i = v.size() - 1, j = d.v.size() - 1; j <= v.size();) {
    int rm = r.get(i + 1) * BIGINT_BASE + r.get(i);
    int m = std::max((int)(db * rm), r.get(i + 1));
    r.sub_mul(b, m, j);
    d.v[j] += m;
    if (!r.get(i + 1))  // Check whether the highest position is already 0 to avoid extreme cases
      --i, --j;
  }
  r.trim();
  // Correct the units digit of the result
  int carry = 0;
  while (!r.absless(b)) {
    r.subtract(b);
    ++carry;
  }
  // Correct the carry at each position
  for (size_t i = 0; i < d.v.size(); ++i) {
    carry += d.v[i];
    d.v[i] = carry % BIGINT_BASE;
    carry /= BIGINT_BASE;
  }
  d.trim();
  d.sign = sign * b.sign;
  return d;
}

BigIntSimple operator/(const BigIntSimple &b) const {
  BigIntSimple r;
  return div_mod(b, r);
}

BigIntSimple operator%(const BigIntSimple &b) const {
  BigIntSimple r;
  div_mod(b, r);
  return r;
}

Karatsuba Multiplication

Let the number of digits in a high-precision number be $n$. Then high-precision-by-high-precision column multiplication takes $O(n^2)$ time. This section introduces a more time-efficient algorithm proposed by the former Soviet (Russian) mathematician Anatoly Karatsuba, which is a divide-and-conquer algorithm.

Consider two decimal integers $x$ and $y$, each containing $n$ digits (leading zeroes are allowed). Choose any $0 < m < n$, and let

where $x_0, y_0, z_0, z_1 < 10^m$. We obtain

Observe that

so to calculate $z_1$, we only need to calculate $(x_1 + x_0) \cdot (y_1 + y_0)$ and then subtract $z_0$ and $z_2$.

The expression above is the core of the Karatsuba algorithm. It transforms a multiplication problem of length $n$ into $3$ smaller subproblems. If $m = \left\lceil \dfrac n 2 \right\rceil$ and the time taken by the Karatsuba algorithm to multiply two $n$-digit integers is $T(n)$, then $T(n) = 3 \cdot T \left(\left\lceil \dfrac n 2 \right\rceil\right) + O(n)$. By the master theorem, $T(n) = \Theta(n^{\log_2 3}) \approx \Theta(n^{1.585})$.

The entire process can be implemented recursively. For clarity, the following code implements polynomial multiplication using the Karatsuba algorithm and handles all carries at the end.

int *karatsuba_polymul(int n, int *a, int *b) {
  if (n <= 32) {
    // Calculate directly for small sizes to avoid the efficiency loss from continuing recursion
    int *r = new int[n * 2 + 1]();
    for (int i = 0; i <= n; ++i)
      for (int j = 0; j <= n; ++j) r[i + j] += a[i] * b[j];
    return r;
  }

  int m = n / 2 + 1;
  int *r = new int[m * 4 + 1]();
  int *z0, *z1, *z2;

  z0 = karatsuba_polymul(m - 1, a, b);
  z2 = karatsuba_polymul(n - m, a + m, b + m);

  // Calculate z1
  // Temporarily modify and restore after calculation
  for (int i = 0; i + m <= n; ++i) a[i] += a[i + m];
  for (int i = 0; i + m <= n; ++i) b[i] += b[i + m];
  z1 = karatsuba_polymul(m - 1, a, b);
  for (int i = 0; i + m <= n; ++i) a[i] -= a[i + m];
  for (int i = 0; i + m <= n; ++i) b[i] -= b[i + m];
  for (int i = 0; i <= (m - 1) * 2; ++i) z1[i] -= z0[i];
  for (int i = 0; i <= (n - m) * 2; ++i) z1[i] -= z2[i];

  // Combine z0, z1, and z2 to obtain the result
  for (int i = 0; i <= (m - 1) * 2; ++i) r[i] += z0[i];
  for (int i = 0; i <= (m - 1) * 2; ++i) r[i + m] += z1[i];
  for (int i = 0; i <= (n - m) * 2; ++i) r[i + m * 2] += z2[i];

  delete[] z0;
  delete[] z1;
  delete[] z2;
  return r;
}

void karatsuba_mul(int a[], int b[], int c[]) {
  int *r = karatsuba_polymul(LEN - 1, a, b);
  memcpy(c, r, sizeof(int) * LEN);
  for (int i = 0; i < LEN - 1; ++i)
    if (c[i] >= 10) {
      c[i + 1] += c[i] / 10;
      c[i] %= 10;
    }
  delete[] r;
}

However, this implementation has a problem: in base $b$, each coefficient of the polynomial may reach the order of $n \cdot b^2$, which may cause integer overflow in a packed-digit high-precision implementation. If carries are handled during polynomial multiplication, however, the results of $x_1 + x_0$ and $y_1 + y_0$ may reach $2 \cdot b^m$, adding one digit (if $x_1 - x_0$ is used instead, negative values must be handled specially). Therefore, the implementation should be chosen according to the actual application scenario.

Efficient Large-Integer Multiplication Based on Polynomials

If the data size reaches $10^{10^5}$ or larger, ordinary high-precision multiplication may time out. This section introduces polynomial-based methods for optimizing this type of multiplication.

For an $n$-digit decimal integer $a$, regard it as a polynomial whose coefficients are integers no greater than $10$: $A=a_{0} 10^0+a_{1} 10^1+\cdots+a_{n-1} 10^{n-1}$. In this way, multiplication of two integers is transformed into multiplication of two polynomials.

The time complexity of ordinary polynomial multiplication is still $O(n^2)$, but it can be optimized using algorithms such as the Fast Fourier Transform and the Number-Theoretic Transform from the polynomial section. The optimized time complexity is $O(n\log n)$.

Wrapper Classes

Here is a packaged high-precision integer class, and here is an ultra-mini implementation class supporting dynamic length and the four arithmetic operations.

constexpr int MAXN = 9999;
// MAXN is the largest digit in one position
constexpr int MAXSIZE = 10024;
// MAXSIZE is the number of positions
constexpr int DLEN = 4;

// DLEN records how many digits are packed
struct Big {
  int a[MAXSIZE], len;
  bool flag;  // Mark the '-' sign

  Big() {
    len = 1;
    memset(a, 0, sizeof a);
    flag = false;
  }

  Big(const int);
  Big(const char*);
  Big(const Big&);
  Big& operator=(const Big&);
  Big operator+(const Big&) const;
  Big operator-(const Big&) const;
  Big operator*(const Big&) const;
  Big operator/(const int&) const;
  // TODO: Big / Big;
  Big operator^(const int&) const;
  // TODO: Big ^ Big;

  // TODO: Big bit operations;

  int operator%(const int&) const;
  // TODO: Big ^ Big;
  bool operator<(const Big&) const;
  bool operator<(const int& t) const;
  void print() const;
};

Big::Big(const int b) {
  int c, d = b;
  len = 0;
  // memset(a,0,sizeof a);
  CLR(a);
  while (d > MAXN) {
    c = d - (d / (MAXN + 1) * (MAXN + 1));
    d = d / (MAXN + 1);
    a[len++] = c;
  }
  a[len++] = d;
}

Big::Big(const char* s) {
  int t, k, index, l;
  CLR(a);
  l = strlen(s);
  len = l / DLEN;
  if (l % DLEN) ++len;
  index = 0;
  for (int i = l - 1; i >= 0; i -= DLEN) {
    t = 0;
    k = i - DLEN + 1;
    if (k < 0) k = 0;
    g(j, k, i) t = t * 10 + s[j] - '0';
    a[index++] = t;
  }
}

Big::Big(const Big& T) : len(T.len) {
  CLR(a);
  f(i, 0, len) a[i] = T.a[i];
  // TODO: Overload here?
}

Big& Big::operator=(const Big& T) {
  CLR(a);
  len = T.len;
  f(i, 0, len) a[i] = T.a[i];
  return *this;
}

Big Big::operator+(const Big& T) const {
  Big t(*this);
  int big = len;
  if (T.len > len) big = T.len;
  f(i, 0, big) {
    t.a[i] += T.a[i];
    if (t.a[i] > MAXN) {
      ++t.a[i + 1];
      t.a[i] -= MAXN + 1;
    }
  }
  if (t.a[big])
    t.len = big + 1;
  else
    t.len = big;
  return t;
}

Big Big::operator-(const Big& T) const {
  int big;
  bool ctf;
  Big t1, t2;
  if (*this < T) {
    t1 = T;
    t2 = *this;
    ctf = true;
  } else {
    t1 = *this;
    t2 = T;
    ctf = false;
  }
  big = t1.len;
  int j = 0;
  f(i, 0, big) {
    if (t1.a[i] < t2.a[i]) {
      j = i + 1;
      while (t1.a[j] == 0) ++j;
      --t1.a[j--];
      // WTF?
      while (j > i) t1.a[j--] += MAXN;
      t1.a[i] += MAXN + 1 - t2.a[i];
    } else
      t1.a[i] -= t2.a[i];
  }
  t1.len = big;
  while (t1.len > 1 && t1.a[t1.len - 1] == 0) {
    --t1.len;
    --big;
  }
  if (ctf) t1.a[big - 1] = -t1.a[big - 1];
  return t1;
}

Big Big::operator*(const Big& T) const {
  Big res;
  int up;
  int te, tee;
  f(i, 0, len) {
    up = 0;
    f(j, 0, T.len) {
      te = a[i] * T.a[j] + res.a[i + j] + up;
      if (te > MAXN) {
        tee = te - te / (MAXN + 1) * (MAXN + 1);
        up = te / (MAXN + 1);
        res.a[i + j] = tee;
      } else {
        up = 0;
        res.a[i + j] = te;
      }
    }
    if (up) res.a[i + T.len] = up;
  }
  res.len = len + T.len;
  while (res.len > 1 && res.a[res.len - 1] == 0) --res.len;
  return res;
}

Big Big::operator/(const int& b) const {
  Big res;
  int down = 0;
  gd(i, len - 1, 0) {
    res.a[i] = (a[i] + down * (MAXN + 1)) / b;
    down = a[i] + down * (MAXN + 1) - res.a[i] * b;
  }
  res.len = len;
  while (res.len > 1 && res.a[res.len - 1] == 0) --res.len;
  return res;
}

int Big::operator%(const int& b) const {
  int d = 0;
  gd(i, len - 1, 0) d = (d * (MAXN + 1) % b + a[i]) % b;
  return d;
}

Big Big::operator^(const int& n) const {
  Big t(n), res(1);
  int y = n;
  while (y) {
    if (y & 1) res = res * t;
    t = t * t;
    y >>= 1;
  }
  return res;
}

bool Big::operator<(const Big& T) const {
  int ln;
  if (len < T.len) return true;
  if (len == T.len) {
    ln = len - 1;
    while (ln >= 0 && a[ln] == T.a[ln]) --ln;
    if (ln >= 0 && a[ln] < T.a[ln]) return true;
    return false;
  }
  return false;
}

bool Big::operator<(const int& t) const {
  Big tee(t);
  return *this < tee;
}

void Big::print() const {
  printf("%d", a[len - 1]);
  gd(i, len - 2, 0) { printf("%04d", a[i]); }
}

void print(const Big& s) {
  int len = s.len;
  printf("%d", s.a[len - 1]);
  gd(i, len - 2, 0) { printf("%04d", s.a[i]); }
}

char s[100024];

Exercises

• NOIP 2012 King Game
• SPOJ - Fast Multiplication
• SPOJ - GCD2
• UVa - Division
• UVa - Fibonacci Freeze
• Codeforces - Notepad

References and Links

1. Karatsuba algorithm - Wikipedia