search

rtoy / maxima

A Clone of Maxima's repo
Other
0 stars 0 forks source link

a cosine of arcsin limit that is evaluated incorrectly #3910

Open rtoy opened 1 month ago

rtoy commented 1 month ago

Imported from SourceForge on 2024-07-08 21:15:01 Created by davewittemorris on 2022-07-09 18:45:30 Original: https://sourceforge.net/p/maxima/bugs/4004

It is not difficult to see that the value of the following limit is 4*m + 1. However, maxima says the value is 0: 

(%i1) declare(m, integer)$ 
(%i2) assume(equal(cos((4 * %pi * m + %pi)/2), 0))$
(%i3) limit(cos((4*m + 1) * asin(1/sqrt(x^2 + 1)))/abs(x), x, 0);
(%o3)                                  0

This is from sagemath trac ticket #34140.

Maxima version: "5.43.0" Maxima build date: "2019-06-05 13:14:43" Host type: "x86_64-apple-darwin13.4.0" Lisp implementation type: "SBCL" Lisp implementation version: "1.5.3"

rtoy commented 1 month ago

Imported from SourceForge on 2024-07-08 21:15:02 Created by robert_dodier on 2022-07-10 18:43:53 Original: https://sourceforge.net/p/maxima/bugs/4004/#2193

rtoy commented 1 month ago

Imported from SourceForge on 2024-07-08 21:15:05 Created by robert_dodier on 2022-07-10 18:43:54 Original: https://sourceforge.net/p/maxima/bugs/4004/#5c5e

Observed in current version from Git. build_info reports in part:

Maxima version: "branch_5_46_base_242_g503e1ac61_dirty"
Maxima build date: "2022-07-09 07:39:59"
Host type: "i686-pc-linux-gnu"
Lisp implementation type: "SBCL"
Lisp implementation version: "2.1.6"
rtoy commented 1 month ago

Imported from SourceForge on 2024-07-08 21:15:09 Created by robert_dodier on 2022-07-10 19:22:56 Original: https://sourceforge.net/p/maxima/bugs/4004/#b4a1

Looks like the result is coming out of TAYLIM (Lisp function for Taylor limit). taylor is returning something that looks like stuff/x + constant + (more stuff times powers of x) where stuff is the term which asksign was asking about. If you answer zero, then TAYLIM seems to forget about the constant, which, at this point, I believe is the part which gives the correct result.

For the record, I didn't investigate whether TAYLIM does the right thing if you respond something other than zero to asksign.

Here's what I'm seeing. Current Git version as noted above.

(%i2) trace (?taylim, taylor, tlimit);
(%o2) [taylim,taylor,tlimit]
(%i3) limit(cos((4*m + 1) * asin(1/sqrt(x^2 + 1)))/abs(x), x, 0);

1 Enter taylor 
[cos(4*m*asin(1/sqrt(1/(1/W633)^2+1))+asin(1/sqrt(1/(1/W633)^2+1)))/W633,W633,
 0,4]
1 Exit  taylor 
cos((4*%pi*m+%pi)/2)/W633+(4*sin((4*%pi*m+%pi)/2)*m+sin((4*%pi*m+%pi)/2))
                         -((16*cos((4*%pi*m+%pi)/2)*m^2
                          +8*cos((4*%pi*m+%pi)/2)*m+cos((4*%pi*m+%pi)/2))
                          *W633)
                          /2
                         -((64*sin((4*%pi*m+%pi)/2)*m^3
                          +48*sin((4*%pi*m+%pi)/2)*m^2
                          +20*sin((4*%pi*m+%pi)/2)*m+3*sin((4*%pi*m+%pi)/2))
                          *W633^2)
                          /6
                         +((256*cos((4*%pi*m+%pi)/2)*m^4
                          +256*cos((4*%pi*m+%pi)/2)*m^3
                          +224*cos((4*%pi*m+%pi)/2)*m^2
                          +80*cos((4*%pi*m+%pi)/2)*m+9*cos((4*%pi*m+%pi)/2))
                          *W633^3)
                          /24
                         +((1024*sin((4*%pi*m+%pi)/2)*m^5
                          +1280*sin((4*%pi*m+%pi)/2)*m^4
                          +1920*sin((4*%pi*m+%pi)/2)*m^3
                          +1120*sin((4*%pi*m+%pi)/2)*m^2
                          +356*sin((4*%pi*m+%pi)/2)*m+45*sin((4*%pi*m+%pi)/2))
                          *W633^4)
                          /120
