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
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
|
/* ocamlgsl - OCaml interface to GSL */
/* Copyright () 2002 - Olivier Andrieu */
/* distributed under the terms of the GPL version 2 */
#include <caml/mlvalues.h>
#include <caml/alloc.h>
#include <caml/memory.h>
#include <caml/fail.h>
#include <gsl/gsl_poly.h>
#include "wrappers.h"
value ml_gsl_poly_eval(value c, value x)
{
int len = Double_array_length(c);
return copy_double(gsl_poly_eval(Double_array_val(c), len, Double_val(x)));
}
value ml_gsl_poly_solve_quadratic(value a, value b, value c)
{
double x0, x1;
int n ;
n = gsl_poly_solve_quadratic(Double_val(a), Double_val(b),
Double_val(c), &x0, &x1);
{
CAMLparam0();
CAMLlocal1(r);
if(n == 0)
r = Val_int(0);
else{
r = alloc_small(2, 0);
Store_field(r, 0, copy_double(x0));
Store_field(r, 1, copy_double(x1));
} ;
CAMLreturn(r);
}
}
value ml_gsl_poly_complex_solve_quadratic(value a, value b, value c)
{
gsl_complex z0, z1;
gsl_poly_complex_solve_quadratic(Double_val(a), Double_val(b),
Double_val(c), &z0, &z1);
{
CAMLparam0();
CAMLlocal3(r,rz0,rz1);
rz0 = alloc_small(2 * Double_wosize, Double_array_tag);
Store_double_field(rz0, 0, GSL_REAL(z0));
Store_double_field(rz0, 1, GSL_IMAG(z0));
rz1 = alloc_small(2 * Double_wosize, Double_array_tag);
Store_double_field(rz1, 0, GSL_REAL(z1));
Store_double_field(rz1, 1, GSL_IMAG(z1));
r = alloc_small(2, 0);
Field(r,0) = rz0 ;
Field(r,1) = rz1 ;
CAMLreturn(r);
}
}
value ml_gsl_poly_solve_cubic(value a, value b, value c)
{
double x0, x1, x2;
int n ;
n = gsl_poly_solve_cubic(Double_val(a), Double_val(b),
Double_val(c), &x0, &x1, &x2);
{
CAMLparam0();
CAMLlocal1(r);
r = Val_int(0); /* to silence compiler warnings */
switch(n){
case 0:
break;
case 1:
r = alloc_small(1, 0);
Store_field(r, 0, copy_double(x0));
break;
case 3:
r = alloc_small(3, 1);
Store_field(r, 0, copy_double(x0));
Store_field(r, 1, copy_double(x1));
Store_field(r, 2, copy_double(x2));
} ;
CAMLreturn(r);
};
}
value ml_gsl_poly_complex_solve_cubic(value a, value b, value c)
{
gsl_complex z0, z1, z2;
gsl_poly_complex_solve_cubic(Double_val(a), Double_val(b),
Double_val(c), &z0, &z1, &z2);
{
CAMLparam0();
CAMLlocal4(r,rz0, rz1, rz2);
rz0 = alloc_small(2 * Double_wosize, Double_array_tag);
Store_double_field(rz0, 0, GSL_REAL(z0));
Store_double_field(rz0, 1, GSL_IMAG(z0));
rz1 = alloc_small(2 * Double_wosize, Double_array_tag);
Store_double_field(rz1, 0, GSL_REAL(z1));
Store_double_field(rz1, 1, GSL_IMAG(z1));
rz2 = alloc_small(2 * Double_wosize, Double_array_tag);
Store_double_field(rz2, 0, GSL_REAL(z2));
Store_double_field(rz2, 1, GSL_IMAG(z2));
r = alloc_small(3, 0);
Field(r,0) = rz0 ;
Field(r,1) = rz1 ;
Field(r,2) = rz2 ;
CAMLreturn(r);
}
}
#define POLY_WS(v) (gsl_poly_complex_workspace *)Field((v), 0)
ML1_alloc(gsl_poly_complex_workspace_alloc, Int_val, Abstract_ptr)
ML1(gsl_poly_complex_workspace_free, POLY_WS, Unit)
value ml_gsl_poly_complex_solve(value a, value ws, value r)
{
gsl_poly_complex_solve(Double_array_val(a), Double_array_length(a),
POLY_WS(ws), (gsl_complex_packed_ptr) r);
return Val_unit;
}
|
