goodstein(S, B, N, [S|Ls]) :-
	goodstein_(S, B, N, Ls).
goodstein_(_S, _B, 1, []).
goodstein_(S0, B, N, [S|Ss]) :-
	N > 1,
	goodstep(S0, B, B1, S),
	N1 is N - 1,
	goodstein_(S, B1, N1, Ss).

good(S, B, N, L) :-
	goodstein(S, B, N, L).

goodstep(N0, B, B1, N) :-
	hbn(N0, B, L),
	B1 is B + 1,
	substitute(B1, L, L0),
	eval(L0, N1),
	N is N1 - 1.

eval([], 0).
eval(N, N) :-
	number(N).
eval([C*B^Exp|Rest], X) :-
	eval(Exp, ExpE),
	eval(Rest, RestE),
	X is C * B ^ ExpE + RestE.

substitute(_B, [], []).
substitute(_B, C, C) :-
	number(C).
substitute(B, [X|Xs], [Y|Ys]) :-
	X = C * _ ^ Exp,
	substitute(B, Exp, ExpS),
	Y = C * B ^ ExpS,
	substitute(B, Xs, Ys).

%%%

hbn(N, B, L) :-
	coefs(N, B, Cs),
	length(Cs, Exp0),
	Exp is Exp0 - 1,
	notation(B, Exp, Cs, L).

notation(_B, _Exp, [], []).
notation(B, Exp, [0|Cs], Ls) :-
	Exp1 is Exp - 1,
	notation(B, Exp1, Cs, Ls).
notation(B, Exp, [C|Cs], [L|Ls]) :-
	Exp1 is Exp - 1,
	C > 0,
	notation(B, Exp1, Cs, Ls),
	notation_element(C, B, Exp, L).

notation_element(C, B, Exp, X) :-
	Exp < B,
	X = C * B ^ Exp.
notation_element(C, B, Exp, X) :-
	Exp >= B,
	hbn(Exp, B, Exp0),
	X = C * B ^ Exp0.

coefs(N, B, Cs) :-
	max_exponent(N, B, Exp),
	base_exp_num_coefs(B, Exp, N, Cs).

max_exponent(N, B, Max) :-
	max_exponent(N, B, 0, Max).
max_exponent(N, B, E0, E) :-
	succ(E0, E1),
	(  B ^ E1 > N
	-> E = E1
	;  max_exponent(N, B, E1, E)
	).

base_exp_num_coefs(_Base, N, _I, []) :-
	N < 0.
base_exp_num_coefs(Base, N, I, [C|Coefs]) :-
	N >= 0,
	NN is Base ^ N,
	C is I // NN,
	N1 is N - 1,
	I1 is I - (C*NN),
	base_exp_num_coefs(Base, N1, I1, Coefs).