Bug in series expansion

[peter....@gmail.com]

restart;

egf := (1 + x*exp(x*z+z))/(x + 1);
egf1 := (1 + x*exp(x*z)*exp(z))/(x + 1);

# Maple believes that both expressions are equal:

simplify(egf - egf1);

# If you expand the exponential gf, it works with 'egf':

c := n -> expand(n!*coeff(series(egf, z, 9), z, n)):
for n from 0 to 6 do seq(coeff(c(n), x, k), k = 0..n) od;

# but not if you use 'egf1' instead of 'egf' in c:
# Error, unable to compute coeff

[acer]

Your use of `expand` is inadequate in that case. Using `normal` instead of `expand`, both cases work.

egf := (1 + x*exp(x*z+z))/(x + 1):

egf1 := (1 + x*exp(x*z)*exp(z))/(x + 1):

simplify(egf - egf1);
0

c := n -> normal(n!*coeff(series(egf, z, 9), z, n)):

for n from 0 to 6 do seq(coeff(c(n), x, k), k = 0..n) od;

1

0, 1

0, 1, 1

0, 1, 2, 1

0, 1, 3, 3, 1

0, 1, 4, 6, 4, 1

0, 1, 5, 10, 10, 5, 1

c := n -> normal(n!*coeff(series(egf1, z, 9), z, n)):

for n from 0 to 6 do seq(coeff(c(n), x, k), k = 0..n) od;

1

0, 1

0, 1, 1

0, 1, 2, 1

0, 1, 3, 3, 1

0, 1, 4, 6, 4, 1

0, 1, 5, 10, 10, 5, 1

[peter....@gmail.com]

> Your use of `expand` is inadequate in that case.

Dear acer,
could you please explain this in more detail?

[acer]

The result series(egf1, z, 9) is of type `polynom(anything,z)`, but its coefficients are not of type `polynom(anything,x)`. They are rational polynomials with `x+1` as denominator.

egf1 := (1 + x*exp(x*z)*exp(z))/(x + 1):
series(egf1, z, 9);
series(1+(x^2+x)/(x+1)*z+(1/2*x+x^2+1/2*x^3)/(x+1)*z^2+(1/6*x+1/2*x^2+1/2*x^3+1
/6*x^4)/(x+1)*z^3+(1/24*x+1/6*x^2+1/4*x^3+1/6*x^4+1/24*x^5)/(x+1)*z^4+(1/120*x+
1/24*x^2+1/12*x^3+1/12*x^4+1/24*x^5+1/120*x^6)/(x+1)*z^5+(1/720*x+1/120*x^2+1/
48*x^3+1/36*x^4+1/48*x^5+1/120*x^6+1/720*x^7)/(x+1)*z^6+(1/5040*x+1/720*x^2+1/
240*x^3+1/144*x^4+1/144*x^5+1/240*x^6+1/720*x^7+1/5040*x^8)/(x+1)*z^7+(1/40320*
x+1/5040*x^2+1/1440*x^3+1/720*x^4+1/576*x^5+1/720*x^6+1/1440*x^7+1/5040*x^8+1/
40320*x^9)/(x+1)*z^8+O(z^9),z,9)

However, if `normal` is applied to them then the results are of type `polynom(anything,x)` (due to factorization). The calls `coeff(c(n), x, k)` will error out if those coefficients have not been thusly simplified/normalized, because those `coeff` calls expect a form that is of type `polynom(anything,x)`.

normal(series(egf1, z, 9));
series(1+x*z+1/2*x*(x+1)*z^2+1/6*(x^2+2*x+1)*x*z^3+1/24*x*(x^3+3*x^2+3*x+1)*z^4
+1/120*x*(x^4+4*x^3+6*x^2+4*x+1)*z^5+1/720*x*(x^5+5*x^4+10*x^3+10*x^2+5*x+1)*z^\
6+1/5040*x*(x^6+6*x^5+15*x^4+20*x^3+15*x^2+6*x+1)*z^7+1/40320*x*(x^7+7*x^6+21*x
^5+35*x^4+35*x^3+21*x^2+7*x+1)*z^8+O(z^9),z,9)

The same thing occurs if you take the coefficients one at a time. For example,

