[Maxima] Re: [Maxima-bugs] [ maxima-Bugs-1281737 ] limit(atan(x)/(1/exp(1)-exp(-(1+x)^2)),x,inf,plus) - wrong

Raymond Toy raymond.toy@ericsson.com
Thu, 15 Sep 2005 09:30:49 -0400


>>>>> ">" == SourceForge net <SourceForge.net> writes:

    >> Bugs item #1281737, was opened at 2005-09-04 15:30

[snip]
    >> Summary: limit(atan(x)/(1/exp(1)-exp(-(1+x)^2)),x,inf,plus) - wrong

    >> Initial Comment:

    >> (%i1) limit(atan(x)/(1/exp(1)-exp(-(1+x)^2)),x,inf,plus);
    >>                                 %e %i log(- 1)
    >> (%o1)                           --------------
    >>                                       2
    >> Correct result is  %e*%pi/2

[snip]

    >> Comment By: Raymond Toy (rtoy)
    >> Date: 2005-09-14 17:36

    >> Message:
    >> Logged In: YES 
    >> user_id=28849

    >> I think this is very likely caused by the fact that maxima
    >> thinks

    >> limit(log((1-%i*x)/(1+%i*x)),x,inf,plus)

    >> is log(-1) = %i*%pi/2.  The correct answer is -%i*%pi/2,
    >> because argument of the log is always in the third quadrant
    >> (x < 0, y < 0).

    >> I think it's because maxima thinks limit log(x) = log(limit
    >> x) in this case since the limit of the arg is a finite number.

I've looked a bit more at this.  A simpler example is 

limit(log(-1-%i*x),x,0,plus);

Maxima returns log(-1) = %i*%pi.  However, the correct answer is
-%i*%pi because we're always in the third quadrant.

Here is a potential patch for this.  What happens is we add %log to
one of the special cases in simplimit (clause 3).  This allows mlogp
case to happen.  We modify simplimln very simply by checking to see if
the imagpart is identically zero or not.  If it is zero, we use the
original simplimln.  If not, we take the limit of both the real part
and the imaginary part of the log expression and combine the results.

With this fix we get the expected limits for the original bug report,
and the two simplified cases.  All of the tests pass, but that's
hardly conclusive.

Perhaps someone who knows more about limit could comment on this.

Ray

Index: limit.lisp
===================================================================
RCS file: /cvsroot/maxima/maxima/src/limit.lisp,v
retrieving revision 1.14
diff -u -r1.14 limit.lisp
--- limit.lisp	9 Aug 2005 03:48:23 -0000	1.14
+++ limit.lisp	15 Sep 2005 13:28:04 -0000
@@ -1549,7 +1549,7 @@
     ((eq var exp) val)
     ((or (atom exp) (mnump exp)) exp)
     ((and (not (infinityp val))
-	  (not (amongl '(%sin %cos %atanh %cosh %sinh %tanh mfactorial)
+	  (not (amongl '(%sin %cos %atanh %cosh %sinh %tanh mfactorial %log)
 		       exp))
 	  (not (inf-typep exp))
 	  (simplimsubst val exp)))
@@ -2408,6 +2408,7 @@
 	   (t (return ($radcan (ridofab (subin val e))))))
      (return (simplimtimes (list n1 d1)))))
 
+#+nil
 (defun simplimln (arg)
   (let* ((arglim (limit arg var val 'think))
 	 (real-lim (ridofab arglim)))
@@ -2426,6 +2427,36 @@
 		 (if (equal dir 1) '$zeroa 0)))
 	      (t (simplify `((%log) ,real-lim)))))))
 
+(defun simplimln (arg)
+  ;; We need to be careful with log because of the branch cut on the
+  ;; negative real axis.  So we look at the imagpart of the log.  If
+  ;; it's not identically zero, we compute the limit of the real and
+  ;; imaginary parts and combine them.  Otherwise, we can use the
+  ;; original method for real limits.
+  (let* ((log-form `((%log) ,arg))
+	 (rp ($realpart log-form))
+	 (ip (simplify ($imagpart log-form))))
+    (cond ((and (numberp ip) (zerop ip))
+	   (let* ((arglim (limit arg var val 'think))
+		  (real-lim (ridofab arglim)))
+	     (if (=0 real-lim)
+		 (cond ((eq arglim '$zeroa)  '$minf)
+		       ((eq arglim '$zerob)  '$infinity)
+		       (t (let ((dir (behavior arg var val)))
+			    (cond ((equal dir 1) '$minf)
+				  ((equal dir -1) '$infinity)
+				  (t (throw 'limit t))))))
+		 (cond ((eq arglim '$inf) '$inf)
+		       ((memq arglim '($minf $infinity)) '$infinity)
+		       ((memq arglim '($ind $und)) '$und)
+		       ((equal arglim 1)
+			(let ((dir (behavior arg var val)))
+			  (if (equal dir 1) '$zeroa 0)))
+		       (t (simplify `((%log) ,real-lim)))))))
+	  (t
+	   (add (limit rp var val 'think)
+		(mul '$%i (limit ip var val 'think)))))))
+
 (defun simplimfact (exp var val arg)
   (cond ((eq arg '$inf) '$inf)
 	((memq arg '($minf $infinity $und $ind)) '$und)