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