W. F. Lunnon, Multi-dimensional map-folding, "The Computer Journal", Volume 14, Issue 1, 1971, Pages 75–80, at page 79, https://doi.org/10.1093/comjnl/14.1.75

procedure foldings (p, job); integer array p; procedure job;
  begin comment perform job (A, B) on each folding of a 
    p[1] x ... x p[d] map, where A and B are the above 
    and below vectors. p[d + 1] < 0 terminates p;

  integer d, n, j, i, m, l, g, gg, dd;
  n := 1; i := d := 0; for i := i + 1 while p[i] ≥ 0 do 
    begin d := i; n := n*p[i] end
  comment d dimensions and n leaves;

  begin integer array A, B, count, gapter [0:n], gap [0:n*n];
  comment B[m] is the leaf below leaf m in the current 
    folding, A[m] the leaf above. count[m] is the no. of 
    sections in which there is a gap for the new leaf l below 
    leaf m, gap[gapter[l - 1] + j] is the j-th (possible or 
    actual) gap for leaf l, and later gap [gapter[l]] is the 
    gap where leaf l is currently inserted;

  integer array P [0:d], C [0:d, 0:n], D [0:d, 0:n, 0:n];
  P[0] := 1; for i := 1 step 1 until d do P[i] := 
    P[i - 1]*p[i];
  for i := 1 step 1 until d do for m := 1 step 1 until n do
      C[i, m] := (m - 1)÷P[i - 1] - 
        (m - 1)÷P[i]*p[i] + 1;
  for i := 1 step 1 until d do for l := 1 step 1 until n do
    for m := 0 step 1 until l do
      D[i,l,m] := if m = 0 then 0 else
        if C[i,l] - C[i,m] = (C[i,l] - C[i,m])÷2*2 
          then

Page 80

            (if C[i,m] = 1 then m else m - P[i - 1]) 
              else
            (if C[i,m] = p[i] > m + P[i - 1] > 1 then 
              m else m + P[i - 1])
      comment P[i] = p[1] x ... x p[i], C[i, m] = i-th co-
        ordinate of leaf m, 
          D[i,l,m] = leaf connected to m in section i when 
            inserting l;

      for m := 0 step 1 until n do count[m] := 0;
      A[0] := B[0] := g := l := 0; goto entry;
      comment kick off with null folding;

down: add := 0; gg := g := gapter[l - 1];
      comment dd is the no. of sections in which l is 
                unconstrained, 
              gg the no. of possible and g the no. of 
                actual gaps for l, + gapter[l - 1];

      comment find the possible gaps for leaf l in each 
                section, then discard those not common 
                to all. All possible if dd = d;
      for i := 1 step 1 until d do if D[i,l,l] = l then 
        dd := dd + 1 else

Page 80 column 2

          for m := D[i,l,l], D[i,l,B[m]] while 
            m ≠ l do
              begin gap[gg] := m; if count[m] = 0 
                then gg := gg + 1;
              count[m] := count[m] + 1 end;

      if dd = d then for m := 0 step 1 until l - 1 do
              begin gap[gg] := m; gg := gg + 1; end;
      for j := g step 1 until gg - 1 do
              begin gap[g] := gap[j]; if 
                count[gap[j]] = d - dd then
                  g := g + 1;
              count[gap[j]] := 0 end;

      comment for each gap insert leaf l, call self recursively, 
                remove leaf l;
along: if g = gapter[l - 1] then goto up; g := g - 1;
      A[l] := gap[g]; B[l] := B[A[l]]; 
        B[A[l]] := A[B[l]] := l;
entry: gapter[l] := g; l := l + 1; if l ≤ m then goto 
        down else job(A, B);

up:   l := l - 1; B[A[l]] := B[l]; A[B[l]] := A[l];
      if l > 0 then goto along;
      end; end of foldings;
