Index: trunk/omega/tests/omega_unit.ml =================================================================== --- trunk/omega/tests/omega_unit.ml (revision 8861) +++ trunk/omega/tests/omega_unit.ml (revision 8862) @@ -1,229 +1,230 @@ (* omega_unit.ml -- Copyright (C) 1999-2023 by Wolfgang Kilian Thorsten Ohl Juergen Reuter Christian Speckner WHIZARD is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2, or (at your option) any later version. WHIZARD is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. *) open OUnit let unattended = ref true let skip_if_unattended () = skip_if !unattended "not suitable for unattended tests" let trivial_test = "trivial" >:: (bracket (fun () -> true) (fun b -> assert_bool "always true" b) (fun b -> ())) let short_random_list n = let l = ref [] in for i = 1 to n do l := Random.int 1024 :: !l done; !l let allowed_recursion_depth () = let rec allowed_recursion_depth' n = try allowed_recursion_depth' (succ n) with | Stack_overflow -> n in allowed_recursion_depth' 0 let long_random_list factor = let n = factor * allowed_recursion_depth () in let l = ref [] in for i = 1 to n do l := Random.int n :: !l done; !l module Integer = struct type t = int let compare = compare let pp_printer = Format.pp_print_int let pp_print_sep = OUnitDiff.pp_comma_separator end module Integer_List = OUnitDiff.ListSimpleMake(Integer) module ThoList_Unit_Tests = struct let inner_list = ThoList.range 1 5 let outer_list = List.map (( * ) 10) (ThoList.range 1 4) let f n = List.map ((+) n) inner_list let flatmap = "flatmap" >:: (fun () -> let result = ThoList.flatmap f outer_list and expected = List.flatten (List.map f outer_list) in assert_equal expected result) let rev_flatmap = "rev_flatmap" >:: (fun () -> let result = ThoList.rev_flatmap f outer_list and expected = List.rev (ThoList.flatmap f outer_list) in Integer_List.assert_equal expected result) let flatmap_stack_overflow = "flatmap_stack_overflow" >:: (fun () -> skip_if !unattended "memory limits not suitable for unattended tests"; let l = long_random_list 2 in let f n = List.map ((+) n) (short_random_list 2) in assert_raises Stack_overflow (fun () -> ThoList.flatmap f l)) let rev_flatmap_no_stack_overflow = "rev_flatmap_no_stack_overflow" >:: (fun () -> skip_if !unattended "memory limits not suitable for unattended tests"; let l = long_random_list 10 in let f n = List.map ((+) n) (short_random_list 10) in ignore (ThoList.rev_flatmap f l); assert_bool "always true" true) let suite = "ThoList" >::: [flatmap; flatmap_stack_overflow; rev_flatmap; rev_flatmap_no_stack_overflow ] end module IListSet = Set.Make (struct type t = int list let compare = compare end) let list_elements_unique l = let rec list_elements_unique' set = function | [] -> true | x :: rest -> if IListSet.mem x set then false else list_elements_unique' (IListSet.add x set) rest in list_elements_unique' IListSet.empty l let ilistset_test = "IListSet" >:: (fun () -> assert_bool "true" (list_elements_unique [[1];[2]]); assert_bool "false" (not (list_elements_unique [[1];[1]]))) module Combinatorics_Unit_Tests = struct let permute = "permute" >:: (fun () -> let n = 8 in let l = ThoList.range 1 n in let result = Combinatorics.permute l in assert_equal (Combinatorics.factorial n) (List.length result); assert_bool "unique" (list_elements_unique result)) let permute_no_stack_overflow = "permute_no_stack_overflow" >:: (fun () -> skip_if !unattended "memory limits not suitable for unattended tests"; let n = 10 in (* n = 10 needs 1 GB, n = 11 needs 7.3 GB *) let l = ThoList.range 1 n in let result = Combinatorics.permute l in assert_equal (Combinatorics.factorial n) (List.length result)) let suite = "Combinatorics" >::: [permute; permute_no_stack_overflow] end let selftest_suite = "testsuite" >::: [trivial_test; ilistset_test] module Permutation_Test_Using_Lists = Permutation.Test(Permutation.Using_Lists) module Permutation_Test_Using_Arrays = Permutation.Test(Permutation.Using_Arrays) let suite = "omega" >::: [selftest_suite; ThoList_Unit_Tests.suite; ThoList.Test.suite; ThoArray.Test.suite; ThoString.Test.suite; Partial.Test.suite; Permutation_Test_Using_Lists.suite; Permutation_Test_Using_Arrays.suite; Combinatorics_Unit_Tests.suite; Combinatorics.Test.suite; Young.Test.suite; Algebra.Q.Test.suite; Algebra.QC.Test.suite; Algebra.Laurent.Test.suite; Color.Flow.Test.suite; Color.Arrow.Test.suite; Color.Birdtracks.Test.suite; Color.SU3.Test.suite; (* [Color.U3.Test.suite;] *) UFO_targets.Fortran.Test.suite; UFO_Lorentz.Test.suite; + UFOx.Test.suite; UFO.Test.suite; Format_Fortran.Test.suite; Dirac.Chiral.test_suite; Dirac.Dirac.test_suite; Dirac.Majorana.test_suite] let suite_long = "omega long" >::: [Young.Test.suite_long; Color.Flow.Test.suite_long; Color.Arrow.Test.suite_long; Color.Birdtracks.Test.suite_long; Color.SU3.Test.suite_long; (* [Color.U3.Test.suite_long] *) ] let run_suite_long = ref false let _ = ignore (run_test_tt_main ~arg_specs:[("-attended", Arg.Clear unattended, " run tests that depend on the environment"); ("-unattended", Arg.Set unattended, " don't run tests depend on the environment"); ("-long", Arg.Set run_suite_long, " also run the very long tests")] suite); if !run_suite_long then ignore (run_test_tt suite_long); exit 0 Index: trunk/omega/src/UFOx.mli =================================================================== --- trunk/omega/src/UFOx.mli (revision 8861) +++ trunk/omega/src/UFOx.mli (revision 8862) @@ -1,250 +1,250 @@ (* vertex.mli -- Copyright (C) 1999-2023 by Wolfgang Kilian Thorsten Ohl Juergen Reuter with contributions from Christian Speckner WHIZARD is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2, or (at your option) any later version. WHIZARD is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. *) module Expr : sig type t val of_string : string -> t val of_strings : string list -> t val substitute : string -> t -> t -> t val rename : (string * string) list -> t -> t val half : string -> t val variables : t -> Sets.String_Caseless.t val functions : t -> Sets.String_Caseless.t end module Value : sig type t val of_expr : Expr.t -> t val to_string : t -> string val to_coupling : (string -> 'b) -> t -> 'b Coupling.expr end (* \begin{dubious} UFO represents rank-2 indices $(i,j)$ as $1000\cdot j + i$. This should be replaced by a proper union type eventually. Unfortunately, this requires many changes in the [Atom]s in [UFOx]. Therefore, we try a quick'n'dirty proof of principle first. \end{dubious} *) module type Index = sig type t = int val position : t -> int val factor : t -> int val unpack : t -> int * int val pack : int -> int -> t val map_position : (int -> int) -> t -> t val to_string : t -> string val list_to_string : t list -> string (* Indices are represented by a pair [int * 'r], where ['r] denotes the representation the index belongs to. *) (* [free indices] returns all free indices in the list [indices], i.\,e.~all positive indices. *) val free : (t * 'r) list -> (t * 'r) list (* [summation indices] returns all summation indices in the list [indices], i.\,e.~all negative indices. *) val summation : (t * 'r) list -> (t * 'r) list val classes_to_string : ('r -> string) -> (t * 'r) list -> string (* Generate summation indices, starting from~$-1001$. TODO: check that there are no clashes with explicitely named indices. *) val fresh_summation : unit -> t val named_summation : string -> unit -> t end module Index : Index module type Tensor = sig type atom (* A tensor is a linear combination of products of [atom]s with rational coefficients. The following could be refined by introducing [scalar] atoms and restricting the denominators to [(scalar list * Algebra.QC.t) list]. At the moment, this restriction is implemented dynamically by [of_expr] and not statically in the type system. Polymorphic variants appear to be the right tool, either directly or as phantom types. However, this is certainly only \textit{nice-to-have} and is not essential. *) type 'a linear = ('a list * Algebra.QC.t) list type t = | Linear of atom linear | Ratios of (atom linear * atom linear) list (* We might need to replace atoms if the syntax is not context free. *) val map_atoms : (atom -> atom) -> t -> t (* We need to rename indices to implement permutations \ldots *) val map_indices : (int -> int) -> t -> t (* \ldots{} but in order to to clean up inconsistencies in the syntax of \texttt{lorentz.py} and \texttt{propagators.py} we also need to rename indices without touching the second argument of \texttt{P}, the argument of \texttt{Mass} etc. *) val rename_indices : (int -> int) -> t -> t (* We need scale coefficients. *) val map_coeff : (Algebra.QC.t -> Algebra.QC.t) -> t -> t (* Try to contract adjacent pairs of [atoms] as allowed but [Atom.contract_pair]. This is not exhaustive, but helps a lot with invariant squares of momenta in applications of [Lorentz]. *) val contract_pairs : t -> t (* The list of variable referenced in the tensor expression, that will need to be imported by the numerical code. *) val variables : t -> string list (* Parsing and unparsing. Lists of [string]s are interpreted as sums. *) val of_expr : UFOx_syntax.expr -> t val of_string : string -> t val of_strings : string list -> t val to_string : t -> string (* The supported representations. *) type r val classify_indices : t -> (int * r) list val rep_to_string : r -> string val rep_to_string_whizard : r -> string val rep_of_int : bool -> int -> r val rep_conjugate : r -> r val rep_trivial : r -> bool (* There is not a 1-to-1 mapping between the representations in the model files and the representations used by O'Mega, e.\,g.~in [Coupling.lorentz]. We might need to use heuristics. *) type r_omega val omega : r -> r_omega end module type Atom = sig type t val map_indices : (int -> int) -> t -> t val rename_indices : (int -> int) -> t -> t val contract_pair : t -> t -> t option val variable : t -> string option val scalar : t -> bool val is_unit : t -> bool val invertible : t -> bool val invert : t -> t val of_expr : string -> UFOx_syntax.expr list -> t list val to_string : t -> string type r val classify_indices : t list -> (int * r) list val disambiguate_indices : t list -> t list val rep_to_string : r -> string val rep_to_string_whizard : r -> string val rep_of_int : bool -> int -> r val rep_conjugate : r -> r val rep_trivial : r -> bool type r_omega val omega : r -> r_omega end module type Lorentz_Atom = sig type dirac = private | C of int * int | Gamma of int * int * int | Gamma5 of int * int | Identity of int * int | ProjP of int * int | ProjM of int * int | Sigma of int * int * int * int type vector = (* private *) | Epsilon of int * int * int * int | Metric of int * int | P of int * int type scalar = (* private *) | Mass of int | Width of int | P2 of int | P12 of int * int | Variable of string | Coeff of Value.t type t = (* private *) | Dirac of dirac | Vector of vector | Scalar of scalar | Inverse of scalar val map_indices_scalar : (int -> int) -> scalar -> scalar val map_indices_vector : (int -> int) -> vector -> vector val rename_indices_vector : (int -> int) -> vector -> vector end module Lorentz_Atom : Lorentz_Atom module Lorentz : Tensor with type atom = Lorentz_Atom.t and type r_omega = Coupling.lorentz module type Color_Atom = sig type t = (* private *) | Identity of int * int | Identity8 of int * int | T of int * int * int | F of int * int * int | D of int * int * int | Epsilon of int * int * int | EpsilonBar of int * int * int | T6 of int * int * int | K6 of int * int * int | K6Bar of int * int * int end module Color_Atom : Color_Atom module Color : Tensor with type atom = Color_Atom.t and type r_omega = Color.t module type Test = sig - val example : unit -> unit val suite : OUnit.test end +module Test : Test Index: trunk/omega/src/UFOx.ml =================================================================== --- trunk/omega/src/UFOx.ml (revision 8861) +++ trunk/omega/src/UFOx.ml (revision 8862) @@ -1,1522 +1,1556 @@ (* vertex.ml -- Copyright (C) 1999-2023 by Wolfgang Kilian Thorsten Ohl Juergen Reuter with contributions from Christian Speckner WHIZARD is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2, or (at your option) any later version. WHIZARD is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. *) let error_in_string text start_pos end_pos = let i = max 0 start_pos.Lexing.pos_cnum in let j = min (String.length text) (max (i + 1) end_pos.Lexing.pos_cnum) in String.sub text i (j - i) let error_in_file name start_pos end_pos = Printf.sprintf "%s:%d.%d-%d.%d" name start_pos.Lexing.pos_lnum (start_pos.Lexing.pos_cnum - start_pos.Lexing.pos_bol) end_pos.Lexing.pos_lnum (end_pos.Lexing.pos_cnum - end_pos.Lexing.pos_bol) module SMap = Map.Make (struct type t = string let compare = compare end) module Expr = struct type t = UFOx_syntax.expr let of_string text = try UFOx_parser.input UFOx_lexer.token (UFOx_lexer.init_position "" (Lexing.from_string text)) with | UFO_tools.Lexical_Error (msg, start_pos, end_pos) -> invalid_arg (Printf.sprintf "lexical error (%s) at: `%s'" msg (error_in_string text start_pos end_pos)) | UFOx_syntax.Syntax_Error (msg, start_pos, end_pos) -> invalid_arg (Printf.sprintf "syntax error (%s) at: `%s'" msg (error_in_string text start_pos end_pos)) | Parsing.Parse_error -> invalid_arg ("parse error: " ^ text) let of_strings = function | [] -> UFOx_syntax.integer 0 | string :: strings -> List.fold_right (fun s acc -> UFOx_syntax.add (of_string s) acc) strings (of_string string) open UFOx_syntax let rec map f = function | Integer _ | Float _ | Quoted _ as e -> e | Variable s as e -> begin match f s with | Some value -> value | None -> e end | Sum (e1, e2) -> Sum (map f e1, map f e2) | Difference (e1, e2) -> Difference (map f e1, map f e2) | Product (e1, e2) -> Product (map f e1, map f e2) | Quotient (e1, e2) -> Quotient (map f e1, map f e2) | Power (e1, e2) -> Power (map f e1, map f e2) | Application (s, el) -> Application (s, List.map (map f) el) let substitute name value expr = map (fun s -> if s = name then Some value else None) expr let rename1 name_map name = try Some (Variable (SMap.find name name_map)) with Not_found -> None let rename alist_names value = let name_map = List.fold_left (fun acc (name, name') -> SMap.add name name' acc) SMap.empty alist_names in map (rename1 name_map) value let half name = Quotient (Variable name, Integer 2) let variables = UFOx_syntax.variables let functions = UFOx_syntax.functions end module Value = struct module S = UFOx_syntax module Q = Algebra.Q type builtin = | Sqrt | Exp | Log | Log10 | Sin | Asin | Cos | Acos | Tan | Atan | Sinh | Asinh | Cosh | Acosh | Tanh | Atanh | Sec | Asec | Csc | Acsc | Conj | Abs let builtin_to_string = function | Sqrt -> "sqrt" | Exp -> "exp" | Log -> "log" | Log10 -> "log10" | Sin -> "sin" | Cos -> "cos" | Tan -> "tan" | Asin -> "asin" | Acos -> "acos" | Atan -> "atan" | Sinh -> "sinh" | Cosh -> "cosh" | Tanh -> "tanh" | Asinh -> "asinh" | Acosh -> "acosh" | Atanh -> "atanh" | Sec -> "sec" | Csc -> "csc" | Asec -> "asec" | Acsc -> "acsc" | Conj -> "conjg" | Abs -> "abs" let builtin_of_string = function | "cmath.sqrt" -> Sqrt | "cmath.exp" -> Exp | "cmath.log" -> Log | "cmath.log10" -> Log10 | "cmath.sin" -> Sin | "cmath.cos" -> Cos | "cmath.tan" -> Tan | "cmath.asin" -> Asin | "cmath.acos" -> Acos | "cmath.atan" -> Atan | "cmath.sinh" -> Sinh | "cmath.cosh" -> Cosh | "cmath.tanh" -> Tanh | "cmath.asinh" -> Asinh | "cmath.acosh" -> Acosh | "cmath.atanh" -> Atanh | "sec" -> Sec | "csc" -> Csc | "asec" -> Asec | "acsc" -> Acsc | "complexconjugate" -> Conj | "abs" -> Abs | name -> failwith ("UFOx.Value: unsupported function: " ^ name) type t = | Integer of int | Rational of Q.t | Real of float | Complex of float * float | Variable of string | Sum of t list | Difference of t * t | Product of t list | Quotient of t * t | Power of t * t | Application of builtin * t list let rec to_string = function | Integer i -> string_of_int i | Rational q -> Q.to_string q | Real x -> string_of_float x | Complex (0.0, 1.0) -> "I" | Complex (0.0, -1.0) -> "-I" | Complex (0.0, i) -> string_of_float i ^ "*I" | Complex (r, 1.0) -> string_of_float r ^ "+I" | Complex (r, -1.0) -> string_of_float r ^ "-I" | Complex (r, i) -> string_of_float r ^ (if i < 0.0 then "-" else "+") ^ string_of_float (abs_float i) ^ "*I" | Variable s -> s | Sum [] -> "0" | Sum [e] -> to_string e | Sum es -> "(" ^ String.concat "+" (List.map maybe_parentheses es) ^ ")" | Difference (e1, e2) -> to_string e1 ^ "-" ^ maybe_parentheses e2 | Product [] -> "1" | Product ((Integer (-1) | Real (-1.)) :: es) -> "-" ^ maybe_parentheses (Product es) | Product es -> String.concat "*" (List.map maybe_parentheses es) - | Quotient (e1, e2) -> to_string e1 ^ "/" ^ maybe_parentheses e2 + | Quotient (e1, e2) -> maybe_parentheses e1 ^ "/" ^ maybe_parentheses e2 | Power ((Integer i as e), Integer p) -> if p < 0 then maybe_parentheses (Real (float_of_int i)) ^ "^(" ^ string_of_int p ^ ")" else if p = 0 then "1" else if p <= 4 then maybe_parentheses e ^ "^" ^ string_of_int p else maybe_parentheses (Real (float_of_int i)) ^ "^" ^ string_of_int p | Power (e1, e2) -> maybe_parentheses e1 ^ "^" ^ maybe_parentheses e2 | Application (f, [Integer i]) -> to_string (Application (f, [Real (float i)])) | Application (f, es) -> builtin_to_string f ^ "(" ^ String.concat "," (List.map to_string es) ^ ")" and maybe_parentheses = function | Integer i as e -> if i < 0 then "(" ^ to_string e ^ ")" else to_string e | Real x as e -> if x < 0.0 then "(" ^ to_string e ^ ")" else to_string e | Complex (x, 0.0) -> to_string (Real x) | Complex (0.0, 1.0) -> "I" | Variable _ | Power (_, _) | Application (_, _) as e -> to_string e | Sum [e] -> to_string e | Product [e] -> maybe_parentheses e | e -> "(" ^ to_string e ^ ")" let rec to_coupling atom = function | Integer i -> Coupling.Integer i | Rational q -> let n, d = Q.to_ratio q in Coupling.Quot (Coupling.Integer n, Coupling.Integer d) | Real x -> Coupling.Float x | Product es -> Coupling.Prod (List.map (to_coupling atom) es) | Variable s -> Coupling.Atom (atom s) | Complex (r, 0.0) -> Coupling.Float r | Complex (0.0, 1.0) -> Coupling.I | Complex (0.0, -1.0) -> Coupling.Prod [Coupling.I; Coupling.Integer (-1)] | Complex (0.0, i) -> Coupling.Prod [Coupling.I; Coupling.Float i] | Complex (r, 1.0) -> Coupling.Sum [Coupling.Float r; Coupling.I] | Complex (r, -1.0) -> Coupling.Diff (Coupling.Float r, Coupling.I) | Complex (r, i) -> Coupling.Sum [Coupling.Float r; Coupling.Prod [Coupling.I; Coupling.Float i]] | Sum es -> Coupling.Sum (List.map (to_coupling atom) es) | Difference (e1, e2) -> Coupling.Diff (to_coupling atom e1, to_coupling atom e2) | Quotient (e1, e2) -> Coupling.Quot (to_coupling atom e1, to_coupling atom e2) | Power (e1, Integer e2) -> Coupling.Pow (to_coupling atom e1, e2) | Power (e1, e2) -> Coupling.PowX (to_coupling atom e1, to_coupling atom e2) | Application (f, [e]) -> apply1 (to_coupling atom e) f | Application (f, []) -> failwith ("UFOx.Value.to_coupling: " ^ builtin_to_string f ^ ": empty argument list") | Application (f, _::_::_) -> failwith ("UFOx.Value.to_coupling: " ^ builtin_to_string f ^ ": more than one argument in list") and apply1 e = function | Sqrt -> Coupling.Sqrt e | Exp -> Coupling.Exp e | Log -> Coupling.Log e | Log10 -> Coupling.Log10 e | Sin -> Coupling.Sin e | Cos -> Coupling.Cos e | Tan -> Coupling.Tan e | Asin -> Coupling.Asin e | Acos -> Coupling.Acos e | Atan -> Coupling.Atan e | Sinh -> Coupling.Sinh e | Cosh -> Coupling.Cosh e | Tanh -> Coupling.Tanh e | Sec -> Coupling.Quot (Coupling.Integer 1, Coupling.Cos e) | Csc -> Coupling.Quot (Coupling.Integer 1, Coupling.Sin e) | Asec -> Coupling.Acos (Coupling.Quot (Coupling.Integer 1, e)) | Acsc -> Coupling.Asin (Coupling.Quot (Coupling.Integer 1, e)) | Conj -> Coupling.Conj e | Abs -> Coupling.Abs e | (Asinh | Acosh | Atanh as f) -> failwith ("UFOx.Value.to_coupling: function `" ^ builtin_to_string f ^ "' not supported yet!") let compress terms = terms let rec of_expr e = compress (of_expr' e) and of_expr' = function | S.Integer i -> Integer i | S.Float x -> Real x | S.Variable "cmath.pi" -> Variable "pi" | S.Quoted name -> invalid_arg ("UFOx.Value.of_expr: unexpected quoted variable '" ^ name ^ "'") | S.Variable name -> Variable name | S.Sum (e1, e2) -> begin match of_expr e1, of_expr e2 with | (Integer 0 | Real 0.), e -> e | e, (Integer 0 | Real 0.) -> e | Sum e1, Sum e2 -> Sum (e1 @ e2) | e1, Sum e2 -> Sum (e1 :: e2) | Sum e1, e2 -> Sum (e2 :: e1) | e1, e2 -> Sum [e1; e2] end | S.Difference (e1, e2) -> begin match of_expr e1, of_expr e2 with | e1, (Integer 0 | Real 0.) -> e1 | e1, e2 -> Difference (e1, e2) end | S.Product (e1, e2) -> begin match of_expr e1, of_expr e2 with | (Integer 0 | Real 0.), _ -> Integer 0 | _, (Integer 0 | Real 0.) -> Integer 0 | (Integer 1 | Real 1.), e -> e | e, (Integer 1 | Real 1.) -> e | Product e1, Product e2 -> Product (e1 @ e2) | e1, Product e2 -> Product (e1 :: e2) | Product e1, e2 -> Product (e2 :: e1) | e1, e2 -> Product [e1; e2] end | S.Quotient (e1, e2) -> begin match of_expr e1, of_expr e2 with | e1, (Integer 0 | Real 0.) -> invalid_arg "UFOx.Value: divide by 0" | e1, (Integer 1 | Real 1.) -> e1 | e1, e2 -> Quotient (e1, e2) end | S.Power (e, p) -> begin match of_expr e, of_expr p with | (Integer 0 | Real 0.), (Integer 0 | Real 0.) -> invalid_arg "UFOx.Value: 0^0" | _, (Integer 0 | Real 0.) -> Integer 1 | e, (Integer 1 | Real 1.) -> e | Integer e, Integer p -> if p < 0 then Power (Real (float_of_int e), Integer p) else if p = 0 then Integer 1 else if p <= 4 then Power (Integer e, Integer p) else Power (Real (float_of_int e), Integer p) | e, p -> Power (e, p) end | S.Application ("complex", [r; i]) -> begin match of_expr r, of_expr i with | r, (Integer 0 | Real 0.0) -> r | Real r, Real i -> Complex (r, i) | Integer r, Real i -> Complex (float_of_int r, i) | Real r, Integer i -> Complex (r, float_of_int i) | Integer r, Integer i -> Complex (float_of_int r, float_of_int i) | _ -> invalid_arg "UFOx.Value: complex expects two numeric arguments" end | S.Application ("complex", _) -> invalid_arg "UFOx.Value: complex expects two arguments" | S.Application ("complexconjugate", [e]) -> Application (Conj, [of_expr e]) | S.Application ("complexconjugate", _) -> invalid_arg "UFOx.Value: complexconjugate expects single argument" | S.Application ("cmath.sqrt", [e]) -> Application (Sqrt, [of_expr e]) | S.Application ("cmath.sqrt", _) -> invalid_arg "UFOx.Value: sqrt expects single argument" | S.Application (name, args) -> Application (builtin_of_string name, List.map of_expr args) end let positive integers = List.filter (fun (i, _) -> i > 0) integers let not_positive integers = List.filter (fun (i, _) -> i <= 0) integers module type Index = sig type t = int val position : t -> int val factor : t -> int val unpack : t -> int * int val pack : int -> int -> t val map_position : (int -> int) -> t -> t val to_string : t -> string val list_to_string : t list -> string val free : (t * 'r) list -> (t * 'r) list val summation : (t * 'r) list -> (t * 'r) list val classes_to_string : ('r -> string) -> (t * 'r) list -> string val fresh_summation : unit -> t val named_summation : string -> unit -> t end module Index : Index = struct type t = int let free i = positive i let summation i = not_positive i let position i = if i > 0 then i mod 1000 else i let factor i = if i > 0 then i / 1000 else invalid_arg "UFOx.Index.factor: argument not positive" let unpack i = if i > 0 then (position i, factor i) else (i, 0) let pack i j = if j > 0 then if i > 0 then 1000 * j + i else invalid_arg "UFOx.Index.pack: position not positive" else if j = 0 then i else invalid_arg "UFOx.Index.pack: factor negative" let map_position f i = let pos, fac = unpack i in pack (f pos) fac let to_string i = let pos, fac = unpack i in if fac = 0 then Printf.sprintf "%d" pos else Printf.sprintf "%d.%d" pos fac let to_string' = string_of_int let list_to_string is = "[" ^ String.concat ", " (List.map to_string is) ^ "]" let classes_to_string rep_to_string index_classes = let reps = ThoList.uniq (List.sort compare (List.map snd index_classes)) in "[" ^ String.concat ", " (List.map (fun r -> (rep_to_string r) ^ "=" ^ (list_to_string (List.map fst (List.filter (fun (_, r') -> r = r') index_classes)))) reps) ^ "]" type factory = { mutable named : int SMap.t; mutable used : Sets.Int.t } let factory = { named = SMap.empty; used = Sets.Int.empty } let first_anonymous = -1001 let fresh_summation () = let next_anonymous = try pred (Sets.Int.min_elt factory.used) with | Not_found -> first_anonymous in factory.used <- Sets.Int.add next_anonymous factory.used; next_anonymous let named_summation name () = try SMap.find name factory.named with | Not_found -> begin let next_named = fresh_summation () in factory.named <- SMap.add name next_named factory.named; next_named end end module type Atom = sig type t val map_indices : (int -> int) -> t -> t val rename_indices : (int -> int) -> t -> t val contract_pair : t -> t -> t option val variable : t -> string option val scalar : t -> bool val is_unit : t -> bool val invertible : t -> bool val invert : t -> t val of_expr : string -> UFOx_syntax.expr list -> t list val to_string : t -> string type r val classify_indices : t list -> (Index.t * r) list val disambiguate_indices : t list -> t list val rep_to_string : r -> string val rep_to_string_whizard : r -> string val rep_of_int : bool -> int -> r val rep_conjugate : r -> r val rep_trivial : r -> bool type r_omega val omega : r -> r_omega end module type Tensor = sig type atom type 'a linear = ('a list * Algebra.QC.t) list type t = | Linear of atom linear | Ratios of (atom linear * atom linear) list val map_atoms : (atom -> atom) -> t -> t val map_indices : (int -> int) -> t -> t val rename_indices : (int -> int) -> t -> t val map_coeff : (Algebra.QC.t -> Algebra.QC.t) -> t -> t val contract_pairs : t -> t val variables : t -> string list val of_expr : UFOx_syntax.expr -> t val of_string : string -> t val of_strings : string list -> t val to_string : t -> string type r val classify_indices : t -> (Index.t * r) list val rep_to_string : r -> string val rep_to_string_whizard : r -> string val rep_of_int : bool -> int -> r val rep_conjugate : r -> r val rep_trivial : r -> bool type r_omega val omega : r -> r_omega end module Tensor (A : Atom) : Tensor with type atom = A.t and type r = A.r and type r_omega = A.r_omega = struct module S = UFOx_syntax (* TODO: we have to switch to [Algebra.QC] to support complex coefficients, as used in custom propagators. *) module Q = Algebra.Q module QC = Algebra.QC type atom = A.t type 'a linear = ('a list * Algebra.QC.t) list type t = | Linear of atom linear | Ratios of (atom linear * atom linear) list let term_to_string (tensors, c) = if QC.is_null c then "" else match tensors with | [] -> QC.to_string c | tensors -> String.concat "*" ((if QC.is_unit c then [] else [QC.to_string c]) @ List.map A.to_string tensors) let linear_to_string terms = String.concat "" (List.map term_to_string terms) let to_string = function | Linear terms -> linear_to_string terms | Ratios ratios -> String.concat " + " (List.map (fun (n, d) -> Printf.sprintf "(%s)/(%s)" (linear_to_string n) (linear_to_string d)) ratios) let variables_of_atoms atoms = List.fold_left (fun acc a -> match A.variable a with | None -> acc | Some name -> Sets.String.add name acc) Sets.String.empty atoms let variables_of_linear linear = List.fold_left (fun acc (atoms, _) -> Sets.String.union (variables_of_atoms atoms) acc) Sets.String.empty linear let variables_set = function | Linear linear -> variables_of_linear linear | Ratios ratios -> List.fold_left (fun acc (numerator, denominator) -> Sets.String.union (variables_of_linear numerator) (Sets.String.union (variables_of_linear denominator) acc)) Sets.String.empty ratios let variables t = Sets.String.elements (variables_set t) let map_ratios f = function | Linear n -> Linear (f n) | Ratios ratios -> Ratios (List.map (fun (n, d) -> (f n, f d)) ratios) let map_summands f t = map_ratios (List.map f) t let map_numerators f = function | Linear n -> Linear (List.map f n) | Ratios ratios -> Ratios (List.map (fun (n, d) -> (List.map f n, d)) ratios) let map_atoms f t = map_summands (fun (atoms, q) -> (List.map f atoms, q)) t let map_indices f t = map_atoms (A.map_indices f) t let rename_indices f t = map_atoms (A.rename_indices f) t let map_coeff f t = map_numerators (fun (atoms, q) -> (atoms, f q)) t type result = | Matched of atom list | Unmatched of atom list (* [contract_pair a rev_prefix suffix] returns [Unmatched (a :: List.rev_append rev_prefix suffix] if there is no match (as defined by [A.contract_pair]) and [Matched] with the reduced list otherwise. *) let rec contract_pair a rev_prefix = function | [] -> Unmatched (a :: List.rev rev_prefix) | a' :: suffix -> begin match A.contract_pair a a' with | None -> contract_pair a (a' :: rev_prefix) suffix | Some a'' -> if A.is_unit a'' then Matched (List.rev_append rev_prefix suffix) else Matched (List.rev_append rev_prefix (a'' :: suffix)) end (* Use [contract_pair] to find all pairs that match according to [A.contract_pair]. *) let rec contract_pairs1 = function | ([] | [_] as t) -> t | a :: t -> begin match contract_pair a [] t with | Unmatched ([]) -> [] | Unmatched (a' :: t') -> a' :: contract_pairs1 t' | Matched t' -> contract_pairs1 t' end let contract_pairs t = map_summands (fun (t', c) -> (contract_pairs1 t', c)) t let add t1 t2 = match t1, t2 with | Linear l1, Linear l2 -> Linear (l1 @ l2) | Ratios r, Linear l | Linear l, Ratios r -> Ratios ((l, [([], QC.unit)]) :: r) | Ratios r1, Ratios r2 -> Ratios (r1 @ r2) let multiply1 (t1, c1) (t2, c2) = (List.sort compare (t1 @ t2), QC.mul c1 c2) let multiply2 t1 t2 = Product.list2 multiply1 t1 t2 let multiply t1 t2 = match t1, t2 with | Linear l1, Linear l2 -> Linear (multiply2 l1 l2) | Ratios r, Linear l | Linear l, Ratios r -> Ratios (List.map (fun (n, d) -> (multiply2 l n, d)) r) | Ratios r1, Ratios r2 -> Ratios (Product.list2 (fun (n1, d1) (n2, d2) -> (multiply2 n1 n2, multiply2 d1 d2)) r1 r2) let rec power n t = if n < 0 then invalid_arg "UFOx.Tensor.power: n < 0" else if n = 0 then Linear [([], QC.unit)] else if n = 1 then t else multiply t (power (pred n) t) let compress ratios = map_ratios (fun terms -> List.map (fun (t, cs) -> (t, QC.sum cs)) (ThoList.factorize terms)) ratios let rec of_expr e = contract_pairs (compress (of_expr' e)) and of_expr' = function | S.Integer i -> Linear [([], QC.make (Q.make i 1) Q.null)] | S.Float _ -> invalid_arg "UFOx.Tensor.of_expr: unexpected float" | S.Quoted name -> invalid_arg ("UFOx.Tensor.of_expr: unexpected quoted variable '" ^ name ^ "'") | S.Variable name -> (* There should be a gatekeeper here or in [A.of_expr]: *) Linear [(A.of_expr name [], QC.unit)] | S.Application ("complex", [re; im]) -> begin match of_expr re, of_expr im with | Linear [([], re)], Linear [([], im)] -> if QC.is_real re && QC.is_real im then Linear [([], QC.make (QC.real re) (QC.real im))] else invalid_arg ("UFOx.Tensor.of_expr: argument of complex is complex") | _ -> invalid_arg "UFOx.Tensor.of_expr: unexpected argument of complex" end | S.Application (name, args) -> Linear [(A.of_expr name args, QC.unit)] | S.Sum (e1, e2) -> add (of_expr e1) (of_expr e2) | S.Difference (e1, e2) -> add (of_expr e1) (of_expr (S.Product (S.Integer (-1), e2))) | S.Product (e1, e2) -> multiply (of_expr e1) (of_expr e2) | S.Quotient (n, d) -> begin match of_expr n, of_expr d with | n, Linear [] -> invalid_arg "UFOx.Tensor.of_expr: zero denominator" | n, Linear [([], q)] -> map_coeff (fun c -> QC.div c q) n | n, Linear ([(invertibles, q)] as d) -> if List.for_all A.invertible invertibles then let inverses = List.map A.invert invertibles in multiply (Linear [(inverses, QC.inv q)]) n else multiply (Ratios [[([], QC.unit)], d]) n | n, (Linear d as d')-> if List.for_all (fun (t, _) -> List.for_all A.scalar t) d then multiply (Ratios [[([], QC.unit)], d]) n else invalid_arg ("UFOx.Tensor.of_expr: non scalar denominator: " ^ to_string d') | n, (Ratios _ as d) -> invalid_arg ("UFOx.Tensor.of_expr: illegal denominator: " ^ to_string d) end | S.Power (e, p) -> begin match of_expr e, of_expr p with | Linear [([], q)], Linear [([], p)] -> if QC.is_real p then let re_p = QC.real p in if Q.is_integer re_p then Linear [([], QC.pow q (Q.to_integer re_p))] else invalid_arg "UFOx.Tensor.of_expr: rational power of number" else invalid_arg "UFOx.Tensor.of_expr: complex power of number" | Linear [([], q)], _ -> invalid_arg "UFOx.Tensor.of_expr: non-numeric power of number" | t, Linear [([], p)] -> if QC.is_integer p then power (Q.to_integer (QC.real p)) t else invalid_arg "UFOx.Tensor.of_expr: non integer power of tensor" | _ -> invalid_arg "UFOx.Tensor.of_expr: non numeric power of tensor" end type r = A.r let rep_to_string = A.rep_to_string let rep_to_string_whizard = A.rep_to_string_whizard let rep_of_int = A.rep_of_int let rep_conjugate = A.rep_conjugate let rep_trivial = A.rep_trivial let numerators = function | Linear tensors -> tensors | Ratios ratios -> ThoList.flatmap fst ratios let classify_indices' filter tensors = ThoList.uniq (List.sort compare (List.map (fun (t, c) -> filter (A.classify_indices t)) (numerators tensors))) (* NB: the number of summation indices is not guarateed to be the same! Therefore it was foolish to try to check for uniqueness \ldots *) let classify_indices tensors = match classify_indices' Index.free tensors with | [] -> (* There's always at least an empty list! *) failwith "UFOx.Tensor.classify_indices: can't happen!" | [f] -> f | _ -> invalid_arg "UFOx.Tensor.classify_indices: incompatible free indices!" let disambiguate_indices1 (atoms, q) = (A.disambiguate_indices atoms, q) let disambiguate_indices tensors = map_ratios (List.map disambiguate_indices1) tensors let check_indices t = ignore (classify_indices t) let of_expr e = let t = disambiguate_indices (of_expr e) in check_indices t; t let of_string s = of_expr (Expr.of_string s) let of_strings s = of_expr (Expr.of_strings s) type r_omega = A.r_omega let omega = A.omega end module type Lorentz_Atom = sig type dirac = private | C of int * int | Gamma of int * int * int | Gamma5 of int * int | Identity of int * int | ProjP of int * int | ProjM of int * int | Sigma of int * int * int * int type vector = (* private *) | Epsilon of int * int * int * int | Metric of int * int | P of int * int type scalar = (* private *) | Mass of int | Width of int | P2 of int | P12 of int * int | Variable of string | Coeff of Value.t type t = (* private *) | Dirac of dirac | Vector of vector | Scalar of scalar | Inverse of scalar val map_indices_scalar : (int -> int) -> scalar -> scalar val map_indices_vector : (int -> int) -> vector -> vector val rename_indices_vector : (int -> int) -> vector -> vector end module Lorentz_Atom = struct type dirac = | C of int * int | Gamma of int * int * int | Gamma5 of int * int | Identity of int * int | ProjP of int * int | ProjM of int * int | Sigma of int * int * int * int type vector = | Epsilon of int * int * int * int | Metric of int * int | P of int * int type scalar = | Mass of int | Width of int | P2 of int | P12 of int * int | Variable of string | Coeff of Value.t type t = | Dirac of dirac | Vector of vector | Scalar of scalar | Inverse of scalar let map_indices_scalar f = function | Mass i -> Mass (f i) | Width i -> Width (f i) | P2 i -> P2 (f i) | P12 (i, j) -> P12 (f i, f j) | (Variable _ | Coeff _ as s) -> s let map_indices_vector f = function | Epsilon (mu, nu, ka, la) -> Epsilon (f mu, f nu, f ka, f la) | Metric (mu, nu) -> Metric (f mu, f nu) | P (mu, n) -> P (f mu, f n) let rename_indices_vector f = function | Epsilon (mu, nu, ka, la) -> Epsilon (f mu, f nu, f ka, f la) | Metric (mu, nu) -> Metric (f mu, f nu) | P (mu, n) -> P (f mu, n) end module Lorentz_Atom' : Atom with type t = Lorentz_Atom.t and type r_omega = Coupling.lorentz = struct type t = Lorentz_Atom.t open Lorentz_Atom let map_indices_dirac f = function | C (i, j) -> C (f i, f j) | Gamma (mu, i, j) -> Gamma (f mu, f i, f j) | Gamma5 (i, j) -> Gamma5 (f i, f j) | Identity (i, j) -> Identity (f i, f j) | ProjP (i, j) -> ProjP (f i, f j) | ProjM (i, j) -> ProjM (f i, f j) | Sigma (mu, nu, i, j) -> Sigma (f mu, f nu, f i, f j) let rename_indices_dirac = map_indices_dirac let map_indices_scalar f = function | Mass i -> Mass (f i) | Width i -> Width (f i) | P2 i -> P2 (f i) | P12 (i, j) -> P12 (f i, f j) | Variable s -> Variable s | Coeff c -> Coeff c let map_indices f = function | Dirac d -> Dirac (map_indices_dirac f d) | Vector v -> Vector (map_indices_vector f v) | Scalar s -> Scalar (map_indices_scalar f s) | Inverse s -> Inverse (map_indices_scalar f s) let rename_indices2 fd fv = function | Dirac d -> Dirac (rename_indices_dirac fd d) | Vector v -> Vector (rename_indices_vector fv v) | Scalar s -> Scalar s | Inverse s -> Inverse s let rename_indices f atom = rename_indices2 f f atom let contract_pair a1 a2 = match a1, a2 with | Vector (P (mu1, i1)), Vector (P (mu2, i2)) -> if mu1 <= 0 && mu1 = mu2 then if i1 = i2 then Some (Scalar (P2 i1)) else Some (Scalar (P12 (i1, i2))) else None | Scalar s, Inverse s' | Inverse s, Scalar s' -> if s = s' then Some (Scalar (Coeff (Value.Integer 1))) else None | _ -> None let variable = function | Scalar (Variable s) | Inverse (Variable s) -> Some s | _ -> None let scalar = function | Dirac _ | Vector _ -> false | Scalar _ | Inverse _ -> true let is_unit = function | Scalar (Coeff c) | Inverse (Coeff c) -> begin match c with | Value.Integer 1 -> true | Value.Rational q -> Algebra.Q.is_unit q | _ -> false end | _ -> false let invertible = scalar let invert = function | Dirac _ -> invalid_arg "UFOx.Lorentz_Atom.invert Dirac" | Vector _ -> invalid_arg "UFOx.Lorentz_Atom.invert Vector" | Scalar s -> Inverse s | Inverse s -> Scalar s let i2s = Index.to_string let dirac_to_string = function | C (i, j) -> Printf.sprintf "C(%s,%s)" (i2s i) (i2s j) | Gamma (mu, i, j) -> Printf.sprintf "Gamma(%s,%s,%s)" (i2s mu) (i2s i) (i2s j) | Gamma5 (i, j) -> Printf.sprintf "Gamma5(%s,%s)" (i2s i) (i2s j) | Identity (i, j) -> Printf.sprintf "Identity(%s,%s)" (i2s i) (i2s j) | ProjP (i, j) -> Printf.sprintf "ProjP(%s,%s)" (i2s i) (i2s j) | ProjM (i, j) -> Printf.sprintf "ProjM(%s,%s)" (i2s i) (i2s j) | Sigma (mu, nu, i, j) -> Printf.sprintf "Sigma(%s,%s,%s,%s)" (i2s mu) (i2s nu) (i2s i) (i2s j) let vector_to_string = function | Epsilon (mu, nu, ka, la) -> Printf.sprintf "Epsilon(%s,%s,%s,%s)" (i2s mu) (i2s nu) (i2s ka) (i2s la) | Metric (mu, nu) -> Printf.sprintf "Metric(%s,%s)" (i2s mu) (i2s nu) | P (mu, n) -> Printf.sprintf "P(%s,%d)" (i2s mu) n let scalar_to_string = function | Mass id -> Printf.sprintf "Mass(%d)" id | Width id -> Printf.sprintf "Width(%d)" id | P2 id -> Printf.sprintf "P(%d)**2" id | P12 (id1, id2) -> Printf.sprintf "P(%d)*P(%d)" id1 id2 | Variable s -> s | Coeff c -> Value.to_string c let to_string = function | Dirac d -> dirac_to_string d | Vector v -> vector_to_string v | Scalar s -> scalar_to_string s | Inverse s -> "1/" ^ scalar_to_string s module S = UFOx_syntax (* \begin{dubious} Here we handle some special cases in order to be able to parse propagators. This needs to be made more general, but unfortunately the syntax for the propagator extension is not well documented and appears to be a bit chaotic! \end{dubious} *) let quoted_index s = Index.named_summation s () let integer_or_id = function | S.Integer n -> n | S.Variable "id" -> 1 | _ -> failwith "UFOx.Lorentz_Atom.integer_or_id: impossible" let vector_index = function | S.Integer n -> n | S.Quoted mu -> quoted_index mu | S.Variable id -> let l = String.length id in if l > 1 then if id.[0] = 'l' then int_of_string (String.sub id 1 (pred l)) else invalid_arg ("UFOx.Lorentz_Atom.vector_index: " ^ id) else invalid_arg "UFOx.Lorentz_Atom.vector_index: empty variable" | _ -> invalid_arg "UFOx.Lorentz_Atom.vector_index" let spinor_index = function | S.Integer n -> n | S.Variable id -> let l = String.length id in if l > 1 then if id.[0] = 's' then int_of_string (String.sub id 1 (pred l)) else invalid_arg ("UFOx.Lorentz_Atom.spinor_index: " ^ id) else invalid_arg "UFOx.Lorentz_Atom.spinor_index: empty variable" | _ -> invalid_arg "UFOx.Lorentz_Atom.spinor_index" let of_expr name args = match name, args with | "C", [i; j] -> [Dirac (C (spinor_index i, spinor_index j))] | "C", _ -> invalid_arg "UFOx.Lorentz.of_expr: invalid arguments to C()" | "Epsilon", [mu; nu; ka; la] -> [Vector (Epsilon (vector_index mu, vector_index nu, vector_index ka, vector_index la))] | "Epsilon", _ -> invalid_arg "UFOx.Lorentz.of_expr: invalid arguments to Epsilon()" | "Gamma", [mu; i; j] -> [Dirac (Gamma (vector_index mu, spinor_index i, spinor_index j))] | "Gamma", _ -> invalid_arg "UFOx.Lorentz.of_expr: invalid arguments to Gamma()" | "Gamma5", [i; j] -> [Dirac (Gamma5 (spinor_index i, spinor_index j))] | "Gamma5", _ -> invalid_arg "UFOx.Lorentz.of_expr: invalid arguments to Gamma5()" | "Identity", [i; j] -> [Dirac (Identity (spinor_index i, spinor_index j))] | "Identity", _ -> invalid_arg "UFOx.Lorentz.of_expr: invalid arguments to Identity()" | "Metric", [mu; nu] -> [Vector (Metric (vector_index mu, vector_index nu))] | "Metric", _ -> invalid_arg "UFOx.Lorentz.of_expr: invalid arguments to Metric()" | "P", [mu; id] -> [Vector (P (vector_index mu, integer_or_id id))] | "P", _ -> invalid_arg "UFOx.Lorentz.of_expr: invalid arguments to P()" | "ProjP", [i; j] -> [Dirac (ProjP (spinor_index i, spinor_index j))] | "ProjP", _ -> invalid_arg "UFOx.Lorentz.of_expr: invalid arguments to ProjP()" | "ProjM", [i; j] -> [Dirac (ProjM (spinor_index i, spinor_index j))] | "ProjM", _ -> invalid_arg "UFOx.Lorentz.of_expr: invalid arguments to ProjM()" | "Sigma", [mu; nu; i; j] -> if mu <> nu then [Dirac (Sigma (vector_index mu, vector_index nu, spinor_index i, spinor_index j))] else invalid_arg "UFOx.Lorentz.of_expr: implausible arguments to Sigma()" | "Sigma", _ -> invalid_arg "UFOx.Lorentz.of_expr: invalid arguments to Sigma()" | "PSlash", [i; j; id] -> let mu = Index.fresh_summation () in [Dirac (Gamma (mu, spinor_index i, spinor_index j)); Vector (P (mu, integer_or_id id))] | "PSlash", _ -> invalid_arg "UFOx.Lorentz.of_expr: invalid arguments to PSlash()" | "Mass", [id] -> [Scalar (Mass (integer_or_id id))] | "Mass", _ -> invalid_arg "UFOx.Lorentz.of_expr: invalid arguments to Mass()" | "Width", [id] -> [Scalar (Width (integer_or_id id))] | "Width", _ -> invalid_arg "UFOx.Lorentz.of_expr: invalid arguments to Width()" | name, [] -> [Scalar (Variable name)] | name, _ -> invalid_arg ("UFOx.Lorentz.of_expr: invalid tensor '" ^ name ^ "'") type r = S | V | T | Sp | CSp | Maj | VSp | CVSp | VMaj | Ghost let rep_trivial = function | S | Ghost -> true | V | T | Sp | CSp | Maj | VSp | CVSp | VMaj -> false let rep_to_string = function | S -> "0" | V -> "1" | T -> "2" | Sp -> "1/2" | CSp-> "1/2bar" | Maj -> "1/2M" | VSp -> "3/2" | CVSp -> "3/2bar" | VMaj -> "3/2M" | Ghost -> "Ghost" let rep_to_string_whizard = function | S -> "0" | V -> "1" | T -> "2" | Sp | CSp | Maj -> "1/2" | VSp | CVSp | VMaj -> "3/2" | Ghost -> "Ghost" let rep_of_int neutral = function | -1 -> Ghost | 1 -> S | 2 -> if neutral then Maj else Sp | -2 -> if neutral then Maj else CSp (* used by [UFO.Particle.force_conjspinor] *) | 3 -> V | 4 -> if neutral then VMaj else VSp | -4 -> if neutral then VMaj else CVSp (* used by [UFO.Particle.force_conjspinor] *) | 5 -> T | s when s > 0 -> failwith "UFOx.Lorentz: spin > 2 not supported!" | _ -> invalid_arg "UFOx.Lorentz: invalid non-positive spin value" let rep_conjugate = function | S -> S | V -> V | T -> T | Sp -> CSp (* ??? *) | CSp -> Sp (* ??? *) | Maj -> Maj | VSp -> CVSp | CVSp -> VSp | VMaj -> VMaj | Ghost -> Ghost let classify_vector_indices1 = function | Epsilon (mu, nu, ka, la) -> [(mu, V); (nu, V); (ka, V); (la, V)] | Metric (mu, nu) -> [(mu, V); (nu, V)] | P (mu, n) -> [(mu, V)] let classify_dirac_indices1 = function | C (i, j) -> [(i, CSp); (j, Sp)] (* ??? *) | Gamma5 (i, j) | Identity (i, j) | ProjP (i, j) | ProjM (i, j) -> [(i, CSp); (j, Sp)] | Gamma (mu, i, j) -> [(mu, V); (i, CSp); (j, Sp)] | Sigma (mu, nu, i, j) -> [(mu, V); (nu, V); (i, CSp); (j, Sp)] let classify_indices1 = function | Dirac d -> classify_dirac_indices1 d | Vector v -> classify_vector_indices1 v | Scalar _ | Inverse _ -> [] module IMap = Map.Make (struct type t = int let compare = compare end) exception Incompatible_factors of r * r let product rep1 rep2 = match rep1, rep2 with | V, V -> T | V, Sp -> VSp | V, CSp -> CVSp | V, Maj -> VMaj | Sp, V -> VSp | CSp, V -> CVSp | Maj, V -> VMaj | _, _ -> raise (Incompatible_factors (rep1, rep2)) let combine_or_add_index (i, rep) map = let pos, fac = Index.unpack i in try let fac', rep' = IMap.find pos map in if pos < 0 then IMap.add pos (fac, rep) map else if fac <> fac' then IMap.add pos (0, product rep rep') map else if rep <> rep' then (* Can be disambiguated! *) IMap.add pos (0, product rep rep') map else invalid_arg (Printf.sprintf "UFO: duplicate subindex %d" pos) with | Not_found -> IMap.add pos (fac, rep) map | Incompatible_factors (rep1, rep2) -> invalid_arg (Printf.sprintf "UFO: incompatible factors (%s,%s) at %d" (rep_to_string rep1) (rep_to_string rep2) pos) let combine_or_add_indices atom map = List.fold_right combine_or_add_index (classify_indices1 atom) map let project_factors (pos, (fac, rep)) = if fac = 0 then (pos, rep) else invalid_arg (Printf.sprintf "UFO: leftover subindex %d.%d" pos fac) let classify_indices atoms = List.map project_factors (IMap.bindings (List.fold_right combine_or_add_indices atoms IMap.empty)) let add_factor fac indices pos = if pos > 0 then if Sets.Int.mem pos indices then Index.pack pos fac else pos else pos let disambiguate_indices1 indices atom = rename_indices2 (add_factor 1 indices) (add_factor 2 indices) atom let vectorspinors atoms = List.fold_left (fun acc (i, r) -> match r with | S | V | T | Sp | CSp | Maj | Ghost -> acc | VSp | CVSp | VMaj -> Sets.Int.add i acc) Sets.Int.empty (classify_indices atoms) let disambiguate_indices atoms = let vectorspinor_indices = vectorspinors atoms in List.map (disambiguate_indices1 vectorspinor_indices) atoms type r_omega = Coupling.lorentz let omega = function | S -> Coupling.Scalar | V -> Coupling.Vector | T -> Coupling.Tensor_2 | Sp -> Coupling.Spinor | CSp -> Coupling.ConjSpinor | Maj -> Coupling.Majorana | VSp -> Coupling.Vectorspinor | CVSp -> Coupling.Vectorspinor (* TODO: not really! *) | VMaj -> Coupling.Vectorspinor (* TODO: not really! *) | Ghost -> Coupling.Scalar end module Lorentz = Tensor(Lorentz_Atom') module type Color_Atom = sig type t = (* private *) | Identity of int * int | Identity8 of int * int | T of int * int * int | F of int * int * int | D of int * int * int | Epsilon of int * int * int | EpsilonBar of int * int * int | T6 of int * int * int | K6 of int * int * int | K6Bar of int * int * int end module Color_Atom = struct type t = | Identity of int * int | Identity8 of int * int | T of int * int * int | F of int * int * int | D of int * int * int | Epsilon of int * int * int | EpsilonBar of int * int * int | T6 of int * int * int | K6 of int * int * int | K6Bar of int * int * int end module Color_Atom' : Atom with type t = Color_Atom.t and type r_omega = Color.t = struct type t = Color_Atom.t module S = UFOx_syntax open Color_Atom let map_indices f = function | Identity (i, j) -> Identity (f i, f j) | Identity8 (a, b) -> Identity8 (f a, f b) | T (a, i, j) -> T (f a, f i, f j) | F (a, i, j) -> F (f a, f i, f j) | D (a, i, j) -> D (f a, f i, f j) | Epsilon (i, j, k) -> Epsilon (f i, f j, f k) | EpsilonBar (i, j, k) -> EpsilonBar (f i, f j, f k) | T6 (a, i', j') -> T6 (f a, f i', f j') | K6 (i', j, k) -> K6 (f i', f j, f k) | K6Bar (i', j, k) -> K6Bar (f i', f j, f k) let rename_indices = map_indices let contract_pair _ _ = None let variable _ = None let scalar _ = false let invertible _ = false let is_unit _ = false let invert _ = invalid_arg "UFOx.Color_Atom.invert" let of_expr1 name args = match name, args with | "Identity", [S.Integer i; S.Integer j] -> Identity (i, j) | "Identity", _ -> invalid_arg "UFOx.Color.of_expr: invalid arguments to Identity()" | "T", [S.Integer a; S.Integer i; S.Integer j] -> T (a, i, j) | "T", _ -> invalid_arg "UFOx.Color.of_expr: invalid arguments to T()" | "f", [S.Integer a; S.Integer b; S.Integer c] -> F (a, b, c) | "f", _ -> invalid_arg "UFOx.Color.of_expr: invalid arguments to f()" | "d", [S.Integer a; S.Integer b; S.Integer c] -> D (a, b, c) | "d", _ -> invalid_arg "UFOx.Color.of_expr: invalid arguments to d()" | "Epsilon", [S.Integer i; S.Integer j; S.Integer k] -> Epsilon (i, j, k) | "Epsilon", _ -> invalid_arg "UFOx.Color.of_expr: invalid arguments to Epsilon()" | "EpsilonBar", [S.Integer i; S.Integer j; S.Integer k] -> EpsilonBar (i, j, k) | "EpsilonBar", _ -> invalid_arg "UFOx.Color.of_expr: invalid arguments to EpsilonBar()" | "T6", [S.Integer a; S.Integer i'; S.Integer j'] -> T6 (a, i', j') | "T6", _ -> invalid_arg "UFOx.Color.of_expr: invalid arguments to T6()" | "K6", [S.Integer i'; S.Integer j; S.Integer k] -> K6 (i', j, k) | "K6", _ -> invalid_arg "UFOx.Color.of_expr: invalid arguments to K6()" | "K6Bar", [S.Integer i'; S.Integer j; S.Integer k] -> K6Bar (i', j, k) | "K6Bar", _ -> invalid_arg "UFOx.Color.of_expr: invalid arguments to K6Bar()" | name, _ -> invalid_arg ("UFOx.Color.of_expr: invalid tensor '" ^ name ^ "'") let of_expr name args = [of_expr1 name args] let to_string = function | Identity (i, j) -> Printf.sprintf "Identity(%d,%d)" i j | Identity8 (a, b) -> Printf.sprintf "Identity8(%d,%d)" a b | T (a, i, j) -> Printf.sprintf "T(%d,%d,%d)" a i j | F (a, b, c) -> Printf.sprintf "f(%d,%d,%d)" a b c | D (a, b, c) -> Printf.sprintf "d(%d,%d,%d)" a b c | Epsilon (i, j, k) -> Printf.sprintf "Epsilon(%d,%d,%d)" i j k | EpsilonBar (i, j, k) -> Printf.sprintf "EpsilonBar(%d,%d,%d)" i j k | T6 (a, i', j') -> Printf.sprintf "T6(%d,%d,%d)" a i' j' | K6 (i', j, k) -> Printf.sprintf "K6(%d,%d,%d)" i' j k | K6Bar (i', j, k) -> Printf.sprintf "K6Bar(%d,%d,%d)" i' j k type r = S | F | C | A let rep_trivial = function | S -> true | F | C | A -> false let rep_to_string = function | S -> "1" | F -> "3" | C -> "3bar" | A-> "8" let rep_to_string_whizard = function | S -> "1" | F -> "3" | C -> "-3" | A-> "8" let rep_of_int neutral = function | 1 -> S | 3 -> F | -3 -> C | 8 -> A | 6 | -6 -> failwith "UFOx.Color: sextets not supported yet!" | 10 | -10 -> failwith "UFOx.Color: decuplets not supported yet!" | n -> invalid_arg (Printf.sprintf "UFOx.Color: impossible representation color = %d!" n) let rep_conjugate = function | S -> S | C -> F | F -> C | A -> A let classify_indices1 = function | Identity (i, j) -> [(i, C); (j, F)] | Identity8 (a, b) -> [(a, A); (b, A)] | T (a, i, j) -> [(i, F); (j, C); (a, A)] | Color_Atom.F (a, b, c) | D (a, b, c) -> [(a, A); (b, A); (c, A)] | Epsilon (i, j, k) -> [(i, F); (j, F); (k, F)] | EpsilonBar (i, j, k) -> [(i, C); (j, C); (k, C)] | T6 (a, i', j') -> failwith "UFOx.Color: sextets not supported yet!" | K6 (i', j, k) -> failwith "UFOx.Color: sextets not supported yet!" | K6Bar (i', j, k) -> failwith "UFOx.Color: sextets not supported yet!" let classify_indices tensors = List.sort compare (List.fold_right (fun v acc -> classify_indices1 v @ acc) tensors []) let disambiguate_indices atoms = atoms type r_omega = Color.t (* FIXME: $N_C=3$ should not be hardcoded! *) let omega = function | S -> Color.Singlet | F -> Color.SUN (3) | C -> Color.SUN (-3) | A -> Color.AdjSUN (3) end module Color = Tensor(Color_Atom') module type Test = sig - val example : unit -> unit val suite : OUnit.test end + +module Test : Test = + struct + + open OUnit + + let parse_unparse s = + Value.to_string (Value.of_expr (Expr.of_string s)) + + let assert_parse_unparse unparsed expr = + assert_equal ~printer:(fun s -> s) unparsed (parse_unparse expr) + + let suite_expr = + "unparse/parse" >::: + [ "a + b" >:: + (fun () -> assert_parse_unparse "(a+b)" "a+b"); + + "(a - b) / c" >:: + (fun () -> assert_parse_unparse "(a-b)/c" "(a-b)/c"); + + "(a + b - c) / d" >:: + (fun () -> assert_parse_unparse "((a+b)-c)/d" "(a+b-c)/d"); + + "S2HDMIV:lam1" >:: + (fun () -> + assert_parse_unparse + "(((Mh3^2*RA3x1^2)+(Mh1^2*RA1x1^2)+(Mh2^2*RA2x1^2))-(musq*SB^2))/(CB^2*vH^2)" + "(Mh1**2*RA1x1**2 + Mh2**2*RA2x1**2 + Mh3**2*RA3x1**2 - musq*SB**2)/(CB**2*vH**2)") ] + + let suite = + "UFOx" >::: + [suite_expr] + + end +