substitution in denominator is skipped

[dkrenn]

The following comes very unexpected

sage: ((1+x^2)/x^2).subs({x^2: 42})
43/x^2

The problem is the internal representation as

sage: ((1+x^2)/x^2).operands()
[x^2 + 1, x^(-2)]

Reported upstream as ​https://github.com/pynac/pynac/issues/186

[rws]

Note that there is the walkaround of substituting twice, i.e. as given followed by a substitution of the inverse pattern with the inverse replacement, e.g., ((1+x^2)/(x^2)).subs({x^2: 42}).subs({x^-2: 1/42})

[rws]

Remaining fail, hopefully the author(s) can see where they worked around this bug so that their workaround can be removed:

File "src/sage/rings/asymptotic/asymptotics_multivariate_generating_functions.py", line 1717, in sage.rings.asymptotic.asymptotics_multivariate_generating_functions.FractionWithFactoredDenominator.?
Failed example:
    asy # long time
Expected:
    (4/3*sqrt(3)*sqrt(r)/sqrt(pi) + 47/216*sqrt(3)/(sqrt(pi)*sqrt(r)),
     1, 4/3*sqrt(3)*sqrt(r)/sqrt(pi) + 47/216*sqrt(3)/(sqrt(pi)*sqrt(r)))
Got:
    (4/3*sqrt(3)*sqrt(r)/sqrt(pi) + 587/216*sqrt(3)/(sqrt(pi)*sqrt(r)),
     1,
     4/3*sqrt(3)*sqrt(r)/sqrt(pi) + 587/216*sqrt(3)/(sqrt(pi)*sqrt(r)))
**********************************************************************
File "src/sage/rings/asymptotic/asymptotics_multivariate_generating_functions.py", line 1720, in sage.rings.asymptotic.asymptotics_multivariate_generating_functions.FractionWithFactoredDenominator.?
Failed example:
    F.relative_error(asy[0], alpha, [1, 2, 4, 8], asy[1]) # long time
Expected:
    [((3, 3, 2), 0.9812164307, [1.515572606], [-0.54458543...]),
     ((6, 6, 4), 1.576181132, [1.992989399], [-0.26444185...]),
     ((12, 12, 8), 2.485286378, [2.712196351], [-0.091301338...]),
     ((24, 24, 16), 3.700576827, [3.760447895], [-0.016178847...])]
Got:
    [((3, 3, 2), 0.9812164307, [3.958585166], [-3.034364939]),
     ((6, 6, 4), 1.576181132, [3.720460147], [-1.360426775]),
     ((12, 12, 8), 2.485286378, [3.933702631], [-0.5827965202]),
     ((24, 24, 16), 3.700576827, [4.624183269], [-0.2495844526])]

---

New commits:

​67f4651

21071: substitution in denominator is skipped

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​deeea0a

21071: refine algorithm, fixes fail

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​5e39b7b

Merge branch 'develop' into t/21071/substitution_in_denominator_is_skipped

[tscrim]

I'm not completely sure I agree with this:

     sage: eq.substitute(a=x, x=1)
-    x + 1 == sin(1/x)
+    x + 1 == sin(1)

On the LHS, you simplify x then a (or simultaneously where one substitution does not affect the other), whereas on the RHS, you currently do a, then x (with the new x from the a). So I see this as an inconsistency. IMO, they should be done simultaneously, and the RHS should be sin(1/x).

[rws]

I agree, but the inconsistency already exists in develop:

sage: sin(x/a).subs({1/x : 1, a : x, x : 1})
sin(1)
sage: sin(x/a).subs({a^-1 : x^-1, x^-1 : 1})
sin(1)

(EDIT: simplified)

[rws]

However, the following shows that it is not due to sequential application but premature evaluation (x/x --> 1).

sage: sin(x/a).subs({a^-1:x^-1,x^-1:pi})
sin(1)

[nbruin]

Quoting from ​http://www.ginac.de/reference/basic_8cpp_source.html#l00606, it looks like this is pretty fundamental:

  606 ex basic::subs(const exmap & m, unsigned options) const
  607 {
  608     size_t num = nops();
  609     if (num) {
  610 
  611         // Substitute in subexpressions
  612         for (size_t i=0; i<num; i++) {
  613             const ex & orig_op = op(i);
  614             const ex & subsed_op = orig_op.subs(m, options);
  615             if (!are_ex_trivially_equal(orig_op, subsed_op)) {
  616 
  617                 // Something changed, clone the object
  618                 basic *copy = duplicate();
  619                 copy->clearflag(status_flags::hash_calculated | status_flags::expanded);
  620 
  621                 // Substitute the changed operand
  622                 copy->let_op(i++) = subsed_op;
  623 
  624                 // Substitute the other operands
  625                 for (; i<num; i++)
  626                     copy->let_op(i) = op(i).subs(m, options);
  627 
  628                 // Perform substitutions on the new object as a whole
  629                 return copy->subs_one_level(m, options);
  630             }
  631         }
  632     }
  633 
  634     // Nothing changed or no subexpressions
  635     return subs_one_level(m, options);
  636 }

The problem is that substitutions are done on the subexpressions and *then* on the entire expression. So, we have:

sage: (x/a).subs({a:x,x:1})
1/x
sage: (x/a).subs({a:x,x:1,1/x:17})
17

I'm pretty sure this is because of line 629. This is not because of premature evaluation (which, as you show, can also be an issue)

[rws]

I would like to solve this outside Pynac via an intermediate substitution step. It would depend on #22102.

[rws]

There is an ambiguity. What is the expected result of

sage: f = piecewise([((0,2), x), ([-2,0], -x)]
sage: (x+f).subs(x==1)

piecewise(x|-->1 on (0, 2), x|-->-1 on [-2, 0]; x) + 1
or
2?

Plot expects this to be 2.

[rws]

I think subs should do nothing inside piecewise and only evaluate if the main variable is substituted. If one wishes to change the piece definition one should use another method that is still to be implemented, like e.g. substitute_piece. Otherwise I don't see how to untangle this issue.

[tscrim]

I think what it should do is return the constant function 3 as this would be in line with how subs works on general functions/expressions:

sage: f(x) = x^2
sage: f.subs(x=2)
x |--> 4
sage: f
x |--> x^2
sage: (f+2).subs(x=2)
x |--> 6
sage: (f(x)+2).subs(x=2)
6

[nbruin]

Substitution in computer algebra systems doesn't tend to be semantically defined. It's an operation on syntax trees. So I think piecewise functions can be treated in two ways:

Either you regard them as "atomic" objects. Then substitution doesn't do anything within them. That's not unreasonable. A piecewise function cannot be evaluated at a symbolic input either ...

Or you regard them as a syntax tree (probably a list of (a,b,expression) tuples or so) and you do a substitution on that (hopefully raising an error when the substitution leads to something illegal).