Opened 7 years ago

Closed 4 years ago

#14326 closed defect (duplicate)

Substituting numeric one in symbolic expression gives symbolic one

Reported by: zimmerma Owned by: AlexGhitza
Priority: major Milestone: sage-duplicate/invalid/wontfix
Component: basic arithmetic Keywords:
Cc: mmezzarobba, burcin Merged in:
Authors: Reviewers:
Report Upstream: N/A Work issues:
Branch: Commit:
Dependencies: Stopgaps:

Description

consider the following in Sage 5.8:

sage: u(n) = n^100 / 100^n; u(1.)
1/100

This is inconsistent with:

sage: n=1.; n^100 / 100^n
0.0100000000000000

and with:

sage: v = lambda(n): n^100 / 100^n; v(1.)
0.0100000000000000

and:

sage: def w(n): return n^100 / 100^n       
sage: w(1.) 
0.0100000000000000

Change History (17)

comment:1 Changed 7 years ago by zimmerma

  • Cc mmezzarobba added

comment:2 Changed 7 years ago by zimmerma

Jeroen, do you have an idea who to include in cc to help isolate this?

Paul

comment:3 Changed 6 years ago by jdemeyer

  • Milestone changed from sage-5.11 to sage-5.12

comment:4 Changed 6 years ago by vbraun_spam

  • Milestone changed from sage-6.1 to sage-6.2

comment:5 Changed 6 years ago by vbraun_spam

  • Milestone changed from sage-6.2 to sage-6.3

comment:6 Changed 5 years ago by vbraun_spam

  • Milestone changed from sage-6.3 to sage-6.4

comment:7 Changed 5 years ago by kcrisman

  • Cc burcin added
  • Summary changed from function gives symbolic output for numeric input to Substituting numeric one in symbolic expression gives symbolic one

Wow, this is weird. Here is a much simpler example.

sage: w(n) = n
sage: w(1.)
1.00000000000000
sage: w(n) = n^2
sage: w(1.)
1

In fact, even

sage: (x^2).subs(x=1.)
1

works. Yuck.

Somehow the custom power method is not doing its job when you substitute . But I don't see an obvious place in Ginac where this would get screwed up...

Aha.

sage: z = x^2
sage: z.subs(x=1.)
1
sage: z.subs(x=2.)
4.00000000000000

Because one does get treated differently. Though

sage: z = x
sage: z.subs(x=1.)
1.00000000000000

so it also has something to do with the coercion that happens in the power method for symbolic expressions.

comment:8 Changed 5 years ago by kcrisman

Okay, I think this might be a bug in Ginac, or possibly in how we use Ginac. In the Ginac definition of automatic rewriting of power::eval, we have

00399     // ^(x,1) -> x
00400     if (eexponent.is_equal(_ex1))
00401         return basis;
...
00413     // ^(1,x) -> 1
00414     if (ebasis.is_equal(_ex1))
00415         return _ex1;

The other rewriting rules are probably harmless, though if

sage: z
x^n
sage: z.subs(x=0.)
0.000000000000000^n
sage: z.subs(x=0)
0^n
sage: z.subs(n=0)
1
sage: z.subs(n=0.)
1

where the 0^n business is because Ginac checks if the exponent is numerical because it doesn't want to evaluate something that could, in principle, still become 0^0.

Unfortunately, I'm not sure how to monkey-patch Pynac into recognizing this situation, and I certainly don't want to do a catch in the symbolic expression code, that is really the wrong place. Here's hoping someone really intimately familiar with our back-and-forth to Pynac sees an easy fix.

comment:9 Changed 5 years ago by jdemeyer

The problem looks similar to #17130, but that's about __call__.

comment:10 Changed 5 years ago by zimmerma

do the patches from #17130 solve this ticket?

Paul

comment:11 follow-up: Changed 5 years ago by zimmerma

Dear Karl-Dieter,

I think this might be a bug in Ginac, or possibly in how we use Ginac.

please could you ask the Ginac developers if this is a bug in Ginac? And if not, how to replace _ex1 to get the correct "one"?

Paul

comment:12 in reply to: ↑ 11 Changed 5 years ago by kcrisman

I think this might be a bug in Ginac, or possibly in how we use Ginac.

please could you ask the Ginac developers if this is a bug in Ginac? And if not, how to replace _ex1 to get the correct "one"?

I am really not familiar enough with Ginac proper to do either of these with any technical knowledge, unfortunately. And the Ginac developer(s) are not particularly responsive right now to any but the most informed pieces of information, apparently. If someone can figured out how to replicate this in Ginac that would be great, but again it could be us misusing it, I'm not sure.

I would be surprised if #17130 fixed this, based on my experiments above.

comment:13 follow-up: Changed 5 years ago by burcin

Looking at the code snippet from comment:8 only, I don't think this is a bug in Ginac. The intended behavior is just different in Pynac, so we have to patch Pynac.

Ginac wants to keep unique reference counted expression objects for expressions that are equal. That is why they return _ex1 on line 415. In Pynac, we do not have a unique "one", we should just return ebasis in this line to keep the precision/type of the base.

comment:14 in reply to: ↑ 13 Changed 5 years ago by kcrisman

Replying to burcin:

Thanks for replying!

Looking at the code snippet from comment:8 only, I don't think this is a bug in Ginac. The intended behavior is just different in Pynac, so we have to patch Pynac.

Ginac wants to keep unique reference counted expression objects for expressions that are equal. That is why they return _ex1 on line 415. In Pynac, we do not have a unique "one", we should just return ebasis in this line to keep the precision/type of the base.

Of course! That makes perfect sense - once you say it, before it was murky :(

What about for the other case, of x^1.0? We have

sage: (x^2).subs(x=1.)
1
sage: (2^x).subs(x=1.)
2
sage: (x^2).subs(x=2.)
4.00000000000000
sage: (2^x).subs(x=2.)
4.00000000000000

so it would seem that if the exponent is not exact, we want the whole thing to be numerical. I guess we could just strip out that simplification completely, but I don't know if that would give us anything useful either.

comment:15 Changed 5 years ago by rws

  • Milestone changed from sage-6.4 to sage-duplicate/invalid/wontfix
  • Status changed from new to needs_review

This is fixed in pynac-0.3.6 (and a duplicate of #12257), please review #18362 and #12257.

comment:16 Changed 5 years ago by mmezzarobba

  • Status changed from needs_review to positive_review

comment:17 Changed 4 years ago by vbraun

  • Resolution set to duplicate
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.