Source file stdcompat__float.ml

 include Float

(*

module Array = struct
  type t = Stdcompat__init.floatarray

  include Stdcompat__array.Floatarray

  let unsafe_fill array offset len v =
    for i = offset to offset + len - 1 do
      unsafe_set array i v
    done

  let unsafe_blit src src_offset dst dst_offset len =
    for i = 0 to len - 1 do
      unsafe_set dst (dst_offset + i) (unsafe_get src (src_offset + i))
    done

  let make len v =
    let result = create len in
    unsafe_fill result 0 len v;
    result

  let init len f =
    let result = create len in
    for i = 0 to len - 1 do
      unsafe_set result i (f i)
    done;
    result

  let append a1 a2 =
    let l1 = length a1 in
    let l2 = length a2 in
    let result = create (l1 + l2) in
    unsafe_blit a1 0 result 0 l1;
    unsafe_blit a2 0 result l1 l2;
    result

  let concat l =
    let len = List.fold_left (fun acc a -> acc + length a) 0 l in
    let result = create len in
    let rec loop l i =
      match l with
      | [] -> ()
      | hd :: tl ->
          let len = length hd in
          unsafe_blit hd 0 result i len;
          loop tl (i + len) in
    loop l 0;
    result

  let check a offset len msg =
    if offset < 0 || len < 0 || offset + len < 0 || offset + len > length a then
      failwith msg

  let sub a offset len =
    check a offset len "Float.Array.sub";
    let result = create len in
    unsafe_blit a offset result 0 len;
    result

  let copy a =
    let len = length a in
    let result = create len in
    unsafe_blit a 0 result 0 len;
    result

  let fill a offset len v =
    check a offset len "Float.Array.fill";
    unsafe_fill a offset len v

  let blit src src_offset dst dst_offset len =
    check src src_offset len "Float.Array.blit";
    check dst dst_offset len "Float.Array.blit";
    unsafe_blit dst src_offset dst dst_offset len

  let to_list a =
    Stdcompat__list.init (length a) (unsafe_get a)

  let of_list l =
    let result = create (List.length l) in
    let rec fill i l =
      match l with
      | [] -> result
      | hd :: tl ->
          unsafe_set result i hd;
          fill (succ i) tl in
    fill 0 l

  let iter f a =
    for i = 0 to length a - 1 do
      f (unsafe_get a i)
    done

  let iteri f a =
    for i = 0 to length a - 1 do
      f i (unsafe_get a i)
    done

  let map f a =
    let len = length a in
    let result = create len in
    for i = 0 to len - 1 do
      unsafe_set result i (f (unsafe_get a i))
    done;
    result

  let mapi f a =
    let len = length a in
    let result = create len in
    for i = 0 to len - 1 do
      unsafe_set result i (f i (unsafe_get a i))
    done;
    result

  let fold_left f acc a =
    let len = length a in
    let rec aux i acc =
      if i < len then
        aux (succ i) (f acc (unsafe_get a i))
      else
        acc in
    aux 0 acc

  let fold_right f a acc =
    let len = length a in
    let rec aux i acc =
      if i > 0 then
        let j = pred i in
        aux j (f (unsafe_get a j) acc)
      else
        acc in
    aux len acc

  let iter2 f a b =
    let len = length a in
    if len <> length b then
      failwith "Float.Array.iter2";
    for i = 0 to len - 1 do
      f (unsafe_get a i) (unsafe_get b i)
    done

  let map2 f a b =
    let len = length a in
    if len <> length b then
      failwith "Float.Array.iter2";
    let result = create len in
    for i = 0 to len - 1 do
      unsafe_set result i (f (unsafe_get a i) (unsafe_get b i))
    done;
    result

  let for_all f a =
    let len = length a in
    let rec aux i =
      i >= len || f (unsafe_get a i) && aux (succ i) in
    aux 0

  let exists f a =
    let len = length a in
    let rec aux i =
      i < len && (f (unsafe_get a i) || aux (succ i)) in
    aux 0

  let mem x a =
    exists (fun y -> compare x y = 0) a

  let mem_ieee x a =
    exists (( = ) x) a

  let sort p a =
    failwith "unimplemented"

  let stable_sort = sort

  let fast_sort = sort

  let to_seq a =
    let rec aux i () =
      if i < length a then
        let x = unsafe_get a i in
        Stdcompat__seq.Cons (x, aux (i + 1))
      else
        Stdcompat__seq.Nil in
    aux 0

  let to_seqi a =
    let rec aux i () =
      if i < length a then
        let x = unsafe_get a i in
        Stdcompat__seq.Cons ((i, x), aux (i + 1))
      else
        Stdcompat__seq.Nil in
    aux 0

  let of_rev_list l =
    let len = List.length l in
    let result = create len in
    let rec fill i l =
      match l with
      | [] -> result
      | hd :: tl ->
          unsafe_set result i hd;
          fill (pred i) tl in
    fill (len - 1) l

  let of_seq i =
    let l = Stdcompat__seq.fold_left (fun acc x -> x :: acc) [] i in
    of_rev_list l

  let map_to_array f a =
    Array.init (length a) (fun i -> f (unsafe_get a i))

  let map_from_array f a =
    init (Array.length a) (fun i -> f (Array.unsafe_get a i))
end


 include (Float : module type of Float with module Array := Array)

*)
(*
type fpclass = Pervasives.fpclass =
  | FP_normal 
  | FP_subnormal 
  | FP_zero 
  | FP_infinite 
  | FP_nan 
(** @since 4.07.0:
    type fpclass = Pervasives.fpclass =
      | FP_normal 
      | FP_subnormal 
      | FP_zero 
      | FP_infinite 
      | FP_nan 
 *)
type t = float
(** @since 4.07.0: type t = float
 *)
let infinity = Pervasives.infinity
(** @since 4.07.0: val infinity : float
 *)
let neg_infinity = Pervasives.neg_infinity
(** @since 4.07.0: val neg_infinity : float
 *)
let nan = Pervasives.nan
(** @since 4.07.0: val nan : float
 *)

let pi = 0x1.921fb54442d18p+1

(*
let pi = 4. *. atan 1.
*)
(** @since 4.07.0: val pi : float
 *)
let max_float = Pervasives.max_float
(** @since 4.07.0: val max_float : float
 *)
let min_float = Pervasives.min_float
(** @since 4.07.0: val min_float : float
 *)
let epsilon = Pervasives.epsilon_float
(** @since 4.07.0: val epsilon : float
 *)
let of_string_opt = Stdcompat__pervasives.float_of_string_opt
(** @since 4.07.0: val of_string_opt : string -> float option
 *)
let to_string = Pervasives.string_of_float
(** @since 4.07.0: val to_string : float -> string
 *)
external compare : float -> float -> int = "%compare"
(** @since 4.07.0: val compare : t -> t -> int
 *)
let equal x y = compare x y = 0
(** @since 4.07.0: val equal : t -> t -> bool
 *)
let hash = Hashtbl.hash
(** @since 4.07.0: val hash : t -> int
 *)

external hypot : float -> float -> float = "caml_hypot_float" "caml_hypot"
[@@unboxed ][@@noalloc ]

(*

external hypot :
  float -> float -> float = "caml_hypot_float" "caml_hypot" "float"

(*
let hypot = Stdcompat__pervasives.hypot
*)

*)
(** @since 4.03.0:
    external hypot :
      float -> float -> float = "caml_hypot_float" "caml_hypot"[@@unboxed ]
    [@@noalloc ]
@since 4.00.0:
external hypot :
  float -> float -> float = "caml_hypot_float" "caml_hypot" "float"
 *)

external expm1 : float -> float = "caml_expm1_float" "caml_expm1"[@@unboxed ]
[@@noalloc ]

(*

external expm1 : float -> float = "caml_expm1_float" "caml_expm1" "float"

(*
let expm1 = Stdcompat__pervasives.expm1
*)

*)
(** @since 4.03.0:
    external expm1 : float -> float = "caml_expm1_float" "caml_expm1"
    [@@unboxed ][@@noalloc ]
@since 3.12.0:
external expm1 : float -> float = "caml_expm1_float" "caml_expm1" "float"
 *)

external log1p : float -> float = "caml_log1p_float" "caml_log1p"[@@unboxed ]
[@@noalloc ]

(*

external log1p : float -> float = "caml_log1p_float" "caml_log1p" "float"

(*
let log1p = Stdcompat__pervasives.log1p
*)

*)
(** @since 4.03.0:
    external log1p : float -> float = "caml_log1p_float" "caml_log1p"
    [@@unboxed ][@@noalloc ]
@since 3.12.0:
external log1p : float -> float = "caml_log1p_float" "caml_log1p" "float"
 *)
external neg : float -> float = "%negfloat"
(** Alias for {!Float.neg} *)
external add : float -> float -> float = "%addfloat"
(** Alias for {!Float.add} *)
external sub : float -> float -> float = "%subfloat"
(** Alias for {!Float.sub} *)
external mul : float -> float -> float = "%mulfloat"
(** Alias for {!Float.mul} *)
external div : float -> float -> float = "%divfloat"
(** Alias for {!Float.div} *)

external rem : float -> float -> float = "caml_fmod_float" "fmod"[@@unboxed ]
[@@noalloc ]

(*

external rem : float -> float -> float = "caml_fmod_float" "fmod" "float"

(*
external rem : float -> float -> float = "fmod_float" "fmod" "float"
*)
*)
(** @since 4.03.0:
    external rem : float -> float -> float = "caml_fmod_float" "fmod"
    [@@unboxed ][@@noalloc ]
@since 3.11.0:
external rem : float -> float -> float = "caml_fmod_float" "fmod" "float"
 *)
external abs : float -> float = "%absfloat"
(** Alias for {!Float.abs} *)
external of_int : int -> float = "%floatofint"
(** Alias for {!Float.of_int} *)
external to_int : float -> int = "%intoffloat"
(** Alias for {!Float.to_int} *)

external of_string : string -> float = "caml_float_of_string"

(*
external of_string : string -> float = "float_of_string"
*)
(** Alias for {!Float.of_string} *)

external classify_float :
  ((float)[@unboxed ]) -> fpclass = "caml_classify_float"
    "caml_classify_float_unboxed"[@@noalloc ]

(*

external classify_float : float -> fpclass = "caml_classify_float"

(*
external classify_float : float -> fpclass = "classify_float"
*)
*)
(** @since 4.03.0:
    external classify_float :
      ((float)[@unboxed ]) -> fpclass = "caml_classify_float"
        "caml_classify_float_unboxed"[@@noalloc ]
@since 3.11.0:
external classify_float : float -> fpclass = "caml_classify_float"
 *)

external pow : float -> float -> float = "caml_power_float" "pow"[@@unboxed ]
[@@noalloc ]

(*

external pow : float -> float -> float = "caml_power_float" "pow" "float"

(*
external pow : float -> float -> float = "power_float" "pow" "float"
*)
*)
(** @since 4.03.0:
    external pow : float -> float -> float = "caml_power_float" "pow"
    [@@unboxed ][@@noalloc ]
@since 3.11.0:
external pow : float -> float -> float = "caml_power_float" "pow" "float"
 *)

external sqrt : float -> float = "caml_sqrt_float" "sqrt"[@@unboxed ]
[@@noalloc ]

(*

external sqrt : float -> float = "caml_sqrt_float" "sqrt" "float"

(*
external sqrt : float -> float = "sqrt_float" "sqrt" "float"
*)
*)
(** @since 4.03.0:
    external sqrt : float -> float = "caml_sqrt_float" "sqrt"[@@unboxed ]
    [@@noalloc ]
@since 3.11.0:
external sqrt : float -> float = "caml_sqrt_float" "sqrt" "float"
 *)

external exp : float -> float = "caml_exp_float" "exp"[@@unboxed ][@@noalloc
                                                                    ]

(*

external exp : float -> float = "caml_exp_float" "exp" "float"

(*
external exp : float -> float = "exp_float" "exp" "float"
*)
*)
(** @since 4.03.0:
    external exp : float -> float = "caml_exp_float" "exp"[@@unboxed ]
    [@@noalloc ]
@since 3.11.0: external exp : float -> float = "caml_exp_float" "exp" "float"
 *)

external log : float -> float = "caml_log_float" "log"[@@unboxed ][@@noalloc
                                                                    ]

(*

external log : float -> float = "caml_log_float" "log" "float"

(*
external log : float -> float = "log_float" "log" "float"
*)
*)
(** @since 4.03.0:
    external log : float -> float = "caml_log_float" "log"[@@unboxed ]
    [@@noalloc ]
@since 3.11.0: external log : float -> float = "caml_log_float" "log" "float"
 *)

external log10 : float -> float = "caml_log10_float" "log10"[@@unboxed ]
[@@noalloc ]

(*

external log10 : float -> float = "caml_log10_float" "log10" "float"

(*
external log10 : float -> float = "log10_float" "log10" "float"
*)
*)
(** @since 4.03.0:
    external log10 : float -> float = "caml_log10_float" "log10"[@@unboxed ]
    [@@noalloc ]
@since 3.11.0:
external log10 : float -> float = "caml_log10_float" "log10" "float"
 *)

external cos : float -> float = "caml_cos_float" "cos"[@@unboxed ][@@noalloc
                                                                    ]

(*

external cos : float -> float = "caml_cos_float" "cos" "float"

(*
external cos : float -> float = "cos_float" "cos" "float"
*)
*)
(** @since 4.03.0:
    external cos : float -> float = "caml_cos_float" "cos"[@@unboxed ]
    [@@noalloc ]
@since 3.11.0: external cos : float -> float = "caml_cos_float" "cos" "float"
 *)

external sin : float -> float = "caml_sin_float" "sin"[@@unboxed ][@@noalloc
                                                                    ]

(*

external sin : float -> float = "caml_sin_float" "sin" "float"

(*
external sin : float -> float = "sin_float" "sin" "float"
*)
*)
(** @since 4.03.0:
    external sin : float -> float = "caml_sin_float" "sin"[@@unboxed ]
    [@@noalloc ]
@since 3.11.0: external sin : float -> float = "caml_sin_float" "sin" "float"
 *)

external tan : float -> float = "caml_tan_float" "tan"[@@unboxed ][@@noalloc
                                                                    ]

(*

external tan : float -> float = "caml_tan_float" "tan" "float"

(*
external tan : float -> float = "tan_float" "tan" "float"
*)
*)
(** @since 4.03.0:
    external tan : float -> float = "caml_tan_float" "tan"[@@unboxed ]
    [@@noalloc ]
@since 3.11.0: external tan : float -> float = "caml_tan_float" "tan" "float"
 *)

external acos : float -> float = "caml_acos_float" "acos"[@@unboxed ]
[@@noalloc ]

(*

external acos : float -> float = "caml_acos_float" "acos" "float"

(*
external acos : float -> float = "acos_float" "acos" "float"
*)
*)
(** @since 4.03.0:
    external acos : float -> float = "caml_acos_float" "acos"[@@unboxed ]
    [@@noalloc ]
@since 3.11.0:
external acos : float -> float = "caml_acos_float" "acos" "float"
 *)

external asin : float -> float = "caml_asin_float" "asin"[@@unboxed ]
[@@noalloc ]

(*

external asin : float -> float = "caml_asin_float" "asin" "float"

(*
external asin : float -> float = "asin_float" "asin" "float"
*)
*)
(** @since 4.03.0:
    external asin : float -> float = "caml_asin_float" "asin"[@@unboxed ]
    [@@noalloc ]
@since 3.11.0:
external asin : float -> float = "caml_asin_float" "asin" "float"
 *)

external atan : float -> float = "caml_atan_float" "atan"[@@unboxed ]
[@@noalloc ]

(*

external atan : float -> float = "caml_atan_float" "atan" "float"

(*
external atan : float -> float = "atan_float" "atan" "float"
*)
*)
(** @since 4.03.0:
    external atan : float -> float = "caml_atan_float" "atan"[@@unboxed ]
    [@@noalloc ]
@since 3.11.0:
external atan : float -> float = "caml_atan_float" "atan" "float"
 *)

external atan2 : float -> float -> float = "caml_atan2_float" "atan2"
[@@unboxed ][@@noalloc ]

(*

external atan2 : float -> float -> float = "caml_atan2_float" "atan2" "float"

(*
external atan2 : float -> float -> float = "atan2_float" "atan2" "float"
*)
*)
(** @since 4.03.0:
    external atan2 : float -> float -> float = "caml_atan2_float" "atan2"
    [@@unboxed ][@@noalloc ]
@since 3.11.0:
external atan2 : float -> float -> float = "caml_atan2_float" "atan2" "float"
 *)

external cosh : float -> float = "caml_cosh_float" "cosh"[@@unboxed ]
[@@noalloc ]

(*

external cosh : float -> float = "caml_cosh_float" "cosh" "float"

(*
external cosh : float -> float = "cosh_float" "cosh" "float"
*)
*)
(** @since 4.03.0:
    external cosh : float -> float = "caml_cosh_float" "cosh"[@@unboxed ]
    [@@noalloc ]
@since 3.11.0:
external cosh : float -> float = "caml_cosh_float" "cosh" "float"
 *)

external sinh : float -> float = "caml_sinh_float" "sinh"[@@unboxed ]
[@@noalloc ]

(*

external sinh : float -> float = "caml_sinh_float" "sinh" "float"

(*
external sinh : float -> float = "sinh_float" "sinh" "float"
*)
*)
(** @since 4.03.0:
    external sinh : float -> float = "caml_sinh_float" "sinh"[@@unboxed ]
    [@@noalloc ]
@since 3.11.0:
external sinh : float -> float = "caml_sinh_float" "sinh" "float"
 *)

external tanh : float -> float = "caml_tanh_float" "tanh"[@@unboxed ]
[@@noalloc ]

(*

external tanh : float -> float = "caml_tanh_float" "tanh" "float"

(*
external tanh : float -> float = "tanh_float" "tanh" "float"
*)
*)
(** @since 4.03.0:
    external tanh : float -> float = "caml_tanh_float" "tanh"[@@unboxed ]
    [@@noalloc ]
@since 3.11.0:
external tanh : float -> float = "caml_tanh_float" "tanh" "float"
 *)

external ceil : float -> float = "caml_ceil_float" "ceil"[@@unboxed ]
[@@noalloc ]

(*

external ceil : float -> float = "caml_ceil_float" "ceil" "float"

(*
external ceil : float -> float = "ceil_float" "ceil" "float"
*)
*)
(** @since 4.03.0:
    external ceil : float -> float = "caml_ceil_float" "ceil"[@@unboxed ]
    [@@noalloc ]
@since 3.11.0:
external ceil : float -> float = "caml_ceil_float" "ceil" "float"
 *)

external floor : float -> float = "caml_floor_float" "floor"[@@unboxed ]
[@@noalloc ]

(*

external floor : float -> float = "caml_floor_float" "floor" "float"

(*
external floor : float -> float = "floor_float" "floor" "float"
*)
*)
(** @since 4.03.0:
    external floor : float -> float = "caml_floor_float" "floor"[@@unboxed ]
    [@@noalloc ]
@since 3.11.0:
external floor : float -> float = "caml_floor_float" "floor" "float"
 *)

external frexp : float -> (float * int) = "caml_frexp_float"

(*
external frexp : float -> (float * int) = "frexp_float"
*)
(** Alias for {!Float.frexp} *)

external ldexp :
  ((float)[@unboxed ]) -> ((int)[@untagged ]) -> ((float)[@unboxed ]) =
    "caml_ldexp_float" "caml_ldexp_float_unboxed"[@@noalloc ]

(*

external ldexp : float -> int -> float = "caml_ldexp_float"

(*
external ldexp : float -> int -> float = "ldexp_float"
*)
*)
(** @since 4.03.0:
    external ldexp :
      ((float)[@unboxed ]) -> ((int)[@untagged ]) -> ((float)[@unboxed ]) =
        "caml_ldexp_float" "caml_ldexp_float_unboxed"[@@noalloc ]
@since 3.11.0: external ldexp : float -> int -> float = "caml_ldexp_float"
 *)

external modf : float -> (float * float) = "caml_modf_float"

(*
external modf : float -> (float * float) = "modf_float"
*)
(** Alias for {!Float.modf} *)
*)

