1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
|
(* ocamlgsl - OCaml interface to GSL *)
(* Copyright () 2002 - Olivier Andrieu *)
(* distributed under the terms of the GPL version 2 *)
open Gsl_complex
type poly = float array
external eval : poly -> float -> float
= "ml_gsl_poly_eval"
type quad_sol =
| Quad_0
| Quad_2 of float * float
external solve_quadratic : a:float -> b:float -> c:float -> quad_sol
= "ml_gsl_poly_solve_quadratic"
external complex_solve_quadratic : a:float -> b:float -> c:float -> complex * complex
= "ml_gsl_poly_complex_solve_quadratic"
type cubic_sol =
| Cubic_0
| Cubic_1 of float
| Cubic_3 of float * float * float
external solve_cubic : a:float -> b:float -> c:float -> cubic_sol
= "ml_gsl_poly_solve_cubic"
external complex_solve_cubic : a:float -> b:float -> c:float -> complex * complex * complex
= "ml_gsl_poly_complex_solve_cubic"
type ws
external _alloc_ws : int -> ws = "ml_gsl_poly_complex_workspace_alloc"
external _free_ws : ws -> unit= "ml_gsl_poly_complex_workspace_free"
external _solve : poly -> ws -> complex_array -> unit
= "ml_gsl_poly_complex_solve"
let solve poly =
let n = Array.length poly in
let ws = _alloc_ws n in
let sol = Array.make (2*(n-1)) 0. in
_solve poly ws sol ;
_free_ws ws ;
sol
|
