1 % (c) 2020-2026 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 construct_real_number/2,
7 is_real/1, is_real/2, is_real_wf/2,
8 is_float/1, is_float_wf/2,
9 is_not_real/1, is_not_float/1,
10 is_ground_real/1,
11 is_largest_positive_float/1, is_smallest_positive_float/1,
12 is_next_larger_float/2, is_next_smaller_float/2,
13 convert_int_to_real/2,
14 real_floor/2, real_ceiling/2, real_truncate/2,
15 real_addition_wf/4, real_subtraction_wf/4,
16 real_multiplication_wf/4, real_division_wf/5,
17 real_unary_minus_wf/3,
18 real_absolute_value_wf/3,
19 real_sign_wf/3,
20 real_square_root_wf/4,
21 real_round_wf/3,
22 real_unop_wf/4, real_unop_wf/5,
23 real_binop_wf/6,
24 real_power_of_wf/5,
25 real_less_than_wf/3, real_less_than_equal_wf/3,
26 real_comp_wf/5,
27 real_maximum_of_set/4, real_minimum_of_set/4,
28 enumerate_real_wf/3,
29 use_clpfd_real_solver/0, do_not_double_check_solution/0,
30 post_real_equal_expr_wf/4,
31 post_real_neq/2,
32 post_real_equality/3, % reified equality
33 is_interval_closure_or_realset/3,
34 is_real_interval_closure_with_known_card/4
35 ]).
36
37 % Reals/Floats in ProB are represented by terms of the form term(floating(Number))
38 % in future a proper wrapper such as real(Number) might be created
39
40 :- meta_predicate real_comp_wf(2,-,-,-,-).
41 :- meta_predicate real_comp_float_blocking(2,-,-,-).
42
43 :- use_module(module_information,[module_info/2]).
44 :- module_info(group,kernel).
45 :- module_info(description,'This module provides (external) functions to manipulate B reals and floats.').
46
47 :- use_module(error_manager).
48 :- use_module(self_check).
49 :- use_module(library(lists)).
50 :- use_module(debug).
51
52 :- use_module(kernel_objects,[exhaustive_kernel_check_wf/2,exhaustive_kernel_check_wf/3]).
53 :- use_module(extension('counter/counter'),
54 [next_smaller_abs_float/2, next_smaller_float/2,next_larger_float/2, next_float_twoards/3,
55 smallest_float/1, largest_float/1, smallest_abs_float/1]).
56
57
58 :- if((current_prolog_flag(dialect, sicstus),
59 current_prolog_flag(version_data, sicstus(4,VN,_,_,_)), VN>=10)). % 4.10 or higher
60 can_use_clpfd_real_solver.
61 use_clpfd_real_solver :- get_preference(solver_for_reals,S), (S=precise_float_solver ; S=real_solver).
62 use_clpfd_real_solver(S) :- get_preference(solver_for_reals,S), (S=precise_float_solver ; S=real_solver).
63 :- use_module(library(clpfd),[fd_min/2, fd_max/2, fd_batch/1, maximum/2, minimum/2,
64 '#<=>'/2, '$>='/2, '$=<'/2, '$='/2, '#='/2, labeling/2]).
65
66 post_real_equal_expr_wf(X,Y,Span,WF) :-
67 catch('$='(X,Y),
68 error(evaluation_error(ERR),_),
69 process_evaluation_error(ERR,'='(X,Y),Span,WF)).
70 :- else.
71 can_use_clpfd_real_solver :- fail.
72 use_clpfd_real_solver :- fail.
73 use_clpfd_real_solver(_) :- fail.
74 post_real_equal_expr_wf(X,Y,Span,WF) :- add_internal_error('Not available:',post_real_equal_expr_wf(X,Y,Span,WF)),fail.
75 :- endif.
76
77 % there is currently no $\= in CLP(FD) for reals; use regular dif instead
78 % this call only works for variables or numbers as arguments, no expressions
79 post_real_neq(E1,E2) :- dif(E1,E2).
80
81
82 :- use_module(kernel_dif,[frozen_dif/2]).
83 % reified version of equality between X and Y
84 post_real_equality(X,Y,Res) :- nonvar(X),nonvar(Y),!,
85 (X=Y -> Res=pred_true ; Res=pred_false).
86 post_real_equality(X,Y,Res) :- % check for attached dif co-routines; later we should probably use CHR
87 frozen_dif(X,Y),!, Res=pred_false.
88 post_real_equality(X,Y,R) :- use_clpfd_real_solver(S),!,
89 post_real_solver_equality(S,X,Y,R).
90 post_real_equality(X,Y,Res) :- when((?=(X,Y);nonvar(Res)),eq_real2(X,Y,Res)).
91
92 post_real_solver_equality(Solver,X,Y,R) :- R==pred_true,!, post_equal_real(Solver,X,Y).
93 post_real_solver_equality(Solver,X,Y,R) :-
94 '#<=>'('$='(X,Y),R01),
95 % Note: '#<=>'( '$='(X,1.0), R01), X = 1.0. does not instantiate R01, and allows R01=0 (see SPRM-21561)
96 % Hence, if X=Y we force R01 to be 1:
97 when(?=(X,Y),(X=Y -> R01=1; R01=0)),
98 prop_pred_01(Solver,R,R01,X,Y). % link pred_false/pred_true to 0/1
99
100 post_equal_real(precise_float_solver,X,Y) :- !, X=Y. % force full equality
101 post_equal_real(_,X,Y) :- '$='(X,Y),
102 (frozen_dif(X,Y) -> fail ; true).
103
104 :- block prop_pred_01(?,-,-,?,?).
105 prop_pred_01(_S,A,B,_X,_Y) :- B==1,!,A=pred_true. % cut ok: either pred_true or 1 set
106 prop_pred_01(Solver,pred_true,1,X,Y) :- post_equal_real(Solver,X,Y).
107 prop_pred_01(_,pred_false,0,X,Y) :- dif(X,Y).
108
109
110 eq_real2(X,Y,Res) :- nonvar(Res),!,
111 (Res=pred_true -> X=Y ; post_real_neq(X,Y)).
112 eq_real2(X,Y,Res) :-
113 (X=Y -> Res=pred_true ; Res=pred_false).
114
115 % used to construct a Value from real(Atom) AST node
116 construct_real(Atom,term(floating(Float))) :-
117 atom_codes(Atom,C),
118 number_codes(Nr,C),
119 Float is float(Nr). % make sure we store a float; not required for AST literals
120
121 construct_real_number(Atom,Nr) :-
122 atom_codes(Atom,C),
123 number_codes(Nr,C).
124
125 construct_negative_real(Atom,term(floating(Float))) :-
126 atom_codes(Atom,C),
127 number_codes(Nr,C),
128 Float is -float(Nr).
129
130 is_real_wf(X,_WF) :- is_real(X,_).
131 is_real(X) :- is_real(X,_).
132
133 is_real(term(floating(Nr)),Nr).
134 % TO DO: other formats
135
136 is_float_wf(X,_WF) :- is_float(X).
137
138 is_float(term(floating(_))).
139
140 is_not_real(_) :- fail.
141 is_not_float(_) :- fail.
142
143 is_ground_real(term(floating(Nr))) :- number(Nr).
144
145 is_largest_positive_float(term(floating(Nr))) :- largest_float(Nr).
146 is_smallest_positive_float(term(floating(Nr))) :- smallest_abs_float(Nr).
147
148 % corresponds to RNEXT external function
149 is_next_larger_float(term(floating(Nr)),term(floating(Next))) :-
150 block_next_larger_float(Nr,Next).
151 :- block block_next_larger_float(-,-).
152 block_next_larger_float(Nr,Next) :- nonvar(Nr),!,next_larger_float(Nr,Next).
153 block_next_larger_float(Nr,Next) :- next_smaller_float(Next,Nr).
154
155 % corresponds to RPREV external function
156 is_next_smaller_float(term(floating(Nr)),term(floating(Next))) :-
157 block_next_larger_float(Next,Nr).
158
159
160 % convert an integer to a real, corresponds to the real(.) function in Atelier-B; convert_real as AST
161 convert_int_to_real(int(X),term(floating(R))) :-
162 convert_int_to_real_aux(X,R).
163
164 convert_int_to_real_aux(Int,Real) :- number(Int),!,
165 Real is float(Int).
166 convert_int_to_real_aux(Int,Real) :-
167 use_clpfd_real_solver,!,
168 '$='(Real,float(Int)).
169 convert_int_to_real_aux(Int,Real) :- block_convert_int_to_real_aux(Int,Real).
170
171 :- block block_convert_int_to_real_aux(-,-).
172 block_convert_int_to_real_aux(Int,Real) :- number(Int),!,
173 Real is float(Int).
174 block_convert_int_to_real_aux(Int,Real) :- Int1 is floor(Real),
175 Real is float(Int1),
176 Int1=Int.
177
178
179 % convert_int_floor as AST, Atelier-B floor(.)
180 % The value is the greatest integer less or equal to X
181 real_floor(term(floating(R)),int(X)) :- real_floor_aux(R,X).
182
183 real_floor_aux(Real,Int) :-
184 use_clpfd_real_solver,!, % + check if we use CLP(FD) for integers?
185 '#='(Int,floor(Real)).
186 real_floor_aux(Real,Int) :- block_real_floor_aux(Real,Int).
187
188 :- block block_real_floor_aux(-,?).
189 block_real_floor_aux(Real,Int) :- Int is floor(Real).
190
191
192 % convert_int_ceiling as AST, Atelier-B ceiling(.)
193 % The value is the least integer greater or equal to X.
194 real_ceiling(term(floating(R)),int(X)) :-
195 real_ceiling_aux(R,X).
196 real_ceiling_aux(Real,Int) :-
197 use_clpfd_real_solver,!, % + check if we use CLP(FD) for integers?
198 '#='(Int,ceiling(Real)).
199 real_ceiling_aux(Real,Int) :- block_real_ceiling_aux(Real,Int).
200
201 :- block block_real_ceiling_aux(-,?).
202 block_real_ceiling_aux(Real,Int) :- Int is ceiling(Real).
203
204
205 % The value is the closest integer between X and 0
206 real_truncate(term(floating(R)),int(X)) :-
207 real_truncate_aux(R,X).
208 real_truncate_aux(Real,Int) :-
209 use_clpfd_real_solver,!, % + check if we use CLP(FD) for integers?
210 '#='(Int,truncate(Real)).
211 real_truncate_aux(Real,Int) :- block_real_truncate_aux(Real,Int).
212
213 :- block block_real_truncate_aux(-,?).
214 block_real_truncate_aux(Real,Int) :- Int is truncate(Real).
215
216
217
218 :- 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)).
219 :- 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)).
220 :- 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)).
221
222 real_addition_wf(term(floating(X)),term(floating(Y)),term(floating(R)),WF) :-
223 real_add_wf_aux(X,Y,R,WF).
224
225 %real_add_wf_aux(X,Y,R,WF) :- X==Y,!, '$='(R,2.0*X), % TODO: check if this works reliably, probably not
226 % (do_not_double_check_solution -> true ; real_add_wf_aux2(X,Y,R,WF)).
227 real_add_wf_aux(X,Y,R,WF) :-
228 use_clpfd_real_solver,!,
229 '$='(R,X+Y),
230 (do_not_double_check_solution -> true ; real_add_wf_aux2(X,Y,R,WF)).
231 real_add_wf_aux(X,Y,R,WF) :- block_real_add_wf_aux(X,Y,R,WF).
232
233 :- block block_real_add_wf_aux(-,-,?,?), block_real_add_wf_aux(?,-,-,?), block_real_add_wf_aux(-,?,-,?).
234 block_real_add_wf_aux(X,Y,R,WF) :-
235 (nonvar(X)
236 -> (nonvar(Y)
237 -> safe_is(R,X+Y,WF)
238 ; allow_constraint_propagation(Kind),
239 safe_is(YY,R-X,WF),
240 confirm_solution_for_expr(YY,X+SOL,SOL,R,Kind,Solutions,WF)
241 -> set_solution_for_expr(Y,Solutions,WF)
242 ; real_add_wf_aux2(X,Y,R,WF) % delay until Y is known
243 )
244 ; allow_constraint_propagation(Kind),
245 safe_is(XX,R-Y,WF),
246 confirm_solution_for_expr(XX,SOL+Y,SOL,R,Kind,Solutions,WF)
247 -> set_solution_for_expr(X,Solutions,WF)
248 ; real_add_wf_aux2(X,Y,R,WF) % delay until X is known
249 ).
250 :- block real_add_wf_aux2(-,?,?,?), real_add_wf_aux2(?,-,?,?).
251 real_add_wf_aux2(X,Y,R,WF) :- safe_is(R,X+Y,WF).
252
253
254 % confirm_solution_for_expr(XX,ExprX,X,R : check if XX=X is a solution for R = ExprX, where ExprX contains variable X
255 % depending on solver setting : compute a list of candidate solutions around XX
256 confirm_solution_for_expr(XX,_,_,_,no_checking,Solutions,_) :- !, Solutions=XX.
257 confirm_solution_for_expr(XX,ExprX,X,R,simple_checking,Solutions,WF) :- !,
258 check_sol_direct(XX,ExprX,X,R,WF), Solutions=XX. % single solution
259 confirm_solution_for_expr(XX,ExprX,X,R,check_unique_solution,Solutions,WF) :- !,
260 check_sol_direct(XX,ExprX,X,R,WF), Solutions=XX, % single solution
261 % Now check if there are other solutions:
262 next_larger_float(XX,XN),
263 (X=XN,safe_is(R,ExprX,WF)
264 -> debug_format(19,'Reject ambiguous larger solution {~w,~w} for ~w = ~w~n',[XX,XN,ExprX,R]),fail
265 ; true
266 ),
267 next_smaller_float(XX,XP),
268 (X=XP,safe_is(R,ExprX,WF)
269 -> debug_format(19,'Reject ambiguous smaller solution {~w,~w} for ~w = ~w~n',[XX,XP,ExprX,R]),fail
270 ; true).
271 confirm_solution_for_expr(XX,ExprX,X,R,check_small_interval_solution,Interval,WF) :-
272 get_solutions_around(XX,ExprX,X,R,Interval,WF).
273
274 :- use_module(kernel_waitflags,[get_wait_flag/4]).
275 % set the solution based on result of confirm_solution_for_expr
276 set_solution_for_expr(X,XX,_WF) :- number(XX),!,X=XX. % deterministic constraint propagation found precise solution
277 set_solution_for_expr(X,[XX],_WF) :- !,X=XX. % ditto
278 set_solution_for_expr(_,[],_WF) :- !, fail. % no solutions
279 set_solution_for_expr(X,Interval,WF) :-
280 length(Interval,Len), % we found an interval list of possible solutions
281 get_wait_flag(Len,kernel_reals,WF,LWF),
282 set_sol_aux(X,Interval,LWF).
283 :- block set_sol_aux(-,?,-).
284 set_sol_aux(X,Interval,_) :-
285 ? member(X,Interval).
286
287
288 % get a sorted list of solutions around the candidate XX
289 get_solutions_around(XX,ExprX,X,R,Interval,WF) :-
290 (check_sol_direct(XX,ExprX,X,R,WF) % XX is a solution
291 -> get_preference(solver_strength,SS),
292 Limit is 10+SS, % number of solutions we are willing to find
293 next_larger_float(XX,XN),
294 find_larger_bounds(XN,Limit,ExprX,X,R,WF,LargerSols,[]),
295 next_smaller_float(XX,XP),
296 find_lower_bounds(XP,Limit,ExprX,X,R,WF,Interval,[XX|LargerSols])
297 ; (X='X',debug_format(19,'Reject non-solution ~w for ~w = ~w~n',[XX,ExprX,R]),fail ; true),
298 %next_smaller_float(XX,XP), (check_sol_direct(XP,ExprX,X,R,WF) -> write(xp(XP)),nl ; write(nxp(XP)),nl),
299 %next_larger_float(XX,XN), (check_sol_direct(XN,ExprX,X,R,WF) -> write(xn(XN)),nl ; write(nxp(XN)),nl),
300 Interval = [] % no floating point solutions; TODO: do we need to check smaller/larger values??
301 ).
302
303 % e.g. for x+1.0 = 3.0 we have 2.0 as solution and the previous float 1.9999999999999998, nothing else
304 % e.b. for x+2.0 = 3.0 we have four solutions {0.9999999999999998,0.9999999999999999,1.0,1.0000000000000002}
305 % card({x|x + 2.999 = 3.0}) = 2049
306 % card({x|x + 2.9999 = 3.0}) = 32769
307
308 % increment floats until we hit the first non-solution or we reach the limit
309 % TODO: try and compute this bound, rather then single-stepping
310 % TODO: enumerate the first solution and then generate an enum warning if we reach limit
311 find_larger_bounds(XP,Lim,ExprX,X,R,WF,[XP|Sols1],Sols2) :-
312 check_sol_direct(XP,ExprX,X,R,WF),!, Lim>0, L1 is Lim-1,
313 next_larger_float(XP,XP2),
314 find_larger_bounds(XP2,L1,ExprX,X,R,WF,Sols1,Sols2).
315 find_larger_bounds(_Bound,_Lim,_ExprX,_X,_R,_WF,Sols,Sols).
316
317 % decrement floats until we hit the first non-solution or we reach the limit
318 find_lower_bounds(XP,Lim,ExprX,X,R,WF,[XP|Sols1],Sols2) :-
319 check_sol_direct(XP,ExprX,X,R,WF),!, Lim>0, L1 is Lim-1,
320 next_smaller_float(XP,XP2),
321 find_lower_bounds(XP2,L1,ExprX,X,R,WF,Sols1,Sols2).
322 find_lower_bounds(_Bound,_Lim,_ExprX,_X,_R,_WF,Sols,Sols).
323
324 check_sol_direct(XX,ExprX,X,R,WF) :- !,
325 % (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
326 (\+((X=XX,safe_is(R,ExprX,WF))) -> fail ; true).
327
328
329 :- use_module(probsrc(preferences), [get_preference/2]).
330
331 allow_constraint_propagation(Kind) :-
332 get_preference(solver_for_reals,Solver),
333 get_solver_kind(Solver,Kind).
334
335 get_solver_kind(real_solver,no_checking). % like CLP(R): do not check propagated solution
336 get_solver_kind(float_solver,simple_checking). % check propagated solution, but do not check uniqueness
337 %get_solver_kind(precise_float_solver,check_unique_solution).
338 get_solver_kind(precise_float_solver,check_small_interval_solution). % perform limited constraint propagation, fail directly if no solution exists
339
340 do_not_double_check_solution :- get_preference(solver_for_reals,real_solver).
341
342 precise_float_solver_active :- get_preference(solver_for_reals,precise_float_solver).
343
344 % we could have various other options:
345 % aggressive: do it like CLP(R), which can generate non-solutions
346 % | ?- {X+10.0e10 = 1.0e-9}, write(sol(X)),nl, {X+10.0e10 = 1.0e-9}.
347 % prints sol(-1.0E+11) but then fails
348 % conservative: do the precision check and also check that there are no other solutions (not done yet)
349 % e.g., X+10000000000.0 = 10000000000.0 & (X=0.000000000001 or X=2.0) is FALSE with the float_solver, but TRUE with none
350
351
352 :- 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)).
353 real_subtraction_wf(term(floating(X)),term(floating(Y)),term(floating(R)),WF) :-
354 real_sub_wf_aux(X,Y,R,WF).
355
356 real_sub_wf_aux(X,Y,R,WF) :-
357 use_clpfd_real_solver,!,
358 '$='(R,X-Y),
359 (do_not_double_check_solution -> true ; real_sub_wf_aux2(X,Y,R,WF)).
360 real_sub_wf_aux(X,Y,R,WF) :- block_real_sub_wf_aux(X,Y,R,WF).
361
362 :- block block_real_sub_wf_aux(-,-,?,?), block_real_sub_wf_aux(?,-,-,?), block_real_sub_wf_aux(-,?,-,?).
363 block_real_sub_wf_aux(X,Y,R,WF) :-
364 (nonvar(X)
365 -> (nonvar(Y)
366 -> safe_is(R,X-Y,WF)
367 ; allow_constraint_propagation(Kind),
368 safe_is(YY,X-R,WF),
369 confirm_solution_for_expr(YY,X-SOL,SOL,R,Kind,Solutions,WF)
370 -> set_solution_for_expr(Y,Solutions,WF) % deterministic constraint propagation found precise solution
371 ; real_sub_wf_aux2(X,Y,R,WF) % delay until Y is known
372 )
373 ; allow_constraint_propagation(Kind),
374 safe_is(XX,R+Y,WF),
375 confirm_solution_for_expr(XX,SOL-Y,SOL,R,Kind,Solutions,WF)
376 -> set_solution_for_expr(X,Solutions,WF) % deterministic constraint propagation found precise solution
377 ; real_sub_wf_aux2(X,Y,R,WF) % delay until X is known
378 ).
379
380 :- block real_sub_wf_aux2(-,?,?,?), real_sub_wf_aux2(?,-,?,?).
381 real_sub_wf_aux2(X,Y,R,WF) :- safe_is(R,X-Y,WF).
382
383
384 :- 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)).
385 real_multiplication_wf(term(floating(X)),term(floating(Y)),term(floating(R)),WF) :-
386 real_mul_wf_aux0(X,Y,R,WF).
387
388 real_mul_wf_aux0(X,Y,R,WF) :-
389 use_clpfd_real_solver,!,
390 (X==Y -> '$='(R,X*X) ; % causes segmentation faults up to SICS 4.10.2-beta2
391 '$='(R,X*Y),propagate_square_inverse(X,Y,R,WF)), %tools_printing:print_term_summary(mul(R,X,Y)),nl,
392 (do_not_double_check_solution -> true ; real_mul_wf_aux2(X,Y,R,WF)).
393 real_mul_wf_aux0(X,Y,R,WF) :- block_real_mul_wf_aux0(X,Y,R,WF).
394
395 % | ?- 2.0 $= X*0.0. -> clpfd:(X in fsup..fsup) ? i.e., CLPFD does not detect multiplication by 0.0
396 % CLPFD also does not take square root
397 :- block propagate_square_inverse(-,-,-,?).
398 propagate_square_inverse(X,Y,R,_WF) :- (X==0.0 ; Y==0.0), !, R=0.0.
399 propagate_square_inverse(X,Y,R,WF) :- var(X),X==Y, !, % X*X = R -> compute sqrt
400 %'$='(R,X^2.0). % makes test 2430 fail for c=card({y|y*y=x}) & x = 2.0E-309 & c=14 -> c=1
401 %'$='(R,X*X). % x * x = 1.0 -> segmentation fault
402 kernel_square_inverse(X,R,WF).
403 propagate_square_inverse(X,Y,R,_WF) :- nonvar(R), R \= 0.0, !,
404 dif(X,0.0), dif(Y,0.0).
405 propagate_square_inverse(_,_,_,_).
406
407 :- block block_real_mul_wf_aux0(-,-,-,?).
408 block_real_mul_wf_aux0(X,Y,R,_WF) :- (X==0.0 ; Y==0.0), !, R=0.0.
409 block_real_mul_wf_aux0(X,Y,R,WF) :- var(X),X==Y, !, % X*X = R -> compute sqrt
410 kernel_square_inverse(X,R,WF).
411 block_real_mul_wf_aux0(X,Y,R,WF) :- real_mul_wf_aux(X,Y,R,WF).
412
413
414 :- block real_mul_wf_aux(-,-,?,?), real_mul_wf_aux(?,-,-,?), real_mul_wf_aux(-,?,-,?).
415 real_mul_wf_aux(X,Y,R,_WF) :- (X==0.0 ; Y==0.0), !, R=0.0.
416 real_mul_wf_aux(X,Y,R,WF) :-
417 (nonvar(X)
418 -> (nonvar(Y)
419 -> safe_is(R,X*Y,WF)
420 ; allow_constraint_propagation(Kind),
421 safe_is(YY,R/X,WF),
422 confirm_solution_for_expr(YY,X*SOL,SOL,R,Kind,Solutions,WF)
423 -> set_solution_for_expr(Y,Solutions,WF)
424 ; real_mul_wf_aux2(X,Y,R,WF) % delay until Y is known
425 )
426 ; allow_constraint_propagation(Kind),
427 safe_is(XX,R/Y,WF),
428 confirm_solution_for_expr(XX,SOL*Y,SOL,R,Kind,Solutions,WF)
429 -> set_solution_for_expr(X,Solutions,WF)
430 ; real_mul_wf_aux2(X,Y,R,WF) % delay until X is known
431 ).
432
433 :- block real_mul_wf_aux2(-,?,?,?), real_mul_wf_aux2(?,-,?,?).
434 real_mul_wf_aux2(X,Y,R,WF) :- safe_is(R,X*Y,WF).
435
436 % compute sqrt, R is known
437 kernel_square_inverse(X,R,WF) :-
438 (R < 0.0 -> fail
439 ; allow_constraint_propagation(Kind),
440 safe_is(XX,sqrt(R),WF),
441 confirm_solution_for_expr(XX,SOL*SOL,SOL,R,Kind,SolutionsPos,WF),
442 (XX=0.0 -> Solutions = SolutionsPos
443 ; NX is -XX,
444 confirm_solution_for_expr(NX,SOL*SOL,SOL,R,Kind,SolutionsNeg,WF),
445 append(SolutionsPos,SolutionsNeg,Solutions)
446 )
447 -> set_solution_for_expr(X,Solutions,WF)
448 ; real_mul_wf_aux(X,X,R,WF)).
449
450
451 % used for 'RSQRT' external function
452 :- block 'real_square_root_wf'(-,-,?,?).
453 real_square_root_wf(term(floating(X)),term(floating(R)),Span,WF) :-
454 real_sqrt_wf_aux(X,R,Span,WF).
455
456 real_sqrt_wf_aux(X,R,Span,WF) :-
457 use_clpfd_real_solver,!,
458 '$='(R,sqrt(X)),
459 (do_not_double_check_solution -> true ; real_unop_wf_aux('sqrt',X,R,Span,WF)).
460 real_sqrt_wf_aux(X,R,Span,WF) :- block_real_sqrt_wf_aux(X,R,Span,WF).
461
462 :- block block_real_sqrt_wf_aux(-,-,?,?).
463 block_real_sqrt_wf_aux(X,R,_Span,WF) :- var(X),
464 allow_constraint_propagation(Kind),
465 get_preference(find_abort_values,false),
466 (R >= 0.0
467 -> safe_is(XX,R*R,WF),
468 confirm_solution_for_expr(XX,sqrt(SOL),SOL,R,Kind,Solutions,WF)
469 ; Solutions = []), % RSQRT always returns a positive number
470 !,
471 set_solution_for_expr(X,Solutions,WF).
472 block_real_sqrt_wf_aux(X,R,Span,WF) :-
473 real_unop_wf_aux('sqrt',X,R,Span,WF).
474
475
476 :- 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)).
477 real_division_wf(term(floating(X)),term(floating(Y)),term(floating(R)),Span,WF) :-
478 real_div_wf_aux(X,Y,R,Span,WF).
479
480 real_div_wf_aux(X,Y,R,Span,WF) :-
481 use_clpfd_real_solver,!,
482 '$='(R,X / Y),
483 (do_not_double_check_solution -> true ; check_div(X,Y,R,Span,WF)).
484 real_div_wf_aux(X,Y,R,Span,WF) :- block_real_div_wf_aux(X,Y,R,Span,WF).
485
486 :- block block_real_div_wf_aux(-,?,-,?,?), block_real_div_wf_aux(?,-,?,?,?).
487 block_real_div_wf_aux(X,Y,R,Span,WF) :- nonvar(X),!,safe_is(R,X/Y,Span,WF).
488 block_real_div_wf_aux(X,Y,R,Span,WF) :- Y \= 0.0, allow_constraint_propagation(Kind),
489 safe_is(XX,R*Y,Span,WF),
490 confirm_solution_for_expr(XX,SOL/Y,SOL,R,Kind,Solutions,WF),
491 !,
492 set_solution_for_expr(X,Solutions,WF).
493 block_real_div_wf_aux(X,Y,R,Span,WF) :- check_div(X,Y,R,Span,WF).
494
495 :- block check_div(-,?,?,?,?), check_div(?,-,?,?,?).
496 check_div(X,Y,R,Span,WF) :- safe_is(R,X/Y,Span,WF).
497
498 :- block 'real_unary_minus_wf'(-,-,?).
499 real_unary_minus_wf(term(floating(X)),term(floating(R)),WF) :-
500 real_um_wf_aux(X,R,WF).
501
502 real_um_wf_aux(X,R,WF) :-
503 use_clpfd_real_solver,!,
504 '$='(R,-X),
505 (do_not_double_check_solution -> true ; block_real_um_wf_aux2(X,R,WF)).
506 real_um_wf_aux(X,R,WF) :- block_real_um_wf_aux2(X,R,WF).
507
508
509 :- block block_real_um_wf_aux2(-,-,?).
510 block_real_um_wf_aux2(X,R,_) :- nonvar(X),!,R is -X.
511 block_real_um_wf_aux2(X,R,_) :- X is -R.
512
513 % used for 'RABS' external function
514 :- block 'real_absolute_value_wf'(-,-,?).
515 real_absolute_value_wf(term(floating(X)),term(floating(R)),WF) :-
516 real_abs_wf_aux(X,R,WF).
517 real_abs_wf_aux(X,R,WF) :-
518 use_clpfd_real_solver,!,
519 '$='(R,abs(X)),
520 (do_not_double_check_solution -> true ; when(nonvar(X), safe_is(R,abs(X),WF))).
521 real_abs_wf_aux(X,R,WF) :- block_real_abs_wf_aux(X,R,WF).
522
523 :- block block_real_abs_wf_aux(-,-,?).
524 block_real_abs_wf_aux(X,R,_) :- nonvar(X),!, R is abs(X).
525 block_real_abs_wf_aux(X,R,_) :- R = 0.0, !, X=0.0.
526 block_real_abs_wf_aux(_,R,_) :- R < 0.0, !, fail.
527 block_real_abs_wf_aux(X,R,WF) :- MR is -R,
528 Solutions = [R,MR],
529 set_solution_for_expr(X,Solutions,WF).
530
531
532 % used for 'RSIGN' external function
533 :- block 'real_sign_wf'(-,-,?).
534 real_sign_wf(term(floating(X)),term(floating(R)),WF) :-
535 real_sign_wf_aux(X,R,WF).
536
537 real_sign_wf_aux(X,R,WF) :-
538 use_clpfd_real_solver,!,
539 '#='(IR,sign(X)), % Note: it is important to use sign(.) directly here with #=,
540 % otherwise we get an expected 'numeric argument' error
541 convert_int_to_real_aux(IR,R),
542 (do_not_double_check_solution -> true ; block_real_unop_wf(sign,X,R,unknown,WF)).
543 real_sign_wf_aux(X,R,WF) :- block_real_unop_wf(sign,X,R,unknown,WF).
544
545
546 % used for 'ROUND' external function
547 real_round_wf(term(floating(X)),int(R),_WF) :- real_round_aux(X,R).
548
549 real_round_aux(Real,Int) :-
550 use_clpfd_real_solver,!,
551 '#='(Int,round(Real)).
552 real_round_aux(Real,Int) :- block_real_round_aux(Real,Int).
553
554 :- block block_real_round_aux(-,?).
555 block_real_round_aux(Real,Int) :- Int is round(Real).
556
557
558 :- use_module(kernel_waitflags,[add_wd_error/3, add_wd_error_span/4]).
559 % a version of is/2 which catches overflows
560 safe_is(R,Expr,WF) :-
561 safe_is(R,Expr,unknown,WF).
562 safe_is(R,Expr,Span,WF) :-
563 catch(R is Expr,
564 error(evaluation_error(ERR),_),
565 process_evaluation_error(ERR,Expr,Span,WF)).
566
567 process_evaluation_error(float_overflow,Expr,Span,WF) :- !,
568 add_wd_error_span('Float Overflow while computing:',Expr,Span,WF).
569 process_evaluation_error(zero_divisor,Expr,Span,WF) :- !,
570 add_wd_error_span('Division by zero while computing:',Expr,Span,WF).
571 process_evaluation_error(undefined,Expr,Span,WF) :- !,
572 add_wd_error_span('Arithmetic operator undefined while computing:',Expr,Span,WF).
573 process_evaluation_error(_,Expr,Span,WF) :-
574 add_wd_error_span('Unknown evaluation error while computing:',Expr,Span,WF).
575
576
577 % call a Prolog unary artihmetic operator
578 real_unop_wf(OP,X,R,WF) :-
579 real_unop_wf(OP,X,R,unknown,WF).
580
581 real_unop_wf(OP,RX,RR,Span,WF) :-
582 get_real(RX,X,OP,Span,WF), get_real(RR,R,OP,Span,WF),
583 real_unop_wf_aux(OP,X,R,Span,WF).
584
585 unsupported_unary_operator(float_fractional_part). % RFRACTION
586
587 real_unop_wf_aux(OP,X,R,Span,WF) :-
588 use_clpfd_real_solver,
589 \+ unsupported_unary_operator(OP),!,
590 Expr =.. [OP,X],
591 post_real_equal_expr_wf(R,Expr,Span,WF), % TODO: catch WD issues (RASIN(>1.0), ...)
592 (do_not_double_check_solution -> true
593 ; block_real_unop_wf(OP,X,R,Span,WF)). % relevant e.g. for {y|x = RSQRT(y)}=r & x=4.0
594 real_unop_wf_aux(OP,X,R,Span,WF) :- block_real_unop_wf(OP,X,R,Span,WF).
595
596 :- block block_real_unop_wf(-,?,?,?,?), block_real_unop_wf(?,-,?,?,?).
597 block_real_unop_wf(OP,X,R,Span,WF) :-
598 Expr =.. [OP,X],
599 safe_is(R,Expr,Span,WF).
600
601
602 :- use_module(probsrc(tools_strings),[ajoin/2]).
603 :- use_module(kernel_waitflags,[add_error_wf/5]).
604 get_real(Var,Real,_,_,_) :- var(Var),!, Var=term(floating(Real)).
605 get_real(term(F),Real,_,_,_) :- !, F=floating(Real).
606 get_real(Other,_,OP,Span,WF) :-
607 ajoin(['Argument for ',OP,' is not a real number:'],Msg),
608 add_error_wf(kernel_reals,Msg,Other,Span,WF).
609
610 % when called for ** operator (power_of_real) the exponent Y is a natural number represented as a float;
611 % for RPOW Y can be any float
612 real_power_of_wf(RX,RY,RR,Span,WF) :-
613 get_real(RX,X,'RPOW',Span,WF), get_real(RY,Y,'RPOW',Span,WF), get_real(RR,R,'RPOW',Span,WF),
614 (Y==2.0 % often we use built-in ** with integer 2; maybe we should distinguish ** and RPOW
615 -> real_mul_wf_aux0(X,X,R,WF)
616 ; real_power_of_aux(X,Y,R,Span,WF)
617 ).
618
619 :- if(can_use_clpfd_real_solver). %can_use_clpfd_real_solver). % % 4.10 or higher
620 % code still deactivated as i:10..20 & 10.0**i + 1.0 = 10.0**i from test 2430 crashes SICStus, see also test 2562
621 real_power_of_aux(X,Y,R,Span,WF) :- % write(rpow(X,Y,R)),nl,flush_output,
622 real_binop_wf_aux('^',X,Y,R,Span,WF).
623 :- else.
624 :- block real_power_of_aux(-,-,?,?,?), real_power_of_aux(?,-,-,?,?), real_power_of_aux(-,?,-,?,?).
625 % we currently require to know at least the exponent, TODO: in future we could also compute the log
626 real_power_of_aux(X,Y,R,Span,WF) :-
627 var(X),!,
628 real_power_inverse(X,Y,R,Span,WF).
629 real_power_of_aux(X,Y,Res,Span,WF) :-
630 (Y==2.0 % typically we use built-in ** with integer 2; maybe we should distinguish ** and RPOW
631 -> real_mul_wf_aux0(X,X,Res,WF)
632 ; real_power_direct(X,Y,Res,Span,WF)).
633
634 % X is variable, Y and R are known
635 real_power_inverse(X,Y,R,_Span,WF) :- Y=2.0,!,
636 kernel_square_inverse(X,R,WF).
637 real_power_inverse(X,Y,R,_Span,_WF) :- Y=1.0,!, X=R.
638 real_power_inverse(X,Y,R,Span,WF) :-
639 real_power_direct(X,Y,R,Span,WF).
640 % TODO: add more exponents and compute root
641 :- endif.
642
643 :- block real_power_direct(-,?,?,?,?), real_power_direct(?,-,?,?,?).
644 real_power_direct(X,Y,R,Span,WF) :-
645 safe_is(R,'**'(X,Y),Span,WF). % is the same as exp(X,Y) in SICStus, but SWI does not support exp/2
646
647 % call a Prolog binary artihmetic operator: min, max, atan2, log
648 real_binop_wf(OP,RX,RY,RR,Span,WF) :-
649 get_real(RX,X,OP,Span,WF), get_real(RY,Y,OP,Span,WF), get_real(RR,R,OP,Span,WF),
650 real_binop_wf_aux(OP,X,Y,R,Span,WF).
651
652 % atan2, log are not supported by CLP(FD)
653 supported_binary_operator(min).
654 supported_binary_operator(max).
655 supported_binary_operator('^').
656
657 real_binop_wf_aux(OP,X,Y,R,Span,WF) :- supported_binary_operator(OP),
658 use_clpfd_real_solver,!,
659 Expr =.. [OP,X,Y],
660 post_real_equal_expr_wf(R,Expr,Span,WF), % TODO: double check necessary??
661 (do_not_double_check_solution -> true
662 ; block_real_binop_wf(OP,X,Y,R,Span,WF)). % necessary for min/max ?
663 real_binop_wf_aux(OP,X,Y,R,Span,WF) :- block_real_binop_wf(OP,X,Y,R,Span,WF).
664
665 :- block block_real_binop_wf(-,?,?,?,?,?), block_real_binop_wf(?,-,?,?,?,?), block_real_binop_wf(?,?,-,?,?,?).
666 block_real_binop_wf(OP,X,Y,R,Span,WF) :-
667 Expr =.. [OP,X,Y],
668 safe_is(R,Expr,Span,WF).
669
670 real_less_than_wf(term(floating(X)),term(floating(Y)),_WF) :-
671 clpfd_less_real(X,Y).
672 real_less_than_equal_wf(term(floating(X)),term(floating(Y)),_WF) :-
673 clpfd_leq_real(X,Y).
674
675 % call a Prolog binary artihmetic operator, reified version.
676 real_comp_wf(OP,term(floating(X)),term(floating(Y)),R,_WF) :-
677 real_comp_float(OP,X,Y,R).
678
679 clpfd_leq_real(X,Y) :- use_clpfd_real_solver,!,'$=<'(X,Y).
680 clpfd_leq_real(X,Y) :- real_comp_float_blocking('=<',X,Y,pred_true).
681
682 clpfd_less_real(X,Y) :- use_clpfd_real_solver(S),!,
683 ( number(Y), S=precise_float_solver, next_smaller_float(Y,Y1) -> '$=<'(X,Y1)
684 ; number(X), S=precise_float_solver, next_larger_float(X,X1) -> '$=<'(X1,Y)
685 ; '$=<'(X,Y), dif(X,Y)
686 ).
687 clpfd_less_real(X,Y) :- real_comp_float_blocking('<',X,Y,pred_true).
688
689 real_comp_float('=<',X,Y,R) :- use_clpfd_real_solver, !,
690 '#<=>'('$=<'(X,Y),R01),
691 when(?=(X,Y),(X=Y -> R01=1; true)), % overcome limitation of CLP(FD) reification SPRM-21561
692 b_interpreter_check:prop_pred_01(R,R01).
693 real_comp_float('<',X,Y,R) :- use_clpfd_real_solver, !,
694 '#<=>'( '#/\\'('$=<'(X,Y), '#\\'('$='(Y,X))), R01), % does not work due to SPRM-21561
695 when(?=(X,Y),(X=Y -> R01=0; true)), % overcome limitation of CLP(FD) reification SPRM-21561
696 b_interpreter_check:prop_pred_01(R,R01).
697 % TODO: reify =:=, =\\= using post_real_equality
698 real_comp_float(OP,X,Y,R) :- real_comp_float2(OP,X,Y,R).
699
700 :- block real_comp_float2(?,-,?,-), real_comp_float2(?,?,-,-).
701 real_comp_float2('<',X,Y,R) :- R == pred_true,!, clpfd_less_real(X,Y).
702 real_comp_float2('<',X,Y,R) :- R == pred_false,!, clpfd_leq_real(Y,X).
703 %real_comp_float2('=<',X,Y,R) :- R == pred_true,!, clpfd_leq_real(X,Y). % not necessary, as reified
704 %real_comp_float2('=<',X,Y,R) :- R == pred_false,!, clpfd_less_real(Y,X).
705 real_comp_float2(OP,X,Y,R) :- real_comp_float_blocking(OP,X,Y,R).
706
707 :- block real_comp_float_blocking(-,?,?,?), real_comp_float_blocking(?,-,?,?), real_comp_float_blocking(?,?,-,?).
708 real_comp_float_blocking(OP,X,Y,R) :-
709 (call(OP,X,Y) -> R=pred_true ; R=pred_false).
710
711 % -----------------
712
713 is_real_interval_closure_with_known_card([Par],[real],b(Body,_,_),Card) :-
714 get_geq_leq_real_bounds(Body,Par,Low,Up),
715 known_card(Low,Up,Card).
716
717
718 % a version of is_interval_closure_or_integerset that only accepts real intervals
719 is_interval_closure_or_realset(Var,_,_) :- var(Var),!,fail.
720 is_interval_closure_or_realset(global_set('REAL'),minus_inf,inf).
721 is_interval_closure_or_realset(closure([Par],[real],b(Body,_,_)),Low,Up) :-
722 get_geq_leq_real_bounds(Body,Par,Low,Up).
723
724 known_card(Low,_Up,Card) :- Low == minus_inf,!,Card=inf.
725 known_card(_Low,Up,Card) :- Up == inf,!,Card=inf.
726 known_card(Low,Up,Card) :- Low==Up,!,Card=1.
727
728
729
730 :- use_module(custom_explicit_sets,[intersect_intervals_with_inf/6]).
731 % TODO: member
732 get_geq_leq_real_bounds(conjunct(b(LEFT,pred,_),b(RIGHT,pred,_)), Par,Low,Up) :-
733 get_geq_leq_real_bounds(LEFT,Par,From1,To1),
734 get_geq_leq_real_bounds(RIGHT,Par,From2,To2),
735 intersect_intervals_with_inf(From1,To1,From2,To2,Low,Up).
736 get_geq_leq_real_bounds(less_equal_real(b(A,_,_),b(B,_,_)),Par,Low,Up) :- get_rbounds2_leq(A,B,Par,Low,Up).
737 get_geq_leq_real_bounds( less_real(b(A,_,_),b(B,_,_)),Par,Low,Up) :-
738 precise_float_solver_active, %TODO: process strict bounds also in real_solver mode
739 % Note in reals min/max may not be defined for such strict intervals !
740 get_rbounds2_lt(A,B,Par,Low,Up).
741
742 get_rbounds2_leq(identifier(Par),V,Par,minus_inf,X) :- real_value(V,X).
743 get_rbounds2_leq(V,identifier(Par),Par,X,inf) :- real_value(V,X).
744
745 get_rbounds2_lt(identifier(Par),V,Par,minus_inf,VXMinus1) :- real_value(V,VX),block_next_larger_float(VXMinus1,VX).
746 get_rbounds2_lt(V,identifier(Par),Par,VXPlus1,inf) :- real_value(V,VX),block_next_larger_float(VX,VXPlus1).
747
748
749 % try and get a real value from an AST node
750 % see integer_value/2 in custom_explicit_sets
751 real_value(V,_) :- var(V),!, print(var_real_value(V)),nl,fail.
752 real_value(real(X),R) :- !, construct_real_number(X,R).
753 % TODO: unary_minus, ...
754 real_value(value(V),R) :- is_real(V,R).
755 % -----------------
756
757
758 :- use_module(probsrc(kernel_tools),[ground_value_check/2]).
759 :- use_module(probsrc(custom_explicit_sets),[expand_and_convert_to_avl_set/4, expand_custom_set_to_list_wf/5,
760 max_of_explicit_set_wf/3, min_of_explicit_set_wf/3]).
761
762 :- block real_maximum_of_set(-,?,?,?).
763 real_maximum_of_set([],_Res,Span,WF) :- !,
764 add_wd_error_span('max applied to empty set of reals:',[],Span,WF).
765 real_maximum_of_set(avl_set(A),Res,_Span,WF) :- !, max_of_explicit_set_wf(avl_set(A),Res,WF).
766 real_maximum_of_set([RX|T],RRes,Span,WF) :- use_clpfd_real_solver, !,
767 get_real(RX,X,maximum,Span,WF),
768 get_real(RRes,Res,maximum,Span,WF),
769 clpfd_leq_real(X,Res),
770 expand_custom_set_to_list_wf(T,ET,_,real_maximum_of_set,WF),
771 real_max_list(ET,[X],Res,Span,WF).
772 real_maximum_of_set(Set,Res,Span,WF) :-
773 is_interval_closure_or_realset(Set,Low,Up),!,
774 (Up==inf
775 -> add_wd_error_span('maximum of unbounded infinite real set not defined:',Set,Span,WF)
776 ; inf_greater_than_equal_real(Up,Low), % TODO: reify comparison and generate WD error if not? ditto for min and integers
777 is_real(Res,Up)).
778 real_maximum_of_set(Set,Res,Span,WF) :-
779 ground_value_check(Set,Gr),
780 rmax_ground_set(Gr,Set,Res,Span,WF).
781
782 % a comparison also taking inf/minus_inf values into account:
783 inf_greater_than_equal_real(X,Y) :-
784 ((X==inf;Y==minus_inf) -> true ; clpfd_leq_real(Y,X)).
785
786 :- block real_max_list(-,?,?,?,?).
787 real_max_list([],List,Res,_Span,_WF) :- clpfd:maximum(Res,List). % requires new CLP(FD) solver in 4.10
788 real_max_list([RX|T],Acc,Res,Span,WF) :-
789 get_real(RX,X,maximum,Span,WF),
790 clpfd_leq_real(X,Res), % not useful if list skeleton known
791 real_max_list(T,[X|Acc],Res,Span,WF).
792
793 :- block rmax_ground_set(-,?,?,?,?).
794 rmax_ground_set(_,Set,Res,Span,WF) :-
795 expand_and_convert_to_avl_set(Set,ESet,'RMAXIMUM','RMAXIMUM'),
796 (ESet=empty
797 -> add_wd_error_span('max applied to empty set of reals:',Set,Span,WF)
798 ; max_of_explicit_set_wf(avl_set(ESet),Res,WF)
799 ).
800
801 :- block real_minimum_of_set(-,?,?,?).
802 real_minimum_of_set([],_Res,Span,WF) :- !,
803 add_wd_error_span('min applied to empty set of reals:',[],Span,WF).
804 real_minimum_of_set(avl_set(A),Res,_Span,WF) :- !, min_of_explicit_set_wf(avl_set(A),Res,WF).
805 real_minimum_of_set([RX|T],RRes,Span,WF) :- use_clpfd_real_solver, !,
806 get_real(RX,X,minumum,Span,WF),
807 get_real(RRes,Res,minumum,Span,WF),
808 clpfd_leq_real(Res,X),
809 expand_custom_set_to_list_wf(T,ET,_,real_minimum_of_set,WF),
810 real_min_list(ET,[X],Res,Span,WF).
811 real_minimum_of_set(Set,Res,Span,WF) :-
812 is_interval_closure_or_realset(Set,Low,Up),!,
813 (Low == minus_inf
814 -> add_wd_error_span('minimum of unbounded infinite real set not defined:',Set,Span,WF)
815 ; inf_greater_than_equal_real(Up,Low),
816 is_real(Res,Low)).
817 real_minimum_of_set(Set,Res,Span,WF) :-
818 ground_value_check(Set,Gr),
819 rmin_ground_set(Gr,Set,Res,Span,WF).
820
821 :- block real_min_list(-,?,?,?,?).
822 real_min_list([],List,Res,_Span,_WF) :- clpfd:minimum(Res,List). % requires new CLP(FD) solver in 4.10
823 real_min_list([RX|T],Acc,Res,Span,WF) :-
824 get_real(RX,X,minimum,Span,WF),
825 clpfd_leq_real(Res,X), % not useful if list skeleton known
826 real_min_list(T,[X|Acc],Res,Span,WF).
827
828 :- block rmin_ground_set(-,?,?,?,?).
829 rmin_ground_set(_,Set,Res,Span,WF) :-
830 expand_and_convert_to_avl_set(Set,ESet,'RMINIMUM','RMINIMUM'),
831 (ESet=empty
832 -> add_wd_error_span('min applied to empty set of reals:',Set,Span,WF)
833 ; min_of_explicit_set_wf(avl_set(ESet),Res,WF)
834 ).
835
836
837 % not used yet: useful for kernel_strings ?
838 %:- block real_to_string(-,?).
839 %real_to_string(term(floating(I)),S) :- real_to_string2(I,S).
840 %
841 %:- block real_to_string2(-,?).
842 %real_to_string2(Num,Res) :-
843 % number_codes(Num,C),
844 % atom_codes(S,C), Res=string(S).
845
846 % -------------------------
847
848 :- use_module(kernel_objects,[gen_enum_warning_wf/6]).
849 :- use_module(library(random),[random/3]).
850
851
852 enumerate_real_wf(term(floating(F)),EnumWarning,WF) :-
853 enum_real(F,EnumWarning,WF).
854
855 enum_real(F,_,_) :- number(F),!.
856 enum_real(F,EnumWarning,WF) :- use_clpfd_real_solver,!,
857 clpfd_enum_real(F,EnumWarning,WF).
858 enum_real(F,EnumWarning,WF) :-
859 gen_enum_warning_wf('REAL',inf,'"0.0","1.0",...',EnumWarning,unknown,WF),
860 ( F = 0.0 ; F = 1.0 ; F = -1.0
861 ; random(0.0,1.0,F)
862 ; random(-1.0,0.0,F)
863 ; preferences:preference(maxint,MaxInt), random(1.0,MaxInt,F)
864 ; preferences:preference(minint,MinInt), random(MinInt,-1.0,F)
865 ).
866
867 :- if(can_use_clpfd_real_solver). % 4.10 or higher
868 clpfd_enum_real(F,EnumWarning,WF) :-
869 fd_min(F,Min), fd_max(F,Max), % finf .. fsup
870 (number(Min),number(Max) ->
871 (preferences:preference(real_solver_precision,Prec), Prec > 0.0
872 -> gen_enum_warning_wf('REAL',Min..Max,precision(Prec),EnumWarning,unknown,WF),
873 labeling([random, precision(Prec)],[F])
874 ; Min \= Max, preferences:preference(solver_for_reals,real_solver)
875 -> gen_enum_warning_wf('REAL',Min..Max,precision(0),EnumWarning,unknown,WF),
876 labeling([random, precision(5.0E-324)],[F]) % there are infinitely many real numbers
877 ; labeling([random],[F])
878 )
879 ; number(Min), preferences:preference(maxint,MaxInt),
880 MaxR is max(float(MaxInt),Min+1.0) % should we use next float instead of + 1.0? or precision(.) preference?
881 -> gen_enum_warning_wf('REAL',Min..fsup,Min..MaxR,EnumWarning,unknown,WF),
882 '$=<'(F,MaxR), labeling([random, precision(0.5)],[F]) % no need to do exhaustive enumeration; use precision
883 ; number(Max), preferences:preference(minint,MinInt),
884 MinR is min(float(MinInt),Max-1.0) % maybe use previous float or use precision(.) preference
885 -> gen_enum_warning_wf('REAL',finf..Max,MinR..Max,EnumWarning,unknown,WF),
886 '$>='(F,MinR), labeling([random, precision(0.5)],[F]) % no need to do exhaustive enumeration; use precision
887 ; preferences:preference(minint,MinInt), MinR is float(MinInt),
888 preferences:preference(maxint,MaxInt), MaxR is float(MaxInt)
889 -> gen_enum_warning_wf('REAL',finf..fsup,MinR..MaxR,EnumWarning,unknown,WF),
890 fd_batch(['$>='(F,MinR), '$=<'(F,MaxR)]),
891 labeling([random, precision(0.5)],[F]) % no need to do exhaustive enumeration; enum warning anyway
892 ).
893 :- else.
894 clpfd_enum_real(F,_,_) :- add_internal_error('Requires SICStus 4.10:',clpfd_enum_real(F,_,_)).
895 :- endif.
896
897 /*
898
899 A few reference tests with SICS 4.10 real support in CLP(FD)
900
901 | ?- X $>= 1.0, fd_dom(X,D), fd_min(X,Min), fd_max(X,Max).
902 D = 1.0 .. fsup,
903 Min = 1.0,
904 Max = fsup,
905 clpfd:(X in1.0 .. fsup) ?
906
907 | ?- X in 1.0 .. 1.000000000000001, labeling([random],[X]).
908 X = 1.0000000000000004 ? ;
909 X = 1.0000000000000009 ? ;
910 X = 1.000000000000001 ? ;
911 X = 1.0000000000000007 ? ;
912 X = 1.0 ? ;
913 X = 1.0000000000000002 ? ;
914 no
915
916 | ?- X in 1.0 .. 2.0, X + 0.5 $= 2.0, labeling([precision(0.0)],[X]), Delta is (X + 0.5) - 2.0, X + 0.5 =:= 2.0.
917 X = 1.5,
918 Delta = 0.0 ? ;
919 X = 1.5000000000000002,
920 Delta = 0.0 ? ;
921 no
922
923 % enumeration always starts with 1.0 here and not very random:
924 | ?- X $>= -1.0, X $=< 3.0, labeling([random],[X]).
925 X = 1.0 ? ;
926 X = -5.551115123125783E-17 ? ;
927 X = 0.49999999999999994 ? ;
928 X = 0.24999999999999992 ? ;
929 X = 0.3749999999999999 ? ;
930 X = 0.4374999999999999 ? ;
931 X = 0.4062499999999999 ? ;
932 X = 0.3906249999999999 ? ;
933 X = 0.3984374999999999 ? ;
934 X = 0.4023437499999999 ? ;
935 X = 0.4042968749999999 ? ;
936 X = 0.4033203124999999 ? ;
937 X = 0.4038085937499999 ? ;
938 X = 0.4035644531249999 ? ;
939 X = 0.4034423828124999 ? ;
940 X = 0.4035034179687499 ? ;
941 X = 0.4035339355468749 ? ;
942 X = 0.4035186767578124 ? ;
943 X = 0.40351104736328114 ? ;
944
945 | ?- X $= 1.0 + Y, X = 2.0.
946 X = 2.0,
947 clpfd:(Y in0.9999999999999996 .. 1.0000000000000004) ?
948 yes
949
950 | ?- X $= 1.0 + Y, X = 2.0, labeling([random],[Y]), write(Y),nl, (X is 1.0 + Y -> true ; write(ko),nl,fail).
951 1.0
952 X = 2.0,
953 Y = 1.0 ? ;
954 0.9999999999999998
955 ko
956 0.9999999999999997
957 ko
958 0.9999999999999996
959 ko
960 0.9999999999999999
961 X = 2.0,
962 Y = 0.9999999999999999 ? ;
963 1.0000000000000004
964 ko
965 1.0000000000000002
966 X = 2.0,
967 Y = 1.0000000000000002 ? ;
968 no
969 | ?- (X $>= 1.0) #<=> R, X $>= 1.1.
970 R = 1,
971 clpfd:(X in1.1 .. fsup) ?
972 yes
973
974 | ?- X in 11..12, Y $>= float(X).
975 clpfd:(Y in 11.0..fsup),
976 X in 11..12 ?
977 yes
978 | ?- Y in 11.1 .. 12.0, X #>= floor(Y).
979 clpfd:(Y in 11.1..12.0),
980 X in 11..sup ?
981 yes
982 | ?- Y in 11.1 .. 12.0, X #>= ceiling(Y).
983 clpfd:(Y in 11.1..12.0),
984 X in 12..sup ?
985 | ?- X in 1.0 .. 2.5, Y $= sin(X).
986 clpfd:(X in 1.0..2.5),
987 clpfd:(Y in 0.5984721441039559..1.0) ?
988 yes
989 | ?- X in 1.0 .. 2.5, Y $= cos(X).
990 clpfd:(X in 1.0..2.5),
991 clpfd:(Y in-0.8011436155469343..0.54030230586814) ?
992 yes
993 | ?- X in 1.0 .. 2.5, Y $= sinh(X).
994 clpfd:(X in 1.0..2.5),
995 clpfd:(Y in 1.1752011936438012..6.050204481039793) ?
996 yes
997 | ?- X in 1.0 .. 2.5, Y $= tanh(X).
998 clpfd:(X in 1.0..2.5),
999 clpfd:(Y in 0.7615941559557647..0.9866142981514303) ?
1000 yes
1001
1002 | ?- I in 1.0 .. 2.0, J in 1.5 .. 2.1, maximum(Max,[I,J]).
1003 clpfd:(I in 1.0..2.0),
1004 clpfd:(J in 1.5..2.1),
1005 clpfd:(Max in 1.5..2.1) ?
1006 yes
1007 */