From 35a0dfcafecaa3c9f3a8dc852be02865bc265332 Mon Sep 17 00:00:00 2001
From: Matthew Flatt <mflatt@racket-lang.org>
Date: Fri, 28 Jun 2019 09:01:06 -0600
Subject: [PATCH] improve precision of complex-number divide

original commit: 4c9a7f6abb1258158d48fcdb656de300902cf3c7
---
 mats/5_3.ms  |  8 ++++++++
 s/5_3.ss     | 35 ++++++++++++++++++++++++++++++++---
 s/library.ss | 32 ++++++++++++++++++++++++++++++--
 3 files changed, 70 insertions(+), 5 deletions(-)

diff --git a/mats/5_3.ms b/mats/5_3.ms
index 09325a2963..ace926700c 100644
--- a/mats/5_3.ms
+++ b/mats/5_3.ms
@@ -1771,6 +1771,14 @@
     (eqv? (/ 0 2/3) 0)
     (eqv? (/ 0 2.0+3.0i) 0)
     (eqv? (/ 0 +nan.0+nan.0i) 0)
+    (eqv? (/ 1e300+1e300i (* 4 1e300+1e300i)) 0.25+0.0i)
+    (eqv? (/ #e1e300+1e300i (* 4 1e300+1e300i)) 0.25+0.0i)
+    (eqv? (/ 1e300+1e300i (* 4 #e1e300+1e300i)) 0.25+0.0i)
+    (eqv? (/ 1e-300+1e-300i (* 4 1e-300+1e-300i)) 0.25+0.0i)
+    (eqv? (/ #e1e-300+1e-300i (* 4 1e-300+1e-300i)) 0.25+0.0i)
+    (eqv? (/ 1e-300+1e-300i (* 4 #e1e-300+1e-300i)) 0.25+0.0i)
+    (eqv? (/ 0.0+0.0i 1+1e-320i) 0.0+0.0i)
+    (eqv? (/ 0.0+0.0i #e1+1e-320i) 0.0+0.0i)
     (test-cp0-expansion eqv? '(/ 1 2) 1/2)
     (test-cp0-expansion eqv? '(/ 1 -2) -1/2)
     (test-cp0-expansion eqv? '(/ 1/2 -2) -1/4)
diff --git a/s/5_3.ss b/s/5_3.ss
index 10a19b0ee5..13288beee9 100644
--- a/s/5_3.ss
+++ b/s/5_3.ss
@@ -2243,12 +2243,41 @@
                     (let ([t (/ x (+ (* c c) (* d d)))])
                       (make-rectangular (* c t) (- (* d t))))))]
              [($exactnum? $inexactnum?)
-              ;; a+bi / c+di => (ac+bd)/(cc+dd) + ((bc-ad)/(cc+dd))i
               (let ([a (real-part x)] [b (imag-part x)]
                     [c (real-part y)] [d (imag-part y)])
-                 (let ([t (+ (* c c) (* d d))])
+                ;; a+bi / c+di => (ac+bd)/(cc+dd) + ((bc-ad)/(cc+dd))i
+                (define (simpler-divide a b c d)
+                  ;; Direct calculuation does not work as well for complex numbers with
+                  ;; large parts, such as `(/ 1e+300+1e+300i 4e+300+4e+300i)`, but it
+                  ;; works better for small parts, as in `(/ 0.0+0.0i 1+1e-320i)`
+                  (let ([t (+ (* c c) (* d d))])
                     (make-rectangular (/ (+ (* a c) (* b d)) t)
-                                      (/ (- (* b c) (* a d)) t))))]
+                                      (/ (- (* b c) (* a d)) t))))
+                ;; Let r = c/d or d/c, depending on which is larger
+                (cond
+                 [(and ($exactnum? x) ($exactnum? y))
+                  (simpler-divide a b c d)]
+                 [(< (abs c) (abs d))
+                  (let ([r (/ d c)])
+                    (if (infinite? r)
+                        ;; Too large; try form that works better with small c or d
+                        (simpler-divide a b c d)
+                        ;; a+bi / c+di => 
+                        (let ([x (+ c (* d r))]) ; x = c+dd/c = (cc+dd)/c
+                          ;; (a+br)/x + ((b-ar)/x)i = (a+bd/c)c/(cc+dd) + ((b-ad/c)c/(cc+dd))i
+                          ;; = (ac+bd)/(cc+dd) + ((bc-ad)/(cc+dd))i
+                          (make-rectangular (/ (+ a (* b r)) x)
+                                            (/ (- b (* a r)) x)))))]
+                 [else
+                  (let ([r (/ c d)])
+                    (if (infinite? r)
+                        ;; Too large; try form that works better with small c or d
+                        (simpler-divide a b c d)
+                        (let ([x (+ d (* c r))]) ; x = d+cc/d = (cc+dd)/d
+                          ;; (b+ar)/x + ((br-a)/x)i = (b+ac/d)d/(cc+dd) + ((bc/d-a)d/(cc+dd))i
+                          ;; = (bd+ac)/(cc+dd) + ((bc-ad)/(cc+dd))i
+                          (make-rectangular (/ (+ b (* a r)) x)
+                                            (/ (- (* b r) a) x)))))]))]
              [else (nonnumber-error who x)])]
          [else (nonnumber-error who y)])))
 
diff --git a/s/library.ss b/s/library.ss
index f9dd164955..7f3356e21a 100644
--- a/s/library.ss
+++ b/s/library.ss
@@ -232,9 +232,37 @@
        ;; a+bi / c+di => (ac+bd)/(cc+dd) + ((bc-ad)/(cc+dd))i
        (let ([a ($inexactnum-real-part x)] [b ($inexactnum-imag-part x)]
              [c ($inexactnum-real-part y)] [d ($inexactnum-imag-part y)])
-          (let ([t (fl+ (fl* c c) (fl* d d))])
+         ;; a+bi / c+di => (ac+bd)/(cc+dd) + ((bc-ad)/(cc+dd))i
+         (define (simpler-divide a b c d)
+           ;; Direct calculuation does not work as well for complex numbers with
+           ;; large parts, such as `(/ 1e+300+1e+300i 4e+300+4e+300i)`, but it
+           ;; works better for small parts, as in `(/ 0.0+0.0i 1+1e-320i)`
+           (let ([t (fl+ (fl* c c) (fl* d d))])
              (fl-make-rectangular (fl/ (fl+ (fl* a c) (fl* b d)) t)
-                                  (fl/ (fl- (fl* b c) (fl* a d)) t))))]))
+                                  (fl/ (fl- (fl* b c) (fl* a d)) t))))
+         ;; Let r = c/d or d/c, depending on which is larger
+         (cond
+          [(fl< (flabs c) (flabs d))
+           (let ([r (fl/ d c)])
+             (if (infinity? r)
+                 ;; Too large; try form that works better with small c or d
+                 (simpler-divide a b c d)
+                 ;; a+bi / c+di => 
+                 (let ([x (fl+ c (fl* d r))]) ; x = c+dd/c = (cc+dd)/c
+                   ;; (a+br)/x + ((b-ar)/x)i = (a+bd/c)c/(cc+dd) + ((b-ad/c)c/(cc+dd))i
+                   ;; = (ac+bd)/(cc+dd) + ((bc-ad)/(cc+dd))i
+                   (fl-make-rectangular (fl/ (fl+ a (fl* b r)) x)
+                                        (fl/ (fl- b (fl* a r)) x)))))]
+          [else
+           (let ([r (fl/ c d)])
+             (if (infinity? r)
+                 ;; Too large; try form that works better with small c or d
+                 (simpler-divide a b c d)
+                 (let ([x (fl+ d (fl* c r))]) ; x = d+cc/d = (cc+dd)/d
+                   ;; (b+ar)/x + ((br-a)/x)i = (b+ac/d)d/(cc+dd) + ((bc/d-a)d/(cc+dd))i
+                   ;; = (bd+ac)/(cc+dd) + ((bc-ad)/(cc+dd))i
+                   (fl-make-rectangular (fl/ (fl+ b (fl* a r)) x)
+                                        (fl/ (fl- (fl* b r) a) x)))))]))]))
 
 (let ()
   (define char-oops