Evaluating Real numbers to decimal representation

Using Coq
1 
Require Import Reals.
(* Euler's formula in Horner's form *)
Definition pi := (nat_rect _ 0 (fun n p => 2 + p * INR n / (2 * INR n + 1)))%R.

What is the easiest way to evaluate pi 100 and print the first 15 decimals?

2 
From Coq Require Import Reals.
From Interval Require Import Tactic.
Definition pi := (nat_rect _ 0 (fun n p => 2 + p * IZR (Z_of_nat n) / (2 * IZR (Z_of_nat n) + 1)))%R.
Goal True.
let v := eval simpl in (pi 100) in
interval_intro v.

gives

H : (8961703559183965 / 2251799813685248 <=
     2 + (2 + (2 + (2 + (2 + (2 + (2 + (2 +  (2 +
             (2 + (2 + (...) * 89 / (... + 1)) * 90 / (2 * 90 + 1)) * 91 / (2 * 91 + 1)) *
            92 / (2 * 92 + 1)) * 93 / (2 * 93 + 1)) * 94 / (2 * 94 + 1)) * 95 /
         (2 * 95 + 1)) * 96 / (2 * 96 + 1)) * 97 / (2 * 97 + 1)) * 98 / 
      (2 * 98 + 1)) * 99 / (2 * 99 + 1) <= 8961703559183979 / 2251799813685248)%R

Now, you need to use some other tool to convert the bounds to decimal, which gives you 3.979795852508501 <= ... <= 3.979795852508507, hence your fifteen digits.

I had to change INR to IZR in your definition of pi, otherwise Coq kind of dies while displaying H, but in theory it should also work with INR for smaller values of n.

If you do not want to use an external tool to get the decimal digits, you can use a slightly more complicated formula:

let v := eval simpl in (IZR (Raux.Zfloor (pi 100 * powerRZ 10 14))) in
interval_intro v with (i_prec 90).
3

Thanks!

How utterly embarrassing to compute pi to 3.97…

Require Import Reals.
Open Scope R.
From Interval Require Import Tactic. (* 'interval' tactic *)

Require Import List. Import ListNotations.
Definition Dec := fold_right (fun a b => a + b/10) 0.

Definition pi m := (nat_rect _ 0 (fun n p => 2 + p * IZR (Z_of_nat (m-n)) / (2 * IZR (Z_of_nat (m-n)) + 1)) m)%R.

Definition x := Dec [3;1;4;1;5;9;2;6;5;3;5;8;9;7;9].
Definition eps := 1/100000000000000.

Goal x < pi 100 < x + eps.
  split; unfold x, eps, pi; simpl; interval.
Qed.

Goal True.

  let X := (eval simpl in (pi 20)) in
  let x := (eval unfold pi in X) in
  let y := (eval simpl     in x) in
  interval_intro y with (i_prec 100) as EQ;
    replace y with X in * by reflexivity.

(*
  EQ : 995610223197390062935504709050 / 316912650057057350374175801344 <= 
       pi 20 <= 995610223197390062935504709052 / 316912650057057350374175801344
  ============================
  True
 *)
1 Like
4

One can use the exact real numbers from math-classes/corn. I can try to give more instructions later when I find the time.
https://math-classes.github.io/

@VincentSe may also be interested

1 Like
5

@larsr Using the CoRN library you can do this to get the decimal expansion of π to 50 digits:

From CoRN Require Import CRtrans.

Local Open Scope Q_scope.

(* This file illustrates how to use the computational reals CR *)

Definition answer (n:positive) (r:CR) : Z :=
 let m := (iter_pos _ (Pmult 10) 1%positive n) in
 let (a,b) := (approximate r (Qpos2QposInf (1#m)))*(Zpos m#1) in
 Z.div a (Zpos b).

(* approximate pi to 50 digits *)
Time Eval vm_compute in answer 50 (CRpi)%CR.
(*      = 314159265358979323846264338327950288419716939937510%Z *)
(*      : Z                                                     *)
(* Finished transaction in 0.105 secs (0.104u,0.s) (successful) *)