1:- module(interval, [interval/2, interval/3, op(150, xfx, ...)]).    2
    3:- multifile int_hook/1.    4:- multifile int_hook/3.    5:- multifile eval_hook/2.    6:- multifile mono/2.    7
    8:- discontiguous interval/3.    9
   10:- set_prolog_flag(float_overflow, infinity).   11:- set_prolog_flag(float_undefined, nan).   12:- set_prolog_flag(float_zero_div, infinity).   13
   14%
   15% User-level
   16%
   17interval(A, Res) :-
   18    interval(A, Res, []).
   19
   20% For maplist
   21interval_(Opt, A, Res) :-
   22    interval(A, Res, Opt).
   23
   24%
   25% If 1st argument already is an interval, do not do anything
   26%
   27interval(L...U, Res, _)
   28 => Res = L...U.
   29
   30%
   31% Hook for custom interval functions
   32%
   33% 1. Declare function with interval:int_hook(Name/Arity)
   34% 2. Calculate result with interval:int_hook(Expr, Res)
   35%
   36% see below example for (/)/2.
   37interval(Expr, Res, Opt),
   38    compound(Expr),
   39    compound_name_arity(Expr, Name, Arity),
   40    int_hook(Name/Arity)
   41 => compound_name_arguments(Expr, Name, Args),
   42    maplist(interval_(Opt), Args, Args1),
   43    compound_name_arguments(Expr1, Name, Args1),
   44    int_hook(Expr1, Res, Opt).
   45
   46%
   47% Default behavior for atoms
   48%
   49% force result to interval (deprecated)
   50interval(A, L...U, _),
   51    atomic(A)
   52 => eval(A, R),
   53    L = R,
   54    U = R.
   55
   56% work with atom
   57interval(A, Res, _),
   58    atomic(A)
   59 => eval(A, Res).
   60
   61%
   62% Monotonically behaving functions
   63%
   64% +: increasing
   65% -: decreasing
   66% *: increasing or decreasing
   67% **: commutative, all *
   68% /: symbolic argument or flag
   69%
   70mono((+)/1, [+]).
   71mono((+)/2, [+, +]).
   72mono((-)/1, [-]).
   73mono((-)/2, [+, -]).
   74mono((*)/2, **).
   75
   76% special case: multiplication ([*, *], commutative)
   77interval(Expr, Res, Opt),
   78    compound(Expr),
   79    compound_name_arity(Expr, Name, Arity),
   80    mono(Name/Arity, **)
   81 => Expr =.. [ Name | Args],
   82    maplist(interval_(Opt), Args, Args1),
   83    findall(R, both(Name, Args1, R), Bounds),
   84    min_list(Bounds, L),
   85    max_list(Bounds, U),
   86    Res = L...U.
   87
   88both(Name, Args, Res) :-
   89    maplist(lower(*), Args, Lower),
   90    Expr =.. [Name | Lower],
   91    eval(Expr, Res).
   92
   93% general case
   94interval(Expr, Res, Opt),
   95    compound(Expr),
   96    compound_name_arity(Expr, Name, Arity),
   97    mono(Name/Arity, Dir)
   98 => Expr =.. [ Name | Args],
   99    maplist(interval_(Opt), Args, Args1),
  100    findall(R, lower(Dir, Name, Args1, R), Lower),
  101    min_list(Lower, L),
  102    findall(R, upper(Dir, Name, Args1, R), Upper),
  103    max_list(Upper, U),
  104    Res = L...U.
  105
  106lower(Dir, Name, Args, Res) :-
  107    maplist(lower, Dir, Args, Lower),
  108    Expr =.. [Name | Lower],
  109    eval(Expr, Res).
  110
  111upper(Dir, Name, Args, Res) :-
  112    maplist(upper, Dir, Args, Upper),
  113    Expr =.. [Name | Upper],
  114    eval(Expr, Res).
  115
  116% Obtain lower and upper bounds
  117lower(+, A..._, L)
  118 => L = A.
  119
  120lower(-, _...A, L)
  121 => L = A.
  122
  123lower(*, A...B, L)
  124 => L = A ; L = B.
  125
  126lower(_, A, L) % either / or A not interval
  127 => L = A.
  128
  129upper(+, _...B, U)
  130 => U = B.
  131
  132upper(-, A..._, U)
  133 => U = A.
  134
  135upper(*, A...B, U)
  136 => U = A ; U = B.
  137
  138upper(_, A, U)
  139 => U = A.
  140
  141%
  142% Arithmetic functions for single numbers
  143%
  144% Define hooks for R functions etc.
  145%
  146eval(Expr, Res),
  147    eval_hook(Expr, R)
  148 => Res = R.
  149
  150eval(X, Res)
  151 => Res is X.
  152
  153%
  154% Comparison
  155%
  156int_hook((<)/2).
  157int_hook(_...A2 < B1..._, Res, _) :-
  158    A2 < B1,
  159    !,
  160    Res = true.
  161
  162int_hook(_..._ < _..._, false, _).
  163
  164int_hook((=<)/2).
  165int_hook(A1..._ =< _...B2, Res, _) :-
  166    A1 =< B2,
  167    !,
  168    Res = true.
  169
  170int_hook(_..._ =< _..._, false, _).
  171
  172int_hook((>)/2).
  173int_hook(A1..._ > _...B2, Res, _) :-
  174    A1 > B2,
  175    !,
  176    Res = true.
  177
  178int_hook(_..._ > _..._, false, _).
  179
  180int_hook((>=)/2).
  181int_hook(_...A2 >= B1..._, Res, _) :-
  182    A2 >= B1,
  183    !,
  184    Res = true.
  185
  186int_hook(_..._ >= _..._, false, _).
  187
  188int_hook(A =\= B, Res, Opt) :-
  189    (   interval(A < B, true, Opt)
  190    ;   interval(A > B, true, Opt)
  191    ), !,
  192    Res = true.
  193
  194int_hook(_..._ =\= _..._, false, _).
  195
  196int_hook(A =:= B, Res, Opt) :-
  197    interval(A =< B, true, Opt),
  198    interval(A >= B, true, Opt),
  199    !,
  200    Res = true.
  201
  202int_hook(_..._ =:= _..._, false, _).
  203
  204%
  205% Division
  206%
  207int_hook((/)/2).
  208int_hook(A1...A2 / B1...B2, Res, _) :-
  209    div(A1...A2, B1...B2, Res).
  210
  211int_hook(A1...A2 / B, Res, _) :-
  212    div(A1...A2, B...B, Res).
  213
  214int_hook(A / B1...B2, Res, _) :-
  215    div(A...A, B1...B2, Res).
  216
  217int_hook(A / B, Res, _) :-
  218    Res is A / B.
  219
  220% Hickey Figure 1
  221mixed(L, U) :-
  222    L < 0,
  223    U > 0.
  224
  225positive(L, U) :-
  226    L >= 0,
  227    U > 0.
  228
  229zeropos(L, U) :-
  230    L =:= 0,
  231    U > 0.
  232
  233strictpos(L, _) :-
  234    L > 0.
  235
  236negative(L, U) :-
  237    L < 0,
  238    U =< 0.
  239
  240zeroneg(L, U) :-
  241    L < 0,
  242    U =:= 0.
  243
  244strictneg(_, U) :-
  245    U < 0.
  246
  247%
  248% Hickey Theorem 8 and Figure 4
  249%
  250% P1 / P (special case, then general case)
  251div(A...B, C...D, Res),
  252    strictpos(A, B),
  253    zeropos(C, D)
  254 => eval(A / D, L),
  255    Res = L...1.0Inf.
  256
  257div(A...B, C...D, Res),
  258    strictpos(A, B),
  259    positive(C, D)
  260 => eval(A / D, L),
  261    eval(B / C, U),
  262    Res = L...U.
  263
  264% P0 / P
  265div(A...B, 0.0...D, Res),
  266    zeropos(A, B),
  267    positive(0.0, D)
  268 => Res = 0.0...1.0Inf.
  269
  270div(A...B, C...D, Res),
  271    zeropos(A, B),
  272    positive(C, D)
  273 => eval(B / C, U),
  274    Res = 0.0...U.
  275
  276% M / P
  277div(A...B, 0.0...D, Res),
  278    mixed(A, B),
  279    positive(0.0, D)
  280 => Res = -1.0Inf...1.0Inf.
  281
  282div(A...B, C...D, Res),
  283    mixed(A, B),
  284    positive(C, D)
  285 => eval(A / C, L),
  286    eval(B / C, U),
  287    Res = L...U.
  288
  289% N0 / P
  290div(A...B, 0.0...D, Res),
  291    zeroneg(A, B),
  292    positive(0.0, D)
  293 => Res = -1.0Inf...0.0.
  294
  295div(A...B, C...D, Res),
  296    zeroneg(A, B),
  297    positive(C, D)
  298 => eval(A / C, L),
  299    Res = L...0.0.
  300
  301% N1 / P
  302div(A...B, 0.0...D, Res),
  303    strictneg(A, B),
  304    positive(0.0, D)
  305 => eval(B / D, U),
  306    Res = -1.0Inf...U.
  307
  308div(A...B, C...D, Res),
  309    strictneg(A, B),
  310    positive(C, D)
  311 => eval(A / C, L),
  312    eval(B / D, U),
  313    Res = L...U.
  314
  315% P1 / M (2 solutions)
  316div(A...B, C...D, Res),
  317    strictpos(A, B),
  318    mixed(C, D)
  319 => (   eval(A / C, U),
  320        Res = -1.0Inf...U
  321    ;   eval(A / D, L),
  322        Res = L...1.0Inf
  323    ).
  324
  325% P0 / M
  326div(A...B, C...D, Res),
  327    zeropos(A, B),
  328    mixed(C, D)
  329 => Res = -1.0Inf...1.0Inf.
  330
  331% M / M
  332div(A...B, C...D, Res),
  333    mixed(A, B),
  334    mixed(C, D)
  335 => Res = -1.0Inf...1.0Inf.
  336
  337% N0 / M
  338div(A...B, C...D, Res),
  339    zeroneg(A, B),
  340    mixed(C, D)
  341 => Res = -1.0Inf...1.0Inf.
  342
  343% N1 / M (2 solutions)
  344div(A...B, C...D, Res),
  345    strictneg(A, B),
  346    mixed(C, D)
  347 => (   eval(B / D, U),
  348        Res = -1.0Inf...U
  349    ;   eval(B / C, L),
  350        Res = L...1.0Inf
  351    ).
  352
  353% P1 / N
  354div(A...B, C...0.0, Res),
  355    strictpos(A, B),
  356    negative(C, 0.0)
  357 => eval(A / C, U),
  358    Res = -1.0Inf...U.
  359
  360div(A...B, C...D, Res),
  361    strictpos(A, B),
  362    negative(C, D)
  363 => eval(B / D, L),
  364    eval(A / C, U),
  365    Res = L...U.
  366
  367% P0 / N
  368div(A...B, C...0.0, Res),
  369    zeropos(A, B),
  370    negative(C, 0.0)
  371 => Res = -1.0Inf...0.0.
  372
  373div(A...B, C...D, Res),
  374    zeropos(A, B),
  375    negative(C, D)
  376 => eval(B / D, L),
  377    Res = L...0.0.
  378
  379% M / N
  380div(A...B, C...0.0, Res),
  381    mixed(A, B),
  382    negative(C, 0.0)
  383 => Res = -1.0Inf...1.0Inf.
  384
  385div(A...B, C...D, Res),
  386    mixed(A, B),
  387    negative(C, D)
  388 => eval(B / D, L),
  389    eval(A / D, U),
  390    Res = L...U.
  391
  392% N0 / N
  393div(A...B, C...0.0, Res),
  394    zeroneg(A, B),
  395    negative(C, 0.0)
  396 => Res = 0.0...1.0Inf.
  397
  398div(A...B, C...D, Res),
  399    zeroneg(A, B),
  400    negative(C, D)
  401 => eval(A / D, U),
  402    Res = 0.0...U.
  403
  404% N1 / N
  405div(A...B, C...0.0, Res),
  406    strictneg(A, B),
  407    negative(C, 0.0)
  408 => eval(B / C, L),
  409    Res = L...1.0Inf.
  410
  411div(A...B, C...D, Res),
  412    strictneg(A, B),
  413    negative(C, D)
  414 => eval(B / C, L),
  415    eval(A / D, U),
  416    Res = L...U.
  417
  418%
  419% Square root
  420%
  421% sqrt/1: "normal" behavior, returns nan for negative argument
  422% sqrt1/1: crops negative part of interval at 0
  423%
  424mono(sqrt/1, [+]).
  425
  426int_hook(sqrt1/1).
  427int_hook(sqrt1(A...B), Res, _) :-
  428    strictneg(A, B),
  429    !,
  430    Res = 1.5NaN.
  431
  432int_hook(sqrt1(A...B), Res, _) :-
  433    zeroneg(A, B),
  434    !,
  435    Res = 0.0.
  436
  437int_hook(sqrt1(A...B), Res, Opt) :-
  438    mixed(A, B),
  439    !,
  440    interval(sqrt(0...B), Res, Opt).
  441
  442int_hook(sqrt1(X), Res, Opt) :-
  443    interval(sqrt(X), Res, Opt).
  444
  445%
  446% Absolute value
  447%
  448int_hook(abs/1).
  449int_hook(abs(A...B), Res, _) :-
  450    positive(A, B),
  451    !,
  452    Res = A...B.
  453
  454int_hook(abs(A...B), Res, _) :-
  455    negative(A, B),
  456    !,
  457    eval(abs(A), U),
  458    eval(abs(B), L),
  459    Res = L...U.
  460
  461% mixed
  462int_hook(abs(A...B), Res, _) :-
  463    !,
  464    L = 0.0,
  465    U is max(abs(A), abs(B)),
  466    Res = L...U.
  467
  468%
  469% round interval
  470%
  471int_hook(round/1).
  472int_hook(round(A...B), Res, Opt) :-
  473    option(digit(Dig), Opt, 2),
  474    eval(floor(A, Dig), A1),
  475    eval(ceiling(B, Dig), B1),
  476    Res = A1...B1.
  477
  478eval_hook(floor(A, Dig), Res) :-
  479    Mul is 10^Dig,
  480    Res is floor(A * Mul) / Mul.
  481
  482eval_hook(ceiling(A, Dig), Res) :-
  483    Mul is 10^Dig,
  484    Res is ceiling(A * Mul) / Mul