1## Copyright 2014-2016 Oliver Heimlich 2## 3## This program is free software; you can redistribute it and/or modify 4## it under the terms of the GNU General Public License as published by 5## the Free Software Foundation; either version 3 of the License, or 6## (at your option) any later version. 7## 8## This program is distributed in the hope that it will be useful, 9## but WITHOUT ANY WARRANTY; without even the implied warranty of 10## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 11## GNU General Public License for more details. 12## 13## You should have received a copy of the GNU General Public License 14## along with this program; if not, see <http://www.gnu.org/licenses/>. 15 16## -*- texinfo -*- 17## @documentencoding UTF-8 18## @deftypemethod {@@infsup} {@var{X} =} mulrev (@var{B}, @var{C}, @var{X}) 19## @deftypemethodx {@@infsup} {@var{X} =} mulrev (@var{B}, @var{C}) 20## @deftypemethodx {@@infsup} {[@var{U}, @var{V}] =} mulrev (@var{B}, @var{C}) 21## @deftypemethodx {@@infsup} {[@var{U}, @var{V}] =} mulrev (@var{B}, @var{C}, @var{X}) 22## 23## Compute the reverse multiplication function or the two-output division. 24## 25## That is, an enclosure of all @code{x ∈ @var{X}} where 26## @code{x .* b ∈ @var{C}} for any @code{b ∈ @var{B}}. 27## 28## This function is similar to interval division @code{@var{C} ./ @var{B}}. 29## However, it treats the case 0/0 as “any real number” instead of “undefined”. 30## 31## Interval division, considered as a set, can have zero, one or two disjoint 32## connected components as a result. If called with two output parameters, 33## this function returns the components separately. @var{U} contains the 34## negative or unique component, whereas @var{V} contains the positive 35## component in cases with two components. 36## 37## Accuracy: The result is a tight enclosure. 38## 39## @example 40## @group 41## c = infsup (1); 42## b = infsup (-inf, inf); 43## [u, v] = mulrev (b, c) 44## @result{} 45## u = [-inf, 0] 46## v = [0, inf] 47## @end group 48## @end example 49## @seealso{@@infsup/times} 50## @end deftypemethod 51 52## Author: Oliver Heimlich 53## Keywords: interval 54## Created: 2014-10-13 55 56function [u, v] = mulrev (b, c, x) 57 58 if (nargin < 2 || nargin > 3) 59 print_usage (); 60 return 61 endif 62 if (nargin < 3) 63 x = infsup (-inf, inf); 64 endif 65 if (not (isa (b, "infsup"))) 66 b = infsup (b); 67 endif 68 if (not (isa (c, "infsup"))) 69 c = infsup (c); 70 endif 71 if (not (isa (x, "infsup"))) 72 x = infsup (x); 73 endif 74 75 ## Resize, if broadcasting is needed 76 if (not (size_equal (b.inf, c.inf))) 77 b.inf = ones (size (c.inf)) .* b.inf; 78 b.sup = ones (size (c.inf)) .* b.sup; 79 c.inf = ones (size (b.inf)) .* c.inf; 80 c.sup = ones (size (b.inf)) .* c.sup; 81 endif 82 83 u = v = x; 84 u.inf = v.inf = inf (size (b.inf)); 85 u.sup = v.sup = -inf (size (b.inf)); 86 87 emptyresult = b.inf == inf | c.inf == inf ... 88 | (b.inf == 0 & b.sup == 0 & (c.sup < 0 | c.inf > 0)); # x * 0 ~= 0 89 twocomponents = b.inf < 0 & b.sup > 0 & not (emptyresult) & ... 90 (c.sup < 0 | c.inf > 0); 91 onecomponent = not (twocomponents) & not (emptyresult); 92 93 u.inf(twocomponents) = -inf; 94 v.sup(twocomponents) = inf; 95 dom = twocomponents & c.inf <= 0 & c.sup >= 0; 96 u.sup(dom) = v.inf(dom) = 0; 97 dom = twocomponents & c.inf > 0; 98 if (not (isempty (dom))) 99 u.sup(dom) = mpfr_function_d ('rdivide', +inf, c.inf(dom), b.inf(dom)); 100 v.inf(dom) = mpfr_function_d ('rdivide', -inf, c.inf(dom), b.sup(dom)); 101 endif 102 dom = twocomponents & c.sup < 0; 103 if (not (isempty (dom))) 104 u.sup(dom) = mpfr_function_d ('rdivide', +inf, c.sup(dom), b.sup(dom)); 105 v.inf(dom) = mpfr_function_d ('rdivide', -inf, c.sup(dom), b.inf(dom)); 106 endif 107 108 dom = onecomponent & b.inf >= 0 & c.inf >= 0; 109 if (not (isempty (dom))) 110 b.inf(dom & b.inf == 0) = +0; 111 c.inf(dom & c.inf == 0) = +0; 112 u.inf(dom) = mpfr_function_d ('rdivide', -inf, c.inf(dom), b.sup(dom)); 113 u.sup(dom) = mpfr_function_d ('rdivide', +inf, c.sup(dom), b.inf(dom)); 114 endif 115 dom = onecomponent & b.sup <= 0 & c.inf >= 0; 116 if (not (isempty (dom))) 117 b.sup(dom & b.sup == 0) = -0; 118 c.inf(dom & c.inf == 0) = +0; 119 u.inf(dom) = mpfr_function_d ('rdivide', -inf, c.sup(dom), b.sup(dom)); 120 u.sup(dom) = mpfr_function_d ('rdivide', +inf, c.inf(dom), b.inf(dom)); 121 endif 122 dom = onecomponent & b.inf >= 0 & c.sup <= 0; 123 if (not (isempty (dom))) 124 b.inf(dom & b.inf == 0) = +0; 125 c.sup(dom & c.sup == 0) = -0; 126 u.inf(dom) = mpfr_function_d ('rdivide', -inf, c.inf(dom), b.inf(dom)); 127 u.sup(dom) = mpfr_function_d ('rdivide', +inf, c.sup(dom), b.sup(dom)); 128 endif 129 dom = onecomponent & b.sup <= 0 & c.sup <= 0; 130 if (not (isempty (dom))) 131 b.sup(dom & b.sup == 0) = -0; 132 c.sup(dom & c.sup == 0) = -0; 133 u.inf(dom) = mpfr_function_d ('rdivide', -inf, c.sup(dom), b.inf(dom)); 134 u.sup(dom) = mpfr_function_d ('rdivide', +inf, c.inf(dom), b.sup(dom)); 135 endif 136 dom = onecomponent & c.inf < 0 & c.sup > 0 & b.inf > 0; 137 if (not (isempty (dom))) 138 u.inf(dom) = mpfr_function_d ('rdivide', -inf, c.inf(dom), b.inf(dom)); 139 u.sup(dom) = mpfr_function_d ('rdivide', +inf, c.sup(dom), b.inf(dom)); 140 endif 141 dom = onecomponent & c.inf < 0 & c.sup > 0 & b.sup < 0; 142 if (not (isempty (dom))) 143 u.inf(dom) = mpfr_function_d ('rdivide', -inf, c.sup(dom), b.sup(dom)); 144 u.sup(dom) = mpfr_function_d ('rdivide', +inf, c.inf(dom), b.sup(dom)); 145 endif 146 dom = onecomponent & b.inf <= 0 & b.sup >= 0 & c.inf <= 0 & c.sup >= 0; 147 # x * 0 == 0 148 u.inf(dom) = -inf; 149 u.sup(dom) = inf; 150 151 u.inf(u.inf == 0) = -0; 152 u.sup(u.sup == 0) = +0; 153 v.inf(v.inf == 0) = -0; 154 v.sup(v.sup == 0) = +0; 155 156 ## Intersect u and v with x 157 u.inf = max (u.inf, x.inf); 158 u.sup = min (u.sup, x.sup); 159 v.inf = max (v.inf, x.inf); 160 v.sup = min (v.sup, x.sup); 161 162 if (nargout < 2) 163 u.inf(twocomponents) = min (u.inf(twocomponents), v.inf(twocomponents)); 164 u.sup(twocomponents) = max (u.sup(twocomponents), v.sup(twocomponents)); 165 emptyresult = u.inf > u.sup; 166 u.inf(emptyresult) = inf; 167 u.sup(emptyresult) = -inf; 168 else 169 empty_u = u.inf > u.sup; 170 u.inf(empty_u) = inf; 171 u.sup(empty_u) = -inf; 172 empty_v = v.inf > v.sup; 173 v.inf(empty_v) = inf; 174 v.sup(empty_v) = -inf; 175 ## It can happen that the twocomponents result has only one component, 176 ## because x is positive for example. Then, only one component shall be 177 ## returned 178 swap = twocomponents & isempty (u) & not (isempty (v)); 179 [u.inf(swap), u.sup(swap), v.inf(swap), ... 180 v.sup(swap)] = deal (v.inf(swap), v.sup(swap), u.inf(swap), u.sup(swap)); 181 endif 182 183endfunction 184 185%!#IEEE Std 1788-2015 mulRevToPair examples 186%!test 187%! [u, v] = mulrev (infsup (0), infsup (1, 2)); 188%! assert (isempty (u) & isempty (v)); 189%!test 190%! [u, v] = mulrev (infsup (0), infsup (0, 1)); 191%! assert (isentire (u) & isempty (v)); 192%!test 193%! [u, v] = mulrev (infsup (1), infsup (1, 2)); 194%! assert (eq (u, infsup (1, 2)) & isempty (v)); 195%!test 196%! [u, v] = mulrev (infsup (1, inf), infsup (1)); 197%! assert (eq (u, infsup (0, 1)) & isempty (v)); 198%!test 199%! [u, v] = mulrev (infsup (-1, 1), infsup (1, 2)); 200%! assert (eq (u, infsup (-inf, -1)) & eq (v, infsup (1, inf))); 201%!test 202%! [u, v] = mulrev (infsup (-inf, inf), infsup (1)); 203%! assert (eq (u, infsup (-inf, 0)) & eq (v, infsup (0, inf))); 204 205%!shared testdata 206%! # Load compiled test data (from src/test/*.itl) 207%! testdata = load (file_in_loadpath ("test/itl.mat")); 208 209%!test 210%! # Scalar evaluation 211%! testcases = testdata.NoSignal.infsup.mulRevToPair1; 212%! for testcase = [testcases]' 213%! assert (isequaln (... 214%! nthargout (1, 2, @mulrev, testcase.in{1}, testcase.in{2}), ... 215%! testcase.out)); 216%! endfor 217 218%!test 219%! # Vector evaluation 220%! testcases = testdata.NoSignal.infsup.mulRevToPair1; 221%! in1 = vertcat (vertcat (testcases.in){:, 1}); 222%! in2 = vertcat (vertcat (testcases.in){:, 2}); 223%! out = vertcat (testcases.out); 224%! assert (isequaln (nthargout (1, 2, @mulrev, in1, in2), out)); 225 226%!test 227%! # N-dimensional array evaluation 228%! testcases = testdata.NoSignal.infsup.mulRevToPair1; 229%! in1 = vertcat (vertcat (testcases.in){:, 1}); 230%! in2 = vertcat (vertcat (testcases.in){:, 2}); 231%! out = vertcat (testcases.out); 232%! # Reshape data 233%! i = -1; 234%! do 235%! i = i + 1; 236%! testsize = factor (numel (in1) + i); 237%! until (numel (testsize) > 2) 238%! in1 = reshape ([in1; in1(1:i)], testsize); 239%! in2 = reshape ([in2; in2(1:i)], testsize); 240%! out = reshape ([out; out(1:i)], testsize); 241%! assert (isequaln (nthargout (1, 2, @mulrev, in1, in2), out)); 242 243%!test 244%! # Scalar evaluation 245%! testcases = testdata.NoSignal.infsup.mulRevToPair2; 246%! for testcase = [testcases]' 247%! assert (isequaln (... 248%! nthargout (2, @mulrev, testcase.in{1}, testcase.in{2}), ... 249%! testcase.out)); 250%! endfor 251 252%!test 253%! # Vector evaluation 254%! testcases = testdata.NoSignal.infsup.mulRevToPair2; 255%! in1 = vertcat (vertcat (testcases.in){:, 1}); 256%! in2 = vertcat (vertcat (testcases.in){:, 2}); 257%! out = vertcat (testcases.out); 258%! assert (isequaln (nthargout (2, @mulrev, in1, in2), out)); 259 260%!test 261%! # N-dimensional array evaluation 262%! testcases = testdata.NoSignal.infsup.mulRevToPair2; 263%! in1 = vertcat (vertcat (testcases.in){:, 1}); 264%! in2 = vertcat (vertcat (testcases.in){:, 2}); 265%! out = vertcat (testcases.out); 266%! # Reshape data 267%! i = -1; 268%! do 269%! i = i + 1; 270%! testsize = factor (numel (in1) + i); 271%! until (numel (testsize) > 2) 272%! in1 = reshape ([in1; in1(1:i)], testsize); 273%! in2 = reshape ([in2; in2(1:i)], testsize); 274%! out = reshape ([out; out(1:i)], testsize); 275%! assert (isequaln (nthargout (2, @mulrev, in1, in2), out)); 276 277%!test 278%! # Scalar evaluation 279%! testcases = testdata.NoSignal.infsup.mulRev; 280%! for testcase = [testcases]' 281%! assert (isequaln (... 282%! mulrev (testcase.in{1}, testcase.in{2}), ... 283%! testcase.out)); 284%! endfor 285 286%!test 287%! # Vector evaluation 288%! testcases = testdata.NoSignal.infsup.mulRev; 289%! in1 = vertcat (vertcat (testcases.in){:, 1}); 290%! in2 = vertcat (vertcat (testcases.in){:, 2}); 291%! out = vertcat (testcases.out); 292%! assert (isequaln (mulrev (in1, in2), out)); 293 294%!test 295%! # N-dimensional array evaluation 296%! testcases = testdata.NoSignal.infsup.mulRev; 297%! in1 = vertcat (vertcat (testcases.in){:, 1}); 298%! in2 = vertcat (vertcat (testcases.in){:, 2}); 299%! out = vertcat (testcases.out); 300%! # Reshape data 301%! i = -1; 302%! do 303%! i = i + 1; 304%! testsize = factor (numel (in1) + i); 305%! until (numel (testsize) > 2) 306%! in1 = reshape ([in1; in1(1:i)], testsize); 307%! in2 = reshape ([in2; in2(1:i)], testsize); 308%! out = reshape ([out; out(1:i)], testsize); 309%! assert (isequaln (mulrev (in1, in2), out)); 310 311%!test 312%! # Scalar evaluation 313%! testcases = testdata.NoSignal.infsup.mulRevTen; 314%! for testcase = [testcases]' 315%! assert (isequaln (... 316%! mulrev (testcase.in{1}, testcase.in{2}, testcase.in{3}), ... 317%! testcase.out)); 318%! endfor 319 320%!test 321%! # Vector evaluation 322%! testcases = testdata.NoSignal.infsup.mulRevTen; 323%! in1 = vertcat (vertcat (testcases.in){:, 1}); 324%! in2 = vertcat (vertcat (testcases.in){:, 2}); 325%! in3 = vertcat (vertcat (testcases.in){:, 3}); 326%! out = vertcat (testcases.out); 327%! assert (isequaln (mulrev (in1, in2, in3), out)); 328 329%!test 330%! # N-dimensional array evaluation 331%! testcases = testdata.NoSignal.infsup.mulRevTen; 332%! in1 = vertcat (vertcat (testcases.in){:, 1}); 333%! in2 = vertcat (vertcat (testcases.in){:, 2}); 334%! in3 = vertcat (vertcat (testcases.in){:, 3}); 335%! out = vertcat (testcases.out); 336%! # Reshape data 337%! i = -1; 338%! do 339%! i = i + 1; 340%! testsize = factor (numel (in1) + i); 341%! until (numel (testsize) > 2) 342%! in1 = reshape ([in1; in1(1:i)], testsize); 343%! in2 = reshape ([in2; in2(1:i)], testsize); 344%! in3 = reshape ([in3; in3(1:i)], testsize); 345%! out = reshape ([out; out(1:i)], testsize); 346%! assert (isequaln (mulrev (in1, in2, in3), out)); 347