1 Enter taylim 
[cos(4*m*asin(1/sqrt(x^2+1))+asin(1/sqrt(x^2+1)))/x,x,zeroa,think]
 1 Enter taylor [cos(4*m*asin(1/sqrt(x^2+1))+asin(1/sqrt(x^2+1)))/x,x,0,4]
 1 Exit  taylor 
 cos((4*%pi*m+%pi)/2)/x+(4*sin((4*%pi*m+%pi)/2)*m+sin((4*%pi*m+%pi)/2))
                       -((16*cos((4*%pi*m+%pi)/2)*m^2
                        +8*cos((4*%pi*m+%pi)/2)*m+cos((4*%pi*m+%pi)/2))
                        *x)
                        /2
                       -((64*sin((4*%pi*m+%pi)/2)*m^3
                        +48*sin((4*%pi*m+%pi)/2)*m^2+20*sin((4*%pi*m+%pi)/2)*m
                        +3*sin((4*%pi*m+%pi)/2))
                        *x^2)
                        /6
                       +((256*cos((4*%pi*m+%pi)/2)*m^4
                        +256*cos((4*%pi*m+%pi)/2)*m^3
                        +224*cos((4*%pi*m+%pi)/2)*m^2
                        +80*cos((4*%pi*m+%pi)/2)*m+9*cos((4*%pi*m+%pi)/2))
                        *x^3)
                        /24
                       +((1024*sin((4*%pi*m+%pi)/2)*m^5
                        +1280*sin((4*%pi*m+%pi)/2)*m^4
                        +1920*sin((4*%pi*m+%pi)/2)*m^3
                        +1120*sin((4*%pi*m+%pi)/2)*m^2
                        +356*sin((4*%pi*m+%pi)/2)*m+45*sin((4*%pi*m+%pi)/2))
                        *x^4)
                        /120
1 Enter taylor 
[cos(4*m*asin(1/sqrt(1/(1/W884)^2+1))+asin(1/sqrt(1/(1/W884)^2+1)))/W884,W884,
 0,4]
1 Exit  taylor 
cos((4*%pi*m+%pi)/2)/W884+(4*sin((4*%pi*m+%pi)/2)*m+sin((4*%pi*m+%pi)/2))
                         -((16*cos((4*%pi*m+%pi)/2)*m^2
                          +8*cos((4*%pi*m+%pi)/2)*m+cos((4*%pi*m+%pi)/2))
                          *W884)
                          /2
                         -((64*sin((4*%pi*m+%pi)/2)*m^3
                          +48*sin((4*%pi*m+%pi)/2)*m^2
                          +20*sin((4*%pi*m+%pi)/2)*m+3*sin((4*%pi*m+%pi)/2))
                          *W884^2)
                          /6
                         +((256*cos((4*%pi*m+%pi)/2)*m^4
                          +256*cos((4*%pi*m+%pi)/2)*m^3
                          +224*cos((4*%pi*m+%pi)/2)*m^2
                          +80*cos((4*%pi*m+%pi)/2)*m+9*cos((4*%pi*m+%pi)/2))
                          *W884^3)
                          /24
                         +((1024*sin((4*%pi*m+%pi)/2)*m^5
                          +1280*sin((4*%pi*m+%pi)/2)*m^4
                          +1920*sin((4*%pi*m+%pi)/2)*m^3
                          +1120*sin((4*%pi*m+%pi)/2)*m^2
                          +356*sin((4*%pi*m+%pi)/2)*m+45*sin((4*%pi*m+%pi)/2))
                          *W884^4)
                          /120
