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,
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_square_root_wf/4,
20		real_round_wf/3,
21		real_unop_wf/4, real_unop_wf/5,
22		real_binop_wf/6,
23		real_power_of_wf/5,
24		real_less_than_wf/3, real_less_than_equal_wf/3,
25		real_comp_wf/5,
26		real_maximum_of_set/4, real_minimum_of_set/4
27		]).
28
29		% Reals/Floats in ProB are represented by terms of the form term(floating(Number))
30		% in future a proper wrapper such as real(Number) might be created
31
32		:- meta_predicate real_comp_wf(2,-,-,-,-).
33		:- meta_predicate real_comp_wf_aux(2,-,-,-,-).
34
35		:- use_module(module_information,[module_info/2]).
36		:- module_info(group,kernel).
37		:- module_info(description,'This module provides (external) functions to manipulate B reals and floats.').
38
39		:- use_module(error_manager).
40		:- use_module(self_check).
41		:- use_module(library(lists)).
42		:- use_module(debug).
43
44		:- use_module(kernel_objects,[exhaustive_kernel_check_wf/2,exhaustive_kernel_check_wf/3]).
45		:- use_module(extension('counter/counter'),
46		[next_smaller_abs_float/2, next_smaller_float/2,next_larger_float/2, next_float_twoards/3,
47		smallest_float/1, largest_float/1, smallest_abs_float/1]).
48
49		% used to construct a Value from real(Atom) AST node
50		construct_real(Atom,term(floating(Float))) :-
51		atom_codes(Atom,C),
52		number_codes(Nr,C),
53		Float is float(Nr). % make sure we store a float; not required for AST literals
54
55		construct_real_number(Atom,Nr) :-
56		atom_codes(Atom,C),
57		number_codes(Nr,C).
58
59		construct_negative_real(Atom,term(floating(Float))) :-
60		atom_codes(Atom,C),
61		number_codes(Nr,C),
62		Float is -float(Nr).
63
64		is_real_wf(X,_WF) :- is_real(X,_).
65		is_real(X) :- is_real(X,_).
66
67		is_real(term(floating(Nr)),Nr).
68		% TO DO: other formats
69
70		is_float_wf(X,_WF) :- is_float(X).
71
72		is_float(term(floating(_))).
73
74		is_not_real(_) :- fail.
75		is_not_float(_) :- fail.
76
77		is_ground_real(term(floating(Nr))) :- number(Nr).
78
79		is_largest_positive_float(term(floating(Nr))) :- largest_float(Nr).
80		is_smallest_positive_float(term(floating(Nr))) :- smallest_abs_float(Nr).
81
82		is_next_larger_float(term(floating(Nr)),term(floating(Next))) :-
83		block_next_larger_float(Nr,Next).
84		:- block block_next_larger_float(-,-).
85		block_next_larger_float(Nr,Next) :- nonvar(Nr),!,next_larger_float(Nr,Next).
86		block_next_larger_float(Nr,Next) :- next_smaller_float(Next,Nr).
87
88		is_next_smaller_float(term(floating(Nr)),term(floating(Next))) :-
89		block_next_larger_float(Next,Nr).
90
91
92		% convert an integer to a real, corresponds to the real(.) function in Atelier-B; convert_real as AST
93		convert_int_to_real(int(X),term(floating(R))) :-
94		convert_int_to_real_aux(X,R).
95
96		:- block convert_int_to_real_aux(-,-).
97		convert_int_to_real_aux(Int,Real) :- number(Int),!,
98		Real is float(Int).
99		convert_int_to_real_aux(Int,Real) :- Int1 is floor(Real),
100		Real is float(Int1),
101		Int1=Int.
102
103
104		% convert_int_floor as AST, Atelier-B floor(.)
105		real_floor(term(floating(R)),int(X)) :- real_floor_aux(R,X).
106		:- block real_floor_aux(-,?).
107		real_floor_aux(Real,Int) :- Int is floor(Real).
108
109		% convert_int_ceiling as AST, Atelier-B ceiling(.)
110		real_ceiling(term(floating(R)),int(X)) :- real_ceiling_aux(R,X).
111		:- block real_ceiling_aux(-,?).
112		real_ceiling_aux(Real,Int) :- Int is ceiling(Real).
113
114		:- 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)).
115		:- 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)).
116		:- 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)).
117
118		real_addition_wf(term(floating(X)),term(floating(Y)),term(floating(R)),WF) :-
119		real_add_wf_aux(X,Y,R,WF).
120		:- block real_add_wf_aux(-,-,?,?), real_add_wf_aux(?,-,-,?), real_add_wf_aux(-,?,-,?).
121		real_add_wf_aux(X,Y,R,WF) :-
122		(nonvar(X)
123		-> (nonvar(Y)
124		-> safe_is(R,X+Y,WF)
125		; allow_constraint_propagation(Kind),
126		safe_is(YY,R-X,WF),
127		confirm_solution_for_expr(YY,X+SOL,SOL,R,Kind,Solutions,WF)
128		-> set_solution_for_expr(Y,Solutions,WF)
129		; real_add_wf_aux2(X,Y,R,WF) % delay until Y is known
130		)
131		; allow_constraint_propagation(Kind),
132		safe_is(XX,R-Y,WF),
133		confirm_solution_for_expr(XX,SOL+Y,SOL,R,Kind,Solutions,WF)
134		-> set_solution_for_expr(X,Solutions,WF)
135		; real_add_wf_aux2(X,Y,R,WF) % delay until X is known
136		).
137
138		:- block real_add_wf_aux2(-,?,?,?), real_add_wf_aux2(?,-,?,?).
139		real_add_wf_aux2(X,Y,R,WF) :- safe_is(R,X+Y,WF).
140
141
142		% confirm_solution_for_expr(XX,ExprX,X,R : check if XX=X is a solution for R = ExprX, where ExprX contains variable X
143		% depending on solver setting : compute a list of candidate solutions around XX
144		confirm_solution_for_expr(XX,_,_,_,no_checking,Solutions,_) :- !, Solutions=XX.
145		confirm_solution_for_expr(XX,ExprX,X,R,simple_checking,Solutions,WF) :- !,
146		check_sol_direct(XX,ExprX,X,R,WF), Solutions=XX. % single solution
147		confirm_solution_for_expr(XX,ExprX,X,R,check_unique_solution,Solutions,WF) :- !,
148		check_sol_direct(XX,ExprX,X,R,WF), Solutions=XX, % single solution
149		% Now check if there are other solutions:
150		next_larger_float(XX,XN),
151		(X=XN,safe_is(R,ExprX,WF)
152		-> debug_format(19,'Reject ambiguous larger solution {~w,~w} for ~w = ~w~n',[XX,XN,ExprX,R]),fail
153		; true
154		),
155		next_smaller_float(XX,XP),
156		(X=XP,safe_is(R,ExprX,WF)
157		-> debug_format(19,'Reject ambiguous smaller solution {~w,~w} for ~w = ~w~n',[XX,XP,ExprX,R]),fail
158		; true).
159		confirm_solution_for_expr(XX,ExprX,X,R,check_small_interval_solution,Interval,WF) :-
160		get_solutions_around(XX,ExprX,X,R,Interval,WF).
161
162		:- use_module(kernel_waitflags,[get_wait_flag/4]).
163		% set the solution based on result of confirm_solution_for_expr
164		set_solution_for_expr(X,XX,_WF) :- number(XX),!,X=XX. % deterministic constraint propagation found precise solution
165		set_solution_for_expr(X,[XX],_WF) :- !,X=XX. % ditto
166		set_solution_for_expr(_,[],_WF) :- !, fail. % no solutions
167		set_solution_for_expr(X,Interval,WF) :-
168		length(Interval,Len), % we found an interval list of possible solutions
169		get_wait_flag(Len,kernel_reals,WF,LWF),
170		set_sol_aux(X,Interval,LWF).
171		:- block set_sol_aux(-,?,-).
172		set_sol_aux(X,Interval,_) :-
173	?	member(X,Interval).
174
175
176		% get a sorted list of solutions around the candidate XX
177		get_solutions_around(XX,ExprX,X,R,Interval,WF) :-
178		(check_sol_direct(XX,ExprX,X,R,WF) % XX is a solution
179		-> get_preference(solver_strength,SS),
180		Limit is 10+SS, % number of solutions we are willing to find
181		next_larger_float(XX,XN),
182		find_larger_bounds(XN,Limit,ExprX,X,R,WF,LargerSols,[]),
183		next_smaller_float(XX,XP),
184		find_lower_bounds(XP,Limit,ExprX,X,R,WF,Interval,[XX|LargerSols])
185		; (X='X',debug_format(19,'Reject non-solution ~w for ~w = ~w~n',[XX,ExprX,R]),fail ; true),
186		%next_smaller_float(XX,XP), (check_sol_direct(XP,ExprX,X,R,WF) -> write(xp(XP)),nl ; write(nxp(XP)),nl),
187		%next_larger_float(XX,XN), (check_sol_direct(XN,ExprX,X,R,WF) -> write(xn(XN)),nl ; write(nxp(XN)),nl),
188		Interval = [] % no floating point solutions; TODO: do we need to check smaller/larger values??
189		).
190
191		% e.g. for x+1.0 = 3.0 we have 2.0 as solution and the previous float 1.9999999999999998, nothing else
192		% e.b. for x+2.0 = 3.0 we have four solutions {0.9999999999999998,0.9999999999999999,1.0,1.0000000000000002}
193		% card({x|x + 2.999 = 3.0}) = 2049
194		% card({x|x + 2.9999 = 3.0}) = 32769
195
196		% increment floats until we hit the first non-solution or we reach the limit
197		% TODO: try and compute this bound, rather then single-stepping
198		% TODO: enumerate the first solution and then generate an enum warning if we reach limit
199		find_larger_bounds(XP,Lim,ExprX,X,R,WF,[XP|Sols1],Sols2) :-
200		check_sol_direct(XP,ExprX,X,R,WF),!, Lim>0, L1 is Lim-1,
201		next_larger_float(XP,XP2),
202		find_larger_bounds(XP2,L1,ExprX,X,R,WF,Sols1,Sols2).
203		find_larger_bounds(_Bound,_Lim,_ExprX,_X,_R,_WF,Sols,Sols).
204
205		% decrement floats until we hit the first non-solution or we reach the limit
206		find_lower_bounds(XP,Lim,ExprX,X,R,WF,[XP|Sols1],Sols2) :-
207		check_sol_direct(XP,ExprX,X,R,WF),!, Lim>0, L1 is Lim-1,
208		next_smaller_float(XP,XP2),
209		find_lower_bounds(XP2,L1,ExprX,X,R,WF,Sols1,Sols2).
210		find_lower_bounds(_Bound,_Lim,_ExprX,_X,_R,_WF,Sols,Sols).
211
212		check_sol_direct(XX,ExprX,X,R,WF) :- !,
213		% (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
214		(\+((X=XX,safe_is(R,ExprX,WF))) -> fail ; true).
215
216
217		:- use_module(probsrc(preferences), [get_preference/2]).
218
219		allow_constraint_propagation(Kind) :-
220		get_preference(solver_for_reals,Solver),
221		get_solver_kind(Solver,Kind).
222
223		get_solver_kind(aggressive_float_solver,no_checking). % like CLP(R): do not check propagated solution
224		get_solver_kind(float_solver,simple_checking). % check propagated solution, but do not check uniqueness
225		%get_solver_kind(precise_float_solver,check_unique_solution).
226		get_solver_kind(precise_float_solver,check_small_interval_solution). % perform limited constraint propagation, fail directly if no solution exists
227
228		% we could have various other options:
229		% aggressive: do it like CLP(R), which can generate non-solutions
230		% | ?- {X+10.0e10 = 1.0e-9}, write(sol(X)),nl, {X+10.0e10 = 1.0e-9}.
231		% prints sol(-1.0E+11) but then fails
232		% conservative: do the precision check and also check that there are no other solutions (not done yet)
233		% e.g., X+10000000000.0 = 10000000000.0 & (X=0.000000000001 or X=2.0) is FALSE with the float_solver, but TRUE with none
234
235
236		:- 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)).
237		real_subtraction_wf(term(floating(X)),term(floating(Y)),term(floating(R)),WF) :-
238		real_sub_wf_aux(X,Y,R,WF).
239		:- block real_sub_wf_aux(-,-,?,?), real_sub_wf_aux(?,-,-,?), real_sub_wf_aux(-,?,-,?).
240		real_sub_wf_aux(X,Y,R,WF) :-
241		(nonvar(X)
242		-> (nonvar(Y)
243		-> safe_is(R,X-Y,WF)
244		; allow_constraint_propagation(Kind),
245		safe_is(YY,X-R,WF),
246		confirm_solution_for_expr(YY,X-SOL,SOL,R,Kind,Solutions,WF)
247		-> set_solution_for_expr(Y,Solutions,WF) % deterministic constraint propagation found precise solution
248		; real_sub_wf_aux2(X,Y,R,WF) % delay until Y is known
249		)
250		; allow_constraint_propagation(Kind),
251		safe_is(XX,R+Y,WF),
252		confirm_solution_for_expr(XX,SOL-Y,SOL,R,Kind,Solutions,WF)
253		-> set_solution_for_expr(X,Solutions,WF) % deterministic constraint propagation found precise solution
254		; real_sub_wf_aux2(X,Y,R,WF) % delay until X is known
255		).
256
257		:- block real_sub_wf_aux2(-,?,?,?), real_sub_wf_aux2(?,-,?,?).
258		real_sub_wf_aux2(X,Y,R,WF) :- safe_is(R,X-Y,WF).
259
260
261		:- 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)).
262		real_multiplication_wf(term(floating(X)),term(floating(Y)),term(floating(R)),WF) :-
263		real_mul_wf_aux0(X,Y,R,WF).
264
265		:- block real_mul_wf_aux0(-,-,-,?).
266		real_mul_wf_aux0(X,Y,R,_WF) :- (X==0.0 ; Y==0.0), !, R=0.0.
267		real_mul_wf_aux0(X,Y,R,WF) :- var(X),X==Y, !, % X*X = R -> compute sqrt
268		kernel_square_inverse(X,R,WF).
269		real_mul_wf_aux0(X,Y,R,WF) :- real_mul_wf_aux(X,Y,R,WF).
270
271		:- block real_mul_wf_aux(-,-,?,?), real_mul_wf_aux(?,-,-,?), real_mul_wf_aux(-,?,-,?).
272		real_mul_wf_aux(X,Y,R,_WF) :- (X==0.0 ; Y==0.0), !, R=0.0.
273		real_mul_wf_aux(X,Y,R,WF) :-
274		(nonvar(X)
275		-> (nonvar(Y)
276		-> safe_is(R,X*Y,WF)
277		; allow_constraint_propagation(Kind),
278		safe_is(YY,R/X,WF),
279		confirm_solution_for_expr(YY,X*SOL,SOL,R,Kind,Solutions,WF)
280		-> set_solution_for_expr(Y,Solutions,WF)
281		; real_mul_wf_aux2(X,Y,R,WF) % delay until Y is known
282		)
283		; allow_constraint_propagation(Kind),
284		safe_is(XX,R/Y,WF),
285		confirm_solution_for_expr(XX,SOL*Y,SOL,R,Kind,Solutions,WF)
286		-> set_solution_for_expr(X,Solutions,WF)
287		; real_mul_wf_aux2(X,Y,R,WF) % delay until X is known
288		).
289
290		:- block real_mul_wf_aux2(-,?,?,?), real_mul_wf_aux2(?,-,?,?).
291		real_mul_wf_aux2(X,Y,R,WF) :- safe_is(R,X*Y,WF).
292
293		% compute sqrt, R is known
294		kernel_square_inverse(X,R,WF) :-
295		(R < 0.0 -> fail
296		; allow_constraint_propagation(Kind),
297		safe_is(XX,sqrt(R),WF),
298		confirm_solution_for_expr(XX,SOL*SOL,SOL,R,Kind,SolutionsPos,WF),
299		(XX=0.0 -> Solutions = SolutionsPos
300		; NX is -XX,
301		confirm_solution_for_expr(NX,SOL*SOL,SOL,R,Kind,SolutionsNeg,WF),
302		append(SolutionsPos,SolutionsNeg,Solutions)
303		)
304		-> set_solution_for_expr(X,Solutions,WF)
305		; real_mul_wf_aux(X,X,R,WF)).
306
307
308		:- block 'real_square_root_wf'(-,-,?,?).
309		real_square_root_wf(term(floating(X)),term(floating(R)),Span,WF) :-
310		real_sqrt_wf_aux(X,R,Span,WF).
311		:- block real_sqrt_wf_aux(-,-,?,?).
312		real_sqrt_wf_aux(X,R,_Span,WF) :- var(X),
313		allow_constraint_propagation(Kind),
314		get_preference(find_abort_values,false),
315		(R >= 0.0
316		-> safe_is(XX,R*R,WF),
317		confirm_solution_for_expr(XX,sqrt(SOL),SOL,R,Kind,Solutions,WF)
318		; Solutions = []), % RSQRT always returns a positive number
319		!,
320		set_solution_for_expr(X,Solutions,WF).
321		real_sqrt_wf_aux(X,R,Span,WF) :-
322		real_unop_wf_aux('sqrt',X,R,Span,WF).
323
324		:- 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)).
325		real_division_wf(term(floating(X)),term(floating(Y)),term(floating(R)),Span,WF) :-
326		real_div_wf_aux(X,Y,R,Span,WF).
327		:- block real_div_wf_aux(-,?,-,?,?), real_div_wf_aux(?,-,?,?,?).
328		real_div_wf_aux(X,Y,R,Span,WF) :- nonvar(X),!,safe_is(R,X/Y,Span,WF).
329		real_div_wf_aux(X,Y,R,Span,WF) :- Y \= 0.0, allow_constraint_propagation(Kind),
330		safe_is(XX,R*Y,Span,WF),
331		confirm_solution_for_expr(XX,SOL/Y,SOL,R,Kind,Solutions,WF),
332		!,
333		set_solution_for_expr(X,Solutions,WF).
334		real_div_wf_aux(X,Y,R,Span,WF) :- when(nonvar(X),safe_is(R,X/Y,Span,WF)).
335
336		:- block 'real_unary_minus_wf'(-,-,?).
337		real_unary_minus_wf(term(floating(X)),term(floating(R)),WF) :-
338		real_um_wf_aux(X,R,WF).
339		:- block real_um_wf_aux(-,-,?).
340		real_um_wf_aux(X,R,_) :- nonvar(X),!,R is -X.
341		real_um_wf_aux(X,R,_) :- X is -R.
342
343		:- block 'real_absolute_value_wf'(-,-,?).
344		real_absolute_value_wf(term(floating(X)),term(floating(R)),WF) :-
345		real_abs_wf_aux(X,R,WF).
346		:- block real_abs_wf_aux(-,-,?).
347		real_abs_wf_aux(X,R,_) :- nonvar(X),!, R is abs(X).
348		real_abs_wf_aux(X,R,_) :- R = 0.0, !, X=0.0.
349		real_abs_wf_aux(_,R,_) :- R < 0.0, !, fail.
350		real_abs_wf_aux(X,R,WF) :- MR is -R,
351		Solutions = [R,MR],
352		set_solution_for_expr(X,Solutions,WF).
353
354		real_round_wf(term(floating(X)),int(R),_WF) :-
355		when(nonvar(X), R is round(X)).
356
357		:- use_module(kernel_waitflags,[add_wd_error/3, add_wd_error_span/4]).
358		% a version of is/2 which catches overflows
359		safe_is(R,Expr,WF) :-
360		safe_is(R,Expr,unknown,WF).
361		safe_is(R,Expr,Span,WF) :-
362		catch(R is Expr,
363		error(evaluation_error(ERR),_),
364		process_evaluation_error(ERR,Expr,Span,WF)).
365
366		process_evaluation_error(float_overflow,Expr,Span,WF) :- !,
367		add_wd_error_span('Float Overflow while computing:',Expr,Span,WF).
368		process_evaluation_error(zero_divisor,Expr,Span,WF) :- !,
369		add_wd_error_span('Division by zero while computing:',Expr,Span,WF).
370		process_evaluation_error(undefined,Expr,Span,WF) :- !,
371		add_wd_error_span('Arithmetic operator undefined while computing:',Expr,Span,WF).
372		process_evaluation_error(_,Expr,Span,WF) :-
373		add_wd_error_span('Unknown evaluation error while computing:',Expr,Span,WF).
374
375
376		% call a Prolog unary artihmetic operator
377		real_unop_wf(OP,X,R,WF) :-
378		real_unop_wf(OP,X,R,unknown,WF).
379
380		real_unop_wf(OP,RX,RR,Span,WF) :-
381		get_real(RX,X,OP,Span,WF), get_real(RR,R,OP,Span,WF),
382		real_unop_wf_aux(OP,X,R,Span,WF).
383		:- block real_unop_wf_aux(-,?,?,?,?), real_unop_wf_aux(?,-,?,?,?).
384		real_unop_wf_aux(OP,X,R,Span,WF) :-
385		Expr =.. [OP,X],
386		safe_is(R,Expr,Span,WF).
387
388		:- use_module(probsrc(tools_strings),[ajoin/2]).
389		:- use_module(kernel_waitflags,[add_error_wf/5]).
390		get_real(Var,Real,_,_,_) :- var(Var),!, Var=term(floating(Real)).
391		get_real(term(F),Real,_,_,_) :- !, F=floating(Real).
392		get_real(Other,_,OP,Span,WF) :-
393		ajoin(['Argument for ',OP,' is not a real number:'],Msg),
394		add_error_wf(kernel_reals,Msg,Other,Span,WF).
395
396		% when called for ** operator the exponent Y is a natural number represented as a float; for RPOW Y can be any float
397		real_power_of_wf(RX,RY,RR,Span,WF) :-
398		get_real(RX,X,'exp',Span,WF), get_real(RY,Y,'exp',Span,WF), get_real(RR,R,'exp',Span,WF),
399		real_power_of_aux(X,Y,R,Span,WF).
400
401		:- block real_power_of_aux(-,-,?,?,?), real_power_of_aux(?,-,-,?,?), real_power_of_aux(-,?,-,?,?).
402		% we currently require to know at least the exponent, TODO: in future we could also compute the log
403		real_power_of_aux(X,Y,R,Span,WF) :-
404		var(X),!,
405		real_power_inverse(X,Y,R,Span,WF).
406		real_power_of_aux(X,Y,Res,Span,WF) :-
407		real_power_direct(X,Y,Res,Span,WF).
408
409		% X is variable, Y and R are known
410		real_power_inverse(X,Y,R,_Span,WF) :- Y=2.0,!,
411		kernel_square_inverse(X,R,WF).
412		real_power_inverse(X,Y,R,_Span,_WF) :- Y=1.0,!, X=R.
413		real_power_inverse(X,Y,R,Span,WF) :-
414		real_power_direct(X,Y,R,Span,WF).
415		% TODO: add more exponents and compute root
416
417		:- block real_power_direct(-,?,?,?,?), real_power_direct(?,-,?,?,?).
418		real_power_direct(X,Y,R,Span,WF) :-
419		safe_is(R,exp(X,Y),Span,WF).
420
421		% call a Prolog binary artihmetic operator
422		real_binop_wf(OP,RX,RY,RR,Span,WF) :-
423		get_real(RX,X,OP,Span,WF), get_real(RY,Y,OP,Span,WF), get_real(RR,R,OP,Span,WF),
424		real_binop_wf_aux(OP,X,Y,R,Span,WF).
425		:- block real_binop_wf_aux(-,?,?,?,?,?), real_binop_wf_aux(?,-,?,?,?,?), real_binop_wf_aux(?,?,-,?,?,?).
426		real_binop_wf_aux(OP,X,Y,R,Span,WF) :-
427		Expr =.. [OP,X,Y],
428		safe_is(R,Expr,Span,WF).
429
430		real_less_than_wf(X,Y,WF) :- real_comp_wf('<',X,Y,pred_true,WF).
431		real_less_than_equal_wf(X,Y,WF) :- real_comp_wf('=<',X,Y,pred_true,WF).
432
433
434		% call a Prolog bianry artihmetic operator
435		real_comp_wf(OP,term(floating(X)),term(floating(Y)),R,WF) :-
436		real_comp_wf_aux(OP,X,Y,R,WF).
437		:- block real_comp_wf_aux(-,?,?,?,?), real_comp_wf_aux(?,-,?,?,?), real_comp_wf_aux(?,?,-,?,?).
438		real_comp_wf_aux(OP,X,Y,R,_) :-
439		(call(OP,X,Y) -> R=pred_true ; R=pred_false).
440
441		% -----------------
442
443		:- use_module(probsrc(kernel_tools),[ground_value_check/2]).
444		real_maximum_of_set(Set,Res,Span,WF) :-
445		ground_value_check(Set,Gr),
446		rmax_set(Gr,Set,Res,Span,WF).
447
448		:- block rmax_set(-,?,?,?,?).
449		rmax_set(_,Set,Res,Span,WF) :-
450		custom_explicit_sets:expand_and_convert_to_avl_set(Set,ESet,'RMAXIMUM','RMAXIMUM'),
451		(ESet=empty
452		-> add_wd_error_span('max applied to empty set of reals:',Set,Span,WF)
453		; custom_explicit_sets:max_of_explicit_set_wf(avl_set(ESet),Res,WF)
454		).
455
456		real_minimum_of_set(Set,Res,Span,WF) :-
457		ground_value_check(Set,Gr),
458		rmin_set(Gr,Set,Res,Span,WF).
459
460		:- block rmin_set(-,?,?,?,?).
461		rmin_set(_,Set,Res,Span,WF) :-
462		custom_explicit_sets:expand_and_convert_to_avl_set(Set,ESet,'RMINIMUM','RMINIMUM'),
463		(ESet=empty
464		-> add_wd_error_span('min applied to empty set of reals:',Set,Span,WF)
465		; custom_explicit_sets:min_of_explicit_set_wf(avl_set(ESet),Res,WF)
466		).
467
468
469		% not used yet: useful for kernel_strings ?
470		%:- block real_to_string(-,?).
471		%real_to_string(term(floating(I)),S) :- real_to_string2(I,S).
472		%
473		%:- block real_to_string2(-,?).
474		%real_to_string2(Num,Res) :-
475		% number_codes(Num,C),
476		% atom_codes(S,C), Res=string(S).
477