{************************************************************
 Parallel composition og two relations
************************************************************}

parcomp(R,S) 
  DECL XX = PROD(R*R^,S*S^);
       YY = PROD(R^*R,S^*S);
       Q1, Q2
  BEG  Q1 = p-1(XX) * R * p-1(YY)^;
       Q2 = p-2(XX) * S * p-2(YY)^
       RETURN Q1 & Q2
  END.

{************************************************************
 Solution of a chess problem: Given R as attack relation
 of a chess piece, how many possibilities exist to place k 
 such chess pieces on the chessboard such that no pair of 
 pieces attacks one another? The number k is given as the
 length of a vector with the same name
************************************************************}

card(v,k) = cardfilter(v,k+FALSE()) & -cardfilter(v,k).

indset(R,k)
  DECL H, M
  BEG  M = epsi(On1(R));
       H = -(M & R * M) & Ln1(R) * card(L1n(M)^,k)^
       RETURN -(L1n(R) * -H)^ 
  END.
       
{************************************************************
 Solution of a chess problem: Given R as attack relation
 of a chess piece, how many possibilities exist to place k 
 such chess pieces on the chessboard such that each square 
 is attacked. The number k is given as length of a vector 
 with the same name
************************************************************}

domset(R,k)
  DECL H, M
  BEG  M = epsi(On1(R));
       H = (M | R * M) & Ln1(R) * card(L1n(M)^,k)^
       RETURN -(L1n(R) * -H)^ 
  END.

{************************************************************
 Computation of the attack relations for the chess pieces
 rook, bishop, queen, king, knight on rectangular chess-
 boards. Input relations are the successor relations on 
 the sets of rows and columns
************************************************************}

up(Sx) = 
  parcomp(Sx,I(Sx)).

right(Sy) =
   parcomp(I(Sy),Sy).

pdiag(Sx,Sy) = 
  parcomp(Sx,Sy).

ndiag(Sx,Sy) = 
  parcomp(Sx^,Sy).

rook(Sx,Sy) =
  symm(trans(up(Sx))) | symm(trans(right(Sy))).

bishop(Sx,Sy) =
  symm(trans(pdiag(Sx,Sy))) | symm(trans(ndiag(Sx,Sy))).

queen(Sx,Sy) =
  rook(Sx,Sy) | bishop(Sx,Sy).

king(Sx,Sy) =
  symm(up(Sx) | right(Sy) | pdiag(Sx,Sy) | ndiag(Sx,Sy)).

knight(Sx,Sy) =
  symm(parcomp(Sx*Sx,Sy) | 
       parcomp(Sx,Sy*Sy) | 
       parcomp((Sx*Sx)^,Sy) | 
       parcomp(Sx^,Sy*Sy)).

{************************************************************
 Computation of the permanent of a relation which is viewed 
 as a quadratic 0/1-matrix. The permanent is the number of 
 all perfect matchings in the undirected and loop-free bi-
 partite graph with relation B(R,R^). If the permanent is 
 0, then the routine yields O : 1 <-> 1, otherwise it yields 
 a vector mv with the permanent as number of entries. (See
 Ph.D. thesis of B. Leoniuk, Kiel University, 2002)
************************************************************}

LineGraph(R)                       
  DECL Prod = PROD(R,R);
       VecDirR, M, H
  BEG  VecDirR = dom(p-1(Prod) * direct(R) & p-2(Prod));
       M = inj(VecDirR) * (p-1(Prod) | p-2(Prod));
       H  = M * M^
       RETURN H & -I(H)
  END.
 
Permanent(R)
  DECL LG, H, Epsi, L, leN, geN, eqN, s, res
  BEG  IF universal(dom(R)) & universal(ran(R))
         THEN LG = LineGraph(B(R,R^));
              Epsi = epsi(On1(LG));
              L = L1n(Epsi)^;
              leN = cardfilter(L, On1(R) + unit);
              geN = -cardfilter(L, On1(R));
              eqN = leN & geN;
              H = Epsi & LG * Epsi;
              s = -(L1n(LG) * H)^ & eqN;
              IF -empty(s) THEN res = s
                           ELSE res = false FI
         ELSE res = false FI
       RETURN res
  END.

{************************************************************
 Definition of some arithmetic functions on natural numbers 
 from N = [0..n], based on the  partial successor function 
 succ : N <-> N
************************************************************}

iszero(succ) = -ran(succ).

if(b,R,S) = (b * L1n(R) & R) | (-b * L1n(S) & S).

add(succ)
  DECL Pdom = PROD(succ,succ);
       pi, rho, X, Y, H
  BEG  pi  = p-1(Pdom);
       rho = p-2(Pdom);
       H   = [pi*succ^, rho*succ];
       X   = O(pi);
       Y   = if(pi*iszero(succ), rho, H*X);
       WHILE -eq(X,Y) DO
         X = Y;
         Y = if(pi*iszero(succ), rho, H*X)
         OD
       RETURN X
  END.

mult(succ)
  DECL Pdom = PROD(succ,succ);
       pi, rho, X, Y, add, H, zero
  BEG  pi   = p-1(Pdom);
       rho  = p-2(Pdom);
       add  = add(succ);
       zero = Ln1(pi)*iszero(succ)^;
       H    = [pi*succ^, rho];
       X    = O(pi);
       Y    = if(pi*iszero(succ), zero, [H*X, rho] * add); 
       WHILE -eq(X,Y) DO
         X = Y;
         Y = if(pi*iszero(succ), zero, [H*X, rho] * add)
         OD
       RETURN X
  END.

quad(succ)
  DECL Pdom = PROD(succ,succ);
       pi, rho, mult, I
  BEG  pi   = p-1(Pdom);
       rho  = p-2(Pdom);
       mult = mult(succ);
       I    = I(succ)
       RETURN [I, I] * mult
  END.

{************************************************************
 Transformation of R : X <-> Y to a vector v : X*Y <-> 1 and 
 vice versa. In VecToRel the second argument S defines the 
 direct product X*Y which is the argument domain of v. See 
 e.g., proceedings WG 86 (LNCS 246), p. 72
************************************************************}

RelToVec(R)
  DECL Prod = PROD(R*R^,R^*R);
       pi, rho, v
  BEG  pi  = p-1(Prod);
       rho = p-2(Prod);
       v   = dom(pi * R & rho)
       RETURN v
  END.

VecToRel(v,S)
  DECL Prod = PROD(S*S^,S^*S);
       pi, rho, R
  BEG  pi  = p-1(Prod);
       rho = p-2(Prod);
       R   = pi^ * (rho & v * L1n(S))
       RETURN R
  END.

{************************************************************
 Transformation of v : X <-> 1 to a point p : 2^X <-> 1 and 
 vice versa. In PoiToVec the second argument w defines the 
 set X which is the argument domain of v. See e.g., proceed-
 ings WG 86 (LNCS 246), p. 72
************************************************************}

VecToPoi(v)
  DECL Epsi, p
  BEG  Epsi = epsi(v);
       p    = syq(Epsi,v)
       RETURN p
  END.

PoiToVec(p,w)
  DECL Epsi, v
  BEG  Epsi = epsi(w);
       v    = Epsi * p
       RETURN v
  END.