Merge pull request #5 from bobot/qqbar

Add qqbar from_roots and from_enclosure

.ocamlformat

@@ -1 +1 @@
-version=0.24.1
+version=0.25.1

antic.opam

@@ -31,3 +31,4 @@ build: [
   ]
 ]
 dev-repo: "git+https://github.com/bobot/ocaml-flint.git"
+conflicts: [ "ocaml-option-bytecode-only" ]

arb.opam

@@ -31,3 +31,4 @@ build: [
   ]
 ]
 dev-repo: "git+https://github.com/bobot/ocaml-flint.git"
+conflicts: [ "ocaml-option-bytecode-only" ]

arb/arb.ml

@@ -29,6 +29,14 @@ module ARF = struct
     let z = Flint.FMPZ.C.mk_fmpz () in
     ignore (arf_get_fmpz_fixed_si z x (Signed.Long.of_int n));
     Flint.FMPZ.to_z z
+
+  let of_fmpz_2exp ~exp m =
+    let arf = C.mk_arf () in
+    arf_set_fmpz_2exp arf m exp;
+    arf
+
+  let of_2exp ~exp m =
+    of_fmpz_2exp ~exp:(Flint.FMPZ.of_z exp) (Flint.FMPZ.of_z m)
 end
 
 module MAG = struct
@@ -79,12 +87,34 @@ module ARB = struct
   let pp fmt x = Format.pp_print_string fmt (External.to_string x)
   let mid x = x |-> ARB.mid
   let rad x = x |-> ARB.rad
+
+  let of_round_fmpz_2exp ?(prec = 0) ~exp base =
+    let arb = C.mk_arb () in
+    arb_set_round_fmpz_2exp arb base exp (Signed.Long.of_int prec);
+    arb
+
+  let of_round_2exp ?(prec = 0) ~exp base =
+    of_round_fmpz_2exp ~prec ~exp:(Flint.FMPZ.of_z exp) (Flint.FMPZ.of_z base)
+
+  let of_interval ?(prec = 0) min max =
+    let arb = C.mk_arb () in
+    arb_set_interval_arf arb min max (Signed.Long.of_int prec);
+    arb
+
+  let zero () =
+    let arb = C.mk_arb () in
+    arb_zero arb;
+    arb
 end
 
 module ACB = struct
   type t = ACB.t
 
   module C = struct
+    type acb = C.Type.ACB.a
+
+    let acb_struct = C.Type.ACB.s
+    let convert = Fun.id
     let acb_t = C.Type.ACB.t
 
     let mk_acb () : ACB.t =
@@ -106,4 +136,9 @@ module ACB = struct
   let rel_accuracy_bits t = Signed.Long.to_int @@ acb_rel_accuracy_bits t
   let real x = x |-> ACB.real
   let imag x = x |-> ACB.imag
+
+  let make ~real ~imag =
+    let acb = C.mk_acb () in
+    acb_set_arb_arb acb real imag;
+    acb
 end

arb/arb.mli

@@ -8,6 +8,8 @@ module ARF : sig
 
   val pp : Format.formatter -> t -> unit
   val get_fmpz_fixed_si : t -> int -> Z.t
+  val of_fmpz_2exp : exp:Flint.FMPZ.t -> Flint.FMPZ.t -> t
+  val of_2exp : exp:Z.t -> Z.t -> t
 end
 
 module MAG : sig
@@ -32,12 +34,20 @@ module ARB : sig
   val pp : Format.formatter -> t -> unit
   val mid : t -> ARF.t
   val rad : t -> MAG.t
+  val of_round_fmpz_2exp : ?prec:int -> exp:Flint.FMPZ.t -> Flint.FMPZ.t -> t
+  val of_round_2exp : ?prec:int -> exp:Z.t -> Z.t -> t
+  val of_interval : ?prec:int -> ARF.t -> ARF.t -> t
+  val zero : unit -> t
 end
 
 module ACB : sig
   type t
 
   module C : sig
+    type acb
+
+    val acb_struct : acb Ctypes.structure Ctypes.typ
+    val convert : acb Ctypes.structure Ctypes.ptr -> t
     val acb_t : t Ctypes.typ
     val mk_acb : unit -> t
   end
@@ -46,4 +56,5 @@ module ACB : sig
   val rel_accuracy_bits : t -> int
   val real : t -> ARB.t
   val imag : t -> ARB.t
+  val make : real:ARB.t -> imag:ARB.t -> t
 end

