# QMake : a generic efficient implementation of rational numbers

Initial authors : Benjamin Gregoire, Laurent Thery, INRIA, 2007

Require Import BigNumPrelude ROmega.
Require Import QArith Qcanon Qpower Qminmax.
Require Import NSig ZSig QSig.

We will build rationals out of an implementation of integers ZType for numerators and an implementation of natural numbers NType for denominators. But first we will need some glue between NType and ZType.

Module Type NType_ZType (N:NType)(Z:ZType).
Parameter Z_of_N : N.t -> Z.t.
Parameter spec_Z_of_N : forall n, Z.to_Z (Z_of_N n) = N.to_Z n.
Parameter Zabs_N : Z.t -> N.t.
Parameter spec_Zabs_N : forall z, N.to_Z (Zabs_N z) = Zabs (Z.to_Z z).
End NType_ZType.

Module Make (N:NType)(Z:ZType)(Import NZ:NType_ZType N Z) <: QType.

The notation of a rational number is either an integer x, interpreted as itself or a pair (x,y) of an integer x and a natural number y interpreted as x/y. The pairs (x,0) and (0,y) are all interpreted as 0.

Inductive t_ :=
| Qz : Z.t -> t_
| Qq : Z.t -> N.t -> t_.

Definition t := t_.

Specification with respect to QArith

Local Open Scope Q_scope.

Definition of_Z x: t := Qz (Z.of_Z x).

Definition of_Q (q:Q) : t :=
let (x,y) := q in
match y with
| 1%positive => Qz (Z.of_Z x)
| _ => Qq (Z.of_Z x) (N.of_N (Npos y))
end.

Definition to_Q (q: t) :=
match q with
| Qz x => Z.to_Z x # 1
| Qq x y => if N.eq_bool y N.zero then 0
else Z.to_Z x # Z2P (N.to_Z y)
end.

Notation "[ x ]" := (to_Q x).

Lemma N_to_Z_pos :
forall x, (N.to_Z x <> N.to_Z N.zero)%Z -> (0 < N.to_Z x)%Z.
Ltac destr_eqb :=
match goal with
| |- context [Z.eq_bool ?x ?y] =>
rewrite (Z.spec_eq_bool x y);
generalize (Zeq_bool_if (Z.to_Z x) (Z.to_Z y));
case (Zeq_bool (Z.to_Z x) (Z.to_Z y));
destr_eqb
| |- context [N.eq_bool ?x ?y] =>
rewrite (N.spec_eq_bool x y);
generalize (Zeq_bool_if (N.to_Z x) (N.to_Z y));
case (Zeq_bool (N.to_Z x) (N.to_Z y));
[ | let H:=fresh "H" in
try (intro H;generalize (N_to_Z_pos _ H); clear H)];
destr_eqb
| _ => idtac
end.

Hint Rewrite
Zplus_0_r Zplus_0_l Zmult_0_r Zmult_0_l Zmult_1_r Zmult_1_l
Z.spec_0 N.spec_0 Z.spec_1 N.spec_1 Z.spec_m1 Z.spec_opp
Z.spec_compare N.spec_compare
Z.spec_gcd N.spec_gcd Zgcd_Zabs Zgcd_1
spec_Z_of_N spec_Zabs_N
: nz.
Ltac nzsimpl := autorewrite with nz in *.

Ltac qsimpl := try red; unfold to_Q; simpl; intros;
destr_eqb; simpl; nzsimpl; intros;
rewrite ?Z2P_correct by auto;
auto.

Theorem strong_spec_of_Q: forall q: Q, [of_Q q] = q.

Theorem spec_of_Q: forall q: Q, [of_Q q] == q.

Definition eq x y := [x] == [y].

Definition zero: t := Qz Z.zero.
Definition one: t := Qz Z.one.
Definition minus_one: t := Qz Z.minus_one.

Lemma spec_0: [zero] == 0.

Lemma spec_1: [one] == 1.

Lemma spec_m1: [minus_one] == -(1).

