Opened 3 years ago

Last modified 3 years ago

#21071 needs_work defect

substitution in denominator is skipped

Reported by: dkrenn Owned by:
Priority: major Milestone: sage-7.4
Component: symbolics Keywords:
Cc: tscrim, dkrenn Merged in:
Authors: Ralf Stephan Reviewers:
Report Upstream: N/A Work issues:
Branch: u/rws/substitution_in_denominator_is_skipped (Commits) Commit: 5e39b7b65f7f3ae6e778adacf9514f137c930900
Dependencies: #22102 Stopgaps:

Description (last modified by rws)

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

Change History (17)

comment:1 Changed 3 years ago by rws

  • Description modified (diff)
  • Report Upstream changed from N/A to Reported upstream. Developers acknowledge bug.
  • Summary changed from substitution in denomintor is skipped to substitution in denominator is skipped

comment:2 Changed 3 years ago by 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})

comment:3 Changed 3 years ago by rws

  • Branch set to u/rws/substitution_in_denominator_is_skipped

comment:4 Changed 3 years ago by rws

  • Authors set to Ralf Stephan
  • Cc tscrim dkrenn added
  • Commit set to 67f46511fb82b1c96182c212cb6112bceda240e0
  • Milestone changed from sage-7.3 to sage-7.4
  • Report Upstream changed from Reported upstream. Developers acknowledge bug. to N/A

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:

67f465121071: substitution in denominator is skipped

comment:5 Changed 3 years ago by git

  • Commit changed from 67f46511fb82b1c96182c212cb6112bceda240e0 to deeea0ac1ab28e246944059bbd576c554a6f4579

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

deeea0a21071: refine algorithm, fixes fail

comment:6 Changed 3 years ago by rws

  • Status changed from new to needs_review

comment:7 Changed 3 years ago by git

  • Commit changed from deeea0ac1ab28e246944059bbd576c554a6f4579 to 5e39b7b65f7f3ae6e778adacf9514f137c930900

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

5e39b7bMerge branch 'develop' into t/21071/substitution_in_denominator_is_skipped

comment:8 Changed 3 years ago by 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).

comment:9 Changed 3 years ago by 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)

Last edited 3 years ago by rws (previous) (diff)

comment:10 Changed 3 years ago by 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)
Last edited 3 years ago by rws (previous) (diff)

comment:11 Changed 3 years ago by 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)

comment:12 Changed 3 years ago by rws

  • Dependencies set to #22102

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

comment:13 Changed 3 years ago by 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.

Last edited 3 years ago by rws (previous) (diff)

comment:14 Changed 3 years ago by 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.

comment:15 Changed 3 years ago by rws

  • Status changed from needs_review to needs_work

comment:16 Changed 3 years ago by 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

comment:17 Changed 3 years ago by 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).

Note: See TracTickets for help on using tickets.