1%
    2% Hook for custom interval functions
    3%
    4% 1. Declare function
    5%    with int_hook(Name, pred_name(Args), Ret, [Options]).
    6%    'Args' specify the type of the arguments, '...' for intervals and 'atomic'
    7%    for numbers and logical values. Same for 'Ret'.
    8%    Add 'evaluate(false)' to 'Options' list to skip evaluation of arguments.
    9% 2. Calculate result with interval:pred_name(Args, Res)
   10%    'Args' are the argument types as defined in the hook with variable names,
   11%    e.g., 'L...U', 'atomic(A)'.
   12%
   13% See example below for (/)/2.
   14%
   15
   16%
   17% Monotonically behaving functions
   18%
   19% +: increasing
   20% -: decreasing
   21% *: increasing or decreasing
   22% **: commutative, all *
   23% /: symbolic argument or flag
   24mono((+)/1, [+]).
   25mono((+)/2, [+, +]).
   26mono((-)/1, [-]).
   27mono((-)/2, [+, -]).
   28mono((*)/2, **).
   29mono((**)/2, [*, /]). % power
   30mono(exp/1, [+]).
   31mono(min/2, [+, +]).
   32mono(max/2, [+, +]).
   33
   34%
   35% Comparison
   36%
   37int_hook(<, less1(atomic, atomic), _, []).
   38
   39less1(atomic(A), atomic(B), Res, _Flags) :-
   40    A < B,
   41    !,
   42    Res = true.
   43
   44less1(atomic(_) < atomic(_), Res, _Flags) :-
   45    !,
   46    Res = false.
   47
   48int_hook(<, less2(..., ...), _, []).
   49
   50less2(_...A2, B1..._, Res, _Flags) :-
   51    A2 < B1,
   52    !,
   53    Res = true.
   54
   55less2(_..._, _..._, false2, _Flags).
   56
   57int_hook(=<, less_eq(_, _), _, []).
   58less_eq(A, B, Res, Flags) :-
   59    interval_(A, A1..._, Flags),
   60    interval_(B, _...B2, Flags),
   61    A1 =< B2,
   62    !,
   63    Res = true.
   64
   65less_eq(_, _, false, _).
   66
   67int_hook(>, great(..., ...), _, []).
   68great(A1..._, _...B2, Res, _Flags) :-
   69    A1 > B2,
   70    !,
   71    Res = true.
   72
   73great(_..._, _..._, false, _Flags).
   74
   75int_hook(>=, great_eq(_, _), _, []).
   76great_eq(A, B, Res, Flags) :-
   77    interval_(A, _...A2, Flags),
   78    interval_(B, B1..._, Flags),
   79    A2 >= B1,
   80    !,
   81    Res = true.
   82
   83great_eq(_, _, false, _).
   84
   85int_hook(=\=, not_eq(..., ...), _, []).
   86not_eq(A...B, C...D, Res, Flags) :-
   87    (   less2(A...B, C...D, true, Flags)
   88    ;   great(A...B, C...D, true, Flags)
   89    ), !,
   90    Res = true.
   91
   92not_eq(_..._, _..._, false, _Flags).
   93
   94int_hook(=:=, eq0(_, _), _, []).
   95eq0(A, B, Res, Flags) :-
   96    less_eq(A, B, true, Flags),
   97    great_eq(A, B, true, Flags),
   98    !,
   99    Res = true.
  100
  101eq0(_, _, false, _Flags). 
  102
  103%
  104% Division
  105%
  106int_hook(/, div1(atomic, atomic), atomic, []).
  107div1(atomic(A), atomic(B), atomic(Res), _Flags) :-
  108    Res is A / B.
  109
  110int_hook(/, div2(..., ...), ..., []).
  111div2(A...B, C...D, Res, Flags) :-
  112    !,
  113    div(A...B, C...D, Res, Flags).
  114
  115int_hook(/, div3(..., atomic), ..., []).
  116div3(L...U, atomic(A), Res, Flags) :-
  117    !,
  118    div(L...U, A...A, Res, Flags).
  119
  120int_hook(/, div4(atomic, ...), ..., []).
  121div4(atomic(A), L...U, Res, Flags) :-
  122    !,
  123    div(A...A, L...U, Res, Flags).
  124
  125% Hickey Figure 1
  126mixed(L, U) :-
  127    L < 0,
  128    U > 0.
  129
  130positive(L, U) :-
  131    L >= 0,
  132    U > 0.
  133
  134zeropos(L, U) :-
  135    L =:= 0,
  136    U > 0.
  137
  138strictpos(L, _) :-
  139    L > 0.
  140
  141negative(L, U) :-
  142    L < 0,
  143    U =< 0.
  144
  145zeroneg(L, U) :-
  146    L < 0,
  147    U =:= 0.
  148
  149strictneg(_, U) :-
  150    U < 0.
  151
  152zero(L, U) :-
  153    L =:= 0,
  154    U =:= 0.
  155
  156%
  157% atomic 0 in numerator or denominator
  158%
  159div(A...B, C...D, Res, _Flags),
  160    zero(A, B),
  161    (   negative(C, D)
  162    ;   mixed(C, D) 
  163    ;   positive(C, D)
  164    )
  165 => Res = atomic(0).
  166
  167div(A...B, C...D, Res, _Flags),
  168    zero(C, D),
  169    zeropos(A, B)
  170 => Res = atomic(1.0Inf).
  171
  172div(A...B, C...D, Res, _Flags),
  173    zero(C, D),
  174    zeroneg(A, B)
  175 => Res = atomic(-1.0Inf).
  176
  177div(A...B, C...D, Res, _Flags),
  178    zero(C, D),
  179    mixed(A, B)
  180 => (   Res = atomic(-1.0Inf)
  181    ;   Res = atomic(1.0Inf)
  182    ).
  183
  184div(A...B, C...D, Res, _Flags),
  185    zero(A, B),
  186    zero(C, D)
  187 => Res = atomic(1.5NaN).
  188
  189%
  190% Hickey Theorem 8 and Figure 4
  191%
  192% P1 / P (special case, then general case)
  193div(A...B, C...D, Res, _Flags),
  194    strictpos(A, B),
  195    zeropos(C, D)
  196 => eval(A / D, L),
  197    Res = L...1.0Inf.
  198
  199div(A...B, C...D, Res, _Flags),
  200    strictpos(A, B),
  201    positive(C, D)
  202 => eval(A / D, L),
  203    eval(B / C, U),
  204    Res = L...U.
  205
  206% P0 / P
  207div(A...B, 0.0...D, Res, _Flags),
  208    zeropos(A, B),
  209    positive(0.0, D)
  210 => Res = 0.0...1.0Inf.
  211
  212div(A...B, C...D, Res, _Flags),
  213    zeropos(A, B),
  214    positive(C, D)
  215 => eval(B / C, U),
  216    Res = 0.0...U.
  217
  218% M / P
  219div(A...B, 0.0...D, Res, _Flags),
  220    mixed(A, B),
  221    positive(0.0, D)
  222 => Res = -1.0Inf...1.0Inf.
  223
  224div(A...B, C...D, Res, _Flags),
  225    mixed(A, B),
  226    positive(C, D)
  227 => eval(A / C, L),
  228    eval(B / C, U),
  229    Res = L...U.
  230
  231% N0 / P
  232div(A...B, 0.0...D, Res, _Flags),
  233    zeroneg(A, B),
  234    positive(0.0, D)
  235 => Res = -1.0Inf...0.0.
  236
  237div(A...B, C...D, Res, _Flags),
  238    zeroneg(A, B),
  239    positive(C, D)
  240 => eval(A / C, L),
  241    Res = L...0.0.
  242
  243% N1 / P
  244div(A...B, 0.0...D, Res, _Flags),
  245    strictneg(A, B),
  246    positive(0.0, D)
  247 => eval(B / D, U),
  248    Res = -1.0Inf...U.
  249
  250div(A...B, C...D, Res, _Flags),
  251    strictneg(A, B),
  252    positive(C, D)
  253 => eval(A / C, L),
  254    eval(B / D, U),
  255    Res = L...U.
  256
  257% P1 / M (2 solutions)
  258div(A...B, C...D, Res, _Flags),
  259    strictpos(A, B),
  260    mixed(C, D)
  261 => (   eval(A / C, U),
  262        Res = -1.0Inf...U
  263    ;   eval(A / D, L),
  264        Res = L...1.0Inf
  265    ).
  266
  267% P0 / M
  268div(A...B, C...D, Res, _Flags),
  269    zeropos(A, B),
  270    mixed(C, D)
  271 => Res = -1.0Inf...1.0Inf.
  272
  273% M / M
  274div(A...B, C...D, Res, _Flags),
  275    mixed(A, B),
  276    mixed(C, D)
  277 => Res = -1.0Inf...1.0Inf.
  278
  279% N0 / M
  280div(A...B, C...D, Res, _Flags),
  281    zeroneg(A, B),
  282    mixed(C, D)
  283 => Res = -1.0Inf...1.0Inf.
  284
  285% N1 / M (2 solutions)
  286div(A...B, C...D, Res, _Flags),
  287    strictneg(A, B),
  288    mixed(C, D)
  289 => (   eval(B / D, U),
  290        Res = -1.0Inf...U
  291    ;   eval(B / C, L),
  292        Res = L...1.0Inf
  293    ).
  294
  295% P1 / N
  296div(A...B, C...D, Res, _Flags),
  297    strictpos(A, B),
  298    zeroneg(C, D)
  299 => eval(A / C, U),
  300    Res = -1.0Inf...U.
  301
  302div(A...B, C...D, Res, _Flags),
  303    strictpos(A, B),
  304    negative(C, D)
  305 => eval(B / D, L),
  306    eval(A / C, U),
  307    Res = L...U.
  308
  309% P0 / N
  310div(A...B, C...D, Res, _Flags),
  311    zeropos(A, B),
  312    zeroneg(C, D)
  313 => Res = -1.0Inf...0.0.
  314
  315div(A...B, C...D, Res, _Flags),
  316    zeropos(A, B),
  317    negative(C, D)
  318 => eval(B / D, L),
  319    Res = L...0.0.
  320
  321% M / N
  322div(A...B, C...D, Res, _Flags),
  323    mixed(A, B),
  324    zeroneg(C, D)
  325 => Res = -1.0Inf...1.0Inf.
  326
  327div(A...B, C...D, Res, _Flags),
  328    mixed(A, B),
  329    negative(C, D)
  330 => eval(B / D, L),
  331    eval(A / D, U),
  332    Res = L...U.
  333
  334% N0 / N
  335div(A...B, C...D, Res, _Flags),
  336    zeroneg(A, B),
  337    zeroneg(C, D)
  338 => Res = 0.0...1.0Inf.
  339
  340div(A...B, C...D, Res, _Flags),
  341    zeroneg(A, B),
  342    negative(C, D)
  343 => eval(A / D, U),
  344    Res = 0.0...U.
  345
  346% N1 / N
  347div(A...B, C...D, Res, _Flags),
  348    strictneg(A, B),
  349    zeroneg(C, D)
  350 => eval(B / C, L),
  351    Res = L...1.0Inf.
  352
  353div(A...B, C...D, Res, _Flags),
  354    strictneg(A, B),
  355    negative(C, D)
  356 => eval(B / C, L),
  357    eval(A / D, U),
  358    Res = L...U.
  359
  360%
  361% Square root
  362%
  363% sqrt/1: "normal" behavior, returns nan for negative argument
  364% sqrt1/1: crops negative part of interval at 0
  365%
  366mono(sqrt/1, [+]).
  367
  368int_hook(sqrt1, sqrt1(...), _, []).
  369sqrt1(A...B, Res, _Flags) :-
  370    strictneg(A, B),
  371    !,
  372    Res = 1.5NaN.
  373
  374sqrt1(A...B, Res, _Flags) :-
  375    zeroneg(A, B),
  376    !,
  377    Res = 0.0.
  378
  379sqrt1(A...B, Res, _Flags) :-
  380    mixed(A, B),
  381    !,
  382    eval(sqrt(B), U),
  383    Res = 0.0...U.
  384
  385%
  386% Power
  387%
  388int_hook(^, pow0(atomic, atomic), atomic, []).
  389pow0(atomic(X), atomic(Y), atomic(Res), _Flags) :-
  390    eval(X^Y, Res).
  391
  392% Even exponent with negative base
  393int_hook(^, pow(..., atomic), ..., []).
  394pow(L...U, atomic(Exp), Res, _Flags),
  395    negative(L, U),
  396    natural(Exp),
  397    even(Exp)
  398 => eval(U^Exp, L^Exp, Res).
  399
  400% Even exponent with mixed base
  401pow(A...B, atomic(Exp), Res, _Flags),
  402    mixed(A, B),
  403    natural(Exp),
  404    even(Exp)
  405 => eval(max(A^Exp, B^Exp), U),
  406    Res = 0...U.
  407
  408% Positive also works with negative exponents
  409pow(L...U, atomic(Exp), Res, _Flags),
  410    positive(L, U),
  411    Exp < 0
  412 => eval(U^Exp, L^Exp, Res).
  413
  414% General case
  415pow(L...U, atomic(Exp), Res, _Flags),
  416    natural(Exp)
  417 => eval(L^Exp, U^Exp, Res).
  418
  419pow(L...U, atomic(Exp), Res, _Flags),
  420    positive(L, U),
  421    Exp > 0,
  422    Exp < 1
  423 => eval(L^Exp, U^Exp, Res).
  424
  425% Utility
  426even(A) :-
  427    A mod 2 =:= 0.
  428
  429natural(A) :-
  430    A >= 0,
  431    integer(A).
  432
  433%
  434% Absolute value
  435%
  436int_hook(abs, abs1(...), ..., []).
  437abs1(A...B, Res, _Flags) :-
  438    positive(A, B),
  439    !,
  440    Res = A...B.
  441
  442abs1(A...B, Res, _Flags) :-
  443    negative(A, B),
  444    !,
  445    eval(abs(A), U),
  446    eval(abs(B), L),
  447    Res = L...U.
  448
  449% mixed
  450abs1(A...B, Res, _Flags) :-
  451    !,
  452    L = 0.0,
  453    U is max(abs(A), abs(B)),
  454    Res = L...U.
  455
  456%
  457% Round
  458%
  459int_hook(round, round0(_), _, []).
  460round0(A, Res, Flags) :-
  461    interval_(round(A, atomic(0)), Res, Flags).
  462
  463int_hook(round, round1(atomic, atomic), ..., []).
  464round1(atomic(A), Dig, Res, Flags) :-
  465    round2(A...A, Dig, Res, Flags).
  466    
  467int_hook(round, round2(..., atomic), ..., []).
  468round2(A...B, atomic(Dig), Res, _Flags) :-
  469    eval(floor(A, Dig), A1),
  470    eval(ceiling(B, Dig), B1),
  471    Res = A1...B1.
  472
  473int_hook(round, round3(_, atomic), _, []).
  474round3(A, _Dig, A, _Flags).
  475
  476eval_hook(floor(A, Dig), Res) :-
  477    Mul is 10^Dig,
  478    Res is floor(A * Mul) / Mul.
  479
  480eval_hook(ceiling(A, Dig), Res) :-
  481    Mul is 10^Dig,
  482    Res is ceiling(A * Mul) / Mul.
  483
  484%
  485% sine
  486%
  487int_hook(sin, sin(...), ..., []).
  488
  489% interval extends over more than 2 max/mins
  490sin(A...B, Res, _Flags) :-
  491    A1 is A/pi - 1/2,
  492    B1 is B/pi - 1/2,
  493    B1 >= ceiling(A1) + 1,
  494    !,
  495    Res = -1...1.
  496
  497% interval extends over 1 max
  498sin(A...B, Res, _Flags) :-
  499    A1 is A / (2*pi) - 1/4,
  500    B1 is B / (2*pi) - 1/4,
  501    B1 >= ceiling(A1),
  502    !,
  503    L is min(sin(A), sin(B)),
  504    Res = L...1.
  505
  506% interval extends over 1 min
  507sin(A...B, Res, _Flags) :-
  508    A1 is A / (2*pi) + 1/4,
  509    B1 is B / (2*pi) + 1/4,
  510    B1 >= ceiling(A1),
  511    !,
  512    U is max(sin(A), sin(B)),
  513    Res = -1...U.
  514
  515% default rising
  516sin(A