(*
* Exact real arithmetic (Constructive reals).
* Copyright (C) 2000 Jean-Christophe FILLIATRE
*
* This software is free software; you can redistribute it and/or
* modify it under the terms of the GNU Library General Public
*
* This software is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
*
* See the GNU Library General Public License version 2 for more details
* (enclosed in the file LGPL).
*)

(*i $Id: creal.ml,v 1.28 2002/07/23 09:13:23 filliatr Exp$ i*)

(*i*)
open Gmp
(*i*)

(*s This module implements constructive reals (exact real numbers),
following the algorithms given in Valérie Ménissier-Morain's
thesis (\small\verb!http://www-calfor.lip6.fr/~vmm/!).
In the following, pages refer to this document. *)

(*s {\bf Representation.} A constructive real is represented by an
approximation function (field [approximate]).  If $x$ is a real
number, its approximation function applied to an integer $n$ (in
type [int]) returns an integer $\ap{x}{n}$ (in type [Z.t]) such
that $|4^n\cdot x - \ap{x}{n}| < 1$.

For efficiency, we add a field [cache] to keep the best
approximation computed so far.  (Notice that it is safe to use
type [int] for the number of digits, since an integer with a
number of digits exceeding the capacity of machine integers would
exceed the memory capacity.)

The field [msd] is a cache for the most significant digit
(see section~\ref{msd} below). *)

type t = {
mutable cache : (int * Z.t) option;
mutable msd : int option;
approximate : int -> Z.t }

let create f = { cache = None; msd = None; approximate = f }

(*s Computing the approximation of [x] to precision [n] is easy:
either we have already computed a better approximation and the
result is just a shift'' of that value (Property 6 page 46), or
we compute [x.approximate n] and we cache its result before
returning it. *)

let fdiv_Bexp z n =
if n == 0 then
z
else if n > 0 then
Z.fdiv_q_2exp z (n + n)
else
Z.mul2exp z (-(n + n))

(*i
let max_prec = ref 0
let _ = at_exit (fun () -> Printf.printf "max_prec=%d\n" !max_prec)
i*)

let approx x n =
let compute () =
let z = x.approximate n in x.cache <- Some (n,z); z
in
match x.cache with
| None -> compute ()
| Some (m,a) -> if n <= m then fdiv_Bexp a (m - n) else compute ()

(*s Some useful constants in [Z.t] and [Q.t]. *)

let z_zero = Z.from_int 0
let z_one = Z.from_int 1
let z_two = Z.from_int 2
let z_three = Z.from_int 3
let z_four = Z.from_int 4

let q_half = Q.from_ints 1 2
let q_zero = Q.from_ints 0 1
let q_one = Q.from_ints 1 1
let q_two = Q.from_ints 2 1
let q_four = Q.from_ints 4 1

(*s Utility functions over [Z.t] and [Q.t]. *)

let z_gt x y = Z.cmp x y > 0
let z_le x y = Z.cmp x y <= 0

let z_between lb x up = z_le lb x && z_le x up

let z_even x = (Z.cmp (Z.cdiv_r_ui x 2) z_zero) == 0

let q_max q1 q2 = if Q.cmp q1 q2 >= 0 then q1 else q2

let q_abs q = if Q.sgn q < 0 then Q.neg q else q

(*s Roundings. Floor, ceil and Gau\ss\ rounding over [Q.t]. The Gau\ss\
rounding of $x$, written $\gauss{x}$, is the (only) integer such that
$\gauss{x} - \half \le x < \gauss{x} + \half$. *)

let q_floor q = Z.fdiv_q (Q.get_num q) (Q.get_den q)

let q_ceil q = Z.cdiv_q (Q.get_num q) (Q.get_den q)

let gauss_round q =
let q' = Q.add q q_half in
Z.fdiv_q (Q.get_num q') (Q.get_den q')

let gauss_round_z_over_4 z = Z.fdiv_q_2exp (Z.add_ui z 2) 2

(*s Addition (Algorithm 2 page 50).
We have $\ap{(x+y)}{n} = \lfloor(\ap{x}{n+1}+\ap{y}{n+1})/4\rceil$.
We do not try to cache a value for [x+y] given the cached values for
[x] and [y], if any, since it may require some computation (some shifts).
Moreover, this is exactly what will be done by the first call to
[approx] on [x+y] if the precision asked is less than $min(x,y)-2$. *)

create
(function n ->
let sn = succ n in
gauss_round_z_over_4 (Z.add (approx x sn) (approx y sn)))

(*s Negation is immediate and subtraction is the composition of
addition and negation (Algorithm 3 page 51). The cached value for
[x] is immediatly cached in [-x] (at no cost). *)

let cache_neg = function
| None -> None
| Some (n,a) -> Some (n, Z.neg a)

let neg x =
{ cache = cache_neg x.cache;
msd = x.msd;
approximate = function n -> Z.neg (approx x n) }

let sub x y = x +! (neg y)

let (-!) = sub

(*s Absolute value. *)

let abs x = create (function n -> Z.abs (approx x n))

(*s Most significant digit ([msd], Definition 9 page 47). \label{msd}
It is defined by $$\msd{x} = \min\ \{n\in Z ~|~ |x_n|>1 \}$$
Note that it does not terminate in 0. *)

let compute_msd x =
let rec look_up n =
(* $|\ap{x}{n-1}| \le 1$ *)
let xn = Z.abs (approx x n) in
if z_gt xn z_one then n else look_up (succ n)
and look_down n =
(* $|\ap{x}{n+1}| > 1$ *)
let xn = Z.abs (approx x n) in
if z_gt xn z_one then look_down (pred n) else succ n
in
let x0 = Z.abs (approx x 0) in
if z_gt x0 z_one then look_down (-1) else look_up 1

let msd x = match x.msd with
| None -> let m = compute_msd x in x.msd <- Some m; m
| Some m -> m

(*s Version of [msd] with a maximal bound on the iteration process
(used in function [mul] to avoid non-termination when
multiplicating by 0). *)

let msd_with_max m x =
let rec look_up n =
if n >= m then n else
let xn = Z.abs (approx x n) in
if z_gt xn z_one then n else look_up (succ n)
and look_down n =
let xn = Z.abs (approx x n) in
if z_gt xn z_one then look_down (pred n) else succ n
in
let x0 = Z.abs (approx x 0) in
if z_gt x0 z_one then look_down (-1) else look_up 1

(*s [mul_Bexp] and [div_Bexp] respectively multiplies and divides
an integer by $B^n$ (works whatever the sign of [n] is).
The result is a rational. *)

let mul_Bexp z n =
if n == 0 then
Q.from_z z
else if n > 0 then
Q.from_z (Z.mul2exp z (n + n))
else
Q.from_zs z (Z.mul2exp z (-(n + n)))

let bexp n = mul_Bexp z_one n

let div_Bexp z n =
if n == 0 then
Q.from_z z
else if n > 0 then
Q.from_zs z (Z.mul2exp z_one (n + n))
else
Q.from_z (Z.mul2exp z (-(n + n)))

(*s Multiplication (Algorithm 4 page 51). *)

let mul x y =
create
(function n ->
let d = (n + 2) / 2 in
let msd' = msd_with_max (n + 3 - d) in
let px = max (n - (msd' y) + 3) d
and py = max (n - (msd' x) + 3) d in
let xpx = approx x px
and ypy = approx y py in
let z = gauss_round (div_Bexp (Z.add_ui (Z.abs (Z.mul xpx ypy)) 1)
(px + py - n)) in
if Z.sgn xpx = Z.sgn ypy then z else Z.neg z)

let ( *! ) = mul

(*s Inverse (Algorithm 5 page 53) and division. *)

let inv x =
create
(function n ->
let msdx = msd x in
if n <= -msdx then
z_zero
else
let k = n + 2 * msdx + 1 in
let xk = approx x k in
let q = Q.div (bexp (k + n)) (Q.from_z xk) in
if z_gt xk z_one then q_ceil q else q_floor q)

let div x y = x *! (inv y)

let (/!) = div

(*s Square root (Algorithm 6 page 56). *)

let sqrt x =
create
(function n ->
let x2n = approx x (n + n) in
if Z.sgn x2n < 0 then invalid_arg "Creal.sqrt";
Z.sqrt x2n)

(*s Coercions from integers and rationals (Algorithm 1 page 49)
and coercion to rationals. *)

let fmul_Bexp q n =
if n == 0 then
q_floor q
else if n > 0 then
Z.fdiv_q (Z.mul2exp (Q.get_num q) (n + n)) (Q.get_den q)
else
q_floor (Q.div q (Q.from_z (Z.mul2exp z_one (-(n + n)))))

let of_z z =
{ cache = Some (0,z);
msd = None;
approximate = function n -> fmul_Bexp (Q.from_z z) n }

let of_q q = create (fmul_Bexp q)

let to_q x n =
let xn = approx x n in
Q.div (Q.from_z xn) (bexp n)

let of_int n = of_z (Z.from_int n)

let zero = of_int 0
let one = of_int 1
let two = of_int 2
let four = of_int 4

(*s Power of a real to a small integer. *)

let rec pow_int x n =
if n == 0 then
one
else if n < 0 then
inv (pow_int x (-n))
else
let y = pow_int (mul x x) (n / 2) in
if n mod 2 == 0 then y else mul y x

let rec pow_z x n =
let c = Z.cmp_si n 0 in
if c == 0 then
one
else if c < 0 then
inv (pow_z x (Z.neg n))
else
let y = pow_z (mul x x) (Z.fdiv_q_2exp n 1) in
if Z.cmp_si (Z.dmod_ui n 2) 0 == 0 then y else mul y x

(*s Alternate power series. The following function
[alternate_powerserie_] computes $B^p S$ where $S$ is the partial
sum of an alternate power serie such that the remainder is less
than $B^{-p}$, that is $S = \sum_{i=0}^{i=n}(-1)^ia_i$ with
$a_{n+1} < B^{-p}$. The alternate power serie is given by its
first term $a_0$ and a function [next] such that $a_{n+1} = \textit{next} ~ n ~ a_n$. *)

let alternate_powerserie_ a0 next p =
let eps = bexp (-p) in
let rec sum s n an =
(* [s] is already the sum up to $a_n$ *)
let asn = next n an in
if Q.cmp (q_abs asn) eps < 0 then
s
else
sum (if n mod 2 == 0 then Q.sub s asn else Q.add s asn) (n + 1) asn
in
Q.div (sum a0 0 a0) eps

(*s A specialized function to compute $atan(1/m)$ where [m] is a small
integer. *)

let arctan_reciproqual m =
let m_inverse = Q.from_ints 1 m in
let m_inverse_square = Q.mul m_inverse m_inverse in
create
(fun n ->
let eps = bexp (-n) in
let rec sum s sign k p =
(* [s] is already the sum up to $a_k$ *)
let p' = Q.mul p m_inverse_square in
let t = Q.mul p' (Q.from_ints 1 (k + 2)) in
if Q.cmp t eps < 0 then
s
else
sum (if sign then Q.add s t else Q.sub s t) (not sign) (k + 2) p'
in
fmul_Bexp (sum m_inverse false 1 m_inverse) n)

(*s $\pi$ is defined using [arctan], with the well-known formula (Algorithm
13 page 68)
$$\frac{\pi}{4} = 12 \arctan\left(\frac{1}{18}\right) + 8 \arctan\left(\frac{1}{57}\right) - 5 \arctan\left(\frac{1}{239}\right)$$ *)

let pi =
(of_int 48 *! arctan_reciproqual 18)
+! (of_int 32 *! arctan_reciproqual 57)
-! (of_int 20 *! arctan_reciproqual 239)

(*i
let pi =  (of_int 16 *! arctan_reciproqual 5)
-! (of_int 4  *! arctan_reciproqual 239)
i*)

let pi_over_2 = pi /! two

(*s Arctangent (Algorithm 12 page 64). *)

let arctan_ x =
let square_x = Q.mul x x in
let next n an =
Q.mul (Q.mul an square_x) (Q.from_ints (2 * n + 1) (2 * n + 3))
in
alternate_powerserie_ x next

let arctan_def x =
create
(function n ->
let k = max 0 (n + 1) in
let xk = approx x k in
if Z.cmp_si xk 0 == 0 then
z_zero
else
let q = Q.from_zs xk (Z.pow_ui_ui 4 k) in
(Q.div (Q.add (arctan_ q (n + 1)) q_one) q_four)
(Q.div
(bexp (n + k))
(Q.add (bexp (2 * n + 2))
(Q.from_z (Z.add (Z.mul xk xk) xk))))))

(*s The above definition of [arctan] converges very slowly when $|x|\ge 1$.
The convergence is accelerated using the following identities:
\begin{displaymath}
\begin{array}{lll}
\arctan(x)
& = -\pi/2 - \arctan(1/x) & \mbox{ when }x<-1 \\
& = -\pi/4 - \arctan((1-x^2)/(2x))/2 & \mbox{ when }x\approx-1 \\
& = +\pi/4 - \arctan((1-x^2)/(2x))/2 & \mbox{ when }x\approx1 \\
& = +\pi/2 - \arctan(1/x) & \mbox{ when }x>1
\end{array}
\end{displaymath}
We use the approximation of $x$ at order 1 to discriminate between the
cases. *)

let arctan x =
let x1 = approx x 1 in
let cmp_x1_neg_4 = Z.cmp_si x1 (-4) in
let cmp_x1_4 = Z.cmp_si x1 4 in
if cmp_x1_neg_4 < 0 then
(* $x < -1$ *)
neg (pi_over_2 +! arctan_def (inv x))
else if cmp_x1_neg_4 == 0 then
(* $x$ close to $-1$ *)
neg (pi_over_2 +! arctan_def ((one -! x *! x) /! (two *! x))) /! two
else if cmp_x1_4 == 0 then
(* $x$ close to 1 *)
(pi_over_2 -! arctan_def ((one -! x *! x) /! (two *! x))) /! two
else if cmp_x1_4 > 0 then
(* $x > 1$ *)
pi_over_2 -! arctan_def (inv x)
else
(* $x$ close to 0 *)
arctan_def x

(*s Arcsinus and arccosinus are derived from arctangent (Algorithm 14
page 69). We use $\arcsin(x)+\arccos(x)=\pi/2$ to avoid
non-termination of $\arcsin(1)$ and $\arccos(0)$. *)

let arcsin_def x = arctan (x /! (sqrt (one -! (x *! x))))

let arccos_def x = arctan ((sqrt (one -! (x *! x))) /! x)

let arcsin x =
let x1 = approx x 1 in
if z_le (Z.abs x1) z_two then arcsin_def x else pi_over_2 -! arccos_def x

let arccos x =
let x1 = approx x 1 in

if z_le (Z.abs x1) z_two then pi_over_2 -! arcsin_def x else arccos_def x

(*s Sinus (Algorithm 15 page 69). *)

let rec sin_ x p =
if Q.cmp x q_zero >= 0 then
let square_x = Q.mul x x in
let next n an =
Q.mul (Q.mul (Q.mul an square_x) (Q.from_ints 1 (2 * n + 2)))
(Q.from_ints 1 (2 * n + 3))
in
alternate_powerserie_ x next p
else
Q.neg (sin_ (Q.neg x) p)

let sin x =
let p = Z.sub_ui (approx (x /! pi) 0) 1 in
let theta = if Z.cmp_si p 0 == 0 then x else x -! ((of_z p) *! pi) in
let z = pi_over_2 in
create
(fun n ->
let k = max 2 (n + 2) in
let zk = approx z k in
let twozk = Z.mul2exp zk 1 in
let threezk = Z.mul_ui zk 3 in
let fourzk = Z.mul2exp zk 2 in
let thetak = approx theta k in
if (z_between z_zero thetak z_one) ||
(z_between (Z.sub_ui fourzk 4) thetak (Z.add_ui fourzk 4)) ||
(z_between (Z.sub_ui twozk 2) thetak (Z.add_ui twozk 2)) then
z_zero
else if z_between (Z.sub_ui zk 1) thetak (Z.add_ui zk 1) then
let bn = Z.mul2exp z_one (n + n) in
if z_even p then bn else Z.neg bn
else if z_between (Z.sub_ui threezk 3) thetak (Z.add_ui threezk 3) then
let bn = Z.mul2exp z_one (n + n) in
if z_even p then Z.neg bn else bn
else
let q = Q.from_zs thetak (Z.pow_ui_ui 4 k) in
let s = sin_ q (n + 2) in
let bw = Q.from_ints 16 1 in
let bn_k = bexp (n - k) in
let r =
if (z_between z_two thetak (Z.sub_ui zk 2)) ||
(z_between (Z.add_ui zk 2) thetak (Z.sub_ui twozk 3)) then
else
q_ceil (Q.sub (Q.div (Q.sub s q_one) bw) bn_k)
in
if z_even p then r else Z.neg r)

(*s Cosinus and tangent are derived from sinus (Algorithm 16 page 78). *)

let cos x = sin (pi_over_2 -! x)

let tan x = (sin x) /! (cos x)

(*s Euler constant [e]. *)

type sum_cache = {
mutable order : int;
mutable sum : Q.t; (* sum up to [order] *)
mutable term : Q.t; (* last term $a_{order}$ *)
mutable prec : int
}

let e =
let e_cache = { order = 1; sum = q_two; term = q_one; prec = 0 } in
create
(fun p ->
if p <= e_cache.prec then
fmul_Bexp e_cache.sum p
else
let eps = bexp (-p) in
let rec sum s n an =
let rn = Q.mul (Q.from_ints 1 n) an in
if Q.cmp rn eps <= 0 then begin
e_cache.order <- n;
e_cache.sum <- s;
e_cache.term <- an;
e_cache.prec <- p;
fmul_Bexp s p
end else
let asn = Q.mul (Q.from_ints 1 (n + 1)) an in
sum (Q.add s asn) (n + 1) asn
in
sum e_cache.sum e_cache.order e_cache.term)

(*s Natural logarithm (Algorithm 9 page 62). *)

let ln_above_1 r =
let y = Q.div (Q.sub r q_one) (Q.add r q_one) in
let y_square = Q.mul y y in
let one_minus_y_square = Q.sub q_one y_square in
fun n ->
let eps = bexp (-n) in
let rec sum s k p =
(* [s] is already the sum up to $a_k$ *)
let p' = Q.mul p y_square in
let t = Q.mul p' (Q.from_ints 1 (k + 2)) in
if Q.cmp (Q.div t one_minus_y_square) eps < 0 then
Q.mul q_two s
else
sum (Q.add s t) (k + 2) p'
in
Q.div (sum y 1 y) eps

let rec ln_ r =
if Q.cmp r q_zero <= 0 then invalid_arg "Creal.ln";
let cmp1 = Q.cmp r q_one in
if cmp1 < 0 then
(* $r < 1$ *)
let ln_inverse_r = ln_ (Q.inv r) in
(fun n -> Q.neg (ln_inverse_r n))
else if cmp1 == 0 then
(* $r = 1$ *)
(fun n -> q_zero)
else
(* $r > 1$ *)
ln_above_1 r

let ln_4 = let f = ln_above_1 q_four in create (fun n -> q_floor (f n))

let rec ln x =
let msd_x = msd x in
let k = -msd_x + 1 in
if k != 0 then
ln (x /! (of_q (bexp k))) +! (of_int k) *! ln_4
else
create
(fun n ->
let w = 2 - min 0 n in
let k = n + msd_x + w in
let xk = Q.from_z (approx x k) in
let q = Q.div xk (bexp k) in
q_floor (Q.add (Q.div (Q.add (ln_ q (n + w)) q_one) (bexp w))
(Q.div (bexp n) xk)))

let log x y = ln y /! ln x

(*s Inverses of hyperbolic functions. *)

let arcsinh x = ln (x +! sqrt (x *! x +! one))
let arccosh x = ln (x +! sqrt (x *! x -! one))
let arctanh x = ln ((one +! x) /! (one -! x)) /! two

(*s Exponential (Algorithm 7 page 57). *)

let exp_neg_ r =
(* $-1 \le r < 0$ *)
let r = q_abs r in
let next n an = Q.mul (Q.mul an r) (Q.from_ints 1 (n + 1)) in
create (fun n -> q_floor (alternate_powerserie_ q_one next n))

let exp_ r =
if Q.cmp r q_zero == 0 then
one
else
let s_floor_r = Z.add_ui (q_floor r) 1 in
mul (pow_z e s_floor_r) (exp_neg_ (Q.sub r (Q.from_z s_floor_r)))

let exp x =
create
(fun n ->
let qbn = bexp n in
let bn = of_q qbn in
let invqbn = Q.inv qbn in
let one_plus_invqbn = Q.add q_one invqbn in
let test1 () =
let lsup = log four (of_int 7 /! ln ((bn +! one) /! (bn -! one))) in
let l = Z.int_from (approx lsup 0) + 1 in
let xl = approx x l in
let log1 = q_floor (ln_ (Q.sub q_one invqbn) l) in
let log2 = q_floor (ln_ one_plus_invqbn l) in
(Z.cmp (Z.add log1 z_two) xl < 0) &&
(Z.cmp xl (Z.sub log2 z_two) < 0)
in
let test2 () =
let x0 = approx x 0 in
let m = Z.sub (q_floor (ln_ one_plus_invqbn 0)) z_two in
Z.cmp x0 m <= 0
in
if (n > 0 && test1 ()) || (n <= 0 && test2 ()) then
fmul_Bexp q_one n
else
let msd_x = msd x in
let clogBe =
if Z.cmp (approx x msd_x) z_one >= 0 then
Q.from_ints 577080 100000
else
Q.from_ints (-72134) 100000
in
let d2 = Q.div clogBe (bexp msd_x) in
let p = max 0 (n + 1) in
let d = q_max (Q.from_ints (-p) 1) d2 in
let k2 = q_ceil (Q.add d (Q.from_ints 44732 100000)) in
let k = max 1 (max msd_x (p + 1 + Z.int_from k2)) in
let bk = bexp k in
let xk = approx x k in
let xkBk = Q.div (Q.from_z xk) bk in
let exp_xkBk_p = approx (exp_ xkBk) p in
if Z.cmp exp_xkBk_p z_zero <= 0 then
z_zero
else
q_ceil (Q.mul (Q.sub (Q.div (Q.from_z exp_xkBk_p) q_four) q_one)
(Q.sub q_one (Q.inv bk))))

let pow x y = exp (y *! ln x)

(*s Hyperbolic functions. *)

let sinh x = (exp x -! exp (neg x)) /! two

let cosh x = (exp x +! exp (neg x)) /! two

let tanh x = sinh x /! cosh x

(*s Comparisons. [cmp] is absolute comparison and [rel_cmp] is comparison
up to $4^{-k}$. *)

let cmp x y =
let rec cmp_rec n =
let xn = approx x n in
let yn = approx y n in
if z_gt (Z.add_ui xn 1) (Z.sub_ui yn 1) &&
z_gt (Z.add_ui yn 1) (Z.sub_ui xn 1)
then
cmp_rec (succ n)
else
if z_le (Z.add_ui xn 1) (Z.sub_ui yn 1) then -1 else 1
in
cmp_rec 0

let rel_cmp k x y =
let rec cmp_rec n =
let xn = approx x n in
let yn = approx y n in
if z_gt (Z.add_ui xn 1) (Z.sub_ui yn 1) &&
z_gt (Z.add_ui yn 1) (Z.sub_ui xn 1) && n <= k + 2
then
cmp_rec (succ n)
else if z_le (Z.add_ui xn 1) (Z.sub_ui yn 1) then
-1
else if z_le (Z.add_ui yn 1) (Z.sub_ui xn 1) then
1
else
0
in
cmp_rec 0

(*s Coercions to and from type [float]. *)

let of_float f = of_q (Q.from_float f)

let to_float x n = Q.float_from (to_q x n)

(*s Coercion to and from type [string]. *)

let of_string s base =
try
begin
try
let n = String.length s in
let p = String.index s '.' in
let dec = n - p - 1 in
let s' = (String.sub s 0 p) ^ (String.sub s (p + 1) dec) in
of_q (Q.from_zs (Z.from_string_base base s') (Z.pow_ui_ui base dec))
with Not_found ->
of_z (Z.from_string_base base s)
end
with Invalid_argument _ -> invalid_arg "Creal.of_string"

(*s Decimal approximation of [x] at order [p]. We look for an integer [n]
such that $|10^px - n| < 1/2$ i.e. the integer closest to
$10^px$. There is sometimes no such integer but then we can find a
decimal approximation at order [p+1].

We first compute $y = 10^px$ and [approx y 3] i.e. an approximation
$y_3/64$ of $y$. Let $q$ and $r$ be the quotient and remainder of the
division $y_3/64$ such that $y_3 = 64q+r$ and $0\le r<63$.
If $r\le 31$ then $n$ is $q$; If $r\ge 33$ then $n$ is $q+1$;
Otherwise $10q+5$ is a decimal approximation of $x$ at order $p+1$. *)

let to_string_aux x p =
if p < 0 then invalid_arg "Creal.to_string";
let tenp = Z.pow_ui_ui 10 p in
let y = mul (of_z tenp) x in
let y3 = approx y 3 in
let q,r = Z.fdiv_qr_ui y3 64 in
let r = Z.int_from r in
let n,p =
if r <= 31 then q, p
else if r >= 33 then Z.add_ui q 1, p
else Z.add_ui (Z.mul_ui q 10) 5, succ p
in
let ns = Z.string_from (Z.abs n) in
let lns = String.length ns in
let ins,dns =
if lns >= p+1 then
String.sub ns 0 (lns - p), String.sub ns (lns - p) p
else
"0", String.make (p - lns) '0' ^ ns
in
Z.sgn n, ins, dns

let to_string x p =
let sgn,i,f = to_string_aux x p in
(if sgn < 0 then "-" else "") ^ i ^ "." ^ f

(*s Coercion to type [string] with digits packed 5 by 5. *)

let string_concat = String.concat ""

let beautiful s =
let n = String.length s in
let eol i = if (i + 5) mod 65 == 0 then "\n" else " " in
let rec cut i =
String.sub s i (min 5 (n - i)) ::
if i < n - 5 then eol i :: cut (i + 5) else []
in
string_concat (cut 0)

let to_beautiful_string x p =
let sgn,i,f = to_string_aux x p in
let nl = if String.length i + String.length f > 75 then "\n" else "" in
(if sgn < 0 then "-" else "") ^ i ^ "." ^ nl ^ beautiful f

module Infixes = struct