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} =} cosrev (@var{C}, @var{X}) 19## @deftypemethodx {@@infsup} {@var{X} =} cosrev (@var{C}) 20## 21## Compute the reverse cosine function. 22## 23## That is, an enclosure of all @code{x ∈ @var{X}} where 24## @code{cos (x) ∈ @var{C}}. 25## 26## Accuracy: The result is a valid enclosure. 27## 28## @example 29## @group 30## cosrev (infsup (0), infsup (6, 9)) 31## @result{} ans ⊂ [7.8539, 7.854] 32## @end group 33## @end example 34## @seealso{@@infsup/cos} 35## @end deftypemethod 36 37## Author: Oliver Heimlich 38## Keywords: interval 39## Created: 2014-10-19 40 41function result = cosrev (c, x) 42 43 if (nargin > 2) 44 print_usage (); 45 return 46 endif 47 48 if (nargin < 2) 49 x = infsup (-inf, inf); 50 endif 51 if (not (isa (c, "infsup"))) 52 c = infsup (c); 53 endif 54 if (not (isa (x, "infsup"))) 55 x = infsup (x); 56 endif 57 58 arccosine = acos (c); 59 result = x; 60 61 ## Resize, if broadcasting is needed 62 if (not (size_equal (arccosine.inf, result.inf))) 63 arccosine.inf = ones (size (result.inf)) .* arccosine.inf; 64 arccosine.sup = ones (size (result.inf)) .* arccosine.sup; 65 result.inf = ones (size (arccosine.inf)) .* result.inf; 66 result.sup = ones (size (arccosine.inf)) .* result.sup; 67 endif 68 69 result.inf(isempty (arccosine)) = inf; 70 result.sup(isempty (arccosine)) = -inf; 71 72 idx.type = '()'; 73 74 persistent pi = infsup ("pi"); 75 select = not (isempty (result)) ... 76 & not (subset (infsup (0, pi.sup), arccosine)); 77 78 if (any (select(:))) 79 ## Find a smaller upper bound for x, if the restriction from c allows it 80 u = inf (size (result.inf)); 81 select_u = select & result.sup < inf; 82 ## Find n, such that result.sup is within a distance of pi/2 83 ## around (n + 1/2) * pi. 84 n = result.sup; 85 n(select_u) = floor (sup (n(select_u) ./ pi)); 86 arccosineshifted = arccosine; 87 idx.subs = {(select_u & rem (n, 2) == 0)}; 88 arccosineshifted = subsasgn (arccosineshifted, idx, ... 89 subsref (arccosine, idx) + ... 90 subsref (n, idx) .* pi); 91 idx.subs = {(select_u & rem (n, 2) ~= 0)}; 92 arccosineshifted = subsasgn (arccosineshifted, idx, ... 93 (infsup (subsref (n, idx)) + 1) .* pi - ... 94 subsref (arccosine, idx)); 95 overlapping = not (isempty (intersect (result, arccosineshifted))); 96 u(select_u & overlapping) = ... 97 min (result.sup(select_u & overlapping), ... 98 arccosineshifted.sup(select_u & overlapping)); 99 m = n; 100 m(select_u & ~overlapping) = ... 101 mpfr_function_d ('minus', +inf, m(select_u & ~overlapping), 1); 102 idx.subs = {(select_u & ~overlapping & rem (n, 2) == 0)}; 103 u(idx.subs {1}) = ... 104 sup (subsref (n, idx) .* pi - subsref (arccosine, idx)); 105 idx.subs = {(select_u & ~overlapping & rem (n, 2) ~= 0)}; 106 u(idx.subs {1}) = ... 107 sup (subsref (arccosine, idx) + subsref (m, idx) .* pi); 108 109 ## Find a larger lower bound for x, if the restriction from c allows it 110 l = -inf (size (result.inf)); 111 select_l = select & result.inf > -inf; 112 ## Find n, such that result.inf is within a distance of pi/2 113 ## around (n + 1/2) * pi. 114 n = result.inf; 115 n(select_l) = floor (inf (n(select_l) ./ pi)); 116 arccosineshifted = arccosine; 117 idx.subs = {(select_l & rem (n, 2) == 0)}; 118 arccosineshifted = subsasgn (arccosineshifted, idx, ... 119 subsref (arccosine, idx) + ... 120 subsref (n, idx) .* pi); 121 idx.subs = {(select_l & rem (n, 2) ~= 0)}; 122 arccosineshifted = subsasgn (arccosineshifted, idx, ... 123 (infsup (subsref (n, idx)) + 1) .* pi - ... 124 subsref (arccosine, idx)); 125 overlapping = not (isempty (intersect (result, arccosineshifted))); 126 l(select_l & overlapping) = ... 127 max (result.inf(select_l & overlapping), ... 128 arccosineshifted.inf(select_l & overlapping)); 129 m = n; 130 m(select_l & ~overlapping) = ... 131 mpfr_function_d ('plus', -inf, m(select_l & ~overlapping), 1); 132 idx.subs = {(select_l & ~overlapping & rem (n, 2) == 0)}; 133 l(idx.subs {1}) = ... 134 inf ((infsup (subsref (m, idx)) + 1) .* pi - subsref (arccosine, idx)); 135 idx.subs = {(select_l & ~overlapping & rem (n, 2) ~= 0)}; 136 l(idx.subs {1}) = ... 137 inf (subsref (arccosine, idx) + subsref (m, idx) .* pi); 138 139 result.inf(select) = max (l(select), result.inf(select)); 140 result.sup(select) = min (u(select), result.sup(select)); 141 142 result.inf(result.inf > result.sup) = inf; 143 result.sup(result.inf > result.sup) = -inf; 144 endif 145 146endfunction 147 148%!# from the documentation string 149%!assert (cosrev (0, infsup (6, 9)) == "[0x1.F6A7A2955385Ep2, 0x1.F6A7A2955386p2]"); 150 151%!# correct use of signed zeros 152%!test 153%! x = cosrev (infsup (1), infsup (0)); 154%! assert (signbit (inf (x))); 155%! assert (not (signbit (sup (x)))); 156 157%!shared testdata 158%! # Load compiled test data (from src/test/*.itl) 159%! testdata = load (file_in_loadpath ("test/itl.mat")); 160 161%!test 162%! # Scalar evaluation 163%! testcases = testdata.NoSignal.infsup.cosRev; 164%! for testcase = [testcases]' 165%! assert (isequaln (... 166%! cosrev (testcase.in{1}), ... 167%! testcase.out)); 168%! endfor 169 170%!test 171%! # Vector evaluation 172%! testcases = testdata.NoSignal.infsup.cosRev; 173%! in1 = vertcat (vertcat (testcases.in){:, 1}); 174%! out = vertcat (testcases.out); 175%! assert (isequaln (cosrev (in1), out)); 176 177%!test 178%! # N-dimensional array evaluation 179%! testcases = testdata.NoSignal.infsup.cosRev; 180%! in1 = vertcat (vertcat (testcases.in){:, 1}); 181%! out = vertcat (testcases.out); 182%! # Reshape data 183%! i = -1; 184%! do 185%! i = i + 1; 186%! testsize = factor (numel (in1) + i); 187%! until (numel (testsize) > 2) 188%! in1 = reshape ([in1; in1(1:i)], testsize); 189%! out = reshape ([out; out(1:i)], testsize); 190%! assert (isequaln (cosrev (in1), out)); 191 192%!test 193%! # Scalar evaluation 194%! testcases = testdata.NoSignal.infsup.cosRevBin; 195%! for testcase = [testcases]' 196%! assert (isequaln (... 197%! cosrev (testcase.in{1}, testcase.in{2}), ... 198%! testcase.out)); 199%! endfor 200 201%!test 202%! # Vector evaluation 203%! testcases = testdata.NoSignal.infsup.cosRevBin; 204%! in1 = vertcat (vertcat (testcases.in){:, 1}); 205%! in2 = vertcat (vertcat (testcases.in){:, 2}); 206%! out = vertcat (testcases.out); 207%! assert (isequaln (cosrev (in1, in2), out)); 208 209%!test 210%! # N-dimensional array evaluation 211%! testcases = testdata.NoSignal.infsup.cosRevBin; 212%! in1 = vertcat (vertcat (testcases.in){:, 1}); 213%! in2 = vertcat (vertcat (testcases.in){:, 2}); 214%! out = vertcat (testcases.out); 215%! # Reshape data 216%! i = -1; 217%! do 218%! i = i + 1; 219%! testsize = factor (numel (in1) + i); 220%! until (numel (testsize) > 2) 221%! in1 = reshape ([in1; in1(1:i)], testsize); 222%! in2 = reshape ([in2; in2(1:i)], testsize); 223%! out = reshape ([out; out(1:i)], testsize); 224%! assert (isequaln (cosrev (in1, in2), out)); 225