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