The imaginary parts of reals are ignored when doing complex division.

This commit is contained in:
Vincent St-Amour 2010-07-14 12:59:37 -04:00
parent 443d8b9f91
commit c587038b33
5 changed files with 130 additions and 99 deletions

View File

@ -8,4 +8,6 @@
(/ 1.0+2.0i 2.0)
(/ 1.0 2.0+4.0i 3.0+6.0i)
(/ 1.0+2.0i 2.0 3.0+6.0i)
(/ 1.0+2.0i 2.0+4.0i 3.0)))
(/ 1.0+2.0i 2.0+4.0i 3.0)
(/ 1.0+2.0i 2.0 3.0)
(/ 1.0 2.0 3.0+6.0i)))

View File

@ -9,28 +9,28 @@
(unboxed-gensym-7 3.0+6.0i)
(unboxed-gensym-8 (unsafe-flreal-part unboxed-gensym-7))
(unboxed-gensym-9 (unsafe-flimag-part unboxed-gensym-7))
(unboxed-gensym-12 (unsafe-fl+ (unsafe-fl* unboxed-gensym-5 unboxed-gensym-5)
(unboxed-gensym-14 (unsafe-fl+ (unsafe-fl* unboxed-gensym-5 unboxed-gensym-5)
(unsafe-fl* unboxed-gensym-6 unboxed-gensym-6)))
(unboxed-gensym-13 (unsafe-fl+ (unsafe-fl* unboxed-gensym-8 unboxed-gensym-8)
(unsafe-fl* unboxed-gensym-9 unboxed-gensym-9)))
(unboxed-gensym-14 (unsafe-fl/ (unsafe-fl+ (unsafe-fl* unboxed-gensym-2
(unboxed-gensym-12 (unsafe-fl/ (unsafe-fl+ (unsafe-fl* unboxed-gensym-2
unboxed-gensym-5)
(unsafe-fl* unboxed-gensym-3
unboxed-gensym-6))
unboxed-gensym-12))
(unboxed-gensym-15 (unsafe-fl/ (unsafe-fl- (unsafe-fl* unboxed-gensym-3
unboxed-gensym-14))
(unboxed-gensym-13 (unsafe-fl/ (unsafe-fl- (unsafe-fl* unboxed-gensym-3
unboxed-gensym-5)
(unsafe-fl* unboxed-gensym-2
unboxed-gensym-6))
unboxed-gensym-12))
(unboxed-gensym-10 (unsafe-fl/ (unsafe-fl+ (unsafe-fl* unboxed-gensym-14
unboxed-gensym-14))
(unboxed-gensym-15 (unsafe-fl+ (unsafe-fl* unboxed-gensym-8 unboxed-gensym-8)
(unsafe-fl* unboxed-gensym-9 unboxed-gensym-9)))
(unboxed-gensym-10 (unsafe-fl/ (unsafe-fl+ (unsafe-fl* unboxed-gensym-12
unboxed-gensym-8)
(unsafe-fl* unboxed-gensym-15
(unsafe-fl* unboxed-gensym-13
unboxed-gensym-9))
unboxed-gensym-13))
(unboxed-gensym-11 (unsafe-fl/ (unsafe-fl- (unsafe-fl* unboxed-gensym-15
unboxed-gensym-15))
(unboxed-gensym-11 (unsafe-fl/ (unsafe-fl- (unsafe-fl* unboxed-gensym-13
unboxed-gensym-8)
(unsafe-fl* unboxed-gensym-14
(unsafe-fl* unboxed-gensym-12
unboxed-gensym-9))
unboxed-gensym-13)))
unboxed-gensym-15)))
(unsafe-make-flrectangular unboxed-gensym-10 unboxed-gensym-11)))

View File

@ -12,11 +12,10 @@
unboxed-gensym-3)
(unsafe-fl* unboxed-gensym-4
unboxed-gensym-4)))
(unboxed-gensym-5 (unsafe-fl/ (unsafe-fl+ (unsafe-fl* unboxed-gensym-1
unboxed-gensym-3)
(unsafe-fl* 0.0 unboxed-gensym-4))
(unboxed-gensym-5 (unsafe-fl/ (unsafe-fl* unboxed-gensym-1
unboxed-gensym-3)
unboxed-gensym-7))
(unboxed-gensym-6 (unsafe-fl/ (unsafe-fl- (unsafe-fl* 0.0 unboxed-gensym-3)
(unboxed-gensym-6 (unsafe-fl/ (unsafe-fl- 0.0
(unsafe-fl* unboxed-gensym-1
unboxed-gensym-4))
unboxed-gensym-7)))
@ -26,17 +25,8 @@
(unboxed-gensym-2 (unsafe-flreal-part unboxed-gensym-1))
(unboxed-gensym-3 (unsafe-flimag-part unboxed-gensym-1))
(unboxed-gensym-4 2.0)
(unboxed-gensym-7 (unsafe-fl+ (unsafe-fl* unboxed-gensym-4
unboxed-gensym-4)
(unsafe-fl* 0.0 0.0)))
(unboxed-gensym-5 (unsafe-fl/ (unsafe-fl+ (unsafe-fl* unboxed-gensym-2
unboxed-gensym-4)
(unsafe-fl* unboxed-gensym-3 0.0))
unboxed-gensym-7))
(unboxed-gensym-6 (unsafe-fl/ (unsafe-fl- (unsafe-fl* unboxed-gensym-3
unboxed-gensym-4)
(unsafe-fl* unboxed-gensym-2 0.0))
unboxed-gensym-7)))
(unboxed-gensym-5 (unsafe-fl/ unboxed-gensym-2 unboxed-gensym-4))
(unboxed-gensym-6 (unsafe-fl/ unboxed-gensym-3 unboxed-gensym-4)))
(unsafe-make-flrectangular unboxed-gensym-5 unboxed-gensym-6))
(let* ((unboxed-gensym-1 1.0)
@ -46,32 +36,31 @@
(unboxed-gensym-5 3.0+6.0i)
(unboxed-gensym-6 (unsafe-flreal-part unboxed-gensym-5))
(unboxed-gensym-7 (unsafe-flimag-part unboxed-gensym-5))
(unboxed-gensym-10 (unsafe-fl+ (unsafe-fl* unboxed-gensym-3
(unboxed-gensym-12 (unsafe-fl+ (unsafe-fl* unboxed-gensym-3
unboxed-gensym-3)
(unsafe-fl* unboxed-gensym-4
unboxed-gensym-4)))
(unboxed-gensym-11 (unsafe-fl+ (unsafe-fl* unboxed-gensym-6
(unboxed-gensym-10 (unsafe-fl/ (unsafe-fl* unboxed-gensym-1
unboxed-gensym-3)
unboxed-gensym-12))
(unboxed-gensym-11 (unsafe-fl/ (unsafe-fl- 0.0
(unsafe-fl* unboxed-gensym-1
unboxed-gensym-4))
unboxed-gensym-12))
(unboxed-gensym-13 (unsafe-fl+ (unsafe-fl* unboxed-gensym-6
unboxed-gensym-6)
(unsafe-fl* unboxed-gensym-7
unboxed-gensym-7)))
(unboxed-gensym-12 (unsafe-fl/ (unsafe-fl+ (unsafe-fl* unboxed-gensym-1
unboxed-gensym-3)
(unsafe-fl* 0.0 unboxed-gensym-4))
unboxed-gensym-10))
(unboxed-gensym-13 (unsafe-fl/ (unsafe-fl- (unsafe-fl* 0.0 unboxed-gensym-3)
(unsafe-fl* unboxed-gensym-1
unboxed-gensym-4))
unboxed-gensym-10))
(unboxed-gensym-8 (unsafe-fl/ (unsafe-fl+ (unsafe-fl* unboxed-gensym-12
(unboxed-gensym-8 (unsafe-fl/ (unsafe-fl+ (unsafe-fl* unboxed-gensym-10
unboxed-gensym-6)
(unsafe-fl* unboxed-gensym-13
(unsafe-fl* unboxed-gensym-11
unboxed-gensym-7))
unboxed-gensym-11))
(unboxed-gensym-9 (unsafe-fl/ (unsafe-fl- (unsafe-fl* unboxed-gensym-13
unboxed-gensym-13))
(unboxed-gensym-9 (unsafe-fl/ (unsafe-fl- (unsafe-fl* unboxed-gensym-11
unboxed-gensym-6)
(unsafe-fl* unboxed-gensym-12
(unsafe-fl* unboxed-gensym-10
unboxed-gensym-7))
unboxed-gensym-11)))
unboxed-gensym-13)))
(unsafe-make-flrectangular unboxed-gensym-8 unboxed-gensym-9))
(let* ((unboxed-gensym-1 1.0+2.0i)
@ -81,31 +70,22 @@
(unboxed-gensym-5 3.0+6.0i)
(unboxed-gensym-6 (unsafe-flreal-part unboxed-gensym-5))
(unboxed-gensym-7 (unsafe-flimag-part unboxed-gensym-5))
(unboxed-gensym-10 (unsafe-fl+ (unsafe-fl* unboxed-gensym-4
unboxed-gensym-4)
(unsafe-fl* 0.0 0.0)))
(unboxed-gensym-11 (unsafe-fl+ (unsafe-fl* unboxed-gensym-6
(unboxed-gensym-10 (unsafe-fl/ unboxed-gensym-2 unboxed-gensym-4))
(unboxed-gensym-11 (unsafe-fl/ unboxed-gensym-3 unboxed-gensym-4))
(unboxed-gensym-13 (unsafe-fl+ (unsafe-fl* unboxed-gensym-6
unboxed-gensym-6)
(unsafe-fl* unboxed-gensym-7
unboxed-gensym-7)))
(unboxed-gensym-12 (unsafe-fl/ (unsafe-fl+ (unsafe-fl* unboxed-gensym-2
unboxed-gensym-4)
(unsafe-fl* unboxed-gensym-3 0.0))
unboxed-gensym-10))
(unboxed-gensym-13 (unsafe-fl/ (unsafe-fl- (unsafe-fl* unboxed-gensym-3
unboxed-gensym-4)
(unsafe-fl* unboxed-gensym-2 0.0))
unboxed-gensym-10))
(unboxed-gensym-8 (unsafe-fl/ (unsafe-fl+ (unsafe-fl* unboxed-gensym-12
(unboxed-gensym-8 (unsafe-fl/ (unsafe-fl+ (unsafe-fl* unboxed-gensym-10
unboxed-gensym-6)
(unsafe-fl* unboxed-gensym-13
(unsafe-fl* unboxed-gensym-11
unboxed-gensym-7))
unboxed-gensym-11))
(unboxed-gensym-9 (unsafe-fl/ (unsafe-fl- (unsafe-fl* unboxed-gensym-13
unboxed-gensym-13))
(unboxed-gensym-9 (unsafe-fl/ (unsafe-fl- (unsafe-fl* unboxed-gensym-11
unboxed-gensym-6)
(unsafe-fl* unboxed-gensym-12
(unsafe-fl* unboxed-gensym-10
unboxed-gensym-7))
unboxed-gensym-11)))
unboxed-gensym-13)))
(unsafe-make-flrectangular unboxed-gensym-8 unboxed-gensym-9))
(let* ((unboxed-gensym-1 1.0+2.0i)
@ -115,29 +95,54 @@
(unboxed-gensym-5 (unsafe-flreal-part unboxed-gensym-4))
(unboxed-gensym-6 (unsafe-flimag-part unboxed-gensym-4))
(unboxed-gensym-7 3.0)
(unboxed-gensym-10 (unsafe-fl+ (unsafe-fl* unboxed-gensym-5
(unboxed-gensym-12 (unsafe-fl+ (unsafe-fl* unboxed-gensym-5
unboxed-gensym-5)
(unsafe-fl* unboxed-gensym-6
unboxed-gensym-6)))
(unboxed-gensym-11 (unsafe-fl+ (unsafe-fl* unboxed-gensym-7
unboxed-gensym-7)
(unsafe-fl* 0.0 0.0)))
(unboxed-gensym-12 (unsafe-fl/ (unsafe-fl+ (unsafe-fl* unboxed-gensym-2
(unboxed-gensym-10 (unsafe-fl/ (unsafe-fl+ (unsafe-fl* unboxed-gensym-2
unboxed-gensym-5)
(unsafe-fl* unboxed-gensym-3
unboxed-gensym-6))
unboxed-gensym-10))
(unboxed-gensym-13 (unsafe-fl/ (unsafe-fl- (unsafe-fl* unboxed-gensym-3
unboxed-gensym-12))
(unboxed-gensym-11 (unsafe-fl/ (unsafe-fl- (unsafe-fl* unboxed-gensym-3
unboxed-gensym-5)
(unsafe-fl* unboxed-gensym-2
unboxed-gensym-6))
unboxed-gensym-10))
(unboxed-gensym-8 (unsafe-fl/ (unsafe-fl+ (unsafe-fl* unboxed-gensym-12
unboxed-gensym-7)
(unsafe-fl* unboxed-gensym-13 0.0))
unboxed-gensym-12))
(unboxed-gensym-8 (unsafe-fl/ unboxed-gensym-10
unboxed-gensym-7))
(unboxed-gensym-9 (unsafe-fl/ unboxed-gensym-11
unboxed-gensym-7)))
(unsafe-make-flrectangular unboxed-gensym-8 unboxed-gensym-9))
(let* ((unboxed-gensym-1 1.0+2.0i)
(unboxed-gensym-2 (unsafe-flreal-part unboxed-gensym-1))
(unboxed-gensym-3 (unsafe-flimag-part unboxed-gensym-1))
(unboxed-gensym-4 2.0)
(unboxed-gensym-5 3.0)
(unboxed-gensym-8 (unsafe-fl/ unboxed-gensym-2 unboxed-gensym-4))
(unboxed-gensym-9 (unsafe-fl/ unboxed-gensym-3 unboxed-gensym-4))
(unboxed-gensym-6 (unsafe-fl/ unboxed-gensym-8 unboxed-gensym-5))
(unboxed-gensym-7 (unsafe-fl/ unboxed-gensym-9 unboxed-gensym-5)))
(unsafe-make-flrectangular unboxed-gensym-6 unboxed-gensym-7))
(let* ((unboxed-gensym-1 1.0)
(unboxed-gensym-2 2.0)
(unboxed-gensym-3 3.0+6.0i)
(unboxed-gensym-4 (unsafe-flreal-part unboxed-gensym-3))
(unboxed-gensym-5 (unsafe-flimag-part unboxed-gensym-3))
(unboxed-gensym-8 (unsafe-fl/ unboxed-gensym-1 unboxed-gensym-2))
(unboxed-gensym-9 0.0)
(unboxed-gensym-11 (unsafe-fl+ (unsafe-fl* unboxed-gensym-4
unboxed-gensym-4)
(unsafe-fl* unboxed-gensym-5
unboxed-gensym-5)))
(unboxed-gensym-6 (unsafe-fl/ (unsafe-fl+ (unsafe-fl* unboxed-gensym-8
unboxed-gensym-4)
(unsafe-fl* unboxed-gensym-9
unboxed-gensym-5))
unboxed-gensym-11))
(unboxed-gensym-9 (unsafe-fl/ (unsafe-fl- (unsafe-fl* unboxed-gensym-13
unboxed-gensym-7)
(unsafe-fl* unboxed-gensym-12 0.0))
(unboxed-gensym-7 (unsafe-fl/ (unsafe-fl- (unsafe-fl* unboxed-gensym-9
unboxed-gensym-4)
(unsafe-fl* unboxed-gensym-8
unboxed-gensym-5))
unboxed-gensym-11)))
(unsafe-make-flrectangular unboxed-gensym-8 unboxed-gensym-9))))
(unsafe-make-flrectangular unboxed-gensym-6 unboxed-gensym-7))))

View File

@ -8,4 +8,6 @@
(/ 1.0+2.0i 2.0)
(/ 1.0 2.0+4.0i 3.0+6.0i)
(/ 1.0+2.0i 2.0 3.0+6.0i)
(/ 1.0+2.0i 2.0+4.0i 3.0)))
(/ 1.0+2.0i 2.0+4.0i 3.0)
(/ 1.0+2.0i 2.0 3.0)
(/ 1.0 2.0 3.0+6.0i)))

View File

@ -114,45 +114,67 @@
#:with real-part (unboxed-gensym)
#:with imag-part (unboxed-gensym)
#:with reals (syntax->list #'(c1.real-part c2.real-part cs.real-part ...))
;; we currently don't skip imaginary parts of reals
#:with imags (map (lambda (x) (if (syntax->datum x) x #'0.0))
(syntax->list #'(c1.imag-part c2.imag-part cs.imag-part ...)))
#:with (denominators ...)
(for/list
([e1 (cdr (syntax->list #'reals))]
[e2 (cdr (syntax->list #'imags))])
#`(#,(unboxed-gensym) (unsafe-fl+ (unsafe-fl* #,e1 #,e1) (unsafe-fl* #,e2 #,e2))))
#:with (bindings ...)
(begin (log-optimization "unboxed binary inexact complex" #'op)
#`(c1.bindings ... c2.bindings ... cs.bindings ... ... denominators ...
#`(c1.bindings ... c2.bindings ... cs.bindings ... ...
;; we want to bind the intermediate results to reuse them
;; the final results are bound to real-part and imag-part
#,@(let loop ([o1 (car (syntax->list #'reals))]
[o2 (car (syntax->list #'imags))]
[e1 (cdr (syntax->list #'reals))]
[e2 (cdr (syntax->list #'imags))]
[d (map (lambda (x) (car (syntax-e x)))
(syntax->list #'(denominators ...)))]
[rs (append (map (lambda (x) (unboxed-gensym))
(syntax->list #'(cs.real-part ...)))
(list #'real-part))]
[is (append (map (lambda (x) (unboxed-gensym))
(syntax->list #'(cs.imag-part ...)))
(list #'imag-part))]
[ds (map (lambda (x) (unboxed-gensym))
(syntax->list #'(c2.real-part cs.real-part ...)))]
[res '()])
(if (null? e1)
(reverse res)
(loop (car rs) (car is) (cdr e1) (cdr e2) (cdr d) (cdr rs) (cdr is)
;; complex division, imag part, then real part (reverse)
(list* #`(#,(car is)
(unsafe-fl/ (unsafe-fl- (unsafe-fl* #,o2 #,(car e1))
(unsafe-fl* #,o1 #,(car e2)))
#,(car d)))
#`(#,(car rs)
(unsafe-fl/ (unsafe-fl+ (unsafe-fl* #,o1 #,(car e1))
(unsafe-fl* #,o2 #,(car e2)))
#,(car d)))
res)))))))
(loop (car rs) (car is) (cdr e1) (cdr e2) (cdr rs) (cdr is) (cdr ds)
;; complex division, imag part, real part, then denominator (reverse)
(let ((o-real? (equal? (syntax->datum o2) 0.0))
(e-real? (equal? (syntax->datum (car e2)) 0.0)))
(cond [(and o-real? e-real?)
(list*
#`(#,(car is) 0.0) ; currently not propagated
#`(#,(car rs) (unsafe-fl/ #,o1 #,(car e1)))
res)]
[o-real?
(list*
#`(#,(car is)
(unsafe-fl/ (unsafe-fl- 0.0
(unsafe-fl* #,o1 #,(car e2)))
#,(car ds)))
#`(#,(car rs) (unsafe-fl/ (unsafe-fl* #,o1 #,(car e1))
#,(car ds)))
#`(#,(car ds) (unsafe-fl+ (unsafe-fl* #,(car e1) #,(car e1))
(unsafe-fl* #,(car e2) #,(car e2))))
res)]
[e-real?
(list*
#`(#,(car is) (unsafe-fl/ #,o2 #,(car e1)))
#`(#,(car rs) (unsafe-fl/ #,o1 #,(car e1)))
res)]
[else
(list*
#`(#,(car is)
(unsafe-fl/ (unsafe-fl- (unsafe-fl* #,o2 #,(car e1))
(unsafe-fl* #,o1 #,(car e2)))
#,(car ds)))
#`(#,(car rs)
(unsafe-fl/ (unsafe-fl+ (unsafe-fl* #,o1 #,(car e1))
(unsafe-fl* #,o2 #,(car e2)))
#,(car ds)))
#`(#,(car ds)
(unsafe-fl+ (unsafe-fl* #,(car e1) #,(car e1))
(unsafe-fl* #,(car e2) #,(car e2))))
res)]))))))))
(pattern (#%plain-app (~and op (~literal conjugate)) c:unboxed-inexact-complex-opt-expr)
#:with real-part #'c.real-part
#:with imag-part (unboxed-gensym)