```// ************************************************************************** //
//                                                                            //
//    eses                   eses                                             //
//   eses                     eses                                            //
//  eses    eseses  esesese    eses   Embedded Systems Group                  //
//  ese    ese  ese ese         ese                                           //
//  ese    eseseses eseseses    ese   Department of Computer Science          //
//  eses   eses          ese   eses                                           //
//   eses   eseses  eseseses  eses    University of Kaiserslautern            //
//    eses                   eses                                             //
//                                                                            //
// ************************************************************************** //
// This file contains functions for radix-B arithmetic. A radix-B number is   //
// thereby represented as a []nat array where B=2^n, e.g., B=2^16 on Abacus   //
// with 16-bit words.                                                         //
// ************************************************************************** //
bool correct;

// -----------------------------------------------------------------------------
// Radix-B Addition y[0..len] = x1[0..len-1] + x2[0..len-1]
// -----------------------------------------------------------------------------
// For the correctness of the computation, note that we should compute c' and
// y[i] so that c'*B+y[i] = x1[i] + x2[i] + c holds. It is clear by the code
// that g*B + s = x1[i] + x2[i] and p*B + y[i] = s + c hold, so that we also
// have (g+p)*B + y[i] = x1[i] + x2[i] + c. Thus, we could determine c' simply
// as c' = g+p, but we can alternatively also define it as bitwise disjunction
// or bitwise xor, since only one of p,g can be 1 as proved below.
//  We prove that g,p,c∈{0,1} and that g=1 implies p=0:
// Since 0 <= x1[i]+x2[i] <= B+(B-2) holds, and B+(B-2) is the machine double
// word |0..01|1..101|, it follows that g is one of the machine words |0..0| or
// |0..01|, and if g=|0..01|, then s!=|1..1|. If c is also one of |0..0| or
// |0..01|, it follows that 0 <= s+c <= B, and the latter is the machine double
// word |0..01|0..0|. Hence, p is also either one of the machine words |0..0| or
// |0..01|, and the same follows now for c which is obtained as bitwise
// disjunction of words that are either |0..0| or |0..01|.
//   Finally, if g=|0..01|, then s!=|1..1| (as noted above), and therefore when
// adding c which is either |0..0| or |0..01|, the sum will not require
// additional bits so that p=|0..0| holds in that case. So, g and p cannot be
// both |1..1|, and thus all assignments c' = g+p, c' = g|p or c' = !(q<->p)
// are equivalent to each other.
// -----------------------------------------------------------------------------

procedure NatAddSameLen([]nat x1,x2,y,nat len) {
nat c,g,p,s,i;
c = 0;                      // initial carry digit
for(i=0..len-1) {           // ripple from least to highest digit
g,s = x1[i] + x2[i];    // g:generate carry; s:preliminary sum
p,y[i] = s + c;         // p:propagate carry; add previous carry
c = g+p;                // define next carry digit 0 or 1
}
y[len] = c;
return;
}

// -----------------------------------------------------------------------------
// Radix-B Subtraction y[0..len] = x1[0..len-1] - x2[0..len-1]
// -----------------------------------------------------------------------------
// For the correctness of the computation, note that we should compute c' and
// y[i] so that -c'*B+y[i] = x1[i] - x2[i] + c holds with c',c∈{0,1}. It is
// clear by the code that g*B + s = x1[i] - x2[i] and p*B + y[i] = s - c hold,
// so that we also have (g+p)*B + y[i] = x1[i] - x2[i] - c. Thus, we could
// determine c' simply as c' = -(g+p), but we have to consider that all this
// will be correct in case of overflows as well.
//  To this end, we will see that g,p are either |0..0| or |1..1| to be read as
// two-complement numbers 0 and -1, respectively. To compute x1[i] - x2[i] with
// an unsigned subtraction, both numbers x1[i] and x2[i] are extended by |0..0|
// to perform the subtraction in full precision.
//   Since -(B-1) <= x1[i]-x2[i] <= +(B-1) holds, and -(B-1) and +(B-1) are the
// machine double words |1..1|0..01| and |0..0|1..1|, respectively, it follows
// that g is one of the machine words |1..1| or |0..0|, and if g=|1..1|, then
// s!=|0..0|. Note that we should read the words here as 2-complement numbers,
// so that g is either -1 or 0, which means that we have to subtract B or not,
// respectively, in the next step via c.
//   If c∈{0,1}, it follows that -1 <= s-c <= B-1, and the latter is the machine
// double word |0..0|1..1| while -1 is |1..1|1..1|. Hence, p is also either one
// of the machine words |0..0| or |1..1|. Finally, note that since g=|1..1|,
// then s!=|0..0|, and therefore, 0<=s-c<=B-1, so that p=|0..0| follows. Hence,
// g and p cannot be both |1..1|, so that g+p, g|p and !(g<->p) are all
// equivalent and will either produce |0..0| or |1..1|, meaning the integers 0
// and -1, respectively. Changing their signs by subtracting from 0 yields
// finally c'∈{0,1}.
// -----------------------------------------------------------------------------

procedure NatSubSameLen([]nat x1,x2,y,nat len) {
nat c,g,p,s,i;
c = 0;                      // initial carry digit ∈{0,1}
for(i=0..len-1) {           // ripple from least to highest digit
g,s = x1[i] - x2[i];    // g:generate carry; p∈{|0..0|,|1..1|}
p,y[i] = s - c;         // p:propagate carry; p∈{|0..0|,|1..1|}
c = 0-(g+p);            // g+p∈{|0..0|,|1..1|}, thus c∈{0,1}
}
y[len] = g+p;               // non-negative iff y[len]=|0..0|
return;                     // negative iff y[len]=|1..1|
}

// -----------------------------------------------------------------------------
// Radix-B Multiplication y[0..len1+len2-1] = x1[0..len1-1] * x2[0..len2-1]
// -----------------------------------------------------------------------------
// For the correctness, note that the maximum of pr[i+j] + x1[i] * x2[j] + c is
// (B-1)*B + (B-1) so that the maximum of the sum of the carry digits computed
// in the algorithm is B-1 and therefore their addition does not overflow.
// -----------------------------------------------------------------------------

procedure NatMul([]nat x1,x2,pr,nat len1,len2) {
nat c1,c,p,i,j,k;
for(j=0..len2-1) {
// compute pr[j+len1..j] = pr[j+len1-1..j] + x1[len1-1..0] * x2[j]
c = 0;
for(i=0..len1-1) {
c1,p = x1[i] * x2[j];
c,p = p + c;
c = c + c1;
k = i+j;
c1,pr[k] = p + pr[k];
c = c + c1;
}
pr[len1+j] = c;
}
return;
}

// -----------------------------------------------------------------------------
// Radix-B Equality checking y = x1[0..len1-1] = x2[0..len2-1]
// -----------------------------------------------------------------------------

function NatEquSameLen([]nat x1,x2,nat len) : bool {
nat i;
bool eq;
i = len-1;
eq = true;
while(i>0 & eq) {
eq = eq & x1[i]==x2[i];
i = i-1;
}
return eq;
}

// -----------------------------------------------------------------------------
// Radix-B Comparison y = x1[0..len1-1] ? x2[0..len2-1] yields -1 if x1<x2,
// 0 if x1==x2 and +1 if x1>x2
// -----------------------------------------------------------------------------

function NatCmpSameLen([]nat x1,x2,nat len) : int {
nat i;
bool eq;
i = len;
eq = true;
while(i>0 & eq) {
i = i-1;
eq = eq & x1[i]==x2[i];
}
if(x1[i]==x2[i]) return +0;
else if(x1[i]<x2[i]) return -1;
else return +1;
}

// -----------------------------------------------------------------------------
// compute sum digits y[0..len] = x1[0..len-1] + x2[0..len-1]
// -----------------------------------------------------------------------------

function NatAddCLA ([]nat x1,x2,y,c, []bool g,p, nat len): []nat {
nat i,j,s,s1,len1;  // loop variables
len1 = len-1;       // compute this constant once and forall

// first row of generate and propagate signals
for(i=0..len1) {
c[i],y[i] = x1[i] + x2[i];  // preliminary sum and carry digits
g[i] = c[i]>0;              // generate signal at level 0
p[i] = y[i]==(nat)(-1);     // propagate signal at level 0
}

// compute generate/propagate signals by forward traversal on lookahead tree
s = 1;
while (s<len) {
i = s-1;                        // middle index of considered subtree
s = s+s;                        // width of current subtrees s=2^l
j = s-1;                        // leftmost index of considered subtree
if(j>len1)
j = len1;
while(i<len1) {
g[j] = g[j] | g[i] & p[j];
p[j] = p[j] & p[i];
i = i+s;
j = j+s;
if(j>len1) j = len1;
}
}

// compute carry signals c[i] by backward traversal of the lookahead tree
c[0] = 0;                                           // rightmost carry signal
c[len] = (nat) (g[len1] | p[len1] & (bool)c[0]);    // leftmost carry signal
while(s>1) {
s1 = s/2;
i = s1;
j = 0;
while(i<len) {
c[i] = (nat) (g[i-1] | p[i-1] & (bool)c[j]);
i = i+s;
j = j+s;
}
s = s1;
}

// compute the final sum digits with the carry and preliminary sum digits
// note that we should have c[i] as either 1 or 0, but due to the casts from
// boolean values to type nat, it becomes MaxNat and subtracting this is
// the same as subtracting -1 (if we ignore overflows)
for(i=1..len) {
y[i] = y[i] - c[i];
}

return y;
}

// ----------------------------------------------------------------------------
// NatAdd implements carry lookahead addition of radix-B numbers x1 and x2 which
// both have length len. The sum is stored in array y whose length is len+1. The
// carry digits are computed with a prefix tree with a depth of log(len) so that
// the procedure could be computed in time log(len) using O(len) processing
// elements.
// ----------------------------------------------------------------------------

function NatAddCLA ([]nat x1,x2,[]nat y,c, []bool g,p, nat len): []nat {
nat i,j,s,s1,len1;  // loop variables
len1 = len-1;       // compute this constant once and forall

// first row of generate and propagate signals
for(i=0..len1) {
c[i],y[i] = x1[i] + x2[i];  // preliminary sum and carry digits
g[i] = c[i]>0;              // generate signal at level 0
p[i] = y[i]==(nat)(-1);     // propagate signal at level 0
}

// compute generate/propagate signals by forward traversal on lookahead tree
s = 1;
while (s<len) {
i = s-1;                        // middle index of considered subtree
s = s+s;                        // width of current subtrees s=2^l
j = s-1;                        // leftmost index of considered subtree
if(j>len1)
j = len1;
while(i<len1) {
g[j] = g[j] | g[i] & p[j];
p[j] = p[j] & p[i];
i = i+s;
j = j+s;
if(j>len1) j = len1;
}
}

// compute carry signals c[i] by backward traversal of the lookahead tree
c[0] = 0;                                           // rightmost carry signal
c[len] = (nat) (g[len1] | p[len1] & (bool)c[0]);    // leftmost carry signal
while(s>1) {
s1 = s/2;
i = s1;
j = 0;
while(i<len) {
c[i] = (nat) (g[i-1] | p[i-1] & (bool)c[j]);
i = i+s;
j = j+s;
}
s = s1;
}

// compute the final sum digits with the carry and preliminary sum digits
// note that we should have c[i] as either 1 or 0, but due to the casts from
// boolean values to type nat, it becomes MaxNat and subtracting this is
// the same as subtracting -1 (if we ignore overflows)
for(i=1..len) {
y[i] = y[i] - c[i];
}

return y;
}

// -----------------------------------------------------------------------------
// compute sum digits y[0..len] = x1[0..len1-1] + x2[0..len2-1]
// -----------------------------------------------------------------------------

procedure NatAdd([]nat x1,x2,sum,nat len1,len2) {
nat c,g,p,s,i,lmin;
// determine set of common digits
if(len1<len2) lmin = len1; else lmin = len2;
// carry ripple addition for common digits
c = 0;
for(i=0..lmin-1) {
g,s = x1[i] + x2[i];
p,sum[i] = s + c;
c = (nat) ((bool) g | (bool) p);
}
// propagate c through the final digits of the longer number
if(len1<len2) {
for(i=lmin..len2-1)
c,sum[i] = x2[i] + c;
sum[len2] = c;
} else {
for(i=lmin..len1-1)
c,sum[i] = x1[i] + c;
sum[len1] = c;
}
return;
}

thread test {
[8]nat a;
[6]nat b;
[9]nat sm;
nat i;
for(i=0..5)
a[i] = (nat) -1;
b[0] = 1;
NatAdd(a,b,sm,8,6);
}
```