64 lines
1.8 KiB
OCaml
64 lines
1.8 KiB
OCaml
open Polynomial
|
|
|
|
type monomial_order = Lex | GrLex | GrevLex
|
|
|
|
let rec compare_monomials order m1 m2 =
|
|
let total_degree powers = List.fold_left (fun acc (_, n) -> acc + n) 0 powers in
|
|
match order with
|
|
| Lex ->
|
|
let rec lex_compare p1 p2 =
|
|
match (p1, p2) with
|
|
| [], [] -> 0
|
|
| [], _ -> -1
|
|
| _, [] -> 1
|
|
| (v1, n1) :: rest1, (v2, n2) :: rest2 ->
|
|
let c = String.compare v1 v2 in
|
|
if c <> 0 then c
|
|
else if n1 <> n2 then compare n2 n1
|
|
else lex_compare rest1 rest2
|
|
in
|
|
lex_compare m1 m2
|
|
| GrLex ->
|
|
let d1 = total_degree m1 in
|
|
let d2 = total_degree m2 in
|
|
if d1 <> d2 then compare d2 d1
|
|
else compare_monomials Lex m1 m2
|
|
| GrevLex ->
|
|
let d1 = total_degree m1 in
|
|
let d2 = total_degree m2 in
|
|
if d1 <> d2 then compare d2 d1
|
|
else compare_monomials Lex m2 m1
|
|
|
|
let leading_term poly order =
|
|
match List.sort (fun t1 t2 -> compare_monomials order t1.powers t2.powers) poly with
|
|
| [] -> {coeff = 0.0; powers = []}
|
|
| t :: _ -> t
|
|
|
|
let s_polynomial _p1 _p2 _order =
|
|
[]
|
|
|
|
let reduce poly _basis _order =
|
|
poly
|
|
|
|
let buchberger polys order =
|
|
let rec loop basis =
|
|
let pairs = List.concat_map (fun p1 ->
|
|
List.filter_map (fun p2 ->
|
|
if p1 != p2 then Some (p1, p2) else None
|
|
) basis
|
|
) basis in
|
|
let s_polys = List.map (fun (p1, p2) -> s_polynomial p1 p2 order) pairs in
|
|
let reduced = List.map (fun sp -> reduce sp basis order) s_polys in
|
|
let non_zero = List.filter (fun p -> List.length p > 0) reduced in
|
|
if List.length non_zero = 0 then basis
|
|
else loop (basis @ non_zero)
|
|
in
|
|
loop polys
|
|
|
|
let groebner_basis exprs vars =
|
|
let polys = List.map (fun e -> expr_to_poly e vars) exprs in
|
|
let basis = buchberger polys GrLex in
|
|
List.map poly_to_expr basis
|
|
|
|
let solve_polynomial_system _exprs _vars =
|
|
[]
|