% rational arithmetic package from Porto's programming examples

op(700, xfx, [isr, norms_to, denorms_to])!

X isr Y :- rat(Y, X), !.	% rational equivalent to "is" predicate.

rat(U + V, W) :-
	integer(U),
	integer(V), !,
	W is U + V.
rat(U + V, W) :- !,
	rat(U, X),
	rat(V, Y),
	X norms_to X1/X2,
	Y norms_to Y1/Y2,
	lcm(X2, Y2, P, X12, Y12),
	Z is X1 * X12 + Y1 * Y12,
	gcd(Z, P, _, W1, W2),
	W1/W2 denorms_to W.
rat(U * V, W) :-
	integer(U),
	integer(V), !,
	W is U * V.
rat(U * V, W) :- !,
	rat(U, X),
	rat(V, Y),
	X norms_to X1/X2,
	Y norms_to Y1/Y2,
	Z1 is X1 * Y1,
	Z2 is X2 * Y2,
	gcd(Z1, Z2, _, W1, W2),
	W1/W2 denorms_to W.
rat(-U, W) :-
	integer(W), !,
	V is - W.
rat(-U, W) :- !,
	rat(U, V),
	V norms_to U1/U2,
	V1 is -U1,
	V1/U2 norms_to W.
rat(U - V, W) :- !,
	rat(-V, V1),
	rat(U + V1, W).
rat(U/1, V) :- !,
	rat(U, V).
rat(U/V, W) :- !,
	rat(U, X),
	rat(V, Y),
	X norms_to X1/X2,
	Y norms_to Y1/Y2,
	sign(Y1, S, A),
	Z1 is X1 * Y2 * S,
	Z2 is X2 * A,
	gcd(Z1, Z2, _, W1, W2),
	W1/W2 denorms_to W.
rat(U^V, W) :-
	rat(U, X),
	X norms_to X1/X2,
	sign(Y, _, E),
	Z1 is X1 ^ E,
	Z2 is X2 ^ E,
	(Y < 0 ->
		sign(Z1, S, W1),
		W2 is S * Z2,
		W2/W1 denorms_to W
	| Z1/Z2 denorms_to W).
rat(U, U) :- integer(U).

X norms_to X/1 :-
	integer(X), !.
X/Y norms_to X/Y :-
	integer(X),
	integer(Y),
	Y > 0.

X/1 denorms_to X :- !.
X denorms_to X.

sign(X, -1, A) :-
	X < 0, !,
	A is -X.
sign(X, 1, X).

lcm(A, B, P, C, E) :-
	gcd(A, B, D, E, C),
	P is A * B / D.

gcd(A, B, C, D, E) :-
	sign(A, _, A1),
	sign(B, _, B1),
	(B1 < A1 -> gcd1(A1, B1, C) | gcd1(B1, A1, C)),
	D is A / C,
	E is B / C.

gcd1(A, 0, A) :- !.
gcd1(A, B, D) :-
	C is A mod B,
	gcd1(B, C, D).
