#4004 a cosine of arcsin limit that is evaluated incorrectly

[Dave Morris]

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"

[Robert Dodier]

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"

[Robert Dodier]

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.

[Robert Dodier]

A somewhat simpler example. I've substituted n for 4*m + 1. I call ratdisrep in between taylor and limit because that's what TAYLIM does.

(%i1) assume (x > 0);
(%o1)                               [x > 0]
(%i2) taylor (cos(n*asin(1/sqrt(x^2+1)))/x, x, 0, 1);
                 %pi n                        %pi n   2
             cos(-----)                  (cos(-----) n ) x
                   2          %pi n             2
(%o2)/T/     ---------- + sin(-----) n - ----------------- + . . .
                 x              2                2
(%i3) ratdisrep (%);
                    2     %pi n           %pi n
                   n  cos(-----) x    cos(-----)
                            2               2            %pi n
(%o3)           (- ---------------) + ---------- + n sin(-----)
                          2               x                2
(%i4) limit (%, x, 0, plus);
       %pi n
Is cos(-----) positive, negative or zero?
         2

z;
(%o4)                                  0

[Dave Morris]

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

[Robert Dodier]

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.

[Dave Morris]

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.

[Robert Dodier]

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.

[Dave Morris]

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.

[Stavros Macrakis]

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

[Robert Dodier]

Dave, excellent observation about LOCOEF. You're right, the substitution approach doesn't get the right result for the example you showed above. I'll take another look at it -- @macrakis can I ask you to look at LOCOEF as Dave M. was suggesting yesterday (July 13)?

[Barton Willis]

We need to follow through on several observations in this thread.

But at the risk of derailing following through, a simple way to make this bug disappear is to insert trigexpand into gruntz1 just before it calls limitinf.

There might be a few more tweaks: insert code into gruntz1 for assumptions about the limit variable, for one. The gruntz code changes the limit point, so the limit code that does assumptions on the limit variable are out of date.

Here I've pasted in a call to trigexpand into gruntz1:

(%i1) load("limit.lisp")$
(%i2)  display2d : false$
(%i3)  declare(m, integer)$

(%i4) assume(equal(cos((4 * %pi * m + %pi)/2), 0))$

(%i5) limit(cos((4*m + 1) * asin(1/sqrt(x^2 + 1)))/abs(x), x, 0);
(%o5) 4*m+1

After inserting trigexpand into gruntz1, the function gruntz is not able to compute this limit, but it does allow gruntz1 to find some other limit that allows limit to correctly return 4m + 1.

Inserting trigexpand into gruntz1 does not fix other bugs mentioned in this thread:

(%i6) assume(equal(a, 0));

(%o6) [equal(a,0)]
(%i7) limit((2*x + a)/(x + a), x, 0);