1 | % (c) 2020-2024 Lehrstuhl fuer Softwaretechnik und Programmiersprachen, | |
2 | % Heinrich Heine Universitaet Duesseldorf | |
3 | % This software is licenced under EPL 1.0 (http://www.eclipse.org/org/documents/epl-v10.html) | |
4 | ||
5 | :- module(kernel_reals, [construct_real/2, construct_negative_real/2, | |
6 | is_real/1, is_real/2, is_real_wf/2, | |
7 | is_float/1, is_float_wf/2, | |
8 | is_not_real/1, is_not_float/1, | |
9 | is_ground_real/1, | |
10 | is_largest_positive_float/1, is_smallest_positive_float/1, is_next_larger_float/2, | |
11 | convert_int_to_real/2, | |
12 | real_floor/2, real_ceiling/2, | |
13 | real_addition_wf/4, real_subtraction_wf/4, | |
14 | real_multiplication_wf/4, real_division_wf/5, | |
15 | real_unary_minus_wf/3, | |
16 | real_round_wf/3, | |
17 | real_unop_wf/4, real_unop_wf/5, | |
18 | real_binop_wf/6, | |
19 | real_power_of_wf/5, | |
20 | real_less_than_wf/3, real_less_than_equal_wf/3, | |
21 | real_comp_wf/5, | |
22 | real_maximum_of_set/4, real_minimum_of_set/4 | |
23 | ]). | |
24 | ||
25 | % Reals/Floats in ProB are represented by terms of the form term(floating(Number)) | |
26 | % in future a proper wrapper such as real(Number) might be created | |
27 | ||
28 | :- meta_predicate real_comp_wf(2,-,-,-,-). | |
29 | :- meta_predicate real_comp_wf_aux(2,-,-,-,-). | |
30 | ||
31 | :- use_module(module_information,[module_info/2]). | |
32 | :- module_info(group,kernel). | |
33 | :- module_info(description,'This module provides (external) functions to manipulate B reals and floats.'). | |
34 | ||
35 | :- use_module(error_manager). | |
36 | :- use_module(self_check). | |
37 | :- use_module(library(lists)). | |
38 | ||
39 | :- use_module(kernel_objects,[exhaustive_kernel_check_wf/2,exhaustive_kernel_check_wf/3]). | |
40 | :- use_module(extension('counter/counter'), | |
41 | [next_smaller_abs_float/2, next_smaller_float/2,next_larger_float/2, next_float_twoards/3, | |
42 | smallest_float/1, largest_float/1, smallest_abs_float/1]). | |
43 | ||
44 | % used to construct a Value from real(Atom) AST node | |
45 | construct_real(Atom,term(floating(Float))) :- | |
46 | atom_codes(Atom,C), | |
47 | number_codes(Nr,C), | |
48 | Float is float(Nr). % make sure we store a float; not required for AST literals | |
49 | ||
50 | construct_negative_real(Atom,term(floating(Float))) :- | |
51 | atom_codes(Atom,C), | |
52 | number_codes(Nr,C), | |
53 | Float is -float(Nr). | |
54 | ||
55 | is_real_wf(X,_WF) :- is_real(X,_). | |
56 | is_real(X) :- is_real(X,_). | |
57 | ||
58 | is_real(term(floating(Nr)),Nr). | |
59 | % TO DO: other formats | |
60 | ||
61 | is_float_wf(X,_WF) :- is_float(X). | |
62 | ||
63 | is_float(term(floating(_))). | |
64 | ||
65 | is_not_real(_) :- fail. | |
66 | is_not_float(_) :- fail. | |
67 | ||
68 | is_ground_real(term(floating(Nr))) :- number(Nr). | |
69 | ||
70 | is_largest_positive_float(term(floating(Nr))) :- largest_float(Nr). | |
71 | is_smallest_positive_float(term(floating(Nr))) :- smallest_abs_float(Nr). | |
72 | ||
73 | is_next_larger_float(term(floating(Nr)),term(floating(Next))) :- | |
74 | block_next_larger_float(Nr,Next). | |
75 | ||
76 | :- block block_next_larger_float(-,-). | |
77 | block_next_larger_float(Nr,Next) :- nonvar(Nr),!,next_larger_float(Nr,Next). | |
78 | block_next_larger_float(Nr,Next) :- next_smaller_float(Next,Nr). | |
79 | ||
80 | % convert an integer to a real, corresponds to the real(.) function in Atelier-B; convert_real as AST | |
81 | convert_int_to_real(int(X),term(floating(R))) :- | |
82 | convert_int_to_real_aux(X,R). | |
83 | ||
84 | :- block convert_int_to_real_aux(-,-). | |
85 | convert_int_to_real_aux(Int,Real) :- number(Int),!, | |
86 | Real is float(Int). | |
87 | convert_int_to_real_aux(Int,Real) :- Int1 is floor(Real), | |
88 | Real is float(Int1), | |
89 | Int1=Int. | |
90 | ||
91 | ||
92 | % convert_int_floor as AST, Atelier-B floor(.) | |
93 | real_floor(term(floating(R)),int(X)) :- real_floor_aux(R,X). | |
94 | :- block real_floor_aux(-,?). | |
95 | real_floor_aux(Real,Int) :- Int is floor(Real). | |
96 | ||
97 | % convert_int_ceiling as AST, Atelier-B ceiling(.) | |
98 | real_ceiling(term(floating(R)),int(X)) :- real_ceiling_aux(R,X). | |
99 | :- block real_ceiling_aux(-,?). | |
100 | real_ceiling_aux(Real,Int) :- Int is ceiling(Real). | |
101 | ||
102 | :- assert_must_succeed(exhaustive_kernel_check_wf([commutative],kernel_reals:real_addition_wf(term(floating(1.0)),term(floating(2.0)),term(floating(3.0)),WF),WF)). | |
103 | :- assert_must_succeed(exhaustive_kernel_check_wf([commutative],kernel_reals:real_addition_wf(term(floating(0.0)),term(floating(2.0)),term(floating(2.0)),WF),WF)). | |
104 | :- assert_must_succeed(exhaustive_kernel_check_wf([commutative],kernel_reals:real_addition_wf(term(floating(-1.0)),term(floating(1.0)),term(floating(0.0)),WF),WF)). | |
105 | ||
106 | real_addition_wf(term(floating(X)),term(floating(Y)),term(floating(R)),WF) :- | |
107 | real_add_wf_aux(X,Y,R,WF). | |
108 | :- block real_add_wf_aux(-,-,?,?), real_add_wf_aux(?,-,-,?), real_add_wf_aux(-,?,-,?). | |
109 | real_add_wf_aux(X,Y,R,WF) :- | |
110 | (nonvar(X) | |
111 | -> (nonvar(Y) | |
112 | -> safe_is(R,X+Y,WF) | |
113 | ; allow_constraint_propagation(Kind), | |
114 | safe_is(YY,R-X,WF), | |
115 | check_solution_for_add(YY,X,R,Kind,WF) | |
116 | -> Y=YY % deterministic constraint propagation found solution | |
117 | ; real_add_wf_aux2(X,Y,R,WF) % delay until Y is known | |
118 | ) | |
119 | ; allow_constraint_propagation(Kind), | |
120 | safe_is(XX,R-Y,WF), | |
121 | check_solution_for_add(XX,Y,R,Kind,WF) -> X=XX % deterministic constraint propagation found solution | |
122 | ; real_add_wf_aux2(X,Y,R,WF) % delay until X is known | |
123 | ). | |
124 | ||
125 | :- block real_add_wf_aux2(-,?,?,?), real_add_wf_aux2(?,-,?,?). | |
126 | real_add_wf_aux2(X,Y,R,WF) :- safe_is(R,X+Y,WF). | |
127 | ||
128 | check_solution_for_add(_,_,_,no_checking,_) :- !. | |
129 | check_solution_for_add(XX,Y,R,simple_checking,WF) :- !, | |
130 | % (x + y = z <=> x = z - y) does not hold for all floats, e.g., when x very small and y very large (x+y)-y = 0 | |
131 | (safe_is(R,XX+Y,WF) -> true ; format('Reject non-solution ~w + ~w = ~w~n',[XX,Y,R]),fail). | |
132 | check_solution_for_add(XX,Y,R,_,WF) :- | |
133 | (safe_is(R,XX+Y,WF) -> true ; format('Reject non-solution ~w + ~w = ~w~n',[XX,Y,R]),fail), | |
134 | next_larger_float(XX,XN), | |
135 | (safe_is(R,XN+Y,WF) -> format('Reject ambiguous larger solution {~w,~w} + ~w = ~w~n',[XX,XN,Y,R]),fail ; true), | |
136 | next_smaller_float(XX,XP), | |
137 | (safe_is(R,XP+Y,WF) -> format('Reject ambiguous smaller solution {~w,~w} + ~w = ~w~n',[XX,XP,Y,R]),fail ; true). | |
138 | % TODO: create choice point if only few (two solutions) | |
139 | % e.g. for x+1.0 = 3.0 we have 2.0 as solution and the previous float 1.9999999999999998, nothing else | |
140 | ||
141 | ||
142 | :- use_module(probsrc(preferences), [get_preference/2]). | |
143 | ||
144 | allow_constraint_propagation(Kind) :- | |
145 | get_preference(solver_for_reals,Solver), | |
146 | get_solver_kind(Solver,Kind). | |
147 | ||
148 | get_solver_kind(aggressive_float_solver,no_checking). % like CLP(R): do not check propagated solution | |
149 | get_solver_kind(float_solver,simple_checking). % check propagated solution, but do not check uniqueness | |
150 | get_solver_kind(precise_float_solver,check_unique_solution). | |
151 | ||
152 | % we could have various other options: | |
153 | % aggressive: do it like CLP(R), which can generate non-solutions | |
154 | % | ?- {X+10.0e10 = 1.0e-9}, write(sol(X)),nl, {X+10.0e10 = 1.0e-9}. | |
155 | % prints sol(-1.0E+11) but then fails | |
156 | % conservative: do the precision check and also check that there are no other solutions (not done yet) | |
157 | % e.g., X+10000000000.0 = 10000000000.0 & (X=0.000000000001 or X=2.0) is FALSE with the float_solver, but TRUE with none | |
158 | ||
159 | ||
160 | :- assert_must_succeed(exhaustive_kernel_check_wf(kernel_reals:real_subtraction_wf(term(floating(3.0)),term(floating(2.0)),term(floating(1.0)),WF),WF)). | |
161 | real_subtraction_wf(term(floating(X)),term(floating(Y)),term(floating(R)),WF) :- | |
162 | real_sub_wf_aux(X,Y,R,WF). | |
163 | :- block real_sub_wf_aux(-,-,?,?), real_sub_wf_aux(?,-,-,?), real_sub_wf_aux(-,?,-,?). | |
164 | real_sub_wf_aux(X,Y,R,WF) :- | |
165 | (nonvar(X) | |
166 | -> (nonvar(Y) | |
167 | -> safe_is(R,X-Y,WF) | |
168 | ; allow_constraint_propagation(_), | |
169 | safe_is(YY,X-R,WF), | |
170 | safe_is(R,X-YY,WF) -> Y=YY % deterministic constraint propagation found precise solution | |
171 | ; real_sub_wf_aux2(X,Y,R,WF) % delay until Y is known | |
172 | ) | |
173 | ; allow_constraint_propagation(_), | |
174 | safe_is(XX,R+Y,WF), | |
175 | safe_is(R,XX-Y,WF) -> X=XX % deterministic constraint propagation found precise solution | |
176 | ; real_sub_wf_aux2(X,Y,R,WF) % delay until X is known | |
177 | ). | |
178 | ||
179 | :- block real_sub_wf_aux2(-,?,?,?), real_sub_wf_aux2(?,-,?,?). | |
180 | real_sub_wf_aux2(X,Y,R,WF) :- safe_is(R,X-Y,WF). | |
181 | ||
182 | % TODO: check solutions | |
183 | ||
184 | :- assert_must_succeed(exhaustive_kernel_check_wf([commutative],kernel_reals:real_multiplication_wf(term(floating(3.0)),term(floating(2.0)),term(floating(6.0)),WF),WF)). | |
185 | real_multiplication_wf(term(floating(X)),term(floating(Y)),term(floating(R)),WF) :- | |
186 | real_mul_wf_aux(X,Y,R,WF). | |
187 | :- block real_mul_wf_aux(-,?,?,?), real_mul_wf_aux(?,-,?,?). | |
188 | real_mul_wf_aux(X,Y,R,WF) :- safe_is(R,X*Y,WF). | |
189 | ||
190 | :- assert_must_succeed(exhaustive_kernel_check_wf(kernel_reals:real_division_wf(term(floating(6.0)),term(floating(2.0)),term(floating(3.0)),unknown,WF),WF)). | |
191 | real_division_wf(term(floating(X)),term(floating(Y)),term(floating(R)),Span,WF) :- | |
192 | real_div_wf_aux(X,Y,R,Span,WF). | |
193 | :- block real_div_wf_aux(-,?,?,?,?), real_div_wf_aux(?,-,?,?,?). | |
194 | real_div_wf_aux(X,Y,R,Span,WF) :- safe_is(R,X/Y,Span,WF). | |
195 | ||
196 | ||
197 | real_unary_minus_wf(term(floating(X)),term(floating(R)),WF) :- | |
198 | real_um_wf_aux(X,R,WF). | |
199 | :- block real_um_wf_aux(-,?,?). | |
200 | real_um_wf_aux(X,R,_) :- R is -X. | |
201 | ||
202 | real_round_wf(term(floating(X)),int(R),_WF) :- | |
203 | when(nonvar(X), R is round(X)). | |
204 | ||
205 | :- use_module(kernel_waitflags,[add_wd_error/3, add_wd_error_span/4]). | |
206 | % a version of is/2 which catches overflows | |
207 | safe_is(R,Expr,WF) :- | |
208 | safe_is(R,Expr,unknown,WF). | |
209 | safe_is(R,Expr,Span,WF) :- | |
210 | catch(R is Expr, | |
211 | error(evaluation_error(ERR),_), | |
212 | process_evaluation_error(ERR,Expr,Span,WF)). | |
213 | ||
214 | process_evaluation_error(float_overflow,Expr,Span,WF) :- !, | |
215 | add_wd_error_span('Float Overflow while computing:',Expr,Span,WF). | |
216 | process_evaluation_error(zero_divisor,Expr,Span,WF) :- !, | |
217 | add_wd_error_span('Division by zero while computing:',Expr,Span,WF). | |
218 | process_evaluation_error(undefined,Expr,Span,WF) :- !, | |
219 | add_wd_error_span('Arithmetic operator undefined while computing:',Expr,Span,WF). | |
220 | process_evaluation_error(_,Expr,Span,WF) :- | |
221 | add_wd_error_span('Unknown evaluation error while computing:',Expr,Span,WF). | |
222 | ||
223 | ||
224 | % call a Prolog unary artihmetic operator | |
225 | real_unop_wf(OP,X,R,WF) :- | |
226 | real_unop_wf(OP,X,R,unknown,WF). | |
227 | ||
228 | real_unop_wf(OP,RX,RR,Span,WF) :- | |
229 | get_real(RX,X,OP,Span,WF), get_real(RR,R,OP,Span,WF), | |
230 | real_unop_wf_aux(OP,X,R,Span,WF). | |
231 | :- block real_unop_wf_aux(-,?,?,?,?), real_unop_wf_aux(?,-,?,?,?). | |
232 | real_unop_wf_aux(OP,X,R,Span,WF) :- | |
233 | Expr =.. [OP,X], | |
234 | safe_is(R,Expr,Span,WF). | |
235 | ||
236 | :- use_module(probsrc(tools_strings),[ajoin/2]). | |
237 | :- use_module(kernel_waitflags,[add_error_wf/5]). | |
238 | get_real(Var,Real,_,_,_) :- var(Var),!, Var=term(floating(Real)). | |
239 | get_real(term(F),Real,_,_,_) :- !, F=floating(Real). | |
240 | get_real(Other,_,OP,Span,WF) :- | |
241 | ajoin(['Argument for ',OP,' is not a real number:'],Msg), | |
242 | add_error_wf(kernel_reals,Msg,Other,Span,WF). | |
243 | ||
244 | ||
245 | real_power_of_wf(X,Y,Res,Span,WF) :- | |
246 | real_binop_wf('exp',X,Y,Res,Span,WF). | |
247 | ||
248 | % call a Prolog binary artihmetic operator | |
249 | real_binop_wf(OP,RX,RY,RR,Span,WF) :- | |
250 | get_real(RX,X,OP,Span,WF), get_real(RY,Y,OP,Span,WF), get_real(RR,R,OP,Span,WF), | |
251 | real_binnop_wf_aux(OP,X,Y,R,Span,WF). | |
252 | :- block real_binnop_wf_aux(-,?,?,?,?,?), real_binnop_wf_aux(?,-,?,?,?,?), real_binnop_wf_aux(?,?,-,?,?,?). | |
253 | real_binnop_wf_aux(OP,X,Y,R,Span,WF) :- | |
254 | Expr =.. [OP,X,Y], | |
255 | safe_is(R,Expr,Span,WF). | |
256 | ||
257 | real_less_than_wf(X,Y,WF) :- real_comp_wf('<',X,Y,pred_true,WF). | |
258 | real_less_than_equal_wf(X,Y,WF) :- real_comp_wf('=<',X,Y,pred_true,WF). | |
259 | ||
260 | ||
261 | % call a Prolog bianry artihmetic operator | |
262 | real_comp_wf(OP,term(floating(X)),term(floating(Y)),R,WF) :- | |
263 | real_comp_wf_aux(OP,X,Y,R,WF). | |
264 | :- block real_comp_wf_aux(-,?,?,?,?), real_comp_wf_aux(?,-,?,?,?), real_comp_wf_aux(?,?,-,?,?). | |
265 | real_comp_wf_aux(OP,X,Y,R,_) :- | |
266 | (call(OP,X,Y) -> R=pred_true ; R=pred_false). | |
267 | ||
268 | % ----------------- | |
269 | ||
270 | :- use_module(probsrc(kernel_tools),[ground_value_check/2]). | |
271 | real_maximum_of_set(Set,Res,Span,WF) :- | |
272 | ground_value_check(Set,Gr), | |
273 | rmax_set(Gr,Set,Res,Span,WF). | |
274 | ||
275 | :- block rmax_set(-,?,?,?,?). | |
276 | rmax_set(_,Set,Res,Span,WF) :- | |
277 | custom_explicit_sets:expand_and_convert_to_avl_set(Set,ESet,'RMAXIMUM','RMAXIMUM'), | |
278 | (ESet=empty | |
279 | -> add_wd_error_span('max applied to empty set of reals:',Set,Span,WF) | |
280 | ; custom_explicit_sets:max_of_explicit_set_wf(avl_set(ESet),Res,WF) | |
281 | ). | |
282 | ||
283 | real_minimum_of_set(Set,Res,Span,WF) :- | |
284 | ground_value_check(Set,Gr), | |
285 | rmin_set(Gr,Set,Res,Span,WF). | |
286 | ||
287 | :- block rmin_set(-,?,?,?,?). | |
288 | rmin_set(_,Set,Res,Span,WF) :- | |
289 | custom_explicit_sets:expand_and_convert_to_avl_set(Set,ESet,'RMINIMUM','RMINIMUM'), | |
290 | (ESet=empty | |
291 | -> add_wd_error_span('min applied to empty set of reals:',Set,Span,WF) | |
292 | ; custom_explicit_sets:min_of_explicit_set_wf(avl_set(ESet),Res,WF) | |
293 | ). | |
294 | ||
295 | ||
296 | % not used yet: useful for kernel_strings ? | |
297 | %:- block real_to_string(-,?). | |
298 | %real_to_string(term(floating(I)),S) :- real_to_string2(I,S). | |
299 | % | |
300 | %:- block real_to_string2(-,?). | |
301 | %real_to_string2(Num,Res) :- | |
302 | % number_codes(Num,C), | |
303 | % atom_codes(S,C), Res=string(S). | |
304 | ||
305 |