1%% Copyright (C) 2014-2017, 2019 Colin B. Macdonald
2%%
3%% This file is part of OctSymPy.
4%%
5%% OctSymPy is free software; you can redistribute it and/or modify
6%% it under the terms of the GNU General Public License as published
7%% by the Free Software Foundation; either version 3 of the License,
8%% or (at your option) any later version.
9%%
10%% This software is distributed in the hope that it will be useful,
11%% but WITHOUT ANY WARRANTY; without even the implied warranty
12%% of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See
13%% the GNU General Public License for more details.
14%%
15%% You should have received a copy of the GNU General Public
16%% License along with this software; see the file COPYING.
17%% If not, see <http://www.gnu.org/licenses/>.
18
19%% -*- texinfo -*-
20%% @documentencoding UTF-8
21%% @defmethod @@sym double (@var{x})
22%% Convert symbolic to doubles.
23%%
24%% Example:
25%% @example
26%% @group
27%% x = sym(1) / 3
28%%   @result{} x = (sym) 1/3
29%% @c doctest: +SKIP_IF(compare_versions (OCTAVE_VERSION(), '6.0.0', '<'))
30%% double (x)
31%%   @result{} ans =  0.3333
32%% @end group
33%% @end example
34%%
35%% Despite the name, this is one way to convert a complex sym to
36%% floating point:
37%% @example
38%% @group
39%% z = sym(4i) - 3;
40%% double (z)
41%%   @result{} ans = -3 + 4i
42%% @end group
43%% @end example
44%%
45%% If conversion fails, you get an error:
46%% @example
47%% @group
48%% syms x
49%% double (x)
50%%   @print{} ??? ... can't convert expression ...
51%% @end group
52%% @end example
53%%
54%% @seealso{sym, vpa}
55%% @end defmethod
56
57
58function y = double(x)
59
60  % FIXME: port to uniop?
61
62  if ~(isscalar(x))
63    % sympy N() works fine on matrices but it gives objects like "Matrix([[1.0,2.0]])"
64    y = zeros(size(x));
65    for j = 1:numel(x)
66      % temp = x(j)  (Issue #17)
67      idx.type = '()';
68      idx.subs = {j};
69      temp = double(subsref(x,idx));
70      if (isempty(temp))
71        y = [];
72        return
73      end
74      y(j) = temp;
75    end
76    return
77  end
78
79  cmd = { '(x,) = _ins'
80          'if x == zoo:'  % zoo -> Inf + Infi
81          '    return (float(sp.oo), float(sp.oo))'
82          'if x == nan:'
83          '    return (float(nan), 0.0)'
84          'x = complex(x)'
85          'return (x.real, x.imag)'
86        };
87
88  [A, B] = pycall_sympy__ (cmd, x);
89
90  %y = A + B*i;  % not quite the same for Inf + InFi
91  if (B == 0.0)
92    y = A;
93  else
94    y = complex(A, B);
95  end
96
97end
98
99
100%!test
101%! % numeric scalar
102%! a = double(sym(10));
103%! assert (a == 10)
104%! assert (isa (a, 'double'))
105
106%!test
107%! % numeric vectors
108%! a = double(sym([10 12]));
109%! assert (isequal (a, [10 12]))
110%! assert (isa (a, 'double'))
111
112%!test
113%! % complex
114%! a = 3 + 4i;
115%! b = sym(a);
116%! assert (isequal (double (b), a))
117
118%!xtest
119%! % unexpected, precisely same floating point
120%! a = 3 + 4i;
121%! b = sym(a);
122%! assert (isequal (double (b/pi), a/pi))
123
124%!test
125%! % floating point
126%! x = sqrt(sym(2));
127%! assert( abs(double(x) - sqrt(2)) < 2*eps)
128%! x = sym(pi);
129%! assert( abs(double(x) - pi) < 2*eps)
130
131%!test
132%! oo = sym(inf);
133%! assert( double(oo) == inf )
134%! assert( double(-oo) == -inf )
135%! assert( isnan(double(0*oo)) )
136
137%!test
138%! zoo = sym('zoo');
139%! assert (double(zoo) == complex(inf, inf))
140
141%!test
142%! zoo = sym('zoo');
143%! assert (double(-zoo) == double(zoo) )
144%! assert( isnan(double(0*zoo)) )
145
146%!test
147%! % nan
148%! snan = sym(nan);
149%! assert( isnan(double(snan)))
150
151%!test
152%! % don't want NaN+NaNi
153%! snan = sym(nan);
154%! assert (isreal (double (snan)))
155
156%!test
157%! % arrays
158%! a = [1 2; 3 4];
159%! assert( isequal(  double(sym(a)), a  ))
160%! assert( isequal(  double(sym(a)), a  ))
161
162%! % should fail with error for non-double
163%!error <TypeError> syms x; double(x)
164%!error <TypeError> syms x; double([1 2 x])
165