extended_precision_arith.m

% extending integer precision beyond 32-bit. Maximum integer value in
% Matlab is 'intmax = 2147483647'. This is a 'int' type in C. So if using
% 'intmax', double() must be used to convert it to double type in order to
% perform further computation. Others are by default double so no need of
% conversion.
% Reference: C Interfaces and Implmentations, David Hanson.

clear all

% must turn intmax to double to further computation. 
% intmax is int32 class by default. All others are double by default
% u and base define how a number is represented. For example,
% 2147483647_10 = [255 255 255 127]_256, where "_" stands for radix or
% base.
u = double(intmax);
base = 2^8;


% Step 1. convert integer to array form in 'base'
temp = 0;
i = 1;

% store u into an array
while ( u > 0)
    temp = mod(u, base);
    a(i) = temp;
    u = (u - temp)/base;
    i = i + 1;
end

% Step 2. return an array into an integer form
len = length(a);
v = 0;
for i = 1 : len
    v = v + base ^ (i - 1) * a(i);
end

% just for sanity check
diff = intmax - v;

% Arith1: Addition
carry = 0;

% addition. All operators should have the same length. Padding 0's if
% needed. In this example, we have 2147483647 + 255 = 2147483902, which is
% [254 0 0 128] = 254 + 128*256^3.
b = [255 0 0 0];

for i = 1 : length(a)
    tmp_add = carry + a(i) + b(i);
    c(i) = mod(tmp_add, base);
    carry = floor(tmp_add / base);
end

% Arith2: substraction
borrow = 0;
c = [255 0 0 0];

for i = 1 : length(a)
    tmp_sub = (a(i) + base) - borrow - c(i);
    d(i) = mod(tmp_sub, base);
    borrow = floor(1 - d(i)/base);
end

% Arith3: multiplication
e = [255 255 0 0];
carry = 0;

z = zeros(1, length(a)+length(e));

for i = 1 : length(a)
    for j = 1 : length(e)
        carry = carry + a(i) * e(j) + z(i+j);
        z(i+j) = mod(carry, base);
        carry = floor(carry/base);
    end
    for j = length(e) : length(a) + length(e) - i
        carry = carry + z(i+j);
        z(i+j) = mod(carry, base);
        carry = floor(carry/base);
    end
end

z = z(2:end);