From d2e7848de2cc4b94a8a0edcc3d605ab9ca3aaf1a Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 28 Dec 2020 01:08:55 +0100 Subject: [PATCH] Finished common in org-mode --- common/README.org | 4 +- common/lib/dune | 2 +- common/lib/math_functions.c | 85 --- common/lib/powers.ml | 57 +- common/lib/powers.mli | 66 +- common/lib/qcaml.ml | 12 +- common/lib/qcaml.mli | 14 +- common/lib/range.ml | 23 +- common/lib/range.mli | 34 +- common/lib/spin.ml | 67 +- common/lib/spin.mli | 19 +- common/lib/util.c | 97 +++ common/lib/util.ml | 435 +++++++------ common/lib/util.mli | 194 ++---- common/lib/zkey.ml | 236 ++++--- common/lib/zkey.mli | 84 +-- common/lib/zmap.ml | 5 +- common/lib/zmap.mli | 7 +- common/powers.org | 128 ++++ common/qcaml.org | 44 ++ common/range.org | 96 +++ common/spin.org | 63 ++ common/test/{math_functions.ml => util.ml} | 50 +- common/util.org | 683 +++++++++++++++++++++ common/zkey.org | 346 +++++++++++ common/zmap.org | 29 + docs/config.el | 1 + test/run_tests.ml | 2 +- 28 files changed, 2155 insertions(+), 728 deletions(-) delete mode 100644 common/lib/math_functions.c create mode 100644 common/lib/util.c create mode 100644 common/powers.org create mode 100644 common/qcaml.org create mode 100644 common/range.org create mode 100644 common/spin.org rename common/test/{math_functions.ml => util.ml} (69%) create mode 100644 common/util.org create mode 100644 common/zkey.org create mode 100644 common/zmap.org diff --git a/common/README.org b/common/README.org index 4a103a8..bf6559d 100644 --- a/common/README.org +++ b/common/README.org @@ -52,12 +52,12 @@ This directory contains many utility functions used by all the other directories *** Extra C files - The ~math_functions~ file contains small C snippets to add missing + The ~util.c~ file contains small C snippets to add missing functionalities to OCaml, such as support for the ~popcnt~ instruction. #+begin_src elisp :tangle (org-entry-get nil "dune" t) (c_names - math_functions + util ) (c_flags (:standard) -Ofast -march=native -fPIC diff --git a/common/lib/dune b/common/lib/dune index 66db8a1..3068c2d 100644 --- a/common/lib/dune +++ b/common/lib/dune @@ -11,7 +11,7 @@ ) (c_names - math_functions + util ) (c_flags (:standard) -Ofast -march=native -fPIC diff --git a/common/lib/math_functions.c b/common/lib/math_functions.c deleted file mode 100644 index 353141f..0000000 --- a/common/lib/math_functions.c +++ /dev/null @@ -1,85 +0,0 @@ -#include -#include -#include - - -CAMLprim value erf_float_bytecode(value x) -{ - return copy_double(erf(Double_val(x))); -} - -CAMLprim double erf_float(double x) -{ - return erf(x); -} - - -CAMLprim value erfc_float_bytecode(value x) -{ - return copy_double(erfc(Double_val(x))); -} - -CAMLprim double erfc_float(double x) -{ - return erfc(x); -} - - -CAMLprim value gamma_float_bytecode(value x) -{ - return copy_double(tgamma(Double_val(x))); -} - - -CAMLprim double gamma_float(double x) -{ - return tgamma(x); -} - - -#include -CAMLprim int32_t popcnt(int64_t i) -{ - return __builtin_popcountll (i); -} - - -CAMLprim value popcnt_bytecode(value i) -{ - return caml_copy_int32(__builtin_popcountll (Int64_val(i))); -} - - -CAMLprim int32_t trailz(int64_t i) -{ - return __builtin_ctzll (i); -} - - -CAMLprim value trailz_bytecode(value i) -{ - return caml_copy_int32(__builtin_ctzll (Int64_val(i))); -} - -CAMLprim int32_t leadz(int64_t i) -{ - return __builtin_clzll(i); -} - - -CAMLprim value leadz_bytecode(value i) -{ - return caml_copy_int32(__builtin_clzll (Int64_val(i))); -} - - - -#include - -CAMLprim value unix_vfork(value unit) -{ - int ret; - ret = vfork(); - return Val_int(ret); -} - diff --git a/common/lib/powers.ml b/common/lib/powers.ml index 62fc1ee..6e34578 100644 --- a/common/lib/powers.ml +++ b/common/lib/powers.ml @@ -1,5 +1,29 @@ -type t = { x: int ; y : int ; z : int ; tot : int } + +(* ~tot~ always contains ~x+y+z~. *) + + +(* [[file:../powers.org::*Type][Type:2]] *) +type t = { + x : int ; + y : int ; + z : int ; + tot : int ; +} +(* Type:2 ends here *) + + + +(* #+begin_example + * Powers.of_int_tuple (2,3,1);; + * - : Powers.t = {Qcaml.Common.Powers.x = 2; y = 3; z = 1; tot = 6} + * + * Powers.(to_int_tuple (of_int_tuple (2,3,1)));; + * - : int * int * int = (2, 3, 1) + * #+end_example *) + + +(* [[file:../powers.org::*Conversions][Conversions:2]] *) let of_int_tuple t = let result = match t with @@ -12,8 +36,30 @@ let of_int_tuple t = invalid_arg (__FILE__^": of_int_tuple"); result -let to_int_tuple { x ; y ; z ; _ } = (x,y,z) +let to_int_tuple { x ; y ; z ; _ } = (x,y,z) +(* Conversions:2 ends here *) + + + +(* | ~get~ | Returns the value of the power for $x$, $y$ or $z$ + * | ~incr~ | Returns a new ~Powers.t~ with the power on the given axis incremented | + * | ~decr~ | Returns a new ~Powers.t~ with the power on the given axis decremented. As opposed to ~of_int_tuple~, the values may become negative| + * + * #+begin_example + * Powers.get Coordinate.Y (Powers.of_int_tuple (2,3,1));; + * - : int = 3 + * + * Powers.incr Coordinate.Y (Powers.of_int_tuple (2,3,1));; + * - : Powers.t = {Qcaml.Common.Powers.x = 2; y = 4; z = 1; tot = 7} + * + * Powers.decr Coordinate.Y (Powers.of_int_tuple (2,3,1));; + * - : Powers.t = {Qcaml.Common.Powers.x = 2; y = 2; z = 1; tot = 5} + * #+end_example *) + + + +(* [[file:../powers.org::*Operations][Operations:2]] *) let get coord t = match coord with | Coordinate.X -> t.x @@ -31,6 +77,9 @@ let decr coord t = | Coordinate.X -> let r = t.x-1 in { t with x = r ; tot = t.tot-1 } | Coordinate.Y -> let r = t.y-1 in { t with y = r ; tot = t.tot-1 } | Coordinate.Z -> let r = t.z-1 in { t with z = r ; tot = t.tot-1 } +(* Operations:2 ends here *) - - +(* [[file:../powers.org::*Printers][Printers:2]] *) +let pp ppf t = + Format.fprintf ppf "@[x^%d + y^%d + z^%d@]" t.x t.y t.z +(* Printers:2 ends here *) diff --git a/common/lib/powers.mli b/common/lib/powers.mli index bcc9732..69cdcdf 100644 --- a/common/lib/powers.mli +++ b/common/lib/powers.mli @@ -1,57 +1,35 @@ -(** Contains powers of x, y and z describing the polynomials in atomic basis sets. *) +(* Type *) + +(* [[file:../powers.org::*Type][Type:1]] *) type t = private { - x : int ; - y : int ; - z : int ; - tot : int ; (* x + y + z *) - } + x : int ; + y : int ; + z : int ; + tot : int ; +} +(* Type:1 ends here *) + +(* Conversions *) +(* [[file:../powers.org::*Conversions][Conversions:1]] *) val of_int_tuple : int * int * int -> t -(** Example: - [of_int_tuple (2,3,1) -> { x=2 ; y=3 ; z=1 ; tot=6 }] - @raise Invalid_argument if x, y or z < 0. - *) - - -val to_int_tuple : t -> int * int * int -(** Example: - [to_int_tuple { x=2 ; y=3 ; z=1 ; tot=6 } -> (2,3,1) ] - *) +val to_int_tuple : t -> int * int * int +(* Conversions:1 ends here *) + +(* Operations *) +(* [[file:../powers.org::*Operations][Operations:1]] *) val get : Coordinate.axis -> t -> int -(** Example: - - [Powers.get Coordinate.Y { x=2 ; y=3 ; z=1 ; tot=6 } -> 3] - - *) - - val incr : Coordinate.axis -> t -> t -(** Returns a new {!Powers.t} with the power on the given axis incremented. - - Example: - - {[ - Powers.incr Coordinate.Y { x=2 ; y=3 ; z=1 ; tot=6 } -> - { x=2 ; y=4 ; z=1 ; tot=7 } - ]} - - *) - - val decr : Coordinate.axis -> t -> t -(** Returns a new {!Powers.t} with the power on the given axis decremented. - As opposed to {!of_int_tuple}, the values may become negative. +(* Operations:1 ends here *) - Example: +(* Printers *) - {[ - Powers.incr Coordinate.Y { x=2 ; y=3 ; z=1 ; tot=6 } -> - { x=2 ; y=2 ; z=1 ; tot=5 } - ]} - - *) +(* [[file:../powers.org::*Printers][Printers:1]] *) +val pp : Format.formatter -> t -> unit +(* Printers:1 ends here *) diff --git a/common/lib/qcaml.ml b/common/lib/qcaml.ml index 38aa96b..66cc4e7 100644 --- a/common/lib/qcaml.ml +++ b/common/lib/qcaml.ml @@ -1,6 +1,12 @@ -let name = "QCaml" +(* | ~root~ | Path to the QCaml source directory | + * | ~name~ | ~"QCaml"~ | *) + + +(* [[file:../qcaml.org::*QCaml][QCaml:2]] *) +let name = "QCaml" + let root = let rec chop = function | [] -> [] @@ -12,6 +18,4 @@ let root = |> chop |> List.rev |> String.concat Filename.dir_sep - - - +(* QCaml:2 ends here *) diff --git a/common/lib/qcaml.mli b/common/lib/qcaml.mli index 025fa67..fc21df1 100644 --- a/common/lib/qcaml.mli +++ b/common/lib/qcaml.mli @@ -1,6 +1,12 @@ +(* QCaml + * :PROPERTIES: + * :header-args: :noweb yes :comments both + * :END: + * + * QCaml-specific parameters *) + -val name : string -(** Name of the QCaml source directory on the file system. *) - +(* [[file:../qcaml.org::*QCaml][QCaml:1]] *) val root : string -(** Path to the QCaml source directory. *) +val name : string +(* QCaml:1 ends here *) diff --git a/common/lib/range.ml b/common/lib/range.ml index f424717..91a1671 100644 --- a/common/lib/range.ml +++ b/common/lib/range.ml @@ -1,5 +1,8 @@ +(* [[file:../range.org::*Type][Type:2]] *) type t = int list +(* Type:2 ends here *) +(* [[file:../range.org::*Conversion][Conversion:2]] *) let to_int_list r = r let expand_range r = @@ -24,13 +27,13 @@ let of_string s = match s.[0] with | '0' .. '9' -> [ int_of_string s ] | _ -> - assert (s.[0] = '[') ; - assert (s.[(String.length s)-1] = ']') ; - let s = String.sub s 1 ((String.length s) - 2) in - let l = String.split_on_char ',' s in - let l = List.map expand_range l in - List.concat l - |> List.sort_uniq compare + assert (s.[0] = '[') ; + assert (s.[(String.length s)-1] = ']') ; + let s = String.sub s 1 ((String.length s) - 2) in + let l = String.split_on_char ',' s in + let l = List.map expand_range l in + List.concat l + |> List.sort_uniq compare let to_string l = @@ -38,9 +41,9 @@ let to_string l = (List.map string_of_int l |> String.concat ",") ^ "]" +(* Conversion:2 ends here *) - +(* [[file:../range.org::*Printers][Printers:2]] *) let pp ppf t = Format.fprintf ppf "@[%s@]" (to_string t) - - +(* Printers:2 ends here *) diff --git a/common/lib/range.mli b/common/lib/range.mli index 8a18174..e216da8 100644 --- a/common/lib/range.mli +++ b/common/lib/range.mli @@ -1,30 +1,22 @@ -(** A range is a sorted list of integers in an interval. - - {[ "[36-53,72-107,126-131]" ]} - represents the list of integers - {[ [ 37 ; 37 ; 38 ; ... ; 52 ; 53 ; 72 ; 73 ; ... ; 106 ; 107 ; 126 ; 127 ; ... - ; 130 ; 131 ] ]} -*) +(* Type *) + +(* [[file:../range.org::*Type][Type:1]] *) type t +(* Type:1 ends here *) -val of_string : string -> t -(** Create from a string: - - - "[a-b]" : range between a and b (included) - - "[a]" : the list with only one integer a - - "a" : equivalent to "[a]" -*) - - -val to_string : t -> string -(** String representation. *) +(* Conversion *) +(* [[file:../range.org::*Conversion][Conversion:1]] *) +val of_string : string -> t +val to_string : t -> string val to_int_list : t -> int list -(** Transform into a list of ints. *) +(* Conversion:1 ends here *) -(** {2 Printers} *) +(* Printers *) + +(* [[file:../range.org::*Printers][Printers:1]] *) val pp : Format.formatter -> t -> unit - +(* Printers:1 ends here *) diff --git a/common/lib/spin.ml b/common/lib/spin.ml index 63ffa8b..cfc8062 100644 --- a/common/lib/spin.ml +++ b/common/lib/spin.ml @@ -1,47 +1,32 @@ -(** Electron spin *) + +(* Note : + * ~Alfa~ if written with an 'f' instead of 'ph' because it has the same number of + * letters as ~Beta~, so the alignment of the code is nicer. *) + + +(* [[file:../spin.org::*Type][Type:2]] *) type t = (* m_s *) -| Alfa (* {% $m_s = +1/2$ %} *) -| Beta (* {% $m_s = -1/2$ %} *) + | Alfa (* {% $m_s = +1/2$ %} *) + | Beta (* {% $m_s = -1/2$ %} *) +(* Type:2 ends here *) + + +(* Returns the opposite spin *) + + +(* [[file:../spin.org::*Functions][Functions:2]] *) let other = function -| Alfa -> Beta -| Beta -> Alfa + | Alfa -> Beta + | Beta -> Alfa -(* -let half = 1. /. 2. +let to_string = function + | Alfa -> "Alpha" + | Beta -> "Beta " +(* Functions:2 ends here *) - -(** {% $\alpha(m_s)$ %} *) -let alfa = function -| n, Alfa -> n -| _, Beta -> 0. - - -(** {% $\beta(m_s)$ %} *) -let beta = function -| _, Alfa -> 0. -| n, Beta -> n - - -(** {% $S_z(m_s)$ %} *) -let s_z = function -| n, Alfa -> half *. n, Alfa -| n, Beta -> -. half *. n, Beta - -let s_plus = function -| n, Beta -> n , Alfa -| _, Alfa -> 0., Alfa - -let s_minus = function -| n, Alfa -> n , Beta -| _, Beta -> 0., Beta - -let ( ++ ) (n, t) (m,t) = - (m. +n.) - -let s2 s = - s_minus @@ s_plus s +. - s_z s +. - s_z @@ s_z s - *) +(* [[file:../spin.org::*Printers][Printers:2]] *) +let pp ppf t = + Format.fprintf ppf "@[%s@]" (to_string t) +(* Printers:2 ends here *) diff --git a/common/lib/spin.mli b/common/lib/spin.mli index 89e3a82..27b2720 100644 --- a/common/lib/spin.mli +++ b/common/lib/spin.mli @@ -1,13 +1,20 @@ -(** Electron spin. +(* Type *) -Note : -[Alfa] if written with an 'f' instead of 'ph' because it has the same number of -letters as [Beta], so the alignment of the code is nicer. - -*) +(* [[file:../spin.org::*Type][Type:1]] *) type t = Alfa | Beta +(* Type:1 ends here *) +(* Functions *) + + +(* [[file:../spin.org::*Functions][Functions:1]] *) val other : t -> t +(* Functions:1 ends here *) + +(* Printers *) +(* [[file:../spin.org::*Printers][Printers:1]] *) +val pp : Format.formatter -> t -> unit +(* Printers:1 ends here *) diff --git a/common/lib/util.c b/common/lib/util.c new file mode 100644 index 0000000..87cf953 --- /dev/null +++ b/common/lib/util.c @@ -0,0 +1,97 @@ +/* External C functions */ + +/* | ~erf_float~ | Error function ~erf~ from =libm= | */ +/* | ~erfc_float~ | Complementary error function ~erfc~ from =libm= | */ +/* | ~gamma_float~ | Gamma function ~gamma~ from =libm= | */ +/* | ~popcnt~ | ~popcnt~ instruction | */ +/* | ~trailz~ | ~ctz~ instruction | */ +/* | ~leadz~ | ~bsf~ instruction | */ + + +/* [[file:../util.org::*External C functions][External C functions:1]] */ +#include +#include +#include +/* External C functions:1 ends here */ + +/* Erf */ + + +/* [[file:../util.org::*Erf][Erf:1]] */ +CAMLprim value erf_float_bytecode(value x) { + return copy_double(erf(Double_val(x))); +} + +CAMLprim double erf_float(double x) { + return erf(x); +} +/* Erf:1 ends here */ + +/* Erfc */ + + +/* [[file:../util.org::*Erfc][Erfc:1]] */ +CAMLprim value erfc_float_bytecode(value x) { + return copy_double(erfc(Double_val(x))); +} + +CAMLprim double erfc_float(double x) { + return erfc(x); +} +/* Erfc:1 ends here */ + +/* Gamma */ + + +/* [[file:../util.org::*Gamma][Gamma:1]] */ +CAMLprim value gamma_float_bytecode(value x) { + return copy_double(tgamma(Double_val(x))); +} + + +CAMLprim double gamma_float(double x) { + return tgamma(x); +} +/* Gamma:1 ends here */ + +/* Popcnt */ + + +/* [[file:../util.org::*Popcnt][Popcnt:1]] */ +CAMLprim int32_t popcnt(int64_t i) { + return __builtin_popcountll (i); +} + + +CAMLprim value popcnt_bytecode(value i) { + return caml_copy_int32(__builtin_popcountll (Int64_val(i))); +} +/* Popcnt:1 ends here */ + +/* Trailz */ + + +/* [[file:../util.org::*Trailz][Trailz:1]] */ +CAMLprim int32_t trailz(int64_t i) { + return __builtin_ctzll (i); +} + + +CAMLprim value trailz_bytecode(value i) { + return caml_copy_int32(__builtin_ctzll (Int64_val(i))); +} +/* Trailz:1 ends here */ + +/* Leadz */ + + +/* [[file:../util.org::*Leadz][Leadz:1]] */ +CAMLprim int32_t leadz(int64_t i) { + return __builtin_clzll(i); +} + + +CAMLprim value leadz_bytecode(value i) { + return caml_copy_int32(__builtin_clzll (Int64_val(i))); +} +/* Leadz:1 ends here */ diff --git a/common/lib/util.ml b/common/lib/util.ml index 103501b..f02c738 100644 --- a/common/lib/util.ml +++ b/common/lib/util.ml @@ -1,72 +1,155 @@ -(** All utilities which should be included in all source files are defined here *) +(* [[file:../util.org::*Erf][Erf:3]] *) +external erf_float : float -> float = "erf_float_bytecode" "erf_float" [@@unboxed] [@@noalloc] +(* Erf:3 ends here *) -(** {1 Functions from libm} *) +(* [[file:../util.org::*Erfc][Erfc:3]] *) +external erfc_float : float -> float = "erfc_float_bytecode" "erfc_float" [@@unboxed] [@@noalloc] +(* Erfc:3 ends here *) -open Constants - - - -external erf_float : float -> float = "erf_float_bytecode" "erf_float" - [@@unboxed] [@@noalloc] - -external erfc_float : float -> float = "erfc_float_bytecode" "erfc_float" - [@@unboxed] [@@noalloc] - -external gamma_float : float -> float = "gamma_float_bytecode" "gamma_float" - [@@unboxed] [@@noalloc] +(* [[file:../util.org::*Gamma][Gamma:3]] *) +external gamma_float : float -> float = "gamma_float_bytecode" "gamma_float" [@@unboxed] [@@noalloc] +(* Gamma:3 ends here *) +(* [[file:../util.org::*Popcnt][Popcnt:3]] *) external popcnt : int64 -> int32 = "popcnt_bytecode" "popcnt" - [@@unboxed] [@@noalloc] - (** popcnt instruction *) +[@@unboxed] [@@noalloc] -let popcnt i = (popcnt [@inlined] ) i |> Int32.to_int +let popcnt i = (popcnt [@inlined] ) i |> Int32.to_int +(* Popcnt:3 ends here *) +(* [[file:../util.org::*Trailz][Trailz:3]] *) external trailz : int64 -> int32 = "trailz_bytecode" "trailz" "int" - [@@unboxed] [@@noalloc] - (** ctz instruction *) +[@@unboxed] [@@noalloc] -let trailz i = trailz i |> Int32.to_int +let trailz i = trailz i |> Int32.to_int +(* Trailz:3 ends here *) +(* [[file:../util.org::*Leadz][Leadz:3]] *) external leadz : int64 -> int32 = "leadz_bytecode" "leadz" "int" - [@@unboxed] [@@noalloc] - (** bsf instruction *) +[@@unboxed] [@@noalloc] -external vfork : unit -> int = "unix_vfork" "unix_vfork" - -let leadz i = leadz i |> Int32.to_int +let leadz i = leadz i |> Int32.to_int +(* Leadz:3 ends here *) -exception SIGTERM -let () = - let f _ = raise SIGTERM in - Sys.set_signal Sys.sigint (Sys.Signal_handle f) - ;; +(* | ~fact~ | Factorial function. | + * | ~binom~ | Binomial coefficient. ~binom n k~ = $C_n^k$ | + * | ~binom_float~ | float variant of ~binom~ | + * | ~pow~ | Fast implementation of the power function for small integer powers | + * | ~chop~ | In ~chop a f~, evaluate ~f~ only if the absolute value of ~a~ is larger than ~Constants.epsilon~, and return ~a *. f ()~. | + * | ~float_of_int_fast~ | Faster implementation of float_of_int for small positive ints | + * | ~not_implemented~ | Fails with error message if some functionality is not implemented | + * | ~of_some~ | Extracts the value of an option | *) -let not_implemented () = - failwith "Not implemented" +(* [[file:../util.org::*General functions][General functions:2]] *) let memo_float_of_int = Array.init 64 float_of_int let float_of_int_fast i = - if Int.logand i 63 = i then - memo_float_of_int.(i) - else - float_of_int i + if Int.logand i 63 = i then + memo_float_of_int.(i) + else + float_of_int i let factmax = 150 -(* Incomplete gamma function : - {% $\gamma(\alpha,x) = \int_0^x e^{-t} t^{\alpha-1} dt$ %} - {% $p: \frac{1}{\Gamma(\alpha)} \int_0^x e^{-t} t^{\alpha-1} dt$ %} - {% $q: \frac{1}{\Gamma(\alpha)} \int_x^\infty e^{-t} t^{\alpha-1} dt$ %} +let fact_memo = + let rec aux accu_l accu = function + | 0 -> (aux [@tailcall]) [1.] 1. 1 + | i when (i = factmax) -> + let x = (float_of_int factmax) *. accu in + List.rev (x::accu_l) + | i -> let x = (float_of_int i) *. accu in + (aux [@tailcall]) (x::accu_l) x (i+1) + in + aux [] 0. 0 + |> Array.of_list - reference - Haruhiko Okumura: C-gengo niyoru saishin algorithm jiten - (New Algorithm handbook in C language) (Gijyutsu hyouron - sha, Tokyo, 1991) p.227 [in Japanese] *) +let fact = function + | i when (i < 0) -> + raise (Invalid_argument "Argument of factorial should be non-negative") + | i when (i > 150) -> + raise (Invalid_argument "Result of factorial is infinite") + | i -> fact_memo.(i) + +let binom = + let memo = + let m = Array.make_matrix 64 64 0 in + for n=0 to Array.length m - 1 do + m.(n).(0) <- 1; + m.(n).(n) <- 1; + for k=1 to (n - 1) do + m.(n).(k) <- m.(n-1).(k-1) + m.(n-1).(k) + done + done; + m + in + let rec f n k = + assert (k >= 0); + assert (n >= k); + if k = 0 || k = n then + 1 + else if n < 64 then + memo.(n).(k) + else + f (n-1) (k-1) + f (n-1) k + in f + + +let binom_float n k = + binom n k + |> float_of_int_fast + + +let rec pow a = function + | 0 -> 1. + | 1 -> a + | 2 -> a *. a + | 3 -> a *. a *. a + | -1 -> 1. /. a + | n when n > 0 -> + let b = pow a (n / 2) in + b *. b *. (if n mod 2 = 0 then 1. else a) + | n when n < 0 -> (pow [@tailcall]) (1./.a) (-n) + | _ -> assert false + + +let chop f g = + if (abs_float f) < Constants.epsilon then 0. + else f *. (g ()) + + +let not_implemented () = + failwith "Not implemented" + + +let of_some = function + | Some a -> a + | None -> assert false +(* General functions:2 ends here *) + + + +(* The lower [[https://en.wikipedia.org/wiki/Incomplete_gamma_function][Incomplete Gamma function]] is implemented : + * \[ + * \gamma(\alpha,x) = \int_0^x e^{-t} t^{\alpha-1} dt + * \] + * + * p: $\frac{1}{\Gamma(\alpha)} \int_0^x e^{-t} t^{\alpha-1} dt$ + * + * q: $\frac{1}{\Gamma(\alpha)} \int_x^\infty e^{-t} t^{\alpha-1} dt$ + * + * Reference : Haruhiko Okumura: C-gengo niyoru saishin algorithm jiten + * (New Algorithm handbook in C language) (Gijyutsu hyouron sha, + * Tokyo, 1991) p.227 [in Japanese] *) + + + +(* [[file:../util.org::*Functions related to the Boys function][Functions related to the Boys function:2]] *) let incomplete_gamma ~alpha x = assert (alpha >= 0.); assert (x >= 0.); @@ -108,124 +191,54 @@ let incomplete_gamma ~alpha x = let lb = (1. +. x -. a) in qg_loop min_float (w /. lb) 1. lb w 2.0 in - gf *. p_gamma x + gf *. p_gamma x +(* Functions related to the Boys function:2 ends here *) +(* The [[https://link.springer.com/article/10.1007/s10910-005-9023-3][Generalized Boys function]] is implemented, + * ~maxm~ is the maximum total angular momentum. + * + * \[ + * F_m(x) = \frac{\gamma(m+1/2,x)}{2x^{m+1/2}} + * \] + * where $\gamma$ is the incomplete gamma function. + * + * - $F_0(0.) = 1$ + * - $F_0(t) = \frac{\sqrt{\pi}}{2\sqrt{t}} \text{erf} ( \sqrt{t} )$ + * - $F_m(0.) = \frac{1}{2m+1}$ + * - $F_m(t) = \frac{\gamma{m+1/2,t}}{2t^{m+1/2}}$ + * - $F_m(t) = \frac{ 2t\, F_{m+1}(t) + e^{-t} }{2m+1}$ *) -let fact_memo = - let rec aux accu_l accu = function - | 0 -> (aux [@tailcall]) [1.] 1. 1 - | i when (i = factmax) -> - let x = (float_of_int factmax) *. accu in - List.rev (x::accu_l) - | i -> let x = (float_of_int i) *. accu in - (aux [@tailcall]) (x::accu_l) x (i+1) - in - aux [] 0. 0 - |> Array.of_list - - - -let fact = function - | i when (i < 0) -> - raise (Invalid_argument "Argument of factorial should be non-negative") - | i when (i > 150) -> - raise (Invalid_argument "Result of factorial is infinite") - | i -> fact_memo.(i) - - -let binom = - let memo = - let m = - Array.make_matrix 64 64 0 - in - for n=0 to Array.length m - 1 do - m.(n).(0) <- 1; - m.(n).(n) <- 1; - for k=1 to (n - 1) do - m.(n).(k) <- m.(n-1).(k-1) + m.(n-1).(k) - done - done; - m - in - let rec f n k = - assert (k >= 0); - assert (n >= k); - if k = 0 || k = n then - 1 - else if n < 64 then - memo.(n).(k) - else - f (n-1) (k-1) + f (n-1) k - in f - - -let binom_float n k = - binom n k - |> float_of_int_fast - - -let rec pow a = function - | 0 -> 1. - | 1 -> a - | 2 -> a *. a - | 3 -> a *. a *. a - | -1 -> 1. /. a - | n when n > 0 -> - let b = pow a (n / 2) in - b *. b *. (if n mod 2 = 0 then 1. else a) - | n when n < 0 -> (pow [@tailcall]) (1./.a) (-n) - | _ -> assert false - - - - -let chop f g = - if (abs_float f) < Constants.epsilon then 0. - else f *. (g ()) - - - -(** Generalized Boys function. - maxm : Maximum total angular momentum - {% $F_m(x) = \frac{\gamma(m+1/2,x)}{2x^{m+1/2}}$ %} - where %{ $\gamma$ %} is the incomplete gamma function. - - {% $F_0(0.) = 1$ %} - {% $F_0(t) = \frac{\sqrt{\pi}}{2\sqrt{t}} \text{erf} ( \sqrt{t} )$ %} - {% $F_m(0.) = \frac{1}{2m+1}$ %} - {% $F_m(t) = \frac{\gamma{m+1/2,t}}{2t^{m+1/2}} - {% $F_m(t) = \frac{ 2t\, F_{m+1}(t) + e^{-t} }{2m+1}$ %} -*) +(* [[file:../util.org::*Functions related to the Boys function][Functions related to the Boys function:4]] *) let boys_function ~maxm t = assert (t >= 0.); match maxm with | 0 -> - begin - if t = 0. then [| 1. |] else - let sq_t = sqrt t in - [| (sq_pi_over_two /. sq_t) *. erf_float sq_t |] - end + begin + if t = 0. then [| 1. |] else + let sq_t = sqrt t in + [| (Constants.sq_pi_over_two /. sq_t) *. erf_float sq_t |] + end | _ -> - begin - assert (maxm > 0); - let result = - Array.init (maxm+1) (fun m -> 1. /. float_of_int (2*m+1)) - in - let power_t_inv = (maxm+maxm+1) in - try + begin + assert (maxm > 0); + let result = + Array.init (maxm+1) (fun m -> 1. /. float_of_int (2*m+1)) + in + let power_t_inv = (maxm+maxm+1) in + try let fmax = let t_inv = sqrt (1. /. t) in let n = float_of_int maxm in let dm = 0.5 +. n in let f = (pow t_inv power_t_inv ) in match classify_float f with - | FP_normal -> (incomplete_gamma ~alpha:dm t) *. 0.5 *. f - | FP_zero - | FP_subnormal -> 0. - | _ -> raise Exit + | FP_normal -> (incomplete_gamma ~alpha:dm t) *. 0.5 *. f + | FP_zero + | FP_subnormal -> 0. + | _ -> raise Exit in let emt = exp (-. t) in result.(maxm) <- fmax; @@ -233,17 +246,18 @@ let boys_function ~maxm t = result.(n) <- ( (t+.t) *. result.(n+1) +. emt) *. result.(n) done; result - with Exit -> result - end + with Exit -> result + end +(* Functions related to the Boys function:4 ends here *) -let of_some = function - | Some a -> a - | None -> assert false + +(* | ~list_some~ | Filters out all ~None~ elements of the list, and returns the elements without the ~Some~ | + * | ~list_range~ | Creates a list of consecutive integers | + * | ~list_pack~ | ~list_pack n l~ Creates a list of ~n~-elements lists | *) -(** {2 List functions} *) - +(* [[file:../util.org::*List functions][List functions:2]] *) let list_some l = List.filter (function None -> false | _ -> true) l |> List.rev_map (function Some x -> x | _ -> assert false) @@ -252,30 +266,57 @@ let list_some l = let list_range first last = if last < first then [] else - let rec aux accu = function - | 0 -> first :: accu - | i -> (aux [@tailcall]) ( (first+i)::accu ) (i-1) - in - aux [] (last-first) + let rec aux accu = function + | 0 -> first :: accu + | i -> (aux [@tailcall]) ( (first+i)::accu ) (i-1) + in + aux [] (last-first) let list_pack n l = assert (n>=0); let rec aux i accu1 accu2 = function - | [] -> if accu1 = [] then - List.rev accu2 - else - List.rev ((List.rev accu1) :: accu2) - | a :: rest -> - match i with - | 0 -> (aux [@tailcall]) (n-1) [] ((List.rev (a::accu1)) :: accu2) rest - | _ -> (aux [@tailcall]) (i-1) (a::accu1) accu2 rest + | [] -> if accu1 = [] then + List.rev accu2 + else + List.rev ((List.rev accu1) :: accu2) + | a :: rest -> + match i with + | 0 -> (aux [@tailcall]) (n-1) [] ((List.rev (a::accu1)) :: accu2) rest + | _ -> (aux [@tailcall]) (i-1) (a::accu1) accu2 rest in aux (n-1) [] [] l +(* List functions:2 ends here *) -(** {2 Stream functions} *) +(* | ~array_range~ | Creates an array of consecutive integers | + * | ~array_sum~ | Returns the sum of all the elements of the array | + * | ~array_product~ | Returns the product of all the elements of the array | *) + + +(* [[file:../util.org::*Array functions][Array functions:2]] *) +let array_range first last = + if last < first then [| |] else + Array.init (last-first+1) (fun i -> i+first) + + +let array_sum a = + Array.fold_left ( +. ) 0. a + + +let array_product a = + Array.fold_left ( *. ) 1. a +(* Array functions:2 ends here *) + + + +(* | ~stream_range~ | Creates a stream returning consecutive integers | + * | ~stream_to_list~ | Read a stream and put items in a list | + * | ~stream_fold~ | Apply a fold to the elements of the stream | *) + + +(* [[file:../util.org::*Stream functions][Stream functions:2]] *) let stream_range first last = Stream.from (fun i -> let result = i+first in @@ -303,41 +344,55 @@ let stream_fold f init stream = try let element = Stream.next stream in Some (f accu element) - with Stream.Failure -> None + with Stream.Failure -> None in match new_accu with | Some new_accu -> (aux [@tailcall]) new_accu | None -> accu in aux init - -(** {2 Array functions} *) - -let array_range first last = - if last < first then [| |] else - Array.init (last-first+1) (fun i -> i+first) - - -let array_sum a = - Array.fold_left ( +. ) 0. a - -let array_product a = - Array.fold_left ( *. ) 1. a +(* Stream functions:2 ends here *) -(** {2 Printers} *) +(* | ~pp_float_array~ | Printer for float arrays | + * | ~pp_float_array_size~ | Printer for float arrays with size | + * | ~pp_float_2darray~ | Printer for matrices | + * | ~pp_float_2darray_size~ | Printer for matrices with size | + * | ~pp_bitstring~ | Printer for bit strings (used by ~Bitstring~ module) | + * + * #+begin_example + * pp_float_array_size: + * [ 6: 1.000000 1.732051 1.732051 1.000000 1.732051 1.000000 ] + * + * pp_float_array: + * [ 1.000000 1.732051 1.732051 1.000000 1.732051 1.000000 ] + * + * pp_float_2darray_size + * [ + * 2:[ 6: 1.000000 1.732051 1.732051 1.000000 1.732051 1.000000 ] + * [ 4: 1.000000 2.000000 3.000000 4.000000 ] ] + * + * pp_float_2darray: + * [ [ 1.000000 1.732051 1.732051 1.000000 1.732051 1.000000 ] + * [ 1.000000 2.000000 3.000000 4.000000 ] ] + * + * pp_bitstring 14: + * +++++------+-- + * #+end_example *) + + +(* [[file:../util.org::*Printers][Printers:2]] *) +let pp_float_array ppf a = + Format.fprintf ppf "@[<2>[@ "; + Array.iter (fun f -> Format.fprintf ppf "@[%10f@]@ " f) a; + Format.fprintf ppf "]@]" let pp_float_array_size ppf a = Format.fprintf ppf "@[<2>@[ %d:@[<2>" (Array.length a); Array.iter (fun f -> Format.fprintf ppf "@[%10f@]@ " f) a; Format.fprintf ppf "]@]@]" -let pp_float_array ppf a = - Format.fprintf ppf "@[<2>[@ "; - Array.iter (fun f -> Format.fprintf ppf "@[%10f@]@ " f) a; - Format.fprintf ppf "]@]" - let pp_float_2darray ppf a = Format.fprintf ppf "@[<2>[@ "; Array.iter (fun f -> Format.fprintf ppf "@[%a@]@ " pp_float_array f) a; @@ -348,11 +403,7 @@ let pp_float_2darray_size ppf a = Array.iter (fun f -> Format.fprintf ppf "@[%a@]@ " pp_float_array_size f) a; Format.fprintf ppf "]@]@]" - let pp_bitstring n ppf bs = String.init n (fun i -> if (Z.testbit bs i) then '+' else '-') - |> Format.fprintf ppf "@[%s@]" - - - - + |> Format.fprintf ppf "@[%s@]" +(* Printers:2 ends here *) diff --git a/common/lib/util.mli b/common/lib/util.mli index fb02c2b..c9746ac 100644 --- a/common/lib/util.mli +++ b/common/lib/util.mli @@ -1,164 +1,90 @@ -(** All utilities which should be included in all source files are defined here *) +(* [[file:../util.org::*Erf][Erf:2]] *) +external erf_float : float -> float = "erf_float_bytecode" "erf_float" [@@unboxed] [@@noalloc] +(* Erf:2 ends here *) +(* [[file:../util.org::*Erfc][Erfc:2]] *) +external erfc_float : float -> float = "erfc_float_bytecode" "erfc_float" [@@unboxed] [@@noalloc] +(* Erfc:2 ends here *) +(* [[file:../util.org::*Gamma][Gamma:2]] *) +external gamma_float : float -> float = "gamma_float_bytecode" "gamma_float" [@@unboxed] [@@noalloc] +(* Gamma:2 ends here *) -(** {2 Functions from libm} *) - -external erf_float : float -> float = "erf_float_bytecode" "erf_float" - [@@unboxed] [@@noalloc] - (** Error function [erf] from [libm] *) - -external erfc_float : float -> float = "erfc_float_bytecode" "erfc_float" - [@@unboxed] [@@noalloc] - (** Complementary error function [erfc] from [libm] *) - -external gamma_float : float -> float = "gamma_float_bytecode" "gamma_float" - [@@unboxed] [@@noalloc] - (** Gamma function [gamma] from [libm] *) - -external vfork : unit -> int = "unix_vfork" "unix_vfork" - -(* -external popcnt : int64 -> int32 = "popcnt_bytecode" "popcnt" - [@@unboxed] [@@noalloc] - (** popcnt instruction *) - -external trailz : int64 -> int32 = "trailz_bytecode" "trailz" - [@@unboxed] [@@noalloc] - (** ctz instruction *) -*) - +(* [[file:../util.org::*Popcnt][Popcnt:2]] *) val popcnt : int64 -> int -(** popcnt instruction *) +(* Popcnt:2 ends here *) +(* [[file:../util.org::*Trailz][Trailz:2]] *) val trailz : int64 -> int -(** ctz instruction *) +(* Trailz:2 ends here *) +(* [[file:../util.org::*Leadz][Leadz:2]] *) val leadz : int64 -> int -(** ctz instruction *) +(* Leadz:2 ends here *) + +(* General functions *) +(* [[file:../util.org::*General functions][General functions:1]] *) +val fact : int -> float +(* @raise Invalid_argument for negative arguments or arguments >100. *) +val binom : int -> int -> int +val binom_float : int -> int -> float -(** {2 General functions} *) +val chop : float -> (unit -> float) -> float +val pow : float -> int -> float +val float_of_int_fast : int -> float val not_implemented : unit -> 'a -(** Fails with error message if some functionality is not implemented. *) - -val fact : int -> float -(** Factorial function. - @raise Invalid_argument for negative arguments or arguments >100. - *) - -val binom : int -> int -> int -(** Binomial coefficient. [binom n k] {% $= C_n^k$ %}. *) - -val binom_float : int -> int -> float -(** Binomial coefficient. [binom n k] {% $= C_n^k$ %}. *) - -val pow : float -> int -> float -(** Fast implementation of the power function for small integer powers *) - -val chop : float -> (unit -> float) -> float -(** In [chop a f], evaluate [f] only if the absolute value of [a] is larger - than {!Constants.epsilon}, and return [a *. f ()]. - *) - -val float_of_int_fast : int -> float -(* Faster implementation of float_of_int for small positive ints *) - val of_some : 'a option -> 'a +(* General functions:1 ends here *) -(** {2 Functions related to the Boys function} *) +(* Functions related to the Boys function *) + +(* [[file:../util.org::*Functions related to the Boys function][Functions related to the Boys function:1]] *) val incomplete_gamma : alpha:float -> float -> float -(** {{:https://en.wikipedia.org/wiki/Incomplete_gamma_function} - Lower incomplete gamma function} - @raise Failure when the calculation doesn't converge. - *) +(* @raise Failure when the calculation doesn't converge. *) +(* Functions related to the Boys function:1 ends here *) +(* [[file:../util.org::*Functions related to the Boys function][Functions related to the Boys function:3]] *) val boys_function : maxm:int -> float -> float array -(** {{:https://link.springer.com/article/10.1007/s10910-005-9023-3} - Generalized Boys function}. - @param maxm Maximum total angular momentum. - *) +(* Functions related to the Boys function:3 ends here *) +(* List functions *) + -(** {2 Extension of the Array module} *) - -val array_range : int -> int -> int array -(** [array_range first last] returns an array [| first; first+1 ; ... ; last-1 ; last |]. *) - -val array_sum : float array -> float -(** Returns the sum of all the elements of the array *) - -val array_product : float array -> float -(** Returns the product of all the elements of the array *) - - -(** {2 Extension of the List module} *) - -val list_some : 'a option list -> 'a list -(** Filters out all [None] elements of the list, and returns the elements without - the [Some]. *) - +(* [[file:../util.org::*List functions][List functions:1]] *) +val list_some : 'a option list -> 'a list val list_range : int -> int -> int list -(** [list_range first last] returns a list [first; first+1 ; ... ; last-1 ; last ]. *) +val list_pack : int -> 'a list -> 'a list list +(* List functions:1 ends here *) -val list_pack : int -> 'a list -> 'a list list -(** Example: -{[ - list_pack 3 [ 1; 2; 3; ...; 18; 19; 20 ] = - [[1; 2; 3]; [4; 5; 6]; [7; 8; 9]; [10; 11; 12]; [13; 14; 15]; - [16; 17; 18]; [19; 20]] -]} -*) +(* Array functions *) + -(** {2 Useful streams} *) -val stream_range : int -> int -> int Stream.t -(** [stream_range first last] returns a stream . *) +(* [[file:../util.org::*Array functions][Array functions:1]] *) +val array_range : int -> int -> int array +val array_sum : float array -> float +val array_product : float array -> float +(* Array functions:1 ends here *) +(* Stream functions *) + + +(* [[file:../util.org::*Stream functions][Stream functions:1]] *) +val stream_range : int -> int -> int Stream.t val stream_to_list : 'a Stream.t -> 'a list -(** Read a stream and put items in a list. *) +val stream_fold : ('a -> 'b -> 'a) -> 'a -> 'b Stream.t -> 'a +(* Stream functions:1 ends here *) -val stream_fold : ('a -> 'b -> 'a) -> 'a -> 'b Stream.t -> 'a -(** Apply a fold to the elements of the stream. *) +(* Printers *) -(** {2 Printers} *) -val pp_float_array_size : Format.formatter -> float array -> unit -(** Example: -{[ -[ 6: 1.000000 1.732051 1.732051 1.000000 1.732051 1.000000 - ] -]} -*) - -val pp_float_array : Format.formatter -> float array -> unit -(** Example: -{[ -[ 1.000000 1.732051 1.732051 1.000000 1.732051 1.000000 - ] -]} -*) - +(* [[file:../util.org::*Printers][Printers:1]] *) +val pp_float_array_size : Format.formatter -> float array -> unit +val pp_float_array : Format.formatter -> float array -> unit val pp_float_2darray_size : Format.formatter -> float array array -> unit -(** Example: -{[ -[ - 2:[ 6: 1.000000 1.732051 1.732051 1.000000 1.732051 1.000000 ] - [ 4: 1.000000 2.000000 3.000000 4.000000 ] ] -]} -*) - -val pp_float_2darray : Format.formatter -> float array array -> unit -(** Example: -{[ -[ [ 1.000000 1.732051 1.732051 1.000000 1.732051 1.000000 ] - [ 1.000000 2.000000 3.000000 4.000000 ] ] -]} -*) - -val pp_bitstring : int -> Format.formatter -> Z.t -> unit -(** Example: [ pp_bitstring 14 ppf z -> +++++------+-- ] *) - - +val pp_float_2darray : Format.formatter -> float array array -> unit +val pp_bitstring : int -> Format.formatter -> Z.t -> unit +(* Printers:1 ends here *) diff --git a/common/lib/zkey.ml b/common/lib/zkey.ml index c044192..79a6f8f 100644 --- a/common/lib/zkey.ml +++ b/common/lib/zkey.ml @@ -1,5 +1,4 @@ -open Powers - +(* [[file:../zkey.org::*Types][Types:2]] *) type t = { mutable left : int; @@ -7,7 +6,31 @@ type t = kind : int ; } +open Powers +type kind = + | Three of Powers.t + | Four of (int * int * int * int) + | Six of (Powers.t * Powers.t) + | Nine of (Powers.t * Powers.t * Powers.t) + | Twelve of (Powers.t * Powers.t * Powers.t * Powers.t) +(* Types:2 ends here *) + + + +(* | ~of_powers_three | Create from a ~Powers.t~ | + * | ~of_powers_six | Create from two ~Powers.t~ | + * | ~of_powers_nine | Create from three ~Powers.t~ | + * | ~of_powers_twelve | Create from four ~Powers.t~ | + * | ~of_powers | Create using the ~kind~ type | + * | ~of_int_array | Convert from an ~int~ array | + * | ~of_int_four | Create from four ~ints~ | + * | ~to_int_array | Convert to an ~int~ array | + * | ~to_powers | Convert to an ~Powers.t~ array | + * | ~to_string | Pretty printing | *) + + +(* [[file:../zkey.org::*Conversions][Conversions:2]] *) (** Creates a Zkey. *) let make ~kind right = { left = 0 ; right ; kind } @@ -29,13 +52,6 @@ let (<+) z x = z -type kind = - | Three of Powers.t - | Four of (int * int * int * int) - | Six of (Powers.t * Powers.t) - | Nine of (Powers.t * Powers.t * Powers.t) - | Twelve of (Powers.t * Powers.t * Powers.t * Powers.t) - let of_powers_three { x=a ; y=b ; z=c ; _ } = assert ( let alpha = a lor b lor c in @@ -43,6 +59,7 @@ let of_powers_three { x=a ; y=b ; z=c ; _ } = ); make ~kind:3 a <+ b <+ c + let of_int_four i j k l = assert ( let alpha = i lor j lor k lor l in @@ -50,6 +67,7 @@ let of_int_four i j k l = ); make ~kind:4 i <+ j <+ k <+ l + let of_powers_six { x=a ; y=b ; z=c ; _ } { x=d ; y=e ; z=f ; _ } = assert ( let alpha = a lor b lor c lor d lor e lor f in @@ -57,25 +75,27 @@ let of_powers_six { x=a ; y=b ; z=c ; _ } { x=d ; y=e ; z=f ; _ } = ); make ~kind:6 a << b << c << d << e << f + let of_powers_nine { x=a ; y=b ; z=c ; _ } { x=d ; y=e ; z=f ; _ } - { x=g ; y=h ; z=i ; _ } = + { x=g ; y=h ; z=i ; _ } = assert ( let alpha = a lor b lor c lor d lor e lor f lor g lor h lor i in alpha >= 0 && alpha < (1 lsl 10) ); make ~kind:9 a << b << c << d << e << f - <| g << h << i + <| g << h << i + let of_powers_twelve { x=a ; y=b ; z=c ; _ } { x=d ; y=e ; z=f ; _ } - { x=g ; y=h ; z=i ; _ } { x=j ; y=k ; z=l ; _ } = + { x=g ; y=h ; z=i ; _ } { x=j ; y=k ; z=l ; _ } = assert ( let alpha = a lor b lor c lor d lor e lor f - lor g lor h lor i lor j lor k lor l + lor g lor h lor i lor j lor k lor l in alpha >= 0 && alpha < (1 lsl 10) ); make ~kind:12 a << b << c << d << e << f - <| g << h << i << j << k << l + <| g << h << i << j << k << l let of_powers a = @@ -92,61 +112,61 @@ and mask15 = 0x7fff let of_int_array = function -| [| a ; b ; c ; d |] -> of_int_four a b c d -| _ -> invalid_arg "of_int_array" + | [| a ; b ; c ; d |] -> of_int_four a b c d + | _ -> invalid_arg "of_int_array" (** Transform the Zkey into an int array *) let to_int_array { left ; right ; kind } = match kind with | 3 -> [| - mask15 land (right lsr 30) ; - mask15 land (right lsr 15) ; - mask15 land right - |] + mask15 land (right lsr 30) ; + mask15 land (right lsr 15) ; + mask15 land right + |] | 4 -> [| - mask15 land (right lsr 45) ; - mask15 land (right lsr 30) ; - mask15 land (right lsr 15) ; - mask15 land right - |] + mask15 land (right lsr 45) ; + mask15 land (right lsr 30) ; + mask15 land (right lsr 15) ; + mask15 land right + |] | 6 -> [| - mask10 land (right lsr 50) ; - mask10 land (right lsr 40) ; - mask10 land (right lsr 30) ; - mask10 land (right lsr 20) ; - mask10 land (right lsr 10) ; - mask10 land right - |] - + mask10 land (right lsr 50) ; + mask10 land (right lsr 40) ; + mask10 land (right lsr 30) ; + mask10 land (right lsr 20) ; + mask10 land (right lsr 10) ; + mask10 land right + |] + | 12 -> [| - mask10 land (left lsr 50) ; - mask10 land (left lsr 40) ; - mask10 land (left lsr 30) ; - mask10 land (left lsr 20) ; - mask10 land (left lsr 10) ; - mask10 land left ; - mask10 land (right lsr 50) ; - mask10 land (right lsr 40) ; - mask10 land (right lsr 30) ; - mask10 land (right lsr 20) ; - mask10 land (right lsr 10) ; - mask10 land right - |] + mask10 land (left lsr 50) ; + mask10 land (left lsr 40) ; + mask10 land (left lsr 30) ; + mask10 land (left lsr 20) ; + mask10 land (left lsr 10) ; + mask10 land left ; + mask10 land (right lsr 50) ; + mask10 land (right lsr 40) ; + mask10 land (right lsr 30) ; + mask10 land (right lsr 20) ; + mask10 land (right lsr 10) ; + mask10 land right + |] | 9 -> [| - mask10 land (left lsr 20) ; - mask10 land (left lsr 10) ; - mask10 land left ; - mask10 land (right lsr 50) ; - mask10 land (right lsr 40) ; - mask10 land (right lsr 30) ; - mask10 land (right lsr 20) ; - mask10 land (right lsr 10) ; - mask10 land right - |] + mask10 land (left lsr 20) ; + mask10 land (left lsr 10) ; + mask10 land left ; + mask10 land (right lsr 50) ; + mask10 land (right lsr 40) ; + mask10 land (right lsr 30) ; + mask10 land (right lsr 20) ; + mask10 land (right lsr 10) ; + mask10 land right + |] | _ -> invalid_arg (__FILE__^": to_int_array") @@ -155,67 +175,75 @@ let to_int_array { left ; right ; kind } = let to_powers { left ; right ; kind } = match kind with | 3 -> Three (Powers.of_int_tuple ( - mask15 land (right lsr 30) , - mask15 land (right lsr 15) , - mask15 land right - )) + mask15 land (right lsr 30) , + mask15 land (right lsr 15) , + mask15 land right + )) | 6 -> Six (Powers.of_int_tuple - ( mask10 land (right lsr 50) , - mask10 land (right lsr 40) , - mask10 land (right lsr 30)), + ( mask10 land (right lsr 50) , + mask10 land (right lsr 40) , + mask10 land (right lsr 30)), Powers.of_int_tuple - ( mask10 land (right lsr 20) , - mask10 land (right lsr 10) , - mask10 land right ) - ) - + ( mask10 land (right lsr 20) , + mask10 land (right lsr 10) , + mask10 land right ) + ) + | 12 -> Twelve (Powers.of_int_tuple - ( mask10 land (left lsr 50) , - mask10 land (left lsr 40) , - mask10 land (left lsr 30)), - Powers.of_int_tuple - ( mask10 land (left lsr 20) , - mask10 land (left lsr 10) , - mask10 land left ) , - Powers.of_int_tuple - ( mask10 land (right lsr 50) , - mask10 land (right lsr 40) , - mask10 land (right lsr 30)), - Powers.of_int_tuple - ( mask10 land (right lsr 20) , - mask10 land (right lsr 10) , - mask10 land right ) - ) + ( mask10 land (left lsr 50) , + mask10 land (left lsr 40) , + mask10 land (left lsr 30)), + Powers.of_int_tuple + ( mask10 land (left lsr 20) , + mask10 land (left lsr 10) , + mask10 land left ) , + Powers.of_int_tuple + ( mask10 land (right lsr 50) , + mask10 land (right lsr 40) , + mask10 land (right lsr 30)), + Powers.of_int_tuple + ( mask10 land (right lsr 20) , + mask10 land (right lsr 10) , + mask10 land right ) + ) | 9 -> Nine (Powers.of_int_tuple - ( mask10 land (left lsr 20) , - mask10 land (left lsr 10) , - mask10 land left ) , + ( mask10 land (left lsr 20) , + mask10 land (left lsr 10) , + mask10 land left ) , Powers.of_int_tuple - ( mask10 land (right lsr 50) , - mask10 land (right lsr 40) , - mask10 land (right lsr 30)), + ( mask10 land (right lsr 50) , + mask10 land (right lsr 40) , + mask10 land (right lsr 30)), Powers.of_int_tuple - ( mask10 land (right lsr 20) , - mask10 land (right lsr 10) , - mask10 land right ) + ( mask10 land (right lsr 20) , + mask10 land (right lsr 10) , + mask10 land right ) ) | _ -> invalid_arg (__FILE__^": to_powers") +(* Conversions:2 ends here *) +(* | ~hash~ | Associates a nonnegative integer to any Zkey | + * | ~equal~ | The equal function. True if two Zkeys are equal | + * | ~compare~ | Comparison function, used for sorting | *) + + +(* [[file:../zkey.org::*Functions for hash tables][Functions for hash tables:2]] *) let hash = Hashtbl.hash let equal - { right = r1 ; left = l1 ; kind = k1 } - { right = r2 ; left = l2 ; kind = k2 } = + { right = r1 ; left = l1 ; kind = k1 } + { right = r2 ; left = l2 ; kind = k2 } = r1 = r2 && l1 = l2 && k1 = k2 + let compare - { right = r1 ; left = l1 ; kind = k1 } - { right = r2 ; left = l2 ; kind = k2 } = + { right = r1 ; left = l1 ; kind = k1 } + { right = r2 ; left = l2 ; kind = k2 } = if k1 <> k2 then invalid_arg (__FILE__^": cmp"); if r1 < r2 then -1 else if r1 > r2 then 1 @@ -223,11 +251,17 @@ let compare else if l1 > l2 then 1 else 0 + let to_string { left ; right ; kind } = "< " ^ string_of_int left ^ string_of_int right ^ " | " ^ ( - to_int_array { left ; right ; kind } - |> Array.map string_of_int - |> Array.to_list - |> String.concat ", " + to_int_array { left ; right ; kind } + |> Array.map string_of_int + |> Array.to_list + |> String.concat ", " ) ^ " >" +(* Functions for hash tables:2 ends here *) +(* [[file:../zkey.org::*Printers][Printers:2]] *) +let pp ppf t = + Format.fprintf ppf "@[%s@]" (to_string t) +(* Printers:2 ends here *) diff --git a/common/lib/zkey.mli b/common/lib/zkey.mli index 27475ad..4808047 100644 --- a/common/lib/zkey.mli +++ b/common/lib/zkey.mli @@ -1,77 +1,45 @@ -(** Encodes the powers of x, y, z in a compact form, suitable for being - used as keys in a hash table. +(* Types *) - Internally, the {Zkey.t} is made of two integers, [left] and [right]. - The small integers x, y and z are stored compactly in this 126-bits - space: - - {[ - Left Right - 3 [--------------------------------------------------------------] [------------------|---------------|---------------|---------------] - x y z - - 6 [--------------------------------------------------------------] [---|----------|----------|----------|----------|----------|---------] - x1 y1 z1 x2 y2 z2 - - 9 [---------------------------------|----------|----------|---------] [---|----------|----------|----------|----------|----------|---------] - x1 y1 z1 x2 y2 z2 x3 y3 z3 - -12 [---|----------|----------|----------|----------|----------|---------] [---|----------|----------|----------|----------|----------|---------] - x1 y1 z1 x2 y2 z2 x3 y3 z3 x4 y4 z4 - ]} - - The values of x,y,z should be positive and should not exceed 32767 for - [kind=3]. For all other kinds kinds the values should not exceed 1023. - -*) +(* [[file:../zkey.org::*Types][Types:1]] *) type t -val to_string : t -> string -(** Pretty printing *) - -val of_powers_three : Powers.t -> t -(** Create from a {!Powers.t}. *) - -val of_int_four : int -> int -> int -> int -> t -(** Create from four integers. *) - -val of_powers_six : Powers.t -> Powers.t -> t -(** Create from two {!Powers.t}. *) - -val of_powers_nine : Powers.t -> Powers.t -> Powers.t -> t -(** Create from three {!Powers.t}. *) - -val of_powers_twelve : Powers.t -> Powers.t -> Powers.t -> Powers.t -> t -(** Create from four {!Powers.t}. *) - type kind = | Three of Powers.t | Four of (int * int * int * int) | Six of (Powers.t * Powers.t) | Nine of (Powers.t * Powers.t * Powers.t) | Twelve of (Powers.t * Powers.t * Powers.t * Powers.t) +(* Types:1 ends here *) -val of_powers : kind -> t -(** Create using the [kind] type *) +(* Conversions *) -val to_int_array : t -> int array -(** Convert to an int array. *) -val of_int_array : int array -> t -(** Convert from an int array. *) +(* [[file:../zkey.org::*Conversions][Conversions:1]] *) +val of_powers_three : Powers.t -> t +val of_powers_six : Powers.t -> Powers.t -> t +val of_powers_nine : Powers.t -> Powers.t -> Powers.t -> t +val of_powers_twelve : Powers.t -> Powers.t -> Powers.t -> Powers.t -> t +val of_powers : kind -> t +val of_int_array : int array -> t +val of_int_four : int -> int -> int -> int -> t +val to_int_array : t -> int array +val to_powers : t -> kind +val to_string : t -> string +(* Conversions:1 ends here *) -val to_powers : t -> kind +(* Functions for hash tables *) -(** {1 Functions for hash tables} *) - -val hash : t -> int -(** Associates a nonnegative integer to any Zkey. *) - -val equal : t -> t -> bool -(** The equal function. True if two Zkeys are equal. *) +(* [[file:../zkey.org::*Functions for hash tables][Functions for hash tables:1]] *) +val hash : t -> int +val equal : t -> t -> bool val compare : t -> t -> int -(** Comparison function, used for sorting. *) +(* Functions for hash tables:1 ends here *) + +(* Printers *) +(* [[file:../zkey.org::*Printers][Printers:1]] *) +val pp : Format.formatter -> t -> unit +(* Printers:1 ends here *) diff --git a/common/lib/zmap.ml b/common/lib/zmap.ml index 96696b1..e62a761 100644 --- a/common/lib/zmap.ml +++ b/common/lib/zmap.ml @@ -1,5 +1,4 @@ -(** Hash table where the keys are of type Zkey.t (tuples of integers) *) - +(* [[file:../zmap.org::*Type][Type:2]] *) module Zmap = Hashtbl.Make(Zkey) include Zmap - +(* Type:2 ends here *) diff --git a/common/lib/zmap.mli b/common/lib/zmap.mli index b2a922a..4625947 100644 --- a/common/lib/zmap.mli +++ b/common/lib/zmap.mli @@ -1,5 +1,6 @@ -(** Hash table where the keys are of type Zkey.t (tuples of integers). *) +(* Type *) + +(* [[file:../zmap.org::*Type][Type:1]] *) include module type of Hashtbl.Make(Zkey) - - +(* Type:1 ends here *) diff --git a/common/powers.org b/common/powers.org new file mode 100644 index 0000000..5319330 --- /dev/null +++ b/common/powers.org @@ -0,0 +1,128 @@ +#+begin_src elisp tangle: no :results none :exports none +(setq pwd (file-name-directory buffer-file-name)) +(setq name (file-name-nondirectory (substring buffer-file-name 0 -4))) +(setq lib (concat pwd "lib/")) +(setq testdir (concat pwd "test/")) +(setq mli (concat lib name ".mli")) +(setq ml (concat lib name ".ml")) +(setq test-ml (concat testdir name ".ml")) +(org-babel-tangle) +#+end_src + +* Powers + :PROPERTIES: + :header-args: :noweb yes :comments both + :END: + + Contains powers of x, y and z describing the polynomials in atomic basis sets. + +** Type + + #+begin_src ocaml :tangle (eval mli) +type t = private { + x : int ; + y : int ; + z : int ; + tot : int ; +} + #+end_src + + ~tot~ always contains ~x+y+z~. + + #+begin_src ocaml :tangle (eval ml) :exports none +type t = { + x : int ; + y : int ; + z : int ; + tot : int ; +} + #+end_src + +** Conversions + + #+begin_src ocaml :tangle (eval mli) +val of_int_tuple : int * int * int -> t +val to_int_tuple : t -> int * int * int + #+end_src + + #+begin_example +Powers.of_int_tuple (2,3,1);; +- : Powers.t = {Qcaml.Common.Powers.x = 2; y = 3; z = 1; tot = 6} + +Powers.(to_int_tuple (of_int_tuple (2,3,1)));; +- : int * int * int = (2, 3, 1) + #+end_example + + #+begin_src ocaml :tangle (eval ml) :exports none +let of_int_tuple t = + let result = + match t with + | (x,y,z) -> { x ; y ; z ; tot=x+y+z } + in + if result.x < 0 || + result.y < 0 || + result.z < 0 || + result.tot < 0 then + invalid_arg (__FILE__^": of_int_tuple"); + result + + +let to_int_tuple { x ; y ; z ; _ } = (x,y,z) + #+end_src + +** Operations + + #+begin_src ocaml :tangle (eval mli) +val get : Coordinate.axis -> t -> int +val incr : Coordinate.axis -> t -> t +val decr : Coordinate.axis -> t -> t + #+end_src + +| ~get~ | Returns the value of the power for $x$, $y$ or $z$ +| ~incr~ | Returns a new ~Powers.t~ with the power on the given axis incremented | +| ~decr~ | Returns a new ~Powers.t~ with the power on the given axis decremented. As opposed to ~of_int_tuple~, the values may become negative| + + #+begin_example +Powers.get Coordinate.Y (Powers.of_int_tuple (2,3,1));; +- : int = 3 + +Powers.incr Coordinate.Y (Powers.of_int_tuple (2,3,1));; +- : Powers.t = {Qcaml.Common.Powers.x = 2; y = 4; z = 1; tot = 7} + +Powers.decr Coordinate.Y (Powers.of_int_tuple (2,3,1));; +- : Powers.t = {Qcaml.Common.Powers.x = 2; y = 2; z = 1; tot = 5} + #+end_example + + + #+begin_src ocaml :tangle (eval ml) :exports none +let get coord t = + match coord with + | Coordinate.X -> t.x + | Coordinate.Y -> t.y + | Coordinate.Z -> t.z + +let incr coord t = + match coord with + | Coordinate.X -> let r = t.x+1 in { t with x = r ; tot = t.tot+1 } + | Coordinate.Y -> let r = t.y+1 in { t with y = r ; tot = t.tot+1 } + | Coordinate.Z -> let r = t.z+1 in { t with z = r ; tot = t.tot+1 } + +let decr coord t = + match coord with + | Coordinate.X -> let r = t.x-1 in { t with x = r ; tot = t.tot-1 } + | Coordinate.Y -> let r = t.y-1 in { t with y = r ; tot = t.tot-1 } + | Coordinate.Z -> let r = t.z-1 in { t with z = r ; tot = t.tot-1 } + #+end_src + + +** Printers + + #+begin_src ocaml :tangle (eval mli) +val pp : Format.formatter -> t -> unit + #+end_src + + #+begin_src ocaml :tangle (eval ml) :exports none +let pp ppf t = + Format.fprintf ppf "@[x^%d + y^%d + z^%d@]" t.x t.y t.z + #+end_src + diff --git a/common/qcaml.org b/common/qcaml.org new file mode 100644 index 0000000..373ade3 --- /dev/null +++ b/common/qcaml.org @@ -0,0 +1,44 @@ +#+begin_src elisp tangle: no :results none :exports none +(setq pwd (file-name-directory buffer-file-name)) +(setq name (file-name-nondirectory (substring buffer-file-name 0 -4))) +(setq lib (concat pwd "lib/")) +(setq testdir (concat pwd "test/")) +(setq mli (concat lib name ".mli")) +(setq ml (concat lib name ".ml")) +(setq test-ml (concat testdir name ".ml")) +(org-babel-tangle) +#+end_src + +* QCaml + :PROPERTIES: + :header-args: :noweb yes :comments both + :END: + + QCaml-specific parameters + + #+begin_src ocaml :tangle (eval mli) +val root : string +val name : string + #+end_src + + | ~root~ | Path to the QCaml source directory | + | ~name~ | ~"QCaml"~ | + + #+begin_src ocaml :tangle (eval ml) :exports none +let name = "QCaml" + +let root = + let rec chop = function + | [] -> [] + | x :: _ as l when x = name -> l + | _ :: rest -> chop rest + in + String.split_on_char Filename.dir_sep.[0] (Sys.getcwd ()) + |> List.rev + |> chop + |> List.rev + |> String.concat Filename.dir_sep + + + #+end_src + diff --git a/common/range.org b/common/range.org new file mode 100644 index 0000000..c8520cd --- /dev/null +++ b/common/range.org @@ -0,0 +1,96 @@ +#+begin_src elisp tangle: no :results none :exports none +(setq pwd (file-name-directory buffer-file-name)) +(setq name (file-name-nondirectory (substring buffer-file-name 0 -4))) +(setq lib (concat pwd "lib/")) +(setq testdir (concat pwd "test/")) +(setq mli (concat lib name ".mli")) +(setq ml (concat lib name ".ml")) +(setq test-ml (concat testdir name ".ml")) +(org-babel-tangle) +#+end_src + +* Range + :PROPERTIES: + :header-args: :noweb yes :comments both + :END: + + A range is a sorted list of integers in an interval. + + - ~"[a-b]"~ : range between a and b (included) + - ~"[a]"~ : the list with only one integer a + - ~"a"~ : equivalent to "[a]" + - ~"[36-53,72-107,126-131]"~ represents the list of integers + [ 37 ; 37 ; 38 ; ... ; 52 ; 53 ; 72 ; 73 ; ... ; 106 ; 107 ; 126 ; 127 ; ... ; 130 ; 131 ]. + +** Type + + #+begin_src ocaml :tangle (eval mli) +type t + #+end_src + + #+begin_src ocaml :tangle (eval ml) :exports none +type t = int list + #+end_src + +** Conversion + + #+begin_src ocaml :tangle (eval mli) +val of_string : string -> t +val to_string : t -> string +val to_int_list : t -> int list + #+end_src + + #+begin_src ocaml :tangle (eval ml) :exports none +let to_int_list r = r + +let expand_range r = + match String.split_on_char '-' r with + | s :: f :: [] -> + begin + let start = int_of_string s + and finish = int_of_string f + in + assert (start <= finish) ; + let rec do_work = function + | i when i=finish -> [ i ] + | i -> i::(do_work (i+1)) + in do_work start + end + | r :: [] -> [int_of_string r] + | [] -> [] + | _ -> invalid_arg "Only one range expected" + + +let of_string s = + match s.[0] with + | '0' .. '9' -> [ int_of_string s ] + | _ -> + assert (s.[0] = '[') ; + assert (s.[(String.length s)-1] = ']') ; + let s = String.sub s 1 ((String.length s) - 2) in + let l = String.split_on_char ',' s in + let l = List.map expand_range l in + List.concat l + |> List.sort_uniq compare + + +let to_string l = + "[" ^ + (List.map string_of_int l + |> String.concat ",") ^ + "]" + + #+end_src + + +** Printers + + #+begin_src ocaml :tangle (eval mli) +val pp : Format.formatter -> t -> unit + #+end_src + + #+begin_src ocaml :tangle (eval ml) :exports none +let pp ppf t = + Format.fprintf ppf "@[%s@]" (to_string t) + #+end_src + diff --git a/common/spin.org b/common/spin.org new file mode 100644 index 0000000..ca18383 --- /dev/null +++ b/common/spin.org @@ -0,0 +1,63 @@ +#+begin_src elisp tangle: no :results none :exports none +(setq pwd (file-name-directory buffer-file-name)) +(setq name (file-name-nondirectory (substring buffer-file-name 0 -4))) +(setq lib (concat pwd "lib/")) +(setq testdir (concat pwd "test/")) +(setq mli (concat lib name ".mli")) +(setq ml (concat lib name ".ml")) +(setq test-ml (concat testdir name ".ml")) +(org-babel-tangle) +#+end_src + +* Spin + :PROPERTIES: + :header-args: :noweb yes :comments both + :END: + + Electron spin + +** Type + + #+begin_src ocaml :tangle (eval mli) +type t = Alfa | Beta + #+end_src + + Note : + ~Alfa~ if written with an 'f' instead of 'ph' because it has the same number of + letters as ~Beta~, so the alignment of the code is nicer. + + #+begin_src ocaml :tangle (eval ml) :exports none +type t = (* m_s *) + | Alfa (* {% $m_s = +1/2$ %} *) + | Beta (* {% $m_s = -1/2$ %} *) + #+end_src + +** Functions + + #+begin_src ocaml :tangle (eval mli) +val other : t -> t + #+end_src + + Returns the opposite spin + + #+begin_src ocaml :tangle (eval ml) :exports none +let other = function + | Alfa -> Beta + | Beta -> Alfa + +let to_string = function + | Alfa -> "Alpha" + | Beta -> "Beta " + #+end_src + +** Printers + + #+begin_src ocaml :tangle (eval mli) +val pp : Format.formatter -> t -> unit + #+end_src + + #+begin_src ocaml :tangle (eval ml) :exports none +let pp ppf t = + Format.fprintf ppf "@[%s@]" (to_string t) + #+end_src + diff --git a/common/test/math_functions.ml b/common/test/util.ml similarity index 69% rename from common/test/math_functions.ml rename to common/test/util.ml index 3e9c8a2..6058546 100644 --- a/common/test/math_functions.ml +++ b/common/test/util.ml @@ -1,6 +1,15 @@ +(* Test header :noexport: *) + + +(* [[file:../util.org::*Test header][Test header:1]] *) open Common.Util open Alcotest +(* Test header:1 ends here *) +(* Test *) + + +(* [[file:../util.org::*Test][Test:1]] *) let test_external () = check (float 1.e-15) "erf" 0.842700792949715 (erf_float 1.0); check (float 1.e-15) "erf" 0.112462916018285 (erf_float 0.1); @@ -24,7 +33,9 @@ let test_external () = check int "leadz" 63 (leadz @@ Int64.of_int 1); check int "leadz" 64 (leadz @@ Int64.of_int 0); () - +(* Test:1 ends here *) + +(* [[file:../util.org::*General functions][General functions:3]] *) let test_general () = check int "of_some_of_int_fast" 1 (of_some (Some 1)) ; check int "binom" 35 (binom 7 4); @@ -33,7 +44,9 @@ let test_general () = check (float 1.e-15) "pow" 729.0 (pow 3.0 6); check (float 1.e-15) "float_of_int_fast" 10.0 (float_of_int_fast 10); () +(* General functions:3 ends here *) +(* [[file:../util.org::*Functions related to the Boys function][Functions related to the Boys function:5]] *) let test_boys () = check (float 1.e-15) "incomplete_gamma" 0.0 (incomplete_gamma ~alpha:0.5 0.); check (float 1.e-15) "incomplete_gamma" 1.114707979049507 (incomplete_gamma ~alpha:0.5 0.4); @@ -48,27 +61,36 @@ let test_boys () = check (float 1.e-15) "boys" 0.14075053682591263 (boys_function ~maxm:2 0.5).(2); check (float 1.e-15) "boys" 0.00012711171070276764 (boys_function ~maxm:3 15.).(3); () +(* Functions related to the Boys function:5 ends here *) +(* [[file:../util.org::*List functions][List functions:3]] *) +let test_list () = + check bool "list_range" true ([ 2; 3; 4 ] = list_range 2 4); + check bool "list_some" true ([ 2; 3; 4 ] = + list_some ([ None ; Some 2 ; None ; Some 3 ; None ; None ; Some 4]) ); + check bool "list_pack" true (list_pack 3 (list_range 1 20) = + [[1; 2; 3]; [4; 5; 6]; [7; 8; 9]; [10; 11; 12]; [13; 14; 15]; + [16; 17; 18]; [19; 20]]); + () +(* List functions:3 ends here *) + +(* [[file:../util.org::*Array functions][Array functions:3]] *) let test_array () = check bool "array_range" true ([| 2; 3; 4 |] = array_range 2 4); check (float 1.e-15) "array_sum" 9. (array_sum [| 2.; 3.; 4. |]); check (float 1.e-15) "array_product" 24. (array_product [| 2.; 3.; 4. |]); () +(* Array functions:3 ends here *) -let test_list () = - check bool "list_range" true ([ 2; 3; 4 ] = list_range 2 4); - check bool "list_some" true ([ 2; 3; 4 ] = - list_some ([ None ; Some 2 ; None ; Some 3 ; None ; None ; Some 4]) ); - check bool "list_pack" true (list_pack 3 (list_range 1 20) = - [[1; 2; 3]; [4; 5; 6]; [7; 8; 9]; [10; 11; 12]; [13; 14; 15]; - [16; 17; 18]; [19; 20]]); - () +(* Test footer :noexport: *) +(* [[file:../util.org::*Test footer][Test footer:1]] *) let tests = [ - "External", `Quick, test_external; - "General", `Quick, test_general; - "Array", `Quick, test_array; - "List", `Quick, test_list; - "Boys", `Quick, test_boys; + "External", `Quick, test_external; + "General" , `Quick, test_general; + "Boys" , `Quick, test_boys; + "List" , `Quick, test_list; + "Array" , `Quick, test_array; ] +(* Test footer:1 ends here *) diff --git a/common/util.org b/common/util.org new file mode 100644 index 0000000..62cbcad --- /dev/null +++ b/common/util.org @@ -0,0 +1,683 @@ +#+begin_src elisp tangle: no :results none :exports none +(setq pwd (file-name-directory buffer-file-name)) +(setq name (file-name-nondirectory (substring buffer-file-name 0 -4))) +(setq lib (concat pwd "lib/")) +(setq testdir (concat pwd "test/")) +(setq mli (concat lib name ".mli")) +(setq ml (concat lib name ".ml")) +(setq c (concat lib name ".c")) +(setq test-ml (concat testdir name ".ml")) +(org-babel-tangle) +#+end_src + +* Util + :PROPERTIES: + :header-args: :noweb yes :comments both + :END: + + Utility functions. + + +** Test header :noexport: + + #+begin_src ocaml :tangle (eval test-ml) :exports none +open Common.Util +open Alcotest + #+end_src + +** External C functions + + | ~erf_float~ | Error function ~erf~ from =libm= | + | ~erfc_float~ | Complementary error function ~erfc~ from =libm= | + | ~gamma_float~ | Gamma function ~gamma~ from =libm= | + | ~popcnt~ | ~popcnt~ instruction | + | ~trailz~ | ~ctz~ instruction | + | ~leadz~ | ~bsf~ instruction | + + #+begin_src c :tangle (eval c) :exports none +#include +#include +#include + #+end_src + +*** Erf + + #+begin_src c :tangle (eval c) :exports none +CAMLprim value erf_float_bytecode(value x) { + return copy_double(erf(Double_val(x))); +} + +CAMLprim double erf_float(double x) { + return erf(x); +} + #+end_src + + #+begin_src ocaml :tangle (eval mli) +external erf_float : float -> float = "erf_float_bytecode" "erf_float" [@@unboxed] [@@noalloc] + #+end_src + + #+begin_src ocaml :tangle (eval ml) :exports none +external erf_float : float -> float = "erf_float_bytecode" "erf_float" [@@unboxed] [@@noalloc] + #+end_src + +*** Erfc + + #+begin_src c :tangle (eval c) :exports none +CAMLprim value erfc_float_bytecode(value x) { + return copy_double(erfc(Double_val(x))); +} + +CAMLprim double erfc_float(double x) { + return erfc(x); +} + #+end_src + + #+begin_src ocaml :tangle (eval mli) +external erfc_float : float -> float = "erfc_float_bytecode" "erfc_float" [@@unboxed] [@@noalloc] + #+end_src + + #+begin_src ocaml :tangle (eval ml) :exports none +external erfc_float : float -> float = "erfc_float_bytecode" "erfc_float" [@@unboxed] [@@noalloc] + #+end_src + +*** Gamma + + #+begin_src c :tangle (eval c) :exports none +CAMLprim value gamma_float_bytecode(value x) { + return copy_double(tgamma(Double_val(x))); +} + + +CAMLprim double gamma_float(double x) { + return tgamma(x); +} + #+end_src + + #+begin_src ocaml :tangle (eval mli) +external gamma_float : float -> float = "gamma_float_bytecode" "gamma_float" [@@unboxed] [@@noalloc] + #+end_src + + #+begin_src ocaml :tangle (eval ml) :exports none +external gamma_float : float -> float = "gamma_float_bytecode" "gamma_float" [@@unboxed] [@@noalloc] + #+end_src + +*** Popcnt + + #+begin_src c :tangle (eval c) :exports none +CAMLprim int32_t popcnt(int64_t i) { + return __builtin_popcountll (i); +} + + +CAMLprim value popcnt_bytecode(value i) { + return caml_copy_int32(__builtin_popcountll (Int64_val(i))); +} + #+end_src + + #+begin_src ocaml :tangle (eval mli) +val popcnt : int64 -> int + #+end_src + + #+begin_src ocaml :tangle (eval ml) :exports none +external popcnt : int64 -> int32 = "popcnt_bytecode" "popcnt" +[@@unboxed] [@@noalloc] + +let popcnt i = (popcnt [@inlined] ) i |> Int32.to_int + #+end_src + +*** Trailz + + #+begin_src c :tangle (eval c) :exports none +CAMLprim int32_t trailz(int64_t i) { + return __builtin_ctzll (i); +} + + +CAMLprim value trailz_bytecode(value i) { + return caml_copy_int32(__builtin_ctzll (Int64_val(i))); +} + #+end_src + + #+begin_src ocaml :tangle (eval mli) +val trailz : int64 -> int + #+end_src + + #+begin_src ocaml :tangle (eval ml) :exports none +external trailz : int64 -> int32 = "trailz_bytecode" "trailz" "int" +[@@unboxed] [@@noalloc] + +let trailz i = trailz i |> Int32.to_int + #+end_src + +*** Leadz + + #+begin_src c :tangle (eval c) :exports none +CAMLprim int32_t leadz(int64_t i) { + return __builtin_clzll(i); +} + + +CAMLprim value leadz_bytecode(value i) { + return caml_copy_int32(__builtin_clzll (Int64_val(i))); +} + #+end_src + + #+begin_src ocaml :tangle (eval mli) +val leadz : int64 -> int + #+end_src + + #+begin_src ocaml :tangle (eval ml) :exports none +external leadz : int64 -> int32 = "leadz_bytecode" "leadz" "int" +[@@unboxed] [@@noalloc] + +let leadz i = leadz i |> Int32.to_int + #+end_src + +*** Test + + #+begin_src ocaml :tangle (eval test-ml) :exports none +let test_external () = + check (float 1.e-15) "erf" 0.842700792949715 (erf_float 1.0); + check (float 1.e-15) "erf" 0.112462916018285 (erf_float 0.1); + check (float 1.e-15) "erf" (-0.112462916018285) (erf_float (-0.1)); + check (float 1.e-15) "erfc" 0.157299207050285 (erfc_float 1.0); + check (float 1.e-15) "erfc" 0.887537083981715 (erfc_float 0.1); + check (float 1.e-15) "erfc" (1.112462916018285) (erfc_float (-0.1)); + check (float 1.e-14) "gamma" (1.77245385090552) (gamma_float 0.5); + check (float 1.e-14) "gamma" (9.51350769866873) (gamma_float (0.1)); + check (float 1.e-14) "gamma" (-3.54490770181103) (gamma_float (-0.5)); + check int "popcnt" 6 (popcnt @@ Int64.of_int 63); + check int "popcnt" 8 (popcnt @@ Int64.of_int 299605); + check int "popcnt" 1 (popcnt @@ Int64.of_int 65536); + check int "popcnt" 0 (popcnt @@ Int64.of_int 0); + check int "trailz" 3 (trailz @@ Int64.of_int 8); + check int "trailz" 2 (trailz @@ Int64.of_int 12); + check int "trailz" 0 (trailz @@ Int64.of_int 1); + check int "trailz" 64 (trailz @@ Int64.of_int 0); + check int "leadz" 60 (leadz @@ Int64.of_int 8); + check int "leadz" 60 (leadz @@ Int64.of_int 12); + check int "leadz" 63 (leadz @@ Int64.of_int 1); + check int "leadz" 64 (leadz @@ Int64.of_int 0); + () + #+end_src + +** General functions + + #+begin_src ocaml :tangle (eval mli) +val fact : int -> float +(* @raise Invalid_argument for negative arguments or arguments >100. *) +val binom : int -> int -> int +val binom_float : int -> int -> float + +val chop : float -> (unit -> float) -> float +val pow : float -> int -> float +val float_of_int_fast : int -> float + +val not_implemented : unit -> 'a +val of_some : 'a option -> 'a + #+end_src + + | ~fact~ | Factorial function. | + | ~binom~ | Binomial coefficient. ~binom n k~ = $C_n^k$ | + | ~binom_float~ | float variant of ~binom~ | + | ~pow~ | Fast implementation of the power function for small integer powers | + | ~chop~ | In ~chop a f~, evaluate ~f~ only if the absolute value of ~a~ is larger than ~Constants.epsilon~, and return ~a *. f ()~. | + | ~float_of_int_fast~ | Faster implementation of float_of_int for small positive ints | + | ~not_implemented~ | Fails with error message if some functionality is not implemented | + | ~of_some~ | Extracts the value of an option | + + #+begin_src ocaml :tangle (eval ml) :exports none +let memo_float_of_int = + Array.init 64 float_of_int + +let float_of_int_fast i = + if Int.logand i 63 = i then + memo_float_of_int.(i) + else + float_of_int i + + +let factmax = 150 + +let fact_memo = + let rec aux accu_l accu = function + | 0 -> (aux [@tailcall]) [1.] 1. 1 + | i when (i = factmax) -> + let x = (float_of_int factmax) *. accu in + List.rev (x::accu_l) + | i -> let x = (float_of_int i) *. accu in + (aux [@tailcall]) (x::accu_l) x (i+1) + in + aux [] 0. 0 + |> Array.of_list + +let fact = function + | i when (i < 0) -> + raise (Invalid_argument "Argument of factorial should be non-negative") + | i when (i > 150) -> + raise (Invalid_argument "Result of factorial is infinite") + | i -> fact_memo.(i) + + +let binom = + let memo = + let m = Array.make_matrix 64 64 0 in + for n=0 to Array.length m - 1 do + m.(n).(0) <- 1; + m.(n).(n) <- 1; + for k=1 to (n - 1) do + m.(n).(k) <- m.(n-1).(k-1) + m.(n-1).(k) + done + done; + m + in + let rec f n k = + assert (k >= 0); + assert (n >= k); + if k = 0 || k = n then + 1 + else if n < 64 then + memo.(n).(k) + else + f (n-1) (k-1) + f (n-1) k + in f + + +let binom_float n k = + binom n k + |> float_of_int_fast + + +let rec pow a = function + | 0 -> 1. + | 1 -> a + | 2 -> a *. a + | 3 -> a *. a *. a + | -1 -> 1. /. a + | n when n > 0 -> + let b = pow a (n / 2) in + b *. b *. (if n mod 2 = 0 then 1. else a) + | n when n < 0 -> (pow [@tailcall]) (1./.a) (-n) + | _ -> assert false + + +let chop f g = + if (abs_float f) < Constants.epsilon then 0. + else f *. (g ()) + + +let not_implemented () = + failwith "Not implemented" + + +let of_some = function + | Some a -> a + | None -> assert false + + #+end_src + + + #+begin_src ocaml :tangle (eval test-ml) :exports none +let test_general () = + check int "of_some_of_int_fast" 1 (of_some (Some 1)) ; + check int "binom" 35 (binom 7 4); + check (float 1.e-15) "fact" 5040. (fact 7); + check (float 1.e-15) "binom_float" 35.0 (binom_float 7 4); + check (float 1.e-15) "pow" 729.0 (pow 3.0 6); + check (float 1.e-15) "float_of_int_fast" 10.0 (float_of_int_fast 10); + () + #+end_src + +** Functions related to the Boys function + + #+begin_src ocaml :tangle (eval mli) +val incomplete_gamma : alpha:float -> float -> float +(* @raise Failure when the calculation doesn't converge. *) + #+end_src + + The lower [[https://en.wikipedia.org/wiki/Incomplete_gamma_function][Incomplete Gamma function]] is implemented : + \[ + \gamma(\alpha,x) = \int_0^x e^{-t} t^{\alpha-1} dt + \] + + p: $\frac{1}{\Gamma(\alpha)} \int_0^x e^{-t} t^{\alpha-1} dt$ + + q: $\frac{1}{\Gamma(\alpha)} \int_x^\infty e^{-t} t^{\alpha-1} dt$ + + Reference : Haruhiko Okumura: C-gengo niyoru saishin algorithm jiten + (New Algorithm handbook in C language) (Gijyutsu hyouron sha, + Tokyo, 1991) p.227 [in Japanese] + + + #+begin_src ocaml :tangle (eval ml) :exports none +let incomplete_gamma ~alpha x = + assert (alpha >= 0.); + assert (x >= 0.); + let a = alpha in + let a_inv = 1./. a in + let gf = gamma_float alpha in + let loggamma_a = log gf in + let rec p_gamma x = + if x >= 1. +. a then 1. -. q_gamma x + else if x = 0. then 0. + else + let rec pg_loop prev res term k = + if k > 1000. then failwith "p_gamma did not converge." + else if prev = res then res + else + let term = term *. x /. (a +. k) in + (pg_loop [@tailcall]) res (res +. term) term (k +. 1.) + in + let r0 = exp (a *. log x -. x -. loggamma_a) *. a_inv in + pg_loop min_float r0 r0 1. + + and q_gamma x = + if x < 1. +. a then 1. -. p_gamma x + else + let rec qg_loop prev res la lb w k = + if k > 1000. then failwith "q_gamma did not converge." + else if prev = res then res + else + let k_inv = 1. /. k in + let kma = (k -. 1. -. a) *. k_inv in + let la, lb = + lb, kma *. (lb -. la) +. (k +. x) *. lb *. k_inv + in + let w = w *. kma in + let prev, res = res, res +. w /. (la *. lb) in + (qg_loop [@tailcall]) prev res la lb w (k +. 1.) + in + let w = exp (a *. log x -. x -. loggamma_a) in + let lb = (1. +. x -. a) in + qg_loop min_float (w /. lb) 1. lb w 2.0 + in + gf *. p_gamma x + #+end_src + + #+begin_src ocaml :tangle (eval mli) +val boys_function : maxm:int -> float -> float array + #+end_src + + The [[https://link.springer.com/article/10.1007/s10910-005-9023-3][Generalized Boys function]] is implemented, + ~maxm~ is the maximum total angular momentum. + + \[ + F_m(x) = \frac{\gamma(m+1/2,x)}{2x^{m+1/2}} + \] + where $\gamma$ is the incomplete gamma function. + + - $F_0(0.) = 1$ + - $F_0(t) = \frac{\sqrt{\pi}}{2\sqrt{t}} \text{erf} ( \sqrt{t} )$ + - $F_m(0.) = \frac{1}{2m+1}$ + - $F_m(t) = \frac{\gamma{m+1/2,t}}{2t^{m+1/2}}$ + - $F_m(t) = \frac{ 2t\, F_{m+1}(t) + e^{-t} }{2m+1}$ + + #+begin_src ocaml :tangle (eval ml) :exports none +let boys_function ~maxm t = + assert (t >= 0.); + match maxm with + | 0 -> + begin + if t = 0. then [| 1. |] else + let sq_t = sqrt t in + [| (Constants.sq_pi_over_two /. sq_t) *. erf_float sq_t |] + end + | _ -> + begin + assert (maxm > 0); + let result = + Array.init (maxm+1) (fun m -> 1. /. float_of_int (2*m+1)) + in + let power_t_inv = (maxm+maxm+1) in + try + let fmax = + let t_inv = sqrt (1. /. t) in + let n = float_of_int maxm in + let dm = 0.5 +. n in + let f = (pow t_inv power_t_inv ) in + match classify_float f with + | FP_normal -> (incomplete_gamma ~alpha:dm t) *. 0.5 *. f + | FP_zero + | FP_subnormal -> 0. + | _ -> raise Exit + in + let emt = exp (-. t) in + result.(maxm) <- fmax; + for n=maxm-1 downto 0 do + result.(n) <- ( (t+.t) *. result.(n+1) +. emt) *. result.(n) + done; + result + with Exit -> result + end + #+end_src + + #+begin_src ocaml :tangle (eval test-ml) :exports none +let test_boys () = + check (float 1.e-15) "incomplete_gamma" 0.0 (incomplete_gamma ~alpha:0.5 0.); + check (float 1.e-15) "incomplete_gamma" 1.114707979049507 (incomplete_gamma ~alpha:0.5 0.4); + check (float 1.e-15) "incomplete_gamma" 1.4936482656248544 (incomplete_gamma ~alpha:0.5 1.); + check (float 1.e-15) "incomplete_gamma" 1.7724401246392805 (incomplete_gamma ~alpha:0.5 10.); + check (float 1.e-15) "incomplete_gamma" 1.7724538509055159 (incomplete_gamma ~alpha:0.5 100.); + + check (float 1.e-15) "boys" 1.0 (boys_function ~maxm:0 0.).(0); + check (float 1.e-15) "boys" 0.2 (boys_function ~maxm:2 0.).(2); + check (float 1.e-15) "boys" (1./.3.) (boys_function ~maxm:2 0.).(1); + check (float 1.e-15) "boys" 0.8556243918921488 (boys_function ~maxm:0 0.5).(0); + check (float 1.e-15) "boys" 0.14075053682591263 (boys_function ~maxm:2 0.5).(2); + check (float 1.e-15) "boys" 0.00012711171070276764 (boys_function ~maxm:3 15.).(3); + () + #+end_src + +** List functions + + #+begin_src ocaml :tangle (eval mli) +val list_some : 'a option list -> 'a list +val list_range : int -> int -> int list +val list_pack : int -> 'a list -> 'a list list + #+end_src + + | ~list_some~ | Filters out all ~None~ elements of the list, and returns the elements without the ~Some~ | + | ~list_range~ | Creates a list of consecutive integers | + | ~list_pack~ | ~list_pack n l~ Creates a list of ~n~-elements lists | + + #+begin_src ocaml :tangle (eval ml) :exports none +let list_some l = + List.filter (function None -> false | _ -> true) l + |> List.rev_map (function Some x -> x | _ -> assert false) + |> List.rev + + +let list_range first last = + if last < first then [] else + let rec aux accu = function + | 0 -> first :: accu + | i -> (aux [@tailcall]) ( (first+i)::accu ) (i-1) + in + aux [] (last-first) + + +let list_pack n l = + assert (n>=0); + let rec aux i accu1 accu2 = function + | [] -> if accu1 = [] then + List.rev accu2 + else + List.rev ((List.rev accu1) :: accu2) + | a :: rest -> + match i with + | 0 -> (aux [@tailcall]) (n-1) [] ((List.rev (a::accu1)) :: accu2) rest + | _ -> (aux [@tailcall]) (i-1) (a::accu1) accu2 rest + in + aux (n-1) [] [] l + + #+end_src + + #+begin_src ocaml :tangle (eval test-ml) :exports none +let test_list () = + check bool "list_range" true ([ 2; 3; 4 ] = list_range 2 4); + check bool "list_some" true ([ 2; 3; 4 ] = + list_some ([ None ; Some 2 ; None ; Some 3 ; None ; None ; Some 4]) ); + check bool "list_pack" true (list_pack 3 (list_range 1 20) = + [[1; 2; 3]; [4; 5; 6]; [7; 8; 9]; [10; 11; 12]; [13; 14; 15]; + [16; 17; 18]; [19; 20]]); + () + #+end_src + +** Array functions + + #+begin_src ocaml :tangle (eval mli) +val array_range : int -> int -> int array +val array_sum : float array -> float +val array_product : float array -> float + #+end_src + + | ~array_range~ | Creates an array of consecutive integers | + | ~array_sum~ | Returns the sum of all the elements of the array | + | ~array_product~ | Returns the product of all the elements of the array | + + #+begin_src ocaml :tangle (eval ml) :exports none +let array_range first last = + if last < first then [| |] else + Array.init (last-first+1) (fun i -> i+first) + + +let array_sum a = + Array.fold_left ( +. ) 0. a + + +let array_product a = + Array.fold_left ( *. ) 1. a + #+end_src + + #+begin_src ocaml :tangle (eval test-ml) :exports none +let test_array () = + check bool "array_range" true ([| 2; 3; 4 |] = array_range 2 4); + check (float 1.e-15) "array_sum" 9. (array_sum [| 2.; 3.; 4. |]); + check (float 1.e-15) "array_product" 24. (array_product [| 2.; 3.; 4. |]); + () + #+end_src + +** Stream functions + + #+begin_src ocaml :tangle (eval mli) +val stream_range : int -> int -> int Stream.t +val stream_to_list : 'a Stream.t -> 'a list +val stream_fold : ('a -> 'b -> 'a) -> 'a -> 'b Stream.t -> 'a + #+end_src + + | ~stream_range~ | Creates a stream returning consecutive integers | + | ~stream_to_list~ | Read a stream and put items in a list | + | ~stream_fold~ | Apply a fold to the elements of the stream | + + #+begin_src ocaml :tangle (eval ml) :exports none +let stream_range first last = + Stream.from (fun i -> + let result = i+first in + if result <= last then + Some result + else None + ) + +let stream_to_list stream = + let rec aux accu = + let new_accu = + try + Some (Stream.next stream :: accu) + with Stream.Failure -> None + in + match new_accu with + | Some new_accu -> (aux [@tailcall]) new_accu + | None -> accu + in List.rev @@ aux [] + + +let stream_fold f init stream = + let rec aux accu = + let new_accu = + try + let element = Stream.next stream in + Some (f accu element) + with Stream.Failure -> None + in + match new_accu with + | Some new_accu -> (aux [@tailcall]) new_accu + | None -> accu + in + aux init + + #+end_src + +** Printers + + #+begin_src ocaml :tangle (eval mli) +val pp_float_array_size : Format.formatter -> float array -> unit +val pp_float_array : Format.formatter -> float array -> unit +val pp_float_2darray_size : Format.formatter -> float array array -> unit +val pp_float_2darray : Format.formatter -> float array array -> unit +val pp_bitstring : int -> Format.formatter -> Z.t -> unit + #+end_src + + | ~pp_float_array~ | Printer for float arrays | + | ~pp_float_array_size~ | Printer for float arrays with size | + | ~pp_float_2darray~ | Printer for matrices | + | ~pp_float_2darray_size~ | Printer for matrices with size | + | ~pp_bitstring~ | Printer for bit strings (used by ~Bitstring~ module) | + + #+begin_example +pp_float_array_size: +[ 6: 1.000000 1.732051 1.732051 1.000000 1.732051 1.000000 ] + +pp_float_array: +[ 1.000000 1.732051 1.732051 1.000000 1.732051 1.000000 ] + +pp_float_2darray_size +[ + 2:[ 6: 1.000000 1.732051 1.732051 1.000000 1.732051 1.000000 ] + [ 4: 1.000000 2.000000 3.000000 4.000000 ] ] + +pp_float_2darray: +[ [ 1.000000 1.732051 1.732051 1.000000 1.732051 1.000000 ] + [ 1.000000 2.000000 3.000000 4.000000 ] ] + +pp_bitstring 14: ++++++------+-- + #+end_example + + #+begin_src ocaml :tangle (eval ml) :exports none +let pp_float_array ppf a = + Format.fprintf ppf "@[<2>[@ "; + Array.iter (fun f -> Format.fprintf ppf "@[%10f@]@ " f) a; + Format.fprintf ppf "]@]" + +let pp_float_array_size ppf a = + Format.fprintf ppf "@[<2>@[ %d:@[<2>" (Array.length a); + Array.iter (fun f -> Format.fprintf ppf "@[%10f@]@ " f) a; + Format.fprintf ppf "]@]@]" + +let pp_float_2darray ppf a = + Format.fprintf ppf "@[<2>[@ "; + Array.iter (fun f -> Format.fprintf ppf "@[%a@]@ " pp_float_array f) a; + Format.fprintf ppf "]@]" + +let pp_float_2darray_size ppf a = + Format.fprintf ppf "@[<2>@[ %d:@[" (Array.length a); + Array.iter (fun f -> Format.fprintf ppf "@[%a@]@ " pp_float_array_size f) a; + Format.fprintf ppf "]@]@]" + +let pp_bitstring n ppf bs = + String.init n (fun i -> if (Z.testbit bs i) then '+' else '-') + |> Format.fprintf ppf "@[%s@]" + #+end_src + + +** Test footer :noexport: + + #+begin_src ocaml :tangle (eval test-ml) :exports none +let tests = [ + "External", `Quick, test_external; + "General" , `Quick, test_general; + "Boys" , `Quick, test_boys; + "List" , `Quick, test_list; + "Array" , `Quick, test_array; +] + #+end_src diff --git a/common/zkey.org b/common/zkey.org new file mode 100644 index 0000000..27c20cf --- /dev/null +++ b/common/zkey.org @@ -0,0 +1,346 @@ +#+begin_src elisp tangle: no :results none :exports none +(setq pwd (file-name-directory buffer-file-name)) +(setq name (file-name-nondirectory (substring buffer-file-name 0 -4))) +(setq lib (concat pwd "lib/")) +(setq testdir (concat pwd "test/")) +(setq mli (concat lib name ".mli")) +(setq ml (concat lib name ".ml")) +(setq test-ml (concat testdir name ".ml")) +(org-babel-tangle) +#+end_src + +* Zkey + :PROPERTIES: + :header-args: :noweb yes :comments both + :END: + + Encodes the powers of x, y, z in a compact form, suitable for being + used as keys in a hash table. + + Internally, the ~Zkey.t~ is made of two integers, ~left~ and ~right~. + The small integers x, y and z are stored compactly in this 126-bits + space: + + #+begin_example + Left Right + 3 [--------------------------------------------------------------] [------------------|---------------|---------------|---------------] + x y z + + 6 [--------------------------------------------------------------] [---|----------|----------|----------|----------|----------|---------] + x1 y1 z1 x2 y2 z2 + + 9 [---------------------------------|----------|----------|---------] [---|----------|----------|----------|----------|----------|---------] + x1 y1 z1 x2 y2 z2 x3 y3 z3 + +12 [---|----------|----------|----------|----------|----------|---------] [---|----------|----------|----------|----------|----------|---------] + x1 y1 z1 x2 y2 z2 x3 y3 z3 x4 y4 z4 + #+end_example + + The values of x,y,z should be positive and should not exceed 32767 for + ~kind=3~. For all other kinds kinds the values should not exceed 1023. + +** Types + + #+begin_src ocaml :tangle (eval mli) +type t + +type kind = + | Three of Powers.t + | Four of (int * int * int * int) + | Six of (Powers.t * Powers.t) + | Nine of (Powers.t * Powers.t * Powers.t) + | Twelve of (Powers.t * Powers.t * Powers.t * Powers.t) + #+end_src + + #+begin_src ocaml :tangle (eval ml) :exports none +type t = + { + mutable left : int; + mutable right : int; + kind : int ; + } + +open Powers + +type kind = + | Three of Powers.t + | Four of (int * int * int * int) + | Six of (Powers.t * Powers.t) + | Nine of (Powers.t * Powers.t * Powers.t) + | Twelve of (Powers.t * Powers.t * Powers.t * Powers.t) + #+end_src + +** Conversions + + #+begin_src ocaml :tangle (eval mli) +val of_powers_three : Powers.t -> t +val of_powers_six : Powers.t -> Powers.t -> t +val of_powers_nine : Powers.t -> Powers.t -> Powers.t -> t +val of_powers_twelve : Powers.t -> Powers.t -> Powers.t -> Powers.t -> t +val of_powers : kind -> t +val of_int_array : int array -> t +val of_int_four : int -> int -> int -> int -> t +val to_int_array : t -> int array +val to_powers : t -> kind +val to_string : t -> string + #+end_src + + | ~of_powers_three | Create from a ~Powers.t~ | + | ~of_powers_six | Create from two ~Powers.t~ | + | ~of_powers_nine | Create from three ~Powers.t~ | + | ~of_powers_twelve | Create from four ~Powers.t~ | + | ~of_powers | Create using the ~kind~ type | + | ~of_int_array | Convert from an ~int~ array | + | ~of_int_four | Create from four ~ints~ | + | ~to_int_array | Convert to an ~int~ array | + | ~to_powers | Convert to an ~Powers.t~ array | + | ~to_string | Pretty printing | + + #+begin_src ocaml :tangle (eval ml) :exports none +(** Creates a Zkey. *) +let make ~kind right = + { left = 0 ; right ; kind } + +(** Move [right] to [left] and set [right = x] *) +let (<|) z x = + z.left <- z.right; + z.right <- x; + z + +(** Shift left [right] by 10 bits, and add [x]. *) +let (<<) z x = + z.right <- (z.right lsl 10) lor x ; + z + +(** Shift left [right] by 10 bits, and add [x]. *) +let (<+) z x = + z.right <- (z.right lsl 15) lor x ; + z + + +let of_powers_three { x=a ; y=b ; z=c ; _ } = + assert ( + let alpha = a lor b lor c in + alpha >= 0 && alpha < (1 lsl 15) + ); + make ~kind:3 a <+ b <+ c + + +let of_int_four i j k l = + assert ( + let alpha = i lor j lor k lor l in + alpha >= 0 && alpha < (1 lsl 15) + ); + make ~kind:4 i <+ j <+ k <+ l + + +let of_powers_six { x=a ; y=b ; z=c ; _ } { x=d ; y=e ; z=f ; _ } = + assert ( + let alpha = a lor b lor c lor d lor e lor f in + alpha >= 0 && alpha < (1 lsl 10) + ); + make ~kind:6 a << b << c << d << e << f + + +let of_powers_nine { x=a ; y=b ; z=c ; _ } { x=d ; y=e ; z=f ; _ } + { x=g ; y=h ; z=i ; _ } = + assert ( + let alpha = a lor b lor c lor d lor e lor f lor g lor h lor i in + alpha >= 0 && alpha < (1 lsl 10) + ); + make ~kind:9 a << b << c << d << e << f + <| g << h << i + + +let of_powers_twelve { x=a ; y=b ; z=c ; _ } { x=d ; y=e ; z=f ; _ } + { x=g ; y=h ; z=i ; _ } { x=j ; y=k ; z=l ; _ } = + assert ( + let alpha = a lor b lor c lor d lor e lor f + lor g lor h lor i lor j lor k lor l + in + alpha >= 0 && alpha < (1 lsl 10) + ); + make ~kind:12 a << b << c << d << e << f + <| g << h << i << j << k << l + + +let of_powers a = + match a with + | Three a -> of_powers_three a + | Six (a,b) -> of_powers_six a b + | Twelve (a,b,c,d) -> of_powers_twelve a b c d + | Nine (a,b,c) -> of_powers_nine a b c + | _ -> invalid_arg "of_powers" + + +let mask10 = 0x3ff +and mask15 = 0x7fff + + +let of_int_array = function + | [| a ; b ; c ; d |] -> of_int_four a b c d + | _ -> invalid_arg "of_int_array" + + +(** Transform the Zkey into an int array *) +let to_int_array { left ; right ; kind } = + match kind with + | 3 -> [| + mask15 land (right lsr 30) ; + mask15 land (right lsr 15) ; + mask15 land right + |] + + | 4 -> [| + mask15 land (right lsr 45) ; + mask15 land (right lsr 30) ; + mask15 land (right lsr 15) ; + mask15 land right + |] + + | 6 -> [| + mask10 land (right lsr 50) ; + mask10 land (right lsr 40) ; + mask10 land (right lsr 30) ; + mask10 land (right lsr 20) ; + mask10 land (right lsr 10) ; + mask10 land right + |] + + | 12 -> [| + mask10 land (left lsr 50) ; + mask10 land (left lsr 40) ; + mask10 land (left lsr 30) ; + mask10 land (left lsr 20) ; + mask10 land (left lsr 10) ; + mask10 land left ; + mask10 land (right lsr 50) ; + mask10 land (right lsr 40) ; + mask10 land (right lsr 30) ; + mask10 land (right lsr 20) ; + mask10 land (right lsr 10) ; + mask10 land right + |] + + | 9 -> [| + mask10 land (left lsr 20) ; + mask10 land (left lsr 10) ; + mask10 land left ; + mask10 land (right lsr 50) ; + mask10 land (right lsr 40) ; + mask10 land (right lsr 30) ; + mask10 land (right lsr 20) ; + mask10 land (right lsr 10) ; + mask10 land right + |] + | _ -> invalid_arg (__FILE__^": to_int_array") + + + +(** Transform the Zkey into an int tuple *) +let to_powers { left ; right ; kind } = + match kind with + | 3 -> Three (Powers.of_int_tuple ( + mask15 land (right lsr 30) , + mask15 land (right lsr 15) , + mask15 land right + )) + + | 6 -> Six (Powers.of_int_tuple + ( mask10 land (right lsr 50) , + mask10 land (right lsr 40) , + mask10 land (right lsr 30)), + Powers.of_int_tuple + ( mask10 land (right lsr 20) , + mask10 land (right lsr 10) , + mask10 land right ) + ) + + | 12 -> Twelve (Powers.of_int_tuple + ( mask10 land (left lsr 50) , + mask10 land (left lsr 40) , + mask10 land (left lsr 30)), + Powers.of_int_tuple + ( mask10 land (left lsr 20) , + mask10 land (left lsr 10) , + mask10 land left ) , + Powers.of_int_tuple + ( mask10 land (right lsr 50) , + mask10 land (right lsr 40) , + mask10 land (right lsr 30)), + Powers.of_int_tuple + ( mask10 land (right lsr 20) , + mask10 land (right lsr 10) , + mask10 land right ) + ) + + | 9 -> Nine (Powers.of_int_tuple + ( mask10 land (left lsr 20) , + mask10 land (left lsr 10) , + mask10 land left ) , + Powers.of_int_tuple + ( mask10 land (right lsr 50) , + mask10 land (right lsr 40) , + mask10 land (right lsr 30)), + Powers.of_int_tuple + ( mask10 land (right lsr 20) , + mask10 land (right lsr 10) , + mask10 land right ) + ) + + | _ -> invalid_arg (__FILE__^": to_powers") + + + #+end_src + +** Functions for hash tables + + #+begin_src ocaml :tangle (eval mli) +val hash : t -> int +val equal : t -> t -> bool +val compare : t -> t -> int + #+end_src + + | ~hash~ | Associates a nonnegative integer to any Zkey | + | ~equal~ | The equal function. True if two Zkeys are equal | + | ~compare~ | Comparison function, used for sorting | + + #+begin_src ocaml :tangle (eval ml) :exports none +let hash = Hashtbl.hash + +let equal + { right = r1 ; left = l1 ; kind = k1 } + { right = r2 ; left = l2 ; kind = k2 } = + r1 = r2 && l1 = l2 && k1 = k2 + + +let compare + { right = r1 ; left = l1 ; kind = k1 } + { right = r2 ; left = l2 ; kind = k2 } = + if k1 <> k2 then invalid_arg (__FILE__^": cmp"); + if r1 < r2 then -1 + else if r1 > r2 then 1 + else if l1 < l2 then -1 + else if l1 > l2 then 1 + else 0 + + +let to_string { left ; right ; kind } = + "< " ^ string_of_int left ^ string_of_int right ^ " | " ^ ( + to_int_array { left ; right ; kind } + |> Array.map string_of_int + |> Array.to_list + |> String.concat ", " + ) ^ " >" + #+end_src + +** Printers + + #+begin_src ocaml :tangle (eval mli) +val pp : Format.formatter -> t -> unit + #+end_src + + #+begin_src ocaml :tangle (eval ml) :exports none +let pp ppf t = + Format.fprintf ppf "@[%s@]" (to_string t) + #+end_src diff --git a/common/zmap.org b/common/zmap.org new file mode 100644 index 0000000..1883da2 --- /dev/null +++ b/common/zmap.org @@ -0,0 +1,29 @@ +#+begin_src elisp tangle: no :results none :exports none +(setq pwd (file-name-directory buffer-file-name)) +(setq name (file-name-nondirectory (substring buffer-file-name 0 -4))) +(setq lib (concat pwd "lib/")) +(setq testdir (concat pwd "test/")) +(setq mli (concat lib name ".mli")) +(setq ml (concat lib name ".ml")) +(setq test-ml (concat testdir name ".ml")) +(org-babel-tangle) +#+end_src + +* Zmap + :PROPERTIES: + :header-args: :noweb yes :comments both + :END: + + A hash table where the keys are ~Zkey~ + +** Type + + #+begin_src ocaml :tangle (eval mli) +include module type of Hashtbl.Make(Zkey) + #+end_src + + #+begin_src ocaml :tangle (eval ml) :exports none +module Zmap = Hashtbl.Make(Zkey) +include Zmap + #+end_src + diff --git a/docs/config.el b/docs/config.el index 6349c11..d1244e4 100755 --- a/docs/config.el +++ b/docs/config.el @@ -72,6 +72,7 @@ with class 'color and highest min-color value." (setq org-html-htmlize-output-type 'css) ; default: 'inline-css (setq org-html-htmlize-font-prefix "org-") ; default: "org-" +(setq c "c") (setq ml "ml") (setq mli "mli") (setq test-ml "test-ml") diff --git a/test/run_tests.ml b/test/run_tests.ml index af83b15..e03aeb1 100644 --- a/test/run_tests.ml +++ b/test/run_tests.ml @@ -4,7 +4,7 @@ let test_suites: unit Alcotest.test list = [ "Common.Bitstring", Test_common.Bitstring.tests; - "Common.Util", Test_common.Math_functions.tests; + "Common.Util", Test_common.Util.tests; "Linear_algebra.Vector", Test_linear_algebra.Vector.tests; "Particles.Nuclei", Test_particles.Nuclei.tests; "Particles.Electrons", Test_particles.Electrons.tests;