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