1 Enter taylim 
[-cos(4*m*asin(1/sqrt(x^2+1))+asin(1/sqrt(x^2+1)))/x,x,zerob,think]
 1 Enter taylor [-cos(4*m*asin(1/sqrt(x^2+1))+asin(1/sqrt(x^2+1)))/x,x,0,4]
 1 Exit  taylor 
 (-cos((4*%pi*m+%pi)/2)/x)+((-4*sin((4*%pi*m+%pi)/2)*m)-sin((4*%pi*m+%pi)/2))
                          +((16*cos((4*%pi*m+%pi)/2)*m^2
                           +8*cos((4*%pi*m+%pi)/2)*m+cos((4*%pi*m+%pi)/2))
                           *x)
                           /2
                          +((64*sin((4*%pi*m+%pi)/2)*m^3
                           +48*sin((4*%pi*m+%pi)/2)*m^2
                           +20*sin((4*%pi*m+%pi)/2)*m+3*sin((4*%pi*m+%pi)/2))
                           *x^2)
                           /6
                          -((256*cos((4*%pi*m+%pi)/2)*m^4
                           +256*cos((4*%pi*m+%pi)/2)*m^3
                           +224*cos((4*%pi*m+%pi)/2)*m^2
                           +80*cos((4*%pi*m+%pi)/2)*m+9*cos((4*%pi*m+%pi)/2))
                           *x^3)
                           /24
                          -((1024*sin((4*%pi*m+%pi)/2)*m^5
                           +1280*sin((4*%pi*m+%pi)/2)*m^4
                           +1920*sin((4*%pi*m+%pi)/2)*m^3
                           +1120*sin((4*%pi*m+%pi)/2)*m^2
                           +356*sin((4*%pi*m+%pi)/2)*m
                           +45*sin((4*%pi*m+%pi)/2))
                           *x^4)
                           /120
1 Enter taylor 
[cos(4*m*asin(1/sqrt(1/(1/W1200)^2+1))+asin(1/sqrt(1/(1/W1200)^2+1)))/W1200,
 W1200,0,4]
1 Exit  taylor 
cos((4*%pi*m+%pi)/2)/W1200+(4*sin((4*%pi*m+%pi)/2)*m+sin((4*%pi*m+%pi)/2))
                          -((16*cos((4*%pi*m+%pi)/2)*m^2
                           +8*cos((4*%pi*m+%pi)/2)*m+cos((4*%pi*m+%pi)/2))
                           *W1200)
                           /2
                          -((64*sin((4*%pi*m+%pi)/2)*m^3
                           +48*sin((4*%pi*m+%pi)/2)*m^2
                           +20*sin((4*%pi*m+%pi)/2)*m+3*sin((4*%pi*m+%pi)/2))
                           *W1200^2)
                           /6
                          +((256*cos((4*%pi*m+%pi)/2)*m^4
                           +256*cos((4*%pi*m+%pi)/2)*m^3
                           +224*cos((4*%pi*m+%pi)/2)*m^2
                           +80*cos((4*%pi*m+%pi)/2)*m+9*cos((4*%pi*m+%pi)/2))
                           *W1200^3)
                           /24
                          +((1024*sin((4*%pi*m+%pi)/2)*m^5
                           +1280*sin((4*%pi*m+%pi)/2)*m^4
                           +1920*sin((4*%pi*m+%pi)/2)*m^3
                           +1120*sin((4*%pi*m+%pi)/2)*m^2
                           +356*sin((4*%pi*m+%pi)/2)*m
                           +45*sin((4*%pi*m+%pi)/2))
                           *W1200^4)
                           /120
1 Enter taylim 
[cos(4*m*asin(1/sqrt(x^2+1))+asin(1/sqrt(x^2+1)))/x,x,zeroa,think]
 1 Enter taylor [cos(4*m*asin(1/sqrt(x^2+1))+asin(1/sqrt(x^2+1)))/x,x,0,4]
 1 Exit  taylor 
 cos((4*%pi*m+%pi)/2)/x+(4*sin((4*%pi*m+%pi)/2)*m+sin((4*%pi*m+%pi)/2))
                       -((16*cos((4*%pi*m+%pi)/2)*m^2
                        +8*cos((4*%pi*m+%pi)/2)*m+cos((4*%pi*m+%pi)/2))
                        *x)
                        /2
                       -((64*sin((4*%pi*m+%pi)/2)*m^3
                        +48*sin((4*%pi*m+%pi)/2)*m^2+20*sin((4*%pi*m+%pi)/2)*m
                        +3*sin((4*%pi*m+%pi)/2))
                        *x^2)
                        /6
                       +((256*cos((4*%pi*m+%pi)/2)*m^4
                        +256*cos((4*%pi*m+%pi)/2)*m^3
                        +224*cos((4*%pi*m+%pi)/2)*m^2
                        +80*cos((4*%pi*m+%pi)/2)*m+9*cos((4*%pi*m+%pi)/2))
                        *x^3)
                        /24
                       +((1024*sin((4*%pi*m+%pi)/2)*m^5
                        +1280*sin((4*%pi*m+%pi)/2)*m^4
                        +1920*sin((4*%pi*m+%pi)/2)*m^3
                        +1120*sin((4*%pi*m+%pi)/2)*m^2
                        +356*sin((4*%pi*m+%pi)/2)*m+45*sin((4*%pi*m+%pi)/2))
                        *x^4)
                        /120
