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
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
(** This module will contain code for 2nd Degree equation resolution *)
open Ast
(** [convert ast] takes the polynome represented as a recursive type and
converts it to a type
[(coefficient:float, (variable:string * power:int) option) list] *)
let convert : polynome -> monome list =
let rec aux sub = function
| Add (a, b) -> aux sub a @ aux false b
| Sub (a, b) -> aux sub a @ aux true b
| Mon (x, v) ->
if sub then
[ (~-.x, v) ]
else
[ (x, v) ]
in
aux false
(** [ast_to_lists (p1, p2)] takes a type [(polynome * polynome)] and calls
[convert ast] on each polynome to return a [(monome list * monome list)] to
facilitate manipulation *)
let ast_to_lists ((p1, p2) : polynome * polynome) = (convert p1, convert p2)
(** [sqrt x] is a reimplementation of classic sqrt function, it's presence is
justified by the specifications of the ComputorV1 42 subject *)
let sqrt (x : float) =
let rec sqrt_aux guess = function
| 0. -> guess
| iter ->
let guess = guess -. (((guess *. guess) -. x) /. (2. *. guess)) in
sqrt_aux guess (iter -. 1.)
in
sqrt_aux 1. 10.
(*
Alternative, iterative version of the custom SQRT function re-implementation
let sqrt x =
let z = ref 1. in
for i = 10 downto 0 do
begin
z := !z -. (!z *. !z -. x) /. (2. *. !z)
end
done;
!z
*)
(** [solve e] takes both monome lists and does all necessary operations that
lead to solution. Invert all monomes of the righthand side polynome and
concatenate to lefthand side one. Check for more than one variable, exit if
so. Create a Hashtable of [(key:power, value:coefficient)]. Check for a
degree > 2, exit if so. Simplify list by combining same power monomes.
Proceed to calculate delta and pretty-print the reduced form along with
polynomial degree and solution(s). *)
let solve fmt (e : monome list * monome list) =
(* Separate both polynomes of the equation left and right from '=' sign *)
let left, right = e in
(* Move right polynome to the left of the '=' sign, invert signs *)
let right' = List.map (fun (coeff, var) -> (~-.coeff, var)) right in
(* Concatenate *)
let p = left @ right' in
let seen = ref None in
(* Scan for presence of more than one variable, exit if there is *)
List.iter
(fun (_coeff, var) ->
match var with
| None -> ()
| Some (var, _n) -> (
match !seen with
| None -> seen := Some var
| Some seen ->
if seen <> var then begin
Format.fprintf fmt
"There are two different variables `%s` and `%s`!@." var seen;
Format.fprintf fmt "Won't compute.@.";
raise Too_many_variables
end ) )
p;
(*
* Produce new list with monomes matching
* (coeff:float * (var:string, 0) option)
* replaced with (1, None)
*)
let poly =
List.map
(fun (coeff, var) ->
match (coeff, var) with
| coeff, None -> (coeff, var)
| coeff, Some (_x, 0) -> (coeff, None)
| coeff, Some (_x, _n) -> (coeff, var) )
p
in
(*
* Create hashtable
* Keys are Power to which a given coefficient(varible option) is raised
* Values are the said coefficients
*
* Add coefficients of same power together to simplify polynome
*)
let tbl = Hashtbl.create 512 in
List.iter
(fun (coeff, var) ->
match var with
| None -> (
let x = Hashtbl.find_opt tbl 0 in
match x with
| None -> Hashtbl.add tbl 0 coeff
| Some x -> Hashtbl.replace tbl 0 (coeff +. x) )
| Some (_var, n) -> (
let x = Hashtbl.find_opt tbl n in
match x with
| None -> Hashtbl.replace tbl n coeff
| Some x -> Hashtbl.replace tbl n (coeff +. x) ) )
poly;
let var = Option.value !seen ~default:"" in
(* Convert all Hashtbl of monomes (key:power, value:coef) to a list of (power, coef)*)
let allterms = List.of_seq (Hashtbl.to_seq tbl) in
(* Order in growing order of (power, _coef) *)
let sorted_terms =
List.sort (fun (pow1, _) (pow2, _) -> Int.compare pow1 pow2) allterms
in
(* Filter out the monomes that have a coef of 0 *)
let filtered_terms = List.filter (fun x -> snd x <> 0.) sorted_terms in
(* Find highest degree with non null coefficient *)
let max_degree =
try fst (List.hd (List.rev filtered_terms)) with
| Failure _ -> 0
in
(* Raise an exception in case of degree > 2 *)
if max_degree > 2 then begin
Format.fprintf fmt "Reduced form: %a@." Pp.equ (filtered_terms, [], var);
Format.fprintf fmt "Polynomial degree: %d@." max_degree;
Format.fprintf fmt
"The polynomial degree is strictly greater than 2, I can't solve.@.";
raise Big_degree
end;
let a = Option.value (Hashtbl.find_opt tbl 2) ~default:0. in
let b = Option.value (Hashtbl.find_opt tbl 1) ~default:0. in
let c = Option.value (Hashtbl.find_opt tbl 0) ~default:0. in
let normalize fmt root =
Format.fprintf fmt "%g"
( if Float.compare (-0.) root = 0 then
Float.abs root
else
root )
in
(* If a is not 0, then it's not a 2nd degree equation *)
if a <> 0. then (
(* Calculating delta for resolution of
* 2nd Degree equation
*)
let delta = (b *. b) -. (4. *. a *. c) in
let normalize fmt root =
Format.fprintf fmt "%g"
( if Float.compare (-0.) root = 0 then
Float.abs root
else
root )
in
match delta with
| 0. ->
(* Delta is null, there is only one possible solution *)
let solution = ~-.b /. (2. *. a) in
Format.fprintf fmt "Reduced form: %a@." Pp.equ (filtered_terms, [], var);
Format.fprintf fmt "Polynomial degree: %d@." max_degree;
Pp.deltap fmt delta;
Format.fprintf fmt "%a@." normalize solution
(* Delta is positive, 2 solutions possible *)
| delta when Float.compare 0. delta < 0 ->
let solution1 = (~-.b -. sqrt delta) /. (2. *. a) in
let solution2 = (~-.b +. sqrt delta) /. (2. *. a) in
Format.fprintf fmt "Reduced form: %a@." Pp.equ (filtered_terms, [], var);
Format.fprintf fmt "Polynomial degree: %d@." max_degree;
Pp.deltap fmt delta;
Format.fprintf fmt "%a@.%a@." normalize solution1 normalize solution2
(* Solutions are found in complex numbers *)
| delta ->
let reduce_zi fmt root =
if Float.compare 1. root = 0 then
Format.fprintf fmt ""
else
Format.fprintf fmt "%g" root
in
Format.fprintf fmt "Reduced form: %a@." Pp.equ (filtered_terms, [], var);
Format.fprintf fmt "Polynomial degree: %d@." max_degree;
Pp.deltap fmt delta;
let root1 = ~-.b /. (2. *. a) in
let root2 = sqrt ~-.delta /. (2. *. a) in
Format.fprintf fmt "%a - %ai@." normalize root1 reduce_zi root2;
Format.fprintf fmt "%a + %ai@." normalize root1 reduce_zi root2
) else
(* Calculating x for resolution of
* 1st or 0th Degree equation
*)
match b with
| 0. -> (
match c with
(* Case 0 = 0
*)
| 0. ->
Format.fprintf fmt "Reduced form: 0 = 0@.";
Format.fprintf fmt "Polynomial degree: %d@." max_degree;
Pp.deltap fmt 0.;
Format.fprintf fmt "All real numbers are solutions.\n"
(* Case 3 = 0
*)
| _c ->
Format.fprintf fmt "Reduced form: %a@." Pp.equ (filtered_terms, [], var);
Format.fprintf fmt "Impossible.@." )
(* Case x = -c/b
*)
| b ->
let solution = ~-.c /. b in
Format.fprintf fmt "Reduced form: %a@." Pp.equ (filtered_terms, [], var);
Format.fprintf fmt "Polynomial degree: %d@." max_degree;
Pp.deltap fmt 0.;
Format.fprintf fmt "%a@." normalize solution