vm.fail
symbolic mathematics engine in OCaml with differentiation, integration, simplification, and numerical methods
1open Expr
2open Canonical
3
4type cse_result = {
5 final_expr: expr;
6 subexpressions: (string * expr) list;
7}
8
9let common_subexpression_elimination expr =
10 let counter = ref 0 in
11 let subexpr_map = Hashtbl.create 100 in
12 let name_map = Hashtbl.create 100 in
13
14 let get_or_create_name e =
15 let h = hash e in
16 match Hashtbl.find_opt name_map h with
17 | Some name -> Var name
18 | None ->
19 incr counter;
20 let name = Printf.sprintf "cse_%d" !counter in
21 Hashtbl.add name_map h name;
22 Hashtbl.add subexpr_map name e;
23 Var name
24 in
25
26 let should_extract = function
27 | Const _ | SymConst _ | Var _ -> false
28 | _ -> true
29 in
30
31 let rec count_occurrences e counts =
32 let h = hash e in
33 Hashtbl.replace counts h (1 + (Hashtbl.find_opt counts h |> Option.value ~default:0));
34 match e with
35 | Add (e1, e2) | Sub (e1, e2) | Mul (e1, e2) | Div (e1, e2) | Pow (e1, e2) ->
36 count_occurrences e1 counts;
37 count_occurrences e2 counts
38 | Neg e | Sin e | Cos e | Tan e | Sinh e | Cosh e | Tanh e
39 | Asin e | Acos e | Atan e | Exp e | Ln e | Sqrt e | Abs e ->
40 count_occurrences e counts
41 | Atan2 (e1, e2) | Log (e1, e2) ->
42 count_occurrences e1 counts;
43 count_occurrences e2 counts
44 | _ -> ()
45 in
46
47 let counts = Hashtbl.create 100 in
48 count_occurrences expr counts;
49
50 let rec extract e =
51 if should_extract e && Hashtbl.find counts (hash e) >= 2 then
52 get_or_create_name e
53 else
54 match e with
55 | Add (e1, e2) -> Add (extract e1, extract e2)
56 | Sub (e1, e2) -> Sub (extract e1, extract e2)
57 | Mul (e1, e2) -> Mul (extract e1, extract e2)
58 | Div (e1, e2) -> Div (extract e1, extract e2)
59 | Pow (e1, e2) -> Pow (extract e1, extract e2)
60 | Neg e -> Neg (extract e)
61 | Sin e -> Sin (extract e)
62 | Cos e -> Cos (extract e)
63 | Tan e -> Tan (extract e)
64 | Sinh e -> Sinh (extract e)
65 | Cosh e -> Cosh (extract e)
66 | Tanh e -> Tanh (extract e)
67 | Asin e -> Asin (extract e)
68 | Acos e -> Acos (extract e)
69 | Atan e -> Atan (extract e)
70 | Atan2 (e1, e2) -> Atan2 (extract e1, extract e2)
71 | Exp e -> Exp (extract e)
72 | Ln e -> Ln (extract e)
73 | Log (e1, e2) -> Log (extract e1, extract e2)
74 | Sqrt e -> Sqrt (extract e)
75 | Abs e -> Abs (extract e)
76 | e -> e
77 in
78
79 let final = extract expr in
80 let subexprs = Hashtbl.fold (fun name e acc -> (name, e) :: acc) subexpr_map [] in
81 {final_expr = final; subexpressions = List.sort (fun (n1, _) (n2, _) -> String.compare n1 n2) subexprs}
82
83let horner_form expr var =
84 let rec collect_poly e =
85 match e with
86 | Const c -> [(0, Const c)]
87 | Var v when v = var -> [(1, Const 1.0)]
88 | Pow (Var v, Const n) when v = var && Float.is_integer n ->
89 [(int_of_float n, Const 1.0)]
90 | Mul (Const c, Pow (Var v, Const n)) when v = var && Float.is_integer n ->
91 [(int_of_float n, Const c)]
92 | Add (e1, e2) -> collect_poly e1 @ collect_poly e2
93 | _ -> [(0, e)]
94 in
95
96 let terms = collect_poly expr in
97 let max_degree = List.fold_left (fun acc (deg, _) -> max acc deg) 0 terms in
98
99 let coeffs = Array.make (max_degree + 1) (Const 0.0) in
100 List.iter (fun (deg, coeff) ->
101 coeffs.(deg) <- Simplify.simplify (Add (coeffs.(deg), coeff))
102 ) terms;
103
104 let rec build_horner deg =
105 if deg < 0 then Const 0.0
106 else if deg = 0 then coeffs.(0)
107 else Add (coeffs.(deg), Mul (Var var, build_horner (deg - 1)))
108 in
109
110 Simplify.simplify (build_horner max_degree)