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