## [Sbcl-devel] Implementation of complex float operations

 [Sbcl-devel] Implementation of complex float operations From: Paul Khuong - 2009-06-20 07:32:19 Good evening. Before committing the SSE complex patch, I'd like your opinion on a few semi-related issues. 1. The current implementation for complex-complex division goes like: (let* ((rx (realpart x)) (ix (imagpart x)) (ry (realpart y)) (iy (imagpart y))) (if (> (abs ry) (abs iy)) (let* ((r (/ iy ry)) (dn (* ry (+ 1 (* r r))))) (complex (/ (+ rx (* ix r)) dn) (/ (- ix (* rx r)) dn))) (let* ((r (/ ry iy)) (dn (* iy (+ 1 (* r r))))) (complex (/ (+ (* rx r) ix) dn) (/ (- (* ix r) rx) dn))))) In his tome on semi-numerical algorithms, Knuth suggests to compute dn as (+ ry (* iy r)) or (+ ir (* ry r)). Indeed, squaring r only to cancel some of the denominator away is fishy. Unless someone can think of a reason for the current implementation, I don't think that change will be contentious. 2. In any case, the code above is far complete. Compared to what the C99 standard mandates, it pretty much fails to handle every edge case. On the other hand, in a performance-oriented application, the user probably expects the above, and not the maze of branches that GCC has to emit to comply. I don't think C99 is noted for the sanity of its complexes either... 3. The (CL) spec says that operations involving real and complex floats should be executed as though the reals had first been upgraded to complexes. That's not what currently happens, since some operations are optimised specially. e.g.: (* 1.0 #c(-1.0 #.sb-ext:single-float-positive-infinity)) => #C(-1.0 #.SB-EXT:SINGLE-FLOAT-POSITIVE-INFINITY) but (* #c(1.0 0.0) #c(-1.0 #.sb-ext:single-float-positive-infinity)) => #C(# #.SB-EXT:SINGLE-FLOAT-POSITIVE- INFINITY) ;; or an FP condition Following that part of the spec to the letter would likely imply emitting slower code that may return NaN-ful or less precise results than the current version, or adding more branches in the complex- complex case to check for zero imaginary parts (and likely other edge cases). The former isn't particularly desirable, and the latter seems lossy to me, since 0.0 does not necessarily mean 0. Should we instead decide to ignore the spec on that point and even code additional real/ complex transforms? Is it time for a fast-math--like declaration? Other thoughts? Paul Khuong