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