Definition compare (x y: t) :=
match x, y with
| Qz zx, Qz zy => Z.compare zx zy
| Qz zx, Qq ny dy =>
if N.eq_bool dy N.zero then Z.compare zx Z.zero
else Z.compare (Z.mul zx (Z_of_N dy)) ny
| Qq nx dx, Qz zy =>
if N.eq_bool dx N.zero then Z.compare Z.zero zy
else Z.compare nx (Z.mul zy (Z_of_N dx))
| Qq nx dx, Qq ny dy =>
match N.eq_bool dx N.zero, N.eq_bool dy N.zero with
| true, true => Eq
| true, false => Z.compare Z.zero ny
| false, true => Z.compare nx Z.zero
| false, false => Z.compare (Z.mul nx (Z_of_N dy))
(Z.mul ny (Z_of_N dx))
end
end.

Theorem spec_compare: forall q1 q2, (compare q1 q2) = ([q1] ?= [q2]).

Definition lt n m := [n] < [m].
Definition le n m := [n] <= [m].

Definition min n m := match compare n m with Gt => m | _ => n end.
Definition max n m := match compare n m with Lt => m | _ => n end.

Lemma spec_min : forall n m, [min n m] == Qmin [n] [m].

Lemma spec_max : forall n m, [max n m] == Qmax [n] [m].

Definition eq_bool n m :=
match compare n m with Eq => true | _ => false end.

Theorem spec_eq_bool: forall x y, eq_bool x y = Qeq_bool [x] [y].

check_int : is a reduced fraction n/d in fact a integer ?

Definition check_int n d :=
match N.compare N.one d with
| Lt => Qq n d
| Eq => Qz n
| Gt => zero
end.

Theorem strong_spec_check_int : forall n d, [check_int n d] = [Qq n d].

Normalisation function

Definition norm n d : t :=
let gcd := N.gcd (Zabs_N n) d in
match N.compare N.one gcd with
| Lt => check_int (Z.div n (Z_of_N gcd)) (N.div d gcd)
| Eq => check_int n d
| Gt => zero
end.

Theorem spec_norm: forall n q, [norm n q] == [Qq n q].

Theorem strong_spec_norm : forall p q, [norm p q] = Qred [Qq p q].

Reduction function : producing irreducible fractions

Definition red (x : t) : t :=
match x with
| Qz z => x
| Qq n d => norm n d
end.

Class Reduced x := is_reduced : [red x] = [x].

Theorem spec_red : forall x, [red x] == [x].

Theorem strong_spec_red : forall x, [red x] = Qred [x].

Definition add (x y: t): t :=
match x with
| Qz zx =>
match y with
| Qz zy => Qz (Z.add zx zy)
| Qq ny dy =>
if N.eq_bool dy N.zero then x
else Qq (Z.add (Z.mul zx (Z_of_N dy)) ny) dy
end
| Qq nx dx =>
if N.eq_bool dx N.zero then y
else match y with
| Qz zy => Qq (Z.add nx (Z.mul zy (Z_of_N dx))) dx
| Qq ny dy =>
if N.eq_bool dy N.zero then x
else
let n := Z.add (Z.mul nx (Z_of_N dy)) (Z.mul ny (Z_of_N dx)) in
let d := N.mul dx dy in
Qq n d
end
end.

Theorem spec_add : forall x y, [add x y] == [x] + [y].

Definition add_norm (x y: t): t :=
match x with
| Qz zx =>
match y with
| Qz zy => Qz (Z.add zx zy)
| Qq ny dy =>
if N.eq_bool dy N.zero then x
else norm (Z.add (Z.mul zx (Z_of_N dy)) ny) dy
end
| Qq nx dx =>
if N.eq_bool dx N.zero then y
else match y with
| Qz zy => norm (Z.add nx (Z.mul zy (Z_of_N dx))) dx
| Qq ny dy =>
if N.eq_bool dy N.zero then x
else
let n := Z.add (Z.mul nx (Z_of_N dy)) (Z.mul ny (Z_of_N dx)) in
let d := N.mul dx dy in
norm n d
end
end.

Theorem spec_add_norm : forall x y, [add_norm x y] == [x] + [y].

`(Reduced x, Reduced y) : Reduced (add_norm x y).

Definition opp (x: t): t :=
match x with
| Qz zx => Qz (Z.opp zx)
| Qq nx dx => Qq (Z.opp nx) dx
end.

Theorem strong_spec_opp: forall q, [opp q] = -[q].

Theorem spec_opp : forall q, [opp q] == -[q].

Instance strong_spec_opp_norm q `(Reduced q) : Reduced (opp q).

Definition sub x y := add x (opp y).

Theorem spec_sub : forall x y, [sub x y] == [x] - [y].

