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