b := n -> n!*coeff(series(egf1, z, 9), z, n):
b(2);
2*(1/2*x+x^2+1/2*x^3)/(x+1)

type(b(2), polynom(anything,x));
false

coeff(b(2), x, 1);

Error, unable to compute coeff

Now, let's repeat that, but with `normal`. (Using `simplify` or `factor` would also work here.)

normal(b(2));
x*(x+1)

type(normal(b(2)), polynom(anything,x));
true

coeff(normal(b(2)), x, 1);
1

Or, as an adjustment to your original `c`,

c := n -> normal(n!*coeff(series(egf1, z, 9), z, n)):

c(2):lprint(%);
x*(x+1)

type(c(2), polynom(anything,x));
true

coeff(c(2), x, 1);
1

Your original `c` tried it using `expand` instead of `normal`. That is not sufficient, ie.

expand(b(2));
1/(x+1)*x+2/(x+1)*x^2+1/(x+1)*x^3

It happens that the coefficients of `series(egf, z, 9)` happen to come out in a form that is of type `polynom(anything,x)`.

egf := (1 + x*exp(x*z+z))/(x + 1):
series(egf, z, 9);
series(1+x*z+1/2*x*(x+1)*z^2+1/6*x*(x+1)^2*z^3+1/24*x*(x+1)^3*z^4+1/120*x*(x+1)
^4*z^5+1/720*x*(x+1)^5*z^6+1/5040*x*(x+1)^6*z^7+1/40320*x*(x+1)^7*z^8+O(z^9),z,\
9)

An alternative is to use `MultiSeries:-series`, producing the desired form of the coefficients for both `egf` and `egf1`.

MultiSeries:-series(egf1, z, 9);
series(1+x*z+1/2*x*(x+1)*z^2+1/6*(x^2+2*x+1)*x*z^3+1/24*x*(x^3+3*x^2+3*x+1)*z^4
+1/120*x*(x^4+4*x^3+6*x^2+4*x+1)*z^5+1/720*x*(x^5+5*x^4+10*x^3+10*x^2+5*x+1)*z^\
6+1/5040*x*(x^6+6*x^5+15*x^4+20*x^3+15*x^2+6*x+1)*z^7+1/40320*x*(x^7+7*x^6+21*x
^5+35*x^4+35*x^3+21*x^2+7*x+1)*z^8+O(z^9),z,9)

for n from 0 to 6 do seq(coeff(d(n), x, k), k = 0..n) od;

[peter....@gmail.com]

acer:

> The result series(egf1, z, 9) is of type `polynom(anything,z)`, but its coefficients are not of type `polynom(anything,x)`. They are rational polynomials with `x+1` as denominator.

Thanks for this instructive explanation.

[peter....@gmail.com]

> The result series(egf1, z, 9) is of type `polynom(anything,z)`,
> but its coefficients are not of type `polynom(anything,x)`.
> They are rational polynomials with `x+1` as denominator.

Your explanation for Maple's behavior is interesting, but of course
you dodge the question of how to evaluate this behavior.

I think we agree that

x = y -> f(x) = f(y) for a function f.

To what extent does this not also apply to

simplify(x - y) = 0 -> f(x) = f(y) ?

Maybe you can't guarantee that in all cases, but in such a
simple example as the one shown here, I think it's essential.
I just expect that as a user. And other CAS deliver that too.

Here's what the example looks like in SageMath:

x, z = var("x, z")

egf = (1 + x*exp(x*z+z))/(x + 1)

egf1 = (1 + x*exp(x*z)*exp(z))/(x + 1)

simplify(egf - egf1)

# If you expand the exponential gfs, it works
# with 'egf' as well as with 'egf1'.

t = taylor(egf, z, 0, 9)

for n in range(7):
print((factorial(n)*t.coefficient(z, n)).list())

[1]

[0, 1]
[0, 1, 1]
[0, 1, 2, 1]
[0, 1, 3, 3, 1]
[0, 1, 4, 6, 4, 1]

[0, 1, 5, 10, 10, 5, 1]

This is the output in both cases; in the exact same Jupyter
environment, Maple would only display the last line, which
I think is another bug.