suzanne.soy/racket
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

improve precision of complex-number divide

Browse Source 
original commit: 4c9a7f6abb1258158d48fcdb656de300902cf3c7
This commit is contained in:
Matthew Flatt 2019-06-28 09:01:06 -06:00
parent 57c997042e
commit 35a0dfcafe
3 changed files with 70 additions and 5 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
mats/5_3.ms
View File

@ -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)

35
s/5_3.ss
View File

@ -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)])))

32
s/library.ss
View File

@ -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