{Author : Hans Riesel
 See "Prime Numbers and Computer Methods for
 Factorization" by Hans Riesel second edition
 pages 349 - 355}
program multipleprecisionpackage(input, output);
{tests the procedures in the package}
uses wincrt;
label 1, 2;
const sqrtb = 10; vsize = 20; lvsize = 40;
type vector = array[0..vsize] of integer;
     alpha = array[1..8] of char;
var base, b1, b101 : integer;
    a, b, c, d : vector;

function sign(x : integer) : integer;
  begin {sign}
    if x = 0 then sign := 0 else sign := abs(x) div x;
  end {sign};

procedure putin(txt : alpha; var a : vector);
{prints a text string of 8 characters and read a multiple-
 precision integer which is stored in a}
  var a0, m, i : integer;

  begin {putin}
    write(txt); read(a0); m := abs(a0); a[0] := a0;
    for i := m downto 1 do read(a[i])
  end {putin};

procedure putout(txt : alpha; var a : vector);
{prints a text string of 8 characters followed by one
 multiple-precision integer a}
  const l = 0.434294481 {const = 1 / ln(10)};
  var a0, i, j, s, sa, t, u, v : integer;

  begin {putout}
    write(txt); a0 := a[0]; sa := sign(a0); a0 := abs(a0);
    if sa < 0 then write('-') else write(' ');
    v := trunc(l * ln(base) + 0.5); {v is the number of
    digits allowed in each component of a}
    if a0 = 0 then write(0 : v)
    else
      begin {else}
        write(a[a0] : v);
        for i := a0 - 1 downto 1 do
          begin {for}
            t := a[i];
            if t = 0 then s := 1
            else if t = b101 then s := v - 1
            else if t >= base div 10 then s:= v
            else s := trunc(1 + l * ln(t));
            u := v - s;
            write(' ');
            for j := 1 to u do write('0');
            write(t : s)
          end {for}
      end {else};
    writeln;
  end {putout};

procedure red(var c : vector);
{reduces a multiple=precision integer to standard form
 in call of procedure addsub}
  var c0, i, j, s : integer;

  begin {red}
    c0 := c[0]; s := sign(c0); c0 := abs(c0); c[c0 + 1] := 0;
    while (c0 > 0) and (c[c0] = 0) do c0 := c0 - 1;
    {here the leading zeros of c are removed}
    if c[c0] < 0 then for i := 0 to c0 do c[i] := - c[i];
    for j := 1 to 2 do
      begin {for}
        for i := 1 to c0 do
          begin {for}
            if c[i] >= base then
              begin {if}
                c[i + 1] := c[i + 1] + 1;
                c[i] := c[i] - base
              end {if}
            else if c[i] < 0 then
              begin {if}
                c[i + 1] := c[i + 1] - 1;
                c[i] := c[i] + base
              end {if}
          end {for}
      end {for};
    c0 := c0 + 1;
    while (c0 > 0) and (c[c0] = 0) do c0 := c0 - 1;
    c[0] := sign(c[0]) * c0
  end {red};

procedure let(var a : vector; b : vector);
{puts a := b for multiple-precision integers a and b}
  var i, s : integer;

  begin {let}
    s := abs(b[0]);
    for i := 0 to s do a[i] := b[i];
  end {let};

procedure addsub(a, b : vector; sgn : integer; var c : vector);
{with a = sign(a[0]) * sum a[i] * B ^ (i - 1),
      b = sign(b[0]) * sum b[j] * B ^ (j - 1),
      and sgn = +1 or - 1,
      c := a + sgn * b, c may be one of a or b}
  var i, max, min, q, r, s, sa, sb : integer;

  begin {addsub}
    sa := a[0]; sb := b[0] * sgn; q := abs(sa); max := q;
    r := abs(sb); min := r;
    if min > max then begin max := r; min := q end;
    s := sign(sa) * sign(sb);
    for i := 1 to min do c[i] := a[i] + s * b[i];
    for i := min + 1 to max do
      if max = q then c[i] := a[i] else c[i] := s * b[i];
    c[0] := sa; c[max + 1] := 0;
    if sa = 0 then let(c, b); {c := b in case a = 0}
    if s <> 0 then c[0] := max * sign(sa);
    {here the components of c may be > B or < 0}
    red(c);
  end {addsub};

