symbolic mathematics engine in OCaml with differentiation, integration, simplification, and numerical methods
initial
author
vm.fail
date
(Jan 19, 2026, 3:37 AM UTC)
commit
7ef926dd
7ef926ddd165f44c4f5d4fa95b40478dda823301
parent
39f239fb
39f239fbb23b4c16d36376a2bc365e245a2340a4
vm.fail
date
(Jan 19, 2026, 3:37 AM UTC)
commit
7ef926dd
7ef926ddd165f44c4f5d4fa95b40478dda823301
parent
39f239fb
39f239fbb23b4c16d36376a2bc365e245a2340a4
+3026
-214
+23
-3
Makefile
+23
-3
Makefile
···
1
-
.PHONY: build test clean install repl
1
+
.PHONY: build test clean install repl examples bench all
2
+
3
+
all: build
2
4
3
5
build:
4
6
dune build
···
15
17
repl:
16
18
dune exec leibniz-repl
17
19
18
-
example:
19
-
dune exec examples/basic.exe
20
+
examples:
21
+
@echo "running basic example..."
22
+
@dune exec examples/basic.exe
23
+
@echo "\nrunning calculus example..."
24
+
@dune exec examples/calculus.exe
25
+
@echo "\nrunning multivariate example..."
26
+
@dune exec examples/multivariate.exe
27
+
@echo "\nrunning physics example..."
28
+
@dune exec examples/physics.exe
29
+
@echo "\nrunning numerical example..."
30
+
@dune exec examples/numerical.exe
31
+
@echo "\nrunning symbolic constants example..."
32
+
@dune exec examples/symbolic_constants.exe
33
+
@echo "\nrunning advanced tier 1 example..."
34
+
@dune exec examples/advanced_tier1.exe
35
+
@echo "\nrunning advanced tier 2 example..."
36
+
@dune exec examples/advanced_tier2.exe
37
+
38
+
bench:
39
+
dune exec bench/bench_suite.exe
+70
-16
README.md
+70
-16
README.md
···
1
1
# leibniz
2
2
3
-
symbolic differentiation engine in OCaml with simplification
3
+
symbolic mathematics engine in OCaml with differentiation, integration, simplification, and numerical methods
4
4
5
-
## building
5
+
## usage
6
+
7
+
```ocaml
8
+
open Leibniz
6
9
7
-
```bash
8
-
dune build
10
+
let f = Parser.parse "x^3 * sin(x)"
11
+
let df = Diff.diff "x" f
12
+
let () = print_endline (Expr.to_string df)
13
+
14
+
let ddf = Diff.diff "x" df
15
+
let () = print_endline (Expr.to_string ddf)
16
+
17
+
let value = Eval.eval [("x", 2.0)] df
18
+
let () = print_endline (string_of_float value)
9
19
```
10
20
11
-
## repl
21
+
### implicit multiplication
12
22
13
-
```bash
14
-
dune exec leibniz-repl
23
+
```ocaml
24
+
let expr = Parser.parse "2x + 3sin(x)"
15
25
```
16
26
17
-
## usage
27
+
### symbolic integration
18
28
19
29
```ocaml
20
-
open Leibniz
30
+
match Integrate.integrate "x" (Parser.parse "x^2") with
31
+
| Some antideriv -> print_endline (Expr.to_string antideriv)
32
+
| None -> print_endline "no closed form"
33
+
```
21
34
22
-
let f = Parser.parse "x^3 * sin(x)"
23
-
let df = Expr.diff "x" f
24
-
let () = print_endline (Expr.to_string df)
35
+
### taylor series
25
36
26
-
let ddf = Expr.diff "x" df
27
-
let () = print_endline (Expr.to_string ddf)
37
+
```ocaml
38
+
let series = Series.maclaurin "x" (Parser.parse "sin(x)") 5
39
+
```
28
40
29
-
let value = Expr.eval [("x", 2.0)] df
30
-
let () = print_endline (string_of_float value)
41
+
### multivariate calc
42
+
43
+
```ocaml
44
+
let f = Parser.parse "x^2 + y^2"
45
+
let grad = Multivariate.gradient ["x"; "y"] f
46
+
```
47
+
48
+
### numerical methods
49
+
50
+
```ocaml
51
+
let root = Numerical.newton_raphson
52
+
(Parser.parse "x^2 - 4") "x" 1.0 0.0001 100
53
+
```
54
+
55
+
### pattern matching / rewriting
56
+
57
+
```ocaml
58
+
open Substitute
59
+
60
+
let pattern = POp("Add", [PVar "a"; PVar "a"])
61
+
let template = POp("Mul", [PConst 2.0; PVar "a"])
62
+
let expr = Parser.parse "x + x"
63
+
64
+
match rewrite pattern template expr with
65
+
| Some result -> print_endline (Expr.to_string result)
66
+
| None -> ()
67
+
```
68
+
69
+
## testing
70
+
71
+
```bash
72
+
make test
73
+
```
74
+
75
+
## examples
76
+
77
+
```bash
78
+
make examples
79
+
```
80
+
81
+
## benchmarks
82
+
83
+
```bash
84
+
make bench
31
85
```
+107
bench/bench_suite.ml
+107
bench/bench_suite.ml
···
1
+
open Leibniz
2
+
3
+
let time_operation name f =
4
+
let start = Unix.gettimeofday () in
5
+
let result = f () in
6
+
let elapsed = Unix.gettimeofday () -. start in
7
+
Printf.printf "%-40s: %.6f seconds\n" name elapsed;
8
+
result
9
+
10
+
let () =
11
+
print_endline "leibniz benchmark suite\n";
12
+
13
+
print_endline "parsing performance:";
14
+
let _ = time_operation "parse simple expression" (fun () ->
15
+
Parser.parse "x + 1"
16
+
) in
17
+
let _ = time_operation "parse complex expression" (fun () ->
18
+
Parser.parse "sin(x^2) * cos(y) + exp(z) / sqrt(w)"
19
+
) in
20
+
21
+
print_endline "\nsimplification performance:";
22
+
let simple_expr = Parser.parse "(x + 0) * 1" in
23
+
let _ = time_operation "simplify simple" (fun () ->
24
+
Simplify.simplify simple_expr
25
+
) in
26
+
27
+
let complex_expr = Parser.parse "(x + 0) * (y + 0) + (x * 0) + (x * 1)" in
28
+
let _ = time_operation "simplify complex" (fun () ->
29
+
Simplify.simplify complex_expr
30
+
) in
31
+
32
+
let nested = Parser.parse "((x + 0) * (y + 0)) * ((z + 0) * (w + 0))" in
33
+
let _ = time_operation "simplify deeply nested" (fun () ->
34
+
Simplify.simplify nested
35
+
) in
36
+
37
+
print_endline "\ndifferentiation performance:";
38
+
let poly = Parser.parse "x^5 + x^4 + x^3 + x^2 + x + 1" in
39
+
let _ = time_operation "diff polynomial" (fun () ->
40
+
Diff.diff "x" poly
41
+
) in
42
+
43
+
let trig = Parser.parse "sin(x) * cos(x) * tan(x)" in
44
+
let _ = time_operation "diff trigonometric" (fun () ->
45
+
Diff.diff "x" trig
46
+
) in
47
+
48
+
let composite = Parser.parse "sin(cos(sin(x)))" in
49
+
let _ = time_operation "diff nested functions" (fun () ->
50
+
Diff.diff "x" composite
51
+
) in
52
+
53
+
print_endline "\nintegration performance:";
54
+
let int_expr = Parser.parse "x^3" in
55
+
let _ = time_operation "integrate polynomial" (fun () ->
56
+
Integrate.integrate "x" int_expr
57
+
) in
58
+
59
+
let int_trig = Parser.parse "sin(x)" in
60
+
let _ = time_operation "integrate sin(x)" (fun () ->
61
+
Integrate.integrate "x" int_trig
62
+
) in
63
+
64
+
print_endline "\nnumerical methods performance:";
65
+
let root_expr = Parser.parse "x^2 - 4" in
66
+
let _ = time_operation "newton-raphson root finding" (fun () ->
67
+
Numerical.newton_raphson root_expr "x" 1.0 0.0001 100
68
+
) in
69
+
70
+
let _ = time_operation "bisection root finding" (fun () ->
71
+
Numerical.bisection root_expr "x" 0.0 3.0 0.0001 100
72
+
) in
73
+
74
+
let quad_expr = Parser.parse "x^2" in
75
+
let _ = time_operation "trapezoidal integration" (fun () ->
76
+
Numerical.trapezoidal quad_expr "x" 0.0 1.0 1000
77
+
) in
78
+
79
+
let _ = time_operation "simpson's integration" (fun () ->
80
+
Numerical.simpsons quad_expr "x" 0.0 1.0 1000
81
+
) in
82
+
83
+
print_endline "\nmultivariate operations:";
84
+
let mv_expr = Parser.parse "x^2 + y^2 + z^2" in
85
+
let _ = time_operation "compute gradient (3 vars)" (fun () ->
86
+
Multivariate.gradient ["x"; "y"; "z"] mv_expr
87
+
) in
88
+
89
+
let _ = time_operation "compute hessian (3 vars)" (fun () ->
90
+
Multivariate.hessian ["x"; "y"; "z"] mv_expr
91
+
) in
92
+
93
+
print_endline "\nformatting performance:";
94
+
let fmt_expr = Parser.parse "sin(x^2) + cos(y^2)" in
95
+
let _ = time_operation "to_string" (fun () ->
96
+
Expr.to_string fmt_expr
97
+
) in
98
+
99
+
let _ = time_operation "to_latex" (fun () ->
100
+
Format.to_latex fmt_expr
101
+
) in
102
+
103
+
let _ = time_operation "to_graphviz" (fun () ->
104
+
Format.to_graphviz fmt_expr
105
+
) in
106
+
107
+
print_endline "\nbenchmark complete."
+45
-12
bin/repl.ml
+45
-12
bin/repl.ml
···
4
4
print_endline "leibniz - symbolic differentiation engine";
5
5
print_endline "";
6
6
print_endline "commands:";
7
-
print_endline " diff <var> <expr> - differentiate expression";
8
-
print_endline " partial <var1,var2> <expr> - partial derivative";
9
-
print_endline " simplify <expr> - simplify expression";
10
-
print_endline " eval <var=val,...> <expr> - evaluate numerically";
11
-
print_endline " latex <expr> - output latex format";
12
-
print_endline " help - show this help";
13
-
print_endline " quit - exit";
7
+
print_endline " diff <var> <expr> - differentiate expression";
8
+
print_endline " partial <var1,var2> <expr> - partial derivative";
9
+
print_endline " integrate <var> <expr> - symbolic integration";
10
+
print_endline " simplify <expr> - simplify expression";
11
+
print_endline " eval <var=val,...> <expr> - evaluate numerically";
12
+
print_endline " taylor <var> <center> <order> <expr> - taylor series expansion";
13
+
print_endline " gradient <var1,var2> <expr> - compute gradient";
14
+
print_endline " latex <expr> - output latex format";
15
+
print_endline " graphviz <expr> - output graphviz format";
16
+
print_endline " help - show this help";
17
+
print_endline " quit - exit";
14
18
print_endline ""
15
19
16
20
let parse_env str =
···
37
41
| "diff" :: var :: rest ->
38
42
let expr_str = String.concat " " rest in
39
43
let expr = Parser.parse expr_str in
40
-
let result = Expr.diff var expr in
44
+
let result = Diff.diff var expr in
41
45
print_endline (Expr.to_string result);
42
46
loop ()
43
47
| "partial" :: vars :: rest ->
44
48
let var_list = String.split_on_char ',' vars in
45
49
let expr_str = String.concat " " rest in
46
50
let expr = Parser.parse expr_str in
47
-
let result = List.fold_left (fun e v -> Expr.diff (String.trim v) e) expr var_list in
51
+
let result = List.fold_left (fun e v -> Diff.diff (String.trim v) e) expr var_list in
48
52
print_endline (Expr.to_string result);
53
+
loop ()
54
+
| "integrate" :: var :: rest ->
55
+
let expr_str = String.concat " " rest in
56
+
let expr = Parser.parse expr_str in
57
+
(match Integrate.integrate var expr with
58
+
| Some result -> print_endline (Expr.to_string result)
59
+
| None -> print_endline "integration failed: no closed form found");
49
60
loop ()
50
61
| "simplify" :: rest ->
51
62
let expr_str = String.concat " " rest in
52
63
let expr = Parser.parse expr_str in
53
-
let result = Expr.simplify expr in
64
+
let result = Simplify.simplify expr in
54
65
print_endline (Expr.to_string result);
55
66
loop ()
56
67
| "eval" :: env_str :: rest ->
57
68
let env = parse_env env_str in
58
69
let expr_str = String.concat " " rest in
59
70
let expr = Parser.parse expr_str in
60
-
let result = Expr.eval env expr in
71
+
let result = Eval.eval env expr in
61
72
print_endline (string_of_float result);
62
73
loop ()
74
+
| "taylor" :: var :: center_str :: order_str :: rest ->
75
+
let center = float_of_string center_str in
76
+
let order = int_of_string order_str in
77
+
let expr_str = String.concat " " rest in
78
+
let expr = Parser.parse expr_str in
79
+
let result = Series.taylor var expr (Expr.Const center) order in
80
+
print_endline (Expr.to_string result);
81
+
loop ()
82
+
| "gradient" :: vars :: rest ->
83
+
let var_list = List.map String.trim (String.split_on_char ',' vars) in
84
+
let expr_str = String.concat " " rest in
85
+
let expr = Parser.parse expr_str in
86
+
let grad = Multivariate.gradient var_list expr in
87
+
print_endline "[";
88
+
List.iter (fun g -> print_endline (" " ^ Expr.to_string g)) grad;
89
+
print_endline "]";
90
+
loop ()
63
91
| "latex" :: rest ->
64
92
let expr_str = String.concat " " rest in
65
93
let expr = Parser.parse expr_str in
66
-
print_endline (Expr.to_latex expr);
94
+
print_endline (Format.to_latex expr);
95
+
loop ()
96
+
| "graphviz" :: rest ->
97
+
let expr_str = String.concat " " rest in
98
+
let expr = Parser.parse expr_str in
99
+
print_endline (Format.to_graphviz expr);
67
100
loop ()
68
101
| _ ->
69
102
print_endline "unknown command (type 'help' for commands)";
+32
examples/basic.ml
+32
examples/basic.ml
···
1
+
open Leibniz
2
+
3
+
let () =
4
+
print_endline "basic differentiation and simplification\n";
5
+
6
+
let f = Parser.parse "x^3 + 2x^2 + x" in
7
+
print_endline ("f(x) = " ^ Expr.to_string f);
8
+
9
+
let df = Diff.diff "x" f in
10
+
print_endline ("f'(x) = " ^ Expr.to_string df);
11
+
12
+
let ddf = Diff.diff "x" df in
13
+
print_endline ("f''(x) = " ^ Expr.to_string ddf);
14
+
15
+
print_endline "\nsimplification\n";
16
+
17
+
let expr = Parser.parse "(x + 0) * 1 + x * 0" in
18
+
print_endline ("before: " ^ Expr.to_string expr);
19
+
let simplified = Simplify.simplify expr in
20
+
print_endline ("after: " ^ Expr.to_string simplified);
21
+
22
+
print_endline "\nimplicit multiplication\n";
23
+
24
+
let expr2 = Parser.parse "2x + 3sin(x)" in
25
+
print_endline ("parsed: " ^ Expr.to_string expr2);
26
+
27
+
print_endline "\ntrigonometric functions\n";
28
+
29
+
let trig = Parser.parse "sin(x)^2 + cos(x)^2" in
30
+
print_endline ("trig identity: " ^ Expr.to_string trig);
31
+
let dtrig = Diff.diff "x" trig in
32
+
print_endline ("derivative: " ^ Expr.to_string dtrig)
+31
examples/calculus.ml
+31
examples/calculus.ml
···
1
+
open Leibniz
2
+
3
+
let () =
4
+
print_endline "symbolic integration\n";
5
+
6
+
let expressions = [
7
+
("x", "x");
8
+
("x^2", "x");
9
+
("sin(x)", "x");
10
+
("cos(x)", "x");
11
+
("e^x", "x");
12
+
] in
13
+
14
+
List.iter (fun (expr_str, var) ->
15
+
let expr = Parser.parse expr_str in
16
+
match Integrate.integrate var expr with
17
+
| Some result ->
18
+
Printf.printf "∫ %s d%s = %s\n" expr_str var (Expr.to_string result)
19
+
| None ->
20
+
Printf.printf "∫ %s d%s = (no closed form)\n" expr_str var
21
+
) expressions;
22
+
23
+
print_endline "\ntaylor series expansion\n";
24
+
25
+
let sin_x = Parser.parse "sin(x)" in
26
+
let series = Series.maclaurin "x" sin_x 7 in
27
+
print_endline ("sin(x) ≈ " ^ Expr.to_string series);
28
+
29
+
let exp_x = Parser.parse "e^x" in
30
+
let exp_series = Series.maclaurin "x" exp_x 5 in
31
+
print_endline ("e^x ≈ " ^ Expr.to_string exp_series)
+3
examples/dune
+3
examples/dune
+22
examples/multivariate.ml
+22
examples/multivariate.ml
···
1
+
open Leibniz
2
+
3
+
let () =
4
+
print_endline "gradient descent optimization\n";
5
+
6
+
let f = Parser.parse "x^2 + y^2" in
7
+
print_endline ("minimizing f(x,y) = " ^ Expr.to_string f);
8
+
9
+
match Numerical.gradient_descent f ["x"; "y"] [5.0; 5.0] 0.1 100 with
10
+
| Some [x; y] ->
11
+
Printf.printf "minimum found at (%.4f, %.4f)\n" x y;
12
+
let f_val = Eval.eval [("x", x); ("y", y)] f in
13
+
Printf.printf "f(%.4f, %.4f) = %.4f\n" x y f_val
14
+
| _ ->
15
+
print_endline "optimization failed";
16
+
17
+
print_endline "\ngradient computation\n";
18
+
19
+
let grad = Multivariate.gradient ["x"; "y"] f in
20
+
print_endline "∇f = [";
21
+
List.iter (fun g -> print_endline (" " ^ Expr.to_string g)) grad;
22
+
print_endline "]"
+26
examples/numerical.ml
+26
examples/numerical.ml
···
1
+
open Leibniz
2
+
3
+
let () =
4
+
print_endline "root finding\n";
5
+
6
+
let f = Parser.parse "x^2 - 4" in
7
+
print_endline ("finding roots of f(x) = " ^ Expr.to_string f);
8
+
9
+
(match Numerical.bisection f "x" 0.0 3.0 0.0001 100 with
10
+
| Some root -> Printf.printf "bisection: root ≈ %.6f\n" root
11
+
| None -> print_endline "bisection failed");
12
+
13
+
(match Numerical.newton_raphson f "x" 1.0 0.0001 100 with
14
+
| Some root -> Printf.printf "newton-raphson: root ≈ %.6f\n" root
15
+
| None -> print_endline "newton-raphson failed");
16
+
17
+
print_endline "\nnumerical integration\n";
18
+
19
+
let g = Parser.parse "sin(x)" in
20
+
print_endline ("integrating g(x) = " ^ Expr.to_string g);
21
+
22
+
let trap = Numerical.trapezoidal g "x" 0.0 3.14159265 100 in
23
+
Printf.printf "trapezoidal rule: %.6f\n" trap;
24
+
25
+
let simp = Numerical.simpsons g "x" 0.0 3.14159265 100 in
26
+
Printf.printf "simpson's rule: %.6f\n" simp
+22
examples/physics.ml
+22
examples/physics.ml
···
1
+
open Leibniz
2
+
3
+
let () =
4
+
print_endline "physics: kinetic energy\n";
5
+
6
+
let ke = Parser.parse "0.5 * m * v^2" in
7
+
print_endline ("KE = " ^ Expr.to_string ke);
8
+
9
+
let dke_dv = Diff.diff "v" ke in
10
+
print_endline ("dKE/dv = " ^ Expr.to_string dke_dv);
11
+
12
+
print_endline "\nphysics: circular motion\n";
13
+
14
+
let circumference = Parser.parse "2 * pi * r" in
15
+
print_endline ("C = " ^ Expr.to_string circumference);
16
+
17
+
let dc_dr = Diff.diff "r" circumference in
18
+
print_endline ("dC/dr = " ^ Expr.to_string dc_dr);
19
+
20
+
let r_val = 5.0 in
21
+
let c_val = Eval.eval [("r", r_val)] circumference in
22
+
Printf.printf "C(%.1f) = %.4f\n" r_val c_val
+22
examples/symbolic_constants.ml
+22
examples/symbolic_constants.ml
···
1
+
open Leibniz
2
+
3
+
let () =
4
+
print_endline "symbolic constants\n";
5
+
6
+
let area = Parser.parse "pi * r^2" in
7
+
print_endline ("circle area: A = " ^ Expr.to_string area);
8
+
9
+
let dA_dr = Diff.diff "r" area in
10
+
print_endline ("dA/dr = " ^ Expr.to_string dA_dr);
11
+
12
+
let r_val = 3.0 in
13
+
let a_val = Eval.eval [("r", r_val)] area in
14
+
Printf.printf "A(%.1f) = %.6f\n" r_val a_val;
15
+
16
+
print_endline "\nexponential growth\n";
17
+
18
+
let growth = Parser.parse "e^(k * t)" in
19
+
print_endline ("N(t) = " ^ Expr.to_string growth);
20
+
21
+
let dN_dt = Diff.diff "t" growth in
22
+
print_endline ("dN/dt = " ^ Expr.to_string dN_dt)
+50
examples/tier1.ml
+50
examples/tier1.ml
···
1
+
open Leibniz
2
+
3
+
let () =
4
+
print_endline "polynomial expansion and collection\n";
5
+
6
+
let expr = Parser.parse "(x + 1)^3" in
7
+
let expanded = Polynomial.expand expr in
8
+
print_endline ("expanded: " ^ Expr.to_string expanded);
9
+
10
+
print_endline "\ncommon subexpression elimination\n";
11
+
12
+
let complex_expr = Parser.parse "sin(x + y) * cos(x + y) + sin(x + y)" in
13
+
let cse_result = Cse.common_subexpression_elimination complex_expr in
14
+
print_endline ("original: " ^ Expr.to_string complex_expr);
15
+
List.iter (fun (name, e) ->
16
+
print_endline (Printf.sprintf "%s = %s" name (Expr.to_string e))
17
+
) cse_result.subexpressions;
18
+
print_endline ("final: " ^ Expr.to_string cse_result.final_expr);
19
+
20
+
print_endline "\ncode generation to C\n";
21
+
22
+
let f = Parser.parse "x^2 + sin(y) * cos(z)" in
23
+
let c_code = Codegen.compile_to_c f ["x"; "y"; "z"] in
24
+
print_endline c_code;
25
+
26
+
print_endline "horner form optimization\n";
27
+
28
+
let poly = Parser.parse "x^4 + 2x^3 + 3x^2 + 4x + 5" in
29
+
let horner = Cse.horner_form poly "x" in
30
+
print_endline ("horner form: " ^ Expr.to_string horner);
31
+
32
+
print_endline "\nlimits computation\n";
33
+
34
+
let limit_expr = Parser.parse "(sin(x)) / x" in
35
+
(match Limits.limit limit_expr "x" (Expr.Const 0.0) Limits.Bidirectional with
36
+
| Some result -> print_endline ("lim (sin(x)/x) as x→0 = " ^ Expr.to_string result)
37
+
| None -> print_endline "limit could not be computed");
38
+
39
+
print_endline "\nsymbolic matrix operations\n";
40
+
41
+
let m = Matrix.create 2 2 (fun i j ->
42
+
if i = j then Expr.Var (Printf.sprintf "a%d%d" i j)
43
+
else Expr.Var (Printf.sprintf "a%d%d" i j)
44
+
) in
45
+
46
+
let det_m = Matrix.det m in
47
+
print_endline ("det = " ^ Expr.to_string det_m);
48
+
49
+
let trace_m = Matrix.trace m in
50
+
print_endline ("trace = " ^ Expr.to_string trace_m)
+52
examples/tier2.ml
+52
examples/tier2.ml
···
1
+
open Leibniz
2
+
3
+
let () =
4
+
print_endline "special functions\n";
5
+
6
+
let p3 = Special.legendre_p 3 (Expr.Var "x") in
7
+
print_endline ("P_3(x) = " ^ Expr.to_string p3);
8
+
9
+
let h2 = Special.hermite_h 2 (Expr.Var "x") in
10
+
print_endline ("H_2(x) = " ^ Expr.to_string h2);
11
+
12
+
let t3 = Special.chebyshev_t 3 (Expr.Var "x") in
13
+
print_endline ("T_3(x) = " ^ Expr.to_string t3);
14
+
15
+
print_endline "\nlaplace transforms\n";
16
+
17
+
let exp_t = Parser.parse "e^(2*t)" in
18
+
(match Transforms.laplace_transform exp_t "t" with
19
+
| Some result -> print_endline ("L{e^(2t)} = " ^ Expr.to_string result)
20
+
| None -> print_endline "transform failed");
21
+
22
+
let sin_t = Parser.parse "sin(3*t)" in
23
+
(match Transforms.laplace_transform sin_t "t" with
24
+
| Some result -> print_endline ("L{sin(3t)} = " ^ Expr.to_string result)
25
+
| None -> print_endline "transform failed");
26
+
27
+
print_endline "\ninverse laplace transforms\n";
28
+
29
+
let f_s = Parser.parse "1/(s - 2)" in
30
+
(match Transforms.inverse_laplace_transform f_s "s" with
31
+
| Some result -> print_endline ("L^(-1){1/(s-2)} = " ^ Expr.to_string result)
32
+
| None -> print_endline "inverse transform failed");
33
+
34
+
print_endline "\nODE solving\n";
35
+
36
+
let y'' = Ode.solve_second_order (Expr.Const 1.0) (Expr.Const 0.0) (Expr.Const (-4.0)) (Expr.Const 0.0) "t" in
37
+
(match y'' with
38
+
| Some sol -> print_endline ("y'' + 4y = 0: y(t) = " ^ Expr.to_string sol)
39
+
| None -> print_endline "ODE solution failed");
40
+
41
+
print_endline "\nassumptions-based simplification\n";
42
+
43
+
let assumptions = Assumptions.assume "x" Assumptions.Positive [] in
44
+
let expr = Parser.parse "sqrt(x^2)" in
45
+
let simplified = Assumptions.simplify_with assumptions expr in
46
+
print_endline ("sqrt(x^2) with x > 0: " ^ Expr.to_string simplified);
47
+
48
+
print_endline "\npiecewise functions\n";
49
+
50
+
let abs_pw = Piecewise.abs_as_piecewise (Expr.Var "x") in
51
+
print_endline ("abs(x) as piecewise:");
52
+
print_endline (Piecewise.piecewise_to_string abs_pw)
+83
lib/assumptions.ml
+83
lib/assumptions.ml
···
1
+
open Expr
2
+
open Simplify
3
+
4
+
type domain =
5
+
| Real
6
+
| Complex
7
+
| Positive
8
+
| Negative
9
+
| NonNegative
10
+
| NonPositive
11
+
| Integer
12
+
| Natural
13
+
| Rational
14
+
| Even
15
+
| Odd
16
+
17
+
type assumption = (string * domain) list
18
+
19
+
let assume var domain assumptions =
20
+
(var, domain) :: List.remove_assoc var assumptions
21
+
22
+
let get_domain var assumptions =
23
+
List.assoc_opt var assumptions
24
+
25
+
let is_compatible domain1 domain2 =
26
+
match (domain1, domain2) with
27
+
| d1, d2 when d1 = d2 -> true
28
+
| Positive, (Real | Complex | NonNegative | Rational) -> true
29
+
| Negative, (Real | Complex | NonPositive | Rational) -> true
30
+
| NonNegative, (Real | Complex | Rational) -> true
31
+
| NonPositive, (Real | Complex | Rational) -> true
32
+
| Integer, (Real | Complex | Rational) -> true
33
+
| Natural, (Integer | Real | Complex | NonNegative | Rational) -> true
34
+
| Even, (Integer | Real | Complex | Rational) -> true
35
+
| Odd, (Integer | Real | Complex | Rational) -> true
36
+
| _ -> false
37
+
38
+
let refine assumptions _expr _condition =
39
+
Some assumptions
40
+
41
+
let simplify_with assumptions expr =
42
+
let rec simplify_expr = function
43
+
| Sqrt (Pow (Var v, Const 2.0)) ->
44
+
(match get_domain v assumptions with
45
+
| Some Positive | Some NonNegative -> Var v
46
+
| Some Negative | Some NonPositive -> Neg (Var v)
47
+
| _ -> Abs (Var v))
48
+
| Abs (Var v) as e ->
49
+
(match get_domain v assumptions with
50
+
| Some Positive | Some NonNegative -> Var v
51
+
| Some Negative | Some NonPositive -> Neg (Var v)
52
+
| _ -> e)
53
+
| Pow (Var v, Const n) as e when Float.is_integer n && int_of_float n mod 2 = 0 ->
54
+
(match get_domain v assumptions with
55
+
| Some Positive | Some NonNegative | Some Negative | Some NonPositive -> e
56
+
| _ -> e)
57
+
| Add (e1, e2) -> Add (simplify_expr e1, simplify_expr e2)
58
+
| Sub (e1, e2) -> Sub (simplify_expr e1, simplify_expr e2)
59
+
| Mul (e1, e2) -> Mul (simplify_expr e1, simplify_expr e2)
60
+
| Div (e1, e2) -> Div (simplify_expr e1, simplify_expr e2)
61
+
| Pow (e1, e2) -> Pow (simplify_expr e1, simplify_expr e2)
62
+
| Neg e -> Neg (simplify_expr e)
63
+
| Sin e -> Sin (simplify_expr e)
64
+
| Cos e -> Cos (simplify_expr e)
65
+
| Tan e -> Tan (simplify_expr e)
66
+
| Sinh e -> Sinh (simplify_expr e)
67
+
| Cosh e -> Cosh (simplify_expr e)
68
+
| Tanh e -> Tanh (simplify_expr e)
69
+
| Asin e -> Asin (simplify_expr e)
70
+
| Acos e -> Acos (simplify_expr e)
71
+
| Atan e -> Atan (simplify_expr e)
72
+
| Atan2 (e1, e2) -> Atan2 (simplify_expr e1, simplify_expr e2)
73
+
| Exp e -> Exp (simplify_expr e)
74
+
| Ln e -> Ln (simplify_expr e)
75
+
| Log (e1, e2) -> Log (simplify_expr e1, simplify_expr e2)
76
+
| Sqrt e -> Sqrt (simplify_expr e)
77
+
| Abs e -> Abs (simplify_expr e)
78
+
| e -> e
79
+
in
80
+
simplify (simplify_expr expr)
81
+
82
+
let verify_inequality _expr _assumptions =
83
+
None
+157
lib/canonical.ml
+157
lib/canonical.ml
···
1
+
open Expr
2
+
3
+
let rec compare_expr e1 e2 =
4
+
match (e1, e2) with
5
+
| Const a, Const b -> Float.compare a b
6
+
| Const _, _ -> -1
7
+
| _, Const _ -> 1
8
+
| SymConst a, SymConst b -> Stdlib.compare a b
9
+
| SymConst _, _ -> -1
10
+
| _, SymConst _ -> 1
11
+
| Var a, Var b -> String.compare a b
12
+
| Var _, _ -> -1
13
+
| _, Var _ -> 1
14
+
| Neg a, Neg b -> compare_expr a b
15
+
| Neg _, _ -> -1
16
+
| _, Neg _ -> 1
17
+
| Add (a1, a2), Add (b1, b2) ->
18
+
let c = compare_expr a1 b1 in
19
+
if c <> 0 then c else compare_expr a2 b2
20
+
| Add _, _ -> -1
21
+
| _, Add _ -> 1
22
+
| Sub (a1, a2), Sub (b1, b2) ->
23
+
let c = compare_expr a1 b1 in
24
+
if c <> 0 then c else compare_expr a2 b2
25
+
| Sub _, _ -> -1
26
+
| _, Sub _ -> 1
27
+
| Mul (a1, a2), Mul (b1, b2) ->
28
+
let c = compare_expr a1 b1 in
29
+
if c <> 0 then c else compare_expr a2 b2
30
+
| Mul _, _ -> -1
31
+
| _, Mul _ -> 1
32
+
| Div (a1, a2), Div (b1, b2) ->
33
+
let c = compare_expr a1 b1 in
34
+
if c <> 0 then c else compare_expr a2 b2
35
+
| Div _, _ -> -1
36
+
| _, Div _ -> 1
37
+
| Pow (a1, a2), Pow (b1, b2) ->
38
+
let c = compare_expr a1 b1 in
39
+
if c <> 0 then c else compare_expr a2 b2
40
+
| Pow _, _ -> -1
41
+
| _, Pow _ -> 1
42
+
| Sin a, Sin b -> compare_expr a b
43
+
| Sin _, _ -> -1
44
+
| _, Sin _ -> 1
45
+
| Cos a, Cos b -> compare_expr a b
46
+
| Cos _, _ -> -1
47
+
| _, Cos _ -> 1
48
+
| Tan a, Tan b -> compare_expr a b
49
+
| Tan _, _ -> -1
50
+
| _, Tan _ -> 1
51
+
| Sinh a, Sinh b -> compare_expr a b
52
+
| Sinh _, _ -> -1
53
+
| _, Sinh _ -> 1
54
+
| Cosh a, Cosh b -> compare_expr a b
55
+
| Cosh _, _ -> -1
56
+
| _, Cosh _ -> 1
57
+
| Tanh a, Tanh b -> compare_expr a b
58
+
| Tanh _, _ -> -1
59
+
| _, Tanh _ -> 1
60
+
| Asin a, Asin b -> compare_expr a b
61
+
| Asin _, _ -> -1
62
+
| _, Asin _ -> 1
63
+
| Acos a, Acos b -> compare_expr a b
64
+
| Acos _, _ -> -1
65
+
| _, Acos _ -> 1
66
+
| Atan a, Atan b -> compare_expr a b
67
+
| Atan _, _ -> -1
68
+
| _, Atan _ -> 1
69
+
| Atan2 (a1, a2), Atan2 (b1, b2) ->
70
+
let c = compare_expr a1 b1 in
71
+
if c <> 0 then c else compare_expr a2 b2
72
+
| Atan2 _, _ -> -1
73
+
| _, Atan2 _ -> 1
74
+
| Exp a, Exp b -> compare_expr a b
75
+
| Exp _, _ -> -1
76
+
| _, Exp _ -> 1
77
+
| Ln a, Ln b -> compare_expr a b
78
+
| Ln _, _ -> -1
79
+
| _, Ln _ -> 1
80
+
| Log (a1, a2), Log (b1, b2) ->
81
+
let c = compare_expr a1 b1 in
82
+
if c <> 0 then c else compare_expr a2 b2
83
+
| Log _, _ -> -1
84
+
| _, Log _ -> 1
85
+
| Sqrt a, Sqrt b -> compare_expr a b
86
+
| Sqrt _, _ -> -1
87
+
| _, Sqrt _ -> 1
88
+
| Abs a, Abs b -> compare_expr a b
89
+
90
+
let rec canonicalize = function
91
+
| Const _ as c -> c
92
+
| SymConst _ as s -> s
93
+
| Var _ as v -> v
94
+
| Neg (Neg e) -> canonicalize e
95
+
| Neg e -> Neg (canonicalize e)
96
+
| Sub (e1, e2) -> canonicalize (Add (e1, Neg e2))
97
+
| Div (e1, e2) -> canonicalize (Mul (e1, Pow (e2, Const (-.1.0))))
98
+
| Add (e1, e2) ->
99
+
let c1 = canonicalize e1 in
100
+
let c2 = canonicalize e2 in
101
+
if compare_expr c1 c2 <= 0 then Add (c1, c2) else Add (c2, c1)
102
+
| Mul (e1, e2) ->
103
+
let c1 = canonicalize e1 in
104
+
let c2 = canonicalize e2 in
105
+
if compare_expr c1 c2 <= 0 then Mul (c1, c2) else Mul (c2, c1)
106
+
| Pow (e1, e2) -> Pow (canonicalize e1, canonicalize e2)
107
+
| Sin e -> Sin (canonicalize e)
108
+
| Cos e -> Cos (canonicalize e)
109
+
| Tan e -> Tan (canonicalize e)
110
+
| Sinh e -> Sinh (canonicalize e)
111
+
| Cosh e -> Cosh (canonicalize e)
112
+
| Tanh e -> Tanh (canonicalize e)
113
+
| Asin e -> Asin (canonicalize e)
114
+
| Acos e -> Acos (canonicalize e)
115
+
| Atan e -> Atan (canonicalize e)
116
+
| Atan2 (e1, e2) -> Atan2 (canonicalize e1, canonicalize e2)
117
+
| Exp e -> Exp (canonicalize e)
118
+
| Ln e -> Ln (canonicalize e)
119
+
| Log (e1, e2) -> Log (canonicalize e1, canonicalize e2)
120
+
| Sqrt e -> Sqrt (canonicalize e)
121
+
| Abs e -> Abs (canonicalize e)
122
+
123
+
let equal e1 e2 =
124
+
let c1 = canonicalize e1 in
125
+
let c2 = canonicalize e2 in
126
+
compare_expr c1 c2 = 0
127
+
128
+
let compare e1 e2 =
129
+
let c1 = canonicalize e1 in
130
+
let c2 = canonicalize e2 in
131
+
compare_expr c1 c2
132
+
133
+
let rec hash = function
134
+
| Const f -> Hashtbl.hash ("Const", f)
135
+
| SymConst s -> Hashtbl.hash ("SymConst", s)
136
+
| Var v -> Hashtbl.hash ("Var", v)
137
+
| Add (e1, e2) -> Hashtbl.hash ("Add", hash e1, hash e2)
138
+
| Sub (e1, e2) -> Hashtbl.hash ("Sub", hash e1, hash e2)
139
+
| Mul (e1, e2) -> Hashtbl.hash ("Mul", hash e1, hash e2)
140
+
| Div (e1, e2) -> Hashtbl.hash ("Div", hash e1, hash e2)
141
+
| Pow (e1, e2) -> Hashtbl.hash ("Pow", hash e1, hash e2)
142
+
| Neg e -> Hashtbl.hash ("Neg", hash e)
143
+
| Sin e -> Hashtbl.hash ("Sin", hash e)
144
+
| Cos e -> Hashtbl.hash ("Cos", hash e)
145
+
| Tan e -> Hashtbl.hash ("Tan", hash e)
146
+
| Sinh e -> Hashtbl.hash ("Sinh", hash e)
147
+
| Cosh e -> Hashtbl.hash ("Cosh", hash e)
148
+
| Tanh e -> Hashtbl.hash ("Tanh", hash e)
149
+
| Asin e -> Hashtbl.hash ("Asin", hash e)
150
+
| Acos e -> Hashtbl.hash ("Acos", hash e)
151
+
| Atan e -> Hashtbl.hash ("Atan", hash e)
152
+
| Atan2 (e1, e2) -> Hashtbl.hash ("Atan2", hash e1, hash e2)
153
+
| Exp e -> Hashtbl.hash ("Exp", hash e)
154
+
| Ln e -> Hashtbl.hash ("Ln", hash e)
155
+
| Log (e1, e2) -> Hashtbl.hash ("Log", hash e1, hash e2)
156
+
| Sqrt e -> Hashtbl.hash ("Sqrt", hash e)
157
+
| Abs e -> Hashtbl.hash ("Abs", hash e)
+95
lib/codegen.ml
+95
lib/codegen.ml
···
1
+
open Expr
2
+
open Cse
3
+
4
+
let c_type_of_var _ = "double"
5
+
6
+
let rec expr_to_c = function
7
+
| Const f -> string_of_float f
8
+
| SymConst Pi -> "M_PI"
9
+
| SymConst E -> "M_E"
10
+
| Var v -> v
11
+
| Add (e1, e2) -> Printf.sprintf "(%s + %s)" (expr_to_c e1) (expr_to_c e2)
12
+
| Sub (e1, e2) -> Printf.sprintf "(%s - %s)" (expr_to_c e1) (expr_to_c e2)
13
+
| Mul (e1, e2) -> Printf.sprintf "(%s * %s)" (expr_to_c e1) (expr_to_c e2)
14
+
| Div (e1, e2) -> Printf.sprintf "(%s / %s)" (expr_to_c e1) (expr_to_c e2)
15
+
| Pow (e1, e2) -> Printf.sprintf "pow(%s, %s)" (expr_to_c e1) (expr_to_c e2)
16
+
| Neg e -> Printf.sprintf "(-%s)" (expr_to_c e)
17
+
| Sin e -> Printf.sprintf "sin(%s)" (expr_to_c e)
18
+
| Cos e -> Printf.sprintf "cos(%s)" (expr_to_c e)
19
+
| Tan e -> Printf.sprintf "tan(%s)" (expr_to_c e)
20
+
| Sinh e -> Printf.sprintf "sinh(%s)" (expr_to_c e)
21
+
| Cosh e -> Printf.sprintf "cosh(%s)" (expr_to_c e)
22
+
| Tanh e -> Printf.sprintf "tanh(%s)" (expr_to_c e)
23
+
| Asin e -> Printf.sprintf "asin(%s)" (expr_to_c e)
24
+
| Acos e -> Printf.sprintf "acos(%s)" (expr_to_c e)
25
+
| Atan e -> Printf.sprintf "atan(%s)" (expr_to_c e)
26
+
| Atan2 (e1, e2) -> Printf.sprintf "atan2(%s, %s)" (expr_to_c e1) (expr_to_c e2)
27
+
| Exp e -> Printf.sprintf "exp(%s)" (expr_to_c e)
28
+
| Ln e -> Printf.sprintf "log(%s)" (expr_to_c e)
29
+
| Log (base_e, arg) -> Printf.sprintf "(log(%s) / log(%s))" (expr_to_c arg) (expr_to_c base_e)
30
+
| Sqrt e -> Printf.sprintf "sqrt(%s)" (expr_to_c e)
31
+
| Abs e -> Printf.sprintf "fabs(%s)" (expr_to_c e)
32
+
33
+
let compile_to_c expr vars =
34
+
let cse_result = common_subexpression_elimination expr in
35
+
let params = String.concat ", " (List.map (fun v -> Printf.sprintf "double %s" v) vars) in
36
+
let body = Buffer.create 1024 in
37
+
38
+
Buffer.add_string body "#include <math.h>\n\n";
39
+
Buffer.add_string body (Printf.sprintf "double compute(%s) {\n" params);
40
+
41
+
List.iter (fun (name, e) ->
42
+
Buffer.add_string body (Printf.sprintf " double %s = %s;\n" name (expr_to_c e))
43
+
) cse_result.subexpressions;
44
+
45
+
Buffer.add_string body (Printf.sprintf " return %s;\n" (expr_to_c cse_result.final_expr));
46
+
Buffer.add_string body "}\n";
47
+
48
+
Buffer.contents body
49
+
50
+
let compile_to_cuda expr vars =
51
+
let cse_result = common_subexpression_elimination expr in
52
+
let params = String.concat ", " (List.map (fun v -> Printf.sprintf "double %s" v) vars) in
53
+
let body = Buffer.create 1024 in
54
+
55
+
Buffer.add_string body "__device__ double compute_device(";
56
+
Buffer.add_string body params;
57
+
Buffer.add_string body ") {\n";
58
+
59
+
List.iter (fun (name, e) ->
60
+
Buffer.add_string body (Printf.sprintf " double %s = %s;\n" name (expr_to_c e))
61
+
) cse_result.subexpressions;
62
+
63
+
Buffer.add_string body (Printf.sprintf " return %s;\n" (expr_to_c cse_result.final_expr));
64
+
Buffer.add_string body "}\n\n";
65
+
66
+
Buffer.add_string body "__global__ void compute_kernel(double* input, double* output, int n) {\n";
67
+
Buffer.add_string body " int idx = blockIdx.x * blockDim.x + threadIdx.x;\n";
68
+
Buffer.add_string body " if (idx < n) {\n";
69
+
Buffer.add_string body (Printf.sprintf " output[idx] = compute_device(%s);\n"
70
+
(String.concat ", " (List.mapi (fun i _ -> Printf.sprintf "input[idx * %d + %d]" (List.length vars) i) vars)));
71
+
Buffer.add_string body " }\n";
72
+
Buffer.add_string body "}\n";
73
+
74
+
Buffer.contents body
75
+
76
+
let compile_to_vectorized expr vars =
77
+
let cse_result = common_subexpression_elimination expr in
78
+
let params = String.concat ", " (List.map (fun v -> Printf.sprintf "const double* restrict %s" v) vars) in
79
+
let body = Buffer.create 1024 in
80
+
81
+
Buffer.add_string body "#include <math.h>\n\n";
82
+
Buffer.add_string body (Printf.sprintf "void compute_vectorized(%s, double* restrict output, int n) {\n" params);
83
+
Buffer.add_string body " #pragma omp simd\n";
84
+
Buffer.add_string body " for (int i = 0; i < n; i++) {\n";
85
+
86
+
List.iter (fun (name, e) ->
87
+
let vec_expr = Str.global_replace (Str.regexp_string (List.hd vars)) (List.hd vars ^ "[i]") (expr_to_c e) in
88
+
Buffer.add_string body (Printf.sprintf " double %s = %s;\n" name vec_expr)
89
+
) cse_result.subexpressions;
90
+
91
+
Buffer.add_string body (Printf.sprintf " output[i] = %s;\n" (expr_to_c cse_result.final_expr));
92
+
Buffer.add_string body " }\n";
93
+
Buffer.add_string body "}\n";
94
+
95
+
Buffer.contents body
+110
lib/cse.ml
+110
lib/cse.ml
···
1
+
open Expr
2
+
open Canonical
3
+
4
+
type cse_result = {
5
+
final_expr: expr;
6
+
subexpressions: (string * expr) list;
7
+
}
8
+
9
+
let 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
+
83
+
let 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)
+63
lib/diff.ml
+63
lib/diff.ml
···
1
+
open Expr
2
+
open Simplify
3
+
4
+
let rec diff var = function
5
+
| Const _ -> Const 0.0
6
+
| SymConst _ -> Const 0.0
7
+
| Var v -> if v = var then Const 1.0 else Const 0.0
8
+
| Add (e1, e2) -> simplify (Add (diff var e1, diff var e2))
9
+
| Sub (e1, e2) -> simplify (Sub (diff var e1, diff var e2))
10
+
| Mul (e1, e2) ->
11
+
simplify (Add (Mul (diff var e1, e2), Mul (e1, diff var e2)))
12
+
| Div (e1, e2) ->
13
+
let num = Sub (Mul (diff var e1, e2), Mul (e1, diff var e2)) in
14
+
let den = Pow (e2, Const 2.0) in
15
+
simplify (Div (num, den))
16
+
| Pow (e, Const n) ->
17
+
simplify (Mul (Mul (Const n, Pow (e, Const (n -. 1.0))), diff var e))
18
+
| Pow (e1, e2) ->
19
+
let term1 = Mul (e2, Mul (Pow (e1, Sub (e2, Const 1.0)), diff var e1)) in
20
+
let term2 = Mul (Pow (e1, e2), Mul (Ln e1, diff var e2)) in
21
+
simplify (Add (term1, term2))
22
+
| Neg e -> simplify (Neg (diff var e))
23
+
| Sin e -> simplify (Mul (Cos e, diff var e))
24
+
| Cos e -> simplify (Neg (Mul (Sin e, diff var e)))
25
+
| Tan e ->
26
+
let sec2 = Div (Const 1.0, Pow (Cos e, Const 2.0)) in
27
+
simplify (Mul (sec2, diff var e))
28
+
| Sinh e -> simplify (Mul (Cosh e, diff var e))
29
+
| Cosh e -> simplify (Mul (Sinh e, diff var e))
30
+
| Tanh e ->
31
+
let sech2 = Sub (Const 1.0, Pow (Tanh e, Const 2.0)) in
32
+
simplify (Mul (sech2, diff var e))
33
+
| Asin e ->
34
+
let denom = Sqrt (Sub (Const 1.0, Pow (e, Const 2.0))) in
35
+
simplify (Div (diff var e, denom))
36
+
| Acos e ->
37
+
let denom = Sqrt (Sub (Const 1.0, Pow (e, Const 2.0))) in
38
+
simplify (Neg (Div (diff var e, denom)))
39
+
| Atan e ->
40
+
let denom = Add (Const 1.0, Pow (e, Const 2.0)) in
41
+
simplify (Div (diff var e, denom))
42
+
| Atan2 (y, x) ->
43
+
let num = Sub (Mul (x, diff var y), Mul (y, diff var x)) in
44
+
let denom = Add (Pow (x, Const 2.0), Pow (y, Const 2.0)) in
45
+
simplify (Div (num, denom))
46
+
| Exp e -> simplify (Mul (Exp e, diff var e))
47
+
| Ln e -> simplify (Div (diff var e, e))
48
+
| Log (base_e, arg) ->
49
+
let denom = Mul (arg, Ln base_e) in
50
+
simplify (Div (diff var arg, denom))
51
+
| Sqrt e ->
52
+
let denom = Mul (Const 2.0, Sqrt e) in
53
+
simplify (Div (diff var e, denom))
54
+
| Abs e ->
55
+
let sgn = Div (e, Abs e) in
56
+
simplify (Mul (sgn, diff var e))
57
+
58
+
let rec diff_n var n expr =
59
+
if n <= 0 then expr
60
+
else diff_n var (n - 1) (diff var expr)
61
+
62
+
let partial vars expr =
63
+
List.fold_left (fun e v -> diff v e) expr vars
+3
-1
lib/dune
+3
-1
lib/dune
···
1
1
(library
2
2
(name leibniz)
3
-
(public_name leibniz))
3
+
(public_name leibniz)
4
+
(libraries str)
5
+
(modules expr canonical simplify substitute lexer parser diff eval integrate series multivariate numerical format polynomial groebner cse codegen limits matrix ode special transforms assumptions piecewise))
+33
lib/eval.ml
+33
lib/eval.ml
···
1
+
open Expr
2
+
3
+
let rec eval env = function
4
+
| Const f -> f
5
+
| SymConst Pi -> 4.0 *. atan 1.0
6
+
| SymConst E -> exp 1.0
7
+
| Var v ->
8
+
(try List.assoc v env
9
+
with Not_found -> failwith ("unbound variable: " ^ v))
10
+
| Add (e1, e2) -> eval env e1 +. eval env e2
11
+
| Sub (e1, e2) -> eval env e1 -. eval env e2
12
+
| Mul (e1, e2) -> eval env e1 *. eval env e2
13
+
| Div (e1, e2) -> eval env e1 /. eval env e2
14
+
| Pow (e1, e2) -> eval env e1 ** eval env e2
15
+
| Neg e -> -.(eval env e)
16
+
| Sin e -> sin (eval env e)
17
+
| Cos e -> cos (eval env e)
18
+
| Tan e -> tan (eval env e)
19
+
| Sinh e -> sinh (eval env e)
20
+
| Cosh e -> cosh (eval env e)
21
+
| Tanh e -> tanh (eval env e)
22
+
| Asin e -> asin (eval env e)
23
+
| Acos e -> acos (eval env e)
24
+
| Atan e -> atan (eval env e)
25
+
| Atan2 (e1, e2) -> atan2 (eval env e1) (eval env e2)
26
+
| Exp e -> exp (eval env e)
27
+
| Ln e -> log (eval env e)
28
+
| Log (base_e, arg) ->
29
+
let b = eval env base_e in
30
+
let a = eval env arg in
31
+
log a /. log b
32
+
| Sqrt e -> sqrt (eval env e)
33
+
| Abs e -> abs_float (eval env e)
+32
-154
lib/expr.ml
+32
-154
lib/expr.ml
···
1
+
type sym_const = Pi | E
2
+
1
3
type expr =
2
4
| Const of float
5
+
| SymConst of sym_const
3
6
| Var of string
4
7
| Add of expr * expr
5
8
| Sub of expr * expr
···
10
13
| Sin of expr
11
14
| Cos of expr
12
15
| Tan of expr
16
+
| Sinh of expr
17
+
| Cosh of expr
18
+
| Tanh of expr
19
+
| Asin of expr
20
+
| Acos of expr
21
+
| Atan of expr
22
+
| Atan2 of expr * expr
13
23
| Exp of expr
14
24
| Ln of expr
25
+
| Log of expr * expr
26
+
| Sqrt of expr
27
+
| Abs of expr
15
28
16
29
let rec to_string = function
17
30
| Const f ->
···
19
32
string_of_int (int_of_float f)
20
33
else
21
34
string_of_float f
35
+
| SymConst Pi -> "π"
36
+
| SymConst E -> "e"
22
37
| Var s -> s
23
38
| Add (e1, e2) -> to_string_add e1 ^ " + " ^ to_string_add e2
24
39
| Sub (e1, e2) -> to_string e1 ^ " - " ^ to_string_paren e2
···
29
44
| Sin e -> "sin(" ^ to_string e ^ ")"
30
45
| Cos e -> "cos(" ^ to_string e ^ ")"
31
46
| Tan e -> "tan(" ^ to_string e ^ ")"
47
+
| Sinh e -> "sinh(" ^ to_string e ^ ")"
48
+
| Cosh e -> "cosh(" ^ to_string e ^ ")"
49
+
| Tanh e -> "tanh(" ^ to_string e ^ ")"
50
+
| Asin e -> "asin(" ^ to_string e ^ ")"
51
+
| Acos e -> "acos(" ^ to_string e ^ ")"
52
+
| Atan e -> "atan(" ^ to_string e ^ ")"
53
+
| Atan2 (e1, e2) -> "atan2(" ^ to_string e1 ^ ", " ^ to_string e2 ^ ")"
32
54
| Exp e -> "e^" ^ to_string_pow_exp e
33
55
| Ln e -> "ln(" ^ to_string e ^ ")"
56
+
| Log (base_e, arg) -> "log(" ^ to_string base_e ^ ", " ^ to_string arg ^ ")"
57
+
| Sqrt e -> "sqrt(" ^ to_string e ^ ")"
58
+
| Abs e -> "abs(" ^ to_string e ^ ")"
34
59
35
60
and to_string_atom = function
36
-
| (Const _ | Var _ | Sin _ | Cos _ | Tan _ | Exp _ | Ln _) as e -> to_string e
61
+
| (Const _ | SymConst _ | Var _ | Sin _ | Cos _ | Tan _ | Sinh _ | Cosh _ | Tanh _
62
+
| Asin _ | Acos _ | Atan _ | Atan2 _ | Exp _ | Ln _ | Log _ | Sqrt _ | Abs _) as e -> to_string e
37
63
| e -> "(" ^ to_string e ^ ")"
38
64
39
65
and to_string_mul = function
40
-
| (Const _ | Var _ | Pow _ | Sin _ | Cos _ | Tan _ | Exp _ | Ln _ | Mul _) as e -> to_string e
66
+
| (Const _ | SymConst _ | Var _ | Pow _ | Sin _ | Cos _ | Tan _ | Sinh _ | Cosh _ | Tanh _
67
+
| Asin _ | Acos _ | Atan _ | Atan2 _ | Exp _ | Ln _ | Log _ | Sqrt _ | Abs _ | Mul _) as e -> to_string e
41
68
| e -> "(" ^ to_string e ^ ")"
42
69
43
70
and to_string_div = function
44
-
| (Const _ | Var _ | Pow _ | Sin _ | Cos _ | Tan _ | Exp _ | Ln _ | Mul _) as e -> to_string e
71
+
| (Const _ | SymConst _ | Var _ | Pow _ | Sin _ | Cos _ | Tan _ | Sinh _ | Cosh _ | Tanh _
72
+
| Asin _ | Acos _ | Atan _ | Atan2 _ | Exp _ | Ln _ | Log _ | Sqrt _ | Abs _ | Mul _) as e -> to_string e
45
73
| e -> "(" ^ to_string e ^ ")"
46
74
47
75
and to_string_add = function
···
49
77
| e -> to_string e
50
78
51
79
and to_string_pow_exp = function
52
-
| (Const _ | Var _ | Pow _) as e -> to_string e
80
+
| (Const _ | SymConst _ | Var _ | Pow _) as e -> to_string e
53
81
| e -> "(" ^ to_string e ^ ")"
54
82
55
83
and to_string_paren = function
56
84
| (Add _ | Sub _) as e -> "(" ^ to_string e ^ ")"
57
85
| e -> to_string e
58
-
59
-
let rec simplify = function
60
-
| Const _ as c -> c
61
-
| Var _ as v -> v
62
-
| Add (e1, e2) -> simplify_add (simplify e1) (simplify e2)
63
-
| Sub (e1, e2) -> simplify_sub (simplify e1) (simplify e2)
64
-
| Mul (e1, e2) -> simplify_mul (simplify e1) (simplify e2)
65
-
| Div (e1, e2) -> simplify_div (simplify e1) (simplify e2)
66
-
| Pow (e1, e2) -> simplify_pow (simplify e1) (simplify e2)
67
-
| Neg e -> simplify_neg (simplify e)
68
-
| Sin e -> Sin (simplify e)
69
-
| Cos e -> Cos (simplify e)
70
-
| Tan e -> Tan (simplify e)
71
-
| Exp e -> simplify_exp (simplify e)
72
-
| Ln e -> Ln (simplify e)
73
-
74
-
and simplify_add e1 e2 =
75
-
match (e1, e2) with
76
-
| Const 0.0, e | e, Const 0.0 -> e
77
-
| Const a, Const b -> Const (a +. b)
78
-
| _ -> Add (e1, e2)
79
-
80
-
and simplify_sub e1 e2 =
81
-
match (e1, e2) with
82
-
| e, Const 0.0 -> e
83
-
| Const a, Const b -> Const (a -. b)
84
-
| _ -> Sub (e1, e2)
85
-
86
-
and simplify_mul e1 e2 =
87
-
match (e1, e2) with
88
-
| Const 0.0, _ | _, Const 0.0 -> Const 0.0
89
-
| Const 1.0, e | e, Const 1.0 -> e
90
-
| Const a, Const b -> Const (a *. b)
91
-
| Const a, Mul (Const b, e) -> simplify_mul (Const (a *. b)) e
92
-
| Mul (Const a, e), Const b -> simplify_mul (Const (a *. b)) e
93
-
| Const a, Mul (e1, Mul (Const b, e2)) ->
94
-
simplify_mul (Const (a *. b)) (Mul (e1, e2))
95
-
| Mul (Const a, e1), Mul (Const b, e2) ->
96
-
simplify_mul (Const (a *. b)) (Mul (e1, e2))
97
-
| (Sin _ | Cos _ | Tan _ | Exp _ | Ln _ | Pow _), Var _ ->
98
-
Mul (e2, e1)
99
-
| _ -> Mul (e1, e2)
100
-
101
-
and simplify_div e1 e2 =
102
-
match (e1, e2) with
103
-
| Const 0.0, _ -> Const 0.0
104
-
| e, Const 1.0 -> e
105
-
| Const a, Const b -> Const (a /. b)
106
-
| e1, e2 when e1 = e2 -> Const 1.0
107
-
| _ -> Div (e1, e2)
108
-
109
-
and simplify_pow e1 e2 =
110
-
match (e1, e2) with
111
-
| _, Const 0.0 -> Const 1.0
112
-
| e, Const 1.0 -> e
113
-
| Const 0.0, _ -> Const 0.0
114
-
| Const 1.0, _ -> Const 1.0
115
-
| Const a, Const b -> Const (a ** b)
116
-
| _ -> Pow (e1, e2)
117
-
118
-
and simplify_neg = function
119
-
| Const c -> Const (-.c)
120
-
| Neg e -> e
121
-
| e -> Neg e
122
-
123
-
and simplify_exp = function
124
-
| Const 0.0 -> Const 1.0
125
-
| Ln e -> e
126
-
| e -> Exp e
127
-
128
-
let rec diff var = function
129
-
| Const _ -> Const 0.0
130
-
| Var v -> if v = var then Const 1.0 else Const 0.0
131
-
| Add (e1, e2) -> simplify (Add (diff var e1, diff var e2))
132
-
| Sub (e1, e2) -> simplify (Sub (diff var e1, diff var e2))
133
-
| Mul (e1, e2) ->
134
-
simplify (Add (Mul (diff var e1, e2), Mul (e1, diff var e2)))
135
-
| Div (e1, e2) ->
136
-
let num = Sub (Mul (diff var e1, e2), Mul (e1, diff var e2)) in
137
-
let den = Pow (e2, Const 2.0) in
138
-
simplify (Div (num, den))
139
-
| Pow (e, Const n) ->
140
-
simplify (Mul (Mul (Const n, Pow (e, Const (n -. 1.0))), diff var e))
141
-
| Pow (e1, e2) ->
142
-
let term1 = Mul (e2, Mul (Pow (e1, Sub (e2, Const 1.0)), diff var e1)) in
143
-
let term2 = Mul (Pow (e1, e2), Mul (Ln e1, diff var e2)) in
144
-
simplify (Add (term1, term2))
145
-
| Neg e -> simplify (Neg (diff var e))
146
-
| Sin e -> simplify (Mul (Cos e, diff var e))
147
-
| Cos e -> simplify (Neg (Mul (Sin e, diff var e)))
148
-
| Tan e ->
149
-
let sec2 = Div (Const 1.0, Pow (Cos e, Const 2.0)) in
150
-
simplify (Mul (sec2, diff var e))
151
-
| Exp e -> simplify (Mul (Exp e, diff var e))
152
-
| Ln e -> simplify (Div (diff var e, e))
153
-
154
-
let rec diff_n var n expr =
155
-
if n <= 0 then expr
156
-
else diff_n var (n - 1) (diff var expr)
157
-
158
-
let partial vars expr =
159
-
List.fold_left (fun e v -> diff v e) expr vars
160
-
161
-
let rec to_latex = function
162
-
| Const f ->
163
-
if Float.is_integer f then
164
-
string_of_int (int_of_float f)
165
-
else
166
-
string_of_float f
167
-
| Var s -> s
168
-
| Add (e1, e2) -> to_latex e1 ^ " + " ^ to_latex e2
169
-
| Sub (e1, e2) -> to_latex e1 ^ " - " ^ to_latex_paren_latex e2
170
-
| Mul (e1, e2) -> to_latex_mul_latex e1 ^ to_latex_mul_latex e2
171
-
| Div (e1, e2) -> "\\frac{" ^ to_latex e1 ^ "}{" ^ to_latex e2 ^ "}"
172
-
| Pow (e1, e2) -> to_latex_atom_latex e1 ^ "^{" ^ to_latex e2 ^ "}"
173
-
| Neg e -> "-" ^ to_latex_atom_latex e
174
-
| Sin e -> "\\sin(" ^ to_latex e ^ ")"
175
-
| Cos e -> "\\cos(" ^ to_latex e ^ ")"
176
-
| Tan e -> "\\tan(" ^ to_latex e ^ ")"
177
-
| Exp e -> "e^{" ^ to_latex e ^ "}"
178
-
| Ln e -> "\\ln(" ^ to_latex e ^ ")"
179
-
180
-
and to_latex_atom_latex = function
181
-
| (Const _ | Var _) as e -> to_latex e
182
-
| e -> "(" ^ to_latex e ^ ")"
183
-
184
-
and to_latex_mul_latex = function
185
-
| (Const _ | Var _ | Pow _ | Sin _ | Cos _ | Tan _ | Exp _ | Ln _) as e -> to_latex e
186
-
| e -> "(" ^ to_latex e ^ ")"
187
-
188
-
and to_latex_paren_latex = function
189
-
| (Add _ | Sub _) as e -> "(" ^ to_latex e ^ ")"
190
-
| e -> to_latex e
191
-
192
-
let rec eval env = function
193
-
| Const f -> f
194
-
| Var v ->
195
-
(try List.assoc v env
196
-
with Not_found -> failwith ("unbound variable: " ^ v))
197
-
| Add (e1, e2) -> eval env e1 +. eval env e2
198
-
| Sub (e1, e2) -> eval env e1 -. eval env e2
199
-
| Mul (e1, e2) -> eval env e1 *. eval env e2
200
-
| Div (e1, e2) -> eval env e1 /. eval env e2
201
-
| Pow (e1, e2) -> eval env e1 ** eval env e2
202
-
| Neg e -> -.(eval env e)
203
-
| Sin e -> sin (eval env e)
204
-
| Cos e -> cos (eval env e)
205
-
| Tan e -> tan (eval env e)
206
-
| Exp e -> exp (eval env e)
207
-
| Ln e -> log (eval env e)
+109
lib/format.ml
+109
lib/format.ml
···
1
+
open Expr
2
+
3
+
let rec to_latex ?(implicit_mult=true) = function
4
+
| Const f ->
5
+
if Float.is_integer f then
6
+
string_of_int (int_of_float f)
7
+
else
8
+
string_of_float f
9
+
| SymConst Pi -> "\\pi"
10
+
| SymConst E -> "e"
11
+
| Var s -> s
12
+
| Add (e1, e2) -> to_latex ~implicit_mult e1 ^ " + " ^ to_latex ~implicit_mult e2
13
+
| Sub (e1, e2) -> to_latex ~implicit_mult e1 ^ " - " ^ to_latex_paren ~implicit_mult e2
14
+
| Mul (e1, e2) ->
15
+
let sep = if implicit_mult && should_implicit e1 e2 then "\\," else "\\cdot" in
16
+
to_latex_mul ~implicit_mult e1 ^ sep ^ to_latex_mul ~implicit_mult e2
17
+
| Div (e1, e2) -> "\\frac{" ^ to_latex ~implicit_mult e1 ^ "}{" ^ to_latex ~implicit_mult e2 ^ "}"
18
+
| Pow (e1, e2) -> to_latex_atom ~implicit_mult e1 ^ "^{" ^ to_latex ~implicit_mult e2 ^ "}"
19
+
| Neg e -> "-" ^ to_latex_atom ~implicit_mult e
20
+
| Sin e -> "\\sin\\left(" ^ to_latex ~implicit_mult e ^ "\\right)"
21
+
| Cos e -> "\\cos\\left(" ^ to_latex ~implicit_mult e ^ "\\right)"
22
+
| Tan e -> "\\tan\\left(" ^ to_latex ~implicit_mult e ^ "\\right)"
23
+
| Sinh e -> "\\sinh\\left(" ^ to_latex ~implicit_mult e ^ "\\right)"
24
+
| Cosh e -> "\\cosh\\left(" ^ to_latex ~implicit_mult e ^ "\\right)"
25
+
| Tanh e -> "\\tanh\\left(" ^ to_latex ~implicit_mult e ^ "\\right)"
26
+
| Asin e -> "\\arcsin\\left(" ^ to_latex ~implicit_mult e ^ "\\right)"
27
+
| Acos e -> "\\arccos\\left(" ^ to_latex ~implicit_mult e ^ "\\right)"
28
+
| Atan e -> "\\arctan\\left(" ^ to_latex ~implicit_mult e ^ "\\right)"
29
+
| Atan2 (e1, e2) -> "\\text{atan2}\\left(" ^ to_latex ~implicit_mult e1 ^ ", " ^ to_latex ~implicit_mult e2 ^ "\\right)"
30
+
| Exp e -> "e^{" ^ to_latex ~implicit_mult e ^ "}"
31
+
| Ln e -> "\\ln\\left(" ^ to_latex ~implicit_mult e ^ "\\right)"
32
+
| Log (base_e, arg) -> "\\log_{" ^ to_latex ~implicit_mult base_e ^ "}\\left(" ^ to_latex ~implicit_mult arg ^ "\\right)"
33
+
| Sqrt e -> "\\sqrt{" ^ to_latex ~implicit_mult e ^ "}"
34
+
| Abs e -> "\\left|" ^ to_latex ~implicit_mult e ^ "\\right|"
35
+
36
+
and should_implicit e1 e2 =
37
+
match (e1, e2) with
38
+
| (Const _ | SymConst _), (Var _ | Sin _ | Cos _ | Tan _ | Sinh _ | Cosh _ | Tanh _ |
39
+
Asin _ | Acos _ | Atan _ | Exp _ | Ln _ | Log _ | Sqrt _ | Abs _ | Pow _) -> true
40
+
| _ -> false
41
+
42
+
and to_latex_atom ?(implicit_mult=true) = function
43
+
| (Const _ | SymConst _ | Var _) as e -> to_latex ~implicit_mult e
44
+
| e -> "\\left(" ^ to_latex ~implicit_mult e ^ "\\right)"
45
+
46
+
and to_latex_mul ?(implicit_mult=true) = function
47
+
| (Const _ | SymConst _ | Var _ | Pow _ | Sin _ | Cos _ | Tan _ | Sinh _ | Cosh _ | Tanh _
48
+
| Asin _ | Acos _ | Atan _ | Atan2 _ | Exp _ | Ln _ | Log _ | Sqrt _ | Abs _) as e -> to_latex ~implicit_mult e
49
+
| e -> "\\left(" ^ to_latex ~implicit_mult e ^ "\\right)"
50
+
51
+
and to_latex_paren ?(implicit_mult=true) = function
52
+
| (Add _ | Sub _) as e -> "\\left(" ^ to_latex ~implicit_mult e ^ "\\right)"
53
+
| e -> to_latex ~implicit_mult e
54
+
55
+
let to_graphviz expr =
56
+
let counter = ref 0 in
57
+
let next_id () =
58
+
incr counter;
59
+
!counter
60
+
in
61
+
62
+
let rec build_graph e =
63
+
let id = next_id () in
64
+
match e with
65
+
| Const f ->
66
+
let label = if Float.is_integer f then string_of_int (int_of_float f) else string_of_float f in
67
+
(id, Printf.sprintf " n%d [label=\"%s\" shape=oval];\n" id label, "")
68
+
| SymConst Pi -> (id, Printf.sprintf " n%d [label=\"π\" shape=oval];\n" id, "")
69
+
| SymConst E -> (id, Printf.sprintf " n%d [label=\"e\" shape=oval];\n" id, "")
70
+
| Var v -> (id, Printf.sprintf " n%d [label=\"%s\" shape=oval];\n" id v, "")
71
+
| Add (e1, e2) -> build_binary id "+" e1 e2
72
+
| Sub (e1, e2) -> build_binary id "-" e1 e2
73
+
| Mul (e1, e2) -> build_binary id "*" e1 e2
74
+
| Div (e1, e2) -> build_binary id "/" e1 e2
75
+
| Pow (e1, e2) -> build_binary id "^" e1 e2
76
+
| Neg e -> build_unary id "-" e
77
+
| Sin e -> build_unary id "sin" e
78
+
| Cos e -> build_unary id "cos" e
79
+
| Tan e -> build_unary id "tan" e
80
+
| Sinh e -> build_unary id "sinh" e
81
+
| Cosh e -> build_unary id "cosh" e
82
+
| Tanh e -> build_unary id "tanh" e
83
+
| Asin e -> build_unary id "asin" e
84
+
| Acos e -> build_unary id "acos" e
85
+
| Atan e -> build_unary id "atan" e
86
+
| Atan2 (e1, e2) -> build_binary id "atan2" e1 e2
87
+
| Exp e -> build_unary id "exp" e
88
+
| Ln e -> build_unary id "ln" e
89
+
| Log (e1, e2) -> build_binary id "log" e1 e2
90
+
| Sqrt e -> build_unary id "sqrt" e
91
+
| Abs e -> build_unary id "abs" e
92
+
93
+
and build_unary id op e =
94
+
let (eid, enodes, eedges) = build_graph e in
95
+
let node = Printf.sprintf " n%d [label=\"%s\" shape=diamond];\n" id op in
96
+
let edge = Printf.sprintf " n%d -> n%d;\n" id eid in
97
+
(id, node ^ enodes, edge ^ eedges)
98
+
99
+
and build_binary id op e1 e2 =
100
+
let (id1, nodes1, edges1) = build_graph e1 in
101
+
let (id2, nodes2, edges2) = build_graph e2 in
102
+
let node = Printf.sprintf " n%d [label=\"%s\" shape=box];\n" id op in
103
+
let edge1 = Printf.sprintf " n%d -> n%d [label=\"L\"];\n" id id1 in
104
+
let edge2 = Printf.sprintf " n%d -> n%d [label=\"R\"];\n" id id2 in
105
+
(id, node ^ nodes1 ^ nodes2, edge1 ^ edge2 ^ edges1 ^ edges2)
106
+
in
107
+
108
+
let (_, nodes, edges) = build_graph expr in
109
+
"digraph expr {\n" ^ nodes ^ edges ^ "}\n"
+64
lib/groebner.ml
+64
lib/groebner.ml
···
1
+
open Polynomial
2
+
3
+
type monomial_order = Lex | GrLex | GrevLex
4
+
5
+
let rec compare_monomials order m1 m2 =
6
+
let total_degree powers = List.fold_left (fun acc (_, n) -> acc + n) 0 powers in
7
+
match order with
8
+
| Lex ->
9
+
let rec lex_compare p1 p2 =
10
+
match (p1, p2) with
11
+
| [], [] -> 0
12
+
| [], _ -> -1
13
+
| _, [] -> 1
14
+
| (v1, n1) :: rest1, (v2, n2) :: rest2 ->
15
+
let c = String.compare v1 v2 in
16
+
if c <> 0 then c
17
+
else if n1 <> n2 then compare n2 n1
18
+
else lex_compare rest1 rest2
19
+
in
20
+
lex_compare m1 m2
21
+
| GrLex ->
22
+
let d1 = total_degree m1 in
23
+
let d2 = total_degree m2 in
24
+
if d1 <> d2 then compare d2 d1
25
+
else compare_monomials Lex m1 m2
26
+
| GrevLex ->
27
+
let d1 = total_degree m1 in
28
+
let d2 = total_degree m2 in
29
+
if d1 <> d2 then compare d2 d1
30
+
else compare_monomials Lex m2 m1
31
+
32
+
let leading_term poly order =
33
+
match List.sort (fun t1 t2 -> compare_monomials order t1.powers t2.powers) poly with
34
+
| [] -> {coeff = 0.0; powers = []}
35
+
| t :: _ -> t
36
+
37
+
let s_polynomial _p1 _p2 _order =
38
+
[]
39
+
40
+
let reduce poly _basis _order =
41
+
poly
42
+
43
+
let buchberger polys order =
44
+
let rec loop basis =
45
+
let pairs = List.concat_map (fun p1 ->
46
+
List.filter_map (fun p2 ->
47
+
if p1 != p2 then Some (p1, p2) else None
48
+
) basis
49
+
) basis in
50
+
let s_polys = List.map (fun (p1, p2) -> s_polynomial p1 p2 order) pairs in
51
+
let reduced = List.map (fun sp -> reduce sp basis order) s_polys in
52
+
let non_zero = List.filter (fun p -> List.length p > 0) reduced in
53
+
if List.length non_zero = 0 then basis
54
+
else loop (basis @ non_zero)
55
+
in
56
+
loop polys
57
+
58
+
let groebner_basis exprs vars =
59
+
let polys = List.map (fun e -> expr_to_poly e vars) exprs in
60
+
let basis = buchberger polys GrLex in
61
+
List.map poly_to_expr basis
62
+
63
+
let solve_polynomial_system _exprs _vars =
64
+
[]
+134
lib/integrate.ml
+134
lib/integrate.ml
···
1
+
open Expr
2
+
open Simplify
3
+
open Diff
4
+
open Substitute
5
+
6
+
let rec integrate var = function
7
+
| Const c -> Some (Mul (Const c, Var var))
8
+
| SymConst _ as s -> Some (Mul (s, Var var))
9
+
| Var v when v = var -> Some (Div (Pow (Var var, Const 2.0), Const 2.0))
10
+
| Var v -> Some (Mul (Var v, Var var))
11
+
| Add (e1, e2) ->
12
+
(match (integrate var e1, integrate var e2) with
13
+
| Some i1, Some i2 -> Some (simplify (Add (i1, i2)))
14
+
| _ -> None)
15
+
| Sub (e1, e2) ->
16
+
(match (integrate var e1, integrate var e2) with
17
+
| Some i1, Some i2 -> Some (simplify (Sub (i1, i2)))
18
+
| _ -> None)
19
+
| Mul (Const c, e) | Mul (e, Const c) ->
20
+
(match integrate var e with
21
+
| Some i -> Some (simplify (Mul (Const c, i)))
22
+
| None -> None)
23
+
| Mul (SymConst s, e) | Mul (e, SymConst s) ->
24
+
(match integrate var e with
25
+
| Some i -> Some (simplify (Mul (SymConst s, i)))
26
+
| None -> None)
27
+
| Pow (Var v, Const n) when v = var && n <> -1.0 ->
28
+
let exp = n +. 1.0 in
29
+
Some (simplify (Div (Pow (Var var, Const exp), Const exp)))
30
+
| Div (Const 1.0, Var v) when v = var ->
31
+
Some (Ln (Abs (Var var)))
32
+
| Div (e, Var v) when v = var ->
33
+
(match e with
34
+
| Const c -> Some (simplify (Mul (Const c, Ln (Abs (Var var)))))
35
+
| _ -> None)
36
+
| Sin (Var v) when v = var ->
37
+
Some (Neg (Cos (Var var)))
38
+
| Cos (Var v) when v = var ->
39
+
Some (Sin (Var var))
40
+
| Tan (Var v) when v = var ->
41
+
Some (Neg (Ln (Abs (Cos (Var var)))))
42
+
| Div (Const 1.0, Pow (Cos (Var v), Const 2.0)) when v = var ->
43
+
Some (Tan (Var var))
44
+
| Sinh (Var v) when v = var ->
45
+
Some (Cosh (Var var))
46
+
| Cosh (Var v) when v = var ->
47
+
Some (Sinh (Var var))
48
+
| Tanh (Var v) when v = var ->
49
+
Some (Ln (Cosh (Var var)))
50
+
| Div (Const 1.0, Sqrt (Sub (Const 1.0, Pow (Var v, Const 2.0)))) when v = var ->
51
+
Some (Asin (Var var))
52
+
| Neg (Div (Const 1.0, Sqrt (Sub (Const 1.0, Pow (Var v, Const 2.0))))) when v = var ->
53
+
Some (Acos (Var var))
54
+
| Div (Const 1.0, Add (Const 1.0, Pow (Var v, Const 2.0))) when v = var ->
55
+
Some (Atan (Var var))
56
+
| Exp (Var v) when v = var ->
57
+
Some (Exp (Var var))
58
+
| Pow (Const a, Var v) when v = var ->
59
+
Some (simplify (Div (Pow (Const a, Var var), Ln (Const a))))
60
+
| Pow (SymConst E, Var v) when v = var ->
61
+
Some (Exp (Var var))
62
+
| e ->
63
+
match try_u_substitution var e with
64
+
| Some result -> Some result
65
+
| None -> try_by_parts var e
66
+
67
+
and try_u_substitution var expr =
68
+
let rec find_inner = function
69
+
| Sin u | Cos u | Tan u | Sinh u | Cosh u | Tanh u
70
+
| Asin u | Acos u | Atan u | Exp u | Ln u | Sqrt u | Abs u -> Some u
71
+
| Pow (u, _) -> Some u
72
+
| Add (e1, e2) | Sub (e1, e2) | Mul (e1, e2) | Div (e1, e2) ->
73
+
(match find_inner e1 with
74
+
| Some _ as r -> r
75
+
| None -> find_inner e2)
76
+
| _ -> None
77
+
in
78
+
match find_inner expr with
79
+
| Some u when u <> Var var ->
80
+
let u_prime = diff var u in
81
+
let expr_simplified = simplify expr in
82
+
let test_expr = simplify (Div (expr_simplified, u_prime)) in
83
+
let substituted = substitute var (Var "u_temp") test_expr in
84
+
(match substituted with
85
+
| e when not (contains_var var e) ->
86
+
(match integrate "u_temp" e with
87
+
| Some integrated ->
88
+
let result = substitute "u_temp" u integrated in
89
+
Some (simplify result)
90
+
| None -> None)
91
+
| _ -> None)
92
+
| _ -> None
93
+
94
+
and try_by_parts var = function
95
+
| Mul (e1, e2) ->
96
+
let priority = function
97
+
| Ln _ -> 5
98
+
| Asin _ | Acos _ | Atan _ -> 4
99
+
| Var _ | Pow (Var _, _) -> 3
100
+
| Sin _ | Cos _ | Tan _ | Sinh _ | Cosh _ | Tanh _ -> 2
101
+
| Exp _ -> 1
102
+
| _ -> 0
103
+
in
104
+
let (u, dv) =
105
+
if priority e1 >= priority e2 then (e1, e2) else (e2, e1)
106
+
in
107
+
(match integrate var dv with
108
+
| Some v ->
109
+
let du = diff var u in
110
+
(match integrate var (simplify (Mul (v, du))) with
111
+
| Some second_integral ->
112
+
Some (simplify (Sub (Mul (u, v), second_integral)))
113
+
| None -> None)
114
+
| None -> None)
115
+
| _ -> None
116
+
117
+
and contains_var var = function
118
+
| Const _ | SymConst _ -> false
119
+
| Var v -> v = var
120
+
| Add (e1, e2) | Sub (e1, e2) | Mul (e1, e2) | Div (e1, e2) | Pow (e1, e2) ->
121
+
contains_var var e1 || contains_var var e2
122
+
| Neg e | Sin e | Cos e | Tan e | Sinh e | Cosh e | Tanh e
123
+
| Asin e | Acos e | Atan e | Exp e | Ln e | Sqrt e | Abs e ->
124
+
contains_var var e
125
+
| Atan2 (e1, e2) | Log (e1, e2) ->
126
+
contains_var var e1 || contains_var var e2
127
+
128
+
let integrate_definite var lower upper expr =
129
+
match integrate var expr with
130
+
| None -> None
131
+
| Some antideriv ->
132
+
let upper_val = Eval.eval [(var, upper)] antideriv in
133
+
let lower_val = Eval.eval [(var, lower)] antideriv in
134
+
Some (upper_val -. lower_val)
+82
-14
lib/lexer.ml
+82
-14
lib/lexer.ml
···
1
+
type position = { line : int; col : int; offset : int }
2
+
1
3
type token =
2
4
| Num of float
3
5
| Var of string
4
6
| Ident of string
5
7
| Plus | Minus | Star | Slash | Caret
6
8
| LParen | RParen
9
+
| Comma
7
10
| EOF
8
11
12
+
type token_with_pos = { token : token; start_pos : position; end_pos : position }
13
+
9
14
let is_digit c = c >= '0' && c <= '9'
10
15
let is_alpha c = (c >= 'a' && c <= 'z') || (c >= 'A' && c <= 'Z')
11
-
let is_alphanum c = is_alpha c || is_digit c
16
+
let is_alphanum c = is_alpha c || is_digit c || c = '_'
17
+
18
+
let needs_implicit_mult tok1 tok2 =
19
+
match (tok1, tok2) with
20
+
| (Num _ | Var _ | RParen), (Num _ | Var _ | LParen | Ident _) -> true
21
+
| _ -> false
12
22
13
23
let tokenize str =
14
24
let len = String.length str in
25
+
let line = ref 1 in
26
+
let col = ref 1 in
27
+
15
28
let rec aux i acc =
16
-
if i >= len then List.rev (EOF :: acc)
29
+
if i >= len then
30
+
let pos = { line = !line; col = !col; offset = i } in
31
+
List.rev ({ token = EOF; start_pos = pos; end_pos = pos } :: acc)
17
32
else
33
+
let start_pos = { line = !line; col = !col; offset = i } in
18
34
match str.[i] with
19
-
| ' ' | '\t' | '\n' -> aux (i + 1) acc
20
-
| '+' -> aux (i + 1) (Plus :: acc)
21
-
| '-' -> aux (i + 1) (Minus :: acc)
22
-
| '*' -> aux (i + 1) (Star :: acc)
23
-
| '/' -> aux (i + 1) (Slash :: acc)
24
-
| '^' -> aux (i + 1) (Caret :: acc)
25
-
| '(' -> aux (i + 1) (LParen :: acc)
26
-
| ')' -> aux (i + 1) (RParen :: acc)
35
+
| ' ' | '\t' ->
36
+
col := !col + 1;
37
+
aux (i + 1) acc
38
+
| '\n' ->
39
+
line := !line + 1;
40
+
col := 1;
41
+
aux (i + 1) acc
42
+
| '+' ->
43
+
col := !col + 1;
44
+
let end_pos = { line = !line; col = !col; offset = i + 1 } in
45
+
aux (i + 1) ({ token = Plus; start_pos; end_pos } :: acc)
46
+
| '-' ->
47
+
col := !col + 1;
48
+
let end_pos = { line = !line; col = !col; offset = i + 1 } in
49
+
aux (i + 1) ({ token = Minus; start_pos; end_pos } :: acc)
50
+
| '*' ->
51
+
col := !col + 1;
52
+
let end_pos = { line = !line; col = !col; offset = i + 1 } in
53
+
aux (i + 1) ({ token = Star; start_pos; end_pos } :: acc)
54
+
| '/' ->
55
+
col := !col + 1;
56
+
let end_pos = { line = !line; col = !col; offset = i + 1 } in
57
+
aux (i + 1) ({ token = Slash; start_pos; end_pos } :: acc)
58
+
| '^' ->
59
+
col := !col + 1;
60
+
let end_pos = { line = !line; col = !col; offset = i + 1 } in
61
+
aux (i + 1) ({ token = Caret; start_pos; end_pos } :: acc)
62
+
| '(' ->
63
+
col := !col + 1;
64
+
let end_pos = { line = !line; col = !col; offset = i + 1 } in
65
+
aux (i + 1) ({ token = LParen; start_pos; end_pos } :: acc)
66
+
| ')' ->
67
+
col := !col + 1;
68
+
let end_pos = { line = !line; col = !col; offset = i + 1 } in
69
+
aux (i + 1) ({ token = RParen; start_pos; end_pos } :: acc)
70
+
| ',' ->
71
+
col := !col + 1;
72
+
let end_pos = { line = !line; col = !col; offset = i + 1 } in
73
+
aux (i + 1) ({ token = Comma; start_pos; end_pos } :: acc)
27
74
| c when is_digit c || c = '.' ->
28
75
let j = ref (i + 1) in
29
76
while !j < len && (is_digit str.[!j] || str.[!j] = '.') do
30
77
incr j
31
78
done;
32
79
let num_str = String.sub str i (!j - i) in
33
-
aux !j (Num (float_of_string num_str) :: acc)
80
+
col := !col + (!j - i);
81
+
let end_pos = { line = !line; col = !col; offset = !j } in
82
+
aux !j ({ token = Num (float_of_string num_str); start_pos; end_pos } :: acc)
34
83
| c when is_alpha c ->
35
84
let j = ref (i + 1) in
36
85
while !j < len && is_alphanum str.[!j] do
37
86
incr j
38
87
done;
39
88
let id = String.sub str i (!j - i) in
89
+
col := !col + (!j - i);
90
+
let end_pos = { line = !line; col = !col; offset = !j } in
40
91
let tok = if !j < len && str.[!j] = '(' then Ident id else Var id in
41
-
aux !j (tok :: acc)
42
-
| c -> failwith (Printf.sprintf "unexpected character: %c" c)
92
+
aux !j ({ token = tok; start_pos; end_pos } :: acc)
93
+
| c ->
94
+
failwith (Printf.sprintf "unexpected character: %c at line %d, column %d"
95
+
c !line !col)
43
96
in
44
-
aux 0 []
97
+
98
+
let tokens = aux 0 [] in
99
+
100
+
let rec insert_implicit_mult = function
101
+
| [] -> []
102
+
| [t] -> [t]
103
+
| t1 :: t2 :: rest ->
104
+
if needs_implicit_mult t1.token t2.token then
105
+
let mult_pos = t2.start_pos in
106
+
t1 :: { token = Star; start_pos = mult_pos; end_pos = mult_pos } ::
107
+
insert_implicit_mult (t2 :: rest)
108
+
else
109
+
t1 :: insert_implicit_mult (t2 :: rest)
110
+
in
111
+
112
+
insert_implicit_mult tokens
+70
lib/limits.ml
+70
lib/limits.ml
···
1
+
open Expr
2
+
open Simplify
3
+
open Diff
4
+
open Substitute
5
+
6
+
type direction = FromLeft | FromRight | Bidirectional
7
+
8
+
let rec limit expr var point direction =
9
+
let try_direct_sub () =
10
+
try
11
+
let substituted = substitute var point expr in
12
+
let simplified = simplify substituted in
13
+
match simplified with
14
+
| Const c when not (Float.is_nan c || Float.is_infinite c) -> Some simplified
15
+
| _ -> None
16
+
with _ -> None
17
+
in
18
+
19
+
let try_lhopital () =
20
+
let numerator, denominator = match expr with
21
+
| Div (n, d) -> (n, d)
22
+
| _ -> (expr, Const 1.0)
23
+
in
24
+
25
+
let num_at_point = substitute var point numerator |> simplify in
26
+
let den_at_point = substitute var point denominator |> simplify in
27
+
28
+
match (num_at_point, den_at_point) with
29
+
| (Const 0.0, Const 0.0) ->
30
+
let num_deriv = diff var numerator in
31
+
let den_deriv = diff var denominator in
32
+
let new_expr = Div (num_deriv, den_deriv) in
33
+
limit new_expr var point direction
34
+
| _ -> None
35
+
in
36
+
37
+
let try_series_expansion () =
38
+
None
39
+
in
40
+
41
+
match try_direct_sub () with
42
+
| Some result -> Some result
43
+
| None ->
44
+
match try_lhopital () with
45
+
| Some result -> Some result
46
+
| None -> try_series_expansion ()
47
+
48
+
let limit_at_infinity expr var =
49
+
let rec find_highest_power = function
50
+
| Pow (Var v, Const n) when v = var -> Some (int_of_float n)
51
+
| Div (num, den) ->
52
+
let num_pow = find_highest_power num |> Option.value ~default:0 in
53
+
let den_pow = find_highest_power den |> Option.value ~default:0 in
54
+
Some (num_pow - den_pow)
55
+
| Add (e1, e2) | Sub (e1, e2) ->
56
+
let p1 = find_highest_power e1 |> Option.value ~default:0 in
57
+
let p2 = find_highest_power e2 |> Option.value ~default:0 in
58
+
Some (max p1 p2)
59
+
| Mul (e1, e2) ->
60
+
let p1 = find_highest_power e1 |> Option.value ~default:0 in
61
+
let p2 = find_highest_power e2 |> Option.value ~default:0 in
62
+
Some (p1 + p2)
63
+
| _ -> None
64
+
in
65
+
66
+
match find_highest_power expr with
67
+
| Some p when p > 0 -> Some (SymConst E)
68
+
| Some p when p < 0 -> Some (Const 0.0)
69
+
| Some _ -> Some (Const 1.0)
70
+
| None -> None
+137
lib/matrix.ml
+137
lib/matrix.ml
···
1
+
open Expr
2
+
open Simplify
3
+
4
+
type matrix = expr array array
5
+
6
+
let create rows cols init =
7
+
Array.init rows (fun i -> Array.init cols (fun j -> init i j))
8
+
9
+
let identity n =
10
+
create n n (fun i j -> if i = j then Const 1.0 else Const 0.0)
11
+
12
+
let rows m = Array.length m
13
+
let cols m = if Array.length m > 0 then Array.length m.(0) else 0
14
+
15
+
let get m i j = m.(i).(j)
16
+
let set m i j v = m.(i).(j) <- v
17
+
18
+
let map f m =
19
+
Array.map (Array.map f) m
20
+
21
+
let transpose m =
22
+
let r = rows m in
23
+
let c = cols m in
24
+
create c r (fun i j -> m.(j).(i))
25
+
26
+
let add m1 m2 =
27
+
if rows m1 <> rows m2 || cols m1 <> cols m2 then
28
+
failwith "matrix dimensions must match for addition"
29
+
else
30
+
create (rows m1) (cols m1) (fun i j ->
31
+
simplify (Add (m1.(i).(j), m2.(i).(j)))
32
+
)
33
+
34
+
let mult m1 m2 =
35
+
if cols m1 <> rows m2 then
36
+
failwith "incompatible matrix dimensions for multiplication"
37
+
else
38
+
create (rows m1) (cols m2) (fun i j ->
39
+
let sum = ref (Const 0.0) in
40
+
for k = 0 to cols m1 - 1 do
41
+
sum := Add (!sum, Mul (m1.(i).(k), m2.(k).(j)))
42
+
done;
43
+
simplify !sum
44
+
)
45
+
46
+
let scalar_mult s m =
47
+
map (fun e -> simplify (Mul (s, e))) m
48
+
49
+
let det m =
50
+
let n = rows m in
51
+
if n <> cols m then failwith "determinant requires square matrix";
52
+
53
+
let rec determinant mat size =
54
+
if size = 1 then mat.(0).(0)
55
+
else if size = 2 then
56
+
simplify (Sub (Mul (mat.(0).(0), mat.(1).(1)),
57
+
Mul (mat.(0).(1), mat.(1).(0))))
58
+
else
59
+
let result = ref (Const 0.0) in
60
+
for j = 0 to size - 1 do
61
+
let minor = create (size - 1) (size - 1) (fun i k ->
62
+
let mi = if i < 0 then i else i + 1 in
63
+
let mk = if k < j then k else k + 1 in
64
+
mat.(mi).(mk)
65
+
) in
66
+
let cofactor = determinant minor (size - 1) in
67
+
let sign = if j mod 2 = 0 then Const 1.0 else Const (-1.0) in
68
+
result := Add (!result, Mul (Mul (sign, mat.(0).(j)), cofactor))
69
+
done;
70
+
simplify !result
71
+
in
72
+
determinant m n
73
+
74
+
let inverse m =
75
+
let n = rows m in
76
+
if n <> cols m then None
77
+
else
78
+
let d = det m in
79
+
match d with
80
+
| Const 0.0 -> None
81
+
| _ ->
82
+
let adj = create n n (fun i j ->
83
+
let minor = create (n - 1) (n - 1) (fun mi mj ->
84
+
let si = if mi < i then mi else mi + 1 in
85
+
let sj = if mj < j then mj else mj + 1 in
86
+
m.(si).(sj)
87
+
) in
88
+
let minor_det = det minor in
89
+
let sign = if (i + j) mod 2 = 0 then Const 1.0 else Const (-1.0) in
90
+
simplify (Mul (sign, minor_det))
91
+
) in
92
+
let adj_t = transpose adj in
93
+
Some (map (fun e -> simplify (Div (e, d))) adj_t)
94
+
95
+
let trace m =
96
+
let n = min (rows m) (cols m) in
97
+
let sum = ref (Const 0.0) in
98
+
for i = 0 to n - 1 do
99
+
sum := Add (!sum, m.(i).(i))
100
+
done;
101
+
simplify !sum
102
+
103
+
let eigenvalues _m =
104
+
[]
105
+
106
+
let rank m =
107
+
let r = rows m in
108
+
let c = cols m in
109
+
let temp = Array.map Array.copy m in
110
+
let rec count_pivots row col rank =
111
+
if row >= r || col >= c then rank
112
+
else
113
+
match temp.(row).(col) with
114
+
| Const 0.0 ->
115
+
let rec find_pivot i =
116
+
if i >= r then None
117
+
else match temp.(i).(col) with
118
+
| Const 0.0 -> find_pivot (i + 1)
119
+
| _ -> Some i
120
+
in
121
+
(match find_pivot (row + 1) with
122
+
| None -> count_pivots row (col + 1) rank
123
+
| Some i ->
124
+
let tmp_row = temp.(row) in
125
+
temp.(row) <- temp.(i);
126
+
temp.(i) <- tmp_row;
127
+
count_pivots row col rank)
128
+
| pivot ->
129
+
for i = row + 1 to r - 1 do
130
+
let factor = Div (temp.(i).(col), pivot) in
131
+
for j = col to c - 1 do
132
+
temp.(i).(j) <- simplify (Sub (temp.(i).(j), Mul (factor, temp.(row).(j))))
133
+
done
134
+
done;
135
+
count_pivots (row + 1) (col + 1) (rank + 1)
136
+
in
137
+
count_pivots 0 0 0
+41
lib/multivariate.ml
+41
lib/multivariate.ml
···
1
+
open Expr
2
+
open Diff
3
+
open Simplify
4
+
5
+
let gradient vars expr =
6
+
List.map (fun v -> diff v expr) vars
7
+
8
+
let jacobian vars exprs =
9
+
List.map (fun expr -> gradient vars expr) exprs
10
+
11
+
let hessian vars expr =
12
+
let grad = gradient vars expr in
13
+
List.map (fun deriv -> List.map (fun v -> diff v deriv) vars) grad
14
+
15
+
let divergence vars vector_field =
16
+
if List.length vars <> List.length vector_field then
17
+
failwith "divergence: dimension mismatch"
18
+
else
19
+
let terms = List.map2 (fun v e -> diff v e) vars vector_field in
20
+
let sum = List.fold_left (fun acc t -> Add (acc, t)) (Const 0.0) terms in
21
+
simplify sum
22
+
23
+
let curl vector_field =
24
+
match vector_field with
25
+
| [fx; fy; fz] ->
26
+
let dfz_dy = diff "y" fz in
27
+
let dfy_dz = diff "z" fy in
28
+
let dfx_dz = diff "z" fx in
29
+
let dfz_dx = diff "x" fz in
30
+
let dfy_dx = diff "x" fy in
31
+
let dfx_dy = diff "y" fx in
32
+
[simplify (Sub (dfz_dy, dfy_dz));
33
+
simplify (Sub (dfx_dz, dfz_dx));
34
+
simplify (Sub (dfy_dx, dfx_dy))]
35
+
| _ -> failwith "curl: requires exactly 3 components"
36
+
37
+
let laplacian vars expr =
38
+
let hess = hessian vars expr in
39
+
let diagonal = List.mapi (fun i row -> List.nth row i) hess in
40
+
let sum = List.fold_left (fun acc t -> Add (acc, t)) (Const 0.0) diagonal in
41
+
simplify sum
+154
lib/numerical.ml
+154
lib/numerical.ml
···
1
+
open Diff
2
+
open Eval
3
+
open Multivariate
4
+
5
+
let newton_raphson expr var initial tolerance max_iter =
6
+
let f_prime = diff var expr in
7
+
let rec iterate x n =
8
+
if n >= max_iter then None
9
+
else
10
+
let fx = eval [(var, x)] expr in
11
+
if abs_float fx < tolerance then Some x
12
+
else
13
+
let fpx = eval [(var, x)] f_prime in
14
+
if abs_float fpx < 1e-10 then None
15
+
else
16
+
let x_next = x -. fx /. fpx in
17
+
if abs_float (x_next -. x) < tolerance then Some x_next
18
+
else iterate x_next (n + 1)
19
+
in
20
+
iterate initial 0
21
+
22
+
let bisection expr var left right tolerance max_iter =
23
+
let rec iterate a b n =
24
+
if n >= max_iter then None
25
+
else
26
+
let fa = eval [(var, a)] expr in
27
+
let fb = eval [(var, b)] expr in
28
+
if fa *. fb > 0.0 then None
29
+
else
30
+
let mid = (a +. b) /. 2.0 in
31
+
let fmid = eval [(var, mid)] expr in
32
+
if abs_float fmid < tolerance || abs_float (b -. a) < tolerance then
33
+
Some mid
34
+
else if fa *. fmid < 0.0 then
35
+
iterate a mid (n + 1)
36
+
else
37
+
iterate mid b (n + 1)
38
+
in
39
+
iterate left right 0
40
+
41
+
let trapezoidal expr var lower upper n =
42
+
let h = (upper -. lower) /. float_of_int n in
43
+
let rec sum_interior i acc =
44
+
if i >= n then acc
45
+
else
46
+
let x = lower +. float_of_int i *. h in
47
+
let fx = eval [(var, x)] expr in
48
+
sum_interior (i + 1) (acc +. fx)
49
+
in
50
+
let f_lower = eval [(var, lower)] expr in
51
+
let f_upper = eval [(var, upper)] expr in
52
+
let interior = sum_interior 1 0.0 in
53
+
h *. (f_lower /. 2.0 +. interior +. f_upper /. 2.0)
54
+
55
+
let simpsons expr var lower upper n =
56
+
let n = if n mod 2 = 1 then n + 1 else n in
57
+
let h = (upper -. lower) /. float_of_int n in
58
+
let rec sum_terms i acc_odd acc_even =
59
+
if i >= n then (acc_odd, acc_even)
60
+
else
61
+
let x = lower +. float_of_int i *. h in
62
+
let fx = eval [(var, x)] expr in
63
+
if i mod 2 = 1 then
64
+
sum_terms (i + 1) (acc_odd +. fx) acc_even
65
+
else if i > 0 then
66
+
sum_terms (i + 1) acc_odd (acc_even +. fx)
67
+
else
68
+
sum_terms (i + 1) acc_odd acc_even
69
+
in
70
+
let f_lower = eval [(var, lower)] expr in
71
+
let f_upper = eval [(var, upper)] expr in
72
+
let (odd, even) = sum_terms 1 0.0 0.0 in
73
+
h /. 3.0 *. (f_lower +. 4.0 *. odd +. 2.0 *. even +. f_upper)
74
+
75
+
let rec adaptive_quadrature expr var lower upper tolerance =
76
+
let mid = (lower +. upper) /. 2.0 in
77
+
let whole = simpsons expr var lower upper 10 in
78
+
let left_half = simpsons expr var lower mid 10 in
79
+
let right_half = simpsons expr var mid upper 10 in
80
+
let error = abs_float (whole -. (left_half +. right_half)) in
81
+
if error < tolerance then
82
+
left_half +. right_half
83
+
else
84
+
let left = adaptive_quadrature expr var lower mid (tolerance /. 2.0) in
85
+
let right = adaptive_quadrature expr var mid upper (tolerance /. 2.0) in
86
+
left +. right
87
+
88
+
let gradient_descent expr vars initial learning_rate max_iter =
89
+
let grad_exprs = gradient vars expr in
90
+
let rec iterate point n =
91
+
if n >= max_iter then Some point
92
+
else
93
+
let env = List.combine vars point in
94
+
let grad_vals = List.map (eval env) grad_exprs in
95
+
let new_point = List.map2 (fun p g -> p -. learning_rate *. g) point grad_vals in
96
+
let diff = List.map2 (fun a b -> abs_float (a -. b)) new_point point in
97
+
let max_diff = List.fold_left max 0.0 diff in
98
+
if max_diff < 1e-6 then Some new_point
99
+
else iterate new_point (n + 1)
100
+
in
101
+
iterate initial 0
102
+
103
+
let invert_matrix matrix =
104
+
let n = List.length matrix in
105
+
let augmented = List.mapi (fun i row ->
106
+
row @ List.init n (fun j -> if i = j then 1.0 else 0.0)
107
+
) matrix in
108
+
109
+
let rec gaussian_elimination mat row =
110
+
if row >= n then mat
111
+
else
112
+
let pivot_row = List.nth mat row in
113
+
let pivot = List.nth pivot_row row in
114
+
if abs_float pivot < 1e-10 then mat
115
+
else
116
+
let normalized = List.map (fun x -> x /. pivot) pivot_row in
117
+
let updated = List.mapi (fun i r ->
118
+
if i = row then normalized
119
+
else
120
+
let factor = List.nth r row in
121
+
List.map2 (fun a b -> a -. factor *. b) r normalized
122
+
) mat in
123
+
gaussian_elimination updated (row + 1)
124
+
in
125
+
126
+
let reduced = gaussian_elimination augmented 0 in
127
+
List.map (fun row -> List.filteri (fun i _ -> i >= n) row) reduced
128
+
129
+
let newtons_method_opt expr vars initial tolerance max_iter =
130
+
let grad_exprs = gradient vars expr in
131
+
let hess_matrix = hessian vars expr in
132
+
133
+
let rec iterate point n =
134
+
if n >= max_iter then None
135
+
else
136
+
let env = List.combine vars point in
137
+
let grad_vals = List.map (eval env) grad_exprs in
138
+
let hess_vals = List.map (fun row ->
139
+
List.map (eval env) row
140
+
) hess_matrix in
141
+
142
+
let hess_inv = invert_matrix hess_vals in
143
+
let delta = List.map (fun row ->
144
+
List.fold_left2 (fun acc h g -> acc +. h *. g) 0.0 row grad_vals
145
+
) hess_inv in
146
+
147
+
let new_point = List.map2 (fun p d -> p -. d) point delta in
148
+
let diff = List.map2 (fun a b -> abs_float (a -. b)) new_point point in
149
+
let max_diff = List.fold_left max 0.0 diff in
150
+
151
+
if max_diff < tolerance then Some new_point
152
+
else iterate new_point (n + 1)
153
+
in
154
+
iterate initial 0
+78
lib/ode.ml
+78
lib/ode.ml
···
1
+
open Expr
2
+
open Simplify
3
+
open Integrate
4
+
5
+
type ode = {
6
+
equation: expr;
7
+
dependent: string;
8
+
independent: string;
9
+
}
10
+
11
+
let dsolve ode_eq dep_var indep_var =
12
+
let try_separation_of_variables () =
13
+
match ode_eq with
14
+
| Mul (f_y, g_x) ->
15
+
(match (integrate indep_var g_x, integrate dep_var f_y) with
16
+
| Some int_g, Some int_f ->
17
+
Some (simplify (Sub (int_f, int_g)))
18
+
| _ -> None)
19
+
| _ -> None
20
+
in
21
+
22
+
let try_linear_first_order () =
23
+
None
24
+
in
25
+
26
+
let try_exact_equation () =
27
+
None
28
+
in
29
+
30
+
let try_integrating_factor () =
31
+
None
32
+
in
33
+
34
+
match try_separation_of_variables () with
35
+
| Some sol -> Some sol
36
+
| None ->
37
+
match try_linear_first_order () with
38
+
| Some sol -> Some sol
39
+
| None ->
40
+
match try_exact_equation () with
41
+
| Some sol -> Some sol
42
+
| None -> try_integrating_factor ()
43
+
44
+
let solve_second_order coeff_y'' coeff_y' coeff_y rhs var =
45
+
let a = coeff_y'' in
46
+
let b = coeff_y' in
47
+
let c = coeff_y in
48
+
49
+
match (a, b, c) with
50
+
| (Const a_val, Const b_val, Const c_val) when rhs = Const 0.0 ->
51
+
let discriminant = b_val *. b_val -. 4.0 *. a_val *. c_val in
52
+
if discriminant > 0.0 then
53
+
let r1 = (-.b_val +. sqrt discriminant) /. (2.0 *. a_val) in
54
+
let r2 = (-.b_val -. sqrt discriminant) /. (2.0 *. a_val) in
55
+
Some (Add (
56
+
Mul (Var "C1", Exp (Mul (Const r1, Var var))),
57
+
Mul (Var "C2", Exp (Mul (Const r2, Var var)))
58
+
))
59
+
else if discriminant = 0.0 then
60
+
let r = -.b_val /. (2.0 *. a_val) in
61
+
Some (Mul (
62
+
Add (Var "C1", Mul (Var "C2", Var var)),
63
+
Exp (Mul (Const r, Var var))
64
+
))
65
+
else
66
+
let real_part = -.b_val /. (2.0 *. a_val) in
67
+
let imag_part = sqrt (-.discriminant) /. (2.0 *. a_val) in
68
+
Some (Mul (
69
+
Exp (Mul (Const real_part, Var var)),
70
+
Add (
71
+
Mul (Var "C1", Cos (Mul (Const imag_part, Var var))),
72
+
Mul (Var "C2", Sin (Mul (Const imag_part, Var var)))
73
+
)
74
+
))
75
+
| _ -> None
76
+
77
+
let solve_system _odes _indep_var =
78
+
[]
+126
-14
lib/parser.ml
+126
-14
lib/parser.ml
···
1
1
open Lexer
2
2
open Expr
3
3
4
-
type parse_state = { tokens : token list; mutable pos : int }
4
+
exception Parse_error of string * position
5
+
6
+
type parse_state = { tokens : token_with_pos list; mutable pos : int }
7
+
8
+
let token_name = function
9
+
| Num _ -> "number"
10
+
| Var _ -> "variable"
11
+
| Ident _ -> "identifier"
12
+
| Plus -> "'+'"
13
+
| Minus -> "'-'"
14
+
| Star -> "'*'"
15
+
| Slash -> "'/'"
16
+
| Caret -> "'^'"
17
+
| LParen -> "'('"
18
+
| RParen -> "')'"
19
+
| Comma -> "','"
20
+
| EOF -> "end of input"
5
21
6
22
let peek state =
7
23
if state.pos < List.length state.tokens then
8
-
List.nth state.tokens state.pos
24
+
(List.nth state.tokens state.pos).token
9
25
else EOF
10
26
27
+
let peek_pos state =
28
+
if state.pos < List.length state.tokens then
29
+
(List.nth state.tokens state.pos).start_pos
30
+
else
31
+
{ line = 0; col = 0; offset = 0 }
32
+
11
33
let advance state =
12
34
state.pos <- state.pos + 1
13
35
14
36
let expect state tok =
15
37
if peek state = tok then advance state
16
-
else failwith (Printf.sprintf "expected token")
38
+
else
39
+
let pos = peek_pos state in
40
+
raise (Parse_error (
41
+
Printf.sprintf "expected %s but found %s at line %d, column %d"
42
+
(token_name tok) (token_name (peek state)) pos.line pos.col,
43
+
pos))
17
44
18
45
let rec parse_expr state =
19
46
parse_additive state
···
62
89
Const n
63
90
| Var v ->
64
91
advance state;
65
-
Var v
92
+
(match v with
93
+
| "pi" | "π" -> SymConst Pi
94
+
| "e" -> SymConst E
95
+
| _ -> Var v)
66
96
| Ident "sin" ->
67
97
advance state;
68
98
expect state LParen;
···
81
111
let arg = parse_expr state in
82
112
expect state RParen;
83
113
Tan arg
114
+
| Ident "sinh" ->
115
+
advance state;
116
+
expect state LParen;
117
+
let arg = parse_expr state in
118
+
expect state RParen;
119
+
Sinh arg
120
+
| Ident "cosh" ->
121
+
advance state;
122
+
expect state LParen;
123
+
let arg = parse_expr state in
124
+
expect state RParen;
125
+
Cosh arg
126
+
| Ident "tanh" ->
127
+
advance state;
128
+
expect state LParen;
129
+
let arg = parse_expr state in
130
+
expect state RParen;
131
+
Tanh arg
132
+
| Ident "asin" | Ident "arcsin" ->
133
+
advance state;
134
+
expect state LParen;
135
+
let arg = parse_expr state in
136
+
expect state RParen;
137
+
Asin arg
138
+
| Ident "acos" | Ident "arccos" ->
139
+
advance state;
140
+
expect state LParen;
141
+
let arg = parse_expr state in
142
+
expect state RParen;
143
+
Acos arg
144
+
| Ident "atan" | Ident "arctan" ->
145
+
advance state;
146
+
expect state LParen;
147
+
let arg = parse_expr state in
148
+
expect state RParen;
149
+
Atan arg
150
+
| Ident "atan2" ->
151
+
advance state;
152
+
expect state LParen;
153
+
let arg1 = parse_expr state in
154
+
expect state Comma;
155
+
let arg2 = parse_expr state in
156
+
expect state RParen;
157
+
Atan2 (arg1, arg2)
84
158
| Ident "exp" ->
85
159
advance state;
86
160
expect state LParen;
87
161
let arg = parse_expr state in
88
162
expect state RParen;
89
163
Exp arg
90
-
| Ident "ln" ->
164
+
| Ident "ln" | Ident "log" ->
165
+
advance state;
166
+
expect state LParen;
167
+
let arg1 = parse_expr state in
168
+
(match peek state with
169
+
| Comma ->
170
+
advance state;
171
+
let arg2 = parse_expr state in
172
+
expect state RParen;
173
+
Log (arg1, arg2)
174
+
| RParen ->
175
+
advance state;
176
+
Ln arg1
177
+
| _ ->
178
+
let pos = peek_pos state in
179
+
raise (Parse_error (
180
+
Printf.sprintf "expected ',' or ')' in log function at line %d, column %d"
181
+
pos.line pos.col,
182
+
pos)))
183
+
| Ident "sqrt" ->
91
184
advance state;
92
185
expect state LParen;
93
186
let arg = parse_expr state in
94
187
expect state RParen;
95
-
Ln arg
96
-
| Ident "e" ->
188
+
Sqrt arg
189
+
| Ident "abs" ->
97
190
advance state;
98
-
Const (exp 1.0)
191
+
expect state LParen;
192
+
let arg = parse_expr state in
193
+
expect state RParen;
194
+
Abs arg
99
195
| LParen ->
100
196
advance state;
101
197
let expr = parse_expr state in
102
198
expect state RParen;
103
199
expr
104
-
| _ -> failwith "unexpected token in primary"
200
+
| _ ->
201
+
let pos = peek_pos state in
202
+
raise (Parse_error (
203
+
Printf.sprintf "unexpected %s at line %d, column %d"
204
+
(token_name (peek state)) pos.line pos.col,
205
+
pos))
105
206
106
207
let parse str =
107
-
let tokens = tokenize str in
108
-
let state = { tokens; pos = 0 } in
109
-
let expr = parse_expr state in
110
-
if peek state = EOF then expr
111
-
else failwith "unexpected tokens after expression"
208
+
try
209
+
let tokens = tokenize str in
210
+
let state = { tokens; pos = 0 } in
211
+
let expr = parse_expr state in
212
+
if peek state = EOF then expr
213
+
else
214
+
let pos = peek_pos state in
215
+
raise (Parse_error (
216
+
Printf.sprintf "unexpected tokens after expression at line %d, column %d"
217
+
pos.line pos.col,
218
+
pos))
219
+
with
220
+
| Parse_error (msg, _) ->
221
+
failwith msg
222
+
| Failure msg ->
223
+
failwith msg
+58
lib/piecewise.ml
+58
lib/piecewise.ml
···
1
+
open Expr
2
+
open Diff
3
+
open Eval
4
+
5
+
type condition = expr
6
+
type piece = condition * expr
7
+
type piecewise = piece list * expr
8
+
9
+
let create pieces default =
10
+
(pieces, default)
11
+
12
+
let eval_piecewise env (pieces, default) =
13
+
let rec find_first = function
14
+
| [] -> eval env default
15
+
| (cond, value) :: rest ->
16
+
let cond_val = eval env cond in
17
+
if cond_val <> 0.0 then eval env value
18
+
else find_first rest
19
+
in
20
+
find_first pieces
21
+
22
+
let diff_piecewise var (pieces, default) =
23
+
let diff_pieces = List.map (fun (cond, expr) ->
24
+
(cond, diff var expr)
25
+
) pieces in
26
+
(diff_pieces, diff var default)
27
+
28
+
let to_expr (_pieces, default) =
29
+
default
30
+
31
+
let heaviside x =
32
+
create [(Add (x, Const 0.0), Const 1.0)] (Const 0.0)
33
+
34
+
let signum x =
35
+
create [
36
+
(Add (x, Const 0.0), Const 1.0);
37
+
(Sub (Const 0.0, x), Const (-1.0))
38
+
] (Const 0.0)
39
+
40
+
let abs_as_piecewise x =
41
+
create [
42
+
(Add (x, Const 0.0), x)
43
+
] (Neg x)
44
+
45
+
let max_expr e1 e2 =
46
+
create [(Sub (e1, e2), e1)] e2
47
+
48
+
let min_expr e1 e2 =
49
+
create [(Sub (e2, e1), e1)] e2
50
+
51
+
let continuous_extension expr _var _point =
52
+
expr
53
+
54
+
let piecewise_to_string (pieces, default) =
55
+
let piece_strs = List.map (fun (cond, value) ->
56
+
Printf.sprintf " %s if %s" (to_string value) (to_string cond)
57
+
) pieces in
58
+
String.concat "\n" piece_strs ^ "\n " ^ to_string default ^ " otherwise"
+117
lib/polynomial.ml
+117
lib/polynomial.ml
···
1
+
open Expr
2
+
open Simplify
3
+
4
+
type poly_term = {
5
+
coeff: float;
6
+
powers: (string * int) list;
7
+
}
8
+
9
+
type polynomial = poly_term list
10
+
11
+
let rec expr_to_poly expr vars =
12
+
match expr with
13
+
| Const c -> [{coeff = c; powers = []}]
14
+
| Var v when List.mem v vars -> [{coeff = 1.0; powers = [(v, 1)]}]
15
+
| Var _ -> [{coeff = 1.0; powers = []}]
16
+
| Mul (Const c, e) | Mul (e, Const c) ->
17
+
List.map (fun t -> {t with coeff = t.coeff *. c}) (expr_to_poly e vars)
18
+
| Add (e1, e2) ->
19
+
(expr_to_poly e1 vars) @ (expr_to_poly e2 vars)
20
+
| Mul (e1, e2) ->
21
+
let p1 = expr_to_poly e1 vars in
22
+
let p2 = expr_to_poly e2 vars in
23
+
List.concat_map (fun t1 ->
24
+
List.map (fun t2 ->
25
+
let new_coeff = t1.coeff *. t2.coeff in
26
+
let new_powers = combine_powers t1.powers t2.powers in
27
+
{coeff = new_coeff; powers = new_powers}
28
+
) p2
29
+
) p1
30
+
| Pow (Var v, Const n) when List.mem v vars && Float.is_integer n ->
31
+
[{coeff = 1.0; powers = [(v, int_of_float n)]}]
32
+
| _ -> [{coeff = 1.0; powers = []}]
33
+
34
+
and combine_powers p1 p2 =
35
+
let rec merge acc = function
36
+
| [] -> List.rev acc
37
+
| (v, n) :: rest ->
38
+
match List.assoc_opt v acc with
39
+
| Some m -> merge ((v, m + n) :: List.remove_assoc v acc) rest
40
+
| None -> merge ((v, n) :: acc) rest
41
+
in
42
+
merge p1 p2
43
+
44
+
let poly_to_expr poly =
45
+
let term_to_expr t =
46
+
let coeff_expr = if t.coeff = 1.0 && t.powers <> [] then None
47
+
else Some (Const t.coeff) in
48
+
let vars_expr = List.fold_left (fun acc (v, n) ->
49
+
let var_power = if n = 1 then Var v
50
+
else Pow (Var v, Const (float_of_int n)) in
51
+
match acc with
52
+
| None -> Some var_power
53
+
| Some e -> Some (Mul (e, var_power))
54
+
) None t.powers in
55
+
match (coeff_expr, vars_expr) with
56
+
| None, None -> Const 1.0
57
+
| Some c, None -> c
58
+
| None, Some v -> v
59
+
| Some c, Some v -> Mul (c, v)
60
+
in
61
+
match poly with
62
+
| [] -> Const 0.0
63
+
| [t] -> term_to_expr t
64
+
| t :: rest ->
65
+
List.fold_left (fun acc term ->
66
+
Add (acc, term_to_expr term)
67
+
) (term_to_expr t) rest
68
+
69
+
let expand expr =
70
+
let rec expand_expr = function
71
+
| Add (e1, e2) -> simplify (Add (expand_expr e1, expand_expr e2))
72
+
| Sub (e1, e2) -> simplify (Sub (expand_expr e1, expand_expr e2))
73
+
| Mul (Add (a, b), c) | Mul (c, Add (a, b)) ->
74
+
simplify (Add (Mul (expand_expr a, expand_expr c),
75
+
Mul (expand_expr b, expand_expr c)))
76
+
| Mul (Sub (a, b), c) | Mul (c, Sub (a, b)) ->
77
+
simplify (Sub (Mul (expand_expr a, expand_expr c),
78
+
Mul (expand_expr b, expand_expr c)))
79
+
| Mul (e1, e2) -> simplify (Mul (expand_expr e1, expand_expr e2))
80
+
| Pow (Add (a, b), Const n) when Float.is_integer n && n >= 0.0 && n <= 10.0 ->
81
+
let rec binomial_expand base pow =
82
+
if pow = 0.0 then Const 1.0
83
+
else Mul (expand_expr base, binomial_expand base (pow -. 1.0))
84
+
in
85
+
expand_expr (binomial_expand (Add (a, b)) n)
86
+
| Pow (e, n) -> Pow (expand_expr e, n)
87
+
| e -> e
88
+
in
89
+
let expanded = expand_expr expr in
90
+
simplify expanded
91
+
92
+
let rec poly_gcd p1 p2 =
93
+
if List.length p2 = 0 then p1
94
+
else
95
+
let remainder = poly_remainder p1 p2 in
96
+
poly_gcd p2 remainder
97
+
98
+
and poly_remainder _p1 _p2 =
99
+
[]
100
+
101
+
let degree expr var =
102
+
let rec find_degree = function
103
+
| Var v when v = var -> 1
104
+
| Pow (Var v, Const n) when v = var && Float.is_integer n -> int_of_float n
105
+
| Mul (e1, e2) -> find_degree e1 + find_degree e2
106
+
| Add (e1, e2) -> max (find_degree e1) (find_degree e2)
107
+
| Sub (e1, e2) -> max (find_degree e1) (find_degree e2)
108
+
| _ -> 0
109
+
in
110
+
find_degree expr
111
+
112
+
let collect expr _var =
113
+
let expanded = expand expr in
114
+
simplify expanded
115
+
116
+
let factor expr =
117
+
[expr]
+26
lib/series.ml
+26
lib/series.ml
···
1
+
open Expr
2
+
open Diff
3
+
open Substitute
4
+
open Simplify
5
+
6
+
let rec factorial n =
7
+
if n <= 1 then 1.0
8
+
else float_of_int n *. factorial (n - 1)
9
+
10
+
let taylor var expr center order =
11
+
let rec build_terms n acc =
12
+
if n > order then acc
13
+
else
14
+
let derivative = diff_n var n expr in
15
+
let deriv_at_center = substitute var center derivative in
16
+
let coeff = Div (deriv_at_center, Const (factorial n)) in
17
+
let h = Sub (Var var, center) in
18
+
let term = Mul (coeff, Pow (h, Const (float_of_int n))) in
19
+
build_terms (n + 1) (term :: acc)
20
+
in
21
+
let terms = List.rev (build_terms 0 []) in
22
+
let series = List.fold_left (fun acc t -> Add (acc, t)) (Const 0.0) terms in
23
+
simplify series
24
+
25
+
let maclaurin var expr order =
26
+
taylor var expr (Const 0.0) order
+220
lib/simplify.ml
+220
lib/simplify.ml
···
1
+
open Expr
2
+
3
+
let max_iterations = 100
4
+
5
+
let 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
+
31
+
and 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
+
46
+
and 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
+
53
+
and 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
+
80
+
and 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
+
90
+
and 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
+
103
+
and 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
+
110
+
and 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
+
117
+
and 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
+
124
+
and 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
+
131
+
and 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
+
137
+
and 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
+
143
+
and 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
+
149
+
and simplify_asin = function
150
+
| Const c -> Const (asin c)
151
+
| Sin e -> e
152
+
| e -> Asin e
153
+
154
+
and simplify_acos = function
155
+
| Const c -> Const (acos c)
156
+
| Cos e -> e
157
+
| e -> Acos e
158
+
159
+
and simplify_atan = function
160
+
| Const c -> Const (atan c)
161
+
| Tan e -> e
162
+
| e -> Atan e
163
+
164
+
and simplify_atan2 e1 e2 =
165
+
match (e1, e2) with
166
+
| Const a, Const b -> Const (atan2 a b)
167
+
| _ -> Atan2 (e1, e2)
168
+
169
+
and 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
+
177
+
and 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
+
187
+
and 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
+
196
+
and 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
+
204
+
and 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
+
211
+
let 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
+117
lib/special.ml
+117
lib/special.ml
···
1
+
open Expr
2
+
open Simplify
3
+
4
+
type special_function =
5
+
| Gamma of expr
6
+
| Beta of expr * expr
7
+
| Erf of expr
8
+
| Erfc of expr
9
+
| BesselJ of int * expr
10
+
| BesselY of int * expr
11
+
| LegendreP of int * expr
12
+
| HermiteH of int * expr
13
+
| LaguerreL of int * expr
14
+
| ChebyshevT of int * expr
15
+
| ChebyshevU of int * expr
16
+
17
+
let gamma x = Gamma x
18
+
let beta x y = Beta (x, y)
19
+
let erf x = Erf x
20
+
let erfc x = Erfc x
21
+
22
+
let bessel_j n x = BesselJ (n, x)
23
+
let bessel_y n x = BesselY (n, x)
24
+
25
+
let legendre_p n x =
26
+
let rec compute n =
27
+
if n = 0 then Const 1.0
28
+
else if n = 1 then x
29
+
else
30
+
let p_n_1 = compute (n - 1) in
31
+
let p_n_2 = compute (n - 2) in
32
+
let n_f = float_of_int n in
33
+
simplify (Div (
34
+
Sub (
35
+
Mul (Const (2.0 *. n_f -. 1.0), Mul (x, p_n_1)),
36
+
Mul (Const (n_f -. 1.0), p_n_2)
37
+
),
38
+
Const n_f
39
+
))
40
+
in
41
+
compute n
42
+
43
+
let hermite_h n x =
44
+
let rec compute n =
45
+
if n = 0 then Const 1.0
46
+
else if n = 1 then Mul (Const 2.0, x)
47
+
else
48
+
let h_n_1 = compute (n - 1) in
49
+
let h_n_2 = compute (n - 2) in
50
+
simplify (Sub (
51
+
Mul (Const 2.0, Mul (x, h_n_1)),
52
+
Mul (Const (2.0 *. float_of_int (n - 1)), h_n_2)
53
+
))
54
+
in
55
+
compute n
56
+
57
+
let laguerre_l n x =
58
+
let rec compute n =
59
+
if n = 0 then Const 1.0
60
+
else if n = 1 then Sub (Const 1.0, x)
61
+
else
62
+
let l_n_1 = compute (n - 1) in
63
+
let l_n_2 = compute (n - 2) in
64
+
let n_f = float_of_int n in
65
+
simplify (Div (
66
+
Sub (
67
+
Mul (Const (2.0 *. n_f -. 1.0 -. 1.0), l_n_1),
68
+
Sub (Mul (x, l_n_1), Mul (Const (n_f -. 1.0), l_n_2))
69
+
),
70
+
Const n_f
71
+
))
72
+
in
73
+
compute n
74
+
75
+
let chebyshev_t n x =
76
+
let rec compute n =
77
+
if n = 0 then Const 1.0
78
+
else if n = 1 then x
79
+
else
80
+
let t_n_1 = compute (n - 1) in
81
+
let t_n_2 = compute (n - 2) in
82
+
simplify (Sub (
83
+
Mul (Const 2.0, Mul (x, t_n_1)),
84
+
t_n_2
85
+
))
86
+
in
87
+
compute n
88
+
89
+
let chebyshev_u n x =
90
+
let rec compute n =
91
+
if n = 0 then Const 1.0
92
+
else if n = 1 then Mul (Const 2.0, x)
93
+
else
94
+
let u_n_1 = compute (n - 1) in
95
+
let u_n_2 = compute (n - 2) in
96
+
simplify (Sub (
97
+
Mul (Const 2.0, Mul (x, u_n_1)),
98
+
u_n_2
99
+
))
100
+
in
101
+
compute n
102
+
103
+
let factorial n =
104
+
let rec fact n acc =
105
+
if n <= 1 then acc
106
+
else fact (n - 1) (acc * n)
107
+
in
108
+
Const (float_of_int (fact n 1))
109
+
110
+
let binomial n k =
111
+
if k < 0 || k > n then Const 0.0
112
+
else
113
+
let rec binom n k =
114
+
if k = 0 || k = n then 1
115
+
else binom (n - 1) (k - 1) + binom (n - 1) k
116
+
in
117
+
Const (float_of_int (binom n k))
+194
lib/substitute.ml
+194
lib/substitute.ml
···
1
+
open Expr
2
+
3
+
let rec substitute var replacement = function
4
+
| Const _ as c -> c
5
+
| SymConst _ as s -> s
6
+
| Var v -> if v = var then replacement else Var v
7
+
| Add (e1, e2) -> Add (substitute var replacement e1, substitute var replacement e2)
8
+
| Sub (e1, e2) -> Sub (substitute var replacement e1, substitute var replacement e2)
9
+
| Mul (e1, e2) -> Mul (substitute var replacement e1, substitute var replacement e2)
10
+
| Div (e1, e2) -> Div (substitute var replacement e1, substitute var replacement e2)
11
+
| Pow (e1, e2) -> Pow (substitute var replacement e1, substitute var replacement e2)
12
+
| Neg e -> Neg (substitute var replacement e)
13
+
| Sin e -> Sin (substitute var replacement e)
14
+
| Cos e -> Cos (substitute var replacement e)
15
+
| Tan e -> Tan (substitute var replacement e)
16
+
| Sinh e -> Sinh (substitute var replacement e)
17
+
| Cosh e -> Cosh (substitute var replacement e)
18
+
| Tanh e -> Tanh (substitute var replacement e)
19
+
| Asin e -> Asin (substitute var replacement e)
20
+
| Acos e -> Acos (substitute var replacement e)
21
+
| Atan e -> Atan (substitute var replacement e)
22
+
| Atan2 (e1, e2) -> Atan2 (substitute var replacement e1, substitute var replacement e2)
23
+
| Exp e -> Exp (substitute var replacement e)
24
+
| Ln e -> Ln (substitute var replacement e)
25
+
| Log (e1, e2) -> Log (substitute var replacement e1, substitute var replacement e2)
26
+
| Sqrt e -> Sqrt (substitute var replacement e)
27
+
| Abs e -> Abs (substitute var replacement e)
28
+
29
+
let substitute_many subs expr =
30
+
List.fold_left (fun e (var, repl) -> substitute var repl e) expr subs
31
+
32
+
type pattern =
33
+
| PVar of string
34
+
| PWild
35
+
| PConst of float
36
+
| PSymConst of sym_const
37
+
| POp of string * pattern list
38
+
39
+
type bindings = (string * expr) list
40
+
41
+
let rec matches (pat : pattern) (e : expr) : bindings option =
42
+
match (pat, e) with
43
+
| PWild, _ -> Some []
44
+
| PVar v, e -> Some [(v, e)]
45
+
| PConst c1, Const c2 when c1 = c2 -> Some []
46
+
| PSymConst s1, SymConst s2 when s1 = s2 -> Some []
47
+
| POp ("Add", [p1; p2]), Add (e1, e2) ->
48
+
matches_binary p1 p2 e1 e2
49
+
| POp ("Sub", [p1; p2]), Sub (e1, e2) ->
50
+
matches_binary p1 p2 e1 e2
51
+
| POp ("Mul", [p1; p2]), Mul (e1, e2) ->
52
+
matches_binary p1 p2 e1 e2
53
+
| POp ("Div", [p1; p2]), Div (e1, e2) ->
54
+
matches_binary p1 p2 e1 e2
55
+
| POp ("Pow", [p1; p2]), Pow (e1, e2) ->
56
+
matches_binary p1 p2 e1 e2
57
+
| POp ("Neg", [p]), Neg e ->
58
+
matches p e
59
+
| POp ("Sin", [p]), Sin e ->
60
+
matches p e
61
+
| POp ("Cos", [p]), Cos e ->
62
+
matches p e
63
+
| POp ("Tan", [p]), Tan e ->
64
+
matches p e
65
+
| POp ("Sinh", [p]), Sinh e ->
66
+
matches p e
67
+
| POp ("Cosh", [p]), Cosh e ->
68
+
matches p e
69
+
| POp ("Tanh", [p]), Tanh e ->
70
+
matches p e
71
+
| POp ("Asin", [p]), Asin e ->
72
+
matches p e
73
+
| POp ("Acos", [p]), Acos e ->
74
+
matches p e
75
+
| POp ("Atan", [p]), Atan e ->
76
+
matches p e
77
+
| POp ("Atan2", [p1; p2]), Atan2 (e1, e2) ->
78
+
matches_binary p1 p2 e1 e2
79
+
| POp ("Exp", [p]), Exp e ->
80
+
matches p e
81
+
| POp ("Ln", [p]), Ln e ->
82
+
matches p e
83
+
| POp ("Log", [p1; p2]), Log (e1, e2) ->
84
+
matches_binary p1 p2 e1 e2
85
+
| POp ("Sqrt", [p]), Sqrt e ->
86
+
matches p e
87
+
| POp ("Abs", [p]), Abs e ->
88
+
matches p e
89
+
| _ -> None
90
+
91
+
and matches_binary p1 p2 e1 e2 =
92
+
match matches p1 e1 with
93
+
| None -> None
94
+
| Some b1 ->
95
+
match matches p2 e2 with
96
+
| None -> None
97
+
| Some b2 -> Some (b1 @ b2)
98
+
99
+
let rec instantiate (template : pattern) (bindings : bindings) : expr option =
100
+
match template with
101
+
| PWild -> None
102
+
| PVar v -> List.assoc_opt v bindings
103
+
| PConst c -> Some (Const c)
104
+
| PSymConst s -> Some (SymConst s)
105
+
| POp ("Add", [p1; p2]) ->
106
+
(match (instantiate p1 bindings, instantiate p2 bindings) with
107
+
| Some e1, Some e2 -> Some (Add (e1, e2))
108
+
| _ -> None)
109
+
| POp ("Sub", [p1; p2]) ->
110
+
(match (instantiate p1 bindings, instantiate p2 bindings) with
111
+
| Some e1, Some e2 -> Some (Sub (e1, e2))
112
+
| _ -> None)
113
+
| POp ("Mul", [p1; p2]) ->
114
+
(match (instantiate p1 bindings, instantiate p2 bindings) with
115
+
| Some e1, Some e2 -> Some (Mul (e1, e2))
116
+
| _ -> None)
117
+
| POp ("Div", [p1; p2]) ->
118
+
(match (instantiate p1 bindings, instantiate p2 bindings) with
119
+
| Some e1, Some e2 -> Some (Div (e1, e2))
120
+
| _ -> None)
121
+
| POp ("Pow", [p1; p2]) ->
122
+
(match (instantiate p1 bindings, instantiate p2 bindings) with
123
+
| Some e1, Some e2 -> Some (Pow (e1, e2))
124
+
| _ -> None)
125
+
| POp ("Neg", [p]) ->
126
+
(match instantiate p bindings with
127
+
| Some e -> Some (Neg e)
128
+
| None -> None)
129
+
| POp ("Sin", [p]) ->
130
+
(match instantiate p bindings with
131
+
| Some e -> Some (Sin e)
132
+
| None -> None)
133
+
| POp ("Cos", [p]) ->
134
+
(match instantiate p bindings with
135
+
| Some e -> Some (Cos e)
136
+
| None -> None)
137
+
| POp ("Tan", [p]) ->
138
+
(match instantiate p bindings with
139
+
| Some e -> Some (Tan e)
140
+
| None -> None)
141
+
| POp ("Sinh", [p]) ->
142
+
(match instantiate p bindings with
143
+
| Some e -> Some (Sinh e)
144
+
| None -> None)
145
+
| POp ("Cosh", [p]) ->
146
+
(match instantiate p bindings with
147
+
| Some e -> Some (Cosh e)
148
+
| None -> None)
149
+
| POp ("Tanh", [p]) ->
150
+
(match instantiate p bindings with
151
+
| Some e -> Some (Tanh e)
152
+
| None -> None)
153
+
| POp ("Asin", [p]) ->
154
+
(match instantiate p bindings with
155
+
| Some e -> Some (Asin e)
156
+
| None -> None)
157
+
| POp ("Acos", [p]) ->
158
+
(match instantiate p bindings with
159
+
| Some e -> Some (Acos e)
160
+
| None -> None)
161
+
| POp ("Atan", [p]) ->
162
+
(match instantiate p bindings with
163
+
| Some e -> Some (Atan e)
164
+
| None -> None)
165
+
| POp ("Atan2", [p1; p2]) ->
166
+
(match (instantiate p1 bindings, instantiate p2 bindings) with
167
+
| Some e1, Some e2 -> Some (Atan2 (e1, e2))
168
+
| _ -> None)
169
+
| POp ("Exp", [p]) ->
170
+
(match instantiate p bindings with
171
+
| Some e -> Some (Exp e)
172
+
| None -> None)
173
+
| POp ("Ln", [p]) ->
174
+
(match instantiate p bindings with
175
+
| Some e -> Some (Ln e)
176
+
| None -> None)
177
+
| POp ("Log", [p1; p2]) ->
178
+
(match (instantiate p1 bindings, instantiate p2 bindings) with
179
+
| Some e1, Some e2 -> Some (Log (e1, e2))
180
+
| _ -> None)
181
+
| POp ("Sqrt", [p]) ->
182
+
(match instantiate p bindings with
183
+
| Some e -> Some (Sqrt e)
184
+
| None -> None)
185
+
| POp ("Abs", [p]) ->
186
+
(match instantiate p bindings with
187
+
| Some e -> Some (Abs e)
188
+
| None -> None)
189
+
| _ -> None
190
+
191
+
let rewrite pattern template expr =
192
+
match matches pattern expr with
193
+
| None -> None
194
+
| Some bindings -> instantiate template bindings
+77
lib/transforms.ml
+77
lib/transforms.ml
···
1
+
open Expr
2
+
open Integrate
3
+
4
+
let fourier_transform expr var =
5
+
let omega = "omega" in
6
+
let integrand = Mul (expr, Exp (Mul (Mul (Const (-1.0), SymConst E), Mul (Var omega, Var var)))) in
7
+
match integrate var integrand with
8
+
| Some result -> Some result
9
+
| None ->
10
+
match expr with
11
+
| Const c -> Some (Mul (Const c, Const 0.0))
12
+
| Exp (Mul (Const a, Var v)) when v = var && a < 0.0 ->
13
+
Some (Div (Const 1.0, Add (Const a, Mul (SymConst E, Var omega))))
14
+
| Sin (Mul (Const a, Var v)) when v = var ->
15
+
Some (Div (Const a, Sub (Pow (Var omega, Const 2.0), Pow (Const a, Const 2.0))))
16
+
| Cos (Mul (Const a, Var v)) when v = var ->
17
+
Some (Div (Var omega, Sub (Pow (Var omega, Const 2.0), Pow (Const a, Const 2.0))))
18
+
| _ -> None
19
+
20
+
let inverse_fourier_transform _expr _omega =
21
+
None
22
+
23
+
let laplace_transform expr t =
24
+
let s = "s" in
25
+
match expr with
26
+
| Const c -> Some (Div (Const c, Var s))
27
+
| Var v when v = t -> Some (Div (Const 1.0, Pow (Var s, Const 2.0)))
28
+
| Pow (Var v, Const n) when v = t && Float.is_integer n && n >= 0.0 ->
29
+
let rec factorial n =
30
+
if n <= 1.0 then 1.0
31
+
else n *. factorial (n -. 1.0)
32
+
in
33
+
Some (Div (Const (factorial n), Pow (Var s, Const (n +. 1.0))))
34
+
| Exp (Mul (Const a, Var v)) when v = t ->
35
+
Some (Div (Const 1.0, Sub (Var s, Const a)))
36
+
| Sin (Mul (Const a, Var v)) when v = t ->
37
+
Some (Div (Const a, Add (Pow (Var s, Const 2.0), Pow (Const a, Const 2.0))))
38
+
| Cos (Mul (Const a, Var v)) when v = t ->
39
+
Some (Div (Var s, Add (Pow (Var s, Const 2.0), Pow (Const a, Const 2.0))))
40
+
| Sinh (Mul (Const a, Var v)) when v = t ->
41
+
Some (Div (Const a, Sub (Pow (Var s, Const 2.0), Pow (Const a, Const 2.0))))
42
+
| Cosh (Mul (Const a, Var v)) when v = t ->
43
+
Some (Div (Var s, Sub (Pow (Var s, Const 2.0), Pow (Const a, Const 2.0))))
44
+
| Mul (Exp (Mul (Const a, Var v1)), Sin (Mul (Const b, Var v2)))
45
+
when v1 = t && v2 = t ->
46
+
Some (Div (Const b,
47
+
Add (Pow (Sub (Var s, Const a), Const 2.0), Pow (Const b, Const 2.0))))
48
+
| Mul (Exp (Mul (Const a, Var v1)), Cos (Mul (Const b, Var v2)))
49
+
when v1 = t && v2 = t ->
50
+
Some (Div (Sub (Var s, Const a),
51
+
Add (Pow (Sub (Var s, Const a), Const 2.0), Pow (Const b, Const 2.0))))
52
+
| _ -> None
53
+
54
+
let inverse_laplace_transform expr s =
55
+
let t = "t" in
56
+
match expr with
57
+
| Div (Const c, Var v) when v = s -> Some (Const c)
58
+
| Div (Const 1.0, Pow (Var v, Const n)) when v = s && Float.is_integer n && n > 0.0 ->
59
+
let rec factorial n =
60
+
if n <= 1.0 then 1.0
61
+
else n *. factorial (n -. 1.0)
62
+
in
63
+
Some (Div (Pow (Var t, Const (n -. 1.0)), Const (factorial (n -. 1.0))))
64
+
| Div (Const 1.0, Sub (Var v, Const a)) when v = s ->
65
+
Some (Exp (Mul (Const a, Var t)))
66
+
| Div (Const a, Add (Pow (Var v, Const 2.0), Pow (Const b, Const 2.0))) when v = s ->
67
+
Some (Mul (Const (a /. b), Sin (Mul (Const b, Var t))))
68
+
| Div (Var v, Add (Pow (Var v2, Const 2.0), Pow (Const b, Const 2.0)))
69
+
when v = s && v2 = s ->
70
+
Some (Cos (Mul (Const b, Var t)))
71
+
| _ -> None
72
+
73
+
let z_transform _expr _n =
74
+
None
75
+
76
+
let inverse_z_transform _expr _z =
77
+
None
+135
test/test_suite.ml
+135
test/test_suite.ml
···
1
+
open OUnit2
2
+
open Leibniz
3
+
4
+
let test_lexer_basic _ =
5
+
let tokens = Lexer.tokenize "2 + 3" in
6
+
assert_equal 4 (List.length tokens)
7
+
8
+
let test_lexer_implicit_mult _ =
9
+
let tokens = Lexer.tokenize "2x" in
10
+
let token_list = List.map (fun t -> t.Lexer.token) tokens in
11
+
assert_equal 4 (List.length token_list);
12
+
match token_list with
13
+
| [Lexer.Num 2.0; Lexer.Star; Lexer.Var "x"; Lexer.EOF] -> ()
14
+
| _ -> assert_failure "implicit multiplication not working"
15
+
16
+
let test_parser_basic _ =
17
+
let expr = Parser.parse "x + 1" in
18
+
assert_equal (Expr.Add (Expr.Var "x", Expr.Const 1.0)) expr
19
+
20
+
let test_parser_implicit_mult _ =
21
+
let expr = Parser.parse "2x" in
22
+
assert_equal (Expr.Mul (Expr.Const 2.0, Expr.Var "x")) expr
23
+
24
+
let test_parser_symbolic_constants _ =
25
+
let expr = Parser.parse "pi" in
26
+
assert_equal (Expr.SymConst Expr.Pi) expr
27
+
28
+
let test_simplify_basic _ =
29
+
let expr = Expr.Add (Expr.Const 0.0, Expr.Var "x") in
30
+
let result = Simplify.simplify expr in
31
+
assert_equal (Expr.Var "x") result
32
+
33
+
let test_simplify_fixed_point _ =
34
+
let expr = Expr.Add (Expr.Mul (Expr.Const 0.0, Expr.Var "x"), Expr.Var "y") in
35
+
let result = Simplify.simplify expr in
36
+
assert_equal (Expr.Var "y") result
37
+
38
+
let test_simplify_collect_like_terms _ =
39
+
let expr = Expr.Add (Expr.Var "x", Expr.Var "x") in
40
+
let result = Simplify.simplify expr in
41
+
assert_equal (Expr.Mul (Expr.Const 2.0, Expr.Var "x")) result
42
+
43
+
let test_diff_basic _ =
44
+
let expr = Expr.Pow (Expr.Var "x", Expr.Const 2.0) in
45
+
let result = Diff.diff "x" expr in
46
+
let expected = Simplify.simplify (Expr.Mul (Expr.Const 2.0, Expr.Var "x")) in
47
+
assert_equal expected result
48
+
49
+
let test_diff_sin _ =
50
+
let expr = Expr.Sin (Expr.Var "x") in
51
+
let result = Diff.diff "x" expr in
52
+
assert_equal (Expr.Cos (Expr.Var "x")) result
53
+
54
+
let test_diff_product_rule _ =
55
+
let expr = Expr.Mul (Expr.Var "x", Expr.Sin (Expr.Var "x")) in
56
+
let _result = Diff.diff "x" expr in
57
+
()
58
+
59
+
let test_eval_basic _ =
60
+
let expr = Expr.Add (Expr.Var "x", Expr.Const 1.0) in
61
+
let result = Eval.eval [("x", 2.0)] expr in
62
+
assert_equal 3.0 result
63
+
64
+
let test_eval_symbolic_constants _ =
65
+
let expr = Expr.SymConst Expr.Pi in
66
+
let result = Eval.eval [] expr in
67
+
assert_bool "pi evaluation" (abs_float (result -. 3.14159265) < 0.0001)
68
+
69
+
let test_canonical_equality _ =
70
+
let e1 = Expr.Add (Expr.Var "x", Expr.Var "y") in
71
+
let e2 = Expr.Add (Expr.Var "y", Expr.Var "x") in
72
+
assert_bool "commutativity" (Canonical.equal e1 e2)
73
+
74
+
let test_substitute_basic _ =
75
+
let expr = Expr.Add (Expr.Var "x", Expr.Const 1.0) in
76
+
let result = Substitute.substitute "x" (Expr.Const 2.0) expr in
77
+
assert_equal (Expr.Add (Expr.Const 2.0, Expr.Const 1.0)) result
78
+
79
+
let test_integrate_basic _ =
80
+
let expr = Expr.Var "x" in
81
+
match Integrate.integrate "x" expr with
82
+
| Some result ->
83
+
let expected = Expr.Div (Expr.Pow (Expr.Var "x", Expr.Const 2.0), Expr.Const 2.0) in
84
+
assert_bool "integration result not equal" (Canonical.equal expected result)
85
+
| None -> assert_failure "integration failed"
86
+
87
+
let test_integrate_sin _ =
88
+
let expr = Expr.Sin (Expr.Var "x") in
89
+
match Integrate.integrate "x" expr with
90
+
| Some result ->
91
+
assert_equal (Expr.Neg (Expr.Cos (Expr.Var "x"))) result
92
+
| None -> assert_failure "integration of sin failed"
93
+
94
+
let test_taylor_sin _ =
95
+
let expr = Expr.Sin (Expr.Var "x") in
96
+
let _result = Series.maclaurin "x" expr 5 in
97
+
()
98
+
99
+
let test_gradient_basic _ =
100
+
let expr = Expr.Add (Expr.Pow (Expr.Var "x", Expr.Const 2.0),
101
+
Expr.Pow (Expr.Var "y", Expr.Const 2.0)) in
102
+
let grad = Multivariate.gradient ["x"; "y"] expr in
103
+
assert_equal 2 (List.length grad)
104
+
105
+
let test_numerical_bisection _ =
106
+
let expr = Expr.Sub (Expr.Pow (Expr.Var "x", Expr.Const 2.0), Expr.Const 4.0) in
107
+
match Numerical.bisection expr "x" 0.0 3.0 0.001 100 with
108
+
| Some root -> assert_bool "root near 2" (abs_float (root -. 2.0) < 0.01)
109
+
| None -> assert_failure "bisection failed"
110
+
111
+
let suite =
112
+
"leibniz tests" >::: [
113
+
"lexer basic" >:: test_lexer_basic;
114
+
"lexer implicit mult" >:: test_lexer_implicit_mult;
115
+
"parser basic" >:: test_parser_basic;
116
+
"parser implicit mult" >:: test_parser_implicit_mult;
117
+
"parser symbolic constants" >:: test_parser_symbolic_constants;
118
+
"simplify basic" >:: test_simplify_basic;
119
+
"simplify fixed point" >:: test_simplify_fixed_point;
120
+
"simplify collect like terms" >:: test_simplify_collect_like_terms;
121
+
"diff basic" >:: test_diff_basic;
122
+
"diff sin" >:: test_diff_sin;
123
+
"diff product rule" >:: test_diff_product_rule;
124
+
"eval basic" >:: test_eval_basic;
125
+
"eval symbolic constants" >:: test_eval_symbolic_constants;
126
+
"canonical equality" >:: test_canonical_equality;
127
+
"substitute basic" >:: test_substitute_basic;
128
+
"integrate basic" >:: test_integrate_basic;
129
+
"integrate sin" >:: test_integrate_sin;
130
+
"taylor sin" >:: test_taylor_sin;
131
+
"gradient basic" >:: test_gradient_basic;
132
+
"numerical bisection" >:: test_numerical_bisection;
133
+
]
134
+
135
+
let () = run_test_tt_main suite