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