1:- module(interval, [interval/2, op(150, xfx, ...)]).    2
    3:- multifile int_hook/2.    4:- multifile eval_hook/2.    5:- multifile mono/2.    6
    7:- discontiguous interval/2.    8:- discontiguous int_hook/2.    9
   10:- set_prolog_flag(float_overflow, infinity).   11:- set_prolog_flag(float_undefined, nan).   12:- set_prolog_flag(float_zero_div, infinity).

This module adds interval arithemtic to Prolog. An interval is represented as L...U, where L stands for the lower bound and U for the upper bound. If the upper bound is a negative number, it has to be written with an additional space, e.g., 1... -2, or in the infix notation, ...(1, -2). The interval/2 parses and evaluates the arithemtic expression with such intervals to a result.

 interval(+A, ?Res)

Evalutes an expression to an interval. If the first argument is already an interval, no evaluation is performed. Supported operations:

• Basic arithemtic: addition '+', subtraction '-', division '/', multiplication '*'
• Square root for positive interval: 'interval(sqrt(X), Res)'
• Square root for negative or mixed interval: 'interval(sqrt1(X), Res)'
• Power: 'interval(X^atomic(N), Res)' with N being a natural number
• Absolute value: 'interval(abs(X), Res)'
• Comparison: '>', '<', '>=', '=<', '=\=', '=:='
• Rounding: 'interval(round(1.356...1.634, atomic(2)), Res)'

Arguments:

