Raymond Toy pushed to branch issue-456-more-accurate-complex-div at cmucl / cmucl

Commits:

2 changed files:

Changes:

  • src/code/numbers.lisp
    ... ... @@ -605,42 +605,85 @@
    605 605 
         (labels
    606 606 
     	((internal-compreal (a b c d r tt)
    607 607 
     	   (declare (double-float a b c d r tt))
    608 
    +	   ;; Compute the real part of the complex division
    609 
    +	   ;; (a+ib)/(c+id), assuming |c| <= |d|.  r = d/c and tt = 1/(c+d*r).
    610 
    +	   ;;
    611 
    +	   ;; The realpart is (a*c+b*d)/(c^2+d^2).
    612 
    +	   ;;
    613 
    +	   ;;   c^2+d^2 = c*(c+d*(d/c)) = c*(c+d*r)
    614 
    +	   ;;
    615 
    +	   ;; Then
    616 
    +	   ;;
    617 
    +	   ;;   (a*c+b*d)/(c^2+d^2) = (a*c+b*d)/(c*(c+d*r))
    618 
    +	   ;;                       = (a + b*d/c)/(c+d*r)
    619 
    +	   ;;                       = (a + b*r)/(c + d*r).
    620 
    +	   ;;
    621 
    +	   ;; Thus tt = (c + d*r).
    622 
    +	   #+nil
    623 
    +	   (progn
    624 
    +	     (format t "a,b,c,d = ~A ~A ~A ~A~%" a b c d)
    625 
    +	     (format t "  r, tt = ~A ~A~%" r tt))
    608 626 
     	   (cond ((/= r 0)
    609 627 
     		  (let ((br (* b r)))
    628 
    +		    #+nil
    629 
    +		    (format t "br = ~A~%" br)
    610 630 
     		    (if (/= br 0)
    611 
    -			(* (+ a br) tt)
    612 
    -			(+ (* a tt)
    613 
    -			   (* (* b tt)
    631 
    +			(/ (+ a br) tt)
    632 
    +			;; b*r underflows.  Instead, compute
    633 
    +			;;
    634 
    +			;; (a + b*r)*tt = a*tt + b*tt*r
    635 
    +			;;              = a*tt + (b*tt)*r
    636 
    +			;; (a + b*r)/tt = a/tt + b/tt*r
    637 
    +			;;              = a*tt + (b*tt)*r
    638 
    +			(+ (/ a tt)
    639 
    +			   (* (/ b tt)
    614 640 
     			      r)))))
    615 641 
     		 (t
    616 
    -		  (* (+ a (* d (/ b c)))
    642 
    +		  ;; r = 0 so d is very tiny compared to c.
    643 
    +		  ;;
    644 
    +		  ;; (a + b*r)/tt = (a + b*(d/c))/tt
    645 
    +		  #+nil
    646 
    +		  (progn
    647 
    +		    (format t "r = 0~%")
    648 
    +		    (format t "a*tt = ~A~%" (* a tt)))
    649 
    +		  (/ (+ a (* d (/ b c)))
    617 650 
     		     tt))))
    618 651 
     	 (robust-subinternal (a b c d)
    619 652 
     	   (declare (double-float a b c d))
    620 653 
     	   (let* ((r (/ d c))
    621 
    -		  (tt (/ (+ c (* d r)))))
    654 
    +		  (tt (+ c (* d r))))
    655 
    +	     ;; e is the real part and f is the imaginary part.  We
    656 
    +	     ;; can use internal-compreal for the imaginary part by
    657 
    +	     ;; noticing that the imaginary part of (a+i*b)/(c+i*d) is
    658 
    +	     ;; the same as the real part of (b-i*a)/(c+i*d).
    622 659 
     	     (let ((e (internal-compreal a b c d r tt))
    623 660 
     		   (f (internal-compreal b (- a) c d r tt)))
    624 
    -	       #+nil(format t "subint: e f = ~A ~A~%" e f)
    625 661 
     	       (values e
    626 662 
     		       f))))
    627 663 
     	 (robust-internal (x y)
    628 664 
     	   (declare (type (complex double-float) x y))
    629 
    -	   #+nil(format t "x = ~A, y = ~A~%" x y)
    630 665 
     	   (let ((a (realpart x))
    631 666 
     		 (b (imagpart x))
    632 667 
     		 (c (realpart y))
    633 668 
     		 (d (imagpart y)))
    634 669 
     	     (cond
    635 670 
     	       ((<= (abs d) (abs c))
    671 
    +		;; |d| <= |c|, so we can use robust-subinternal to
    672 
    +		;; perform the division.
    636 673 
     		(multiple-value-bind (e f)
    637 674 
     		    (robust-subinternal a b c d)
    638 
    -		  #+nil(format t "1: e f = ~A ~A~%" e f)
    639 675 
     		  (complex e f)))
    640 676 
     	       (t
    677 
    +		;; |d| > |c|.  So, instead compute
    678 
    +		;;
    679 
    +		;;   (b + i*a)/(d + i*c) = ((b*d+a*c) + (a*d-b*c)*i)/(d^2+c^2)
    680 
    +		;;
    681 
    +		;; Compare this to (a+i*b)/(c+i*d) and we see that
    682 
    +		;; realpart of the former is the same, but the
    683 
    +		;; imagpart of the former is the negative of the
    684 
    +		;; desired division.
    641 685 
     		(multiple-value-bind (e f)
    642 686 
     		    (robust-subinternal b a d c)
    643 
    -		  #+nil(format t "2: e f = ~A ~A~%" e f)
    644 687 
     		  (complex e (- f))))))))
    645 688 
           (let* ((a (realpart x))
    646 689 
     	     (b (imagpart x))
    ... ... @@ -648,19 +691,23 @@
    648 691 
     	     (d (imagpart y))
    649 692 
     	     (ab (max (abs a) (abs b)))
    650 693 
     	     (cd (max (abs c) (abs d)))
    651 
    -	     (b 2d0)
    694 
    +	     (bb 2d0)
    652 695 
     	     (s 1d0)
    653 
    -	     (be (/ b (* eps eps))))
    696 
    +	     (be (/ bb (* eps eps))))
    697 
    +	;; If a or b is big, scale down a and b.
    654 698 
     	(when (>= ab (/ ov 2))
    655 699 
     	  (setf x (/ x 2)
    656 700 
     		s (* s 2)))
    701 
    +	;; If c or d is big, scale down c and d.
    657 702 
     	(when (>= cd (/ ov 2))
    658 703 
     	  (setf y (/ y 2)
    659 704 
     		s (/ s 2)))
    660 
    -	(when (<= ab (* un (/ b eps)))
    705 
    +	;; If a or b is tiny, scale up a and b.
    706 
    +	(when (<= ab (* un (/ bb eps)))
    661 707 
     	  (setf x (* x be)
    662 708 
     		s (/ s be)))
    663 
    -	(when (<= cd (* un (/ b eps)))
    709 
    +	;; If c or d is tiny, scale up c and d.
    710 
    +	(when (<= cd (* un (/ bb eps)))
    664 711 
     	  (setf y (* y be)
    665 712 
     		s (* s be)))
    666 713 
     	(* s

  • tests/float.lisp
    ... ... @@ -388,7 +388,7 @@
    388 388 
        (list (complex (scale-float 1d0 -1074) (scale-float 1d0 -1074))
    389 389 
     	 (complex (scale-float 1d0 -1073) (scale-float 1d0 -1074))
    390 390 
     	 (complex 0.6d0 0.2d0)
    391 
    -	 52 least-positive-double-float)
    391 
    +	 53 least-positive-double-float)
    392 392 
        ;; 9
    393 393 
        (list (complex (scale-float 1d0 1015) (scale-float 1d0 -989))
    394 394 
     	 (complex (scale-float 1d0 1023) (scale-float 1d0 1023))
    ... ... @@ -399,11 +399,18 @@
    399 399 
     	 (complex (scale-float 1d0 -343) (scale-float 1d0 -798))
    400 400 
     	 (complex 1.02951151789360578d-84 6.97145987515076231d-220)
    401 401 
     	 53 least-positive-double-float)
    402 
    +   ;; 11
    402 403 
        ;; From Maxima
    403 404 
        (list #c(5.43d-10 1.13d-100)
    404 405 
     	 #c(1.2d-311 5.7d-312)
    405 406 
     	 #c(3.691993880674614517999740937026568563794896024143749539711267954d301
    406 407 
     	    -1.753697093319947872394996242210428954266103103602859195409591583d301)
    408 
    +	 52 least-positive-double-float)
    409 
    +   ;; 12
    410 
    +   ;; Found by ansi tests. z/z should be exactly 1.
    411 
    +   (list #c(1.565640716292489d19 0.0d0)
    412 
    +	 #c(1.565640716292489d19 0.0d0)
    413 
    +	 #c(1d0 0)
    407 414 
     	 53 least-positive-double-float)
    408 415 
        ))
    409 416
    ... ... @@ -494,10 +501,10 @@
    494 501
    495 502 
     (define-test complex-division.single
    496 503 
         (:tag :issues)
    497 
    -  (let ((x #c(1 2))
    498 
    -	(y (complex (expt 2 127) (expt 2 127)))
    499 
    -	(expected (coerce (/ x y)
    500 
    -			  '(complex single-float))))
    504 
    +  (let* ((x #c(1 2))
    505 
    +	 (y (complex (expt 2 127) (expt 2 127)))
    506 
    +	 (expected (coerce (/ x y)
    507 
    +			   '(complex single-float))))
    501 508 
         ;; A naive implementation of complex division would cause an
    502 509 
         ;; overflow in computing the denominator.
    503 510 
         (assert-equal expected

    • —
    View it on GitLab.
    You're receiving this email because of your account on gitlab.common-lisp.net. Manage all notifications · Help Notification message regarding https://gitlab.common-lisp.net/cmucl/cmucl/-/compare/bb5c9e106d89594a1bc3020107ac54a006dac70d...1aa7c11f6cf9da0be58523e74a66692564898a95 at 1765424010