Opened 6 years ago

Last modified 6 years ago

#21071 needs_work defect

substitution in denominator is skipped

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

Status badges

Description (last modified by Ralf Stephan)

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 6 years ago by Ralf Stephan

Description: modified (diff)
Report Upstream: N/AReported upstream. Developers acknowledge bug.
Summary: substitution in denomintor is skippedsubstitution in denominator is skipped

comment:2 Changed 6 years ago by Ralf Stephan

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 6 years ago by Ralf Stephan

Branch: u/rws/substitution_in_denominator_is_skipped

comment:4 Changed 6 years ago by Ralf Stephan

Authors: Ralf Stephan
Cc: Travis Scrimshaw Daniel Krenn added
Commit: 67f46511fb82b1c96182c212cb6112bceda240e0
Milestone: sage-7.3sage-7.4
Report Upstream: Reported upstream. Developers acknowledge bug.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 6 years ago by git

Commit: 67f46511fb82b1c96182c212cb6112bceda240e0deeea0ac1ab28e246944059bbd576c554a6f4579

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

deeea0a21071: refine algorithm, fixes fail

comment:6 Changed 6 years ago by Ralf Stephan

Status: newneeds_review

comment:7 Changed 6 years ago by git

Commit: deeea0ac1ab28e246944059bbd576c554a6f45795e39b7b65f7f3ae6e778adacf9514f137c930900

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 6 years ago by Travis Scrimshaw

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 6 years ago by Ralf Stephan

I agree, but the inconsistency already exists in develop:

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

(EDIT: simplified)

Version 1, edited 6 years ago by Ralf Stephan (previous) (next) (diff)

comment:10 Changed 6 years ago by Ralf Stephan

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 6 years ago by Ralf Stephan (previous) (diff)

comment:11 Changed 6 years ago by Nils Bruin

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 6 years ago by Ralf Stephan

Dependencies: #22102

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

comment:13 Changed 6 years ago by Ralf Stephan

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 6 years ago by Ralf Stephan (previous) (diff)

comment:14 Changed 6 years ago by Ralf Stephan

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 6 years ago by Ralf Stephan

Status: needs_reviewneeds_work

comment:16 Changed 6 years ago by Travis Scrimshaw

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 6 years ago by Nils Bruin

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.