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

The imaginary parts of inexact reals are ignored when doing complex

Browse Source 
addition or subtraction.

original commit: 748e9e47ad65c31d653907623b17689b47c76269
This commit is contained in:
Vincent St-Amour 2010-07-13 17:06:07 -04:00
parent 6ec1b3686b
commit 5f878b83bc
2 changed files with 81 additions and 43 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
collects/typed-scheme/optimizer/float.rkt
View File

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

122
collects/typed-scheme/optimizer/inexact-complex.rkt
View File

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