/* Let us count the ways [08aug05] */

/*
		<Example>

		$rule @p(_,_) $if @q(_,_), @r(_,_) .
Slots:		         1 2         3 4      5 6

Mapping:	1,3 -> X ; 4,5 -> T ; 2,6 -> Y

Result:		$rule @p(X,Y) $if @q(X,T), @r(T,Y) .

		How many logically different ways 
		can the six slots in this rule 
		be mapped to variables?

$map !slot -> {$number}

$map !set -> *!slot

$topic @partitions(!set,*!set)
$topic @map_slot(!slot,*!set,*!set)

$para @partitions
  $rule @partitions([], []).
  $rule @partitions([Slot|Slots], NewPartition) $if
     @partitions(Slots, Partition),
     @map_slot(Slot, Partition, NewPartition).
$end_para

$para @map_slot
  $rule @map_slot(Slot, [], [[Slot]]).
  $rule @map_slot(Slot, [Set|RemPartition], [[Slot|Set]|RemPartition]).
  $rule @map_slot(Slot, [Set|RemPartition], [Set|NewRemPartition]) $if
     @map_slot(Slot, RemPartition, NewRemPartition).
$end_para


  See also

  Newsgroup: comp.lang.prolog
  Subject: Re: occupancy (?) problem
  Date: 1996/07/29
  Message-ID: <4tj700$36q@arcadia.informatik.uni-muenchen.de>

*/


partitions([], []).
partitions([Slot|Slots], NewPartition) :-
     partitions(Slots, Partition),
     map_slot(Slot, Partition, NewPartition).

map_slot(Slot, [], [[Slot]]).
map_slot(Slot, [Set|RemPartition], [[Slot|Set]|RemPartition]).
map_slot(Slot, [Set|RemPartition], [Set|NewRemPartition]) :-
     map_slot(Slot, RemPartition, NewRemPartition).


/*

  Definition: Bell Numbers

  [see http://www.mathematik.uni-bielefeld.de/~sillke/SEQUENCES/series005]

  b(0) = 1

  b(n) = sum{k=0..n-1} b(k)*binomial(n-1,k), n > 0

*/

b(0,1).
b(N,BN) :-
  N > 0,
  NM1 is N-1,
  sum_k(0,NM1,0,BN).

sum_k(K,NM1,In,In) :-
  K > NM1.
sum_k(K,NM1,In,Out) :-
  \+ K > NM1,
  b(K,BK),
  binomial(NM1,K,BNM1K),
  In1 is In + BK*BNM1K,
  K1 is K+1,
  sum_k(K1,NM1,In1,Out).

binomial(N,K,BNK) :-
  NMKP1 is N-K+1,
  product(NMKP1,N,1,Numerator),
  product(1,K,1,Denominator),
  BNK is Numerator/Denominator.

product(M,N,In,In) :-
  M > N.
product(M,N,In,Out) :-
  \+ M > N,
  In1 is M*In,
  M1 is M+1,
  product(M1,N,In1,Out).

/*

 Bell's Triangle

      1
      1    2
      2    3    5
      5    7   10   15
     15   20   27   37   52
     52   67   87  114  151  203
    203  ...

 [http://mathworld.wolfram.com/BellTriangle.html]


% ---- Definition: 

  Let

    bell_row(1) = <1>

    bell_row(n) = <b(1),...,b(n-1), b(n)>, n > 1.

  Then

    bell_row(n+1) = <c(1), ..., c(n+1)>,

  where

   c(1) = b(n)

   c(k) = b(k) + c(k-1) for k=2, ..., n 

  and

   c(n+1) = b(n) + c(n).

% ---- Implementation:

/*
  bell_row(*$integer, *$integer) 
  ------------------------------
  input:  <b(1),...,b(n-1),b(n)> 
  ouput: <<c(1),...,c(n-1),c(n)>[c(n+1)>
*/
bell_row(<Bs_1_to_n_minus_1[B_n>, <Cs_1_to_n[C_n_plus_1>) :-
  C_1 = B_n,
  combine(<Bs_1_to_n_n_minus_1[B_n>, C_1, Cs_1_to_n),
  <_[C_n> = Cs_1_to_n,
  C_n_plus_1 is B_n + C_n.

/*
  combine(*$integer, $integer, $integer)
  --------------------------------------
  input:  <b(1),...,b(n-1),b(n)>, c(1) 
  output: <c(1),...,c(n-1),c(n)>
*/
combine(<>, _, <>).
combine(<B_1]Bs_2_to_n>, C_1, <C_1]Cs_2_to_n>) :-
  C_2 is B_1 + C_1,
  combine(Bs_2_to_n, C_2, Cs_2_to_n).
*/

bell_row([B1|Bs_2_to_n], Cs_1_to_n_plus_1) :-
  get_last([B1|Bs_2_to_n], B_n),
  C_1 = B_n,
  combine([B1|Bs_2_to_n], C_1, Cs_1_to_n),
  get_last(Cs_1_to_n, C_n),
  C_n_plus_1 is B_n + C_n,
  put_last(C_n_plus_1, Cs_1_to_n, Cs_1_to_n_plus_1).

combine([], _, []).
combine([B_1|Bs_2_to_n], C_1, [C_1|Cs_2_to_n]) :-
  C_2 is B_1 + C_1,
  combine(Bs_2_to_n, C_2, Cs_2_to_n).

get_last([X], X).
get_last([_,X|L], Y) :-
   get_last([X|L], Y).

put_last(X, [], [X]).
put_last(X, [Y|L1], [Y|L2]) :-
  put_last(X, L1, L2).

bell_triangle(Bs) :-
  write(Bs),
  bell_row(Bs, NewBs),
  ( more -> nl, !, bell_triangle(NewBs) ; nl, !, write('done'), nl ).
 
/* 
  Note: 

   <> = []
   <X]L> = [X]

*/

more :-
  get_single_char(13).

run :-
  bell_triangle([1]).


/*
  Many thanks to Markus Triska 
  [Message ID <87ps5sdmbv.fsf@logic.at>] 
  for pointing the following out to me:

 By introducing an additional argument that is used to report solutions
incrementally, this kind of querying can be delegated to the toplevel
{by utilizing the built-in SWI-Prolog predicate 'last/2'
instead of using the 'get_last/2' and 'put_last/3' defined above}:
*/

   aitken(R, R).
   aitken(R0, R) :-
       last(R0, L),
       row(R0, L, R1),
       aitken([L|R1], R).

   row([], _, []).
   row([A|As], B, [C|Cs]) :-
       C is A + B,
       row(As, C, Cs).
/*
 Example query:

   %?- aitken([1], R).
   %@ R = [1];
   %@ R = [1, 2] ;
   %@ R = [2, 3, 5] ;
   %@ R = [5, 7, 10, 15] ;
   %@ R = [15, 20, 27, 37, 52] a
   %@ Yes
*/

new_bell_row(R, R).
new_bell_row(R0, R) :-
  last(R0, L),	/* SWI-Prolog built-in [see ':- help(last).'] */
  new_combine(R0, L, R1),
  new_bell_row([L|R1], R).

new_combine([], _, []).
new_combine([B_1|Bs_2_to_n], C_0, [C_1|Cs_2_to_n]) :-
  C_1 is B_1 + C_0,
  new_combine(Bs_2_to_n, C_1, Cs_2_to_n).

new_bell_triangle(R) :-
  new_bell_row([1], R).

new_run :-
  new_bell_triangle(R),
  write(R),
  ( more -> nl, fail ; !, nl, write('done'), nl ).

/* (c) Bill Hogan - April 29, 2007 */