23
29evalNode(Op, P, Is, R) :-
30 g_inc('clpBNR:evalNode'), 31 narrowing_op(Op, P, Is, R),
32 !.
33evalNode(_Op, _, _, _):-
34 g_inc('clpBNR:evalNodeFail'), 35 debugging(clpBNR), 36 current_node_(C),
37 debug_clpBNR_("** fail ** ~w.",C),
38 fail.
39
40sandbox:safe_global_variable('clpBNR:evalNode').
41sandbox:safe_global_variable('clpBNR:evalNodeFail').
42
43clpStatistics :-
44 g_assign('clpBNR:evalNode',0),
45 g_assign('clpBNR:evalNodeFail',0),
46 fail. 47
48clpStatistic(narrowingOps(C)) :- g_read('clpBNR:evalNode',C).
49
50clpStatistic(narrowingFails(C)) :- g_read('clpBNR:evalNodeFail',C).
51
55goal_expansion(non_empty(L,H), 1 =\= cmpr(L,H)). 56
57non_empty((L,H)) :- non_empty(L,H).
58non_empty(L,H) :- 1 =\= cmpr(L,H).
59
63intCase(C, (L,H)) :-
64 ( L>=0 -> C=p 65 ; H=<0 -> C=n 66 ; C=s 67 ).
68
72^((Xl,Xh), (Yl,Yh), (Zl,Zh)) :-
73 Zl is maxr(Xl, Yl), Zh is minr(Xh, Yh), 74 non_empty(Zl,Zh).
75
79union_((Xl,Xh),(Yl,Yh),(Zl,Zh)) :-
80 Zl is minr(Xl,Yl), Zh is maxr(Xh,Yh). 81
87to_float((L,H),(FL,FH)) :-
88 ( float(L) -> FL = L ; FL is roundtoward(float(L), to_negative) ),
89 ( float(H) -> FH = H ; FH is roundtoward(float(H), to_positive) ).
90
95:- discontiguous clpBNR:narrowing_op/4. 96:- style_check(-singleton). 99
101
102narrowing_op(integral, P, $(X), $(NewX)) :-
103 integral_(X,NewX).
104
105integral_((Xl,Xh),(NXl,NXh)) :-
106 NXl is ceiling(Xl), NXh is floor(Xh),
107 non_empty(NXl,NXh).
108
110
114
115narrowing_op(int, P, $(Z,X), $(NewZ,NewX)) :-
116 int_(Z,X,P,NewZ,NewX).
117
118int_(Z,X,p,NewZ,X) :- 119 X = (Xl,Xh),
120 (0 is cmpr(Xl,Xh) 121 -> (integer(Xl) -> B = 1 ; B = 0) 122 ; ceiling(Xl) > floor(Xh), 123 B = 0
124 ),
125 !, 126 ^(Z,(B,B),NewZ).
127
128int_(Z,X,_P,NewZ,NewX) :- 129 (Z = (1,1)
130 -> integral_(X,NewX), 131 NewZ = Z
132 ; NewX = X, 133 ^(Z,(0,1),NewZ) 134 ).
135
137
139
140narrowing_op(eq, P, $(Z,X,Y), Out) :-
141 (^(X,Y,New) 142 -> (Z = (1,1)
143 -> Out = $(Z,New,New),
144 New = (L,H),
145 (0 is cmpr(L,H) -> P = p ; true) 146 ; 147 ( X = (Xl,Xh), Y = (Yl,Yh),
148 0 is cmpr(Xl,Xh) \/ cmpr(Yl,Yh) \/ cmpr(Xl,Yl) 149 -> ^(Z,(1,1),NewZ), 150 P = p 151 ; ^(Z,(0,1),NewZ) 152 ),
153 Out = $(NewZ,X,Y)
154 )
155 ; 156 ^(Z, (0,0), NewZ), 157 Out = $(NewZ,X,Y), P = p 158 ).
159
161
163
164narrowing_op(ne, P, $(Z,X,Y), Out) :-
165 (^(X,Y,(L,H))
166 -> (Z = (1,1)
167 -> ne_int_(X,Y,NewX), 168 ne_int_(Y,NewX,NewY), 169 Out = $(Z,NewX,NewY)
170 ; (X = Y, L = H -> B = (0,0), P = p ; B = (0,1)),
171 ^(Z,B,NewZ),
172 Out = $(NewZ,X,Y)
173 )
174 ; 175 ^(Z, (1,1), NewZ), 176 P = p, 177 Out = $(NewZ,X,Y)
178 ).
179
182ne_int_((X,H), (X,X), (NewL,H)) :- !, 183 NewL is X+1, NewL=<H.
184ne_int_((L,X), (X,X), (L,NewH)) :- !, 185 NewH is X-1, L=<NewH.
186ne_int_(Y, _X, Y). 187
189
191
192narrowing_op(le, P, $(Z, X, Y), Out) :-
193 X = (Xl,Xh), Y = (Yl,Yh),
194 (non_empty(Xh,Yl) 195 -> ^(Z,(1,1),NewZ), P = p, Out = $(NewZ,X,Y) 196 ; (-1 is cmpr(Yh,Xl) 197 -> ^(Z,(0,0),NewZ), P = p, Out = $(NewZ,X,Y) 198 ; (Z = (1,1) 199 -> NXh is minr(Xh,Yh), 200 NYl is maxr(Xl,Yl), 201 Out = $(Z,(Xl,NXh),(NYl,Yh))
202 ; ^(Z,(0,1),NewZ), 203 Out = $(NewZ,X,Y)
204 )
205 )
206 ).
207
209
211
212narrowing_op(lt, P, $(Z, X, Y), Out) :-
213 X = (Xl,Xh), Y = (Yl,Yh),
214 (-1 is cmpr(Xh,Yl) 215 -> ^(Z,(1,1),NewZ), P = p, Out = $(NewZ,X,Y) 216 ; (1 =\= cmpr(Yh,Xl) 217 -> ^(Z,(0,0),NewZ), P = p, Out = $(NewZ,X,Y) 218 ; (Z = (1,1) 219 -> next_lt_(Yh,-1,YhD), 220 NXh is minr(Xh,YhD), 221 next_lt_(Xl,1,XlU), 222 NYl is maxr(XlU,Yl), 223 (NXh < NYl -> P=p ; true), 224 Out = $(Z,(Xl,NXh),(NYl,Yh))
225 ; ^(Z,(0,1),NewZ), 226 Out = $(NewZ,X,Y)
227 )
228 )
229 ).
230
233next_lt_( 1.0Inf, _, 1.0Inf) :- !.
234next_lt_(-1.0Inf, _, -1.0Inf) :- !.
235next_lt_(V, -1, NV) :- NV is maxr(V-1,nexttoward(V,-1.0Inf)). 236next_lt_(V, 1, NV) :- NV is minr(V+1,nexttoward(V, 1.0Inf)).
237
239
241
242narrowing_op(in, P, $(Z, X, Y), $(NewZ, NewX, Y)):-
243 244 (^(X,Y,X1) 245 -> (Z == (1,1)
246 -> NewX = X1 247 ; NewX = X, 248 ^(Z,(0,1),NewZ) 249 )
250 ; ^(Z,(0,0),NewZ), 251 NewX = X, 252 P=p 253 ).
254
256
258
259narrowing_op(add, _, $((Zl,Zh), (Xl,Xh), (Yl,Yh)), $((NZl,NZh), (NXl,NXh), (NYl,NYh))) :-
260 NZl is maxr(Zl, roundtoward( Xl+Yl, to_negative)), 261 NZh is minr(Zh, roundtoward( Xh+Yh, to_positive)),
262 non_empty(NZl,NZh),
263 264 265 NXl is maxr(Xl, roundtoward(NZl+(-Yh), to_negative)), 266 NXh is minr(Xh, roundtoward(NZh+(-Yl), to_positive)),
267 non_empty(NXl,NXh),
268 NYl is maxr(Yl, roundtoward(NZl+(-NXh), to_negative)), 269 NYh is minr(Yh, roundtoward(NZh+(-NXl), to_positive)),
270 non_empty(NYl,NYh).
271
273
275
276narrowing_op(mul, _, $(Z,X,Y), $(NewZ, NewX, NewY)) :-
277 intCase(Cx,X),
278 intCase(Cy,Y),
279 multCase(Cx,Cy,X,Y,Z,NewZ), 280 intCase(CNz,NewZ),
281 odivCase(CNz,Cx,NewZ,X,Y,Yp), 282 intCase(Cyp,Yp),
283 odivCase(CNz,Cyp,NewZ,Yp,X,NewX), 284 285 (Y = Yp, \+(X = NewX) 286 -> intCase(CNx,NewX),
287 odivCase(CNz,CNx,NewZ,NewX,Y,NewY) 288 ; NewY = Yp
289 ).
290
298
299multCase(p,p, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)):-
300 NZl is maxr(Zl,roundtoward(Xl*Yl,to_negative)),
301 NZh is minr(Zh,roundtoward(Xh*Yh,to_positive)),
302 non_empty(NZl,NZh).
303multCase(p,n, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)):-
304 NZl is maxr(Zl,roundtoward(Xh*Yl,to_negative)),
305 NZh is minr(Zh,roundtoward(Xl*Yh,to_positive)),
306 non_empty(NZl,NZh).
307multCase(n,p, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)):-
308 NZl is maxr(Zl,roundtoward(Xl*Yh,to_negative)),
309 NZh is minr(Zh,roundtoward(Xh*Yl,to_positive)),
310 non_empty(NZl,NZh).
311multCase(n,n, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)):-
312 NZl is maxr(Zl,roundtoward(-Xh * -Yh,to_negative)), 313 NZh is minr(Zh,roundtoward(-Xl * -Yl,to_positive)),
314 non_empty(NZl,NZh).
315
316multCase(p,s, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)):-
317 NZl is maxr(Zl,roundtoward(Xh*Yl,to_negative)),
318 NZh is minr(Zh,roundtoward(Xh*Yh,to_positive)),
319 non_empty(NZl,NZh).
320multCase(n,s, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)):-
321 NZl is maxr(Zl,roundtoward(Xl*Yh,to_negative)),
322 NZh is minr(Zh,roundtoward(-Xl * -Yl,to_positive)), 323 non_empty(NZl,NZh).
324multCase(s,p, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)):-
325 NZl is maxr(Zl,roundtoward(Xl*Yh,to_negative)),
326 NZh is minr(Zh,roundtoward(Xh*Yh,to_positive)),
327 non_empty(NZl,NZh).
328multCase(s,n, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)):-
329 NZl is maxr(Zl,roundtoward(Xh*Yl,to_negative)),
330 NZh is minr(Zh,roundtoward(-Xl * -Yl,to_positive)), 331 non_empty(NZl,NZh).
332
333multCase(s,s, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)):-
334 NZl is maxr(Zl,minr(roundtoward(Xl*Yh,to_negative),roundtoward(Xh*Yl,to_negative))),
335 NZh is minr(Zh,maxr(roundtoward(-Xl* -Yl,to_positive),roundtoward(Xh*Yh,to_positive))),
336 non_empty(NZl,NZh).
337
344
345odivCase(p,p, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)) :-
346 sign(Yl)+sign(Xl) > 0 347 -> NZl is maxr(Zl,roundtoward(Xl/Yh,to_negative)),
348 NZh is minr(Zh,roundtoward(Xh/max(0.0,Yl),to_positive)),
349 non_empty(NZl,NZh)
350 ; NZl = Zl, NZh = Zh.
351odivCase(p,n, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)) :-
352 sign(Xl)-sign(Yh) > 0 353 -> NZl is maxr(Zl,roundtoward(-Xh / -min(-0.0,Yh),to_negative)), 354 NZh is minr(Zh,roundtoward(-Xl / -Yl ,to_positive)),
355 non_empty(NZl,NZh)
356 ; NZl = Zl, NZh = Zh.
357odivCase(p,s, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)) :-
358 Xl > 0 359 -> ZIh is roundtoward(Xl/Yh,to_negative), 360 ZIl is roundtoward(Xl/Yl,to_positive), 361 ( 1 is cmpr(Zl,ZIl) -> NZl is maxr(Zl,ZIh) ; NZl = Zl), 362 (-1 is cmpr(Zh,ZIh) -> NZh is minr(Zh,ZIl) ; NZh = Zh), 363 non_empty(NZl,NZh)
364 ; NZl = Zl, NZh = Zh.
365
366odivCase(n,p, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)):-
367 sign(Yl)-sign(Xh) > 0 368 -> NZl is maxr(Zl,roundtoward(-Xl / -max(0.0,Yl),to_negative)),
369 NZh is minr(Zh,roundtoward(-Xh / -Yh,to_positive)),
370 non_empty(NZl,NZh)
371 ; NZl = Zl, NZh = Zh.
372odivCase(n,n, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)) :-
373 sign(Yh)+sign(Xh) < 0 374 -> NZl is maxr(Zl,roundtoward(Xh/Yl,to_negative)),
375 NZh is minr(Zh,roundtoward(Xl/min(-0.0,Yh),to_positive)), 376 non_empty(NZl,NZh)
377 ; NZl = Zl, NZh = Zh.
378odivCase(n,s, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)) :-
379 Xh < 0 380 -> ZIh is roundtoward( Xh / Yl,to_negative), 381 ZIl is roundtoward(-Xh / -Yh,to_positive), 382 ( 1 is cmpr(Zl,ZIh) -> NZl is maxr(Zl,ZIl) ; NZl = Zl), 383 (-1 is cmpr(Zh,ZIl) -> NZh is minr(Zh,ZIh) ; NZh = Zh), 384 non_empty(NZl,NZh)
385 ; NZl = Zl, NZh = Zh.
386
387odivCase(s,p, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)) :-
388 Yl > 0 389 -> NZl is maxr(Zl,roundtoward(-Xl / -Yl,to_negative)),
390 NZh is minr(Zh,roundtoward( Xh / Yl,to_positive)),
391 non_empty(NZl,NZh)
392 ; NZl = Zl, NZh = Zh.
393odivCase(s,n, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)) :-
394 Yh < 0 395 -> NZl is maxr(Zl,roundtoward(Xh/Yh,to_negative)),
396 NZh is minr(Zh,roundtoward(Xl/Yh,to_positive)),
397 non_empty(NZl,NZh)
398 ; NZl = Zl, NZh = Zh.
399
400odivCase(s,s, _, _, Z, Z). 401
403
405
406narrowing_op(min, _, $((Zl,Zh),(Xl,Xh),(Yl,Yh)), New) :-
407 Z1l is maxr(Zl,minr(Xl,Yl)), 408 Z1h is minr(Zh,minr(Xh,Yh)),
409 ( Yh < Xl -> Case = 1 410 ; ( Xh < Yl -> Case = 2 411 ; Case = 3 412 )
413 ),
414 minimax(Case, (Zl,1.0Inf), (Z1l,Z1h),(Xl,Xh),(Yl,Yh), New).
415
416narrowing_op(max, _, $((Zl,Zh),(Xl,Xh),(Yl,Yh)), New) :-
417 Z1l is maxr(Zl,maxr(Xl,Yl)), 418 Z1h is minr(Zh,maxr(Xh,Yh)),
419 ( Xh < Yl -> Case = 1 420 ; ( Yh < Xl -> Case = 2 421 ; Case = 3 422 )
423 ),
424 minimax(Case, (-1.0Inf,Zh), (Z1l,Z1h),(Xl,Xh),(Yl,Yh), New).
425
426minimax(1, Zpartial, Z, X, Y, $(Z, X, NewY)) :- 427 ^(Y,Zpartial,NewY). 428minimax(2, Zpartial, Z, X, Y, $(Z, NewX, Y)) :- 429 ^(X,Zpartial,NewX). 430minimax(3, Zpartial, Z, X, Y, $(Z, NewX, NewY)) :- 431 ^(X,Zpartial,NewX), 432 ^(Y,Zpartial,NewY). 433
435
437
438narrowing_op(abs, _, $(Z,X), $(NewZ, NewX)) :-
439 intCase(Cx,X),
440 absCase(Cx,X,Z,NewZ), !,
441 non_empty(NewZ),
442 intersect_regions_(X,NewZ,NewX),
443 non_empty(NewX).
444
448intersect_regions_(Z,X,NewZ) :-
449 (^(Z,X,ZPos) -> true ; empty_interval(ZPos)), 450 (-(X,Z,ZNeg) -> true ; empty_interval(ZNeg)), 451 union_(ZPos,ZNeg,NewZ).
452
454-((Xl,Xh), (Zl,Zh), (NZl,NZh)) :-
455 NZl is maxr(Zl,-Xh), NZh is minr(Zh,-Xl),
456 non_empty(NZl,NZh).
457
461absCase(p, (Xl,Xh), (Zl,Zh), (NZl,NZh)) :- NZl is maxr(Zl,Xl), NZh is minr(Zh,Xh).
462absCase(n, (Xl,Xh), (Zl,Zh), (NZl,NZh)) :- NZl is maxr(Zl,-Xh), NZh is minr(Zh,-Xl).
463absCase(s, (Xl,Xh), (Zl,Zh), (NZl,NZh)) :- NZl is maxr(Zl,0), NZh is minr(Zh,maxr(-Xl,Xh)).
464
466
468
469narrowing_op(minus, _, $((Zl,Zh),(Xl,Xh)), $((NZl,NZh),(NXl,NXh))) :-
470 NZl is maxr(Zl,-Xh), NZh is minr(Zh,-Xl), 471 NXl is maxr(Xl,-NZh), NXh is minr(Xh,-NZl). 472
474
476
477narrowing_op(sqrt, _, $((Zl,Zh),(Xl,Xh)), $((NZl,NZh),(NXl,NXh))) :-
478 Xh>=0, Xl1 is maxr(0,Xl), 479 NZl is maxr(maxr(0,Zl),roundtoward(sqrt(Xl1),to_negative)),
480 NZh is minr(Zh,roundtoward(sqrt(Xh),to_positive)),
481 non_empty(NZl,NZh),
482 NXh is minr(Xh,roundtoward(NZh*NZh,to_positive)),
483 NXl is maxr(Xl,-NXh),
484 non_empty(NXl,NXh).
485
487
489
490narrowing_op(exp, _, $((Zl,Zh),(Xl,Xh)), $((NZl,NZh),(NXl,NXh))) :-
491 NZl is maxr(Zl,roundtoward(exp(Xl),to_negative)),
492 NZh is minr(Zh,roundtoward(exp(Xh),to_positive)),
493 non_empty(NZl,NZh),
494 NXl is maxr(Xl,roundtoward(log(NZl),to_negative)),
495 NXh is minr(Xh,roundtoward(log(NZh),to_positive)),
496 non_empty(NXl,NXh).
497
499
501
502narrowing_op(pow, P, $(Z,X,Y), $(NewZ,NewX,NewY)) :-
503 powr_(P, Z,X,Y, NewZ,NewX,NewY).
504
505powr_(p, Z,X,Y, NewZ,NewX,NewY) :- 506 Z = (Zl,Zh), X = (Xl,Xh), Y = (Yl,Yh),
507 ( equals_one(Xl), equals_one(Xh) -> 508 !, ^(Z,(1,1),NewZ), NewX = X, NewY = Y 509 ; equals_one(Zl), equals_one(Zh) -> 510 ((Xl =< -1,-1 =< Xh ; Xl =< 1,1 =< Xh)
511 -> fail 512 ; !, NewZ = Z, NewX = X, ^(Y,(0,0),NewY) 513 )
514 ; equals_zero(Yl), equals_zero(Yh) -> 515 !, ^(Z,(1,1),NewZ), NewX = X, NewY = Y 516 ).
517
518powr_(_, Z,X,(R,R), NewZ,NewX,NewY) :- rational(R,N,D), !, 519 N1 is N mod 2, D1 is D rem 2, 520 pt_powrCase(N1,D1,Z,X,R,NewZ), non_empty(NewZ),
521 Ri is 1/R,
522 pt_powrCase(D1,N1,X,NewZ,Ri,NewX), non_empty(NewX),
523 NewY = (R,R).
524
525powr_(_, Z,X,Y, NewZ,NewX,NewY) :- 526 intCase(C,X),
527 powr_intY_(C, Z,X,Y, NewZ,NewX,NewY).
528
529equals_zero(0).
530equals_zero(0.0).
531equals_zero(-0.0).
532
533equals_one(1).
534equals_one(1.0).
535
537pt_powrCase(N1,D1,Z,X,R,NewZ) :- R < 0, !, 538 Rp is -R,
539 universal_interval(UI),
540 intCase(Cz,Z),
541 odivCase(p,Cz,(1,1),Z,UI,Zi), 542 pt_powrCase(N1,D1,Zi,X,Rp,NZi),
543 intCase(CNzi,NZi),
544 odivCase(p,CNzi,(1,1),NZi,Z,NewZ). 545
546pt_powrCase(1,0,(Zl,Zh),(Xl,Xh),R,(NZl,NZh)) :- 547 548 Xh >=0, 549 Zl1 is maxr(0, roundtoward(Xl**R,to_negative)), 550 Zh1 is roundtoward(Xh**R,to_positive),
551 Zpl is maxr(Zl,Zl1), Zph is minr(Zh,Zh1), 552 Znl is maxr(Zl,-Zh1), Znh is minr(Zh,-Zl1), 553 ( 1 is cmpr(Znl,Znh) -> NZl is Zpl, NZh is Zph 554 ; 1 is cmpr(Zpl,Zph) -> NZl is Znl, NZh is Znh 555 ; NZl is minr(Zpl,Znl), NZh is maxr(Zph,Znh) 556 ).
557
558pt_powrCase(0,1,(Zl,Zh),(Xl,Xh),R,(NZl,NZh)) :- 559 560 ( Xh < 0 -> Xl1 is -Xh, Xh1 is -Xl 561 ;(Xl > 0 -> Xl1 = Xl, Xh1 = Xh 562 ; Xl1 = 0, Xh1 is maxr(-Xl,Xh) 563 )),
564 565 NZl is maxr(Zl, roundtoward(Xl1**R,to_negative)), 566 NZh is minr(Zh, roundtoward(Xh1**R,to_positive)).
567
568pt_powrCase(1,1,(Zl,Zh),(Xl,Xh),R,(NZl,NZh)) :- 569 570 NZl is maxr(Zl, roundtoward(Xl**R,to_negative)),
571 NZh is minr(Zh, roundtoward(Xh**R,to_positive)).
572
574powr_intY_(p, Z,X,Y, NewZ,NewX,NewY) :- 575 powr_prim_(Z,X,Y,NewZ), non_empty(NewZ), 576 universal_interval(UI),
577 intCase(Cy,Y),
578 odivCase(p,Cy,(1,1),Y,UI,Yi), 579 powr_prim_(X,NewZ,Yi,NewX), non_empty(NewX), 580 ln(NewZ,Ynum), intCase(Cyn,Ynum),
581 ln(NewX,Yden), intCase(Cyd,Yden),
582 odivCase(Cyn,Cyd,Ynum,Yden,Y,NewY). 583
584powr_intY_(n, Z,X,Y, NewZ,NewX,NewY) :- 585 586 587 588 589 X = (Xl,Xh), X1l is -Xh, X1h is -Xl,
590 powr_intY_(p, Z,(X1l,X1h),Y, (_,Z1h),(X2l,X2h),NewY), 591 X3l is -X2h, X3h is -X2l, ^(X,(X3l,X3h),NewX), 592 Z1l is -Z1h, ^(Z,(Z1l,Z1h),NewZ). 593
594powr_intY_(s, Z,X,Y, NewZ,NewX,NewY):- 595 596 X = (Xl,Xh), Xfl is -Xl, 597 powr_intY_(p, Z,(0,Xfl),Y, (_,Znh),NXn,NYn), 598 powr_intY_(p, Z,(0,Xh), Y, (_,Zph),NXp,NYp), 599 600 NZl is -Znh, NZh is maxr(Znh,Zph), ^(Z,(NZl,NZh),NewZ),
601 union_(NXn,NXp,NX), ^(X,NX,NewX),
602 union_(NYn,NYp,NY), ^(Y,NY,NewY).
603
604powr_prim_((Zl,Zh), (Xl,Xh), (Yl,Yh), (NZl,NZh)) :- 605 Zll is roundtoward(Xl**Yl,to_negative),
606 Zlh is roundtoward(Xl**Yh,to_positive),
607 Zhl is roundtoward(Xh**Yl,to_negative),
608 Zhh is roundtoward(Xh**Yh,to_positive),
609 NZl is maxr(Zl,minr(Zll, minr(Zlh, minr(Zhl, Zhh)))), 610 NZh is minr(Zh,maxr(Zll, maxr(Zlh, maxr(Zhl, Zhh)))). 611
612ln((Xl,Xh), (Zl,Zh)) :-
613 Zl is roundtoward(log(Xl),to_negative),
614 Zh is roundtoward(log(Xh),to_positive).
615
617
619
620narrowing_op(sin, _, $(Z,X), $(NewZ, NewX)) :-
621 sin_(X,S,Z1), ^(Z,Z1,NewZ),
622 asin_(NewZ,S,X1), ^(X,X1,NewX).
623
624sin_((X,X),S,NZ) :- 625 abs(X) =:= 1.0Inf, !,
626 S = (-10,10), NZ = (-1.0,1.0).
627sin_((Xl,Xh),(Sl,Sh),(NZl,NZh)) :-
628 Sl is round(Xl/pi), 629 Sh is round(Xh/pi), 630 to_float((Xl,Xh),(FXl,FXh)), 631 sinCase(FXl,FXh,Sl,Sh,NZl,NZh).
632
633sinCase(Xl,Xh,S,S,NZl,NZh) :- !, 634 ( 0 is S mod 2 635 -> NZl is roundtoward(sin(Xl), to_negative),
636 NZh is roundtoward(sin(Xh), to_positive)
637 ; NZl is roundtoward(sin(Xh), to_negative),
638 NZh is roundtoward(sin(Xl), to_positive)
639 ).
640sinCase(Xl,Xh,Sl,Sh,NZl,NZh) :- Sh-Sl =:= 1, !, 641 ( 0 is Sl mod 2 642 -> NZl is minr(roundtoward(sin(Xl), to_negative), 643 roundtoward(sin(Xh), to_negative)),
644 NZh = 1.0
645 ; NZl = -1.0, 646 NZh is maxr(roundtoward(sin(Xl), to_positive),
647 roundtoward(sin(Xh), to_positive))
648 ).
649sinCase(_,_,_,_,-1.0,1.0). 650
651asin_((Xl,Xh),(Sl,Sh),(NZl,NZh)) :-
652 asinCase(Xl,Xh,Sl,Sh,NZl,NZh).
653
656asinCase(Xl,Xh,S,S,NZl,NZh) :- !, 657 Zl is roundtoward(asin(Xl),to_negative),
658 Zh is roundtoward(asin(Xh),to_positive),
659 Offs is S*pi,
660 bounded_number(Offl,Offh,Offs), 661 ( 0 is S mod 2 662 -> NZl is roundtoward(Offl+Zl,to_negative),
663 NZh is roundtoward(Offh+Zh,to_positive)
664 ; NZl is roundtoward(Offl-Zh,to_negative),
665 NZh is roundtoward(Offh-Zl,to_positive)
666 ).
667asinCase(Xl,Xh,Sl,Sh,NZl,NZh) :- Sh-Sl =:= 1, !, 668 Offl is nexttoward(Sl*pi,-inf), Offh is nexttoward(Sh*pi,inf), 669 ( 0 is Sl mod 2 670 -> Z is roundtoward(asin(Xl),to_negative),
671 NZl is roundtoward(Offl+Z,to_negative),
672 NZh is roundtoward(Offh-Z,to_positive)
673 ; Z is roundtoward(asin(Xh),to_negative),
674 NZl is roundtoward(Offl-Z,to_negative),
675 NZh is roundtoward(Offh+Z,to_positive)
676 ).
677asinCase(_,_,_,_,-1.0Inf,1.0Inf). 678
680
682
683narrowing_op(cos, _, $(Z,X), $(NewZ, NewX)) :-
684 cos_(X,S,Z1), ^(Z,Z1,NewZ),
685 acos_(NewZ,S,X1), ^(X,X1,NewX).
686
687cos_((X,X),S,NZ) :- 688 abs(X) =:= 1.0Inf, !,
689 S = (-10,10), NZ = (-1.0,1.0).
690cos_((Xl,Xh),(Sl,Sh),(NZl,NZh)) :-
691 Sl is floor(Xl/pi),
692 Sh is floor(Xh/pi), 693 to_float((Xl,Xh),(FXl,FXh)), 694 cosCase(FXl,FXh,Sl,Sh,NZl,NZh).
695
696cosCase(Xl,Xh,S,S,NZl,NZh) :- !, 697 ( 1 is S mod 2 698 -> NZl is roundtoward(cos(Xl), to_negative),
699 NZh is roundtoward(cos(Xh), to_positive)
700 ; NZl is roundtoward(cos(Xh), to_negative),
701 NZh is roundtoward(cos(Xl), to_positive)
702 ).
703cosCase(Xl,Xh,Sl,Sh,NZl,NZh) :- Sh-Sl =:= 1, !, 704 ( 1 is Sl mod 2 705 -> NZl is minr(roundtoward(cos(Xl), to_negative), 706 roundtoward(cos(Xh), to_negative)),
707 NZh = 1.0
708 ; NZl = -1.0, 709 NZh is maxr(roundtoward(cos(Xl), to_positive),
710 roundtoward(cos(Xh), to_positive))
711 ).
712cosCase(_,_,_,_,-1.0,1.0). 713
714acos_((Xl,Xh),(Sl,Sh),(NZl,NZh)) :-
715 acosCase(Xl,Xh,Sl,Sh,NZl,NZh).
716
719acosCase(Xl,Xh,S,S,NZl,NZh) :- !, 720 Zl is roundtoward(acos(Xh),to_negative),
721 Zh is roundtoward(acos(Xl),to_positive),
722 ( 1 is S mod 2 723 -> Offs is (S+1)*pi,
724 bounded_number(Offl,Offh,Offs), 725 NZl is roundtoward(Offl-Zh,to_negative),
726 NZh is roundtoward(Offh-Zl,to_positive)
727 ; Offs is S*pi,
728 bounded_number(Offl,Offh,Offs), 729 NZl is roundtoward(Offl+Zl,to_negative),
730 NZh is roundtoward(Offh+Zh,to_positive)
731 ).
732acosCase(Xl,Xh,Sl,Sh,NZl,NZh) :- Sh-Sl =:= 1, !, 733 ( 1 is Sl mod 2 734 -> Z is roundtoward(acos(Xl),to_positive), 735 Offl is nexttoward((Sl+1)*pi,-inf), NZl is roundtoward(Offl-Z,to_negative),
736 Offh is nexttoward(Sh*pi,inf) , NZh is roundtoward(Offh+Z,to_positive)
737 ; Z is roundtoward(acos(Xh),to_positive), 738 Offl is nexttoward(Sl*pi,-inf) , NZl is roundtoward(Offl+Z,to_negative),
739 Offh is nexttoward((Sh+1)*pi,inf), NZh is roundtoward(Offh-Z,to_positive)
740 ).
741acosCase(_,_,_,_,-1.0Inf,1.0Inf). 742
744
746
747narrowing_op(tan, _, $(Z,X), $(NewZ, NewX)) :-
748 X = (Xl,Xh),
749 Sl is round(Xl/pi), 750 Sh is round(Xh/pi),
751 tanCase(Z,X,Sl,Sh,NewZ,NewX).
752
755tanCase(Z,(X,X),_,_,Z,(X,X)) :- 756 abs(X) =:= 1.0Inf, !.
757tanCase((Zl,Zh),(Xl,Xh),S,S,(NZl,NZh),(NXl,NXh)) :- !, 758 Z1l is roundtoward(tan(Xl),to_negative),
759 Z1h is roundtoward(tan(Xh),to_positive),
760 ( non_empty(Z1l,Z1h)
761 -> ^((Zl,Zh),(Z1l,Z1h),(NZl,NZh)), 762 Offs is S*pi,
763 bounded_number(Offl,Offh,Offs), 764 X1l is Offl+roundtoward(atan(NZl),to_negative),
765 X1h is Offh+roundtoward(atan(NZh),to_positive),
766 ^((Xl,Xh),(X1l,X1h),(NXl,NXh))
767 ; 768 ((Z1l =< Zh ; Z1h >= Zl) -> true), 769 NZl = Zl, NZh = Zh, NXl = Xl, NXh = Xh 770 ).
771tanCase(Z,(Xl,Xh),Sl,Sh,(NZl,NZh),(NXl,NXh)) :- Sh-Sl =:= 1, !, 772 Sg is pi*(Sl+1r2), 773 bounded_number(XSl,XSh,Sg), 774 (tanCase(Z,(Xl,XSl),Sl,Sl,(NZ1l,NZ1h),(NXl,NX1h))
775 -> (tanCase(Z,(XSh,Xh),Sh,Sh,_NZ2,(_,NXh))
776 -> NZl = -1.0Inf, NZh = 1.0Inf 777 ; NXh = NX1h, NZl = NZ1l, NZh = NZ1h 778 )
779 ; tanCase(Z,(XSh,Xh),Sh,Sh,(NZl,NZh),(NXl,NXh)) 780 ).
781tanCase(Z,X,Sl,Sh,Z,X). 782
784
786
787narrowing_op(not, P, $(Z,X), $(NewZ, NewX)) :-
788 notB_(X,Z,Z1,P), ^(Z,Z1,NewZ),
789 notB_(NewZ,X,X1,P), ^(X,X1,NewX).
790
791notB_((0,0), _, (1,1), p) :- !.
792notB_((1,1), _, (0,0), p) :- !.
793notB_( _, Z, NewZ, _) :- ^(Z,(0,1),NewZ).
794
796
798
799narrowing_op(and, P, $(Z,X,Y), $(NewZ, NewX, NewY)) :-
800 andB_rel_(Z,X,Y, NewZ, NewX, NewY, P).
801
802andB_rel_(Z,(1,1),(1,1), NewZ,(1,1),(1,1), p) :- !, ^(Z,(1,1),NewZ). 803andB_rel_(Z,(0,0),Y, NewZ,(0,0),Y, p) :- !, ^(Z,(0,0),NewZ), ^(Y,(0,1),NewY).
804andB_rel_(Z,X,(0,0), NewZ,X,(0,0), p) :- !, ^(Z,(0,0),NewZ), ^(X,(0,1),NewX).
805andB_rel_((1,1),X,Y, (1,1),NewX,NewY, p) :- !, ^(X,(1,1),NewX), ^(Y,(1,1),NewY).
806andB_rel_((0,0),X,(1,1), (0,0),NewX,(1,1), p) :- !, ^(X,(0,0),NewX).
807andB_rel_((0,0),(1,1),Y, (0,0),(1,1),NewY, p) :- !, ^(Y,(0,0),NewY).
808andB_rel_(Z,X,Y, NewZ,NewX,NewY, P) :- ^(Z,(0,1),NewZ), ^(X,(0,1),NewX), ^(Y,(0,1),NewY). 809
811
813
814narrowing_op(nand, P, $(Z,X,Y), $(NewZ, NewX, NewY)) :-
815 nandB_rel_(Z,X,Y, NewZ, NewX, NewY, P).
816
817nandB_rel_(Z,(1,1),(1,1), NewZ,(1,1),(1,1), p) :- !, ^(Z,(0,0),NewZ). 818nandB_rel_(Z,(0,0),Y, NewZ,(0,0),Y, p) :- !, ^(Z,(1,1),NewZ), ^(Y,(0,1),NewY).
819nandB_rel_(Z,X,(0,0), NewZ,X,(0,0), p) :- !, ^(Z,(1,1),NewZ), ^(X,(0,1),NewX).
820nandB_rel_((0,0),X,Y, (0,0),NewX,NewY, p) :- !, ^(X,(1,1),NewX), ^(Y,(1,1),NewY).
821nandB_rel_((1,1),X,(1,1), (1,1),NewX,(1,1), p) :- !, ^(X,(0,0),NewX).
822nandB_rel_((1,1),(1,1),Y, (1,1),(1,1),NewY, p) :- !, ^(Y,(0,0),NewY).
823nandB_rel_(Z,X,Y, NewZ,NewX,NewY, P) :- ^(Z,(0,1),NewZ), ^(X,(0,1),NewX), ^(Y,(0,1),NewY). 824
826
828
829narrowing_op(or, P, $(Z,X,Y), $(NewZ, NewX, NewY)) :-
830 orB_rel_(Z,X,Y, NewZ, NewX, NewY, P).
831
832orB_rel_(Z,(0,0),(0,0), NewZ,(0,0),(0,0), p) :- !, ^(Z,(0,0),NewZ). 833orB_rel_(Z,(1,1),Y, NewZ,(1,1),Y, p) :- !, ^(Z,(1,1),NewZ), ^(Y,(0,1),NewY).
834orB_rel_(Z,X,(1,1), NewZ,X,(1,1), p) :- !, ^(Z,(1,1),NewZ), ^(X,(0,1),NewX).
835orB_rel_((0,0),X,Y, (0,0),NewX,NewY, p) :- !, ^(X,(0,0),NewX), ^(Y,(0,0),NewY).
836orB_rel_((1,1),X,(0,0), (1,1),NewX,(0,0), p) :- !, ^(X,(1,1),NewX).
837orB_rel_((1,1),(0,0),Y, (1,1),(0,0),NewY, p) :- !, ^(Y,(1,1),NewY).
838orB_rel_(Z,X,Y, NewZ,NewX,NewY, P) :- ^(Z,(0,1),NewZ), ^(X,(0,1),NewX), ^(Y,(0,1),NewY). 839
841
843
844narrowing_op(nor, P, $(Z,X,Y), $(NewZ, NewX, NewY)) :-
845 norB_rel_(Z,X,Y, NewZ, NewX, NewY, P).
846
847norB_rel_(Z,(0,0),(0,0), NewZ,(0,0),(0,0), p) :- !, ^(Z,(1,1),NewZ). 848norB_rel_(Z,(1,1),Y, NewZ,(1,1),Y, p) :- !, ^(Z,(0,0),NewZ), ^(Y,(0,1),NewY).
849norB_rel_(Z,X,(1,1), NewZ,X,(1,1), p) :- !, ^(Z,(0,0),NewZ), ^(X,(0,1),NewX).
850norB_rel_((1,1),X,Y, (1,1),NewX,NewY, p) :- !, ^(X,(0,0),NewX), ^(Y,(0,0),NewY).
851norB_rel_((0,0),X,(0,0), (0,0),NewX,(0,0), p) :- !, ^(X,(1,1),NewX).
852norB_rel_((0,0),(0,0),Y, (0,0),(0,0),NewY, p) :- !, ^(Y,(1,1),NewY).
853norB_rel_(Z,X,Y, NewZ,NewX,NewY, P) :- ^(Z,(0,1),NewZ), ^(X,(0,1),NewX), ^(Y,(0,1),NewY). 854
856
858
859narrowing_op(xor, P, $(Z,X,Y), $(NewZ, NewX, NewY)) :-
860 xorB_rel_(Z,X,Y, NewZ, NewX, NewY, P).
861
862xorB_rel_(Z,(B,B),(B,B), NewZ,(B,B),(B,B), p) :- !, ^(Z,(0,0),NewZ). 863xorB_rel_(Z,(X,X),(Y,Y), NewZ,(X,X),(Y,Y), p) :- !, ^(Z,(1,1),NewZ).
864xorB_rel_((B,B),X,(B,B), (B,B),NewX,(B,B), p) :- !, ^(X,(0,0),NewX).
865xorB_rel_((Z,Z),X,(Y,Y), (Z,Z),NewX,(Y,Y), p) :- !, ^(X,(1,1),NewX).
866xorB_rel_((B,B),(B,B),Y, (B,B),(B,B),NewY, p) :- !, ^(Y,(0,0),NewY).
867xorB_rel_((Z,Z),(X,X),Y, (Z,Z),(X,X),NewY, p) :- !, ^(Y,(1,1),NewY).
868xorB_rel_(Z,X,Y, NewZ,NewX,NewY, P) :- ^(Z,(0,1),NewZ), ^(X,(0,1),NewX), ^(Y,(0,1),NewY). 869
871
873
874narrowing_op(imB, P, $(Z,X,Y), $(NewZ, NewX, NewY)) :-
875 imB_rel_(Z,X,Y, NewZ, NewX, NewY, P).
876
877imB_rel_(Z,(1,1),Y, New,(1,1),New, P) :- !, ^(Z,Y,New), (New=(B,B) -> P=p ; true). 878imB_rel_(Z,(0,0),Y, NewZ,(0,0),Y, p) :- !, ^(Z,(1,1),NewZ), ^(Y,(0,1),NewY). 879imB_rel_(Z,X,(0,0), NewZ,NewX,(0,0), P) :- !, notB_(X,Z,NewZ,P), notB_(NewZ,X,NewX,P). 880imB_rel_(Z,X,(1,1), NewZ,X,(1,1), P) :- !, ^(Z,(1,1),NewZ). 881imB_rel_((0,0),X,Y, (0,0),NewX,NewY, p) :- !, ^(X,(1,1),NewX), ^(Y,(0,0),NewY). 882imB_rel_(Z,X,Y, NewZ,NewX,NewY, P) :- ^(Z,(0,1),NewZ), ^(X,(0,1),NewX), ^(Y,(0,1),NewY). 883
885
886:- style_check(+singleton). 887
889
890narrowing_op(user(Prim), P, InArgs, OutArgs) :-
891 call_user_primitive(Prim, P, InArgs, OutArgs).