procedure mul(a, b : vector; var c : vector);
{computes c := a * b for multiple-precision integers
 c may be one of a or b}
  type longvector = array[0..lvsize] of integer;
  var i, k, m, m2, m21, mn, n, n2, n21, p, p2, s, sa, sb, sc :
      integer;
      aa, bb : vector;
      cc : longvector;

  begin {mul}
    m := a[0]; sa := sign(m); m := abs(m); n := b[0];
    sb := sign(n); n := abs(n); p := m + n; sc := sa * sb;
    m2 := m * 2; n2 := n * 2; p2 := p * 2; m21 := m2 - 1;
    n21 := n2 - 1; mn := m2 + n21;
    for i := 1 to m do
      begin {for}
        s := a[i]; aa[2 * i] := s div sqrtb;
        aa[2 * i - 1] := s mod sqrtb;
      end {for};
    for i := 1 to n do
      begin {for}
        s := b[i]; bb[2 * i] := s div sqrtb;
        bb[2 * i - 1] := s mod sqrtb;
      end {for};
    for i := 0 to lvsize do cc[i] := 0;
    if m <= n then
      begin {if}
        for i := 1 to m21 do
          begin {for}
            s := 0;
            for k := 1 to i do
              begin {for}
                s := s + aa[k] * bb[i + 1 - k];
                if s >= base then
                  begin {if}
                    s := s - base;
                    cc[i + 2] := cc[i + 2] + 1
                  end {if}
              end {for};
            cc[i] := cc[i] + s
          end {for};
        for i := m2 to n2 do
          begin {for}
            s := 0;
            for k := 1 to m2 do
              begin {for}
                s := s + aa[k] * bb[i + 1 - k];
                if s >= base then
                  begin {if}
                    s := s - base;
                    cc[i + 2] := cc[i + 2] + 1
                  end {if}
              end {for};
            cc[i] := cc[i] + s
          end {for};
        for i := n2 + 1 to mn do
          begin {for}
            s := 0;
            for k := i - n21 to m2 do
              begin {for}
                s := s + aa[k] * bb[i + 1 - k];
                if s >= base then
                  begin {if}
                    s := s - base;
                    cc[i + 2] := cc[i + 2] + 1
                  end {if}
              end {for};
            cc[i] := cc[i] + s
          end {for}
      end {if}
    else
      begin {else}
        for i := 1 to n21 do
          begin {for}
            s := 0;
            for k := 1 to i do
              begin {for}
                s := s + bb[k] * aa[i + 1 - k];
                if s >= base then
                  begin {if}
                    s := s - base;
                    cc[i + 2] := cc[i + 2] + 1
                  end {if}
              end {for};
            cc[i] := cc[i] + s
          end {for};
        for i := n2 to m2 do
          begin {for}
            s := 0;
            for k := 1 to n2 do
              begin {for}
                s := s + bb[k] * aa[i + 1 - k];
                if s >= base then
                  begin {if}
                    s := s - base;
                    cc[i + 2] := cc[i + 2] + 1
                  end {if}
              end {for};
            cc[i] := cc[i] + s
          end {for};
        for i := m2 + 1 to mn do
          begin {for}
            s := 0;
            for k := i - m21 to n2 do
              begin {for}
                s := s + bb[k] * aa[i + 1 - k];
                if s >= base then
                  begin {if}
                    s := s - base;
                    cc[i + 2] := cc[i + 2] + 1
                  end {if}
              end {for};
            cc[i] := cc[i] + s
          end {for}
      end {else};
    for i := 1 to p2 do
      begin {for}
        cc[i + 1] := cc[i + 1] + cc[i] div sqrtb;
        cc[i] := cc[i] mod sqrtb
      end {for};
    for i := 1 to p do
      c[i] := cc[2 * i - 1] + sqrtb * cc[2 * i];
    if c[p] = 0 then p := p - 1; c[0] := p * sc
  end {mul};

