/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Written May 26th 2006 Markus Triska triska@metalevel.at Public domain code. - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ :- 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 .