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 (NN:NType)(ZZ:ZType).
Parameter Z_of_N : NN.t -> ZZ.t.
Parameter spec_Z_of_N : forall n, ZZ.to_Z (Z_of_N n) = NN.to_Z n.
Parameter Zabs_N : ZZ.t -> NN.t.
Parameter spec_Zabs_N : forall z, NN.to_Z (Zabs_N z) = Z.abs (ZZ.to_Z z).
End NType_ZType.

Module Make (NN:NType)(ZZ:ZType)(Import NZ:NType_ZType NN ZZ) <: 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 : ZZ.t -> t_
| Qq : ZZ.t -> NN.t -> t_.

Definition t := t_.

Specification with respect to QArith

Local Open Scope Q_scope.

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

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

Definition to_Q (q: t) :=
match q with
| Qz x => ZZ.to_Z x # 1
| Qq x y => if NN.eqb y NN.zero then 0
else ZZ.to_Z x # Z.to_pos (NN.to_Z y)
end.

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

Lemma N_to_Z_pos :
forall x, (NN.to_Z x <> NN.to_Z NN.zero)%Z -> (0 < NN.to_Z x)%Z.

Ltac destr_zcompare := case Z.compare_spec; intros ?H.

Ltac destr_eqb :=
match goal with
| |- context [ZZ.eqb ?x ?y] =>
rewrite (ZZ.spec_eqb x y);
case (Z.eqb_spec (ZZ.to_Z x) (ZZ.to_Z y));
destr_eqb
| |- context [NN.eqb ?x ?y] =>
rewrite (NN.spec_eqb x y);
case (Z.eqb_spec (NN.to_Z x) (NN.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
ZZ.spec_0 NN.spec_0 ZZ.spec_1 NN.spec_1 ZZ.spec_m1 ZZ.spec_opp
ZZ.spec_compare NN.spec_compare
ZZ.spec_gcd NN.spec_gcd Z.gcd_abs_l Z.gcd_1_r
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 ?Z2Pos.id 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 ZZ.zero.
Definition one: t := Qz ZZ.one.
Definition minus_one: t := Qz ZZ.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 => ZZ.compare zx zy
| Qz zx, Qq ny dy =>
if NN.eqb dy NN.zero then ZZ.compare zx ZZ.zero
else ZZ.compare (ZZ.mul zx (Z_of_N dy)) ny
| Qq nx dx, Qz zy =>
if NN.eqb dx NN.zero then ZZ.compare ZZ.zero zy
else ZZ.compare nx (ZZ.mul zy (Z_of_N dx))
| Qq nx dx, Qq ny dy =>
match NN.eqb dx NN.zero, NN.eqb dy NN.zero with
| true, true => Eq
| true, false => ZZ.compare ZZ.zero ny
| false, true => ZZ.compare nx ZZ.zero
| false, false => ZZ.compare (ZZ.mul nx (Z_of_N dy))
(ZZ.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 NN.compare NN.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 := NN.gcd (Zabs_N n) d in
match NN.compare NN.one gcd with
| Lt => check_int (ZZ.div n (Z_of_N gcd)) (NN.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 (ZZ.add zx zy)
| Qq ny dy =>
if NN.eqb dy NN.zero then x
else Qq (ZZ.add (ZZ.mul zx (Z_of_N dy)) ny) dy
end
| Qq nx dx =>
if NN.eqb dx NN.zero then y
else match y with
| Qz zy => Qq (ZZ.add nx (ZZ.mul zy (Z_of_N dx))) dx
| Qq ny dy =>
if NN.eqb dy NN.zero then x
else
let n := ZZ.add (ZZ.mul nx (Z_of_N dy)) (ZZ.mul ny (Z_of_N dx)) in
let d := NN.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 (ZZ.add zx zy)
| Qq ny dy =>
if NN.eqb dy NN.zero then x
else norm (ZZ.add (ZZ.mul zx (Z_of_N dy)) ny) dy
end
| Qq nx dx =>
if NN.eqb dx NN.zero then y
else match y with
| Qz zy => norm (ZZ.add nx (ZZ.mul zy (Z_of_N dx))) dx
| Qq ny dy =>
if NN.eqb dy NN.zero then x
else
let n := ZZ.add (ZZ.mul nx (Z_of_N dy)) (ZZ.mul ny (Z_of_N dx)) in
let d := NN.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 (ZZ.opp zx)
| Qq nx dx => Qq (ZZ.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 (ZZ.mul zx zy)
| Qz zx, Qq ny dy => Qq (ZZ.mul zx ny) dy
| Qq nx dx, Qz zy => Qq (ZZ.mul nx zy) dx
| Qq nx dx, Qq ny dy => Qq (ZZ.mul nx ny) (NN.mul dx dy)
end.

Ltac nsubst :=
match goal with E : NN.to_Z _ = _ |- _ => rewrite E in * end.

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

Definition norm_denum n d :=
if NN.eqb d NN.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 := NN.gcd (Zabs_N n) d in
match NN.compare gcd NN.one with
| Gt => (ZZ.div n (Z_of_N gcd), NN.div d gcd)
| _ => (n, d)
end.

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

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

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

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

Definition mul_norm (x y: t): t :=
match x, y with
| Qz zx, Qz zy => Qz (ZZ.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 (ZZ.mul ny nx) (NN.mul dx dy)
end.

Lemma spec_mul_norm_Qz_Qq : forall z n d,
[mul_norm_Qz_Qq z n d] == [Qq (ZZ.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 ZZ.compare ZZ.zero z with
| Eq => zero
| Lt => Qq ZZ.one (Zabs_N z)
| Gt => Qq ZZ.minus_one (Zabs_N z)
end
| Qq n d =>
match ZZ.compare ZZ.zero n with
| Eq => zero
| Lt => Qq (Z_of_N d) (Zabs_N n)
| Gt => Qq (ZZ.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 ZZ.compare ZZ.zero z with
| Eq => zero
| Lt => Qq ZZ.one (Zabs_N z)
| Gt => Qq ZZ.minus_one (Zabs_N z)
end
| Qq n d =>
if NN.eqb d NN.zero then zero else
match ZZ.compare ZZ.zero n with
| Eq => zero
| Lt =>
match ZZ.compare n ZZ.one with
| Gt => Qq (Z_of_N d) (Zabs_N n)
| _ => Qz (Z_of_N d)
end
| Gt =>
match ZZ.compare n ZZ.minus_one with
| Lt => Qq (ZZ.opp (Z_of_N d)) (Zabs_N n)
| _ => Qz (ZZ.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 (ZZ.square zx)
| Qq nx dx => Qq (ZZ.square nx) (NN.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 (ZZ.pow_pos zx p)
| Qq nx dx => Qq (ZZ.pow_pos nx p) (NN.pow_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]] ^ Pos.to_nat p.

End Make.