symbolic mathematics engine in OCaml with differentiation, integration, simplification, and numerical methods
1open Expr
2open Simplify
3
4type 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
17let gamma x = Gamma x
18let beta x y = Beta (x, y)
19let erf x = Erf x
20let erfc x = Erfc x
21
22let bessel_j n x = BesselJ (n, x)
23let bessel_y n x = BesselY (n, x)
24
25let 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
43let 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
57let 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
75let 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
89let 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
103let 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
110let 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))