More on compiler optimization
"Paul Werkowski (as pw at snoopy dot qozzy dot com)" <[email protected]>
| Newsgroups | gmane.lisp.lispworks.general |
|---|---|
| Message-ID | <[email protected]> |
For the the past month or so I have been searching for the cause of lisp garbage generation where compiling with (declare (:explain :boxing)) was silent. My test case was/is my Lispworks implementation of the Fast Fourier Transform (FFT) algorithm commonly used in Digital Signal Processing. I have now succeeded in creating a zero allocation (garbage) implementation. While my implementation generates versions for both single-float and double-float data types, only the double-float version on 64-bit LW 8.0.1 fits that description. The main thing to get zero garbage (in addition to proper type and optimize declarations) is to ensure a boxed object, like a DOUBLE-FLOAT value, does not escape the compilation environment it exist in. Such an object can exist on the heap (default), as a raw value in a typed array, or on the control stack. They may also be created by compiler internal functions. It is the later case that can only be found by examining a disassembly listing. Finding one of a family of names starting COMPILER::RAW-FAST-BOX- in the listing is a sure sign that an allocation is happening. Such is necessary if the boxed item is to be an argument to an external (to the compilation environment) function which is expecting a pointer an object containing the raw value. What works is to place functions that need to efficiently pass double-float values around in local functions via FLET. Also it helps to avoid generic CL library functions such as cl:cis. (flet ((%%cis (x) (complex (cos x)(sin x)))) ....) Results? On my 2016 Dell XPS 8910 with quad 4 GHZ dual core processors, a 1024 point transform executes in 40 microseconds (down from 120) with 0 bytes allocated (down from 171,728). That is for the double-float version. A single-float version (created from same code but with double-float changed to single-float) takes 3.6 milliseconds with 687,576 bytes of garbage generated. LW 8.0.1 seems to not optimize AREFs into a typed array. My implementation is attached for the readers entertainment. Compiling the file will produce, in the :FFT package, CFFT-DOUBLE-FLOAT & CFFT-SINGLE-FLOAT. Enjoy! Paul
cfft-export.lisp
(text/plain, 6.9 KB)
;;; cfft-basic.lisp -*- Package: fft -*-
;;;
;;;Copyright (c) 1997, 2025 Paul F. Werkowski
;;;
;;;
(defpackage "FFT" (:use "COMMON-LISP")
(:export "CFFT-DOUBLE-FLOAT" "CFFT-SINGLE-FLOAT"))
(in-package :fft)
(defun next-power-of-two (N)
(ash 1 (ceiling (log N 2))))
(defun power-of-two-p (N)
(and (realp N)(plusp N)(= 1 (logcount N))))
(defun power-of-two-length-p (v)
(power-of-two-p (length v)))
(deftype fft-vector (&optional type size)
`(and (array ,type (,size))
(satisfies power-of-two-length-p)))
(defconstant +fft-max-power-of-2+ 17)
(deftype fft-index () '(integer 0 131072)) ; 2^17
;;;
;;; Bit Reversing from The Fast Fourier Transform by James S. Walker
(defun buneman-permutation-vector(P)
"Returns a vector of bit-reversed index values.
P for FFT of length (expt 2 P)"
(declare (type (integer 1) P))
(let* ((N (ash 1 P))
(J (make-array N :element-type 'fixnum)))
(setf (aref j 0) 0
(aref j 1) 1)
(do ((L 2 (+ L L))
(i 2 (1+ i)))
((> i P) j)
(dotimes (k L)
(incf (aref j k)(aref j k))
(setf (aref j (+ k L))(1+ (aref j k)))))))
;;; Memoize vectors used.
(defvar *cached-permutation-vectors*
#+(and cmu (not hash-new)) (make-hash-table :test #'eq :weak-p t)
#-(and cmu (not hash-new)) (make-hash-table :test #'eq)
"Cache of precomputed BR permutation vectors.")
(defun get-permutation-vector (power)
"Fetch a BR permutation vector for FFT length 2^P"
(declare (type (integer 1 20) power)
(values (simple-array fixnum (*))))
(let ((ht *cached-permutation-vectors*))
(multiple-value-bind (pv foundp) (gethash power ht)
(if foundp pv
(setf (gethash power ht)
(buneman-permutation-vector power))))))
(eval-when (:load-toplevel :execute)
(do ((p 1 (1+ p)))
((>= p +fft-max-power-of-2+))
(get-permutation-vector p)))
(defmacro def-bit-reverse (name ftype)
`(defun ,name (Z)
(declare (type (simple-array (complex ,ftype) (*)) Z)
(optimize (float 0)(safety 0)(hcl:fixnum-safety 0)))
(assert (= 1 (logcount (length Z)))) ; power of 2 length
(let* ((R (1- (integer-length (length Z))))
(m (ash 1 r))
(r2 (ash r -1))
(n1 (if (evenp r) m (ash m -1)))
(n9 (isqrt n1))
(n8 (if (evenp r) n9 (* n9 2)))
(j (get-permutation-vector r2)))
(declare (type (fixnum 0 #.array-dimension-limit) R m r2 n1 n9 n8)
(type (simple-array (fixnum 0)(*)) j))
;; Walker's optimal swaps
(do ((m7 0 n9)
(y 0 n9)
(c (if (evenp r) 1 0)
(1+ c)))
((not (< c 2))) ; one or two iterations.
(do ((L 1 (1+ L))
(i1 (+ m7 n8)(+ i1 n8)))
((not (< L n9)))
(let ((i2 (+ (aref j L) y)))
(rotatef (aref Z i1)(aref Z i2))
(do ((K 1 (1+ K))
(j2 (+ i2 n8)(+ j2 n8)))
((not (< K L)))
(let ((j1 (+ i1 (aref j K))))
(rotatef (aref Z j1)(aref Z j2))))))))))
(def-bit-reverse bit-reverse-double-float double-float)
(def-bit-reverse bit-reverse-single-float single-float)
(defun check-fft-args(v direction)
(let* ((n (length V))
(p (1- (integer-length n))))
(unless (= 1 (logcount n))
(error "Argument vector length must be a power of two, not ~d." n))
(unless (member direction '(:forward :inverse))
(error "Bad argument <~a> is not (member :forward :inverse))" direction))
p))
(defmacro %%cfft-guts(V direction type)
;; This is the generic guts of the complex radix two decimation in time
;; FFT algorithm. It is a macro so that a few key declarations can be made.
;; Caller must ensure argument consistency.
;; V a vector of length (ash 1 P)
;; DIRECTION is :forward or :inverse
;; TYPE the underlying float type.
`(flet ((%%bfly (V x y W)
(declare (type fft-index x y)
(type (simple-array (complex ,type) (*)) V)
(type (complex ,type) W))
(let* ((xc (aref V x))
(yc (aref V y))
(u (* yc W))
(b (- xc u))
(a (+ xc u)))
(setf (aref V x) a
(aref V y) b))
(values))
(%%cis (Z) ;; saves 240 bytes @N 1024 - 1/23/25
(declare (,type Z))
(complex (cos Z)(sin Z)))
(%%cis-arg (N direction type)
(let ((value (if (eq direction :forward)
(/ pi (float (- N)))
(/ PI (float N)))))
(coerce value type))))
(let* ((N (length ,V))
(P (1- (integer-length N))))
(declare (type fft-index N)
(type (fixnum 0) P))
(do* ((L 0 (1+ L))
(le 2 (ash le 1)))
((not (< L P)))
(declare (type fft-index L le))
(let* ((le2 (ash le -1))
(U (%%cis (coerce 0 ',type)))
(W (%%cis (%%cis-arg le2 ,direction ',type))))
(declare (type (complex ,type) U W)
(type fft-index le2))
(dotimes (j le2)
(declare (type fft-index j))
(let* ((x j)
(y (+ x le2)))
(loop for i fixnum from j below N by le
do (let ()
(%%bfly V x y U)
(incf x le)
(incf y le)))
(setq U (* W U)))))))))
(macrolet
((frob (type)
(let ((name1 (intern (format nil "%CFFT-KERNEL-~a" type)))
(name2 (intern (format nil "CFFT-~a" type)))
(name3 (intern (format nil "BIT-REVERSE-~a" type))))
`(progn
(defun ,name1 (V direction &key bit-reverse-p)
"CFFT on a bit-reverse permuted complex vector"
(declare (type (simple-array (complex ,type) (*)) V)
(type (member :forward :inverse) direction)
(optimize (space 0)(safety 0)(float 0))) ; !
(check-fft-args v direction)
(when bit-reverse-p
(,name3 V))
(%%cfft-guts V direction ,type))
(defun ,name2 (V direction)
"CFFT on a complex vector."
(declare (type (simple-array (complex ,type)(*)) V)
(optimize (float 0)(safety 0)))
(,name1 V direction :bit-reverse-p t)
(when (eq direction :inverse)
(let* ((N (length V))
(f (/ 1d0 (float N)))
(s (coerce f ',type)))
(declare (double-float f)
(,type s))
(dotimes (i N)
(declare (type fft-index i))
(setf (aref V i)(* (aref V i) s)))))
V)))))
(frob double-float)
(frob single-float))