function z = bigintadd(x,y,b)
% BIGINTADD - adds x and y to give z, in arbitrary-precision ints, base b
%
% function z = bigintadd(x,y,b).  x, y, z are row-vectors, b is integer
% Note, b should be smaller than half the largest possible integer
% exactly represented by Matlab's floating point doubles, ie 2^53.

nx = numel(x); ny = numel(y);
n = max(nx,ny) + 1;    % longest poss output array
x = [zeros(1,n-nx) x]; % zero pad at the start
y = [zeros(1,n-ny) y];
z = x+y;             % do elementwise add
for i=n:-1:2         % loop from LSB to MSB
  c = floor(z(i)/b); % carry over (0 if no carry)
  z(i) = z(i) - c*b;  % correct if carried
  z(i-1) = z(i-1) + c; % add the carry
end
if z(1)==0, z = z(2:end); end  % remove leading zero