Is cos((4*%pi*m+%pi)/2) positive, negative or zero?

z;
1 Exit  taylim 0
1 Enter taylor 
[cos(4*m*asin(1/sqrt(1/(1/W1451)^2+1))+asin(1/sqrt(1/(1/W1451)^2+1)))/W1451,
 W1451,0,4]
1 Exit  taylor 
cos((4*%pi*m+%pi)/2)/W1451+(4*sin((4*%pi*m+%pi)/2)*m+sin((4*%pi*m+%pi)/2))
                          -((16*cos((4*%pi*m+%pi)/2)*m^2
                           +8*cos((4*%pi*m+%pi)/2)*m+cos((4*%pi*m+%pi)/2))
                           *W1451)
                           /2
                          -((64*sin((4*%pi*m+%pi)/2)*m^3
                           +48*sin((4*%pi*m+%pi)/2)*m^2
                           +20*sin((4*%pi*m+%pi)/2)*m+3*sin((4*%pi*m+%pi)/2))
                           *W1451^2)
                           /6
                          +((256*cos((4*%pi*m+%pi)/2)*m^4
                           +256*cos((4*%pi*m+%pi)/2)*m^3
                           +224*cos((4*%pi*m+%pi)/2)*m^2
                           +80*cos((4*%pi*m+%pi)/2)*m+9*cos((4*%pi*m+%pi)/2))
                           *W1451^3)
                           /24
                          +((1024*sin((4*%pi*m+%pi)/2)*m^5
                           +1280*sin((4*%pi*m+%pi)/2)*m^4
                           +1920*sin((4*%pi*m+%pi)/2)*m^3
                           +1120*sin((4*%pi*m+%pi)/2)*m^2
                           +356*sin((4*%pi*m+%pi)/2)*m
                           +45*sin((4*%pi*m+%pi)/2))
                           *W1451^4)
                           /120
1 Enter taylim 
[-cos(4*m*asin(1/sqrt(x^2+1))+asin(1/sqrt(x^2+1)))/x,x,zerob,think]
 1 Enter taylor [-cos(4*m*asin(1/sqrt(x^2+1))+asin(1/sqrt(x^2+1)))/x,x,0,4]
 1 Exit  taylor 
 (-cos((4*%pi*m+%pi)/2)/x)+((-4*sin((4*%pi*m+%pi)/2)*m)-sin((4*%pi*m+%pi)/2))
                          +((16*cos((4*%pi*m+%pi)/2)*m^2
                           +8*cos((4*%pi*m+%pi)/2)*m+cos((4*%pi*m+%pi)/2))
                           *x)
                           /2
                          +((64*sin((4*%pi*m+%pi)/2)*m^3
                           +48*sin((4*%pi*m+%pi)/2)*m^2
                           +20*sin((4*%pi*m+%pi)/2)*m+3*sin((4*%pi*m+%pi)/2))
                           *x^2)
                           /6
                          -((256*cos((4*%pi*m+%pi)/2)*m^4
                           +256*cos((4*%pi*m+%pi)/2)*m^3
                           +224*cos((4*%pi*m+%pi)/2)*m^2
                           +80*cos((4*%pi*m+%pi)/2)*m+9*cos((4*%pi*m+%pi)/2))
                           *x^3)
                           /24
                          -((1024*sin((4*%pi*m+%pi)/2)*m^5
                           +1280*sin((4*%pi*m+%pi)/2)*m^4
                           +1920*sin((4*%pi*m+%pi)/2)*m^3
                           +1120*sin((4*%pi*m+%pi)/2)*m^2
                           +356*sin((4*%pi*m+%pi)/2)*m
                           +45*sin((4*%pi*m+%pi)/2))
                           *x^4)
                           /120
1 Exit  taylim 0
(%o3) 0

