% AllSeriesFiles.txt   by Richard Neidinger, Davidson College,  6/2/04
% @series directory
% Path ->  C:\Documents and Settings\rineidinger\My Documents\RESEARCH\autodiff\matlab\AD_classes\@series
% 
%   asin     - overloads arcsine of a series object argument
%   atan     - overloads arctangent of a series object argument
%   cos      - overloads cosine of a series object argument
%   display  - Command window display of a series object.
%   double   - Convert series object to vector of doubles
%   exp      - overloads exp with a series object argument
%   log      - overloads natural logarithm of a series object argument
%   minus    - overloads subtraction with at least one series object argument
%   mpower   - overloads power with at least one series object argument.
%   mrdivide - overloads division with at least one series object argument
%   mtimes   - overloads multiplication with at least one series object argument
%   plus     - overloads addition with at least one series object argument
%   series   - series class constructor
%   sin      - overloads sine with a series object argument
%   sincos   - overloads exp with a series object argument
%   sqrt     - overloads square root of a series object argument
%   tan      - overloads tangent of a series object argument
%   uminus   - overloads negation with a series object argument

function x = series(a,c,d)
%SERIES series class constructor
%   by Richard Neidinger 5/27/04, class methods complete 6/2/04
%   x = series(a,c,d) creates a series object for a Talyor polynomial with
%   constant term a, following coefficient(s) in c, and zeros thru degree d
%   put into fields:
%       .val for the constant term value of the function
%       .coef for the vector of coefficients from linear term to highest term
%   h = series(a,c) will assume that c contains all desired coefficients
%   x = series(a,1,d) creates a series for a variable x at value a.
if nargin == 0
    x.val = [];
    x.coef = [];
    x = class(x,'series');
elseif isa(a,'series')
    x = a;
else
    x.val = a;
    x.coef = c;
    if nargin >= 3 
        x.coef = [x.coef zeros(1,d-length(c))];
    end
    x = class(x,'series');
end

function coefs = double(f)
%SERIES/DOUBLE Convert series object to vector of doubles
%   containing all the series coefficient values.
coefs = [f.val f.coef];

function display(f)
%SERIES/DISPLAY Command window display of a series object.
double(f)

function h = plus(u,v)
%SERIES/PLUS overloads addition with at least one series object argument
% assumes .coef lists are of same length
if ~isa(u,'series')
    h = series(u + v.val, v.coef);
elseif ~isa(v,'series')
    h = series(v + u.val, u.coef);
else
    h = series(u.val + v.val, u.coef + v.coef);
end

function h = uminus(u)
%SERIES/UMINUS overloads negation with a series object argument
h = series(-u.val, -u.coef);

function h = minus(u,v)
%SERIES/MINUS overloads subtraction with at least one series object argument
% assumes .coef lists are of same length
if ~isa(u,'series')
    h = series(u - v.val, v.coef);
elseif ~isa(v,'series')
    h = series(u.val - v, u.coef);
else
    h = series(u.val - v.val, u.coef - v.coef);
end

function h = mtimes(u,v)
%SERIES/MTIMES overloads multiplication with at least one series object argument
% assumes .coef lists are of same length
if ~isa(u,'series')
    h = series(u*v.val, u*v.coef);
elseif ~isa(v,'series')
    h = series(v*u.val, v*u.coef);
else
    d = length(u.coef);
    uall = [u.val u.coef];
    vall = [v.val v.coef];
    h = series(u.val*v.val,0,d);
    for k = 1:d
        h.coef(k) = dot( fliplr(uall(1:k+1)), vall(1:k+1) );
    end
end

function h = mrdivide(u,v)
%SERIES/MRDIVIDE overloads division with at least one series object argument
% assumes .coef lists are of same length
if ~isa(v,'series')  % assume v is a scalar
    h = series(u.val/v, u.coef/v);
else
    d = length(v.coef);
    if ~isa(u,'series')  % assume u is a scalar
        u = series(u,0,d);
    end
    hAll = u.val/v.val;
    for k = 1:d
        hk = (u.coef(k) - dot( hAll, v.coef(fliplr(1:k)) ))/v.val;
        hAll = [hAll hk];
    end
    h = series(hAll(1), hAll(2:d+1));
end

function h = sqrt(u)
%SERIES/SQRT overloads square root of a series object argument
d = length(u.coef);
th0 = 2*sqrt(u.val);
hcoef = zeros(1,0);
for k = 1:d
    hk = (u.coef(k) - dot( hcoef, fliplr(hcoef) )) / th0;
    hcoef = [hcoef hk];
