# Maple programs for operations on Hurwitz sequences, based on the paper by
# Xing Gao and William F. Keigher, Interlacing of Hurwitz series, 
# Communications in Algebra, 45:5 (2017), 2163-2185;
# DOI: 10.1080/00927872.2016.1226885

# N. J. A. Sloane, April 12 2018

# Hurwitz product of two sequences
# Also known as the exponential convolution.
Hmul:=proc(a,b) local aa,bb,cc,La,Lb,Lc,i,n;
La:=nops(a); Lb:=nops(b);
if La <= 0 then error("first argument cannot be empty"); fi;
if Lb <= 0 then error("second argument cannot be empty"); fi;
Lc:=min(La,Lb);
aa:=Array(0..La-1,0); bb:=Array(0..Lb-1,0); cc:=Array(0..Lc-1,0);
for n from 1 to La do aa[n-1]:=a[n]; od;
for n from 1 to Lb do bb[n-1]:=b[n]; od;
for n from 0 to Lc-1 do cc[n] := add(binomial(n,i)*aa[i]*bb[n-i],i=0..n); od:
[seq(cc[i],i=0..Lc-1)];
end;

# Hurwitz inverse
Hinv:=proc(a) local a0i,aa,bb,La,i,n;
La:=nops(a);
if La <= 0 then error("argument cannot be empty"); fi;
if a[1]=0 then error("first term must be nonzero"); fi;
aa:=Array(0..La-1,0); bb:=Array(0..La-1,0);
for n from 1 to La do aa[n-1]:=a[n]; od;
a0i:=1/aa[0]; bb[0]:=a0i;
for n from 1 to La-1 do
bb[n] := -a0i*add(binomial(n,i)*bb[i]*aa[n-i],i=0..n-1); od:
[seq(bb[i],i=0..La-1)];
end;

# Hurwitz left shift (or partial derivative)
Hleft:=proc(a) local La,i; La:=nops(a);
if La <= 0 then error("argument cannot be empty"); fi;
if La = 1 then []; else [seq(a[i],i=2..La)]; fi;
end;

# Hurwitz right shift (or integral)
Hright:=proc(a) local La,i; La:=nops(a);
if La = 0 then [0]; else [0,seq(a[i],i=1..La)]; fi;
end;

# Hurwitz log
Hlog:=proc(a) local a0i,La,b1,b2,b3; global Hinv,Hleft,Hright;
La:=nops(a);
if La <= 0 then error("argument cannot be empty"); fi;
if abs(a[1]) <> 1 then error("first term must be +-1"); fi;
b1:=Hleft(a); b2:=Hinv(a); b3:=Hmul(b2,b1);
Hright(b3);
end;

# Hurwitz exponential
Hexp:=proc(a) local i,j,La,b1,b2,t1,dp; global Hinv,Hleft,Hright;
La:=nops(a);
if La <= 0 then error("argument cannot be empty"); fi;
if a[1] <> 0 then error("first term must be 0"); fi;
# divided powers
t1:=[1,seq(0,i=1..La-1)];
dp:=t1;
for j from 1 to La-1 do
b1:=Hleft(a); b2:=Hmul(t1,b1); t1:=Hright(b2);
for i from 1 to La do dp[i] := dp[i]+t1[i]; od;
od;
[seq(dp[i],i=1..La)];
end;

# Hurwitz sine, first M terms
Hsin:= proc(M) local a,i;
a:=[seq(0,i=1..M)];
for i from 1 to floor(M/2) do a[2*i]:=(-1)^(i+1); od;
[seq(a[i],i=1..M)];
end;

# Hurwitz cosine, first M terms
Hcos:= proc(M) local a,i;
a:=[seq(0,i=1..M)];
for i from 0 to floor((M-1)/2) do a[2*i+1]:=(-1)^i; od;
[seq(a[i],i=1..M)];
end;

# Hurwitz secant, first M terms
Hsec := M -> Hinv(Hcos(M));

# Hurwitz tangent, first M terms
Htan := M -> Hmul(Hsin(M), Hsec(M));