[Maxima-commits] CVS: maxima/src rpart.lisp,1.32,1.33

[Raymond T. <rt...@us...>]

Update of /cvsroot/maxima/maxima/src
In directory sfp-cvsdas-1.v30.ch3.sourceforge.com:/tmp/cvs-serv21744/src

Modified Files:
	rpart.lisp 
Log Message:
Updaate sprecip to work around the fact that some Lisp's will give a
floating point error when computing (/ (complex x y)).  Use the
workaround for all Lisp except Clisp, CMUCL, and SBCL.


Index: rpart.lisp
===================================================================
RCS file: /cvsroot/maxima/maxima/src/rpart.lisp,v
retrieving revision 1.32
retrieving revision 1.33
diff -u -d -r1.32 -r1.33
--- rpart.lisp	22 Mar 2010 23:43:46 -0000	1.32
+++ rpart.lisp	31 Mar 2010 03:39:12 -0000	1.33
@@ -510,15 +510,38 @@
 ;; rationals.  Return a cons of the real and imaginary part of the
 ;; result.  We count on the underlying Lisp to be able to compute (/
 ;; (complex x y)) accurately and without unnecessary overflow or
-;; underflow..  If not, complain to your Lisp vendor.
+;; underflow..  If not, complain to your Lisp vendor.  Well, it seems
+;; that Clisp, CMUCL, and SBCL do a nice job.  But others apparently
+;; do not.  (I tested ecl 9.12.3 and ccl 1.4, which both fail.)
+;; Workaround those deficiencies.
 (defun sprecip (sp)
+  
   (destructuring-bind (x . y)
       sp
+    #+(or clisp cmu sbcl)
     (let* ((x (bigfloat:to x))
 	   (y (bigfloat:to y))
 	   (q (bigfloat:/ (bigfloat:complex x y))))
       (cons (to (bigfloat:realpart q))
-	    (to (bigfloat:imagpart q))))))
+	    (to (bigfloat:imagpart q))))
+    #-(or clisp cmu sbcl)
+    (let ((x (bigfloat:to x))
+	  (y (bigfloat:to y)))
+      ;; 1/(x+%i*y).
+      ;;
+      ;; Assume abs(x) > abs(y).  Then 1/(x+%i*y) = 1/x/(1+%i*r) where
+      ;; r = y/x.  Thus 1/(x+%i*y) = (1-%i*r)/(x*(1+r*r))
+      ;;
+      ;; The case for abs(x) <= abs(y) is similar.
+      (if (> (bigfloat:abs x) (bigfloat:abs y))
+	  (let* ((r (bigfloat:/ y x))
+		 (dn (bigfloat:* x (bigfloat:+ 1 (bigfloat:* r r)))))
+	    (cons (to (bigfloat:/ dn))
+		  (to (bigfloat:/ (bigfloat:- r) dn))))
+	  (let* ((r (bigfloat:/ x y))
+		 (dn (bigfloat:* y (bigfloat:+ 1 (bigfloat:* r r)))))
+	    (cons (to (bigfloat:/ x dn))
+		  (to (bigfloat:/ (bigfloat:- x) dn))))))))