symbolic mathematics engine in OCaml with differentiation, integration, simplification, and numerical methods
1open Polynomial
2
3type monomial_order = Lex | GrLex | GrevLex
4
5let 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
32let 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
37let s_polynomial _p1 _p2 _order =
38 []
39
40let reduce poly _basis _order =
41 poly
42
43let 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
58let 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
63let solve_polynomial_system _exprs _vars =
64 []