arb/function_description.ml

@@ -29,4 +29,22 @@ module Functions (F : Ctypes.FOREIGN) = struct
   let arf_get_fmpz_fixed_si =
     foreign "arf_get_fmpz_fixed_si"
       (Flint.FMPZ.C.fmpz_t @-> arf_t @-> long @-> returning bool)
+
+  let acb_set_arb_arb =
+    foreign "acb_set_arb_arb" (ACB.t @-> ARB.t @-> ARB.t @-> returning void)
+
+  let arb_set_round_fmpz_2exp =
+    foreign "arb_set_round_fmpz_2exp"
+      (ARB.t @-> Flint.FMPZ.C.fmpz_t @-> Flint.FMPZ.C.fmpz_t @-> long
+     @-> returning void)
+
+  let arf_set_fmpz_2exp =
+    foreign "arf_set_fmpz_2exp"
+      (ARF.t @-> Flint.FMPZ.C.fmpz_t @-> Flint.FMPZ.C.fmpz_t @-> returning void)
+
+  let arb_set_interval_arf =
+    foreign "arb_set_interval_arf"
+      (ARB.t @-> ARF.t @-> ARF.t @-> long @-> returning void)
+
+  let arb_zero = foreign "arb_zero" (ARB.t @-> returning void)
 end

calcium.opam

@@ -36,3 +36,4 @@ build: [
   ]
 ]
 dev-repo: "git+https://github.com/bobot/ocaml-flint.git"
+conflicts: [ "ocaml-option-bytecode-only" ]

calcium/calcium.ml

@@ -14,6 +14,62 @@ module CTX = struct
     ctx
 end
 
+module QQBAR = struct
+  type t = qqbar_t
+
+  module C = struct
+    let mk_qqbar () : t = allocate_n ~count:1 ~finalise:qqbar_clear qqbar_struct
+  end
+
+  let equal q1 q2 = qqbar_equal q1 q2
+  let compare q1 q2 = qqbar_cmp_root_order q1 q2
+  let hash q1 = Unsigned.ULong.to_int @@ qqbar_hash q1
+  let is_real = qqbar_is_real
+  let is_zero = qqbar_is_zero
+  let is_one = qqbar_is_one
+
+  let poly qqbar : Flint.FMPZ_poly.t =
+    Flint.FMPZ_poly.C.convert @@ (qqbar |-> QQBAR.poly)
+
+  let enclosure qqbar : Arb.ACB.t =
+    Arb.ACB.C.convert @@ (qqbar |-> QQBAR.enclosure)
+
+  let debug_print = qqbar_print
+
+  let from_enclosure p e =
+    let qqbar = C.mk_qqbar () in
+    qqbar_init qqbar;
+    if
+      qqbar_validate_existence_uniqueness (enclosure qqbar) p e
+        qqbar_default_prec
+    then (
+      Flint.FMPZ_poly.C.set ~dst:(poly qqbar) ~src:p;
+      Some qqbar)
+    else None
+
+  let from_roots ?(unsorted = false) ?(irreducible = false) p =
+    let deg = max 0 (Flint.FMPZ_poly.length p - 1) in
+    let finalise v =
+      for i = 0 to deg - 1 do
+        qqbar_clear (v +@ i)
+      done
+    in
+    let v = allocate_n ~count:deg ~finalise qqbar_struct in
+    for i = 0 to deg - 1 do
+      qqbar_init (v +@ i)
+    done;
+    let flags =
+      let flags = [] in
+      let flags =
+        if irreducible then QQBAR_ROOTS_IRREDUCIBLE :: flags else flags
+      in
+      let flags = if unsorted then QQBAR_ROOTS_UNSORTED :: flags else flags in
+      flags
+    in
+    qqbar_roots_fmpz_poly v p flags;
+    Array.init deg (fun i -> v +@ i)
+end
+
 module CA = struct
   type t = ca_t
 
@@ -134,4 +190,17 @@ module CA = struct
   let mod_e ~ctx a b = sub ~ctx a (mul ~ctx (of_z ~ctx (div_e ~ctx a b)) b)
   let mod_t ~ctx a b = sub ~ctx a (mul ~ctx (of_z ~ctx (div_t ~ctx a b)) b)
   let mod_f ~ctx a b = sub ~ctx a (mul ~ctx (of_z ~ctx (div_f ~ctx a b)) b)
