/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
   Written May 26th 2006 Markus Triska triska@metalevel.at
   Public domain code. Tested with Scryer Prolog.
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

:- use_module(library(clpz)).
:- use_module(library(lists)).

/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
   See prost.pl for more information about streams.
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

:- initialization(consult(prost)).

/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  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 Gödel
  number of N is (p_1^d_1) * (p_2^d_2) * --- * (p_k^d_k).
  For instance, the Gödel number of 409 is 2^4*3^0*5^9 = 31250000.
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

product_of_powers(N, Prime, Prod0, Prod) :- Prod #= Prod0*(Prime^N).


incr_([], [1], Ps0, Ps) :- seq_next(Ps0, Ps).
incr_([N0|Ns0], [N|Ns], Ps0, Ps) :-
        (   N0 #< 9 ->
            N #= N0 + 1,
            Ns-Ps = Ns0-Ps0
        ;   N = 0,
            incr_(Ns0, Ns, Ps0, 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).

%?- meertens(M).
%@    M = 81312000
%@ ;  ... .