(*
let sign_of_bits bits =
  not (Stdcompat__int64.equal (Int64.shift_right_logical bits 63) Int64.zero)

let sign_bit x = 
  sign_of_bits (Int64.bits_of_float x)

let is_nan (x : float) =
  x <> x

let is_infinite x =
  1. /. x = 0.

let is_finite x =
  x -. x = 0.

let trunc x =
  if sign_bit x then
    ceil x
  else
    floor x

let is_integer x =
  x = trunc x && is_finite x

let round x =
  if is_finite x then
    let ceil = ceil x in
    let floor = floor x in
    let o = compare (x -. floor) (ceil -. x) in
    if o < 0 then
      floor
    else if o > 0 then
      ceil
    else
      if sign_bit x then
        floor
      else
        ceil
  else
    x

let min_max_num x y =
  if is_nan x then (y, y)
  else if is_nan y then (x, x)
  else if x <= y && (sign_bit x || not (sign_bit y)) then (x, y)
  else (y, x)

let max_num x y =
  if is_nan x then y
  else if is_nan y then x
  else if x <= y && (sign_bit x || not (sign_bit y)) then y
  else x

let min_num x y =
  if is_nan x then y
  else if is_nan y then x
  else if x <= y && (sign_bit x || not (sign_bit y)) then x
  else y

let min_max x y =
  if is_nan x then (x, x)
  else if is_nan y then (y, y)
  else if x <= y && (sign_bit x || not (sign_bit y)) then (x, y)
  else (y, x)

let max x y =
  if is_nan x then x
  else if is_nan y then y
  else if x <= y && (sign_bit x || not (sign_bit y)) then y
  else x

let min x y =
  if is_nan x then x
  else if is_nan y then y
  else if x <= y && (sign_bit x || not (sign_bit y)) then x
  else y


external copy_sign :
  float -> float -> float = "caml_copysign_float" "caml_copysign"[@@unboxed ]
[@@noalloc ]

(*

external copy_sign :
  float -> float -> float = "caml_copysign_float" "caml_copysign" "float"

(*
let copy_sign = Stdcompat__pervasives.copysign
*)
*)

let pred_ x =
  let bits = Int64.bits_of_float x in
  if is_finite x then
    let bits' =
      if sign_of_bits bits then
        Int64.succ bits
      else
        Int64.pred bits in
    Int64.float_of_bits bits'
  else
    if sign_of_bits bits then
      x
    else
      Int64.float_of_bits (Int64.pred bits)

let succ_ x =
  let bits = Int64.bits_of_float x in
  if is_finite x then
    let bits' =
      if sign_of_bits bits then
        Int64.pred bits
      else
        Int64.succ bits in
    Int64.float_of_bits bits'
  else
    if sign_of_bits bits then
      Int64.float_of_bits (Int64.pred bits)
    else
      x

let next_after x y =
  if is_nan x || is_nan y then nan
  else
    let o = compare x y in
    if o < 0 then
      succ_ x
    else if o > 0 then
      pred_ x
    else
      y

let pred x =
  if is_nan x then
    x
  else
    pred_ x

let succ x =
  if is_nan x then
    x
  else
    succ_ x

let fma x y z =
  x *. y +. z

let minus_one = -1.

let one = 1.

let zero = 0.

module ArrayLabels = struct
  include Array

  let init i ~f = init i f

  let sub a ~pos ~len = sub a pos len

  let fill a ~pos ~len v = fill a pos len v

  let blit ~src ~src_pos ~dst ~dst_pos ~len = blit src src_pos dst dst_pos len

  let iter ~f a = iter f a

  let iteri ~f a = iteri f a

  let map ~f a = map f a

  let mapi ~f a = mapi f a

  let fold_left ~f ~init a = fold_left f init a

  let fold_right ~f a ~init = fold_right f a init

  let iter2 ~f a = iter2 f a

  let map2 ~f a = map2 f a

  let for_all ~f a = for_all f a

  let exists ~f a = exists f a

  let mem v ~set = mem v set

  let mem_ieee v ~set = mem_ieee v set

  let sort ~cmp a = sort cmp a

  let stable_sort ~cmp a = stable_sort cmp a

  let fast_sort ~cmp a = fast_sort cmp a

  let map_to_array ~f a = map_to_array f a

  let map_from_array ~f a = map_from_array f a
end
*)


let cbrt x =
  let res = exp ((1. /. 3.) *. log (abs x)) in
  if x >= 0. then
    res
  else
    -. res

let log2e = log 2.

let log2 x =
  log x /. log2e

let exp2 x =
  pow 2. x

let acosh x =
  log (x +. sqrt (x *. x -. 1.))

let asinh x =
  log (x +. sqrt (x *. x +. 1.))

let atanh x =
  0.5 *. log ((1. +. x) /. (1. -. x))

let erf x =
  (* from ocaml source floats.c, based on formula 7.1.26
     from the "Handbook of Mathematical Functions" by Abramowitz and Stegun. *)
  let a1 =  0.254829592 and a2 = -0.284496736 and a3 = 1.421413741
  and a4 = -1.453152027 and a5 = 1.061405429 and p = 0.3275911 in
  let x' = abs x in
  let t = 1. /. (1. +. p *. x') in
  let y =
    1. -. ((((a5 *. t +. a4) *. t +. a3) *. t +. a2) *. t +. a1) *.
      t *. exp (-. x' *. x') in
  if x >= 0. then y
  else -. y

let erfc x =
  1. -. erf x