procedure quot(a, b : vector; var q, r : vector);
{computes q and r satisfying a = b * q + r, 0 <= r < b
 for multiple-precision integers a and b, r may be
 either of a or b}
  label 1, 2, 3;
  var i, j, j1, k, l, m, n, p, s, sa, sb, t, t1, u : integer;
      ar, br, qr : real;
      aa, bb, cc, one, qk, qq : vector;
       
  begin {quot}
    for i := 0 to vsize do qq[i] := 0;
    m := a[0]; n := b[0];
    sa := sign(m); sb := sign(n); s := abs(m); t := abs(n);
    l := s - t + 1;
    if l <= 0 then l := 1;
    if s < t then u := t else u := s;
    if n = 0 then
      begin {if}
        writeln('division by zero attempted in quot');
        l := 0;
        goto 3
      end {if};
    one[0] := 1; one[1] := 1; qk[0] := 1; t1 := t + 1;
    b[0] := t;
    let(aa, a);
    aa[0] := s;
    if t = 1 then br := b[1] else br := b[t] + b[t - 1] / base;
    {here br is the leading part of b}
    1 : if s < t then goto 3;
    if s = t then
      begin {if}
        addsub(aa, b, - 1, cc);
        if cc[0] < 0 then goto 3
      end {if};
    j := aa[0] - t;
    if s <= 1 then ar := aa[s] else ar := aa[s] + aa[s- 1] / base;
    {here ar is the leading part of a}
    qr := ar / br;
    {qr is the approximate quotient of a / b}
    if (qr < 1) and (j = 0) then goto 3
    else if qr < 1 then
      begin {if}
        qr := qr * base; j := j - 1;
        if qr < 1 then qr := 1
      end {if};
    k := trunc(qr); qk[1] := k; qq[j + 1] := qq[j + 1] + k;
    2 : mul(qk, b, bb);
    for i := t1 downto 1 do bb[i + j] := bb[i];
    for i := 1 to j do bb[i] := 0;
    bb[0] := bb[0] + j;
    addsub(aa, bb, - 1, cc);
    s := cc[0]; p := s;
    if p < 0 then
      begin {if}
        qk[1] := qk[1] - 1; qq[j + 1] := qq[j + 1] - 1;
        if qk[1] = 0 then
          begin {if}
            qk[1] := b1; qq[j] := qq[j] + qk[1]; j := j - 1
          end {if};
        goto 2
      end {if};
    for i := 0 to p do aa[i] := cc[i];
    goto 1;
    3 : b[0] := n;
    if qq[l] = 0 then l := l - 1;
    qq[0] := l * sa * sb;
    let(q, qq);
    if l = 0 then aa[0] := m else aa[0] := aa[0] * sa * sb;
    if (sa * sb < 0) and (aa[0] <> 0) then
      begin {if}
        addsub(q, one, - 1, q);
        addsub(aa, b, sb, aa)
      end {if};
    let(r, aa);
  end {quot};

procedure euclid(a, b : vector; var d : vector);
{computes d = gcd(a, b) for multiple-precision integers
 with Euclid's algorithm}
  var a0, aa0, b0, bb0, i, r0 : integer;
      aa, bb, q, r : vector;

  begin {euclid}
    a0 := abs(a[0]); b0 := abs(b[0]);
    let(aa, a);
    aa[0] := a0;
    let(bb, b);
    bb[0] := b0;
    while bb[0] > 0 do
      begin {while}
        quot(aa, bb, q, r);
        let(aa, bb);
        let(bb, r)
      end {while};
    let(d, aa)
  end {euclid};

 procedure apowerb(a, b, n : vector; var r : vector);
 {computes r = a ^ b (mod n) for multiple-precision
  integers, b > 0}
   label 1;
   var aa, bb, gb, p, par, s, two : vector;

   begin {apowerb}
     p[0] := 1; p[1] := 1; two[0] := 1; two[1] := 2;
     let(bb, b);
     if bb[0] <= 0 then
       begin {if}
         writeln('b <= 0 attempted in apowerb');
         goto 1
       end {if};
     quot(a, n, gb, aa);
     while bb[0] > 0 do
       begin {while}
         quot(bb, two, gb, par);
         let(bb, gb);
         if par[0] = 1 then
           begin {if}
             mul(aa, p, s);
             quot(s, n, gb, p)
           end {if};
         mul(aa, aa, s);
         quot(s, n, gb, aa)
       end {while};
     let(r, p);
     1:
   end {apowerb};

begin {multipleprecisionpackage}
  base := sqr(sqrtb); b1 := base - 1; b101 := base div 10 - 1;
  writeln('base = 10 ^ ', trunc(0.434294481 * ln(base) + 0.5) : 2);
  1 : putin('give a: ', a);
  putin('give b: ', b);
  if (a[0] = 0) and (b[0] = 0) then goto 2;
  addsub(a, b, 1, c);
  putout('a + b = ', c);
  addsub(a, b, - 1, c);
  putout('a - b = ', c);
  mul(a, b, c);
  putout('a * b = ', c);
  quot(a, b, c, d);
  putout('a / b = ', c);
  putout('a % b = ', d);
  euclid(a, b, c);
  putout('gcd =   ', c);
  putin('give c: ', c);
  apowerb(a, b, c, d);
  putout('a ^ b = ', d);
  goto 1;
  2:
end {multipleprecisionpackage}.