+
+  let from_qqbar ~ctx qqbar =
+    let ca = mk_ca ~ctx () in
+    ca_set_qqbar ca qqbar ctx;
+    (* only algebraic currently *)
+    ca
+
+  let to_qqbar ~ctx ca =
+    let qqbar = QQBAR.C.mk_qqbar () in
+    let b = ca_get_qqbar qqbar ca ctx in
+    assert b;
+    (* only algebraic currently *)
+    qqbar
 end

calcium/calcium.mli

@@ -4,6 +4,27 @@ module CTX : sig
   val mk : unit -> t
 end
 
+module QQBAR : sig
+  type t
+
+  val equal : t -> t -> bool
+  val compare : t -> t -> int
+
+  val debug_print : t -> unit
+  (** to stdout *)
+
+  val is_real : t -> bool
+  val is_one : t -> bool
+  val is_zero : t -> bool
+  val poly : t -> Flint.FMPZ_poly.t
+  val enclosure : t -> Arb.ACB.t
+  val from_enclosure : Flint.FMPZ_poly.t -> Arb.ACB.t -> t option
+
+  val from_roots :
+    ?unsorted:bool -> ?irreducible:bool -> Flint.FMPZ_poly.t -> t array
+  (** default optional value is false *)
+end
+
 module CA : sig
   type t
 
@@ -73,4 +94,6 @@ module CA : sig
   val sqrt : ctx:CTX.t -> t -> t
   val pow_int : ctx:CTX.t -> t -> int -> t
   val pow : ctx:CTX.t -> t -> Q.t -> t
+  val from_qqbar : ctx:CTX.t -> QQBAR.t -> t
+  val to_qqbar : ctx:CTX.t -> t -> QQBAR.t
 end

calcium/function_description.ml

@@ -105,4 +105,53 @@ module Functions (F : Ctypes.FOREIGN) = struct
 
   let ca_check_is_negative_real =
     foreign "ca_check_is_negative_real" (ca_t @-> ca_ctx_t @-> returning truth_t)
+
+  let qqbar_init = foreign "qqbar_init" (qqbar_t @-> returning void)
+  let qqbar_clear = foreign "qqbar_clear" (qqbar_t @-> returning void)
+
+  let ca_get_qqbar =
+    foreign "ca_get_qqbar" (qqbar_t @-> ca_t @-> ca_ctx_t @-> returning bool)
+
+  let ca_set_qqbar =
+    foreign "ca_set_qqbar" (ca_t @-> qqbar_t @-> ca_ctx_t @-> returning void)
+
+  let flag_qqbar_roots =
+    let cons_if i cst cst' l = if Int.logand i cst = 0 then l else cst' :: l in
+    view int
+      ~read:(fun i ->
+        cons_if i qqbar_roots_irreducible QQBAR_ROOTS_IRREDUCIBLE
+        @@ cons_if i qqbar_roots_unsorted QQBAR_ROOTS_UNSORTED
+        @@ [])
+      ~write:(fun l ->
+        List.fold_left
+          (fun acc -> function
+            | QQBAR_ROOTS_IRREDUCIBLE -> Int.logor acc qqbar_roots_irreducible
+            | QQBAR_ROOTS_UNSORTED -> Int.logor acc qqbar_roots_unsorted)
+          0 l)
+
+  let qqbar_roots_fmpz_poly =
+    foreign "qqbar_roots_fmpz_poly"
+      (qqbar_t @-> Flint.FMPZ_poly.C.fmpz_poly_t @-> flag_qqbar_roots
+     @-> returning void)
+
+  let qqbar_validate_existence_uniqueness =
+    foreign "_qqbar_validate_existence_uniqueness"
+      (Arb.ACB.C.acb_t @-> Flint.FMPZ_poly.C.fmpz_poly_t @-> Arb.ACB.C.acb_t
+     @-> long @-> returning bool)
+
+  let qqbar_set = foreign "qqbar_set" (qqbar_t @-> qqbar_t @-> returning void)
+  let qqbar_swap = foreign "qqbar_swap" (qqbar_t @-> qqbar_t @-> returning void)
+
+  let qqbar_equal =
+    foreign "qqbar_equal" (qqbar_t @-> qqbar_t @-> returning bool)
+
+  let qqbar_hash = foreign "qqbar_hash" (qqbar_t @-> returning ulong)
+
+  let qqbar_cmp_root_order =
+    foreign "qqbar_cmp_root_order" (qqbar_t @-> qqbar_t @-> returning int)
+
+  let qqbar_is_real = foreign "qqbar_is_real" (qqbar_t @-> returning bool)
+  let qqbar_is_zero = foreign "qqbar_is_zero" (qqbar_t @-> returning bool)
+  let qqbar_is_one = foreign "qqbar_is_one" (qqbar_t @-> returning bool)
+  let qqbar_print = foreign "qqbar_print" (qqbar_t @-> returning void)
 end