Definition sub_norm x y := add_norm x (opp y).

Theorem spec_sub_norm : forall x y, [sub_norm x y] == [x] - [y].

Instance strong_spec_sub_norm x y
`(Reduced x, Reduced y) : Reduced (sub_norm x y).

Definition mul (x y: t): t :=
match x, y with
| Qz zx, Qz zy => Qz (Z.mul zx zy)
| Qz zx, Qq ny dy => Qq (Z.mul zx ny) dy
| Qq nx dx, Qz zy => Qq (Z.mul nx zy) dx
| Qq nx dx, Qq ny dy => Qq (Z.mul nx ny) (N.mul dx dy)
end.

Theorem spec_mul : forall x y, [mul x y] == [x] * [y].

Definition norm_denum n d :=
if N.eq_bool d N.one then Qz n else Qq n d.

Lemma spec_norm_denum : forall n d,
[norm_denum n d] == [Qq n d].

Definition irred n d :=
let gcd := N.gcd (Zabs_N n) d in
match N.compare gcd N.one with
| Gt => (Z.div n (Z_of_N gcd), N.div d gcd)
| _ => (n, d)
end.

Lemma spec_irred : forall n d, exists g,
let (n',d') := irred n d in
(Z.to_Z n' * g = Z.to_Z n)%Z /\ (N.to_Z d' * g = N.to_Z d)%Z.

Lemma spec_irred_zero : forall n d,
(N.to_Z d = 0)%Z <-> (N.to_Z (snd (irred n d)) = 0)%Z.

Lemma strong_spec_irred : forall n d,
(N.to_Z d <> 0%Z) ->
let (n',d') := irred n d in Zgcd (Z.to_Z n') (N.to_Z d') = 1%Z.

Definition mul_norm_Qz_Qq z n d :=
if Z.eq_bool z Z.zero then zero
else
let gcd := N.gcd (Zabs_N z) d in
match N.compare gcd N.one with
| Gt =>
let z := Z.div z (Z_of_N gcd) in
let d := N.div d gcd in
norm_denum (Z.mul z n) d
| _ => Qq (Z.mul z n) d
end.

Definition mul_norm (x y: t): t :=
match x, y with
| Qz zx, Qz zy => Qz (Z.mul zx zy)
| Qz zx, Qq ny dy => mul_norm_Qz_Qq zx ny dy
| Qq nx dx, Qz zy => mul_norm_Qz_Qq zy nx dx
| Qq nx dx, Qq ny dy =>
let (nx, dy) := irred nx dy in
let (ny, dx) := irred ny dx in
norm_denum (Z.mul ny nx) (N.mul dx dy)
end.

Lemma spec_mul_norm_Qz_Qq : forall z n d,
[mul_norm_Qz_Qq z n d] == [Qq (Z.mul z n) d].

Instance strong_spec_mul_norm_Qz_Qq z n d :
forall `(Reduced (Qq n d)), Reduced (mul_norm_Qz_Qq z n d).

Theorem spec_mul_norm : forall x y, [mul_norm x y] == [x] * [y].

Instance strong_spec_mul_norm x y :
forall `(Reduced x, Reduced y), Reduced (mul_norm x y).

Definition inv (x: t): t :=
match x with
| Qz z =>
match Z.compare Z.zero z with
| Eq => zero
| Lt => Qq Z.one (Zabs_N z)
| Gt => Qq Z.minus_one (Zabs_N z)
end
| Qq n d =>
match Z.compare Z.zero n with
| Eq => zero
| Lt => Qq (Z_of_N d) (Zabs_N n)
| Gt => Qq (Z.opp (Z_of_N d)) (Zabs_N n)
end
end.

Theorem spec_inv : forall x, [inv x] == /[x].

Definition inv_norm (x: t): t :=
match x with
| Qz z =>
match Z.compare Z.zero z with
| Eq => zero
| Lt => Qq Z.one (Zabs_N z)
| Gt => Qq Z.minus_one (Zabs_N z)
end
| Qq n d =>
if N.eq_bool d N.zero then zero else
match Z.compare Z.zero n with
| Eq => zero
| Lt =>
match Z.compare n Z.one with
| Gt => Qq (Z_of_N d) (Zabs_N n)
| _ => Qz (Z_of_N d)
end
| Gt =>
match Z.compare n Z.minus_one with
| Lt => Qq (Z.opp (Z_of_N d)) (Zabs_N n)
| _ => Qz (Z.opp (Z_of_N d))
end
end
end.

Theorem spec_inv_norm : forall x, [inv_norm x] == /[x].

Instance strong_spec_inv_norm x : Reduced x -> Reduced (inv_norm x).

Definition div x y := mul x (inv y).

Theorem spec_div x y: [div x y] == [x] / [y].

Definition div_norm x y := mul_norm x (inv_norm y).

Theorem spec_div_norm x y: [div_norm x y] == [x] / [y].

Instance strong_spec_div_norm x y
`(Reduced x, Reduced y) : Reduced (div_norm x y).