end
h = series(th0/2, hcoef);

function h = exp(u)
%SERIES/EXP overloads exp with a series object argument
d = length(u.coef);
up = (1:d).*u.coef;
hAll = exp(u.val);
for k = 1:d
    hk = dot( hAll, up(fliplr(1:k)) ) / k;
    hAll = [hAll hk];
end
h = series(hAll(1), hAll(2:d+1));

function h = log(u)
%SERIES/LOG overloads natural logarithm of a series object argument
d = length(u.coef);
hp = u.coef(1)/u.val;
for k = 2:d
    hpk = (k*u.coef(k) - dot( hp, u.coef(fliplr(1:k-1)) )) / u.val;
    hp = [hp hpk];
end
h = series( log(u.val), hp./(1:d) );

function h = mpower(u,r)
%SERIES/MPOWER overloads power with at least one series object argument.
% Assumes .coef lists are of same length, if both are series objects.
if isa(r,'series')
    h = exp(r*log(u));  % if u is scalar, only exp is series "convolution"
                        % if u and r are series, exp, log, and * are "convolutions"
else % assume r is a scalar
    if r == 2
        h = u*u;
    elseif r == 1/2
        h = sqrt(u);
    else % about same computation as exp&log but this allows negative base
        d = length(u.coef);
        up = (1:d).*u.coef;
        hAll = u.val^r;
        hp = zeros(1,0);
        for k = 1:d
            hpk = (r*dot(hAll,up(fliplr(1:k))) - dot(hp,u.coef(fliplr(1:k-1))))/u.val;
            hp = [hp hpk];
            hAll = [hAll (hpk/k)];
        end
        h = series(hAll(1), hAll(2:d+1));
    end
end

function [s, c] = sincos(u)
%SERIES/SINCOS overloads exp with a series object argument
d = length(u.coef);
up = (1:d).*u.coef;
sAll = sin(u.val);
cAll = cos(u.val);
for k = 1:d
    sk = dot( cAll, up(fliplr(1:k)) ) / k;
    ck = -dot( sAll, up(fliplr(1:k)) ) / k;
    sAll = [sAll sk];
    cAll = [cAll ck];
end
s = series(sAll(1), sAll(2:d+1));
c = series(cAll(1), cAll(2:d+1));

function h = sin(u)
%SERIES/SIN overloads sine with a series object argument
[h, g] = sincos(u);

function g = cos(u)
%SERIES/COS overloads cosine of a series object argument
[h, g] = sincos(u);

function h = tan(u)
%SERIES/TAN overloads tangent of a series object argument
d = length(u.coef);
up = (1:d).*u.coef;
hAll = tan(u.val);
hpfliplr = zeros(1,0);
vAll = (cos(u.val))^(-2);  % sec here calls series\sec if in the series directory
for k = 1:d
    hpk = dot(vAll,up(fliplr(1:k)));
    hpfliplr = [hpk hpfliplr];
    vk = 2*dot(hAll,hpfliplr)/k;
    hAll = [hAll (hpk/k)];
    vAll = [vAll vk];
end
h = series(hAll(1), hAll(2:d+1));

function h = asin(u)
%SERIES/ASIN overloads arcsine of a series object argument
d = length(u.coef);
h0 = asin(u.val);
v0 = cos(h0);
hp = zeros(1,0);
vcoeffliplr = zeros(1,0);
ufliplr = [fliplr(u.coef) u.val];
for k = 1:d
    hpk = (k*u.coef(k) - dot(hp,vcoeffliplr))/v0;
    hp = [hp hpk];
    vpk = -dot(hp,ufliplr(d+2-k:d+1));
    vcoeffliplr = [vpk/k vcoeffliplr];
end
h = series(h0,hp./(1:d));

function h = atan(u)
%SERIES/ATAN overloads arctangent of a series object argument
d = length(u.coef);
up = (1:d).*u.coef;
h0 = atan(u.val);
v0 = (cos(h0))^(-2);
hp = zeros(1,0);
vcoeffliplr = zeros(1,0);
ufliplr = [fliplr(u.coef) u.val];
for k = 1:d
    hpk = (up(k) - dot(hp,vcoeffliplr))/v0;
    hp = [hp hpk];
    vpk = 2*dot(up(1:k),ufliplr(d+2-k:d+1));
    vcoeffliplr = [vpk/k vcoeffliplr];
end
h = series(h0,hp./(1:d));