vm.fail
symbolic mathematics engine in OCaml with differentiation, integration, simplification, and numerical methods
1open Expr
2
3let max_iterations = 100
4
5let rec simplify_once = function
6 | Const _ as c -> c
7 | SymConst _ as s -> s
8 | Var _ as v -> v
9 | Add (e1, e2) -> simplify_add (simplify_once e1) (simplify_once e2)
10 | Sub (e1, e2) -> simplify_sub (simplify_once e1) (simplify_once e2)
11 | Mul (e1, e2) -> simplify_mul (simplify_once e1) (simplify_once e2)
12 | Div (e1, e2) -> simplify_div (simplify_once e1) (simplify_once e2)
13 | Pow (e1, e2) -> simplify_pow (simplify_once e1) (simplify_once e2)
14 | Neg e -> simplify_neg (simplify_once e)
15 | Sin e -> simplify_sin (simplify_once e)
16 | Cos e -> simplify_cos (simplify_once e)
17 | Tan e -> simplify_tan (simplify_once e)
18 | Sinh e -> simplify_sinh (simplify_once e)
19 | Cosh e -> simplify_cosh (simplify_once e)
20 | Tanh e -> simplify_tanh (simplify_once e)
21 | Asin e -> simplify_asin (simplify_once e)
22 | Acos e -> simplify_acos (simplify_once e)
23 | Atan e -> simplify_atan (simplify_once e)
24 | Atan2 (e1, e2) -> simplify_atan2 (simplify_once e1) (simplify_once e2)
25 | Exp e -> simplify_exp (simplify_once e)
26 | Ln e -> simplify_ln (simplify_once e)
27 | Log (e1, e2) -> simplify_log (simplify_once e1) (simplify_once e2)
28 | Sqrt e -> simplify_sqrt (simplify_once e)
29 | Abs e -> simplify_abs (simplify_once e)
30
31and simplify_add e1 e2 =
32 match (e1, e2) with
33 | Const 0.0, e | e, Const 0.0 -> e
34 | Const a, Const b -> Const (a +. b)
35 | Mul (Const a, e1), Mul (Const b, e2) when Canonical.equal e1 e2 ->
36 simplify_mul (Const (a +. b)) e1
37 | Mul (Const a, e1), e2 when Canonical.equal e1 e2 ->
38 simplify_mul (Const (a +. 1.0)) e1
39 | e1, Mul (Const b, e2) when Canonical.equal e1 e2 ->
40 simplify_mul (Const (1.0 +. b)) e1
41 | e1, e2 when Canonical.equal e1 e2 ->
42 simplify_mul (Const 2.0) e1
43 | Add (e1, Const a), Const b -> simplify_add e1 (Const (a +. b))
44 | _ -> Add (e1, e2)
45
46and simplify_sub e1 e2 =
47 match (e1, e2) with
48 | e, Const 0.0 -> e
49 | Const a, Const b -> Const (a -. b)
50 | e1, e2 when Canonical.equal e1 e2 -> Const 0.0
51 | _ -> Sub (e1, e2)
52
53and simplify_mul e1 e2 =
54 match (e1, e2) with
55 | Const 0.0, _ | _, Const 0.0 -> Const 0.0
56 | Const 1.0, e | e, Const 1.0 -> e
57 | Const (-1.0), e -> simplify_neg e
58 | e, Const (-1.0) -> simplify_neg e
59 | Const a, Const b -> Const (a *. b)
60 | Const a, Mul (Const b, e) -> simplify_mul (Const (a *. b)) e
61 | Mul (Const a, e), Const b -> simplify_mul (Const (a *. b)) e
62 | Const a, Mul (e1, Mul (Const b, e2)) ->
63 simplify_mul (Const (a *. b)) (Mul (e1, e2))
64 | Mul (Const a, e1), Mul (Const b, e2) ->
65 simplify_mul (Const (a *. b)) (Mul (e1, e2))
66 | Pow (e1, a), Pow (e2, b) when Canonical.equal e1 e2 ->
67 simplify_pow e1 (simplify_add a b)
68 | e1, Pow (e2, b) when Canonical.equal e1 e2 ->
69 simplify_pow e1 (simplify_add (Const 1.0) b)
70 | Pow (e1, a), e2 when Canonical.equal e1 e2 ->
71 simplify_pow e1 (simplify_add a (Const 1.0))
72 | e1, e2 when Canonical.equal e1 e2 ->
73 simplify_pow e1 (Const 2.0)
74 | Exp e1, Exp e2 -> Exp (simplify_add e1 e2)
75 | (Sin _ | Cos _ | Tan _ | Sinh _ | Cosh _ | Tanh _ |
76 Asin _ | Acos _ | Atan _ | Exp _ | Ln _ | Log _ | Sqrt _ | Abs _ | Pow _), Var _ ->
77 Mul (e2, e1)
78 | _ -> Mul (e1, e2)
79
80and simplify_div e1 e2 =
81 match (e1, e2) with
82 | Const 0.0, _ -> Const 0.0
83 | e, Const 1.0 -> e
84 | Const a, Const b -> Const (a /. b)
85 | e1, e2 when Canonical.equal e1 e2 -> Const 1.0
86 | Mul (e1, e2), e3 when Canonical.equal e2 e3 -> e1
87 | Mul (e1, e2), e3 when Canonical.equal e1 e3 -> e2
88 | _ -> Div (e1, e2)
89
90and simplify_pow e1 e2 =
91 match (e1, e2) with
92 | _, Const 0.0 -> Const 1.0
93 | e, Const 1.0 -> e
94 | Const 0.0, _ -> Const 0.0
95 | Const 1.0, _ -> Const 1.0
96 | Const a, Const b -> Const (a ** b)
97 | Pow (e, a), b -> simplify_pow e (simplify_mul a b)
98 | Sqrt e, Const 2.0 -> e
99 | e, Const 0.5 -> simplify_sqrt e
100 | SymConst E, Ln e -> e
101 | _ -> Pow (e1, e2)
102
103and simplify_neg = function
104 | Const c -> Const (-.c)
105 | Neg e -> e
106 | Mul (Const c, e) -> simplify_mul (Const (-.c)) e
107 | Mul (e, Const c) -> simplify_mul (Const (-.c)) e
108 | e -> Neg e
109
110and simplify_sin = function
111 | Const 0.0 -> Const 0.0
112 | Const c -> Const (sin c)
113 | Asin e -> e
114 | Neg e -> simplify_neg (Sin e)
115 | e -> Sin e
116
117and simplify_cos = function
118 | Const 0.0 -> Const 1.0
119 | Const c -> Const (cos c)
120 | Acos e -> e
121 | Neg e -> Cos e
122 | e -> Cos e
123
124and simplify_tan = function
125 | Const 0.0 -> Const 0.0
126 | Const c -> Const (tan c)
127 | Atan e -> e
128 | Neg e -> simplify_neg (Tan e)
129 | e -> Tan e
130
131and simplify_sinh = function
132 | Const 0.0 -> Const 0.0
133 | Const c -> Const (sinh c)
134 | Neg e -> simplify_neg (Sinh e)
135 | e -> Sinh e
136
137and simplify_cosh = function
138 | Const 0.0 -> Const 1.0
139 | Const c -> Const (cosh c)
140 | Neg e -> Cosh e
141 | e -> Cosh e
142
143and simplify_tanh = function
144 | Const 0.0 -> Const 0.0
145 | Const c -> Const (tanh c)
146 | Neg e -> simplify_neg (Tanh e)
147 | e -> Tanh e
148
149and simplify_asin = function
150 | Const c -> Const (asin c)
151 | Sin e -> e
152 | e -> Asin e
153
154and simplify_acos = function
155 | Const c -> Const (acos c)
156 | Cos e -> e
157 | e -> Acos e
158
159and simplify_atan = function
160 | Const c -> Const (atan c)
161 | Tan e -> e
162 | e -> Atan e
163
164and simplify_atan2 e1 e2 =
165 match (e1, e2) with
166 | Const a, Const b -> Const (atan2 a b)
167 | _ -> Atan2 (e1, e2)
168
169and simplify_exp = function
170 | Const 0.0 -> Const 1.0
171 | Const c -> Const (exp c)
172 | Ln e -> e
173 | Add (e1, e2) -> simplify_mul (Exp e1) (Exp e2)
174 | Mul (Const c, e) -> simplify_pow (Exp e) (Const c)
175 | e -> Exp e
176
177and simplify_ln = function
178 | Const 1.0 -> Const 0.0
179 | Const c when c > 0.0 -> Const (log c)
180 | Exp e -> e
181 | SymConst E -> Const 1.0
182 | Mul (e1, e2) -> simplify_add (Ln e1) (Ln e2)
183 | Div (e1, e2) -> simplify_sub (Ln e1) (Ln e2)
184 | Pow (e, Const c) -> simplify_mul (Const c) (Ln e)
185 | e -> Ln e
186
187and simplify_log base_e arg =
188 match (base_e, arg) with
189 | Const b, Const a when b > 0.0 && b <> 1.0 && a > 0.0 ->
190 Const (log a /. log b)
191 | b, a when Canonical.equal b a -> Const 1.0
192 | Const b, Pow (e, c) when Canonical.equal (Const b) e ->
193 c
194 | _ -> Log (base_e, arg)
195
196and simplify_sqrt = function
197 | Const 0.0 -> Const 0.0
198 | Const 1.0 -> Const 1.0
199 | Const c when c >= 0.0 -> Const (sqrt c)
200 | Pow (e, Const 2.0) -> simplify_abs e
201 | Mul (e1, e2) -> simplify_mul (Sqrt e1) (Sqrt e2)
202 | e -> Sqrt e
203
204and simplify_abs = function
205 | Const c -> Const (abs_float c)
206 | Abs e -> Abs e
207 | Neg e -> Abs e
208 | Mul (Const c, e) -> simplify_mul (Const (abs_float c)) (Abs e)
209 | e -> Abs e
210
211let simplify expr =
212 let rec fixed_point e count =
213 if count >= max_iterations then e
214 else
215 let canonical = Canonical.canonicalize e in
216 let simplified = simplify_once canonical in
217 if Canonical.equal simplified canonical then simplified
218 else fixed_point simplified (count + 1)
219 in
220 fixed_point expr 0