Definition square (x: t): t :=
match x with
| Qz zx => Qz (Z.square zx)
| Qq nx dx => Qq (Z.square nx) (N.square dx)
end.

Theorem spec_square : forall x, [square x] == [x] ^ 2.

Definition power_pos (x : t) p : t :=
match x with
| Qz zx => Qz (Z.power_pos zx p)
| Qq nx dx => Qq (Z.power_pos nx p) (N.power_pos dx p)
end.

Theorem spec_power_pos : forall x p, [power_pos x p] == [x] ^ Zpos p.

Instance strong_spec_power_pos x p `(Reduced x) : Reduced (power_pos x p).

Definition power (x : t) (z : Z) : t :=
match z with
| Z0 => one
| Zpos p => power_pos x p
| Zneg p => inv (power_pos x p)
end.

Theorem spec_power : forall x z, [power x z] == [x]^z.

Definition power_norm (x : t) (z : Z) : t :=
match z with
| Z0 => one
| Zpos p => power_pos x p
| Zneg p => inv_norm (power_pos x p)
end.

Theorem spec_power_norm : forall x z, [power_norm x z] == [x]^z.

Instance strong_spec_power_norm x z :
Reduced x -> Reduced (power_norm x z).

Interaction with Qcanon.Qc

Open Scope Qc_scope.

Definition of_Qc q := of_Q (this q).

Definition to_Qc q := !! [q].

Notation "[[ x ]]" := (to_Qc x).

Theorem strong_spec_of_Qc : forall q, [of_Qc q] = q.

Instance strong_spec_of_Qc_bis q : Reduced (of_Qc q).

Theorem spec_of_Qc: forall q, [[of_Qc q]] = q.

Theorem spec_oppc: forall q, [[opp q]] = -[[q]].

Theorem spec_oppc_bis : forall q : Qc, [opp (of_Qc q)] = - q.

Theorem spec_comparec: forall q1 q2,
compare q1 q2 = ([[q1]] ?= [[q2]]).

[[add x y]] = [[x]] + [[y]].

[[add_norm x y]] = [[x]] + [[y]].

Theorem spec_add_normc_bis : forall x y : Qc,
[add_norm (of_Qc x) (of_Qc y)] = x+y.

Theorem spec_subc x y: [[sub x y]] = [[x]] - [[y]].

Theorem spec_sub_normc x y:
[[sub_norm x y]] = [[x]] - [[y]].

Theorem spec_sub_normc_bis : forall x y : Qc,
[sub_norm (of_Qc x) (of_Qc y)] = x-y.

Theorem spec_mulc x y:
[[mul x y]] = [[x]] * [[y]].

Theorem spec_mul_normc x y:
[[mul_norm x y]] = [[x]] * [[y]].

Theorem spec_mul_normc_bis : forall x y : Qc,
[mul_norm (of_Qc x) (of_Qc y)] = x*y.

Theorem spec_invc x:
[[inv x]] = /[[x]].

Theorem spec_inv_normc x:
[[inv_norm x]] = /[[x]].

Theorem spec_inv_normc_bis : forall x : Qc,
[inv_norm (of_Qc x)] = /x.

Theorem spec_divc x y: [[div x y]] = [[x]] / [[y]].

Theorem spec_div_normc x y: [[div_norm x y]] = [[x]] / [[y]].

Theorem spec_div_normc_bis : forall x y : Qc,
[div_norm (of_Qc x) (of_Qc y)] = x/y.

Theorem spec_squarec x: [[square x]] = [[x]]^2.

Theorem spec_power_posc x p:
[[power_pos x p]] = [[x]] ^ nat_of_P p.

End Make.