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 (non_empty(NXh,NYl) -> P=p ; true), 202 Out = $(Z,(Xl,NXh),(NYl,Yh))
203 ; ^(Z,(0,1),NewZ), 204 Out = $(NewZ,X,Y)
205 )
206 )
207 ).
208
210
212
213narrowing_op(lt, P, $(Z, X, Y), Out) :-
214 X = (Xl,Xh), Y = (Yl,Yh),
215 (-1 is cmpr(Xh,Yl) 216 -> ^(Z,(1,1),NewZ), P = p, Out = $(NewZ,X,Y) 217 ; (1 =\= cmpr(Yh,Xl) 218 -> ^(Z,(0,0),NewZ), P = p, Out = $(NewZ,X,Y) 219 ; (Z = (1,1) 220 -> next_lt_(Yh,-1,YhD), 221 NXh is minr(Xh,YhD), 222 next_lt_(Xl,1,XlU), 223 NYl is maxr(XlU,Yl), 224 (NXh < NYl -> P=p ; true), 225 Out = $(Z,(Xl,NXh),(NYl,Yh))
226 ; ^(Z,(0,1),NewZ), 227 Out = $(NewZ,X,Y)
228 )
229 )
230 ).
231
234next_lt_( 1.0Inf, _, 1.0Inf) :- !.
235next_lt_(-1.0Inf, _, -1.0Inf) :- !.
236next_lt_(V, -1, NV) :- NV is maxr(V-1,nexttoward(V,-1.0Inf)). 237next_lt_(V, 1, NV) :- NV is minr(V+1,nexttoward(V, 1.0Inf)).
238
240
242
243narrowing_op(in, P, $(Z, X, Y), $(NewZ, NewX, Y)):-
244 245 (^(X,Y,X1) 246 -> (Z == (1,1)
247 -> NewX = X1 248 ; NewX = X, 249 ^(Z,(0,1),NewZ) 250 )
251 ; ^(Z,(0,0),NewZ), 252 NewX = X, 253 P=p 254 ).
255
257
259
260narrowing_op(add, _, $((Zl,Zh), (Xl,Xh), (Yl,Yh)), $((NZl,NZh), (NXl,NXh), (NYl,NYh))) :-
261 NZl is maxr(Zl, roundtoward( Xl+Yl, to_negative)), 262 NZh is minr(Zh, roundtoward( Xh+Yh, to_positive)),
263 non_empty(NZl,NZh),
264 265 266 NXl is maxr(Xl, roundtoward(NZl+(-Yh), to_negative)), 267 NXh is minr(Xh, roundtoward(NZh+(-Yl), to_positive)),
268 non_empty(NXl,NXh),
269 NYl is maxr(Yl, roundtoward(NZl+(-NXh), to_negative)), 270 NYh is minr(Yh, roundtoward(NZh+(-NXl), to_positive)),
271 non_empty(NYl,NYh).
272
274
276
277narrowing_op(mul, _, $(Z,X,Y), $(NewZ, NewX, NewY)) :-
278 intCase(Cx,X),
279 intCase(Cy,Y),
280 multCase(Cx,Cy,X,Y,Z,NewZ), 281 intCase(CNz,NewZ),
282 odivCase(CNz,Cx,NewZ,X,Y,Yp), 283 intCase(Cyp,Yp),
284 odivCase(CNz,Cyp,NewZ,Yp,X,NewX), 285 286 (Y == Yp, X \== NewX 287 -> intCase(CNx,NewX),
288 odivCase(CNz,CNx,NewZ,NewX,Y,NewY) 289 ; NewY = Yp
290 ).
291
299
300multCase(p,p, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)):-
301 NZl is maxr(Zl,roundtoward(Xl*Yl,to_negative)),
302 NZh is minr(Zh,roundtoward(Xh*Yh,to_positive)),
303 non_empty(NZl,NZh).
304multCase(p,n, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)):-
305 NZl is maxr(Zl,roundtoward(Xh*Yl,to_negative)),
306 NZh is minr(Zh,roundtoward(Xl*Yh,to_positive)),
307 non_empty(NZl,NZh).
308multCase(n,p, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)):-
309 NZl is maxr(Zl,roundtoward(Xl*Yh,to_negative)),
310 NZh is minr(Zh,roundtoward(Xh*Yl,to_positive)),
311 non_empty(NZl,NZh).
312multCase(n,n, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)):-
313 NZl is maxr(Zl,roundtoward(-Xh * -Yh,to_negative)), 314 NZh is minr(Zh,roundtoward(-Xl * -Yl,to_positive)),
315 non_empty(NZl,NZh).
316
317multCase(p,s, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)):-
318 NZl is maxr(Zl,roundtoward(Xh*Yl,to_negative)),
319 NZh is minr(Zh,roundtoward(Xh*Yh,to_positive)),
320 non_empty(NZl,NZh).
321multCase(n,s, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)):-
322 NZl is maxr(Zl,roundtoward(Xl*Yh,to_negative)),
323 NZh is minr(Zh,roundtoward(-Xl * -Yl,to_positive)), 324 non_empty(NZl,NZh).
325multCase(s,p, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)):-
326 NZl is maxr(Zl,roundtoward(Xl*Yh,to_negative)),
327 NZh is minr(Zh,roundtoward(Xh*Yh,to_positive)),
328 non_empty(NZl,NZh).
329multCase(s,n, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)):-
330 NZl is maxr(Zl,roundtoward(Xh*Yl,to_negative)),
331 NZh is minr(Zh,roundtoward(-Xl * -Yl,to_positive)), 332 non_empty(NZl,NZh).
333
334multCase(s,s, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)):-
335 NZl is maxr(Zl,minr(roundtoward(Xl*Yh,to_negative),roundtoward(Xh*Yl,to_negative))),
336 NZh is minr(Zh,maxr(roundtoward(-Xl* -Yl,to_positive),roundtoward(Xh*Yh,to_positive))),
337 non_empty(NZl,NZh).
338
345
346odivCase(p,p, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)) :-
347 sign(Yl)+sign(Xl) > 0 348 -> NZl is maxr(Zl,roundtoward(Xl/Yh,to_negative)),
349 NZh is minr(Zh,roundtoward(Xh/max(0.0,Yl),to_positive)),
350 non_empty(NZl,NZh)
351 ; NZl = Zl, NZh = Zh.
352odivCase(p,n, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)) :-
353 sign(Xl)-sign(Yh) > 0 354 -> NZl is maxr(Zl,roundtoward(-Xh / -min(-0.0,Yh),to_negative)), 355 NZh is minr(Zh,roundtoward(-Xl / -Yl ,to_positive)),
356 non_empty(NZl,NZh)
357 ; NZl = Zl, NZh = Zh.
358odivCase(p,s, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)) :-
359 Xl > 0 360 -> ZIh is roundtoward(Xl/Yh,to_negative), 361 ZIl is roundtoward(Xl/Yl,to_positive), 362 ( 1 is cmpr(Zl,ZIl) -> NZl is maxr(Zl,ZIh) ; NZl = Zl), 363 (-1 is cmpr(Zh,ZIh) -> NZh is minr(Zh,ZIl) ; NZh = Zh), 364 non_empty(NZl,NZh)
365 ; NZl = Zl, NZh = Zh.
366
367odivCase(n,p, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)):-
368 sign(Yl)-sign(Xh) > 0 369 -> NZl is maxr(Zl,roundtoward(-Xl / -max(0.0,Yl),to_negative)), 370 NZh is minr(Zh,roundtoward(-Xh / -max(0.0,Yh),to_positive)),
371 non_empty(NZl,NZh)
372 ; NZl = Zl, NZh = Zh.
373odivCase(n,n, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)) :-
374 sign(Yh)+sign(Xh) < 0 375 -> NZl is maxr(Zl,roundtoward(Xh/Yl,to_negative)),
376 NZh is minr(Zh,roundtoward(Xl/min(-0.0,Yh),to_positive)), 377 non_empty(NZl,NZh)
378 ; NZl = Zl, NZh = Zh.
379odivCase(n,s, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)) :-
380 Xh < 0 381 -> ZIh is roundtoward( Xh / Yl,to_negative), 382 ZIl is roundtoward(-Xh / -Yh,to_positive), 383 ( 1 is cmpr(Zl,ZIh) -> NZl is maxr(Zl,ZIl) ; NZl = Zl), 384 (-1 is cmpr(Zh,ZIl) -> NZh is minr(Zh,ZIh) ; NZh = Zh), 385 non_empty(NZl,NZh)
386 ; NZl = Zl, NZh = Zh.
387
388odivCase(s,p, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)) :-
389 Yl > 0 390 -> NZl is maxr(Zl,roundtoward(-Xl / -Yl,to_negative)),
391 NZh is minr(Zh,roundtoward( Xh / Yl,to_positive)),
392 non_empty(NZl,NZh)
393 ; NZl = Zl, NZh = Zh.
394odivCase(s,n, (Xl,Xh), (Yl,Yh), (Zl,Zh), (NZl,NZh)) :-
395 Yh < 0 396 -> NZl is maxr(Zl,roundtoward(Xh/Yh,to_negative)),
397 NZh is minr(Zh,roundtoward(Xl/Yh,to_positive)),
398 non_empty(NZl,NZh)
399 ; NZl = Zl, NZh = Zh.
400
401odivCase(s,s, _, _, Z, Z). 402
404
406
407narrowing_op(min, _, $((Zl,Zh),(Xl,Xh),(Yl,Yh)), New) :-
408 Z1l is maxr(Zl,minr(Xl,Yl)), 409 Z1h is minr(Zh,minr(Xh,Yh)),
410 minimax((Zl,1.0Inf), (Z1l,Z1h),(Xl,Xh),(Yl,Yh), New).
411
412narrowing_op(max, _, $((Zl,Zh),(Xl,Xh),(Yl,Yh)), New) :-
413 Z1l is maxr(Zl,maxr(Xl,Yl)), 414 Z1h is minr(Zh,maxr(Xh,Yh)),
415 minimax((-1.0Inf,Zh), (Z1l,Z1h),(Xl,Xh),(Yl,Yh), New).
416
417minimax(_, Z1,X,Y, NewXYZ) :- 418 \+(^(Z1,X,_)), !, 419 NewXYZ = $(New, X, New),
420 ^(Y,Z1,New). 421
422minimax(_, Z1,X,Y, NewXYZ) :- 423 \+(^(Z1,Y,_)), !, 424 NewXYZ = $(New, New, Y),
425 ^(X,Z1,New). 426
427minimax(Zpartial, Z1,X,Y, $(Z1, NewX, NewY)) :- 428 ^(X,Zpartial,NewX), 429 ^(Y,Zpartial,NewY). 430
432
434
435narrowing_op(abs, _, $(Z,X), $(NewZ, NewX)) :-
436 intCase(Cx,X),
437 absCase(Cx,X,Z,NewZ), !,
438 non_empty(NewZ),
439 intersect_regions_(X,NewZ,NewX),
440 non_empty(NewX).
441
445intersect_regions_(Z,X,NewZ) :-
446 (^(Z,X,ZPos) -> true ; empty_interval(ZPos)), 447 (-(X,Z,ZNeg) -> true ; empty_interval(ZNeg)), 448 union_(ZPos,ZNeg,NewZ).
449
451-((Xl,Xh), (Zl,Zh), (NZl,NZh)) :-
452 NZl is maxr(Zl,-Xh), NZh is minr(Zh,-Xl),
453 non_empty(NZl,NZh).
454
458absCase(p, (Xl,Xh), (Zl,Zh), (NZl,NZh)) :- NZl is maxr(Zl,Xl), NZh is minr(Zh,Xh).
459absCase(n, (Xl,Xh), (Zl,Zh), (NZl,NZh)) :- NZl is maxr(Zl,-Xh), NZh is minr(Zh,-Xl).
460absCase(s, (Xl,Xh), (Zl,Zh), (NZl,NZh)) :- NZl is maxr(Zl,0), NZh is minr(Zh,maxr(-Xl,Xh)).
461
463
465
466narrowing_op(minus, _, $((Zl,Zh),(Xl,Xh)), $((NZl,NZh),(NXl,NXh))) :-
467 NZl is maxr(Zl,-Xh), NZh is minr(Zh,-Xl), 468 NXl is maxr(Xl,-NZh), NXh is minr(Xh,-NZl). 469
471
473
474narrowing_op(sqrt, _, $((Zl,Zh),(Xl,Xh)), $((NZl,NZh),(NXl,NXh))) :-
475 Xh>=0, Xl1 is maxr(0,Xl), 476 NZl is maxr(maxr(0,Zl),roundtoward(sqrt(Xl1),to_negative)),
477 NZh is minr(Zh,roundtoward(sqrt(Xh),to_positive)),
478 non_empty(NZl,NZh),
479 NXh is minr(Xh,roundtoward(NZh*NZh,to_positive)),
480 NXl is maxr(Xl,-NXh),
481 non_empty(NXl,NXh).
482
484
486
487narrowing_op(exp, _, $((Zl,Zh),(Xl,Xh)), $((NZl,NZh),(NXl,NXh))) :-
488 NZl is maxr(Zl,roundtoward(exp(Xl),to_negative)),
489 NZh is minr(Zh,roundtoward(exp(Xh),to_positive)),
490 non_empty(NZl,NZh),
491 NXl is maxr(Xl,roundtoward(log(NZl),to_negative)),
492 NXh is minr(Xh,roundtoward(log(NZh),to_positive)),
493 non_empty(NXl,NXh).
494
496
498
499narrowing_op(pow, P, $(Z,X,Y), $(NewZ,NewX,NewY)) :-
500 powr_(P, Z,X,Y, NewZ,NewX,NewY).
501
502powr_(p, Z,X,Y, NewZ,NewX,NewY) :- 503 Z = (Zl,Zh), X = (Xl,Xh), Y = (Yl,Yh),
504 ( equals_one(Xl), equals_one(Xh) -> 505 !, ^(Z,(1,1),NewZ), NewX = X, NewY = Y 506 ; equals_one(Zl), equals_one(Zh) -> 507 ((Xl =< -1,-1 =< Xh ; Xl =< 1,1 =< Xh)
508 -> fail 509 ; !, NewZ = Z, NewX = X, ^(Y,(0,0),NewY) 510 )
511 ; equals_zero(Yl), equals_zero(Yh) -> 512 !, ^(Z,(1,1),NewZ), NewX = X, NewY = Y 513 ).
514
515powr_(_, Z,X,(R,R), NewZ,NewX,NewY) :- rational(R,N,D), !, 516 N1 is N mod 2, D1 is D rem 2, 517 pt_powrCase(N1,D1,Z,X,R,NewZ), non_empty(NewZ),
518 Ri is 1/R,
519 pt_powrCase(D1,N1,X,NewZ,Ri,NewX), non_empty(NewX),
520 NewY = (R,R).
521
522powr_(_, Z,X,Y, NewZ,NewX,NewY) :- 523 intCase(C,X),
524 powr_intY_(C, Z,X,Y, NewZ,NewX,NewY).
525
526equals_zero(0).
527equals_zero(0.0).
528equals_zero(-0.0).
529
530equals_one(1).
531equals_one(1.0).
532
534pt_powrCase(N1,D1,Z,X,R,NewZ) :- R < 0, !, 535 Rp is -R,
536 universal_interval(UI),
537 intCase(Cz,Z),
538 odivCase(p,Cz,(1,1),Z,UI,Zi), 539 pt_powrCase(N1,D1,Zi,X,Rp,NZi),
540 intCase(CNzi,NZi),
541 odivCase(p,CNzi,(1,1),NZi,Z,NewZ). 542
543pt_powrCase(1,0,(Zl,Zh),(Xl,Xh),R,(NZl,NZh)) :- 544 545 Xh >=0, 546 Zl1 is maxr(0, roundtoward(Xl**R,to_negative)), 547 Zh1 is roundtoward(Xh**R,to_positive),
548 Zpl is maxr(Zl,Zl1), Zph is minr(Zh,Zh1), 549 Znl is maxr(Zl,-Zh1), Znh is minr(Zh,-Zl1), 550 ( 1 is cmpr(Znl,Znh) -> NZl is Zpl, NZh is Zph 551 ; 1 is cmpr(Zpl,Zph) -> NZl is Znl, NZh is Znh 552 ; NZl is minr(Zpl,Znl), NZh is maxr(Zph,Znh) 553 ).
554
555pt_powrCase(0,1,(Zl,Zh),(Xl,Xh),R,(NZl,NZh)) :- 556 557 ( Xh < 0 -> Xl1 is -Xh, Xh1 is -Xl 558 ;(Xl > 0 -> Xl1 = Xl, Xh1 = Xh 559 ; Xl1 = 0, Xh1 is maxr(-Xl,Xh) 560 )),
561 562 NZl is maxr(Zl, roundtoward(Xl1**R,to_negative)), 563 NZh is minr(Zh, roundtoward(Xh1**R,to_positive)).
564
565pt_powrCase(1,1,(Zl,Zh),(Xl,Xh),R,(NZl,NZh)) :- 566 567 NZl is maxr(Zl, roundtoward(Xl**R,to_negative)),
568 NZh is minr(Zh, roundtoward(Xh**R,to_positive)).
569
571powr_intY_(p, Z,X,Y, NewZ,NewX,NewY) :- 572 powr_prim_(Z,X,Y,NewZ), 573 universal_interval(UI),
574 intCase(Cy,Y),
575 odivCase(p,Cy,(1,1),Y,UI,Yi), 576 powr_prim_(X,NewZ,Yi,NewX), 577 ln(NewZ,Ynum), intCase(Cyn,Ynum),
578 ln(NewX,Yden), intCase(Cyd,Yden),
579 odivCase(Cyn,Cyd,Ynum,Yden,Y,NewY). 580
581powr_intY_(n, Z,X,Y, NewZ,NewX,NewY) :- 582 583 584 585 586 X = (Xl,Xh), X1l is -Xh, X1h is -Xl,
587 powr_intY_(p, Z,(X1l,X1h),Y, (_,Z1h),(X2l,X2h),NewY), 588 X3l is -X2h, X3h is -X2l, ^(X,(X3l,X3h),NewX), 589 Z1l is -Z1h, ^(Z,(Z1l,Z1h),NewZ). 590
591powr_intY_(s, Z,X,Y, NewZ,NewX,NewY):- 592 593 X = (Xl,Xh), Xfl is -Xl, 594 powr_intY_(p, Z,(0,Xfl),Y, (_,Znh),NXn,NYn), 595 powr_intY_(p, Z,(0,Xh), Y, (_,Zph),NXp,NYp), 596 597 NZl is -Znh, NZh is maxr(Znh,Zph), ^(Z,(NZl,NZh),NewZ),
598 union_(NXn,NXp,NX), ^(X,NX,NewX),
599 union_(NYn,NYp,NY), ^(Y,NY,NewY).
600
601powr_prim_((Zl,Zh), (Xl,Xh), (Yl,Yh), NewZ) :- 602 Zll is roundtoward(Xl**Yl,to_negative),
603 Zlh is roundtoward(Xl**Yh,to_positive),
604 Zhl is roundtoward(Xh**Yl,to_negative),
605 Zhh is roundtoward(Xh**Yh,to_positive),
606 NZl is maxr(Zl,minr(Zll, minr(Zlh, minr(Zhl, Zhh)))), 607 NZh is minr(Zh,maxr(Zll, maxr(Zlh, maxr(Zhl, Zhh)))), 608 (non_empty(NZl,NZh) -> NewZ = (NZl,NZh) ; NewZ = (Zl,Zh)). 609
610ln((Xl,Xh), (Zl,Zh)) :-
611 Zl is roundtoward(log(Xl),to_negative),
612 Zh is roundtoward(log(Xh),to_positive).
613
615
617
618narrowing_op(sin, _, $(Z,X), $(NewZ, NewX)) :-
619 sin_(X,S,Z1), ^(Z,Z1,NewZ),
620 asin_(NewZ,S,X1), ^(X,X1,NewX).
621
622sin_((X,X),S,NZ) :- 623 abs(X) =:= 1.0Inf, !,
624 S = (-10,10), NZ = (-1.0,1.0).
625sin_((Xl,Xh),(Sl,Sh),(NZl,NZh)) :-
626 Sl is round(Xl/pi), 627 Sh is round(Xh/pi), 628 to_float((Xl,Xh),(FXl,FXh)), 629 sinCase(FXl,FXh,Sl,Sh,NZl,NZh).
630
631sinCase(Xl,Xh,S,S,NZl,NZh) :- !, 632 ( 0 is S mod 2 633 -> NZl is roundtoward(sin(Xl), to_negative),
634 NZh is roundtoward(sin(Xh), to_positive)
635 ; NZl is roundtoward(sin(Xh), to_negative),
636 NZh is roundtoward(sin(Xl), to_positive)
637 ).
638sinCase(Xl,Xh,Sl,Sh,NZl,NZh) :- Sh-Sl =:= 1, !, 639 ( 0 is Sl mod 2 640 -> NZl is minr(roundtoward(sin(Xl), to_negative), 641 roundtoward(sin(Xh), to_negative)),
642 NZh = 1.0
643 ; NZl = -1.0, 644 NZh is maxr(roundtoward(sin(Xl), to_positive),
645 roundtoward(sin(Xh), to_positive))
646 ).
647sinCase(_,_,_,_,-1.0,1.0). 648
649asin_((Xl,Xh),(Sl,Sh),(NZl,NZh)) :-
650 asinCase(Xl,Xh,Sl,Sh,NZl,NZh).
651
654asinCase(Xl,Xh,S,S,NZl,NZh) :- !, 655 Zl is roundtoward(asin(Xl),to_negative),
656 Zh is roundtoward(asin(Xh),to_positive),
657 Offs is S*pi,
658 bounded_number(Offl,Offh,Offs), 659 ( 0 is S mod 2 660 -> NZl is roundtoward(Offl+Zl,to_negative),
661 NZh is roundtoward(Offh+Zh,to_positive)
662 ; NZl is roundtoward(Offl-Zh,to_negative),
663 NZh is roundtoward(Offh-Zl,to_positive)
664 ).
665asinCase(Xl,Xh,Sl,Sh,NZl,NZh) :- Sh-Sl =:= 1, !, 666 Offl is nexttoward(Sl*pi,-inf), Offh is nexttoward(Sh*pi,inf), 667 ( 0 is Sl mod 2 668 -> Z is roundtoward(asin(Xl),to_negative),
669 NZl is roundtoward(Offl+Z,to_negative),
670 NZh is roundtoward(Offh-Z,to_positive)
671 ; Z is roundtoward(asin(Xh),to_negative),
672 NZl is roundtoward(Offl-Z,to_negative),
673 NZh is roundtoward(Offh+Z,to_positive)
674 ).
675asinCase(_,_,_,_,-1.0Inf,1.0Inf). 676
678
680
681narrowing_op(cos, _, $(Z,X), $(NewZ, NewX)) :-
682 cos_(X,S,Z1), ^(Z,Z1,NewZ),
683 acos_(NewZ,S,X1), ^(X,X1,NewX).
684
685cos_((X,X),S,NZ) :- 686 abs(X) =:= 1.0Inf, !,
687 S = (-10,10), NZ = (-1.0,1.0).
688cos_((Xl,Xh),(Sl,Sh),(NZl,NZh)) :-
689 Sl is floor(Xl/pi),
690 Sh is floor(Xh/pi), 691 to_float((Xl,Xh),(FXl,FXh)), 692 cosCase(FXl,FXh,Sl,Sh,NZl,NZh).
693
694cosCase(Xl,Xh,S,S,NZl,NZh) :- !, 695 ( 1 is S mod 2 696 -> NZl is roundtoward(cos(Xl), to_negative),
697 NZh is roundtoward(cos(Xh), to_positive)
698 ; NZl is roundtoward(cos(Xh), to_negative),
699 NZh is roundtoward(cos(Xl), to_positive)
700 ).
701cosCase(Xl,Xh,Sl,Sh,NZl,NZh) :- Sh-Sl =:= 1, !, 702 ( 1 is Sl mod 2 703 -> NZl is minr(roundtoward(cos(Xl), to_negative), 704 roundtoward(cos(Xh), to_negative)),
705 NZh = 1.0
706 ; NZl = -1.0, 707 NZh is maxr(roundtoward(cos(Xl), to_positive),
708 roundtoward(cos(Xh), to_positive))
709 ).
710cosCase(_,_,_,_,-1.0,1.0). 711
712acos_((Xl,Xh),(Sl,Sh),(NZl,NZh)) :-
713 acosCase(Xl,Xh,Sl,Sh,NZl,NZh).
714
717acosCase(Xl,Xh,S,S,NZl,NZh) :- !, 718 Zl is roundtoward(acos(Xh),to_negative),
719 Zh is roundtoward(acos(Xl),to_positive),
720 ( 1 is S mod 2 721 -> Offs is (S+1)*pi,
722 bounded_number(Offl,Offh,Offs), 723 NZl is roundtoward(Offl-Zh,to_negative),
724 NZh is roundtoward(Offh-Zl,to_positive)
725 ; Offs is S*pi,
726 bounded_number(Offl,Offh,Offs), 727 NZl is roundtoward(Offl+Zl,to_negative),
728 NZh is roundtoward(Offh+Zh,to_positive)
729 ).
730acosCase(Xl,Xh,Sl,Sh,NZl,NZh) :- Sh-Sl =:= 1, !, 731 ( 1 is Sl mod 2 732 -> Z is roundtoward(acos(Xl),to_positive), 733 Offl is nexttoward((Sl+1)*pi,-inf), NZl is roundtoward(Offl-Z,to_negative),
734 Offh is nexttoward(Sh*pi,inf) , NZh is roundtoward(Offh+Z,to_positive)
735 ; Z is roundtoward(acos(Xh),to_positive), 736 Offl is nexttoward(Sl*pi,-inf) , NZl is roundtoward(Offl+Z,to_negative),
737 Offh is nexttoward((Sh+1)*pi,inf), NZh is roundtoward(Offh-Z,to_positive)
738 ).
739acosCase(_,_,_,_,-1.0Inf,1.0Inf). 740
742
744
745narrowing_op(tan, _, $(Z,X), $(NewZ, NewX)) :-
746 X = (Xl,Xh),
747 Sl is round(Xl/pi), 748 Sh is round(Xh/pi),
749 tanCase(Z,X,Sl,Sh,NewZ,NewX).
750
753tanCase(Z,(X,X),_,_,Z,(X,X)) :- 754 abs(X) =:= 1.0Inf, !.
755tanCase((Zl,Zh),(Xl,Xh),S,S,(NZl,NZh),(NXl,NXh)) :- !, 756 Z1l is roundtoward(tan(Xl),to_negative),
757 Z1h is roundtoward(tan(Xh),to_positive),
758 ( non_empty(Z1l,Z1h)
759 -> ^((Zl,Zh),(Z1l,Z1h),(NZl,NZh)), 760 Offs is S*pi,
761 bounded_number(Offl,Offh,Offs), 762 X1l is Offl+roundtoward(atan(NZl),to_negative),
763 X1h is Offh+roundtoward(atan(NZh),to_positive),
764 ^((Xl,Xh),(X1l,X1h),(NXl,NXh))
765 ; 766 ((Z1l =< Zh ; Z1h >= Zl) -> true), 767 NZl = Zl, NZh = Zh, NXl = Xl, NXh = Xh 768 ).
769tanCase(Z,(Xl,Xh),Sl,Sh,(NZl,NZh),(NXl,NXh)) :- Sh-Sl =:= 1, !, 770 Sg is pi*(Sl+1r2), 771 bounded_number(XSl,XSh,Sg), 772 (tanCase(Z,(Xl,XSl),Sl,Sl,(NZ1l,NZ1h),(NXl,NX1h))
773 -> (tanCase(Z,(XSh,Xh),Sh,Sh,_NZ2,(_,NXh))
774 -> NZl = -1.0Inf, NZh = 1.0Inf 775 ; NXh = NX1h, NZl = NZ1l, NZh = NZ1h 776 )
777 ; tanCase(Z,(XSh,Xh),Sh,Sh,(NZl,NZh),(NXl,NXh)) 778 ).
779tanCase(Z,X,Sl,Sh,Z,X). 780
782
784
785narrowing_op(not, P, $(Z,X), $(NewZ, NewX)) :-
786 notB_(X,Z,Z1,P), ^(Z,Z1,NewZ),
787 notB_(NewZ,X,X1,P), ^(X,X1,NewX).
788
789notB_((0,0), _, (1,1), p) :- !.
790notB_((1,1), _, (0,0), p) :- !.
791notB_( _, Z, NewZ, _) :- ^(Z,(0,1),NewZ).
792
794
796
797narrowing_op(and, P, $(Z,X,Y), $(NewZ, NewX, NewY)) :-
798 andB_rel_(Z,X,Y, NewZ, NewX, NewY, P).
799
800andB_rel_(Z,(1,1),(1,1), NewZ,(1,1),(1,1), p) :- !, ^(Z,(1,1),NewZ). 801andB_rel_(Z,(0,0),Y, NewZ,(0,0),Y, p) :- !, ^(Z,(0,0),NewZ), ^(Y,(0,1),NewY).
802andB_rel_(Z,X,(0,0), NewZ,X,(0,0), p) :- !, ^(Z,(0,0),NewZ), ^(X,(0,1),NewX).
803andB_rel_((1,1),X,Y, (1,1),NewX,NewY, p) :- !, ^(X,(1,1),NewX), ^(Y,(1,1),NewY).
804andB_rel_((0,0),X,(1,1), (0,0),NewX,(1,1), p) :- !, ^(X,(0,0),NewX).
805andB_rel_((0,0),(1,1),Y, (0,0),(1,1),NewY, p) :- !, ^(Y,(0,0),NewY).
806andB_rel_(Z,X,Y, NewZ,NewX,NewY, P) :- ^(Z,(0,1),NewZ), ^(X,(0,1),NewX), ^(Y,(0,1),NewY). 807
809
811
812narrowing_op(nand, P, $(Z,X,Y), $(NewZ, NewX, NewY)) :-
813 nandB_rel_(Z,X,Y, NewZ, NewX, NewY, P).
814
815nandB_rel_(Z,(1,1),(1,1), NewZ,(1,1),(1,1), p) :- !, ^(Z,(0,0),NewZ). 816nandB_rel_(Z,(0,0),Y, NewZ,(0,0),Y, p) :- !, ^(Z,(1,1),NewZ), ^(Y,(0,1),NewY).
817nandB_rel_(Z,X,(0,0), NewZ,X,(0,0), p) :- !, ^(Z,(1,1),NewZ), ^(X,(0,1),NewX).
818nandB_rel_((0,0),X,Y, (0,0),NewX,NewY, p) :- !, ^(X,(1,1),NewX), ^(Y,(1,1),NewY).
819nandB_rel_((1,1),X,(1,1), (1,1),NewX,(1,1), p) :- !, ^(X,(0,0),NewX).
820nandB_rel_((1,1),(1,1),Y, (1,1),(1,1),NewY, p) :- !, ^(Y,(0,0),NewY).
821nandB_rel_(Z,X,Y, NewZ,NewX,NewY, P) :- ^(Z,(0,1),NewZ), ^(X,(0,1),NewX), ^(Y,(0,1),NewY). 822
824
826
827narrowing_op(or, P, $(Z,X,Y), $(NewZ, NewX, NewY)) :-
828 orB_rel_(Z,X,Y, NewZ, NewX, NewY, P).
829
830orB_rel_(Z,(0,0),(0,0), NewZ,(0,0),(0,0), p) :- !, ^(Z,(0,0),NewZ). 831orB_rel_(Z,(1,1),Y, NewZ,(1,1),Y, p) :- !, ^(Z,(1,1),NewZ), ^(Y,(0,1),NewY).
832orB_rel_(Z,X,(1,1), NewZ,X,(1,1), p) :- !, ^(Z,(1,1),NewZ), ^(X,(0,1),NewX).
833orB_rel_((0,0),X,Y, (0,0),NewX,NewY, p) :- !, ^(X,(0,0),NewX), ^(Y,(0,0),NewY).
834orB_rel_((1,1),X,(0,0), (1,1),NewX,(0,0), p) :- !, ^(X,(1,1),NewX).
835orB_rel_((1,1),(0,0),Y, (1,1),(0,0),NewY, p) :- !, ^(Y,(1,1),NewY).
836orB_rel_(Z,X,Y, NewZ,NewX,NewY, P) :- ^(Z,(0,1),NewZ), ^(X,(0,1),NewX), ^(Y,(0,1),NewY). 837
839
841
842narrowing_op(nor, P, $(Z,X,Y), $(NewZ, NewX, NewY)) :-
843 norB_rel_(Z,X,Y, NewZ, NewX, NewY, P).
844
845norB_rel_(Z,(0,0),(0,0), NewZ,(0,0),(0,0), p) :- !, ^(Z,(1,1),NewZ). 846norB_rel_(Z,(1,1),Y, NewZ,(1,1),Y, p) :- !, ^(Z,(0,0),NewZ), ^(Y,(0,1),NewY).
847norB_rel_(Z,X,(1,1), NewZ,X,(1,1), p) :- !, ^(Z,(0,0),NewZ), ^(X,(0,1),NewX).
848norB_rel_((1,1),X,Y, (1,1),NewX,NewY, p) :- !, ^(X,(0,0),NewX), ^(Y,(0,0),NewY).
849norB_rel_((0,0),X,(0,0), (0,0),NewX,(0,0), p) :- !, ^(X,(1,1),NewX).
850norB_rel_((0,0),(0,0),Y, (0,0),(0,0),NewY, p) :- !, ^(Y,(1,1),NewY).
851norB_rel_(Z,X,Y, NewZ,NewX,NewY, P) :- ^(Z,(0,1),NewZ), ^(X,(0,1),NewX), ^(Y,(0,1),NewY). 852
854
856
857narrowing_op(xor, P, $(Z,X,Y), $(NewZ, NewX, NewY)) :-
858 xorB_rel_(Z,X,Y, NewZ, NewX, NewY, P).
859
860xorB_rel_(Z,(B,B),(B,B), NewZ,(B,B),(B,B), p) :- !, ^(Z,(0,0),NewZ). 861xorB_rel_(Z,(X,X),(Y,Y), NewZ,(X,X),(Y,Y), p) :- !, ^(Z,(1,1),NewZ).
862xorB_rel_((B,B),X,(B,B), (B,B),NewX,(B,B), p) :- !, ^(X,(0,0),NewX).
863xorB_rel_((Z,Z),X,(Y,Y), (Z,Z),NewX,(Y,Y), p) :- !, ^(X,(1,1),NewX).
864xorB_rel_((B,B),(B,B),Y, (B,B),(B,B),NewY, p) :- !, ^(Y,(0,0),NewY).
865xorB_rel_((Z,Z),(X,X),Y, (Z,Z),(X,X),NewY, p) :- !, ^(Y,(1,1),NewY).
866xorB_rel_(Z,X,Y, NewZ,NewX,NewY, P) :- ^(Z,(0,1),NewZ), ^(X,(0,1),NewX), ^(Y,(0,1),NewY). 867
869
871
872narrowing_op(imB, P, $(Z,X,Y), $(NewZ, NewX, NewY)) :-
873 imB_rel_(Z,X,Y, NewZ, NewX, NewY, P).
874
875imB_rel_(Z,(1,1),Y, New,(1,1),New, P) :- !, ^(Z,Y,New), (New=(B,B) -> P=p ; true). 876imB_rel_(Z,(0,0),Y, NewZ,(0,0),Y, p) :- !, ^(Z,(1,1),NewZ), ^(Y,(0,1),NewY). 877imB_rel_(Z,X,(0,0), NewZ,NewX,(0,0), P) :- !, notB_(X,Z,NewZ,P), notB_(NewZ,X,NewX,P). 878imB_rel_(Z,X,(1,1), NewZ,X,(1,1), P) :- !, ^(Z,(1,1),NewZ). 879imB_rel_((0,0),X,Y, (0,0),NewX,NewY, p) :- !, ^(X,(1,1),NewX), ^(Y,(0,0),NewY). 880imB_rel_(Z,X,Y, NewZ,NewX,NewY, P) :- ^(Z,(0,1),NewZ), ^(X,(0,1),NewX), ^(Y,(0,1),NewY). 881
883
884:- style_check(+singleton). 885
887
888narrowing_op(user(Prim), P, InArgs, OutArgs) :-
889 call_user_primitive(Prim, P, InArgs, OutArgs).