Why3 Standard Library index

Unbounded floating-point numbers

See also this report

module UFloat
  use real.RealInfix
  use real.FromInt
  use real.Abs
  use ieee_float.RoundingMode

  type t

  function uround mode real : t
  function to_real t : real
  function of_int int : t
  axiom to_real_of_int : forall x [of_int x]. to_real (of_int x) = from_int x

  constant eps:real
  constant eta:real
  axiom eps_bounds : 0. <. eps <. 1.
  axiom eta_bounds : 0. <. eta <. 1.

  (* To avoid "inline_trivial" to break the forward_propagation strategy *)
  meta "inline:no" function eps
  meta "inline:no" function eta

  let ghost function uadd (x y:t) : t
    (* TODO: Do we want the two first assertions in our context ?
      We only use them to prove the addition lemma *)
    ensures { abs (to_real result -. (to_real x +. to_real y)) <=. abs (to_real x) }
    ensures { abs (to_real result -. (to_real x +. to_real y)) <=. abs (to_real y) }
    ensures {
      abs (to_real result -. (to_real x +. to_real y))
      <=. abs (to_real x +. to_real y) *. eps
    }
  = uround RNE (to_real x +. to_real y)

  let ghost function usub (x y:t) : t
    (* TODO: Do we want the two first assertions in our context ?
      We only use them to prove the addition lemma *)
    ensures { abs (to_real result -. (to_real x -. to_real y)) <=. abs (to_real x) }
    ensures { abs (to_real result -. (to_real x -. to_real y)) <=. abs (to_real y) }
    ensures {
      abs (to_real result -. (to_real x -. to_real y))
      <=. abs (to_real x -. to_real y) *. eps
    }
  = uround RNE (to_real x -. to_real y)

  let ghost function umul (x y:t) : t
    ensures {
      abs (to_real result -. (to_real x *. to_real y))
        <=. abs (to_real x *. to_real y) *. eps +. eta
    }
  = uround RNE (to_real x *. to_real y)

  let ghost function udiv (x y:t) : t
    requires { to_real y <> 0. }
    ensures {
      abs (to_real result -. (to_real x /. to_real y))
        <=. abs (to_real x /. to_real y) *. eps +. eta
    }
  = uround RNE (to_real x /. to_real y)

  let ghost function uminus (x:t) : t
    ensures { to_real result = -. (to_real x) }
  = uround RNE (-. (to_real x))

  let ghost function ( ++. ) (x:t) (y:t) : t = uadd x y
  let ghost function ( --. ) (x:t) (y:t) : t = usub x y
  let ghost function ( **. ) (x:t) (y:t) : t = umul x y
  (* Why3 doesn't support abbreviations so we need to add the requires *)
  let ghost function ( //. ) (x:t) (y:t) : t
    requires { to_real y <> 0. }
  = udiv x y
  let ghost function ( --._ ) (x:t) : t = uminus x

Infix operators

  (* Some constants *)
  constant u0:t
  axiom to_real_u0 : to_real u0 = 0.0
  constant u1:t
  axiom to_real_u1 : to_real u1 = 1.0
  constant u2:t
  axiom to_real_u2 : to_real u2 = 2.0
  constant u4:t
  axiom to_real_u4 : to_real u4 = 4.0
  constant u8:t
  axiom to_real_u8 : to_real u8 = 8.0
  constant u16:t
  axiom to_real_u16 : to_real u16 = 16.0
  constant u32:t
  axiom to_real_u32 : to_real u32 = 32.0
  constant u64:t
  axiom to_real_u64 : to_real u64 = 64.0
  constant u128:t
  axiom to_real_u128 : to_real u128 = 128.0
  constant u256:t
  axiom to_real_u256 : to_real u256 = 256.0
  constant u512:t
  axiom to_real_u512 : to_real u512 = 512.0
  constant u1024:t
  axiom to_real_u1024 : to_real u1024 = 1024.0
  constant u2048:t
  axiom to_real_u2048 : to_real u2048 = 2048.0
  constant u4096:t
  axiom to_real_u4096 : to_real u4096 = 4096.0
  constant u8192:t
  axiom to_real_u8192 : to_real u8192 = 8192.0
  constant u16384:t
  axiom to_real_u16384 : to_real u16384 = 16384.0
  constant u32768:t
  axiom to_real_u32768 : to_real u32768 = 32768.0
  constant u65536:t
  axiom to_real_u65536 : to_real u65536 = 65536.0
end

module USingle
  use real.RealInfix

Single-precision unbounded floats

  type usingle

  constant eps:real = 0x1p-24 /. (1. +. 0x1p-24)
  constant eta:real = 0x1p-150

  clone export UFloat with
    type t = usingle,
    constant eps = eps,
    constant eta = eta,
    axiom.
end

module UDouble
  use real.RealInfix
  type udouble

Double-precision unbounded floats

  constant eps:real = 0x1p-53 /. (1. +. 0x1p-53)
  constant eta:real = 0x1p-1075

  clone export UFloat with
    type t = udouble,
    constant eps = eps,
    constant eta = eta,
    axiom.
end

(* Helper lemmas to help the proof of propagation lemmas *)
module HelperLemmas
  use real.RealInfix
  use real.Abs

  let ghost div_order_compat (x y z:real)
    requires { x <=. y }
    requires { 0. <. z }
    ensures { x /. z <=. y /. z }
    = ()

  let ghost mult_err (x exact_x x' x_rel_err x_cst_err y:real)
    requires { 0. <=. x_rel_err }
    requires { 0. <=. x_cst_err }
    requires { abs exact_x <=. x' }
    requires { abs (x -. exact_x) <=. x_rel_err *. x' +. x_cst_err }
    ensures { abs (x *. y -. exact_x *. y) <=. x_rel_err *. abs (x' *. y) +. x_cst_err *. abs y }
  =
  assert {
    y >=. 0. ->
    abs (x *. y -. exact_x *. y) <=. abs (x_rel_err *. x' *. y) +. x_cst_err *. abs y
    by
      (exact_x -. x_rel_err *. x' -. x_cst_err) *. y <=. x *. y <=. (exact_x +. x_rel_err *. x' +. x_cst_err) *. y
  };
  assert {
    y <. 0. ->
    abs (x *. y -. exact_x *. y) <=. abs (x_rel_err *. x' *. y) +. x_cst_err *. abs y
    by
      (exact_x +. x_rel_err *. x' +. x_cst_err) *. y <=. x *. y <=. (exact_x -. x_rel_err *. x' -. x_cst_err) *. y
  }

  let ghost mult_err_combine (x exact_x x' x_rel_err x_cst_err y exact_y y' y_rel_err y_cst_err:real)
    requires { 0. <=. x_rel_err }
    requires { 0. <=. y_rel_err }
    requires { 0. <=. x_cst_err }
    requires { 0. <=. y_cst_err }
    requires { abs exact_x <=. x' }
    requires { abs exact_y <=. y' }
    requires { abs (x -. exact_x) <=. x_rel_err *. x' +. x_cst_err }
    requires { abs (y -. exact_y) <=. y_rel_err *. y' +. y_cst_err }
    ensures {
      abs (x *. y -. exact_x *. exact_y)
      <=. (x_rel_err +. y_rel_err +. x_rel_err *. y_rel_err) *. (x' *. y')
        +. (y_cst_err +. y_cst_err *. x_rel_err) *. x'
        +. (x_cst_err +. x_cst_err *. y_rel_err) *. y'
        +. x_cst_err *. y_cst_err
    }
  =
  mult_err x exact_x x' x_rel_err x_cst_err y;
  mult_err y exact_y y' y_rel_err y_cst_err exact_x;
  mult_err y exact_y y' y_rel_err y_cst_err x';
  assert {
    abs (x *. y -. exact_x *. exact_y) <=. (y_rel_err *. x' *. y') +. (y_cst_err *. x') +. (x_rel_err *. abs (x' *. y)) +. x_cst_err *. abs y
  };
  assert {
    abs (x *. y -. exact_x *. exact_y) <=. (y_rel_err *. x' *. y') +. (x_rel_err *. (x' *. y' *. (1. +. y_rel_err) +. x' *. y_cst_err)) +. y_cst_err *. x' +. x_cst_err *. abs y
    by
      abs (x' *. y) <=. x' *. y' *. (1. +. y_rel_err) +. x' *. y_cst_err
  };
  assert {
    x_cst_err *. abs y <=. x_cst_err *. (y' *. (1. +. y_rel_err) +. y_cst_err)
  }

  use real.ExpLog

  let ghost exp_approx_err (x x_approx x' a b :real)
    requires { abs (x_approx -. x) <=. x' *. a +. b }
    requires { x <=. x' }
    ensures {
      abs (exp(x_approx) -. exp(x)) <=. exp(x) *. (exp(a *. x' +. b) -. 1.)
    }
  =
  assert {
    exp(x_approx) <=. exp(x) +. exp(x) *. (exp(a *. x' +. b) -. 1.)
    by
      exp (x_approx) <=. exp(x) *. exp (a *. x' +. b)
  };
  assert {
    exp(x_approx) >=. exp(x) -. exp(x) *. (exp(a *. x' +. b) -. 1.)
    by
      exp (x_approx) >=. exp(x) *. exp (-. a *. x' -. b)
    so
      exp(x_approx) -. exp(x) >=. exp(x) *. (exp (-. a *. x' -. b) -. 1.)
    so
      exp(a *. x' +. b) +. exp(-.a *. x' -. b) >=. 2.
    so
      -. exp(a *. x' +. b) +. 1. <=. exp(-.a *. x' -. b) -. 1.
    so
      exp(x) *. ((-. exp(a *. x' +. b)) +. 1.) <=. exp(x) *. (exp(-. a *. x' -. b) -. 1.)
    so
      -. exp(x) *. (exp(a *. x' +. b) -. 1.) <=. exp(x) *. (exp(-. a *. x' -. b) -. 1.)
  }

  let lemma log_1_minus_x (x:real)
    requires { 0. <=. abs x <. 1. }
    ensures { log (1. +. x) <=. -. log (1. -. x) }
  =
    assert { 1. +. x <=. 1. /. (1. -. x) };
    assert { forall x y z. 0. <=. x -> 0. <. y -> 0. <=. z -> x *. y <=. z -> x <=. z /. y };
    assert { exp (-.log (1. -. x)) = 1. /. (1. -. x) }

  let lemma log2_1_minus_x (x:real)
    requires { 0. <=. abs x <. 1. }
    ensures { log2 (1. +. x) <=. -. log2 (1. -. x) }
  =
  div_order_compat (log (1. +. x)) (-. log (1. -. x)) (log 2.);
  log_1_minus_x x

  let lemma log10_1_minus_x (x:real)
    requires { 0. <=. abs x <. 1. }
    ensures { log10 (1. +. x) <=. -. log10 (1. -. x) }
  =
  div_order_compat (log (1. +. x)) (-. log (1. -. x)) (log 10.);
  log_1_minus_x x

  let ghost log_approx_err (x x_approx x' a b :real)
    requires { abs (x_approx -. x) <=. x' *. a +. b }
    requires { 0. <. (x -. a *. x' -. b) }
    requires { 0. <. x <=. x' }
    ensures {
      abs (log x_approx -. log x) <=. -. log(1. -. ((a *. x' +. b) /. x))
    }
  =
    assert { a *. x' +. b  = x *. ((a *. x' +. b) /. x) };
    assert {
      log (x *. (1. -. (a *. x' +. b) /. x))
      <=. log x_approx
      <=. log (x *. (1. +. (a *. x' +. b) /.x))
      by
        0. <.  (x -. (a *. x' +. b)) <=. x_approx
    };
    log_1_minus_x ((a *. x' +. b) /. x)

  let ghost log2_approx_err (x x_approx x' a b :real)
    requires { abs (x_approx -. x) <=. x' *. a +. b }
    requires { 0. <. (x -. a *. x' -. b) }
    requires { 0. <. x <=. x' }
    ensures {
      abs (log2 x_approx -. log2 x) <=. -. log2(1. -. ((a *. x' +. b) /. x))
    }
  =
    assert { a *. x' +. b  = x *. ((a *. x' +. b) /. x) };
    assert {
      log2 (x *. (1. -. (a *. x' +. b) /. x))
      <=. log2 x_approx
      <=. log2 (x *. (1. +. (a *. x' +. b) /.x))
      by
        0. <.  (x -. (a *. x' +. b)) <=. x_approx
    };
    log2_1_minus_x ((a *. x' +. b) /. x)

  let ghost log10_approx_err (x x_approx x' a b :real)
    requires { abs (x_approx -. x) <=. x' *. a +. b }
    requires { 0. <. (x -. a *. x' -. b) }
    requires { 0. <. x <=. x' }
    ensures {
      abs (log10 x_approx -. log10 x) <=. -. log10(1. -. ((a *. x' +. b) /. x))
    }
  =
    assert { a *. x' +. b  = x *. ((a *. x' +. b) /. x) };
    assert {
      log10 (x *. (1. -. (a *. x' +. b) /. x))
      <=. log10 x_approx
      <=. log10 (x *. (1. +. (a *. x' +. b) /.x))
      by
        0. <.  (x -. (a *. x' +. b)) <=. x_approx
    };
    log10_1_minus_x ((a *. x' +. b) /. x)

  use real.Trigonometry

  lemma sin_of_approx : forall x y. abs (sin x -. sin y) <=. abs (x -. y)
  lemma cos_of_approx : forall x y. abs (cos x -. cos y) <=. abs (x -. y)

  use real.Sum
  use int.Int
  use real.FromInt

  let rec ghost sum_approx_err (rel_err cst_err:real) (f exact_f f' : int -> real) (a b:int)
    requires { a <= b }
    requires { forall i. a <= i < b -> abs (f i -. exact_f i) <=. f' i *. rel_err +. cst_err }
    variant { b - a }
    ensures { abs (sum f a b -. sum exact_f a b) <=. rel_err *. sum f' a b +. cst_err *. from_int (b-a) }
  =
  if (a < b) then
    begin
      sum_approx_err rel_err cst_err f exact_f f' a (b - 1)
    end

end

module USingleLemmas
  use real.RealInfix
  use real.FromInt
  use real.Abs
  use USingle

Single propagation lemmas

  let lemma uadd_single_error_propagation (x_uf y_uf r: usingle) (x x' x_rel_err x_cst_err y y' y_rel_err y_cst_err : real)
    requires {
      abs (to_real x_uf -. x) <=. x_rel_err *. x' +. x_cst_err
    }
    requires {
      abs (to_real y_uf -. y) <=. y_rel_err *. y' +. y_cst_err
    }
    requires { abs x <=. x' }
    requires { abs y <=. y' }
    (* TODO: Use (0 <=. x_rel_err \/ (x' = 0 /\ x_cst_err = 0)), same for y. *)
    requires { 0. <=. x_rel_err }
    requires { 0. <=. y_rel_err }
    requires { 0. <=. x_cst_err }
    requires { 0. <=. y_cst_err }
    requires { r = (x_uf ++. y_uf) }
    ensures {
      abs (to_real r -. (x +. y)) <=.
      (x_rel_err +. y_rel_err +. eps) *. (x' +. y')
          +. ((1. +. eps +. y_rel_err) *. x_cst_err +. (1. +. eps +. x_rel_err) *. y_cst_err)
    }
  =
  let ghost delta = abs (to_real (x_uf ++. y_uf) -. (to_real x_uf +. to_real y_uf)) in
  assert {
    0. <=. x_rel_err /\ 0. <=. y_rel_err ->
    delta <=.
      (eps +. y_rel_err) *. x' +. (eps +. x_rel_err) *. y'
      +. (x_rel_err +. eps) *. y_cst_err +. (y_rel_err +. eps) *. x_cst_err
    by
      (delta <=. x' *. x_rel_err +. x_cst_err +. x'
      so
        x' +. x_cst_err <=. eps *. (y' +. y_cst_err) ->
        delta <=. (eps +. x_rel_err) *. y'
        +. (eps +. y_rel_err) *. x'
        +. (y_rel_err +. eps) *. x_cst_err +. (x_rel_err +. eps) *. y_cst_err
      by
        delta <=. eps *. (y' +. y_cst_err) *. x_rel_err
              +. (eps *. (y' +. y_cst_err)))
      /\
      (delta <=. y' *. y_rel_err +. y_cst_err +. y'
      so
      abs y' +. y_cst_err <=. eps *. (x' +. x_cst_err) ->
      delta <=. (eps +. y_rel_err) *. x'
        +. (eps +. x_rel_err) *. y'
        +. (x_rel_err +. eps) *. y_cst_err +. (y_rel_err +. eps) *. x_cst_err
      by
        delta <=. eps *. (x' +. x_cst_err) *. y_rel_err
              +. (eps *. (x' +. x_cst_err)))
      /\
      (
       (eps *. (x' +. x_cst_err) <. abs y' +. y_cst_err /\
       eps *. (y' +. y_cst_err) <. abs x' +. x_cst_err) ->
       (delta <=.
       (eps +. y_rel_err) *. x' +. (eps +. x_rel_err) *. y'
      +. (x_rel_err +. eps) *. y_cst_err +. (y_rel_err +. eps) *. x_cst_err
      by
        abs (to_real x_uf +. to_real y_uf) <=.
        abs (to_real x_uf -. x) +. x' +. (abs (to_real y_uf -. y) +. y')
      so
        x' *. x_rel_err <=. (y' +. y_cst_err) /. eps *. x_rel_err /\
        y' *. y_rel_err <=. (x' +. x_cst_err) /. eps *. y_rel_err))
  }

  let lemma usub_single_error_propagation (x_uf y_uf r : usingle) (x x' x_rel_err x_cst_err y y' y_rel_err y_cst_err : real)
    requires {
      abs (to_real x_uf -. x) <=. x_rel_err *. x' +. x_cst_err
    }
    requires {
      abs (to_real y_uf -. y) <=. y_rel_err *. y' +. y_cst_err
    }
    requires { abs x <=. x' }
    requires { abs y <=. y' }
    requires { 0. <=. x_cst_err }
    requires { 0. <=. y_cst_err }
    requires { 0. <=. x_rel_err }
    requires { 0. <=. y_rel_err }
    requires { r = x_uf --. y_uf }
    ensures {
      abs (to_real r -. (x -. y))
      <=. (x_rel_err +. y_rel_err +. eps) *. (x' +. y')
          +. ((1. +. eps +. y_rel_err) *. x_cst_err +. (1. +. eps +. x_rel_err) *. y_cst_err)
    }
  = uadd_single_error_propagation x_uf (--. y_uf) r x x' x_rel_err x_cst_err (-. y) y' y_rel_err y_cst_err

  use HelperLemmas

  let lemma umul_single_error_propagation (x_uf y_uf r : usingle) (x x' x_rel_err x_cst_err y y' y_rel_err y_cst_err : real)
    requires {
      abs (to_real x_uf -. x) <=. x_rel_err *. x' +. x_cst_err
    }
    requires {
      abs (to_real y_uf -. y) <=. y_rel_err *. y' +. y_cst_err
    }
    requires { abs x <=. x' }
    requires { abs y <=. y' }
    requires { 0. <=. x_rel_err }
    requires { 0. <=. y_rel_err }
    requires { 0. <=. x_cst_err }
    requires { 0. <=. y_cst_err }
    requires { r = x_uf **. y_uf }
    ensures {
      abs (to_real r -. (x *. y)) <=.
        (eps +. (x_rel_err +. y_rel_err +. x_rel_err *. y_rel_err) *. (1. +. eps)) *. (x' *. y')
        +. (((y_cst_err +. y_cst_err *. x_rel_err) *. x'
        +. (x_cst_err +. x_cst_err *. y_rel_err) *. y'
        +. x_cst_err *. y_cst_err) *. (1. +. eps) +. eta)
    }
  =
  assert {
    to_real x_uf *. to_real y_uf -. abs (to_real x_uf *. to_real y_uf) *. eps -. eta
    <=. to_real (x_uf **. y_uf)
    <=. to_real x_uf *. to_real y_uf +. abs (to_real x_uf *. to_real y_uf) *. eps +. eta
  };
    assert { abs (x *. y) <=. x' *. y' by
       abs x *. abs y <=. x' *. abs y = abs y *. x' <=. y' *. x' };
  mult_err_combine (to_real x_uf) x x' x_rel_err x_cst_err (to_real y_uf) y y' y_rel_err y_cst_err

  use real.ExpLog

  let lemma log_single_error_propagation (log_approx : usingle -> usingle) (x_uf : usingle)
        (exact_x x' log_rel_err log_cst_err x_rel_err x_cst_err : real)
    requires { abs (to_real x_uf -. exact_x) <=. x_rel_err *. x' +. x_cst_err }
    requires {
      abs (to_real (log_approx x_uf) -. log(to_real x_uf))
       <=. log_rel_err *. abs (log (to_real x_uf)) +. log_cst_err
    }
    requires { 0. <. exact_x <=. x' }
    requires { 0. <. (exact_x -. x_rel_err *. x' -. x_cst_err) }
    requires { 0. <=. log_rel_err }
    ensures {
      abs (to_real (log_approx x_uf) -. log (exact_x))
        <=. log_rel_err *. abs (log (exact_x)) +.
          (-. log (1. -. ((x_rel_err *. x' +. x_cst_err) /. exact_x)) *. (1. +. log_rel_err)
          +. log_cst_err)
    }
  =
  log_approx_err exact_x (to_real x_uf) x' x_rel_err x_cst_err;
  assert {
   abs (log (to_real x_uf)) *. log_rel_err
    <=. (abs (log (exact_x)) -. log (1.0 -. (((x_rel_err *. x') +. x_cst_err) /. exact_x))) *. log_rel_err
  }

  let lemma log2_single_error_propagation (log2_approx : usingle -> usingle) (x_uf : usingle)
        (exact_x x' log_rel_err log_cst_err x_rel_err x_cst_err : real)
    requires { abs (to_real x_uf -. exact_x) <=. x_rel_err *. x' +. x_cst_err }
    requires {
      abs (to_real (log2_approx x_uf) -. log2(to_real x_uf))
       <=. log_rel_err *. abs (log2 (to_real x_uf)) +. log_cst_err
    }
    requires { 0. <. exact_x <=. x' }
    requires { 0. <. (exact_x -. x_rel_err *. x' -. x_cst_err) }
    requires { 0. <=. log_rel_err }
    ensures {
      abs (to_real (log2_approx x_uf) -. log2 (exact_x))
        <=. log_rel_err *. abs (log2 (exact_x)) +.
          (-. log2 (1. -. ((x_rel_err *. x' +. x_cst_err) /. exact_x)) *. (1. +. log_rel_err)
          +. log_cst_err)
    }
  =
  log2_approx_err exact_x (to_real x_uf) x' x_rel_err x_cst_err;
  assert {
   abs (log2 (to_real x_uf)) *. log_rel_err
    <=. (abs (log2 (exact_x)) -. log2 (1.0 -. (((x_rel_err *. x') +. x_cst_err) /. exact_x))) *. log_rel_err
  }

  let lemma log10_single_error_propagation (log10_approx : usingle -> usingle) (x_uf : usingle)
        (exact_x x' log_rel_err log_cst_err x_rel_err x_cst_err : real)
    requires { abs (to_real x_uf -. exact_x) <=. x_rel_err *. x' +. x_cst_err }
    requires {
      abs (to_real (log10_approx x_uf) -. log10(to_real x_uf))
       <=. log_rel_err *. abs (log10 (to_real x_uf)) +. log_cst_err
    }
    requires { 0. <. exact_x <=. x' }
    requires { 0. <. (exact_x -. x_rel_err *. x' -. x_cst_err) }
    requires { 0. <=. log_rel_err }
    ensures {
      abs (to_real (log10_approx x_uf) -. log10 (exact_x))
        <=. log_rel_err *. abs (log10 (exact_x)) +.
          (-. log10 (1. -. ((x_rel_err *. x' +. x_cst_err) /. exact_x)) *. (1. +. log_rel_err)
          +. log_cst_err)
    }
  =
  log10_approx_err exact_x (to_real x_uf) x' x_rel_err x_cst_err;
  assert {
   abs (log10 (to_real x_uf)) *. log_rel_err
    <=. (abs (log10 (exact_x)) -. log10 (1.0 -. (((x_rel_err *. x') +. x_cst_err) /. exact_x))) *. log_rel_err
  }

  let lemma exp_single_error_propagation (exp_approx : usingle -> usingle) (x_uf : usingle)
        (exact_x x' exp_rel_err exp_cst_err x_rel_err x_cst_err : real)
    requires { abs (to_real x_uf -. exact_x) <=. x_rel_err *. x' +. x_cst_err }
    requires {
      abs (to_real (exp_approx x_uf) -. exp(to_real x_uf))
        <=. exp_rel_err *. exp (to_real x_uf) +. exp_cst_err
    }
    requires { exact_x <=. x' }
    requires { 0. <=. exp_rel_err <=. 1. }
    ensures {
      abs (to_real (exp_approx x_uf) -. exp (exact_x))
      <=. (exp_rel_err +. (exp(x_rel_err *. x' +. x_cst_err) -. 1.) *. (1. +. exp_rel_err)) *. exp(exact_x)
        +. exp_cst_err
    }
  =
    exp_approx_err exact_x (to_real x_uf) x' x_rel_err x_cst_err;
    assert {
      to_real (exp_approx x_uf) <=. (exp(exact_x) +. exp(exact_x)*.(exp(x_rel_err *. x' +. x_cst_err) -. 1.))*. (1. +. exp_rel_err) +. exp_cst_err
      by
        to_real (exp_approx x_uf) <=. exp(to_real x_uf) *. (1. +. exp_rel_err) +. exp_cst_err
    };
    assert {
      to_real (exp_approx x_uf) >=.
       exp(exact_x) -.
        exp(exact_x)*.(exp_rel_err +. ((exp(x_rel_err *. x' +. x_cst_err) -. 1.) *. (1. +. exp_rel_err))) -. exp_cst_err
      by
        to_real (exp_approx x_uf) >=. exp(to_real x_uf) *. (1. -. exp_rel_err) -. exp_cst_err
      so
        exp (to_real x_uf) >=. exp exact_x -.
          (exp exact_x *. (exp ((x_rel_err *. x') +. x_cst_err) -. 1.))
      so
        to_real (exp_approx x_uf) >=.
        (exp exact_x -.
          (exp exact_x *. (exp ((x_rel_err *. x') +. x_cst_err) -. 1.)))
        *. (1. -. exp_rel_err) -. exp_cst_err
    };

  use real.Trigonometry

  let lemma sin_single_error_propagation (sin_approx : usingle -> usingle) (x_uf : usingle)
        (exact_x x' sin_rel_err sin_cst_err x_rel_err x_cst_err : real)
    requires { abs (to_real x_uf -. exact_x) <=. x_rel_err *. x' +. x_cst_err }
    requires {
      abs (to_real (sin_approx x_uf) -. sin(to_real x_uf))
        <=. sin_rel_err *. abs (sin (to_real x_uf)) +. sin_cst_err
    }
    requires { exact_x <=. x' }
    requires { 0. <=. sin_rel_err }
    ensures {
      abs (to_real (sin_approx x_uf) -. sin (exact_x))
      <=. sin_rel_err *. abs(sin(exact_x))
          +. (((x_rel_err *. x' +. x_cst_err) *. (1. +. sin_rel_err)) +. sin_cst_err)
    }
  =
  assert {
  abs (sin (to_real x_uf)) *. sin_rel_err
  <=. (abs (sin exact_x) +. (x_rel_err *. x' +. x_cst_err)) *. sin_rel_err
  }

  let lemma cos_single_error_propagation (cos_approx : usingle -> usingle) (x_uf : usingle)
        (exact_x x' cos_rel_err cos_cst_err x_rel_err x_cst_err : real)
    requires { abs (to_real x_uf -. exact_x) <=. x_rel_err *. x' +. x_cst_err }
    requires {
      abs (to_real (cos_approx x_uf) -. cos(to_real x_uf))
        <=. cos_rel_err *. abs (cos (to_real x_uf)) +. cos_cst_err
    }
    requires { exact_x <=. x' }
    requires { 0. <=. cos_rel_err }
    ensures {
      abs (to_real (cos_approx x_uf) -. cos (exact_x))
      <=. cos_rel_err *. abs(cos(exact_x))
          +. (((x_rel_err *. x' +. x_cst_err) *. (1. +. cos_rel_err)) +. cos_cst_err)
    }
  =
  assert {
  abs (cos (to_real x_uf)) *. cos_rel_err
  <=. (abs (cos exact_x) +. (x_rel_err *. x' +. x_cst_err)) *. cos_rel_err
  }

  use real.Sum
  use int.Int
  use real.FromInt

  function real_fun (f:int -> usingle) : int -> real = fun i -> to_real (f i)

  let lemma sum_single_error_propagation (x : usingle)
                (f : int -> usingle) (exact_f f' f'' : int -> real) (n:int)
                (sum_rel_err sum_cst_err f_rel_err f_cst_err : real)
    requires {
      forall i. 0 <= i < n ->
        abs ((real_fun f) i -. exact_f i) <=. f_rel_err *. f' i +. f_cst_err
    }
    requires {
      forall i. 0 <= i < n ->
      f' i -. f_rel_err *. f' i -. f_cst_err <=. f'' i <=. f' i +. f_rel_err *. f' i +. f_cst_err
    }
    requires {
      abs (to_real x -. (sum (real_fun f) 0 n))
        <=. sum_rel_err *. (sum f'' 0 n) +. sum_cst_err
    }
    requires { 0. <=. sum_rel_err }
    requires { 0 <= n }
    ensures {
      abs (to_real x -. sum exact_f 0 n)
      <=. (f_rel_err +. (sum_rel_err *. (1. +. f_rel_err))) *. sum f' 0 n +.
        ((f_cst_err *. from_int n *.(1. +. sum_rel_err)) +. sum_cst_err)
    }
  =
  sum_approx_err f_rel_err f_cst_err (real_fun f) exact_f f' 0 n;
  sum_approx_err f_rel_err f_cst_err f'' f' f' 0 n;
  assert {
    sum_rel_err *. sum f'' 0 n <=.
    sum_rel_err *. (sum f' 0 n +. ((f_rel_err *. sum f' 0 n) +. (f_cst_err *. from_int n)))
  }

  predicate is_positive_power_of_2 (x:usingle) =
    x = u1 \/ x = u2 || x = u4 || x = u8 || x = u16 || x = u32 || x = u64
    || x = u128 \/ x = u256 || x = u4096 || x = u8192 || x = u16384 || x = u32768
    || x = u65536

  axiom div_by_positive_power_of_2 : forall x y. is_positive_power_of_2 y ->
    abs (to_real (x //. y) -. to_real x /. to_real y) <=. eta

  let lemma udiv_pow_of_2_single_error_propagation (x_uf y_uf : usingle) (x x' x_rel_err x_cst_err : real)
    requires {
      abs (to_real x_uf -. x) <=. x_rel_err *. x' +. x_cst_err
    }
    requires {
      is_positive_power_of_2 y_uf
    }
    requires { abs x <=. x' }
    requires { 0. <=. x_rel_err }
    requires { 0. <=. x_cst_err }
    ensures {
      abs (to_real (x_uf //. y_uf) -. (x /. (to_real y_uf))) <=.
        x_rel_err *. (x' /. to_real y_uf) +. ((x_cst_err /. to_real y_uf) +. eta)
    }
  =
  div_order_compat (to_real x_uf) (x +. x_rel_err *. x' +. x_cst_err) (to_real y_uf);
  div_order_compat (x -. x_rel_err *. x' -. x_cst_err) (to_real x_uf) (to_real y_uf);

end

module UDoubleLemmas
  use real.RealInfix
  use real.FromInt
  use real.Abs
  use UDouble
Double propagation lemmas
  let lemma uadd_double_error_propagation (x_uf y_uf r : udouble) (x x' x_rel_err x_cst_err y y' y_rel_err y_cst_err : real)
    requires {
      abs (to_real x_uf -. x) <=. x_rel_err *. x' +. x_cst_err
    }
    requires {
      abs (to_real y_uf -. y) <=. y_rel_err *. y' +. y_cst_err
    }
    requires { abs x <=. x' }
    requires { abs y <=. y' }
    (* TODO: Use (0 <=. x_rel_err \/ (x' = 0 /\ x_cst_err = 0)), same for y. *)
    requires { 0. <=. x_rel_err }
    requires { 0. <=. y_rel_err }
    requires { 0. <=. x_cst_err }
    requires { 0. <=. y_cst_err }
    requires { r = x_uf ++. y_uf }
    ensures {
      abs (to_real r -. (x +. y)) <=.
      (x_rel_err +. y_rel_err +. eps) *. (x' +. y')
          +. ((1. +. eps +. y_rel_err) *. x_cst_err +. (1. +. eps +. x_rel_err) *. y_cst_err)
    }
  =
  let ghost delta = abs (to_real (x_uf ++. y_uf) -. (to_real x_uf +. to_real y_uf)) in
  assert {
    0. <=. x_rel_err /\ 0. <=. y_rel_err ->
    delta <=.
      (eps +. y_rel_err) *. x' +. (eps +. x_rel_err) *. y'
      +. (x_rel_err +. eps) *. y_cst_err +. (y_rel_err +. eps) *. x_cst_err
    by
      (delta <=. x' *. x_rel_err +. x_cst_err +. x'
      so
        x' +. x_cst_err <=. eps *. (y' +. y_cst_err) ->
        delta <=. (eps +. x_rel_err) *. y'
        +. (eps +. y_rel_err) *. x'
        +. (y_rel_err +. eps) *. x_cst_err +. (x_rel_err +. eps) *. y_cst_err
      by
        delta <=. eps *. (y' +. y_cst_err) *. x_rel_err
              +. (eps *. (y' +. y_cst_err)))
      /\
      (delta <=. y' *. y_rel_err +. y_cst_err +. y'
      so
      abs y' +. y_cst_err <=. eps *. (x' +. x_cst_err) ->
      delta <=. (eps +. y_rel_err) *. x'
        +. (eps +. x_rel_err) *. y'
        +. (x_rel_err +. eps) *. y_cst_err +. (y_rel_err +. eps) *. x_cst_err
      by
        delta <=. eps *. (x' +. x_cst_err) *. y_rel_err
              +. (eps *. (x' +. x_cst_err)))
      /\
      (
       (eps *. (x' +. x_cst_err) <. abs y' +. y_cst_err /\
       eps *. (y' +. y_cst_err) <. abs x' +. x_cst_err) ->
       (delta <=.
       (eps +. y_rel_err) *. x' +. (eps +. x_rel_err) *. y'
      +. (x_rel_err +. eps) *. y_cst_err +. (y_rel_err +. eps) *. x_cst_err
      by
        abs (to_real x_uf +. to_real y_uf) <=.
        abs (to_real x_uf -. x) +. x' +. (abs (to_real y_uf -. y) +. y')
      so
        x' *. x_rel_err <=. (y' +. y_cst_err) /. eps *. x_rel_err /\
        y' *. y_rel_err <=. (x' +. x_cst_err) /. eps *. y_rel_err))
  }

  let lemma usub_double_error_propagation (x_uf y_uf r : udouble) (x x' x_rel_err x_cst_err y y' y_rel_err y_cst_err : real)
    requires {
      abs (to_real x_uf -. x) <=. x_rel_err *. x' +. x_cst_err
    }
    requires {
      abs (to_real y_uf -. y) <=. y_rel_err *. y' +. y_cst_err
    }
    requires { abs x <=. x' }
    requires { abs y <=. y' }
    requires { 0. <=. x_cst_err }
    requires { 0. <=. y_cst_err }
    requires { 0. <=. x_rel_err }
    requires { 0. <=. y_rel_err }
    requires { r = x_uf --. y_uf }
    ensures {
      abs (to_real r -. (x -. y))
      <=. (x_rel_err +. y_rel_err +. eps) *. (x' +. y')
          +. ((1. +. eps +. y_rel_err) *. x_cst_err +. (1. +. eps +. x_rel_err) *. y_cst_err)
    }
  = uadd_double_error_propagation x_uf (--. y_uf) r x x' x_rel_err x_cst_err (-. y) y' y_rel_err y_cst_err

  use HelperLemmas

  let lemma umul_double_error_propagation (x_uf y_uf r : udouble) (x x' x_rel_err x_cst_err y y' y_rel_err y_cst_err : real)
    requires {
      abs (to_real x_uf -. x) <=. x_rel_err *. x' +. x_cst_err
    }
    requires {
      abs (to_real y_uf -. y) <=. y_rel_err *. y' +. y_cst_err
    }
    requires { abs x <=. x' }
    requires { abs y <=. y' }
    requires { 0. <=. x_rel_err }
    requires { 0. <=. y_rel_err }
    requires { 0. <=. x_cst_err }
    requires { 0. <=. y_cst_err }
    requires { r = x_uf **. y_uf }
    ensures {
      abs (to_real r -. (x *. y)) <=.
        (eps +. (x_rel_err +. y_rel_err +. x_rel_err *. y_rel_err) *. (1. +. eps)) *. (x' *. y')
        +. (((y_cst_err +. y_cst_err *. x_rel_err) *. x'
        +. (x_cst_err +. x_cst_err *. y_rel_err) *. y'
        +. x_cst_err *. y_cst_err) *. (1. +. eps) +. eta)
    }
  =
  assert {
    to_real x_uf *. to_real y_uf -. abs (to_real x_uf *. to_real y_uf) *. eps -. eta
    <=. to_real (x_uf **. y_uf)
    <=. to_real x_uf *. to_real y_uf +. abs (to_real x_uf *. to_real y_uf) *. eps +. eta
  };
    assert { abs (x *. y) <=. x' *. y' by
       abs x *. abs y <=. x' *. abs y = abs y *. x' <=. y' *. x' };
  mult_err_combine (to_real x_uf) x x' x_rel_err x_cst_err (to_real y_uf) y y' y_rel_err y_cst_err

  use real.ExpLog

  let lemma log_double_error_propagation (log_approx : udouble -> udouble) (x_uf : udouble)
        (exact_x x' log_rel_err log_cst_err x_rel_err x_cst_err : real)
    requires { abs (to_real x_uf -. exact_x) <=. x_rel_err *. x' +. x_cst_err }
    requires {
      abs (to_real (log_approx x_uf) -. log(to_real x_uf))
       <=. log_rel_err *. abs (log (to_real x_uf)) +. log_cst_err
    }
    requires { 0. <. exact_x <=. x' }
    requires { 0. <. (exact_x -. x_rel_err *. x' -. x_cst_err) }
    requires { 0. <=. log_rel_err }
    ensures {
      abs (to_real (log_approx x_uf) -. log (exact_x))
        <=. log_rel_err *. abs (log (exact_x)) +.
          (-. log (1. -. ((x_rel_err *. x' +. x_cst_err) /. exact_x)) *. (1. +. log_rel_err)
          +. log_cst_err)
    }
  =
  log_approx_err exact_x (to_real x_uf) x' x_rel_err x_cst_err;
  assert {
   abs (log (to_real x_uf)) *. log_rel_err
    <=. (abs (log (exact_x)) -. log (1.0 -. (((x_rel_err *. x') +. x_cst_err) /. exact_x))) *. log_rel_err
  }

  let lemma log2_double_error_propagation (log2_approx : udouble -> udouble) (x_uf : udouble)
        (exact_x x' log_rel_err log_cst_err x_rel_err x_cst_err : real)
    requires { abs (to_real x_uf -. exact_x) <=. x_rel_err *. x' +. x_cst_err }
    requires {
      abs (to_real (log2_approx x_uf) -. log2(to_real x_uf))
       <=. log_rel_err *. abs (log2 (to_real x_uf)) +. log_cst_err
    }
    requires { 0. <. exact_x <=. x' }
    requires { 0. <. (exact_x -. x_rel_err *. x' -. x_cst_err) }
    requires { 0. <=. log_rel_err }
    ensures {
      abs (to_real (log2_approx x_uf) -. log2 (exact_x))
        <=. log_rel_err *. abs (log2 (exact_x)) +.
          (-. log2 (1. -. ((x_rel_err *. x' +. x_cst_err) /. exact_x)) *. (1. +. log_rel_err)
          +. log_cst_err)
    }
  =
  log2_approx_err exact_x (to_real x_uf) x' x_rel_err x_cst_err;
  assert {
   abs (log2 (to_real x_uf)) *. log_rel_err
    <=. (abs (log2 (exact_x)) -. log2 (1.0 -. (((x_rel_err *. x') +. x_cst_err) /. exact_x))) *. log_rel_err
  }

  let lemma log10_double_error_propagation (log10_approx : udouble -> udouble) (x_uf : udouble)
        (exact_x x' log_rel_err log_cst_err x_rel_err x_cst_err : real)
    requires { abs (to_real x_uf -. exact_x) <=. x_rel_err *. x' +. x_cst_err }
    requires {
      abs (to_real (log10_approx x_uf) -. log10(to_real x_uf))
       <=. log_rel_err *. abs (log10 (to_real x_uf)) +. log_cst_err
    }
    requires { 0. <. exact_x <=. x' }
    requires { 0. <. (exact_x -. x_rel_err *. x' -. x_cst_err) }
    requires { 0. <=. log_rel_err }
    ensures {
      abs (to_real (log10_approx x_uf) -. log10 (exact_x))
        <=. log_rel_err *. abs (log10 (exact_x)) +.
          (-. log10 (1. -. ((x_rel_err *. x' +. x_cst_err) /. exact_x)) *. (1. +. log_rel_err)
          +. log_cst_err)
    }
  =
  log10_approx_err exact_x (to_real x_uf) x' x_rel_err x_cst_err;
  assert {
   abs (log10 (to_real x_uf)) *. log_rel_err
    <=. (abs (log10 (exact_x)) -. log10 (1.0 -. (((x_rel_err *. x') +. x_cst_err) /. exact_x))) *. log_rel_err
  }

  let lemma exp_double_error_propagation (exp_approx : udouble -> udouble) (x_uf : udouble)
        (exact_x x' exp_rel_err exp_cst_err x_rel_err x_cst_err : real)
    requires { abs (to_real x_uf -. exact_x) <=. x_rel_err *. x' +. x_cst_err }
    requires {
      abs (to_real (exp_approx x_uf) -. exp(to_real x_uf))
        <=. exp_rel_err *. exp (to_real x_uf) +. exp_cst_err
    }
    requires { exact_x <=. x' }
    requires { 0. <=. exp_rel_err <=. 1. }
    ensures {
      abs (to_real (exp_approx x_uf) -. exp (exact_x))
      <=. (exp_rel_err +. (exp(x_rel_err *. x' +. x_cst_err) -. 1.) *. (1. +. exp_rel_err)) *. exp(exact_x)
        +. exp_cst_err
    }
  =
    exp_approx_err exact_x (to_real x_uf) x' x_rel_err x_cst_err;
    assert {
      to_real (exp_approx x_uf) <=. (exp(exact_x) +. exp(exact_x)*.(exp(x_rel_err *. x' +. x_cst_err) -. 1.))*. (1. +. exp_rel_err) +. exp_cst_err
      by
        to_real (exp_approx x_uf) <=. exp(to_real x_uf) *. (1. +. exp_rel_err) +. exp_cst_err
    };
    assert {
      to_real (exp_approx x_uf) >=.
       exp(exact_x) -.
        exp(exact_x)*.(exp_rel_err +. ((exp(x_rel_err *. x' +. x_cst_err) -. 1.) *. (1. +. exp_rel_err))) -. exp_cst_err
      by
        to_real (exp_approx x_uf) >=. exp(to_real x_uf) *. (1. -. exp_rel_err) -. exp_cst_err
      so
        exp (to_real x_uf) >=. exp exact_x -.
          (exp exact_x *. (exp ((x_rel_err *. x') +. x_cst_err) -. 1.0))
      so
        to_real (exp_approx x_uf) >=.
        (exp exact_x -.
          (exp exact_x *. (exp ((x_rel_err *. x') +. x_cst_err) -. 1.0)))
        *. (1. -. exp_rel_err) -. exp_cst_err
    };

  use real.Trigonometry

  let lemma sin_double_error_propagation (sin_approx : udouble -> udouble) (x_uf : udouble)
        (exact_x x' sin_rel_err sin_cst_err x_rel_err x_cst_err : real)
    requires { abs (to_real x_uf -. exact_x) <=. x_rel_err *. x' +. x_cst_err }
    requires {
      abs (to_real (sin_approx x_uf) -. sin(to_real x_uf))
        <=. sin_rel_err *. abs (sin (to_real x_uf)) +. sin_cst_err
    }
    requires { exact_x <=. x' }
    requires { 0. <=. sin_rel_err }
    ensures {
      abs (to_real (sin_approx x_uf) -. sin (exact_x))
      <=. sin_rel_err *. abs(sin(exact_x))
          +. (((x_rel_err *. x' +. x_cst_err) *. (1. +. sin_rel_err)) +. sin_cst_err)
    }
  =
  assert {
  abs (sin (to_real x_uf)) *. sin_rel_err
  <=. (abs (sin exact_x) +. (x_rel_err *. x' +. x_cst_err)) *. sin_rel_err
  }

  let lemma cos_double_error_propagation (cos_approx : udouble -> udouble) (x_uf : udouble)
        (exact_x x' cos_rel_err cos_cst_err x_rel_err x_cst_err : real)
    requires { abs (to_real x_uf -. exact_x) <=. x_rel_err *. x' +. x_cst_err }
    requires {
      abs (to_real (cos_approx x_uf) -. cos(to_real x_uf))
        <=. cos_rel_err *. abs (cos (to_real x_uf)) +. cos_cst_err
    }
    requires { exact_x <=. x' }
    requires { 0. <=. cos_rel_err }
    ensures {
      abs (to_real (cos_approx x_uf) -. cos (exact_x))
      <=. cos_rel_err *. abs(cos(exact_x))
          +. (((x_rel_err *. x' +. x_cst_err) *. (1. +. cos_rel_err)) +. cos_cst_err)
    }
  =
  assert {
  abs (cos (to_real x_uf)) *. cos_rel_err
  <=. (abs (cos exact_x) +. (x_rel_err *. x' +. x_cst_err)) *. cos_rel_err
  }

  use real.Sum
  use int.Int
  use real.FromInt

  function real_fun (f:int -> udouble) : int -> real = fun i -> to_real (f i)

  let lemma sum_double_error_propagation (x : udouble)
                (f : int -> udouble) (exact_f f' f'' : int -> real) (n:int)
                (sum_rel_err sum_cst_err f_rel_err f_cst_err : real)
    requires {
      forall i. 0 <= i < n ->
        abs ((real_fun f) i -. exact_f i) <=. f_rel_err *. f' i +. f_cst_err
    }
    requires {
      forall i. 0 <= i < n ->
      f' i -. f_rel_err *. f' i -. f_cst_err <=. f'' i <=. f' i +. f_rel_err *. f' i +. f_cst_err
    }
    requires {
      abs (to_real x -. (sum (real_fun f) 0 n))
        <=. sum_rel_err *. (sum f'' 0 n) +. sum_cst_err
    }
    requires { 0. <=. sum_rel_err }
    requires { 0 <= n }
    ensures {
      abs (to_real x -. sum exact_f 0 n)
      <=. (f_rel_err +. (sum_rel_err *. (1. +. f_rel_err))) *. sum f' 0 n +.
        ((f_cst_err *. from_int n *.(1. +. sum_rel_err)) +. sum_cst_err)
    }
  =
  sum_approx_err f_rel_err f_cst_err (real_fun f) exact_f f' 0 n;
  sum_approx_err f_rel_err f_cst_err f'' f' f' 0 n;
  assert {
    sum_rel_err *. sum f'' 0 n <=.
    sum_rel_err *. (sum f' 0 n +. ((f_rel_err *. sum f' 0 n) +. (f_cst_err *. from_int n)))
  }

  predicate is_positive_power_of_2 (x:udouble) =
    x = u1 \/ x = u2 || x = u4 || x = u8 || x = u16 || x = u32 || x = u64
    || x = u128 \/ x = u256 || x = u4096 || x = u8192 || x = u16384 || x = u32768
    || x = u65536

  axiom div_by_positive_power_of_2 : forall x y. is_positive_power_of_2 y ->
    abs (to_real (x //. y) -. to_real x /. to_real y) <=. eta

  let lemma udiv_pow_of_2_double_error_propagation (x_uf y_uf : udouble) (x x' x_rel_err x_cst_err : real)
    requires {
      abs (to_real x_uf -. x) <=. x_rel_err *. x' +. x_cst_err
    }
    requires {
      is_positive_power_of_2 y_uf
    }
    requires { abs x <=. x' }
    requires { 0. <=. x_rel_err }
    requires { 0. <=. x_cst_err }
    ensures {
      abs (to_real (x_uf //. y_uf) -. (x /. (to_real y_uf))) <=.
        x_rel_err *. (x' /. to_real y_uf) +. ((x_cst_err /. to_real y_uf) +. eta)
    }
  =
  div_order_compat (to_real x_uf) (x +. x_rel_err *. x' +. x_cst_err) (to_real y_uf);
  div_order_compat (x -. x_rel_err *. x' -. x_cst_err) (to_real x_uf) (to_real y_uf);

end

Generated by why3doc 1.7.2+git