Source file owl_maths_interpolate.ml
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
# 1 "src/base/maths/owl_maths_interpolate.ml"
(** Interpolation and Extrapolation *)
let polint xs ys x =
let n = Array.length xs in
let m = Array.length ys in
let error () =
let s =
Printf.sprintf
"polint requires that xs and ys have the same length, but xs length is %i \
whereas ys length is %i"
n
m
in
Owl_exception.INVALID_ARGUMENT s
in
Owl_exception.verify (m = n) error;
let c = Array.copy ys in
let d = Array.copy ys in
let ns = ref 0 in
let dy = ref 0. in
let dif = ref (abs_float (x -. xs.(0))) in
for i = 0 to n - 1 do
let dift = abs_float (x -. xs.(i)) in
if dift < !dif
then (
ns := i;
dif := dift)
done;
let y = ref ys.(!ns) in
ns := !ns - 1;
for m = 0 to n - 2 do
for i = 0 to n - m - 2 do
let ho = xs.(i) -. x in
let hp = xs.(i + m + 1) -. x in
let w = c.(i + 1) -. d.(i) in
let den = ho -. hp in
assert (den <> 0.);
let den = w /. den in
c.(i) <- ho *. den;
d.(i) <- hp *. den
done;
if (2 * !ns) + 2 < n - m - 1
then dy := c.(!ns + 1)
else (
dy := d.(!ns);
ns := !ns - 1);
y := !y +. !dy
done;
!y, !dy
let ratint xs ys x =
let n = Array.length xs in
let m = Array.length ys in
let error () =
let s =
Printf.sprintf
"ratint requires that xs and ys have the same length, but xs length is %i \
whereas ys length is %i"
n
m
in
Owl_exception.INVALID_ARGUMENT s
in
Owl_exception.verify (m = n) error;
let c = Array.copy ys in
let d = Array.copy ys in
let hh = ref (abs_float (x -. xs.(0))) in
let y = ref 0. in
let dy = ref 0. in
let ns = ref 0 in
let eps = 1e-25 in
try
for i = 0 to n do
let h = abs_float (x -. xs.(i)) in
if h = 0.
then (
y := ys.(i);
dy := 0.;
raise Owl_exception.FOUND)
else if h < !hh
then (
ns := i;
hh := h;
c.(i) <- ys.(i);
d.(i) <- ys.(i) +. eps)
done;
y := ys.(!ns);
ns := !ns - 1;
for m = 1 to n - 1 do
for i = 1 to n - m do
let w = c.(i) -. d.(i - 1) in
let h = xs.(i + m - 1) -. x in
let t = (xs.(i - 1) -. x) *. d.(i - 1) /. h in
let dd = t -. c.(i) in
if dd = 0. then failwith "Has a pole";
let dd = w /. dd in
d.(i - 1) <- c.(i) *. dd;
c.(i - 1) <- t *. dd
done
done;
!y, !dy
with
| Owl_exception.FOUND -> !y, !dy
| e -> raise e