• Bug in series expansion

    From [email protected]@21:1/5 to All on Tue Apr 18 03:31:59 2023
    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

    --- SoupGate-Win32 v1.05
    * Origin: fsxNet Usenet Gateway (21:1/5)
  • From acer@21:1/5 to [email protected] on Tue Apr 18 13:27:00 2023
    On Tuesday, April 18, 2023 at 6:32:00 AM UTC-4, [email protected] wrote:
    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

    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

    --- SoupGate-Win32 v1.05
    * Origin: fsxNet Usenet Gateway (21:1/5)
  • From [email protected]@21:1/5 to All on Wed Apr 19 04:42:03 2023
    Your use of `expand` is inadequate in that case.

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

    --- SoupGate-Win32 v1.05
    * Origin: fsxNet Usenet Gateway (21:1/5)
  • From acer@21:1/5 to [email protected] on Wed Apr 19 07:03:41 2023
    On Wednesday, April 19, 2023 at 7:42:05 AM UTC-4, [email protected] wrote:
    Your use of `expand` is inadequate in that case.
    Dear acer,
    could you please explain this in more detail?

    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;
    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

    --- SoupGate-Win32 v1.05
    * Origin: fsxNet Usenet Gateway (21:1/5)
  • From [email protected]@21:1/5 to All on Wed Apr 19 12:03:21 2023
    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.

    --- SoupGate-Win32 v1.05
    * Origin: fsxNet Usenet Gateway (21:1/5)
  • From [email protected]@21:1/5 to All on Mon May 1 00:56:01 2023
    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.

    --- SoupGate-Win32 v1.05
    * Origin: fsxNet Usenet Gateway (21:1/5)