limit seems to call TAYLIM twice, dunno why it decided it needed to do that.

rtoy commented 1 month ago

Imported from SourceForge on 2024-07-08 21:15:13 Created by davewittemorris on 2022-07-10 21:17:41 Original: https://sourceforge.net/p/maxima/bugs/4004/#70d0

Thanks for looking at this. Based on your last message, I found a much simpler example:

(%i1) assume(equal(a, 0));
(%o1) [equal(a,0)]
(%i2) limit(a/x + 1, x, 0);
(%o2) 0
rtoy commented 1 month ago

Imported from SourceForge on 2024-07-08 21:15:16 Created by robert_dodier on 2022-07-11 22:56:05 Original: https://sourceforge.net/p/maxima/bugs/4004/#13bf

Dave, thanks a lot for that example. I see that limit(a*x + 1, x, inf) has the same bug -- the result is 0 if you respond z to "Is a positive, negative or zero?".

Looks like the bug is in RATLIM (src/limit.lisp), which appears to be for computing limits for rational functions. A glance seems to show RATLIM has cases for less than zero and greater than zero, but assumes the result is 0 in the case that the sign is zero -- look for the calls to GETSIGNL and the code after those.

rtoy commented 1 month ago

Imported from SourceForge on 2024-07-08 21:15:20 Created by davewittemorris on 2022-07-13 23:17:07 Original: https://sourceforge.net/p/maxima/bugs/4004/#dd70

Thanks! I know very little lisp, but I am able to understand the gist of the RATLIM function: to calculate a limit x -> 0 (or x ->zeroa) of a rational function, we just need to know the coefficient of the lowest-order term in the numerator and denominator.

I think the best option would be to modify the locoef function to check whether the result is actually nonzero, and continue looking if it is not. Another alternative would be for RATLIM to check whether (locoef n g) and (locoef d g) are zero. If one of them is 0, then the corresponding term can be deleted from e, and RATLIM can be applied (recursively) to the modified expression.
Or RATLIM could search through the numerator (and denominator) for a term that is not 0, but I think that such a search should be delegated to locoef.

PS: It seems to me that there is another (unrelated) bug in RATLIM. Look at lines 3-4:

((eq val '$minf)
     (setq e (maxima-substitute (m^t -1 (m^t 'x -1)) var e)))

This is supposed to convert a limit approaching -inf to a limit approaching 0 from the right. The transformation for this is -1/x, so I think m^t -1 should be changed to m-. If this is indeed wrong, then I don't know why nobody has reported errors or incorrect results.

rtoy commented 1 month ago

Imported from SourceForge on 2024-07-08 21:15:23 Created by robert_dodier on 2022-07-14 16:35:47 Original: https://sourceforge.net/p/maxima/bugs/4004/#58eb

I think I've made progress on the original problem by substituting 0 for the expression which was asked about, and then calling RATLIM on the result. That gives correct results for the cases mentioned in this bug report (although for the original problem, limit needs a hint in the form of a direction plus or minus, since otherwise it can't tell the right and left limits are the same; need to look at that more closely). I'll post more later.

About the buggy code at lines 3 and 4, I agree that looks to be incorrect -- I think that should be (m* -1 (m^t 'x -1)). However, it appears that's never executed, so it should either be corrected or expunged. I'll open a separate ticket about that.

rtoy commented 1 month ago

Imported from SourceForge on 2024-07-08 21:15:27 Created by davewittemorris on 2022-07-14 18:36:03 Original: https://sourceforge.net/p/maxima/bugs/4004/#194a

I don't think it's sufficient to try to fix the problem at the point where RATLIM currently thinks the limit is 0 (although it would certainly be an improvement). For example, this limit should be 2 when a is 0, but RATLIM currently says that the limit is 1, whether there is an assumption or not (and will not ask any questions if you omit the assume):

(%i1) assume(equal(a, 0));
(%o1)                            [equal(a, 0)]
(%i2) limit((2*x + a)/(x + a), x, 0);
(%o2)                                  1

So it seems to me that it would be better to fix locoef.

rtoy commented 1 month ago

Imported from SourceForge on 2024-07-08 21:15:30 Created by macrakis on 2022-07-14 20:33:53 Original: https://sourceforge.net/p/maxima/bugs/4004/#a3ec

I've filed that example as https://sourceforge.net/p/maxima/bugs/4006/