# Arithmetic
Here are some simple implementations of arbitary precision arithmetic. For complicated but robust and fast implementations, see libraries like GMP, MPFR, MPC, etc.
The implementations on this page have no bounds checking, a proper version would refuse to proceed if there is not enough space in the output buffers.
# 1 N
Natural numbers can be stored as null-terminated strings of decimal digits in ASCII or similar character encoding. Zero is represented as the empty string. This is not the most efficient use of storage space or CPU time (due to divisions by 10 when handling carry when adding), but does make decimal input and output very simple.
For simplicity of implementation, little endian digit order is chosen.
# 1.1 Print
To print little endian digits in conventional big endian order, the digits have to be reversed:
void print(const char *a)
{
if (*a)
{
const char *b = a;
while (*b)
{
++b;
}
while (b != a)
{
putchar(*(--b));
}
}
else
{
putchar('0');
}
}
# 1.2 Set
Setting a variable to another variable is a string copy.
void set(char *out, const char *a)
{
while (*a)
{
*out++ = *a++;
}
*out = 0;
}
It would be useful to have a function for setting a natural string from an int variable, but I didn’t make it yet. It would do repeated divisions by 10, storing the remainders.
# 1.3 Add
Addition zips the inputs together, propagating carry (exercise: prove correctness, considering the range of values the carry variable can take). The numbers can have different numbers of digits; at most 2 of the 3 while loops are entered. Don’t forget the final carry out.
void add(char *out, const char *a, const char *b)
{
char carry = 0;
while (*a && *b)
{
char sum = (*a++ - '0') + (*b++ - '0') + carry;
carry = sum / 10;
*out++ = (sum % 10) + '0';
}
while (*a)
{
char sum = (*a++ - '0') + carry;
carry = sum / 10;
*out++ = (sum % 10) + '0';
}
while (*b)
{
char sum = (*b++ - '0') + carry;
carry = sum / 10;
*out++ = (sum % 10) + '0';
}
if (carry)
{
*out++ = carry + '0';
}
*out = 0;
}
This addition can work in place, that is,
out
can be the same buffer as a
, for example.
# 1.4 Mul
Multiplying by a single digit (stored as a plain number, not ASCII encoded!) can be done with repeated addition:
void mul1(char *out, const char *a, char b)
{
out[0] = 0;
while (b--)
{
add(out, out, a);
}
}
Using this algorithm for big number multiplication would be very slow. At the cost of more temporary space (which might be avoided with more code complexity, by fusing multiplication with in-place addition) multiplication can be done a digit at a time. This algorithm is quadratic time in digit count, which is far from optimal, but the advantage is simplicity of implementation.
void mul(char *out, const char *a, const char *b, char *tmp)
{
out[0] = 0;
while (*b)
{
mul1(tmp, a, *b++ - '0');
add(out, out, tmp);
++out;
}
}
Unlike addition, this multiply does not work in place!
The out
buffer must not alias the a
or b
buffers.
# 1.5 Pow
Powers can be done by repeated multiplication, for small powers. For large powers algorithms like exponentiation by squaring are much faster.
void pow(char *out, const char *a, char b, char *tmp1, char *tmp2)
{
out[0] = '1';
out[1] = 0;
while (b--)
{
set(tmp1, out);
mul(out, tmp1, a, tmp2);
}
}
# 1.6 Example
For example, to calculate 2375 to the power of 15, noting that the buffers won’t overflow because \(2375 < 10^4\) and \(4 \times 15 = 60 < 255\).
int main()
{
#define DIGITS 255
char result[DIGITS+1], tmp1[DIGITS+1], tmp2[DIGITS+1];
const char base[] = "5732"; // reversed digits of 2375
char power = 15;
pow(result, base, power, tmp1, tmp2);
print(result);
putchar('\n');
return 0;
}
Expected output (matches bc
):
431473581269153734723431625752709805965423583984375
This code successfully cross-compiles for Commodore C64 using cc65
and produces the correct answer in that environment,
taking a few seconds to run.