% usally can be treated as 1x1 class 'sym'
classdef Expr < sympy.Base
    

    methods (Access=public)
        % ===================== 构造函数（通用表达式构造） =====================
        function obj = Expr(x)
            obj@sympy.Base(x);
        end


        % ===================== 一元运算符重载 =====================
        function res = uminus(a)
            % 重载一元负号：-a
            res = sympy.Expr(-getObj(a));
        end

        function res = uplus(a)
            res = sympy.Expr(+getObj(a));
        end

        % ===================== 二元算术运算符重载 =====================
        function res = plus(a, b)
            % 重载加法：a + b
            a_sp = sympy.to_sp_object(a);
            b_sp = sympy.to_sp_object(b);
            res = sympy.Expr((a_sp + b_sp).simplify());
        end

        function res = minus(a, b)
            a_sp = sympy.to_sp_object(a);
            b_sp = sympy.to_sp_object(b);
            res = sympy.Expr((a_sp - b_sp).simplify());
        end

        function res = mtimes(a, b)
            % 重载矩阵乘法：a * b（符号运算中等同于普通乘法）
            a_sp = sympy.to_sp_object(a);
            b_sp = sympy.to_sp_object(b);
            res = sympy.Expr(a_sp * b_sp);
        end

        function res = times(a, b)
            % 重载元素乘法：a .* b（与矩阵乘法逻辑一致）
            a_sp = sympy.to_sp_object(a);
            b_sp = sympy.to_sp_object(b);
            res = sympy.Expr(a_sp * b_sp);
        end

        function res = mrdivide(a, b)
            % 重载矩阵右除：a / b
            a_sp = sympy.to_sp_object(a);
            b_sp = sympy.to_sp_object(b);
            res = sympy.Expr(a_sp / b_sp);
        end

        function res = ldivide(a, b)
            % 重载矩阵左除：a \ b (等价于 b / a)
            a_sp = sympy.to_sp_object(a);
            b_sp = sympy.to_sp_object(b);
            res = sympy.Expr(b_sp / a_sp);
        end

        function res = mpower(a, b)
            a_sp = sympy.to_sp_object(a);
            b_sp = sympy.to_sp_object(b);
            res = sympy.Expr(a_sp^b_sp);
        end

        % ===================== 比较运算符重载 =====================
        function res = eq(a, b)
            % 重载等于：a == b
            a_sp = sympy.to_sp_object(a);
            b_sp = sympy.to_sp_object(b);
            sp = sympy__sp;
            res = sympy.Eq(sp.Eq(a_sp, b_sp));
        end

        % ===================== 转置/共轭转置重载 =====================
        function res = transpose(a)
            % 重载转置：a.'
            sp = sympy__sp;
            res = sympy.Expr(sp.transpose(getObj(a)));
        end

        function res = ctranspose(a)
            % 重载共轭转置：a'
            sp = sympy__sp;
            res = sympy.Expr(sp.conjugate(sp.transpose(getObj(a))));
        end

        function res = vertcat(varargin)
            sp = sympy__sp;
            n = nargin;
            for i = 1:n
                varargin{i} = sympy.to_sp_object(varargin{i}); % todo: support double matrix
            end
            n = nargin;
            res = sympy.Matrix(n, 1, sp.Matrix(n, 1, pylist(varargin{:})));
        end

        function res = horzcat(varargin)
            sp = sympy__sp;
            for i = 1:nargin
                varargin{i} = sympy.to_sp_object(varargin{i}); % todo: support double matrix
            end
            n = nargin;
            res = sympy.Matrix(1, n, sp.Matrix(1, n, pylist(varargin{:})));
        end

        % ====================== 通用函数 =====================
        function f = double(x)
            % 转换为数值类型
            py = pyroot;
            f = py.float(getObj(x).evalf());  % 修复原代码x.evalf()的错误
        end

        function r = vpa(x, d)
            % 数值近似计算
            narginchk(1, 2)
            if nargin == 1; d = 32; end
            sp = sympy__sp;
            py = pyroot;
            r = sympy.Expr(sp.N(getObj(x), d));
        end


        function df = diff(f, x, varargin)
            % 表达式求偏导。symdiff(f, x, varargin) 相当于 f.diff(x, varargin{:})
            % df = symdiff(f, x) 表示f对x的偏导数；
            % df = symdiff(f, x, order) 表示f对x的order阶偏导数，order应该输入整数；
            % df = symdiff(f, x, ...) 后面输入更多的符号，表示混合偏导数（请自行注意输入顺序）。
            % 例子：
            % Input [1] >> [a, x] = syms('a x'); 
            % Input [2] >> symdiff(exp(a*x), x)
            % a*exp(a*x)
            % Input [3] >> symdiff(exp(a*x), x, 2)
            % a**2*exp(a*x)
            % Input [4] >> symdiff(exp(a*x), x, a)
            % (a*x + 1)*exp(a*x)
            if nargin == 1; x = symvar(f, 1); end
            sp = sympy__sp();
            if nargin == 2; 
                df = sp.diff(getObj(f), getObj(x));
            else
                varargin = cellfun(@(o)sympy.to_sp_object(o), varargin, 'UniformOutput', false);
                df = sp.diff(getObj(f), getObj(x), varargin{:});
            end
            df = f.setObj(df);
        end

        function g = gradient(f, vars)
            sp = sympy__sp;
            % 计算表达式f关于所有自由变量的梯度向量，模仿MATLAB gradient
            % 返回结果为与自变量数相等的向量，每个分量为对该变量的偏导
            if nargin < 2
                vars = symvar(f);  % 获取所有符号变量
                if isempty(vars)
                    error("无法对常数或无变量的表达式求梯度");
                end
                n = numel(vars);
            else
                if isa(vars, 'sympy.Matrix')
                    mat_vars = getObj(vars);  % 获取底层的Python矩阵对象
                    n = numel(vars);  % 获取矩阵维度
                    linear_vars = cell(1, n);
                    py = pyroot; isinstance = py.isinstance;
                    for i = 1:n
                        t = mat_vars{i-1};
                        if ~isinstance(t, sp.Symbol); error("期望输入Symbol数组"); end
                        linear_vars{i} = t;  % Py对象矩阵索引从0计；cell下标从1计
                    end
                    vars = linear_vars;
                elseif isa(vars, 'sympy.Expr')
                    g = diff(f, vars);
                    return;
                else
                    error("期望输入Symbol或Symbol数组");
                end
            end
            f = getObj(f);
            grad_list = pylist;
            for k = 1:n
                grad_list.append(f.diff(vars{k}));
            end
            % 将结果组装为sympy.Matrix形式的“梯度列向量”
            gradient_vector = sp.Matrix(n, 1, grad_list);
            g = sympy.Matrix(n, 1, gradient_vector);
        end

        function v = limit(f, var, a, dir)
            % 求表达式的极限
            % v = limit(f, var, a) 求表达式f中的变量var在a点的（左右）极限
            % v = limit(f, var, a, dir) 求表达式f中的变量var在a点的方向极限，dir可以是"+","-","+-"。
            if nargin < 3; error("输入参数不够"); end;
            if nargin < 4; dir = "+-"; end
            if nargin < 3; a = 0; end
            if any(strcmp(["left", "-"], dir))
                dir = "-";
            elseif any(strcmp(["right", "+"], dir))
                dir = "+";
            else
                dir = "+-";
            end
            sp = sympy__sp;
            expr = getObj(f);
            var = getObj(var);
            v = sp.limit(expr, var, a, dir).simplify();
            v = f.setObj(v);
        end 

        function r = int(f, varargin)
            %SYMINT - 符号积分，支持不定积分、定积分.
            %
            % 语法：
            %   r = int(f);                      - 对 symvar 得到的第一个符号变量求不定积分
            %   r = int(f, var);                 - 对符号变量 var 求不定积分
            %   r = int(f, a, b);                - 对 symvar 得到的第一个符号变量求定积分
            %   r = int(f, var, a, b, varargin); - 支持的键值对包括 "Hold"、"PrincipalValue" 和 "IgnoreSpecialCases"
            %       - "IgnoreAnalyticConstraints": 默认为 false, 忽略解析限制, 功能不完全, 不建议使用;
            %       - "IgnoreSpecialCases": 默认为 false, 当 "IgnoreSpecialCases" 为 true 时, 忽略积分的特殊情况;
            %       - "PrincipalValue": 默认为 false, 当 "PrincipalValue" 为 true 时, 计算广义积分的 Cauchy 主值;
            %       - "Hold": 默认为 false, 当 "Hold" 为 true 时, 保持积分表达式.
            %
            % 示例:
            % 不定积分
            % >> syms x y
            % >> int(cos(x))
            % sin(x)
            % >> int(log(x + y), x)
            % x*log(x + y) - x + y*log(x + y)
            %
            % 定积分、重积分
            % >> int(exp(-x), 0, inf)
            % 1
            % >> int(2*x, sin(y), cos(y))
            % cos(2*y)
            % >> int(int(exp(-x^2 - y^2), x, -inf, inf), y, -inf, inf)
            % pi
            %
            % 矩阵元素的积分和忽略特殊情况
            % >> M = sym(zeros(2, 2));
            % >> M{0, 0} = exp(y);
            % >> M{0, 1} = exp(x * y);
            % >> M{1, 0} = sin(y);
            % >> M{1, 1} = cos(x);
            % >> int(M, y)
            % Matrix([[exp(y), Piecewise((exp(x*y)/x, (x > 0) | (x < 0)), (y, True))], [-cos(y), y*cos(x)]])
            % >> int(M, y, "IgnoreSpecialCases", true)
            % Matrix([[exp(y), exp(x*y)/x], [-cos(y), y*cos(x)]])
            %
            % 保持积分号和计算 Cauchy 主值
            % >> int(1/(x - 1), x, 0, 2)
            % nan
            % >> int(1/(x - 1), x, 0, 2, "Hold", true)
            % Integral(1/(x - 1), (x, 0, 2))
            % >> int(1/(x - 1), x, 0, 2, "PrincipalValue", true)
            % 0
            %
            % 忽略解析限制
            % >> int(asin(sin(x)), x, 0, 5)
            % Integral(asin(sin(x)), (x, 0, 5))
            % >> int(asin(sin(x)), x, 0, 5, "IgnoreAnalyticConstraints", true)
            % 25/2
            %
            py = pyroot;
            sp = sympy__sp;
            defaultVal = struct("var", [], "a", [], "b", [], "IgnoreAnalyticConstraints", false, "IgnoreSpecialCases", false, "PrincipalValue", false, "Hold", false);
            Results = sympy.parseParameters(3, defaultVal, @sympy.to_sp_object, varargin);

            expr = getObj(f);
            var = Results.var;   % optional #1
            a = Results.a;       % optional #2
            b = Results.b;       % optional #3

            hasVar = ~isempty(var);
            hasA = ~isempty(a);
            hasB = ~isempty(b);

            if ~hasVar && ~hasA && ~hasB
                % int(expr): use first free symbol as integration variable
                defvar = getObj(symvar(f, 1));
                r = sp.Integral(expr, defvar);
            elseif hasVar && ~hasA && ~hasB
                % int(expr, var)
                r = sp.Integral(expr, var);
            elseif hasVar && hasA && ~hasB
                % int(expr, a, b): var/a store lower/upper bounds here
                defvar = getObj(symvar(f, 1));
                r = sp.Integral(expr, pytuple(defvar, var, a));
            elseif hasVar && hasA && hasB
                % int(expr, var, a, b)
                r = sp.Integral(expr, pytuple(var, a, b));
            else
                error("积分参数格式错误");
            end
            
            if Results.Hold
                r = r;
            elseif Results.PrincipalValue
                r = r.principal_value();
            elseif Results.IgnoreSpecialCases
                r = r.doit(pyargs("conds", "none"));
                r = r.simplify();
            else
                r = r.doit();
                r = r.simplify();
            end
            r = f.setObj(r);
        end

        function s = vpaintegral(expr, varargin)
            narginchk(2, inf);
            sp = sympy__sp;
            exprSp = getObj(expr);

            args = varargin;
            firstNv = numel(args) + 1;
            for i = 1:numel(args)
                if ischar(args{i}) || isstring(args{i})
                    firstNv = i;
                    break;
                end
            end

            posArgs = args(1:firstNv - 1);
            nvArgs = args(firstNv:end);
            if mod(numel(nvArgs), 2) ~= 0
                error("Name-Value 参数必须成对出现");
            end

            waypoints = [];
            for i = 1:2:numel(nvArgs)
                key = lower(string(nvArgs{i}));
                val = nvArgs{i + 1};
                switch key
                    case "waypoints"
                        waypoints = val;
                    otherwise
                        error("不支持的参数名: %s", key);
                end
            end

            if isempty(posArgs)
                error("vpaintegral 缺少积分区间");
            end

            if isa(posArgs{1}, "sympy.Base")
                var = sympy.to_sp_object(posArgs{1});
                boundArgs = posArgs(2:end);
            else
                defvar = symvar(expr, 1);
                if isempty(defvar)
                    error("表达式中没有可用于积分的变量");
                end
                var = getObj(defvar);
                boundArgs = posArgs;
            end

            if numel(boundArgs) == 1
                lim = boundArgs{1};
                if numel(lim) ~= 2
                    error("积分界限错误");
                end
                a = sympy.Expr.toBound(lim(1), sp);
                b = sympy.Expr.toBound(lim(2), sp);
            elseif numel(boundArgs) == 2
                a = sympy.Expr.toBound(boundArgs{1}, sp);
                b = sympy.Expr.toBound(boundArgs{2}, sp);
            else
                error("积分参数格式错误");
            end

            if isempty(waypoints)
                sSp = sp.N(sp.Integral(exprSp, pytuple(var, a, b)));
            else
                if ~isvector(waypoints)
                    error("Waypoints 必须是向量");
                end
                points = cell(1, numel(waypoints) + 2);
                points{1} = a;
                for i = 1:numel(waypoints)
                    points{i + 1} = sympy.Expr.toScalarSp(waypoints(i), sp);
                end
                points{end} = b;

                t = sp.Symbol("_t");
                sSp = sp.Integer(0);
                for i = 1:numel(points) - 1
                    z0 = points{i};
                    z1 = points{i + 1};
                    dz = z1 - z0;
                    zt = z0 + dz * t;
                    segExpr = exprSp.subs(var, zt) * dz;
                    sSp = sSp + sp.N(sp.Integral(segExpr, pytuple(t, 0, 1)));
                end
            end

            s = expr.setObj(sp.simplify(sSp));
        end

        function ret = changeIntegrationVariable(expr, old, new)
            exprSp = getObj(expr);
            oldSp = sympy.to_sp_object(old);
            newSp = sympy.to_sp_object(new);
            py = pyroot;

            % Integral uses SymPy transform; generic expressions use substitution.
            isIntegral = false;
            try
                isIntegral = logical(py.bool(exprSp.is_Integral));
            catch
            end

            if isIntegral
                retSp = exprSp.transform(oldSp, newSp);
            else
                retSp = exprSp.subs(oldSp, newSp);
            end

            ret = expr.setObj(retSp);
        end

        function ret = integrateByParts(expr, dv)
            % 分部积分：Integral(u*dv) -> u*v - Integral(v*du)
            py = pyroot;
            sp = sympy__sp;
            exprSp = getObj(expr);
            dvSp = sympy.to_sp_object(dv);

            isIntegral = false;
            try
                isIntegral = logical(py.isinstance(exprSp, sp.Integral));
            catch
            end
            if ~isIntegral
                try
                    isIntegral = logical(py.bool(exprSp.is_Integral));
                catch
                end
            end

            if isIntegral
                % 单重积分路径：兼容 limits / variables 两种属性风格。
                funcSp = exprSp.function;
                uSp = funcSp / dvSp;
                var = [];
                lower = [];
                upper = [];
                lenLimit = int64(0);
                hasLimits = false;

                try
                    limits = exprSp.limits;
                    if py.len(limits) == 1
                        hasLimits = true;
                        limit = limits{0};
                        lenLimit = int64(py.len(limit));
                        var = limit{0};
                        if lenLimit == 3
                            lower = limit{1};
                            upper = limit{2};
                        end
                    end
                catch
                end

                if ~hasLimits
                    try
                        vars = exprSp.variables;
                        if py.len(vars) == 1
                            var = vars{0};
                            lenLimit = int64(1);
                        else
                            error("integrateByParts 仅支持单重积分");
                        end
                    catch
                        error("integrateByParts 仅支持单重积分");
                    end
                end

                vSp = sp.integrate(dvSp, var);
                duSp = sp.diff(uSp, var);

                if lenLimit == 1
                    retSp = uSp * vSp - sp.Integral(vSp * duSp, var);
                elseif lenLimit == 3
                    uv = uSp * vSp;
                    boundary = uv.subs(var, upper) - uv.subs(var, lower);
                    retSp = boundary - sp.Integral(vSp * duSp, pytuple(var, lower, upper));
                else
                    error("积分上下限格式错误");
                end
            else
                % 回退路径：当 expr 不是 Integral 时，视作被积函数继续做分部积分。
                var = [];
                try
                    dvVars = py.list(dvSp.free_symbols);
                    if py.len(dvVars) > 0
                        var = dvVars{0};
                    end
                catch
                end
                if isempty(var)
                    vars = py.list(exprSp.free_symbols);
                    if py.len(vars) == 0
                        error("integrateByParts 无法从表达式中推断积分变量");
                    end
                    var = vars{0};
                end
                vSp = sp.integrate(dvSp, var);
                uSp = exprSp / dvSp;
                duSp = sp.diff(uSp, var);
                retSp = uSp * vSp - sp.Integral(vSp * duSp, var);
            end

            ret = expr.setObj(retSp);
        end

        function vars = symvar(s, n)
            % SYMVAR 提取符号表达式/方程中的自由符号变量（适配自定义类）
            % 语法（与MATLAB原生一致）：
            %   vars = symvar(s)        - 返回s中的所有自由符号变量，按字母顺序排序
            %   vars = symvar(s, n)     - 返回前n个符号变量（按字母顺序）
            % 输入：
            %   s - 自定义符号类对象：sympy.Expr/sympy.Eq/sympy.Function（单个），或其元胞数组
            %   n - 可选，正整数，指定返回的变量数量（默认返回所有）
            % 输出：
            %   vars - 符号变量数组（sympy.Expr类对象），按字母顺序排序
            % 核心：通过getObj()获取sympy原生对象，提取free_symbols并处理
            
            % ===================== 步骤1：参数校验与初始化 =====================
            narginchk(1, 2);  % 限制输入参数数量为1或2
            py = pyroot;
            sp = sympy__sp;

            % 默认返回所有符号变量
            if nargin == 1
                n = inf;
            else
                if ~isnumeric(n) || n < 1 || mod(n, 1) ~= 0
                    error('第二个参数必须是正整数！');
                end
            end
            
            s_sp = s.getObj();  % 核心：通过getObj获取sympy原生对象
            % 提取自由符号并加入集合
            syms_set = s_sp.free_symbols;

            % ===================== 步骤3：处理无符号变量的情况 =====================
            if py.len(syms_set) == 0
                vars = [];  % 无符号变量返回空数组
                return;
            end

            % ===================== 步骤4：转换为列表并按字母顺序排序（匹配MATLAB行为） =====================
            % 转换为Python列表并排序（按符号名字母顺序）
            syms_list = py.list(syms_set);

            % ===================== 步骤5：限制返回数量并封装为自定义Expr对象 =====================
            % 确定返回的符号数量
            num_syms = min(py.len(syms_list), n);
            
            if num_syms == 1
                vars = sympy.Symbol(syms_list{0});
                return;
            end
            
            vars = cell(1, num_syms);
            % return sympy.Matrix 1xn
            vars = sympy.Matrix(vars);
        end

        function nexpr = collect(expr, v)
            % nexpr = collect(expr)
            % nexpr = collect(expr,v)
            
            % nexpr = collect(expr) 对符号表达式 expr 中相同幂次的项按照默认变量进行系数合并
            % nexpr = collect(expr,v) 对符号表达式expr中相同幂次的项按照指定变量v的多项式形式返回
            
            %% 例子：
            % 按照默认变量进行合并：
            % Input [1] >> x = sym('x');
            % Input [2] >> collect((x+3)*(x+2))
            % x**2 + 5*x + 6
            %% 按照指定变量进行合并：
            % Input [1] >> [x,y] = syms('x y');
            % Input [2] >> collect(x^2*y + 3*x*y - 2*x^2 - 5*x, x)
            % 输出：x**2*(y - 2) + x*(3*y - 5)
            
            py = pyroot;
            expr = getObj(expr).expand();
            if nargin == 1
                vars = py.list(expr.free_symbols);
                nexpr = expr.collect(vars{0});
            end
            
            if nargin ==2
                nexpr = expr.collect(v);
            end
            nexpr = sympy.Expr(expr);
        end

        function [F] = sum(expr, var, a, b)
            %SYMSUM - 符号级数求和.
            %
            % 语法:
            %   F = sum(expr, var, [a, b]);
            %   F = sum(expr, var, a, b);    - 计算级数 expr 关于求和索引 var 从下界 a 到上界 b 的和.
            %   F = sum(expr, var);          - 计算级数 expr 关于求和索引 var 的不定总和 (反差分).
            %
            % 示例:
            % >> sp = sympy__sp;
            % >> [k, n, x] = syms("k n x");
            % >> sum(k, k, 1, n)
            % n**2/2 + n/2
            % >> sum(1/k^2, k, [1, inf])
            % pi**2/6
            % >> sum(x^k/sp.factorial(k), k, 1, inf)
            % x*(exp(x)/x - 1/x)
            % >> sum(k, k)
            % k*(k - 1)/2
            % >> F = sum(k*x^k, k, 1, 8)
            % 8*x**8 + 7*x**7 + 6*x**6 + 5*x**5 + 4*x**4 + 3*x**3 + 2*x**2 + x
            % >> subs(F, 2)
            % 3586
            %
            narginchk(2, 4);
            sp = sympy__sp;
            expr = getObj(expr);
            if nargin == 4
                var = sympy.to_sp_object(var);
                a = sympy.to_sp_object(a);
                b = sympy.to_sp_object(b);
                F = sp.summation(expr, pytuple(var, a, b));
            elseif nargin == 3
                if numel(a) ~= 2
                    error("求和界限错误.");
                end
                [a, b] = deal(a(1), a(2));
                F = sp.summation(expr, pytuple(var, a, b));
            elseif nargin == 2
                py = pyroot;
                var = sympy.to_sp_object(var);
                f = sp.Function('f');
                n = sp.Symbol('n', pyargs('Integer', true));
                C0 = sp.Symbol('C0');
                rexpr = f(n + 1) - f(n) - expr.subs(var, n);
                F = sp.rsolve(rexpr, f(n));
                F = F.subs(n, var);
                F = F.subs(C0, 0);
            end
            F = sympy.Expr(F);
        end

        function s = taylor(expr, var, point, order)
            % 求解符号表达式的talor展开。
            % s = taylor(expr, var, point, order)
            % 求解表达式expr对于变量var在point处到order阶的taylor展开。var默认为0，order默认是5。
            sp = sympy__sp();
            if nargin < 4; order = 6; end
            if nargin < 3; point = 0; end
            expr = getObj(expr); var = getObj(var);
            s = sp.series(expr, var, point, order);
            try
                s = s.removeO();
            catch
            end
            s = sympy.Expr(sp.simplify(s));
        end

        function s = series(expr, varargin)
            %SERIES - 符号表达式级数展开。
            %
            % 语法：
            %   s = series(expr)
            %   s = series(expr, var)
            %   s = series(expr, var, 'Order', n)
            %   s = series(expr, var, 'ExpansionPoint', x0)
            %
            % 说明：
            %   默认展开变量为 symvar(expr, 1) 返回的首个符号；
            %   默认展开点为 0，默认展开阶数为 6。
            sp = sympy__sp();
            py = pyroot;

            var = [];
            point = 0;
            order = 6;
            args = varargin;

            % Optional: series(expr, var, ...)
            if ~isempty(args)
                firstArg = args{1};
                if ~(ischar(firstArg) || isstring(firstArg))
                    var = sympy.to_sp_object(firstArg);
                    args = args(2:end);
                end
            end

            if mod(numel(args), 2) ~= 0
                error("Name-Value 参数必须成对出现");
            end

            for i = 1:2:numel(args)
                key = string(args{i});
                val = args{i + 1};
                switch lower(key)
                    case "order"
                        order = val;
                    case "expansionpoint"
                        point = val;
                    otherwise
                        error("不支持的参数名: %s", key);
                end
            end

            if isempty(var)
                defaultVar = symvar(expr, 1);
                if isempty(defaultVar)
                    error("表达式中没有可用于展开的变量");
                end
                var = getObj(defaultVar);
            end

            if ~isnumeric(order) || ~isscalar(order) || order < 1 || mod(order, 1) ~= 0
                error("参数 'Order' 必须是正整数");
            end
            order = int32(order);

            if isnumeric(point) && isscalar(point) && isinf(point)
                if point > 0
                    point = sp.oo;
                else
                    point = -sp.oo;
                end
            else
                point = sympy.to_sp_object(point);
            end

            s = sp.series(getObj(expr), var, point, py.int(order));
            try
                s = s.removeO();
            catch
            end
            s = sympy.Expr(s);
        end

        function ft = fourier(expr, t, w)
            %FOURIER - 符号表达式傅里叶变换。
            %
            % 语法：
            %   ft = fourier(expr)
            %   ft = fourier(expr, t)
            %   ft = fourier(expr, t, w)
            %
            % 说明：
            %   默认变换变量为 symvar(expr, 1)；
            %   默认频域变量为符号 w。
            narginchk(1, 3);
            sp = sympy__sp;

            if nargin < 2 || isempty(t)
                defaultVar = symvar(expr, 1);
                if isempty(defaultVar)
                    t = sp.Symbol("x");
                else
                    t = getObj(defaultVar);
                end
            else
                t = sympy.to_sp_object(t);
            end

            if nargin < 3 || isempty(w)
                w = sp.Symbol("w");
            else
                w = sympy.to_sp_object(w);
            end

            ftSp = sp.fourier_transform(getObj(expr), t, w);
            ft = expr.setObj(ftSp);
        end

        function ft = ifourier(expr, varargin)
            %IFOURIER - 符号表达式逆傅里叶变换（与 SymPy inverse_fourier_transform 一致）。
            %
            % 语法：
            %   ft = ifourier(expr)
            %   ft = ifourier(expr, t)
            %   ft = ifourier(expr, w, t)
            %
            % 说明：
            %   单参数：频率变量优先取表达式中的 w（若存在），否则取 symvar(expr,1)；时间变量默认为 t。
            %   双参数：第二个参数为时间域变量；频率变量同上。
            %   三参数：显式指定频率变量 w 与时间域变量 t（对应 SymPy 的 (F, w, t)）。
            narginchk(1, 3);
            sp = sympy__sp;
            exprSp = getObj(expr);

            if nargin == 1
                wFreq = sympy.Expr.ifourierDefaultFreqVar(exprSp, sp);
                tTime = sp.Symbol("t");
            elseif nargin == 2
                tTime = sympy.to_sp_object(varargin{1});
                wFreq = sympy.Expr.ifourierDefaultFreqVar(exprSp, sp);
            else
                wFreq = sympy.to_sp_object(varargin{1});
                tTime = sympy.to_sp_object(varargin{2});
            end

            ftSp = sp.inverse_fourier_transform(exprSp, wFreq, tTime);
            ft = expr.setObj(ftSp);
        end

        function Ft = laplace(expr, varargin)
            %LAPLACE - 拉普拉斯变换（与 SymPy laplace_transform 一致）。
            %
            % 语法（勿使用与 MATLAB 工具箱同名的顶层 laplace()；请用本方法）：
            %   Ft = f.laplace()
            %   Ft = f.laplace(t)
            %   Ft = f.laplace(t, s)
            %
            % 说明：
            %   默认时间变量：若表达式含符号 t 则用 t，否则取 symvar 的首个符号；
            %   默认复频域变量为 s；可显式指定 t、s。
            narginchk(1, 3);
            sp = sympy__sp;
            py = pyroot;
            exprSp = getObj(expr);

            if nargin == 1
                tIn = sympy.Expr.laplaceDefaultTimeVar(exprSp, sp);
                sOut = sp.Symbol("s");
            elseif nargin == 2
                tIn = sympy.to_sp_object(varargin{1});
                sOut = sp.Symbol("s");
            else
                tIn = sympy.to_sp_object(varargin{1});
                sOut = sympy.to_sp_object(varargin{2});
            end

            raw = sp.laplace_transform(exprSp, tIn, sOut);
            if py.isinstance(raw, py.tuple)
                FtSp = raw{0};
            else
                FtSp = raw;
            end
            try
                FtSp = FtSp.simplify();
            catch
            end
            Ft = expr.setObj(FtSp);
        end

        function f = ilaplace(expr, varargin)
            %ILAPLACE - 逆拉普拉斯变换（与 SymPy inverse_laplace_transform 一致）。
            %
            % 语法（勿使用与 MATLAB 工具箱同名的顶层 ilaplace()；请用本方法）：
            %   g = F.ilaplace()
            %   g = F.ilaplace(t)
            %   g = F.ilaplace(t, s)
            %
            % 说明：
            %   默认复频域变量：优先取表达式中的符号 s，否则取首个自由符号；
            %   默认时间变量为 t。
            narginchk(1, 3);
            sp = sympy__sp;
            exprSp = getObj(expr);

            if nargin == 1
                tOut = sp.Symbol("t");
                sIn = sympy.Expr.ilaplaceDefaultLaplaceVar(exprSp, sp);
            elseif nargin == 2
                tOut = sympy.to_sp_object(varargin{1});
                sIn = sympy.Expr.ilaplaceDefaultLaplaceVar(exprSp, sp);
            else
                tOut = sympy.to_sp_object(varargin{1});
                sIn = sympy.to_sp_object(varargin{2});
            end

            fSp = sp.inverse_laplace_transform(exprSp, sIn, tOut);
            try
                fSp = fSp.simplify();
            catch
            end
            f = expr.setObj(fSp);
        end

        % ===================== 数学函数重载 =====================
        function res = abs(a)
            % 重载绝对值函数：abs(a)
            sp = sympy__sp;
            res = sympy.Expr(sp.Abs(getObj(a)));
        end

        function res = sin(a)
            sp = sympy__sp;
            res = sympy.Expr(sp.sin(getObj(a)));
        end

        function res = sinh(a)
            sp = sympy__sp;
            res = sympy.Expr(sp.sinh(getObj(a)));
        end

        function res = asin(a)
            sp = sympy__sp;
            res = sympy.Expr(sp.asin(getObj(a)));
        end

        function res = asinh(a)
            sp = sympy__sp;
            res = sympy.Expr(sp.asinh(getObj(a)));
        end

        function res = cos(a)
            sp = sympy__sp;
            res = sympy.Expr(sp.cos(getObj(a)));
        end

        function res = cosh(a)
            sp = sympy__sp;
            res = sympy.Expr(sp.cosh(getObj(a)));
        end

        function res = acos(a)
            sp = sympy__sp;
            res = sympy.Expr(sp.acos(getObj(a)));
        end

        function res = acosh(a)
            sp = sympy__sp;
            res = sympy.Expr(sp.acosh(getObj(a)));
        end

        function res = tan(a)
            sp = sympy__sp;
            res = sympy.Expr(sp.tan(getObj(a)));
        end

        function res = tanh(a)
            sp = sympy__sp;
            res = sympy.Expr(sp.tanh(getObj(a)));
        end

        function res = atan(a)
            sp = sympy__sp;
            res = sympy.Expr(sp.atan(getObj(a)));
        end

        function res = atanh(a)
            sp = sympy__sp;
            res = sympy.Expr(sp.atanh(getObj(a)));
        end

        function res = exp(a)
            sp = sympy__sp;
            res = sympy.Expr(sp.exp(getObj(a)));
        end

        function res = log(a)
            sp = sympy__sp;
            res = sympy.Expr(sp.log(getObj(a)));
        end

        function res = sqrt(a)
            sp = sympy__sp;
            res = sympy.Expr(sp.sqrt(getObj(a)));
        end

        function res = factorial(a)
            sp = sympy__sp;
            res = sympy.Expr(sp.factorial(getObj(a)));
        end

        function y = heaviside(x)
            %HEAVISIDE - Heaviside step function.
            narginchk(1, 1);
            sp = sympy__sp();
            x = sympy.to_sp_object(x);
            y = sympy.Expr(sp.Heaviside(x));
        end
        
    end

    methods (Static)
        function R = cumprod(A, dim)
            if nargin < 2
                R = sympy.Matrix.cumprod(A);
            else
                R = sympy.Matrix.cumprod(A, dim);
            end
        end

        function R = cumsum(A, dim)
            if nargin < 2
                R = sympy.Matrix.cumsum(A);
            else
                R = sympy.Matrix.cumsum(A, dim);
            end
        end

        function ret = symprod(expr, varargin)
            narginchk(2, 4);
            sp = sympy__sp;
            py = pyroot;

            expr = sympy.to_sp_object(expr);
            var = sympy.to_sp_object(varargin{1});

            idxName = char(py.str(var));
            idx = sp.Symbol(idxName, pyargs("integer", true, "positive", true));
            exprIdx = expr.subs(var, idx);

            if nargin == 2
                a = 1;
                b = var;
            elseif nargin == 3
                lim = varargin{2};
                if numel(lim) ~= 2
                    error("求积界限错误.");
                end
                a = sympy.Expr.toBound(lim(1), sp);
                b = sympy.Expr.toBound(lim(2), sp);
            else
                a = sympy.Expr.toBound(varargin{2}, sp);
                b = sympy.Expr.toBound(varargin{3}, sp);
            end

            if sympy.Expr.isPosInf(b, sp, py)
                N = sp.Symbol("N", pyargs("integer", true, "positive", true));
                finiteRet = sp.product(exprIdx, pytuple(idx, a, N)).doit();
                ret = sp.limit(finiteRet, N, sp.oo);
            else
                ret = sp.product(exprIdx, pytuple(idx, a, b)).doit();
            end
            ret = sympy.Expr(ret.simplify());
        end

        function s = vpasum(expr, varargin)
            narginchk(2, 4);
            sp = sympy__sp;
            exprSp = sympy.to_sp_object(expr);

            if nargin == 2
                var = getObj(sympy.Expr(exprSp).symvar(1));
                a = 1;
                b = sympy.to_sp_object(varargin{1});
            elseif nargin == 3
                var = getObj(sympy.Expr(exprSp).symvar(1));
                a = sympy.Expr.toBound(varargin{1}, sp);
                b = sympy.Expr.toBound(varargin{2}, sp);
            else
                var = sympy.to_sp_object(varargin{1});
                if numel(varargin) == 2
                    lim = varargin{2};
                    if numel(lim) ~= 2
                        error("求和界限错误.");
                    end
                    a = sympy.Expr.toBound(lim(1), sp);
                    b = sympy.Expr.toBound(lim(2), sp);
                else
                    a = sympy.Expr.toBound(varargin{2}, sp);
                    b = sympy.Expr.toBound(varargin{3}, sp);
                end
            end

            s = sp.summation(exprSp, pytuple(var, a, b)).evalf();
            s = sympy.Expr(s);
        end

        function y = dirac(varargin)
            %DIRAC - 狄拉克 δ；勿使用顶层 dirac（与 MATLAB Symbolic 冲突），请用 sympy.Expr.dirac。
            narginchk(1, 2);
            sp = sympy__sp();
            if nargin == 1
                x = sympy.to_sp_object(varargin{1});
                y = sympy.Expr(sp.DiracDelta(x));
            else
                n = varargin{1};
                x = sympy.to_sp_object(varargin{2});
                y = sympy.Expr(sp.DiracDelta(x, n));
            end
        end
    end

    methods (Static, Access=private)
        function out = toScalarSp(x, sp)
            if isa(x, "sympy.Base")
                out = getObj(x);
                return;
            end

            if ~(isnumeric(x) && isscalar(x))
                error("输入必须是标量数值或符号对象");
            end

            if isreal(x)
                out = sympy.Expr.toBound(x, sp);
            else
                re = sympy.to_sp_object(real(x));
                im = sympy.to_sp_object(imag(x));
                out = re + im * sp.I;
            end
        end

        function out = toBound(x, sp)
            if isnumeric(x) && isscalar(x) && isinf(x)
                if x > 0
                    out = sp.oo;
                else
                    out = -sp.oo;
                end
            else
                out = sympy.to_sp_object(x);
            end
        end

        function tf = isPosInf(x, sp, py)
            try
                tf = logical(py.bool(x.equals(sp.oo)));
            catch
                tf = false;
            end
        end

        function w = ifourierDefaultFreqVar(exprSp, sp)
            py = pyroot;
            fs = py.list(exprSp.free_symbols);
            n = int64(py.len(fs));
            for i = 0:n - 1
                s = fs{i};
                if strcmp(char(s.name), "w")
                    w = s;
                    return;
                end
            end
            if n == 0
                w = sp.Symbol("w");
            else
                fs_sorted = py.list(fs);
                fs_sorted.sort();
                w = fs_sorted{0};
            end
        end

        function t = laplaceDefaultTimeVar(exprSp, sp)
            % 与 MATLAB laplace 默认一致：若含 t 则对 t 变换，否则取首个符号（字母序）。
            py = pyroot;
            fs = py.list(exprSp.free_symbols);
            n = int64(py.len(fs));
            for i = 0:n - 1
                sym = fs{i};
                if strcmp(char(sym.name), "t")
                    t = sym;
                    return;
                end
            end
            if n == 0
                error("laplace: 表达式中无符号变量");
            end
            fs_sorted = py.list(fs);
            fs_sorted.sort();
            t = fs_sorted{0};
        end

        function s = ilaplaceDefaultLaplaceVar(exprSp, sp)
            py = pyroot;
            fs = py.list(exprSp.free_symbols);
            n = int64(py.len(fs));
            for i = 0:n - 1
                sym = fs{i};
                if strcmp(char(sym.name), "s")
                    s = sym;
                    return;
                end
            end
            if n == 0
                s = sp.Symbol("s");
            else
                fs_sorted = py.list(fs);
                fs_sorted.sort();
                s = fs_sorted{0};
            end
        end
    end

end