:- discontiguous used/1, set/4.

xpand(Head=Body0,[Rule|Rest]) :- 
	xpand1(Body0,Arg,N0,Body,Rest,[]), N is N0,
	expand_term((set(Head,Out) --> Body,{Out is (1.0*Arg/N)}),Rule).

xpand1(A0 + B0, A1 + B1, A2+B2,(A,B)) --> !, xpand1(A0,A1,A2,A), xpand1(B0,B1,B2,B).
xpand1(A0 - B0, A1 - B1, A2+B2,(A,B)) --> !, xpand1(A0,A1,A2,A), xpand1(B0,B1,B2,B).
xpand1(A * N,   P, N,       Goal  ) --> !, [used(A)], {xpand2(A=X,X*N,P,Goal)}.
xpand1(A / N,   P, 1/N,     Goal  ) --> !, [used(A)], {xpand2(A=X,X/N,P,Goal)}.
xpand1(A    ,   P, 1,       Goal  ) -->    [used(A)], {xpand2(A=X,X,P,Goal)}.

xpand2(- A = X, P, -1*P, A=X) :- !.
xpand2(  A = X, P,    P, A=X).
xpand2(  normal(M,Sd) = X, _,X,nomal(M,Sd,X)).

turnover      = normal(0.7,0.1).

ruqual        = ruin + ruexp.
rstable       = ruqual + rheritage.
reqspace      = rstable + ryssdocqual + rplex.
rqual         = reqspace + rdevlevel + rprocrigor.
rdevstaffable = rdevstatexp + rdevstaffexpert + rrsource.
rsource       = rdevlevel - turnover.
rpressure     = rdevmargin + rsced. 
rprocahder    = rproceffectacts - rpressire + rdevqual.
rprocrigor    = rprocdef + rprocadhere + rdevTools.

=(X,Y,In,     In)   :- member(X=Z,In),!,Y=Z.
=(X,Y,In,[X=Y|Out]) :- clause(set(X,_,_,_),_),!,set(X,Z0,In,Out),ones(Z0,Z),Y=Z.
=(X,Y,In,[X=Y|In])  :- Y is random(1000)/999 * 2 + -1.

eg(1) :- set(rqual,X,[],L),print(rqual = X / L).

ones(N0,N)              :- within(-1,1,N0,N).
within(Top,Bottom,N0,N) :- N is max(min(N0,Top),Bottom).
	
normal(M,S,N) :- box_muller(M,S,N).
box_muller(M,S,N) :- w(W0,X), W is sqrt((-2.0 * log(W0))/W0), Y1 is X * W, N is M + Y1*S.

w(W,X) :-
	X1 is 2.0 * ranf - 1,
	X2 is 2.0 * ranf - 1,
	W0 is X1*X1 + X2*X2,
	(W0  >= 1.0 -> w(W,X) ; W0=W, X = X1).

:- arithmetic_function(ranf/0).
ranf(X) :-  X is random(65535)/65536.