symbolic mathematics engine in OCaml with differentiation, integration, simplification, and numerical methods
17

Configure Feed

Select the types of activity you want to include in your feed.

leibniz / lib / groebner.ml
1.9 kB 64 lines
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 []