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

Commits:

2 changed files:

Changes:

  • src/code/numbers.lisp
    ... ... @@ -593,180 +593,18 @@
    593 593 
     	    (build-ratio x y)
    594 594 
     	    (build-ratio (truncate x gcd) (truncate y gcd))))))
    595 595
    596 
    -
    597 
    -#+nil
    598 
    -(defun two-arg-/ (x y)
    599 
    -  (number-dispatch ((x number) (y number))
    600 
    -    (float-contagion / x y (ratio integer))
    601 
    -
    602 
    -    ((complex complex)
    603 
    -     (let* ((rx (realpart x))
    604 
    -	    (ix (imagpart x))
    605 
    -	    (ry (realpart y))
    606 
    -	    (iy (imagpart y)))
    607 
    -       (if (> (abs ry) (abs iy))
    608 
    -	   (let* ((r (/ iy ry))
    609 
    -		  (dn (+ ry (* r iy))))
    610 
    -	     (canonical-complex (/ (+ rx (* ix r)) dn)
    611 
    -				(/ (- ix (* rx r)) dn)))
    612 
    -	   (let* ((r (/ ry iy))
    613 
    -		  (dn (+ iy (* r ry))))
    614 
    -	     (canonical-complex (/ (+ (* rx r) ix) dn)
    615 
    -				(/ (- (* ix r) rx) dn))))))
    616 
    -    (((foreach integer ratio single-float double-float) complex)
    617 
    -     (let* ((ry (realpart y))
    618 
    -	    (iy (imagpart y)))
    619 
    -       (if (> (abs ry) (abs iy))
    620 
    -	   (let* ((r (/ iy ry))
    621 
    -		  (dn (* ry (+ 1 (* r r)))))
    622 
    -	     (canonical-complex (/ x dn)
    623 
    -				(/ (- (* x r)) dn)))
    624 
    -	   (let* ((r (/ ry iy))
    625 
    -		  (dn (* iy (+ 1 (* r r)))))
    626 
    -	     (canonical-complex (/ (* x r) dn)
    627 
    -				(/ (- x) dn))))))
    628 
    -    ((complex (or rational float))
    629 
    -     (canonical-complex (/ (realpart x) y)
    630 
    -			(/ (imagpart x) y)))
    631 
    -
    632 
    -    ((ratio ratio)
    633 
    -     (let* ((nx (numerator x))
    634 
    -	    (dx (denominator x))
    635 
    -	    (ny (numerator y))
    636 
    -	    (dy (denominator y))
    637 
    -	    (g1 (gcd nx ny))
    638 
    -	    (g2 (gcd dx dy)))
    639 
    -       (build-ratio (* (maybe-truncate nx g1) (maybe-truncate dy g2))
    640 
    -		    (* (maybe-truncate dx g2) (maybe-truncate ny g1)))))
    641 
    -
    642 
    -    ((integer integer)
    643 
    -     (integer-/-integer x y))
    644 
    -
    645 
    -    ((integer ratio)
    646 
    -     (if (zerop x)
    647 
    -	 0
    648 
    -	 (let* ((ny (numerator y))
    649 
    -		(dy (denominator y))
    650 
    -		(gcd (gcd x ny)))
    651 
    -	   (build-ratio (* (maybe-truncate x gcd) dy)
    652 
    -			(maybe-truncate ny gcd)))))
    653 
    -
    654 
    -    ((ratio integer)
    655 
    -     (let* ((nx (numerator x))
    656 
    -	    (gcd (gcd nx y)))
    657 
    -       (build-ratio (maybe-truncate nx gcd)
    658 
    -		    (* (maybe-truncate y gcd) (denominator x)))))))
    659 
    -
    660 
    -;; An implementation of Baudin and Smith's robust complex division.
    661 
    -;; This is a pretty straightforward translation of the original in
    662 
    -;; https://arxiv.org/pdf/1210.4539.
    663 
    -#+nil
    664 
    -(eval-when (:compile-toplevel :load-toplevel :execute)
    665 
    -(macrolet
    666 
    -    ((frob (type max-float min-float eps-float)
    667 
    -       (let ((name (symbolicate "CDIV-" type)))
    668 
    -	 `(progn
    669 
    -	    (let ((ov ,max-float)
    670 
    -		  (un ,min-float)
    671 
    -		  (eps ,eps-float))
    672 
    -	      (defun ,name (x y)
    673 
    -		(declare (type (complex ,type) x y))
    674 
    -		(labels
    675 
    -		    ((internal-compreal (a b c d r tt)
    676 
    -		       (declare (double-float a b c d r tt))
    677 
    -		       ;; Iteration 1:  compare r against DBL_MIN instead of 0
    678 
    -		       (cond ((/= r 0)
    679 
    -			      (let ((br (* b r)))
    680 
    -				(if (/= br 0)
    681 
    -				    (* (+ a br) tt)
    682 
    -				    (+ (* a tt)
    683 
    -				       (* (* b tt)
    684 
    -					  r)))))
    685 
    -			     (t
    686 
    -			      (* (+ a (* d (/ b c)))
    687 
    -				 tt))))
    688 
    -
    689 
    -		     (robust-subinternal (a b c d)
    690 
    -		       (declare (double-float a b c d))
    691 
    -		       (let* ((r (/ d c))
    692 
    -			      (tt (/ (+ c (* d r)))))
    693 
    -			 (let ((e (internal-compreal a b c d r tt))
    694 
    -			       (f (internal-compreal b (- a) c d r tt)))
    695 
    -			   #+nil(format t "subint: e f = ~A ~A~%" e f)
    696 
    -			   (values e
    697 
    -				   f))))
    698 
    -
    699 
    -		     (robust-internal (x y)
    700 
    -		       (declare (type (complex ,type) x y))
    701 
    -		       #+nil(format t "x = ~A, y = ~A~%" x y)
    702 
    -		       (let ((a (realpart x))
    703 
    -			     (b (imagpart x))
    704 
    -			     (c (realpart y))
    705 
    -			     (d (imagpart y)))
    706 
    -			 (cond
    707 
    -			   ((<= (abs d) (abs c))
    708 
    -			    (multiple-value-bind (e f)
    709 
    -				(robust-subinternal a b c d)
    710 
    -			      #+nil(format t "1: e f = ~A ~A~%" e f)
    711 
    -			      (complex e f)))
    712 
    -			   (t
    713 
    -			    (multiple-value-bind (e f)
    714 
    -				(robust-subinternal b a d c)
    715 
    -			      #+nil(format t "2: e f = ~A ~A~%" e f)
    716 
    -			      (complex e (- f))))))))
    717 
    -		  (let* ((a (realpart x))
    718 
    -			 (b (imagpart x))
    719 
    -			 (c (realpart y))
    720 
    -			 (d (imagpart y))
    721 
    -			 (ab (max (abs a) (abs b)))
    722 
    -			 (cd (max (abs c) (abs d)))
    723 
    -			 (b 2d0)
    724 
    -			 (s 1d0)
    725 
    -			 (be (/ b (* eps eps))))
    726 
    -		    (when (>= ab (/ ov 2))
    727 
    -		      (setf x (/ x 2)
    728 
    -			    s (* s 2)))
    729 
    -		    (when (>= cd (/ ov 2))
    730 
    -		      (setf y (/ y 2)
    731 
    -			    s (/ s 2)))
    732 
    -		    (when (<= ab (* un (/ b eps)))
    733 
    -		      (setf x (* x be)
    734 
    -			    s (/ s be)))
    735 
    -		    (when (<= cd (* un (/ b eps)))
    736 
    -		      (setf y (* y be)
    737 
    -			    s (* s be)))
    738 
    -		    ;; Iteration 2
    739 
    -		    #+nil
    740 
    -		    (when (< cd eps)
    741 
    -		      (setf x (/ x eps))
    742 
    -		      (setf y (/ y eps)))
    743 
    -
    744 
    -		    ;; Iteration 4
    745 
    -		    #+nil
    746 
    -		    (when (> cd (/ ov 2))
    747 
    -		      (setf x (/ x 2)
    748 
    -			    y (/ y 2)))
    749 
    -
    750 
    -		    (* s
    751 
    -		       (robust-internal x y))))))))))
    752 
    -  ;; The eps value for double-float is determined by looking at the
    753 
    -  ;; value of %eps in Scilab.
    754 
    -  (frob double-float most-positive-double-float least-positive-normalized-double-float
    755 
    -	(scale-float 1d0 -52))
    756 
    -  ;; The eps value here is a guess.
    757 
    -  (frob single-float most-positive-single-float least-positive-normalized-single-float
    758 
    -	(scale-float 1f0 -22)))
    759 
    -)
    760 
    -
    596 
    +;; An implementation of Baudin and Smith's robust complex division for
    597 
    +;; double-floats.  This is a pretty straightforward translation of the
    598 
    +;; original in https://arxiv.org/pdf/1210.4539.
    761 599 
     (let ((ov most-positive-double-float)
    762 600 
           (un least-positive-normalized-double-float)
    601 
    +      ;; This is the value of Scilab's %eps variable.
    763 602 
           (eps (scale-float 1d0 -52)))
    764 603 
       (defun cdiv-double-float (x y)
    765 604 
         (declare (type (complex double-float) x y))
    766 605 
         (labels
    767 606 
     	((internal-compreal (a b c d r tt)
    768 607 
     	   (declare (double-float a b c d r tt))
    769 
    -	   ;; Iteration 1:  compare r against DBL_MIN instead of 0
    770 608 
     	   (cond ((/= r 0)
    771 609 
     		  (let ((br (* b r)))
    772 610 
     		    (if (/= br 0)
    ... ... @@ -777,7 +615,6 @@
    777 615 
     		 (t
    778 616 
     		  (* (+ a (* d (/ b c)))
    779 617 
     		     tt))))
    780 
    -
    781 618 
     	 (robust-subinternal (a b c d)
    782 619 
     	   (declare (double-float a b c d))
    783 620 
     	   (let* ((r (/ d c))
    ... ... @@ -787,7 +624,6 @@
    787 624 
     	       #+nil(format t "subint: e f = ~A ~A~%" e f)
    788 625 
     	       (values e
    789 626 
     		       f))))
    790 
    -
    791 627 
     	 (robust-internal (x y)
    792 628 
     	   (declare (type (complex double-float) x y))
    793 629 
     	   #+nil(format t "x = ~A, y = ~A~%" x y)
    ... ... @@ -827,21 +663,11 @@
    827 663 
     	(when (<= cd (* un (/ b eps)))
    828 664 
     	  (setf y (* y be)
    829 665 
     		s (* s be)))
    830 
    -	;; Iteration 2
    831 
    -	#+nil
    832 
    -	(when (< cd eps)
    833 
    -	  (setf x (/ x eps))
    834 
    -	  (setf y (/ y eps)))
    835 
    -
    836 
    -	;; Iteration 4
    837 
    -	#+nil
    838 
    -	(when (> cd (/ ov 2))
    839 
    -	  (setf x (/ x 2)
    840 
    -		y (/ y 2)))
    841 
    -
    842 666 
     	(* s
    843 667 
     	   (robust-internal x y))))))
    844 668
    669 
    +;; Smith's algorithm for complex division for (complex single-float).
    670 
    +;; We convert the parts to double-floats before computing the result.
    845 671 
     (defun cdiv-single-float (x y)
    846 672 
       (declare (type (complex single-float) x y))
    847 673 
       (let ((a (float (realpart x) 1d0))
    ... ... @@ -877,12 +703,15 @@
    877 703 
          (cdiv-single-float (coerce x '(complex single-float))
    878 704 
     			y))
    879 705
    880 
    -    (((foreach integer ratio single-float double-float (complex rational) (complex single-float))
    706 
    +    (((foreach integer ratio single-float double-float (complex rational)
    707 
    +	       (complex single-float))
    881 708 
           (complex double-float))
    882 709 
          (cdiv-double-float (coerce x '(complex double-float))
    883 710 
     			y))
    884 711 
         (((complex double-double-float)
    885 
    -      (foreach (complex rational) (complex single-float) (complex double-float) (complex double-double-float)))
    712 
    +      (foreach (complex rational) (complex single-float) (complex double-float)
    713 
    +	       (complex double-double-float)))
    714 
    +     ;; We should do something better for double-double floats.
    886 715 
          (let ((a (realpart x))
    887 716 
     	   (b (imagpart x))
    888 717 
     	   (c (realpart y))
    ... ... @@ -935,6 +764,8 @@
    935 764 
     				(/ (- x) dn))))))
    936 765 
         (((complex rational)
    937 766 
           (complex rational))
    767 
    +     ;; We probably don't need to do Smith's algorithm for rationals.
    768 
    +     ;; A naive implementation of coplex division has no issues.
    938 769 
          (let ((a (realpart x))
    939 770 
     	   (b (imagpart x))
    940 771 
     	   (c (realpart y))
    ... ... @@ -952,7 +783,8 @@
    952 783 
     		     (f (/ (- b (* a r)) denom)))
    953 784 
     		(canonical-complex e f))))))
    954 785
    955 
    -    (((foreach (complex rational) (complex single-float) (complex double-float) (complex double-double-float))
    786 
    +    (((foreach (complex rational) (complex single-float) (complex double-float)
    787 
    +	       (complex double-double-float))
    956 788 
           (or rational float))
    957 789 
          (canonical-complex (/ (realpart x) y)
    958 790 
     			(/ (imagpart x) y)))
    ... ... @@ -962,108 +794,6 @@
    962 794 
          (cdiv-double-float (coerce x '(complex double-float))
    963 795 
     			(coerce y '(complex double-float))))
    964 796
    965 
    -    #+nil
    966 
    -    (((foreach integer ratio single-float double-float) complex)
    967 
    -     (let* ((ry (realpart y))
    968 
    -	    (iy (imagpart y)))
    969 
    -       (if (> (abs ry) (abs iy))
    970 
    -	   (let* ((r (/ iy ry))
    971 
    -		  (dn (* ry (+ 1 (* r r)))))
    972 
    -	     (canonical-complex (/ x dn)
    973 
    -				(/ (- (* x r)) dn)))
    974 
    -	   (let* ((r (/ ry iy))
    975 
    -		  (dn (* iy (+ 1 (* r r)))))
    976 
    -	     (canonical-complex (/ (* x r) dn)
    977 
    -				(/ (- x) dn))))))
    978 
    -    #+nil
    979 
    -    ((complex (or rational float))
    980 
    -     (canonical-complex (/ (realpart x) y)
    981 
    -			(/ (imagpart x) y)))
    982 
    -
    983 
    -    ((ratio ratio)
    984 
    -     (let* ((nx (numerator x))
    985 
    -	    (dx (denominator x))
    986 
    -	    (ny (numerator y))
    987 
    -	    (dy (denominator y))
    988 
    -	    (g1 (gcd nx ny))
    989 
    -	    (g2 (gcd dx dy)))
    990 
    -       (build-ratio (* (maybe-truncate nx g1) (maybe-truncate dy g2))
    991 
    -		    (* (maybe-truncate dx g2) (maybe-truncate ny g1)))))
    992 
    -
    993 
    -    ((integer integer)
    994 
    -     (integer-/-integer x y))
    995 
    -
    996 
    -    ((integer ratio)
    997 
    -     (if (zerop x)
    998 
    -	 0
    999 
    -	 (let* ((ny (numerator y))
    1000 
    -		(dy (denominator y))
    1001 
    -		(gcd (gcd x ny)))
    1002 
    -	   (build-ratio (* (maybe-truncate x gcd) dy)
    1003 
    -			(maybe-truncate ny gcd)))))
    1004 
    -
    1005 
    -    ((ratio integer)
    1006 
    -     (let* ((nx (numerator x))
    1007 
    -	    (gcd (gcd nx y)))
    1008 
    -       (build-ratio (maybe-truncate nx gcd)
    1009 
    -		    (* (maybe-truncate y gcd) (denominator x)))))))
    1010 
    -
    1011 
    -
    1012 
    -#+nil
    1013 
    -(defun two-arg-/ (x y)
    1014 
    -  (number-dispatch ((x number) (y number))
    1015 
    -    (float-contagion / x y (ratio integer))
    1016 
    -
    1017 
    -    (((complex double-float) (complex double-float))
    1018 
    -     (cdiv-double-float x y))
    1019 
    -    (((complex double-float) (complex single-float))
    1020 
    -     (cdiv-double-float x (coerce y '(complex double-float))))
    1021 
    -
    1022 
    -    (((complex single-float) (complex double-float))
    1023 
    -     (cdiv-double-float (coerce y '(complex double-float))
    1024 
    -			y))
    1025 
    -    (((complex single-float) (complex single-float))
    1026 
    -     (cdiv-single-float x y))
    1027 
    -
    1028 
    -    #+nil
    1029 
    -    (((foreach (complex single-float) (complex double-float))
    1030 
    -      (complex double-float))
    1031 
    -     (cdiv-double-float (complex (float (realpart x) 1d0)
    1032 
    -				 (float (imagpart x) 1d0))
    1033 
    -			y))
    1034 
    -    (((foreach (complex rational) (complex double-double-float))
    1035 
    -      (foreach (complex rational) (complex double-double-float)))
    1036 
    -     (let* ((rx (realpart x))
    1037 
    -	    (ix (imagpart x))
    1038 
    -	    (ry (realpart y))
    1039 
    -	    (iy (imagpart y)))
    1040 
    -       (if (> (abs ry) (abs iy))
    1041 
    -	   (let* ((r (/ iy ry))
    1042 
    -		  (dn (+ ry (* r iy))))
    1043 
    -	     (canonical-complex (/ (+ rx (* ix r)) dn)
    1044 
    -				(/ (- ix (* rx r)) dn)))
    1045 
    -	   (let* ((r (/ ry iy))
    1046 
    -		  (dn (+ iy (* r ry))))
    1047 
    -	     (canonical-complex (/ (+ (* rx r) ix) dn)
    1048 
    -				(/ (- (* ix r) rx) dn))))))
    1049 
    -    (((foreach integer ratio single-float double-float)
    1050 
    -      (foreach (complex rational) (complex double-double-float)))
    1051 
    -     (let* ((ry (realpart y))
    1052 
    -	    (iy (imagpart y)))
    1053 
    -       (if (> (abs ry) (abs iy))
    1054 
    -	   (let* ((r (/ iy ry))
    1055 
    -		  (dn (* ry (+ 1 (* r r)))))
    1056 
    -	     (canonical-complex (/ x dn)
    1057 
    -				(/ (- (* x r)) dn)))
    1058 
    -	   (let* ((r (/ ry iy))
    1059 
    -		  (dn (* iy (+ 1 (* r r)))))
    1060 
    -	     (canonical-complex (/ (* x r) dn)
    1061 
    -				(/ (- x) dn))))))
    1062 
    -    (((foreach (complex rational) (complex double-double-float))
    1063 
    -      (or rational float))
    1064 
    -     (canonical-complex (/ (realpart x) y)
    1065 
    -			(/ (imagpart x) y)))
    1066 
    -
    1067 797 
         ((ratio ratio)
    1068 798 
          (let* ((nx (numerator x))
    1069 799 
     	    (dx (denominator x))
    ... ... @@ -1092,7 +822,6 @@
    1092 822 
            (build-ratio (maybe-truncate nx gcd)
    1093 823 
     		    (* (maybe-truncate y gcd) (denominator x)))))))
    1094 824
    1095 
    -
    1096 825 
     (defun %negate (n)
    1097 826 
       (number-dispatch ((n number))
    1098 827 
         (((foreach fixnum single-float double-float #+long-float long-float))

  • src/compiler/float-tran.lisp
    ... ... @@ -1804,119 +1804,6 @@
    1804 1804 
       #+double-double
    1805 1805 
       (frob double-double-float))
    1806 1806
    1807 
    -#+(and nil sse2 complex-fp-vops)
    1808 
    -(macrolet
    1809 
    -    ((frob (type one)
    1810 
    -       `(deftransform / ((x y) (,type ,type) *
    1811 
    -			 :policy (> speed space))
    1812 
    -	  ;; Divide a complex by a complex
    1813 
    -
    1814 
    -	  ;; Here's how we do a complex division
    1815 
    -	  ;;
    1816 
    -	  ;; Compute (xr + i*xi)/(yr + i*yi)
    1817 
    -	  ;;
    1818 
    -	  ;; Assume |yi| < |yr|.  Then
    1819 
    -	  ;;
    1820 
    -	  ;; (xr + i*xi)      (xr + i*xi)
    1821 
    -	  ;; ----------- = -----------------
    1822 
    -	  ;; (yr + i*yi)   yr*(1 + i*(yi/yr))
    1823 
    -	  ;;
    1824 
    -	  ;;               (xr + i*xi)*(1 - i*(yi/yr))
    1825 
    -	  ;;             = ---------------------------
    1826 
    -	  ;;                   yr*(1 + (yi/yr)^2)
    1827 
    -	  ;;
    1828 
    -	  ;;               (xr + i*xi)*(1 - i*(yi/yr))
    1829 
    -	  ;;             = ---------------------------
    1830 
    -	  ;;                   yr + (yi/yr)*yi
    1831 
    -	  ;;
    1832 
    -	  ;; This allows us to use a fast complex multiply followed by
    1833 
    -	  ;; a real division.
    1834 
    -	  '(let* ((ry (realpart y))
    1835 
    -		  (iy (imagpart y)))
    1836 
    -	    (if (> (abs ry) (abs iy))
    1837 
    -		(let* ((r (/ iy ry))
    1838 
    -		       (dn (+ ry (* r iy))))
    1839 
    -		  (/ (* x (complex ,one r))
    1840 
    -		     dn))
    1841 
    -		(let* ((r (/ ry iy))
    1842 
    -		       (dn (+ iy (* r ry))))
    1843 
    -		  (/ (* x (complex r ,(- one)))
    1844 
    -		     dn)))))))
    1845 
    -  (frob (complex single-float) 1f0)
    1846 
    -  (frob (complex double-float) 1d0))
    1847 
    -
    1848 
    -#+(and nil sse2 complex-fp-vops)
    1849 
    -(macrolet
    1850 
    -    ((frob (type one)
    1851 
    -       `(deftransform / ((x y) (,type ,type) *
    1852 
    -			 :policy (> speed space))
    1853 
    -	  ;; Divide a complex by a complex
    1854 
    -
    1855 
    -	  ;; Here's how we do a complex division
    1856 
    -	  ;;
    1857 
    -	  ;; Compute (xr + i*xi)/(yr + i*yi)
    1858 
    -	  ;;
    1859 
    -	  ;; Assume |xi| < |xr| and |yi| < |yr|.  Then
    1860 
    -	  ;;
    1861 
    -	  ;; (xr + i*xi)   xr*(1 + i*(xi/xr))
    1862 
    -	  ;; ----------- = -----------------
    1863 
    -	  ;; (yr + i*yi)   yr*(1 + i*(yi/yr))
    1864 
    -	  ;;
    1865 
    -          ;;               xr*(1+i(xi/xr))*(1-i*(yi/yr))
    1866 
    -	  ;;             = -----------------------------
    1867 
    -          ;;                    yr*(1+(yi/yr)^2)
    1868 
    -	  ;;
    1869 
    -          ;;               xr*(1+i(xi/xr))*(1-i*(yi/yr))
    1870 
    -          ;;             = -----------------------------
    1871 
    -          ;;                    yr + (yi/yr)*yi
    1872 
    -          ;;
    1873 
    -	  ;; This allows us to use a fast complex multiply followed by
    1874 
    -	  ;; a real division.
    1875 
    -	  '(let* ((rx (realpart x))
    1876 
    -		  (ix (imagpart x))
    1877 
    -		  (ry (realpart y))
    1878 
    -		  (iy (imagpart y)))
    1879 
    -	    (if (>= (abs rx) (abs ix))
    1880 
    -		(if (>= (abs ry) (abs iy))
    1881 
    -		    ;; |xr| >= |xi| and |yr| >= |yi|
    1882 
    -		    ;;
    1883 
    -		    ;;     xr*(1+i*(xi/xr))*(1-i*(yi/yr))
    1884 
    -		    ;; w = ------------------------------
    1885 
    -		    ;;        yr + (yi/yr)*yi
    1886 
    -		    ;;
    1887 
    -		    ;; (1+i*p)*(1-i*q) = 1+p*q + i*(p-q)
    1888 
    -		    (let* ((x-div (/ ix rx))
    1889 
    -			   (y-div (/ iy ry))
    1890 
    -			   (dn (+ ry (* y-div iy)))
    1891 
    -			   (factor (/ xr dn)))
    1892 
    -		      (* factor
    1893 
    -			 (complex (1+ (* x-div y-div))
    1894 
    -				  (- x-div y-div))))
    1895 
    -		    ;; |xr| >= |xi| and |yr| < |yi|
    1896 
    -		    ;;
    1897 
    -		    ;;     xr*(1+i*(xi/xr))*(yr/yi-i)
    1898 
    -		    ;; w = ------------------------------
    1899 
    -		    ;;         yi*((yr/yi)^2+1)
    1900 
    -		    ;;
    1901 
    -		    ;;     xr*(1+i*(xi/xr))*(yr/yi-i)
    1902 
    -		    ;; w = ------------------------------
    1903 
    -		    ;;         (yr/yi)*yi + yi
    1904 
    -		    ;;
    1905 
    -		    ;; (1+i*p)*(q-i) = p+q + i*(p*q-1)
    1906 
    -		    (let* ((x-div (/ ix rx))
    1907 
    -			   (y-div (/ ry iy))
    1908 
    -			   (dn (+ yi (* y-div iy)))
    1909 
    -			   (factor (/ xr dn)))
    1910 
    -		      (* factor
    1911 
    -			 (complex (+ x-div y-div)
    1912 
    -				  (1- (* x-div y-div))))))
    1913 
    -		    (let* ((r (/ ry iy))
    1914 
    -			   (dn (+ iy (* r ry))))
    1915 
    -		      (/ (* x (complex r ,(- one)))
    1916 
    -			 dn)))))))
    1917 
    -  (frob (complex single-float) 1f0)
    1918 
    -  (frob (complex double-float) 1d0))
    1919 
    -
    1920 1807 
     (macrolet
    1921 1808 
         ((frob (type name)
    1922 1809 
            `(deftransform / ((x y) ((complex ,type) (complex ,type)) *)

    • —
    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/-/commit/bb5c9e106d89594a1bc3020107ac54a006dac70d at 1765409534