/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
   Written May 26th 2006 Markus Triska triska@gmx.at
   Public domain code.
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */


seq_make(Car, Cdr, seq(Car,Cdr)).
seq_car(seq(Car,_), Car).
seq_cdr(seq(_,Cdr), Cdr).

seq_next(Seq0, Seq) :-
    seq_cdr(Seq0, Cdr0),
    call(Cdr0, Seq).


seq_filter(Pred, Seq0, Seq) :-
    seq_car(Seq0, Car),
    seq_next(Seq0, Seq1),
    (	call(Pred, Car) ->
	seq_make(Car, seq_filter(Pred,Seq1), Seq)
    ;	seq_filter(Pred, Seq1, Seq)
    ).


lp(First, Second, Seq) :-
    Diff is Second - First,
    Third is Second + Diff,
    seq_make(First, lp(Second,Third), Seq).


not_divisible_by(By, D) :- D mod By =\= 0.

primes(Ps) :-
    lp(2, 3, Seq),
    sieve(Seq, Ps).

sieve(Seq0, Seq) :-
    seq_car(Seq0, Car),
    seq_filter(not_divisible_by(Car), Seq0, Seq1),
    seq_make(Car, sieve(Seq1), Seq).


/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  We use the following encoding:
  Let d_1, d_2, ..., d_k denote the digits in the decimal expansion
  of N and let p_1, p_2, ... be the sequence of primes. The Goedel
  number of N is (p_1^d_1) * (p_2^d_2) * --- * (p_k^d_k).
  For instance, the Goedel number of 409 is 2^4*3^0*5^9 = 31250000.
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

incr(m(N0,Ns0,Ps0,Primes0,_), m(N,Ns,Ps,Primes,G)) :-
    N is N0 + 1,
    incr(Ns0, Ps0, Ns, Ps),
    (	Ps0 == Ps ->
	Primes = Primes0
    ;	seq_car(Ps, P),
	Primes = [P|Primes0]
    ),
    product_of_powers(Ns, 1, Primes, G).


product_of_powers([], P, _, P).
product_of_powers([N|Ns], Prod0, [Prime|Primes], Prod) :-
    Prod1 is Prod0*(Prime^N),
    product_of_powers(Ns, Prod1, Primes, Prod).


incr([], Ps0, [1], Ps) :- seq_next(Ps0, Ps).
incr([N0|Ns0], Ps0, [N|Ns], Ps) :-
    (	N0 < 9 ->
	N is N0 + 1,
	Ns-Ps = Ns0-Ps0
    ;	N = 0,
	incr(Ns0, Ps0, Ns, Ps)
    ).


meertens(M) :-
    primes(Ps0),
    Start = m(1,[1],Ps0,[2],2),
    meertens(Start, M).

meertens(M0, M) :- M0 = m(M, _, _, _, M).
meertens(M0, M) :-
    incr(M0, M1),
    meertens(M1, M).