sym() and vpa() refactor

NEWS

@@ -26,6 +26,12 @@ octsympy 2.5.1-dev
     using `sym(0.1, 'f')`.  Passing `'r'` gives the current heuristic
     behaviour.  Neither raises a warning.
 
+  * The default floating-point heuristics now check for square roots of
+    integers, so `sym(sqrt(2))` should work.
+
+  * `sym('i')` and `sym('I')` no longer create the imaginary unit: use
+    `sym(i)` or `sym(1i)`.
+
   * Bug fix: isequal(n) for symfun's correctly tests the argnames.
 
 

inst/@sym/private/detect_special_str.m

@@ -0,0 +1,60 @@
+%% Copyright (C) 2015-2017 Colin B. Macdonald
+%% Copyright (C) 2016 Lagu
+%%
+%% This file is part of OctSymPy.
+%%
+%% OctSymPy is free software; you can redistribute it and/or modify
+%% it under the terms of the GNU General Public License as published
+%% by the Free Software Foundation; either version 3 of the License,
+%% or (at your option) any later version.
+%%
+%% This software is distributed in the hope that it will be useful,
+%% but WITHOUT ANY WARRANTY; without even the implied warranty
+%% of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See
+%% the GNU General Public License for more details.
+%%
+%% You should have received a copy of the GNU General Public
+%% License along with this software; see the file COPYING.
+%% If not, see <http://www.gnu.org/licenses/>.
+
+%% -*- texinfo -*-
+%% @deftypefun {@var{s} =} detect_special_str (@var{x})
+%% Recognize special constants, return their sympy string equivalent.
+%%
+%% Private helper function.
+%%
+%% This function should return a string @var{s} which can be instantiated
+%% directly in Python such as @code{['return' @var{s}]}.  If no special
+%% value is found, it returns the empty string.
+%%
+%% @seealso{sym, vpa}
+%% @end deftypefun
+
+function s = detect_special_str (x)
+
+  % Table of special strings.  Format for each row is:
+  %    {list of strings to recognize}, resulting python expr
+  % Note: case sensitive
+  % Note: python expr should be in list for identity "sym(sympy(x)) == x"
+  table = {{'pi'} 'S.Pi'; ...
+           {'inf' 'Inf' 'oo'} 'S.Infinity'; ...
+           {'NaN' 'nan'} 'S.NaN'; ...
+           {'zoo'} 'S.ComplexInfinity'};
+
+  s = '';
+
+  assert (ischar (x))
+
+  for j = 1:length (table)
+    for n = 1:length (table{j, 1})
+      if (strcmp (x, table{j, 1}{n}) || strcmp (x, ['+' table{j, 1}{n}]))
+        s = table{j, 2};
+        return
+      elseif (strcmp (x, ['-' table{j, 1}{n}]))
+        s = ['-' table{j, 2}];
+        return
+      end
+    end
+  end
+
+end

inst/private/magic_str_str.m → inst/@sym/private/double_to_sym_exact.m

@@ -1,4 +1,4 @@
-%% Copyright (C) 2015 Colin B. Macdonald
+%% Copyright (C) 2017 Colin B. Macdonald
 %%
 %% This file is part of OctSymPy.
 %%
@@ -17,29 +17,22 @@
 %% If not, see <http://www.gnu.org/licenses/>.
 
 %% -*- texinfo -*-
-%% @deftypefun {@var{s} =} magic_str_str (@var{x})
-%% Recognize special string values and substitute others.
+%% @deftypefun {@var{y} =} double_to_sym_exact (@var{x})
+%% Convert a double value to the equivalent rational sym
 %%
-%% Private function.
+%% Private helper function.
 %%
-%% @seealso{sym, vpa}
 %% @end deftypefun
 
-function s = magic_str_str(x, varargin)
-
-  if (~ischar(x))
-    error('OctSymPy:magic_str_str:notstring', ...
-          'Expected a string');
-  end
-
-  if (strcmpi(x, 'inf')) || (strcmpi(x, '+inf'))
-    s = 'oo';
-  elseif (strcmpi(x, '-inf'))
-    s = '-oo';
-  elseif (strcmpi(x, 'i'))
-    s = 'I';
+function y = double_to_sym_exact (x)
+  if (isnan (x))
+    y = python_cmd ('return S.NaN');
+  elseif (isinf (x) && x < 0)
+    y = python_cmd ('return -S.Infinity');
+  elseif (isinf (x))
+    y = python_cmd ('return S.Infinity');
   else
-    s = x;
+    %% Rational will exactly convert from a float
+    y = python_cmd ('return Rational(_ins[0])', x);
   end
-
 end

inst/@sym/private/double_to_sym_heuristic.m

@@ -0,0 +1,64 @@
+%% Copyright (C) 2017 Colin B. Macdonald
+%%
+%% This file is part of OctSymPy.
+%%
+%% OctSymPy is free software; you can redistribute it and/or modify
+%% it under the terms of the GNU General Public License as published
+%% by the Free Software Foundation; either version 3 of the License,
+%% or (at your option) any later version.
+%%
+%% This software is distributed in the hope that it will be useful,
+%% but WITHOUT ANY WARRANTY; without even the implied warranty
+%% of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See
+%% the GNU General Public License for more details.
+%%
+%% You should have received a copy of the GNU General Public
+%% License along with this software; see the file COPYING.
+%% If not, see <http://www.gnu.org/licenses/>.
+
+%% -*- texinfo -*-
+%% @deftypefun {@var{y} =} double_to_sym_heuristic (@var{x}, @var{ratwarn}, @var{argnstr})
+%% Convert a double value to a nearby "nice" sym
+%%
+%% Private helper function.
+%%
+%% @end deftypefun
+
+function y = double_to_sym_heuristic (x, ratwarn, argnstr)
+  assert (isempty (argnstr))
+  if (isnan (x))
+    y = python_cmd ('return S.NaN');
+  elseif (isinf (x) && x < 0)
+    y = python_cmd ('return -S.Infinity');
+  elseif (isinf (x))
+    y = python_cmd ('return S.Infinity');
+  elseif (isequal (x, pi))
+    %% Special case for pi
+    y = python_cmd ('return S.Pi');
+  elseif (isequal (x, -pi))
+    y = python_cmd ('return -S.Pi');
+  elseif ((abs (x) < flintmax) && (mod (x, 1) == 0))
+    y = python_cmd ('return Integer(_ins[0])', int64 (x));
+  else
+    %% Find the nearest rational, rational*pi or sqrt(integer)
+    if (ratwarn)
+      warning('OctSymPy:sym:rationalapprox', ...
+              'passing floating-point values to sym is dangerous, see "help sym"');
+    end
+    [N1, D1] = rat (x);
+    [N2, D2] = rat (x / pi);
+    N3 = round (x^2);
+    err1 = abs (N1 / D1 - x);
+    err2 = abs ((N2*pi) / D2 - x);
+    err3 = abs (sqrt (N3) - x);
+    if (err1 <= err3)
+      if (err1 <= err2)
+        y = python_cmd ('return Rational(*_ins)', int64 (N1), int64 (D1));
+      else
+        y = python_cmd ('return Rational(*_ins)*S.Pi', int64 (N2), int64 (D2));
+      end
+    else
+      y = python_cmd ('return sqrt(Integer(*_ins))', int64 (N3));
+    end
+  end
+end