mirror of https://gitlab.com/scemama/QCaml.git
Finished common in org-mode
This commit is contained in:
parent
36e7dbc7bd
commit
d2e7848de2
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -11,7 +11,7 @@
|
|||
)
|
||||
|
||||
(c_names
|
||||
math_functions
|
||||
util
|
||||
)
|
||||
(c_flags (:standard)
|
||||
-Ofast -march=native -fPIC
|
||||
|
|
|
|||
|
|
@ -1,85 +0,0 @@
|
|||
#include <math.h>
|
||||
#include <caml/mlvalues.h>
|
||||
#include <caml/alloc.h>
|
||||
|
||||
|
||||
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 <stdio.h>
|
||||
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 <unistd.h>
|
||||
|
||||
CAMLprim value unix_vfork(value unit)
|
||||
{
|
||||
int ret;
|
||||
ret = vfork();
|
||||
return Val_int(ret);
|
||||
}
|
||||
|
||||
|
|
@ -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 *)
|
||||
|
|
|
|||
|
|
@ -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 *)
|
||||
|
|
|
|||
|
|
@ -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 *)
|
||||
|
|
|
|||
|
|
@ -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 *)
|
||||
|
|
|
|||
|
|
@ -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 *)
|
||||
|
|
|
|||
|
|
@ -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 *)
|
||||
|
|
|
|||
|
|
@ -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 *)
|
||||
|
|
|
|||
|
|
@ -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 *)
|
||||
|
|
|
|||
|
|
@ -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 <math.h>
|
||||
#include <caml/mlvalues.h>
|
||||
#include <caml/alloc.h>
|
||||
/* 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 */
|
||||
|
|
@ -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 "@[<h>%s@]"
|
||||
|
||||
|
||||
|
||||
|
||||
|> Format.fprintf ppf "@[<h>%s@]"
|
||||
(* Printers:2 ends here *)
|
||||
|
|
|
|||
|
|
@ -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 <first ; first+1 ; ... ; last-1 ; last>. *)
|
||||
(* [[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 *)
|
||||
|
|
|
|||
|
|
@ -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 *)
|
||||
|
|
|
|||
|
|
@ -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 *)
|
||||
|
|
|
|||
|
|
@ -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 *)
|
||||
|
|
|
|||
|
|
@ -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 *)
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
||||
|
|
@ -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
|
||||
|
||||
|
|
@ -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
|
||||
|
||||
|
|
@ -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
|
||||
|
||||
|
|
@ -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 *)
|
||||
|
|
@ -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 <math.h>
|
||||
#include <caml/mlvalues.h>
|
||||
#include <caml/alloc.h>
|
||||
#+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 "@[<h>%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
|
||||
|
|
@ -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
|
||||
|
|
@ -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
|
||||
|
||||
|
|
@ -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")
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
|
|
|
|||
Loading…
Reference in New Issue