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