Library CPlouffe
Require Import ssreflect String NPeano ZArith Coquelicot Reals Field Psatz Plouffe.
(******************************************************************************)
(* *)
(* COMPUTING HEXADECIMAL DIGITS OF PI WITH THE PLOUFFE FORMULA *)
(* *)
(******************************************************************************)
Ltac tlra := try lra.
Ltac tlia := try lia.
Notation "a %:R " := (INR a) (at level 10).
Open Scope nat_scope.
(******************************************************************************)
(* Some theorems for nat *)
Lemma le_minus_0 a b : a <= b -> a - b = 0.
Lemma Ndiv_minus a b c : 0 < c -> (a - b * c) / c = a / c - b.
Lemma mod_minus_mult a b c : 0 < c -> b * c <= a -> (a - b * c) mod c = a mod c.
Lemma pow_mod a b n :
(0 < n -> (a ^ b) mod n = (a mod n) ^ b mod n)%nat.
Lemma mod_between (a b c1 c2 m n : nat) :
(0 <= c1 <= c2 -> m < n -> c2 < b ^ m -> 1 < b ->
((a + c2) mod (b ^n) / (b ^ m) = a mod (b ^n) / (b ^ m)) ->
((a + c1) mod (b ^n) / (b ^ m) = a mod (b ^n) / (b ^ m))).
(******************************************************************************)
Open Scope Z_scope.
(* Some theorems from Z *)
Require Import ZArith.
Lemma pow_Zpower a b :
Z.of_nat (a ^ b) = Z.of_nat a ^ Z.of_nat b.
(******************************************************************************)
Open Scope R_scope.
(* Some theorems from R *)
Lemma Rinv_le_0_compat x : 0 < x -> 0 <= / x.
Lemma sum_f_R0_plus_r (f : nat -> R) m n :
sum_f_R0 f (S (m + n)) =
sum_f_R0 f m + sum_f_R0 (fun n => f (S m + n)%nat) n.
Lemma approx_divR a b : (0 < b)%nat -> 0 <= a%:R / b%:R - (a / b)%:R < 1.
Lemma pow_INR a n : (a ^ n)%:R = a%:R ^ n.
Lemma Int_part_IZR z : Int_part (IZR z) = z.
Lemma frac_part_IZR z : frac_part (IZR z) = 0.
Lemma le_Int_part a b : a <= b -> (Int_part a <= Int_part b)%Z.
Lemma le_Int_part_pos a : 0 <= a -> (0 <= Int_part a)%Z.
Lemma Int_part_INR b c : (0 < c)%nat ->
(Int_part (b * c%:R) / (Z.of_nat c) = Int_part b)%Z.
Definition Rdigit (b : nat) (d : nat) (r : R) :=
(Z.to_nat (Int_part ((Rabs r) * (b%:R ^ d))) mod b)%nat.
Lemma RdigitS b d r : Rdigit b d (r * b%:R) = Rdigit b (S d) r.
Lemma Rdigit_shift b d r n :
Rdigit b d (r * (b%:R ^ n)) = Rdigit b (n + d) r.
Lemma Rdigit_mod_div b r k n : (0 < b)%nat ->
Rdigit (b ^ n) 1 r =
((Z.to_nat (Int_part ((Rabs r) * b%:R ^ (k + n))))
mod (b ^ (k + n)) / (b ^ k))%nat.
Lemma Rdigit_mod_div16 r d p : (3 < p)%nat ->
Rdigit 16 d r =
((Z.to_nat (Int_part ((Rabs r) * 16 ^ d / 16 * 2 ^ p )))
mod (2 ^ p) / (2 ^ p / 16))%nat.
(******************************************************************************)
(* *)
(* COMPUTATION in NAT *)
(* *)
(******************************************************************************)
Open Scope nat_scope.
Section NComputePi.
(* Precision *)
Variable p : nat.
Let pS := 2 ^ p.
Let PpS : 0 < 2 ^ p.
(* Digit position *)
Variable d : nat.
Section power.
Variable m : nat.
(* a ^ b mod m *)
Definition NpowerMod a b m := (a ^ b) mod m.
Lemma NpowerModE a b : 0 < m ->
exists u, a ^ b = u * m + NpowerMod a b m.
End power.
Let f j i := (((16 ^ d / 16) * 2 ^ p) / (16 ^ i * (8 * i + j)%:R))%R.
Let g j i := (/ (16 ^ i * (8 * i + j)%:R))%R.
Lemma fgE j : (Series (f j) = (16 ^ d / 16) * 2 ^ p * Series (g j))%R.
Lemma PigE : ((16 ^ d / 16) * 2 ^ p * PI =
4 * (Series (f 1)) - 2 * (Series (f 4))
- (Series (f 5)) - (Series (f 6)))%R.
(* Iterative state : counter kN = k and result *)
Inductive NstateF := NStateF (k : nat) (s : nat).
(* Un pas : k = k + 1, kN = kN + 1,
s = s + (16^(d - 1 - k)) * 2^p / (8 * k + j) *)
Definition NiterF (j : nat) (st : NstateF) :=
let (k, s) := st in
let r := 8 * k + j in
let s := s + (pS * (NpowerMod 16 (d - 1 - k) r)) / r in
let s := if s <? pS then s else s - pS in
NStateF (S k) s.
Lemma NiterF_mod j k s : 0 < j ->
let (k1, s1) := NiterF j (NStateF k s) in s < pS -> s1 < pS.
Lemma base_sum_approxL j k (v := 8 * k + j) (po := NpowerMod 16 (d - 1 - k) v) :
0 < j -> k < d ->
exists u, (0 <= f j k - (pS * po / v)%:R - u%:R * 2 ^ p < 1)%R.
(* Compute \sum_{i = 0 to d - 1} (16^(d - 1 - i)) * 2^p / (8 * i + j) *)
Definition NiterL j :=
nat_iter d (NiterF j) (NStateF 0 0).
Lemma NiterL_mod j : 0 < j ->
let (k1, s1) := NiterL j in s1 < pS.
Lemma sumLE j : 0 < j -> 0 < d ->
let (k1, s1) := NiterL j in
exists u1,
(0 <= sum_f_R0 (f j) (d - 1) - s1%:R - u1%:R * 2 ^ p < d%:R)%R.
(* Iterative state : counter kN = k shift d and result *)
Inductive NstateG := NStateG (k : nat) (d : nat) (t : nat).
(* Un pas : kN = kN + 1, k = k + 1, d = d / 16,
t = t + (d/(8 * k + j)) *)
Definition NiterG (j : nat) (st : NstateG) :=
let (k, d, t) := st in
let r := 8 * k + j in
let t := t + (d / r) in
NStateG (S k) (d / 16) t.
(* Compute \sum_{i = d to infinity} (16^(d - 1 - i)) / (8 * i + j) *)
Definition NiterR j :=
nat_iter (p / 4) (NiterG j) (NStateG d (pS / 16) 0).
Lemma sumRE j : 0 < j -> 0 < p / 4 ->
let (_,_, s1) := NiterR j in
(0 <= sum_f_R0 (fun i => (f j (d + i))) (p / 4 - 1) - s1%:R < (p / 4)%:R)%R.
Lemma NiterR_mod j : 0 < j ->
let (_,_, s) := NiterR j in 15 * s < pS.
Lemma ex_series_f j : 0 < j -> @ex_series R_AbsRing R_NormedModule (f j).
Lemma series_bound j : 0 < j ->
(0 <= Series (fun k : nat => f j (d + (p / 4 + k))) < 2)%R.
Lemma bound j : 0 < j ->
exists u,
let (_, s1) := NiterL j in
let (_, _, s2) := NiterR j in
(0 <= Series (f j) - s1%:R - s2%:R - u %:R * 2 ^ p < (d + p / 4 + 2)%:R)%R.
(* Compute \sum_{i = 0 to infinity} (16^(d - 1 - i)) / (8 * i + j) *)
Definition NsumV j :=
let (_, s) := NiterL j in
let (_, _, t) := NiterR j in s + t.
Lemma NsumV_mod j : 0 < j -> NsumV j < 2 * pS.
Lemma bound_NsumV j : 0 < j ->
exists u,
(0 <= Series (f j) - (NsumV j)%:R - u %:R * 2 ^ p < (d + p / 4 + 2)%:R)%R.
Lemma main_thm (delta := d + p / 4 + 2)
(X := (4 * Series (f 1) - 2 * Series (f 4)
- Series (f 5) - Series (f 6))%R)
(Y := (4 * (NsumV 1) +
(9 * pS - (2 * NsumV 4 + NsumV 5 + NsumV 6 + 4 * delta)))%nat)
: 3 < p -> 8 * delta < 2 ^ (p - 4) ->
(Y + 8 * delta) mod 2 ^ p / 2 ^ (p - 4) = Y mod 2 ^ p / 2 ^ (p - 4)
-> Rdigit 16 d PI = Y mod 2 ^ p / 2 ^ (p - 4).
(* Extra the first digit from Plouffe formula *)
Definition NpiDigit :=
let delta := d + p / 4 + 2 in
if (3 <? p) then
if (8 * delta <? 2 ^ (p - 4)) then
let Y := (4 * (NsumV 1) +
(9 * pS - (2 * NsumV 4 + NsumV 5 + NsumV 6 + 4 * delta))) in
let v1 := (Y + 8 * delta) mod 2 ^ p / 2 ^ (p - 4) in
let v2 := Y mod 2 ^ p / 2 ^ (p - 4) in
if beq_nat v1 v2 then Some v2 else
None
else None
else None.
Lemma NpiDigit_correct k :
NpiDigit = Some k -> Rdigit 16 d PI = k.
End NComputePi.
(* Hexa conversion *)
Definition nToS n :=
let v := n in
match v with
| 0%nat => "0"%string
| 1%nat => "1"%string
| 2%nat => "2"%string
| 3%nat => "3"%string
| 4%nat => "4"%string
| 5%nat => "5"%string
| 6%nat => "6"%string
| 7%nat => "7"%string
| 8%nat => "8"%string
| 9%nat => "9"%string
| 10%nat => "A"%string
| 11%nat => "B"%string
| 12%nat => "C"%string
| 13%nat => "D"%string
| 14%nat => "E"%string
| 15%nat => "F"%string
| _ => "?"%string
end.
(* How many bits of the fixed-point like computation *)
Definition Nprecision := 14.
Definition NpiDigitF k :=
match NpiDigit Nprecision k with Some v => v | _ => 16 end.
Fixpoint rNpi k n :=
match n with
| S n1 => String.append (nToS (NpiDigitF k)) (rNpi (S k) n1)
| _ => (nToS (NpiDigitF k))
end.
(* Compute string of the hexa representation of pi *)
Definition Npi := rNpi 0.
(* Pi in hexa
3243F6A8885A308D313198A2E03707344A4093822299F31D00
82EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0A
C29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310
BA698DFB5AC2FFD72DBD01ADFB7B8E1AFED6A267E96BA7C904
5F12C7F9924A19947B3916CF70801F2E2858EFC16636920D87
1574E69A458FEA3F4933D7E0D95748F728EB658718BCD58821
54AEE7B54A41DC25A59B59C30D5392AF26013C5D1B02328608
5F0CA417918B8DB38EF8E79DCB0603A180E6C9E0E8BB01E8A3
ED71577C1BD314B2778AF2FDA55605C60E65525F3AA55AB945
748986263E8144055CA396A2AAB10B6B4CC5C341141E8CEA15
*)
(* First 5 digits of Pi *)
Definition Ndigit n := nToS (NpiDigitF n).
(* 6^th decimal of Pi *)
(******************************************************************************)
(* Turning from N to BigN *)
(******************************************************************************)
Require Import BigN.
Open Scope bigN_scope.
Notation " [[ n ]] " := (Z.to_nat [n]).
Lemma specN_add m n : [[m + n]] = ([[m]] + [[n]])%nat.
Lemma specN_sub m n : [[m - n]] = ([[m]] - [[n]])%nat.
Lemma specN_mul m n : [[m * n]] = ([[m]] * [[n]])%nat.
Lemma specN_pow m n : ([[m ^ n]] = [[m]] ^ [[n]])%nat.
Lemma specN_mod m n :
(0 < [[n]])%nat -> [[m mod n]] = ([[m]] mod [[n]])%nat.
Lemma specN_cmp m n : m <? n = ([[m]] <? [[n]])%nat.
Lemma specN_if_cmp m n p q :
[[if m <? n then p else q]] =
(if [[m]] <? [[n]] then [[p]] else [[q]])%nat.
Lemma specN_eqb m n : m =? n = (beq_nat [[m]] [[n]])%nat.
Lemma specN_if_eqp m n p q :
[[if m =? n then p else q]] =
(if beq_nat [[m]] [[n]] then [[p]] else [[q]])%nat.
Lemma specN_shiftl m n :
([[BigN.shiftl m n]] = (2 ^ [[n]]) * [[m]])%nat.
Lemma specN_shiftr m n :
[[BigN.shiftr m n]] = ([[m]] / (2 ^ [[n]]))%nat.
Lemma specN_div m n :
(0 < [[n]] -> [[m / n]] = [[m]] / [[n]])%nat.
(******************************************************************************)
Section ComputePi.
(* Precision *)
Variable p : N.
Definition pN := BigN.of_N p.
Definition pS := BigN.shiftl 1 pN.
(* Digit position *)
Variable d : N.
Definition dN := BigN.of_N d.
Section power.
Variable m : bigN.
Fixpoint powerModc a (b : positive) (c : bigN) :=
match b with
xH => (a * c) mod m
| xO b1 => powerModc ((a * a) mod m) b1 c
| xI b1 => powerModc ((a * a) mod m) b1 ((a * c) mod m)
end.
Lemma specN_powerModc a b c : (0 < [[m]])%nat ->
[[powerModc a b c]] = (([[a]] ^ (Pos.to_nat b) * [[c]]) mod [[m]])%nat.
End power.
(* a ^ b mod m *)
Definition powerMod a (b : N) m :=
match b with Npos b1 => powerModc m a b1 1 | _ =>
if m <? 2 then 0 else 1 end.
Lemma specN_powerMod a b m: (0 < [[m]])%nat ->
[[powerMod a b m]] = (([[a]] ^ (N.to_nat b)) mod [[m]])%nat.
(* Iterative state : counter kN = k and result *)
Inductive stateF := StateF (kN : bigN) (k : N) (s : bigN).
(* Un pas : k = k + 1, kN = kN + 1,
s = s + (16^(d - 1 - k)) * 2^p / (8 * k + j) *)
Definition iterF (j : bigN) (st : stateF) :=
let (kN, k, s) := st in
let r := 8 * kN + j in
let s := s + (BigN.shiftl (powerMod 16 (d - 1 - k)%N r) pN) / r in
let s := if s <? pS then s else s - pS in
StateF (kN + 1) (k + 1)%N s.
Definition eq_StateF (s1 : stateF) (s2 : NstateF) :=
let (k1,l1,n1) := s1 in
let (k2,n2) := s2 in
[[k1]] = k2 /\ N.to_nat l1 = k2 /\ [[n1]] = n2.
Lemma specN_iterF s1 s2 j : (0 < [[j]])%nat ->
eq_StateF s1 s2 ->
eq_StateF (iterF j s1) (NiterF (N.to_nat p) (N.to_nat d) [[j]] s2).
(* Compute \sum_{i = 0 to d - 1} (16^(d - 1 - i)) * 2^p / (8 * i + j) *)
Definition iterL j :=
N.iter d (iterF j) (StateF 0 0%N 0).
Lemma specN_iterL j : (0 < [[j]])%nat ->
eq_StateF (iterL j) (NiterL (N.to_nat p) (N.to_nat d) [[j]]).
(* Iterative state : counter kN = k shift d and result *)
Inductive stateG := StateG (kN : bigN) (kN : N) (d : bigN) (t : bigN).
Definition eq_StateG (s1 : stateG) (s2 : NstateG) :=
let (k1,l1,d1,n1) := s1 in
let (k2,d2,n2) := s2 in
[[k1]] = k2 /\ N.to_nat l1 = k2 /\ [[d1]] = d2 /\ [[n1]] = n2.
(* Un pas : kN = kN + 1, k = k + 1, d = d / 16,
t = t + (d/(8 * k + j)) *)
Definition iterG (j : bigN) (st : stateG) :=
let (kN, k, d, t) := st in
let r := 8 * kN + j in
let t := t + (d / r) in
StateG (kN + 1) (k + 1)%N (BigN.shiftr d 4) t.
Lemma specN_iterG s1 s2 j : (0 < [[j]])%nat ->
eq_StateG s1 s2 ->
eq_StateG (iterG j s1) (NiterG [[j]] s2).
(* Compute \sum_{i = d to infinity} (16^(d - 1 - i)) / (8 * i + j) *)
Definition iterR j :=
N.iter (p / 4) (iterG j) (StateG dN d (BigN.shiftr pS 4) 0).
Lemma specN_iterR j : (0 < [[j]])%nat ->
eq_StateG (iterR j) (NiterR (N.to_nat p) (N.to_nat d) [[j]]).
(* Compute \sum_{i = 0 to infinity} (16^(d - 1 - i)) / (8 * i + j) *)
Definition sumV j :=
let (_, _, s) := iterL j in
let (_, _, _, t) := iterR j in s + t.
Lemma specN_sumV j : (0 < [[j]])%nat ->
[[sumV j]] = NsumV (N.to_nat p) (N.to_nat d) [[j]].
(* Extra the first digit from Plouffe formula *)
Definition piDigit :=
let delta := dN + pN / 4 + 2 in
if (3 <? pN) then
if (8 * delta <? 2 ^ (pN - 4)) then
let Y := (4 * (sumV 1) +
(9 * pS - (2 * sumV 4 + sumV 5 + sumV 6 + 4 * delta))) in
let v1 := (Y + 8 * delta) mod 2 ^ pN / 2 ^ (pN - 4) in
let v2 := Y mod 2 ^ pN / 2 ^ (pN - 4) in
if v1 =? v2 then Some v2 else
None
else None
else None.
Lemma specN_piDigit :
match piDigit with
| Some x => NpiDigit (N.to_nat p) (N.to_nat d) = Some [[x]]
| None => NpiDigit (N.to_nat p) (N.to_nat d) = None
end.
Lemma piDigit_correct k :
piDigit = Some k -> Rdigit 16 (N.to_nat d) PI = [[k]].
End ComputePi.
(* How many bits of the fixed-point like computation *)
Definition precision := 40%N.
Definition piDigitF d :=
match piDigit precision d with Some k => k | _ => 16 end.
Definition NToS n := nToS (Z.to_nat (BigN.to_Z n)).
Fixpoint rpi k n :=
let v := NToS (piDigitF (N.of_nat k)) in
match n with
| S n1 => String.append v (rpi (S k) n1)
| _ => v
end.
(* Compute string of the hexa representation of pi *)
Definition pi := rpi 0.
(* Pi in hexa
3243F6A8885A308D313198A2E03707344A4093822299F31D00
82EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0A
C29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310
BA698DFB5AC2FFD72DBD01ADFB7B8E1AFED6A267E96BA7C904
5F12C7F9924A19947B3916CF70801F2E2858EFC16636920D87
1574E69A458FEA3F4933D7E0D95748F728EB658718BCD58821
54AEE7B54A41DC25A59B59C30D5392AF26013C5D1B02328608
5F0CA417918B8DB38EF8E79DCB0603A180E6C9E0E8BB01E8A3
ED71577C1BD314B2778AF2FDA55605C60E65525F3AA55AB945
748986263E8144055CA396A2AAB10B6B4CC5C341141E8CEA15
*)
(* First 500 digits of Pi by blocks of 50 *)
Definition digit n := NToS (piDigitF n).
(* 100000^th decimal of Pi *)
(******************************************************************************)
(* *)
(* COMPUTING HEXADECIMAL DIGITS OF PI WITH THE PLOUFFE FORMULA *)
(* *)
(******************************************************************************)
Ltac tlra := try lra.
Ltac tlia := try lia.
Notation "a %:R " := (INR a) (at level 10).
Open Scope nat_scope.
(******************************************************************************)
(* Some theorems for nat *)
Lemma le_minus_0 a b : a <= b -> a - b = 0.
Lemma Ndiv_minus a b c : 0 < c -> (a - b * c) / c = a / c - b.
Lemma mod_minus_mult a b c : 0 < c -> b * c <= a -> (a - b * c) mod c = a mod c.
Lemma pow_mod a b n :
(0 < n -> (a ^ b) mod n = (a mod n) ^ b mod n)%nat.
Lemma mod_between (a b c1 c2 m n : nat) :
(0 <= c1 <= c2 -> m < n -> c2 < b ^ m -> 1 < b ->
((a + c2) mod (b ^n) / (b ^ m) = a mod (b ^n) / (b ^ m)) ->
((a + c1) mod (b ^n) / (b ^ m) = a mod (b ^n) / (b ^ m))).
(******************************************************************************)
Open Scope Z_scope.
(* Some theorems from Z *)
Require Import ZArith.
Lemma pow_Zpower a b :
Z.of_nat (a ^ b) = Z.of_nat a ^ Z.of_nat b.
(******************************************************************************)
Open Scope R_scope.
(* Some theorems from R *)
Lemma Rinv_le_0_compat x : 0 < x -> 0 <= / x.
Lemma sum_f_R0_plus_r (f : nat -> R) m n :
sum_f_R0 f (S (m + n)) =
sum_f_R0 f m + sum_f_R0 (fun n => f (S m + n)%nat) n.
Lemma approx_divR a b : (0 < b)%nat -> 0 <= a%:R / b%:R - (a / b)%:R < 1.
Lemma pow_INR a n : (a ^ n)%:R = a%:R ^ n.
Lemma Int_part_IZR z : Int_part (IZR z) = z.
Lemma frac_part_IZR z : frac_part (IZR z) = 0.
Lemma le_Int_part a b : a <= b -> (Int_part a <= Int_part b)%Z.
Lemma le_Int_part_pos a : 0 <= a -> (0 <= Int_part a)%Z.
Lemma Int_part_INR b c : (0 < c)%nat ->
(Int_part (b * c%:R) / (Z.of_nat c) = Int_part b)%Z.
Definition Rdigit (b : nat) (d : nat) (r : R) :=
(Z.to_nat (Int_part ((Rabs r) * (b%:R ^ d))) mod b)%nat.
Lemma RdigitS b d r : Rdigit b d (r * b%:R) = Rdigit b (S d) r.
Lemma Rdigit_shift b d r n :
Rdigit b d (r * (b%:R ^ n)) = Rdigit b (n + d) r.
Lemma Rdigit_mod_div b r k n : (0 < b)%nat ->
Rdigit (b ^ n) 1 r =
((Z.to_nat (Int_part ((Rabs r) * b%:R ^ (k + n))))
mod (b ^ (k + n)) / (b ^ k))%nat.
Lemma Rdigit_mod_div16 r d p : (3 < p)%nat ->
Rdigit 16 d r =
((Z.to_nat (Int_part ((Rabs r) * 16 ^ d / 16 * 2 ^ p )))
mod (2 ^ p) / (2 ^ p / 16))%nat.
(******************************************************************************)
(* *)
(* COMPUTATION in NAT *)
(* *)
(******************************************************************************)
Open Scope nat_scope.
Section NComputePi.
(* Precision *)
Variable p : nat.
Let pS := 2 ^ p.
Let PpS : 0 < 2 ^ p.
(* Digit position *)
Variable d : nat.
Section power.
Variable m : nat.
(* a ^ b mod m *)
Definition NpowerMod a b m := (a ^ b) mod m.
Lemma NpowerModE a b : 0 < m ->
exists u, a ^ b = u * m + NpowerMod a b m.
End power.
Let f j i := (((16 ^ d / 16) * 2 ^ p) / (16 ^ i * (8 * i + j)%:R))%R.
Let g j i := (/ (16 ^ i * (8 * i + j)%:R))%R.
Lemma fgE j : (Series (f j) = (16 ^ d / 16) * 2 ^ p * Series (g j))%R.
Lemma PigE : ((16 ^ d / 16) * 2 ^ p * PI =
4 * (Series (f 1)) - 2 * (Series (f 4))
- (Series (f 5)) - (Series (f 6)))%R.
(* Iterative state : counter kN = k and result *)
Inductive NstateF := NStateF (k : nat) (s : nat).
(* Un pas : k = k + 1, kN = kN + 1,
s = s + (16^(d - 1 - k)) * 2^p / (8 * k + j) *)
Definition NiterF (j : nat) (st : NstateF) :=
let (k, s) := st in
let r := 8 * k + j in
let s := s + (pS * (NpowerMod 16 (d - 1 - k) r)) / r in
let s := if s <? pS then s else s - pS in
NStateF (S k) s.
Lemma NiterF_mod j k s : 0 < j ->
let (k1, s1) := NiterF j (NStateF k s) in s < pS -> s1 < pS.
Lemma base_sum_approxL j k (v := 8 * k + j) (po := NpowerMod 16 (d - 1 - k) v) :
0 < j -> k < d ->
exists u, (0 <= f j k - (pS * po / v)%:R - u%:R * 2 ^ p < 1)%R.
(* Compute \sum_{i = 0 to d - 1} (16^(d - 1 - i)) * 2^p / (8 * i + j) *)
Definition NiterL j :=
nat_iter d (NiterF j) (NStateF 0 0).
Lemma NiterL_mod j : 0 < j ->
let (k1, s1) := NiterL j in s1 < pS.
Lemma sumLE j : 0 < j -> 0 < d ->
let (k1, s1) := NiterL j in
exists u1,
(0 <= sum_f_R0 (f j) (d - 1) - s1%:R - u1%:R * 2 ^ p < d%:R)%R.
(* Iterative state : counter kN = k shift d and result *)
Inductive NstateG := NStateG (k : nat) (d : nat) (t : nat).
(* Un pas : kN = kN + 1, k = k + 1, d = d / 16,
t = t + (d/(8 * k + j)) *)
Definition NiterG (j : nat) (st : NstateG) :=
let (k, d, t) := st in
let r := 8 * k + j in
let t := t + (d / r) in
NStateG (S k) (d / 16) t.
(* Compute \sum_{i = d to infinity} (16^(d - 1 - i)) / (8 * i + j) *)
Definition NiterR j :=
nat_iter (p / 4) (NiterG j) (NStateG d (pS / 16) 0).
Lemma sumRE j : 0 < j -> 0 < p / 4 ->
let (_,_, s1) := NiterR j in
(0 <= sum_f_R0 (fun i => (f j (d + i))) (p / 4 - 1) - s1%:R < (p / 4)%:R)%R.
Lemma NiterR_mod j : 0 < j ->
let (_,_, s) := NiterR j in 15 * s < pS.
Lemma ex_series_f j : 0 < j -> @ex_series R_AbsRing R_NormedModule (f j).
Lemma series_bound j : 0 < j ->
(0 <= Series (fun k : nat => f j (d + (p / 4 + k))) < 2)%R.
Lemma bound j : 0 < j ->
exists u,
let (_, s1) := NiterL j in
let (_, _, s2) := NiterR j in
(0 <= Series (f j) - s1%:R - s2%:R - u %:R * 2 ^ p < (d + p / 4 + 2)%:R)%R.
(* Compute \sum_{i = 0 to infinity} (16^(d - 1 - i)) / (8 * i + j) *)
Definition NsumV j :=
let (_, s) := NiterL j in
let (_, _, t) := NiterR j in s + t.
Lemma NsumV_mod j : 0 < j -> NsumV j < 2 * pS.
Lemma bound_NsumV j : 0 < j ->
exists u,
(0 <= Series (f j) - (NsumV j)%:R - u %:R * 2 ^ p < (d + p / 4 + 2)%:R)%R.
Lemma main_thm (delta := d + p / 4 + 2)
(X := (4 * Series (f 1) - 2 * Series (f 4)
- Series (f 5) - Series (f 6))%R)
(Y := (4 * (NsumV 1) +
(9 * pS - (2 * NsumV 4 + NsumV 5 + NsumV 6 + 4 * delta)))%nat)
: 3 < p -> 8 * delta < 2 ^ (p - 4) ->
(Y + 8 * delta) mod 2 ^ p / 2 ^ (p - 4) = Y mod 2 ^ p / 2 ^ (p - 4)
-> Rdigit 16 d PI = Y mod 2 ^ p / 2 ^ (p - 4).
(* Extra the first digit from Plouffe formula *)
Definition NpiDigit :=
let delta := d + p / 4 + 2 in
if (3 <? p) then
if (8 * delta <? 2 ^ (p - 4)) then
let Y := (4 * (NsumV 1) +
(9 * pS - (2 * NsumV 4 + NsumV 5 + NsumV 6 + 4 * delta))) in
let v1 := (Y + 8 * delta) mod 2 ^ p / 2 ^ (p - 4) in
let v2 := Y mod 2 ^ p / 2 ^ (p - 4) in
if beq_nat v1 v2 then Some v2 else
None
else None
else None.
Lemma NpiDigit_correct k :
NpiDigit = Some k -> Rdigit 16 d PI = k.
End NComputePi.
(* Hexa conversion *)
Definition nToS n :=
let v := n in
match v with
| 0%nat => "0"%string
| 1%nat => "1"%string
| 2%nat => "2"%string
| 3%nat => "3"%string
| 4%nat => "4"%string
| 5%nat => "5"%string
| 6%nat => "6"%string
| 7%nat => "7"%string
| 8%nat => "8"%string
| 9%nat => "9"%string
| 10%nat => "A"%string
| 11%nat => "B"%string
| 12%nat => "C"%string
| 13%nat => "D"%string
| 14%nat => "E"%string
| 15%nat => "F"%string
| _ => "?"%string
end.
(* How many bits of the fixed-point like computation *)
Definition Nprecision := 14.
Definition NpiDigitF k :=
match NpiDigit Nprecision k with Some v => v | _ => 16 end.
Fixpoint rNpi k n :=
match n with
| S n1 => String.append (nToS (NpiDigitF k)) (rNpi (S k) n1)
| _ => (nToS (NpiDigitF k))
end.
(* Compute string of the hexa representation of pi *)
Definition Npi := rNpi 0.
(* Pi in hexa
3243F6A8885A308D313198A2E03707344A4093822299F31D00
82EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0A
C29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310
BA698DFB5AC2FFD72DBD01ADFB7B8E1AFED6A267E96BA7C904
5F12C7F9924A19947B3916CF70801F2E2858EFC16636920D87
1574E69A458FEA3F4933D7E0D95748F728EB658718BCD58821
54AEE7B54A41DC25A59B59C30D5392AF26013C5D1B02328608
5F0CA417918B8DB38EF8E79DCB0603A180E6C9E0E8BB01E8A3
ED71577C1BD314B2778AF2FDA55605C60E65525F3AA55AB945
748986263E8144055CA396A2AAB10B6B4CC5C341141E8CEA15
*)
(* First 5 digits of Pi *)
Definition Ndigit n := nToS (NpiDigitF n).
(* 6^th decimal of Pi *)
(******************************************************************************)
(* Turning from N to BigN *)
(******************************************************************************)
Require Import BigN.
Open Scope bigN_scope.
Notation " [[ n ]] " := (Z.to_nat [n]).
Lemma specN_add m n : [[m + n]] = ([[m]] + [[n]])%nat.
Lemma specN_sub m n : [[m - n]] = ([[m]] - [[n]])%nat.
Lemma specN_mul m n : [[m * n]] = ([[m]] * [[n]])%nat.
Lemma specN_pow m n : ([[m ^ n]] = [[m]] ^ [[n]])%nat.
Lemma specN_mod m n :
(0 < [[n]])%nat -> [[m mod n]] = ([[m]] mod [[n]])%nat.
Lemma specN_cmp m n : m <? n = ([[m]] <? [[n]])%nat.
Lemma specN_if_cmp m n p q :
[[if m <? n then p else q]] =
(if [[m]] <? [[n]] then [[p]] else [[q]])%nat.
Lemma specN_eqb m n : m =? n = (beq_nat [[m]] [[n]])%nat.
Lemma specN_if_eqp m n p q :
[[if m =? n then p else q]] =
(if beq_nat [[m]] [[n]] then [[p]] else [[q]])%nat.
Lemma specN_shiftl m n :
([[BigN.shiftl m n]] = (2 ^ [[n]]) * [[m]])%nat.
Lemma specN_shiftr m n :
[[BigN.shiftr m n]] = ([[m]] / (2 ^ [[n]]))%nat.
Lemma specN_div m n :
(0 < [[n]] -> [[m / n]] = [[m]] / [[n]])%nat.
(******************************************************************************)
Section ComputePi.
(* Precision *)
Variable p : N.
Definition pN := BigN.of_N p.
Definition pS := BigN.shiftl 1 pN.
(* Digit position *)
Variable d : N.
Definition dN := BigN.of_N d.
Section power.
Variable m : bigN.
Fixpoint powerModc a (b : positive) (c : bigN) :=
match b with
xH => (a * c) mod m
| xO b1 => powerModc ((a * a) mod m) b1 c
| xI b1 => powerModc ((a * a) mod m) b1 ((a * c) mod m)
end.
Lemma specN_powerModc a b c : (0 < [[m]])%nat ->
[[powerModc a b c]] = (([[a]] ^ (Pos.to_nat b) * [[c]]) mod [[m]])%nat.
End power.
(* a ^ b mod m *)
Definition powerMod a (b : N) m :=
match b with Npos b1 => powerModc m a b1 1 | _ =>
if m <? 2 then 0 else 1 end.
Lemma specN_powerMod a b m: (0 < [[m]])%nat ->
[[powerMod a b m]] = (([[a]] ^ (N.to_nat b)) mod [[m]])%nat.
(* Iterative state : counter kN = k and result *)
Inductive stateF := StateF (kN : bigN) (k : N) (s : bigN).
(* Un pas : k = k + 1, kN = kN + 1,
s = s + (16^(d - 1 - k)) * 2^p / (8 * k + j) *)
Definition iterF (j : bigN) (st : stateF) :=
let (kN, k, s) := st in
let r := 8 * kN + j in
let s := s + (BigN.shiftl (powerMod 16 (d - 1 - k)%N r) pN) / r in
let s := if s <? pS then s else s - pS in
StateF (kN + 1) (k + 1)%N s.
Definition eq_StateF (s1 : stateF) (s2 : NstateF) :=
let (k1,l1,n1) := s1 in
let (k2,n2) := s2 in
[[k1]] = k2 /\ N.to_nat l1 = k2 /\ [[n1]] = n2.
Lemma specN_iterF s1 s2 j : (0 < [[j]])%nat ->
eq_StateF s1 s2 ->
eq_StateF (iterF j s1) (NiterF (N.to_nat p) (N.to_nat d) [[j]] s2).
(* Compute \sum_{i = 0 to d - 1} (16^(d - 1 - i)) * 2^p / (8 * i + j) *)
Definition iterL j :=
N.iter d (iterF j) (StateF 0 0%N 0).
Lemma specN_iterL j : (0 < [[j]])%nat ->
eq_StateF (iterL j) (NiterL (N.to_nat p) (N.to_nat d) [[j]]).
(* Iterative state : counter kN = k shift d and result *)
Inductive stateG := StateG (kN : bigN) (kN : N) (d : bigN) (t : bigN).
Definition eq_StateG (s1 : stateG) (s2 : NstateG) :=
let (k1,l1,d1,n1) := s1 in
let (k2,d2,n2) := s2 in
[[k1]] = k2 /\ N.to_nat l1 = k2 /\ [[d1]] = d2 /\ [[n1]] = n2.
(* Un pas : kN = kN + 1, k = k + 1, d = d / 16,
t = t + (d/(8 * k + j)) *)
Definition iterG (j : bigN) (st : stateG) :=
let (kN, k, d, t) := st in
let r := 8 * kN + j in
let t := t + (d / r) in
StateG (kN + 1) (k + 1)%N (BigN.shiftr d 4) t.
Lemma specN_iterG s1 s2 j : (0 < [[j]])%nat ->
eq_StateG s1 s2 ->
eq_StateG (iterG j s1) (NiterG [[j]] s2).
(* Compute \sum_{i = d to infinity} (16^(d - 1 - i)) / (8 * i + j) *)
Definition iterR j :=
N.iter (p / 4) (iterG j) (StateG dN d (BigN.shiftr pS 4) 0).
Lemma specN_iterR j : (0 < [[j]])%nat ->
eq_StateG (iterR j) (NiterR (N.to_nat p) (N.to_nat d) [[j]]).
(* Compute \sum_{i = 0 to infinity} (16^(d - 1 - i)) / (8 * i + j) *)
Definition sumV j :=
let (_, _, s) := iterL j in
let (_, _, _, t) := iterR j in s + t.
Lemma specN_sumV j : (0 < [[j]])%nat ->
[[sumV j]] = NsumV (N.to_nat p) (N.to_nat d) [[j]].
(* Extra the first digit from Plouffe formula *)
Definition piDigit :=
let delta := dN + pN / 4 + 2 in
if (3 <? pN) then
if (8 * delta <? 2 ^ (pN - 4)) then
let Y := (4 * (sumV 1) +
(9 * pS - (2 * sumV 4 + sumV 5 + sumV 6 + 4 * delta))) in
let v1 := (Y + 8 * delta) mod 2 ^ pN / 2 ^ (pN - 4) in
let v2 := Y mod 2 ^ pN / 2 ^ (pN - 4) in
if v1 =? v2 then Some v2 else
None
else None
else None.
Lemma specN_piDigit :
match piDigit with
| Some x => NpiDigit (N.to_nat p) (N.to_nat d) = Some [[x]]
| None => NpiDigit (N.to_nat p) (N.to_nat d) = None
end.
Lemma piDigit_correct k :
piDigit = Some k -> Rdigit 16 (N.to_nat d) PI = [[k]].
End ComputePi.
(* How many bits of the fixed-point like computation *)
Definition precision := 40%N.
Definition piDigitF d :=
match piDigit precision d with Some k => k | _ => 16 end.
Definition NToS n := nToS (Z.to_nat (BigN.to_Z n)).
Fixpoint rpi k n :=
let v := NToS (piDigitF (N.of_nat k)) in
match n with
| S n1 => String.append v (rpi (S k) n1)
| _ => v
end.
(* Compute string of the hexa representation of pi *)
Definition pi := rpi 0.
(* Pi in hexa
3243F6A8885A308D313198A2E03707344A4093822299F31D00
82EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0A
C29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310
BA698DFB5AC2FFD72DBD01ADFB7B8E1AFED6A267E96BA7C904
5F12C7F9924A19947B3916CF70801F2E2858EFC16636920D87
1574E69A458FEA3F4933D7E0D95748F728EB658718BCD58821
54AEE7B54A41DC25A59B59C30D5392AF26013C5D1B02328608
5F0CA417918B8DB38EF8E79DCB0603A180E6C9E0E8BB01E8A3
ED71577C1BD314B2778AF2FDA55605C60E65525F3AA55AB945
748986263E8144055CA396A2AAB10B6B4CC5C341141E8CEA15
*)
(* First 500 digits of Pi by blocks of 50 *)
Definition digit n := NToS (piDigitF n).
(* 100000^th decimal of Pi *)