The imaginary parts of inexact reals are ignored when doing complex multiplication.

This commit is contained in:
Vincent St-Amour 2010-07-13 19:01:46 -04:00
parent b0788372ea
commit c645aa7ebc
4 changed files with 54 additions and 40 deletions

View File

@ -5,3 +5,4 @@
(* 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 3.0)

View File

@ -4,24 +4,19 @@
(unboxed-gensym-2 2.0+4.0i)
(unboxed-gensym-3 (unsafe-flreal-part unboxed-gensym-2))
(unboxed-gensym-4 (unsafe-flimag-part unboxed-gensym-2))
(unboxed-gensym-5 (unsafe-fl- (unsafe-fl* unboxed-gensym-1
unboxed-gensym-3)
(unsafe-fl* 0.0 unboxed-gensym-4)))
(unboxed-gensym-6 (unsafe-fl+ (unsafe-fl* 0.0 unboxed-gensym-3)
(unsafe-fl* unboxed-gensym-1
unboxed-gensym-4))))
(unboxed-gensym-5 (unsafe-fl* unboxed-gensym-1
unboxed-gensym-3))
(unboxed-gensym-6 (unsafe-fl* unboxed-gensym-1
unboxed-gensym-4)))
(unsafe-make-flrectangular unboxed-gensym-5 unboxed-gensym-6))
(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 (unsafe-fl- (unsafe-fl* unboxed-gensym-2
unboxed-gensym-4)
(unsafe-fl* unboxed-gensym-3 0.0)))
(unboxed-gensym-6 (unsafe-fl+ (unsafe-fl* unboxed-gensym-3
unboxed-gensym-4)
(unsafe-fl* unboxed-gensym-2
0.0))))
(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)
(unboxed-gensym-2 2.0+4.0i)
@ -30,12 +25,10 @@
(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-1
unboxed-gensym-3)
(unsafe-fl* 0.0 unboxed-gensym-4)))
(unboxed-gensym-11 (unsafe-fl+ (unsafe-fl* 0.0 unboxed-gensym-3)
(unsafe-fl* unboxed-gensym-1
unboxed-gensym-4)))
(unboxed-gensym-10 (unsafe-fl* unboxed-gensym-1
unboxed-gensym-3))
(unboxed-gensym-11 (unsafe-fl* unboxed-gensym-1
unboxed-gensym-4))
(unboxed-gensym-8 (unsafe-fl- (unsafe-fl* unboxed-gensym-10
unboxed-gensym-6)
(unsafe-fl* unboxed-gensym-11
@ -52,12 +45,10 @@
(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-2
unboxed-gensym-4)
(unsafe-fl* unboxed-gensym-3 0.0)))
(unboxed-gensym-11 (unsafe-fl+ (unsafe-fl* unboxed-gensym-3
unboxed-gensym-4)
(unsafe-fl* unboxed-gensym-2 0.0)))
(unboxed-gensym-10 (unsafe-fl* unboxed-gensym-2
unboxed-gensym-4))
(unboxed-gensym-11 (unsafe-fl* unboxed-gensym-3
unboxed-gensym-4))
(unboxed-gensym-8 (unsafe-fl- (unsafe-fl* unboxed-gensym-10
unboxed-gensym-6)
(unsafe-fl* unboxed-gensym-11
@ -82,10 +73,22 @@
unboxed-gensym-5)
(unsafe-fl* unboxed-gensym-2
unboxed-gensym-6)))
(unboxed-gensym-8 (unsafe-fl- (unsafe-fl* unboxed-gensym-10
unboxed-gensym-7)
(unsafe-fl* unboxed-gensym-11 0.0)))
(unboxed-gensym-9 (unsafe-fl+ (unsafe-fl* unboxed-gensym-11
unboxed-gensym-7)
(unsafe-fl* unboxed-gensym-10 0.0))))
(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))

View File

@ -1,7 +1,8 @@
#lang typed/scheme
#lang typed/scheme
(require racket/unsafe/ops)
(* 1.0 2.0+4.0i)
(* 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 3.0)

View File

@ -73,7 +73,6 @@
#`(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
;; we currently don't skip imaginary parts of reals
#,@(let ((l (map (lambda (x) (if (syntax->datum x) x #'0.0))
(syntax->list #'(c1.imag-part c2.imag-part cs.imag-part ...)))))
(let loop ([o1 #'c1.real-part]
@ -91,13 +90,23 @@
(reverse res)
(loop (car rs) (car is) (cdr e1) (cdr e2) (cdr rs) (cdr is)
;; complex multiplication, imag part, then real part (reverse)
(list* #`(#,(car is)
(unsafe-fl+ (unsafe-fl* #,o2 #,(car e1))
(unsafe-fl* #,o1 #,(car e2))))
#`(#,(car rs)
(unsafe-fl- (unsafe-fl* #,o1 #,(car e1))
(unsafe-fl* #,o2 #,(car e2))))
res))))))))
;; we eliminate operations on the imaginary parts of reals
(let ((o-real? (equal? (syntax->datum o2) 0.0))
(e-real? (equal? (syntax->datum (car e2)) 0.0)))
(list* #`(#,(car is)
#,(cond ((and o-real? e-real?) #'0.0)
(o-real? #`(unsafe-fl* #,o1 #,(car e2)))
(e-real? #`(unsafe-fl* #,o2 #,(car e1)))
(else
#`(unsafe-fl+ (unsafe-fl* #,o2 #,(car e1))
(unsafe-fl* #,o1 #,(car e2))))))
#`(#,(car rs)
#,(cond ((or o-real? e-real?)
#`(unsafe-fl* #,o1 #,(car e1)))
(else
#`(unsafe-fl- (unsafe-fl* #,o1 #,(car e1))
(unsafe-fl* #,o2 #,(car e2))))))
res)))))))))
(pattern (#%plain-app (~and op (~literal /))
c1:unboxed-inexact-complex-opt-expr
c2:unboxed-inexact-complex-opt-expr