# type'a
cell
=
{
info
:
'a
;
mutable
prev
:
'a
dlist
;
mutable
next
:
'a
dlist
}
and
'a
dlist
=
Empty
|
Cell
of
'a
cell
;;
type 'a cell = { info: 'a; mutable prev: 'a dlist; mutable next: 'a dlist }
type 'a dlist = | Empty | Cell of 'a cell
# letadd
x
=
function
Empty
->
Cell
{
info
=
x;
prev
=
Empty;
next
=
Empty}
|
Cell
c
as
l
->
let
new_cell
=
{
info
=
x;
prev
=
c.
prev;
next
=
l}
in
let
new_dlist
=
Cell
new_cell
in
c
.
prev
<-
new_dlist
;
(
match
new_cell
.
prevwith
Empty
->
()
|
Cell
pl
->
pl
.
next
<-
new_dlist
)
;
new_dlist
;;
val add : 'a -> 'a dlist -> 'a dlist = <fun>
# letremove_cell
=
function
Empty
->
failwith
"Already empty"
|
Cell
c
->
match
(c
.
prev
,
c
.
next)with
Empty
,
Empty
->
Empty
|
Cell
c1
as
l
,
Empty
->
c1
.
next
<-
Empty
;
l
|
Empty
,
((Cell
c2)
as
l)
->
c2
.
prev
<-
Empty
;
l
|
Cell
c1
as
l1
,
(Cell
c2
as
l2)
->
c1
.
next
<-
l2;
c2
.
prev
<-
l1;
l1
;;
val remove_cell : 'a dlist -> 'a dlist = <fun>
# letrec
remove
x
l
=
let
rec
remove_left
=
function
Empty
->
()
|
Cell
c
as
l
->
let
pl
=
c
.
previn
if
c
.
info
=
x
then
ignore
(remove_cell
l)
;
remove_left
pl
and
remove_right
=
function
Empty
->
()
|
Cell
c
as
l
->
let
nl
=
c
.
nextin
if
c
.
info
=
x
then
ignore
(remove_cell
l)
;
remove_right
nl
in
match
l
with
Empty
->
Empty
|
Cell
c
as
l
->
if
c
.
info
=
x
then
remove
x
(remove_cell
l)
else
(remove_left
c
.
prev;
remove_right
c
.
next;
l)
;;
val remove : 'a -> 'a dlist -> 'a dlist = <fun>
# typevect
=
float
array
;;
type vect = float array
# typemat
=
vect
array
;;
type mat = vect array
# typesyst
=
{
m
:
mat;
v
:
vect}
;;
type syst = { m: mat; v: vect }
# letmy_print_float
s
=
let
x
=
string_of_float
s
in
let
y
=
match
String.length
x
with
5
->
x
|
n
when
n
<
5
->
(String.make
(
5
-
n)
' '
)
^
x
|
n
->
String.sub
x
0
5
in
print_string
y
;;
val my_print_float : float -> unit = <fun>
# letprint_syst
s
=
let
l
=
Array.length
s
.
m
in
for
i
=
0
to
l
-
1
do
print_string
" | "
;
for
j
=
0
to
l
-
1
do
my_print_float
s
.
m.
(i).
(j);
print_string
" "
done
;
if
i
=
l/
2
then
print_string
" | * | x"
else
print_string
" | | x"
;
print_int
(i
+
1
);
if
i
=
l/
2
then
print_string
" | = | "
else
print_string
" | | "
;
my_print_float
s
.
v.
(i);
print_string
" |"
;
print_newline
()
done
;;
val print_syst : syst -> unit = <fun>
# letadd_vect
v1
v2
=
let
l
=
Array.length
v1
in
let
res
=
Array.create
l
0
.
0
in
for
i
=
0
to
l
-
1
do
res
.
(i)
<-
v1
.
(i)
+.
v2
.
(i)done
;
res
;;
val add_vect : float array -> float array -> float array = <fun>
# letmult_scal_vect
x
v
=
let
l
=
Array.length
v
in
let
res
=
Array.create
l
0
.
0
in
for
i
=
0
to
l
-
1
do
res
.
(i)
<-
v
.
(i)
*.
x
done
;
res
;;
val mult_scal_vect : float -> float array -> float array = <fun>
# letmult_mat
m1
m2
=
let
l1
=
Array.length
m1
and
l2
=
Array.length
m2
.
(0
)
and
l3
=
Array.length
m2
in
let
res
=
Array.create_matrix
l1
l2
0
.
0
in
for
i
=
0
to
l1
-
1
do
for
j
=
0
to
l2
-
1
do
for
k
=
0
to
l3
-
1
do
res
.
(i).
(j)
<-
res
.
(i).
(j)
+.
m1
.
(i).
(k)
*.
m2
.
(k).
(j)
done
done
done
;
res
;;
val mult_mat : float array array -> float array array -> float array array =
<fun>
# letmult_mat_vect
m
v
=
let
l1
=
Array.length
m
and
l2
=
Array.length
v
in
let
res
=
Array.create
l1
0
.
0
in
for
i
=
0
to
l1
-
1
do
for
j
=
0
to
l2
-
1
do
res
.
(i)
<-
res
.
(i)
+.
m
.
(i).
(j)
*.
v
.
(j)
done
done
;
res
;;
val mult_mat_vect : float array array -> float array -> float array = <fun>
# letdiv_syst
s
i
=
let
p
=
s
.
m.
(i).
(i)
in
s
.
m.
(i).
(i)
<-
1
.
0
;
for
j
=
i+
1
to
(Array.length
s
.
m.
(0
))
-
1
do
s
.
m.
(i).
(j)
<-
s
.
m.
(i).
(j)
/.
p
done
;
s
.
v.
(i)
<-
s
.
v.
(i)
/.
p
;;
val div_syst : syst -> int -> unit = <fun>
# letpermut_syst
s
i
j
=
let
aux1
=
s
.
m.
(i)and
aux2
=
s
.
v.
(i)
in
s
.
m.
(i)
<-
s
.
m.
(j)
;
s
.
v.
(i)
<-
s
.
v.
(j);
s
.
m.
(j)
<-
aux1
;
s
.
v.
(j)
<-
aux2
;;
val permut_syst : syst -> int -> int -> unit = <fun>
# exceptionNot_linear
;;
exception Not_linear
# lettriangulize
s
=
if
s
.
m=[|
|]
||
s
.
v=[|
|]
then
raise
Not_linear
;
let
l
=
Array.length
s
.
m
in
if
l
<>
Array.lengths
.
m.
(0
)
||
l
<>
Array.lengths
.
vthen
raise
Not_linear
;
for
i
=
0
to
l
-
1
do
if
s
.
m.
(i).
(i)=
0
.
0
then
begin
let
j
=
ref
(i
+
1
)
in
while
!
j<
l
&&
s
.
m.
(!
j).
(i)=
0
.
0
do
incr
j
done
;
if
!
j=
lthen
raise
Not_linear
;
permut_syst
s
i
!
j
end
;
div_syst
s
i
;
for
j
=
i+
1
to
l
-
1
do
s
.
v.
(j)
<-
s
.
v.
(j)
-.
s
.
m.
(j).
(i)
*.
s
.
v.
(i);
s
.
m.
(j)
<-
add_vect
s
.
m.
(j)(mult_scal_vect
(
-.
s
.
m.
(j).
(i))s
.
m.
(i))
done
done
;;
val triangulize : syst -> unit = <fun>
# letsolve
s
=
triangulize
s
;
let
l
=
Array.length
s
.
v
in
let
res
=
Array.copy
s
.
v
in
for
i
=
l
-
1
downto
0
do
let
x
=
ref
res
.
(i)
in
for
j
=
i+
1
to
l
-
1
do
x
:=
!
x
-.
s
.
m.
(i).
(j)
*.
res
.
(j)done
;
res
.
(i)
<-
!
x
done
;
res
;;
val solve : syst -> float array = <fun>
AX = |
æ ç ç ç è |
|
ö ÷ ÷ ÷ ø |
* |
æ ç ç ç è |
|
ö ÷ ÷ ÷ ø |
= |
æ ç ç ç è |
|
ö ÷ ÷ ÷ ø |
= Y |
AX = |
æ ç ç ç è |
|
ö ÷ ÷ ÷ ø |
* |
æ ç ç ç è |
|
ö ÷ ÷ ÷ ø |
= |
æ ç ç ç è |
|
ö ÷ ÷ ÷ ø |
= Y |
AX = |
æ ç ç ç è |
|
ö ÷ ÷ ÷ ø |
* |
æ ç ç ç è |
|
ö ÷ ÷ ÷ ø |
= |
æ ç ç ç è |
|
ö ÷ ÷ ÷ ø |
= Y |
# letax1
=
{
m
=
[|
[|
1
0
.
0
;
7
.
0
;
8
.
0
;
7
.
0
|]
;
[|
7
.
0
;
5
.
0
;
6
.
0
;
5
.
0
|]
;
[|
8
.
0
;
6
.
0
;
1
0
.
0
;
9
.
0
|]
;
[|
7
.
0
;
5
.
0
;
9
.
0
;
1
0
.
0
|]
|]
;
v
=
[|
3
2
.
0
;
2
3
.
0
;
3
3
.
0
;
3
1
.
0
|]
}
;;
val ax1 : syst =
{m=[|[|10; 7; 8; 7|]; [|7; 5; 6; 5|]; [|8; 6; 10; ...|]; ...|]; v=...}
# letr1
=
solve
ax1
;;
val r1 : float array = [|1; 1; 1; 1|]
# letax2
=
{
m
=
[|
[|
1
0
.
0
;
7
.
0
;
8
.
0
;
7
.
0
|]
;
[|
7
.
0
;
5
.
0
;
6
.
0
;
5
.
0
|]
;
[|
8
.
0
;
6
.
0
;
1
0
.
0
;
9
.
0
|]
;
[|
7
.
0
;
5
.
0
;
9
.
0
;
1
0
.
0
|]
|]
;
v
=
[|
3
2
.
1
;
2
2
.
9
;
3
3
.
1
;
3
0
.
9
|]
}
;;
val ax2 : syst =
{m=[|[|10; 7; 8; 7|]; [|7; 5; 6; 5|]; [|8; 6; 10; ...|]; ...|]; v=...}
# letr2
=
solve
ax2
;;
val r2 : float array = [|9.2; -12.6; 4.5; -1.1|]
# letax3
=
{
m
=
[|
[|
1
0
.
0
;
7
.
0
;
8
.
1
;
7
.
2
|]
;
[|
7
.
0
8
;
5
.
0
4
;
6
.
0
;
5
.
0
|]
;
[|
8
.
0
;
5
.
9
8
;
9
.
8
9
;
9
.
0
|]
;
[|
6
.
9
9
;
4
.
9
9
;
9
.
0
;
9
.
9
8
|]
|]
;
v
=
[|
3
2
.
0
;
2
3
.
0
;
3
3
.
0
;
3
1
.
0
|]
}
;;
val ax3 : syst =
{m=
[|[|10; 7; 8.1; 7.2|]; [|7.08; 5.04; 6; 5|]; [|8; 5.98; 9.89; ...|];
...|];
v=...}
# letr3
=
solve
ax3
;;
val r3 : float array = [|-80.9999999999; 137; -33.9999999999; 22|]