calcium/tests/units.expected

@@ -8,3 +8,8 @@ acb8:61 -> 92681
 acb0:28 -> 92681
 acb16:53 -> 92681
 h1=h2:true
+x^2-2
+r0:1.41421 {a where a = 1.41421 [a^2-2=0]}
+r1:-1.41421 {-a where a = 1.41421 [a^2-2=0]}
+acb:((5 * 2^-2) +/- (536870913 * 2^-31), (0) +/- (0))
+a:1.41421 {a where a = 1.41421 [a^2-2=0]}

calcium/tests/units.ml

@@ -29,3 +29,25 @@ let () =
   let z2' = Calcium.CA.pow ~ctx z1 (Q.of_string "0.5") in
   let h2 = Calcium.CA.hash ~ctx z2' in
   Format.printf "h1=h2:%b@." (h1 = h2)
+
+let () =
+  let p = Flint.FMPZ_poly.create [| Z.of_int (-2); Z.of_int 0; Z.of_int 1 |] in
+  Format.printf "%a@." Flint.FMPZ_poly.pp p;
+  let roots = Calcium.QQBAR.from_roots p in
+  Array.iteri
+    (fun i a ->
+      Format.printf "r%i:%a@." i (Calcium.CA.pp ~ctx)
+        (Calcium.CA.from_qqbar ~ctx a))
+    roots
+
+let () =
+  let min = Arb.ARF.of_2exp (Z.of_int 1) ~exp:Z.zero in
+  let max = Arb.ARF.of_2exp (Z.of_int 3) ~exp:Z.minus_one in
+  let arb = Arb.ARB.of_interval min max in
+  let acb = Arb.ACB.make ~real:arb ~imag:(Arb.ARB.zero ()) in
+  Format.printf "acb:%a@." Arb.ACB.pp acb;
+  let p = Flint.FMPZ_poly.create [| Z.of_int (-2); Z.of_int 0; Z.of_int 1 |] in
+  let a = Calcium.QQBAR.from_enclosure p acb in
+  match a with
+  | None -> Format.printf "no roots@."
+  | Some a -> pp "a" (Calcium.CA.from_qqbar ~ctx a)

calcium/type_description.ml

@@ -34,4 +34,29 @@ module Types (F : Ctypes.TYPE) = struct
     let t_unknown = constant "T_UNKNOWN" int64_t in
     enum ~typedef:true "truth_t"
       [ (TRUE, t_true); (FALSE, t_false); (UNKNOWN, t_unknown) ]
+
+  type qqbar
+
+  let qqbar_struct : qqbar structure typ =
+    typedef (structure "qqbar_struct_struct") "qqbar_struct"
+
+  module QQBAR = struct
+    let poly =
+      field qqbar_struct "poly" (F.lift_typ Flint.FMPZ_poly.C.fmpz_poly_struct)
+
+    let enclosure =
+      field qqbar_struct "enclosure" (F.lift_typ Arb.ACB.C.acb_struct)
+
+    let () = seal qqbar_struct
+  end
+
+  type qqbar_t = qqbar structure ptr
+
+  let qqbar_t : qqbar_t typ = ptr qqbar_struct
+
+  type flag_qqbar_roots = QQBAR_ROOTS_IRREDUCIBLE | QQBAR_ROOTS_UNSORTED
+
+  let qqbar_roots_irreducible = constant "QQBAR_ROOTS_IRREDUCIBLE" int
+  let qqbar_roots_unsorted = constant "QQBAR_ROOTS_UNSORTED" int
+  let qqbar_default_prec = constant "QQBAR_DEFAULT_PREC" long
 end

dune

@@ -0,0 +1,4 @@
+(rule (copy opam.template flint.opam.template))
+(rule (copy opam.template antic.opam.template))
+(rule (copy opam.template arb.opam.template))
+(rule (copy opam.template calcium.opam.template))

flint.opam

@@ -32,3 +32,4 @@ build: [
   ]
 ]
 dev-repo: "git+https://github.com/bobot/ocaml-flint.git"
+conflicts: [ "ocaml-option-bytecode-only" ]