A	- is the expression to be evaluted.
Res	- is the result.

   37%
   38% If 1st argument already is an interval, do not do anything
   39%
   40interval(L...U, Res)
   41 => Res = L...U.
   42
   43% Force translation of atom to interval
   44/* interval(atomic(A), L...U)
   45 => eval(A, Res),
   46    L = Res,
   47    U = Res. */
   48
   49interval(atomic(A), Res)
   50 => eval(A, R),
   51    Res = atomic(R).
   52
   53%
   54% Hook for custom interval functions
   55%
   56% 1. Declare function with interval:int_hook(Name, pred_name(Args))
   57%    Args specify the type of the arguments, ... for interval and atomic for numbers.
   58% 2. Calculate result with interval:pred_name(Args, Res)
   59%    Args are the argument types as defined in the hook with variable names.
   60% 
   61%
   62% see below example for (/)/2.
   63interval(Expr, Res),
   64    compound(Expr),
   65    compound_name_arguments(Expr, Name, Args),
   66    int_hook(Name, Mask),
   67    compound_name_arguments(Mask, Fun, Args1),
   68    maplist(instantiate, Args1, Args2),
   69    maplist(interval, Args, Args2)
   70 => compound_name_arguments(Goal, Fun, Args2),
   71    call(Goal, Res).
   72
   73instantiate(atomic, atomic(_)).
   74instantiate(..., _..._).
   75instantiate(ci, ci(_, _)).
   76%
   77% Monotonically behaving functions
   78%
   79% +: increasing
   80% -: decreasing
   81% *: increasing or decreasing
   82% **: commutative, all *
   83% /: symbolic argument or flag
   84mono((+)/1, [+]).
   85mono((+)/2, [+, +]).
   86mono((-)/1, [-]).
   87mono((-)/2, [+, -]).
   88mono((*)/2, **).
   89mono((^)/2, [*, /]).
   90
   91% special case: multiplication ([*, *], commutative)
   92interval(Expr, Res),
   93    compound(Expr),
   94    compound_name_arity(Expr, Name, Arity),
   95    mono(Name/Arity, **)
   96 => compound_name_arguments(Expr, Name, Args),
   97    maplist(interval, Args, Args1),
   98    findall(R, both(Name, Args1, R), Bounds),
   99    min_list(Bounds, L),
  100    max_list(Bounds, U),
  101    Res = L...U.
  102
  103% general case
  104interval(Expr, Res),
  105    compound(Expr),
  106    compound_name_arity(Expr, Name, Arity),
  107    mono(Name/Arity, Dir)
  108 => compound_name_arguments(Expr, Name, Args),
  109    maplist(interval, Args, Args1),
  110    findall(R, lower(Dir, Name, Args1, R), Lower),
  111    min_list(Lower, L),
  112    findall(R, upper(Dir, Name, Args1, R), Upper),
  113    max_list(Upper, U),
  114    Res = L...U.
  115
  116lower(Dir, Name, Args, Res) :-
  117    maplist(lower, Dir, Args, Lower),
  118    Expr =.. [Name | Lower],
  119    eval(Expr, Res).
  120
  121upper(Dir, Name, Args, Res) :-
  122    maplist(upper, Dir, Args, Upper),
  123    Expr =.. [Name | Upper],
  124    eval(Expr, Res).
  125
  126both(Name, Args, Res) :-
  127    maplist(lower(*), Args, Lower),
  128    Expr =.. [Name | Lower],
  129    eval(Expr, Res).
  130
  131% Obtain lower and upper bounds
  132lower(+, A..._, L)
  133 => L = A.
  134
  135lower(-, _...A, L)
  136 => L = A.
  137
  138lower(*, A...B, L)
  139 => L = A ; L = B.
  140
  141lower(_, A, L) % either / or A not interval
  142 => L = A.
  143
  144upper(+, _...B, U)
  145 => U = B.
  146
  147upper(-, A..._, U)
  148 => U = A.
  149
  150upper(*, A...B, U)
  151 => U = A ; U = B.
  152
  153upper(_, A, U)
  154 => U = A.
  155
  156%
  157% Arithmetic functions for single numbers
  158%
  159% Define hooks for R functions etc.
  160%
  161eval(Expr, Res),
  162    eval_hook(Expr, R)
  163 => Res = R.
  164
  165eval(X, Res)
  166 => Res is X.
  167
  168%
  169% Comparison
  170%
  171int_hook(<, less1(atomic, atomic)).
  172
  173less1(atomic(A), atomic(B), Res) :-
  174    A < B,
  175    !,
  176    Res = true.
  177
  178less1(atomic(_) < atomic(_), Res) :-
  179    !,
  180    Res = false.
  181
  182int_hook(<, less2(..., ...)).
  183
  184less2(_...A2, B1..._, Res) :-
  185    A2 < B1,
  186    !,
  187    Res = true.
  188
  189less2(_..._, _..._, false2).
  190
  191int_hook(=<, less_eq(..., ...)).
  192
  193less_eq(A1..._, _...B2, Res) :-
  194    A1 =< B2,
  195    !,
  196    Res = true.
  197
  198less_eq(_..._, _..._, false).
  199
  200int_hook(>, great(..., ...)).
  201great(A1..._, _...B2, Res) :-
  202    A1 > B2,
  203    !,
  204    Res = true.
  205
  206great(_..._, _..._, false).
  207
  208int_hook(>=, great_eq(..., ...)).
  209great_eq(_...A2, B1..._, Res) :-
  210    A2 >= B1,
  211    !,
  212    Res = true.
  213
  214great_eq(_..._, _..._, false).
  215
  216
  217int_hook(=\=, not_eq(..., ...)).
  218not_eq(A...B, C...D, Res) :-
  219    (   less2(A...B, C...D, true)
  220    ;   great(A...B, C...D, true)
  221    ), !,
  222    Res = true.
  223
  224not_eq(_..._, _..._, false).
  225
  226
  227int_hook(=:=, eq(..., ...)).
  228eq(A...B, C...D, Res) :-
  229    less_eq(A...B, C...D, true),
  230    great_eq(A...B, C...D, true),
  231    !,
  232    Res = true.
  233
  234eq(_..._, _..._, false). 
  235
  236%
  237% Division
  238%
  239int_hook(/, div1(..., ...)).
  240div1(A...B, C...D, Res) :-
  241    !,
  242    div(A...B, C...D, Res).
  243
  244int_hook(/, div2(..., atomic)).
  245div2(A1...A2, atomic(B), Res) :-
  246    !,
  247    div(A1...A2, B...B, Res).
  248
  249int_hook(/, div3(atomic, ...)).
  250div3(atomic(A), B1...B2, Res) :-
  251    !,
  252    div(A...A, B1...B2, Res).
  253
  254int_hook(/, div4(atomic, atomic)).
  255div4(atomic(A), atomic(B), Res) :-
  256    Res is A / B.
  257
  258% Hickey Figure 1
  259mixed(L, U) :-
  260    L < 0,
  261    U > 0.
  262
  263positive(L, U) :-
  264    L >= 0,
  265    U > 0.
  266
  267zeropos(L, U) :-
  268    L =:= 0,
  269    U > 0.
  270
  271strictpos(L, _) :-
  272    L > 0.
  273
  274negative(L, U) :-
  275    L < 0,
  276    U =< 0.
  277
  278zeroneg(L, U) :-
  279    L < 0,
  280    U =:= 0.
  281
  282strictneg(_, U) :-
  283    U < 0.
  284
  285%
  286% Hickey Theorem 8 and Figure 4
  287%
  288% P1 / P (special case, then general case)
  289div(A...B, C...D, Res),
  290    strictpos(A, B),
  291    zeropos(C, D)
  292 => eval(A / D, L),
  293    Res = L...1.0Inf.
  294
  295div(A...B, C...D, Res),
  296    strictpos(A, B),
  297    positive(C, D)
  298 => eval(A / D, L),
  299    eval(B / C, U),
  300    Res = L...U.
  301
  302% P0 / P
  303div(A...B, 0.0...D, Res),
  304    zeropos(A, B),
  305    positive(0.0, D)
  306 => Res = 0.0...1.0Inf.
  307
  308div(A...B, C...D, Res),
  309    zeropos(A, B),
  310    positive(C, D)
  311 => eval(B / C, U),
  312    Res = 0.0...U.
  313
  314% M / P
  315div(A...B, 0.0...D, Res),
  316    mixed(A, B),
  317    positive(0.0, D)
  318 => Res = -1.0Inf...1.0Inf.
  319
  320div(A...B, C...D, Res),
  321    mixed(A, B),
  322    positive(C, D)
  323 => eval(A / C, L),
  324    eval(B / C, U),
  325    Res = L...U.
  326
  327% N0 / P
  328div(A...B, 0.0...D, Res),
  329    zeroneg(A, B),
  330    positive(0.0, D)
  331 => Res = -1.0Inf...0.0.
  332
  333div(A...B, C...D, Res),
  334    zeroneg(A, B),
  335    positive(C, D)
  336 => eval(A / C, L),
  337    Res = L...0.0.
  338
  339% N1 / P
  340div(A...B, 0.0...D, Res),
  341    strictneg(A, B),
  342    positive(0.0, D)
  343 => eval(B / D, U),
  344    Res = -1.0Inf...U.
  345
  346div(A...B, C...D, Res),
  347    strictneg(A, B),
  348    positive(C, D)
  349 => eval(A / C, L),
  350    eval(B / D, U),
  351    Res = L...U.
  352
  353% P1 / M (2 solutions)
  354div(A...B, C...D, Res),
  355    strictpos(A, B),
  356    mixed(C, D)
  357 => (   eval(A / C, U),
  358        Res = -1.0Inf...U
  359    ;   eval(A / D, L),
  360        Res = L...1.0Inf
  361    ).
  362
  363% P0 / M
  364div(A...B, C...D, Res),
  365    zeropos(A, B),
  366    mixed(C, D)
  367 => Res = -1.0Inf...1.0Inf.
  368
  369% M / M
  370div(A...B, C...D, Res),
  371    mixed(A, B),
  372    mixed(C, D)
  373 => Res = -1.0Inf...1.0Inf.
  374
  375% N0 / M
  376div(A...B, C...D, Res),
  377    zeroneg(A, B),
  378    mixed(C, D)
  379 => Res = -1.0Inf...1.0Inf.
  380
  381% N1 / M (2 solutions)
  382div(A...B, C...D, Res),
  383    strictneg(A, B),
  384    mixed(C, D)
  385 => (   eval(B / D, U),
  386        Res = -1.0Inf...U
  387    ;   eval(B / C, L),
  388        Res = L...1.0Inf
  389    ).
  390
  391% P1 / N
  392div(A...B, C...0.0, Res),
  393    strictpos(A, B),
  394    negative(C, 0.0)
  395 => eval(A / C, U),
  396    Res = -1.0Inf...U.
  397
  398div(A...B, C...D, Res),
  399    strictpos(A, B),
  400    negative(C, D)
  401 => eval(B / D, L),
  402    eval(A / C, U),
  403    Res = L...U.
  404
  405% P0 / N
  406div(A...B, C...0.0, Res),
  407    zeropos(A, B),
  408    negative(C, 0.0)
  409 => Res = -1.0Inf...0.0.
  410
  411div(A...B, C...D, Res),
  412    zeropos(A, B),
  413    negative(C, D)
  414 => eval(B / D, L),
  415    Res = L...0.0.
  416
  417% M / N
  418div(A...B, C...0.0, Res),
  419    mixed(A, B),
  420    negative(C, 0.0)
  421 => Res = -1.0Inf...1.0Inf.
  422
  423div(A...B, C...D, Res),
  424    mixed(A, B),
  425    negative(C, D)
  426 => eval(B / D, L),
  427    eval(A / D, U),
  428    Res = L...U.
  429
  430% N0 / N
  431div(A...B, C...0.0, Res),
  432    zeroneg(A, B),
  433    negative(C, 0.0)
  434 => Res = 0.0...1.0Inf.
  435
  436div(A...B, C...D, Res),
  437    zeroneg(A, B),
  438    negative(C, D)
  439 => eval(A / D, U),
  440    Res = 0.0...U.
  441
  442% N1 / N
  443div(A...B, C...0.0, Res),
  444    strictneg(A, B),
  445    negative(C, 0.0)
  446 => eval(B / C, L),
  447    Res = L...1.0Inf.
  448
  449div(A...B, C...D, Res),
  450    strictneg(A, B),
  451    negative(C, D)
  452 => eval(B / C, L),
  453    eval(A / D, U),
  454    Res = L...U.
  455
  456%
  457% Square root
  458%
  459% sqrt/1: "normal" behavior, returns nan for negative argument
  460% sqrt1/1: crops negative part of interval at 0
  461%
  462mono(sqrt/1, [+]).
  463
  464int_hook(sqrt1, sqrt1(...)).
  465sqrt1(A...B, Res) :-
  466    strictneg(A, B),
  467    !,
  468    Res = 1.5NaN.
  469
  470sqrt1(A...B, Res) :-
  471    zeroneg(A, B),
  472    !,
  473    Res = 0.0.
  474
  475sqrt1(A...B, Res) :-
  476    mixed(A, B),
  477    !,
  478    eval(sqrt(B), U),
  479    Res = 0.0...U.
  480%
  481% Power
  482%
  483% Even exponent with negative base
  484int_hook((^), pow(..., atomic)).
  485pow(L...U, atomic(Exp), Res),
  486    negative(L, U),
  487    even(Exp),
  488    natural(Exp)
  489 => eval(U^Exp, L^Exp, Res).
  490
  491% Even exponent with mixed base
  492pow(L...U, atomic(Exp), Res),
  493    mixed(L, U),
  494    even(Exp),
  495    natural(Exp)
  496 => eval(max(L^Exp, U^Exp), Upper),
  497    Res = 0...Upper.
  498
  499% General case
  500pow(L...U, atomic(Exp), Res),
  501    natural(Exp)
  502 => eval(L^Exp, U^Exp, Res).
  503
  504% Utility
  505even(A) :-
  506    A mod 2 =:= 0.
  507
  508natural(A) :-
  509    A >=0,
  510    integer(A).
  511
  512%
  513% Absolute value
  514%
  515int_hook(abs, abs1(...)).
  516abs1(A...B, Res) :-
  517    positive(A, B),
  518    !,
  519    Res = A...B.
  520
  521abs1(A...B, Res) :-
  522    negative(A, B),
  523    !,
  524    eval(abs(A), U),
  525    eval(abs(B), L),
  526    Res = L...U.
  527
  528% mixed
  529abs1(A...B, Res) :-
  530    !,
  531    L = 0.0,
  532    U is max(abs(A), abs(B)),
  533    Res = L...U.
  534
  535%
  536% round interval
  537%
  538int_hook(round, round1(..., atomic)).
  539round1(A...B, atomic(Dig), Res) :-
  540    eval(floor(A, Dig), A1),
  541    eval(ceiling(B, Dig), B1),
  542    Res = A1...B1.
  543
  544eval_hook(floor(A, Dig), Res) :-
  545    Mul is 10^Dig,
  546    Res is floor(A * Mul) / Mul.
  547
  548eval_hook(ceiling(A, Dig), Res) :-
  549    Mul is 10^Dig,
  550    Res is ceiling(A * Mul) / Mul.
  551
  552% For convenience
  553eval(Expr1, Expr2, L ... U) :-
  554    interval:eval(Expr1, L),
  555    interval:eval(Expr2, U)