mpowerin.m

来自「计算动力学系统的分岔图」· M 代码 · 共 38 行

function ans = mpowerin(x,y)
% X^Y for adtayl objects.
% The case where X is scalar and Y is an integer.
switch y
  case 1, ans = x;
  case 2, x*x;
  case 3, x*x*x;
  case 4, x*x;  ans*ans; 
  otherwise
    %if ~(isscalar(x) && isnumeric(y) && fix(y)==y)
     % error('X^Y for adtayl objects currently requires X scalar and Y integer')
    %end
    xtc = x.tc;
    if y>0
      pow = power0(xtc,y);
    else
      pow = zeros(size(xtc)); pow(1)=1; % vector [1,0,...,0], representing constant 1.
      if y<0
        pow = filter(1,power0(xtc,-y), pow);
      end
    end
    class(struct('tc',pow),'adtayl');
end

%---------------------------
function ans = power0(xtc,y)
% Assumes xtc is a vector and y>=1.
if y==1
  xtc;
else
 %yover2 = fix(y/2); yrem2 = rem(y,2);
  pyover2 = power0(xtc,fix(y/2));
  filter(pyover2,1, pyover2);
  if rem(y,2)~=0
    filter(xtc,1, ans);
  end
end