From 5f878b83bc8c19df5cb8e3c16d57414169125ece Mon Sep 17 00:00:00 2001 From: Vincent St-Amour Date: Tue, 13 Jul 2010 17:06:07 -0400 Subject: [PATCH] The imaginary parts of inexact reals are ignored when doing complex addition or subtraction. original commit: 748e9e47ad65c31d653907623b17689b47c76269 --- collects/typed-scheme/optimizer/float.rkt | 2 +- .../optimizer/inexact-complex.rkt | 122 ++++++++++++------ 2 files changed, 81 insertions(+), 43 deletions(-) diff --git a/collects/typed-scheme/optimizer/float.rkt b/collects/typed-scheme/optimizer/float.rkt index 15675b37..2c9f02cb 100644 --- a/collects/typed-scheme/optimizer/float.rkt +++ b/collects/typed-scheme/optimizer/float.rkt @@ -7,7 +7,7 @@ (types abbrev type-table utils subtype) (optimizer utils fixnum)) -(provide float-opt-expr float-op mk-float-tbl) +(provide float-opt-expr) (define (mk-float-tbl generic) diff --git a/collects/typed-scheme/optimizer/inexact-complex.rkt b/collects/typed-scheme/optimizer/inexact-complex.rkt index 3d65c7d8..53f6848b 100644 --- a/collects/typed-scheme/optimizer/inexact-complex.rkt +++ b/collects/typed-scheme/optimizer/inexact-complex.rkt @@ -4,7 +4,7 @@ "../utils/utils.rkt" (for-template scheme/base scheme/math scheme/flonum scheme/unsafe/ops) (types abbrev type-table utils subtype) - (optimizer utils float)) + (optimizer utils)) (provide inexact-complex-opt-expr) @@ -14,7 +14,7 @@ ;; we keep the real and imaginary parts unboxed as long as we stay within ;; complex operations (define-syntax-class unboxed-inexact-complex-opt-expr - (pattern (#%plain-app (~and (~var op (float-op binary-inexact-complex-ops)) (~or (~literal +) (~literal -))) + (pattern (#%plain-app (~and op (~literal +)) c1:unboxed-inexact-complex-opt-expr c2:unboxed-inexact-complex-opt-expr cs:unboxed-inexact-complex-opt-expr ...) @@ -25,10 +25,43 @@ #`(#,@(append (syntax->list #'(c1.bindings ... c2.bindings ... cs.bindings ... ...)) (list #`(real-part #,(for/fold ((o #'c1.real-part)) ((e (syntax->list #'(c2.real-part cs.real-part ...)))) - #`(op.unsafe #,o #,e))) - #`(imag-part #,(for/fold ((o #'c1.imag-part)) - ((e (syntax->list #'(c2.imag-part cs.imag-part ...)))) - #`(op.unsafe #,o #,e)))))))) + #`(unsafe-fl+ #,o #,e))) + ;; we can skip the imaginary parts of reals (#f) + #`(imag-part + #,(let ((l (filter syntax->datum + (syntax->list #'(c1.imag-part c2.imag-part cs.imag-part ...))))) + (case (length l) + ((0) #'0.0) + ((1) (car l)) + (else + (for/fold ((o (car l))) + ((e (cdr l))) + #`(unsafe-fl+ #,o #,e))))))))))) + (pattern (#%plain-app (~and op (~literal -)) + c1:unboxed-inexact-complex-opt-expr + c2:unboxed-inexact-complex-opt-expr + cs:unboxed-inexact-complex-opt-expr ...) + #:with real-part (unboxed-gensym) + #:with imag-part (unboxed-gensym) + #:with (bindings ...) + (begin (log-optimization "unboxed binary inexact complex" #'op) + #`(#,@(append (syntax->list #'(c1.bindings ... c2.bindings ... cs.bindings ... ...)) + (list #`(real-part #,(for/fold ((o #'c1.real-part)) + ((e (syntax->list #'(c2.real-part cs.real-part ...)))) + #`(unsafe-fl- #,o #,e))) + ;; unlike addition, we simply can't skip imaginary parts of reals + #`(imag-part + #,(let* ((l1 (map (lambda (x) (if (syntax->datum x) x #'0.0)) + (syntax->list #'(c1.imag-part c2.imag-part cs.imag-part ...)))) + ;; but we can skip all but the first 0 + (l2 (filter (lambda (x) (not (equal? (syntax->datum x) 0.0))) + (cdr l1)))) + (case (length l2) + ((0) (car l1)) + (else + (for/fold ((o (car l1))) + ((e l2)) + #`(unsafe-fl- #,o #,e))))))))))) (pattern (#%plain-app (~and op (~literal *)) c1:unboxed-inexact-complex-opt-expr c2:unboxed-inexact-complex-opt-expr @@ -40,48 +73,55 @@ #`(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 #'c1.real-part] - [o2 #'c1.imag-part] - [e1 (syntax->list #'(c2.real-part cs.real-part ...))] - [e2 (syntax->list #'(c2.imag-part cs.imag-part ...))] - [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))] - [res '()]) - (if (null? e1) - (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 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] + [o2 (car l)] + [e1 (syntax->list #'(c2.real-part cs.real-part ...))] + [e2 (cdr l)] + [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))] + [res '()]) + (if (null? e1) + (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)))))))) (pattern (#%plain-app (~and op (~literal /)) c1:unboxed-inexact-complex-opt-expr c2:unboxed-inexact-complex-opt-expr cs:unboxed-inexact-complex-opt-expr ...) #: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 (syntax->list #'(c2.real-part cs.real-part ...))] - [e2 (syntax->list #'(c2.imag-part cs.imag-part ...))]) + ([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 ... ;; we want to bind the intermediate results to reuse them ;; the final results are bound to real-part and imag-part - #,@(let loop ([o1 #'c1.real-part] - [o2 #'c1.imag-part] - [e1 (syntax->list #'(c2.real-part cs.real-part ...))] - [e2 (syntax->list #'(c2.imag-part cs.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)) @@ -115,10 +155,9 @@ ;; special handling of inexact reals #:when (subtypeof? #'e -Flonum) #:with real-part (unboxed-gensym) - #:with imag-part (unboxed-gensym) + #:with imag-part #f #:with (bindings ...) - #`((real-part #,((optimize) #'e)) - (imag-part 0.0))) + #`((real-part #,((optimize) #'e)))) (pattern e:expr ;; can't work on inexact reals, which are a subtype of inexact ;; complexes, so this has to be equality @@ -134,8 +173,9 @@ (define-syntax-class inexact-complex-unary-op (pattern (~or (~literal real-part) (~literal flreal-part)) #:with unsafe #'unsafe-flreal-part) (pattern (~or (~literal imag-part) (~literal flimag-part)) #:with unsafe #'unsafe-flimag-part)) -(define binary-inexact-complex-ops - (mk-float-tbl (list #'+ #'- #'* #'/))) + +(define-syntax-class inexact-complex-binary-op + (pattern (~or (~literal +) (~literal -) (~literal *) (~literal /) (~literal conjugate)))) (define-syntax-class inexact-complex-expr (pattern e:expr @@ -147,9 +187,7 @@ #:with opt (begin (log-optimization "unary inexact complex" #'op) #'(op.unsafe n.opt))) - (pattern (~and exp (#%plain-app (~or (~var op (float-op binary-inexact-complex-ops)) - (~and op (~literal conjugate))) - e:inexact-complex-expr ...)) + (pattern (~and exp (#%plain-app op:inexact-complex-binary-op e:inexact-complex-expr ...)) #:when (isoftype? #'exp -InexactComplex) #:with exp*:unboxed-inexact-complex-opt-expr #'exp #:with opt