;; large structure tests (require libc.scm) (set! (*s7* 'max-vector-length) (ash 1 36)) (set! (*s7* 'max-string-length) (ash 1 36)) (set! (*s7* 'safety) -1) ;; setting heap-size slows us down (define bold-text (format #f "~C[1m" #\escape)) (define unbold-text (format #f "~C[22m" #\escape)) (define-macro (test a b) (format *stderr* "~S~%" a) `(if (not (equal? ,a ,b)) (format *stderr* " ~A~S -> ~S?~A~%" bold-text ',a ,a unbold-text))) (define (clear-and-gc) (let ((x 0+i)) (do ((i 0 (+ i 1))) ((= i 256)) (set! x (complex i i)))) ; clear temps (gc) (gc)) (define total-memory (with-let *libc* (* (sysconf _SC_PHYS_PAGES) (sysconf _SC_PAGESIZE)))) ; not quite what we want, but... ;; there's no point in overallocating and looking for malloc->null -- malloc doesn't work that way in Linux (define big-size 2500000000) (define fft-size (ash 1 17)) (define little-size 1000000) (unless (defined? 'complex-vector) (define complex-vector vector) (define* (make-complex-vector len (init 0.0)) (make-vector len init complex?)) (define complex-vector-ref vector-ref) (define complex-vector-set! vector-set!) (define complex-vector? vector?) (set! *#readers* (cons (cons #\c (lambda (s) (apply vector (read)))) *#readers*))) ;; -------------------------------------------------------------------------------- (format () "complex fft...~%") (define (complex-fft data n dir) (when data (do ((i 0 (+ i 1)) (j 0)) ((= i n)) (if (> j i) (let ((temp (complex-vector-ref data j))) (complex-vector-set! data j (complex-vector-ref data i)) (complex-vector-set! data i temp))) (let ((m (/ n 2))) (do () ((not (<= 2 m j))) (set! j (- j m)) (set! m (/ m 2))) (set! j (+ j m)))) (let ((ipow (floor (log n 2))) (prev 1)) (do ((lg 0 (+ lg 1)) (mmax 2 (* mmax 2)) (pow (/ n 2) (/ pow 2)) (theta (complex 0.0 (* pi dir)) (* theta 0.5))) ((= lg ipow)) (let ((wpc (exp theta)) (wc 1.0) (tc 0.0)) (do ((ii 0 (+ ii 1))) ((= ii prev)) (do ((jj 0 (+ jj 1)) (i ii (+ i mmax)) (j (+ ii prev) (+ j mmax))) ((>= jj pow)) (set! tc (* wc (complex-vector-ref data j))) (complex-vector-set! data j (- (complex-vector-ref data i) tc)) (complex-vector-set! data i (+ (complex-vector-ref data i) tc))) (set! wc (* wc wpc))) (set! prev mmax)))) data)) (define (complex-checker fvr) (when fvr (let ((pk 0.0) (mx 0.0) (mxloc 0)) (do ((i 0 (+ i 1))) ((= i fft-size)) (set! pk (magnitude (complex-vector-ref fvr i))) (when (> pk mx) (set! mx pk) (set! mxloc i))) (format () "~A ~A~%" mxloc (/ (* 2 (sqrt mx)) fft-size)) (format () "~A ~A~%" (magnitude (fvr 4)) (magnitude (fvr (- fft-size 4))))) (complex-vector-set! fvr 4 0.0) (complex-vector-set! fvr (- fft-size 4) 0.0) (let ((mx 0.0)) (do ((i 0 (+ i 1))) ((= i fft-size)) (set! mx (max mx (magnitude (complex-vector-ref fvr i))))) (format () "noise: ~A~%" mx)))) (define (complex-test) (let ((fvr (make-vector fft-size 0.0 complex?)) (scl (/ (* 8.0 pi) fft-size))) (do ((i 0 (+ i 1)) (x 0.0 (+ x scl))) ((= i fft-size)) (complex-vector-set! fvr i (sin x))) (complex-fft fvr fft-size 1) (complex-checker fvr) (complex-fft #f 0 0) (complex-checker #f))) (complex-test) (clear-and-gc) ;; -------------------------------------------------------------------------------- (format () "ratio fft...~%") (define (ratio-fft rl im n dir) (when rl (do ((i 0 (+ i 1)) (j 0)) ((= i n)) (if (> j i) (let ((tempr (vector-ref rl j)) (tempi (vector-ref im j))) (vector-set! rl j (vector-ref rl i)) (vector-set! im j (vector-ref im i)) (vector-set! rl i tempr) (vector-set! im i tempi))) (let ((m (/ n 2))) (do () ((not (<= 2 m j))) (set! j (- j m)) (set! m (/ m 2))) (set! j (+ j m)))) (let ((ipow (floor (log n 2))) (prev 1)) (do ((lg 0 (+ lg 1)) (mmax 2 (* mmax 2)) (pow (/ n 2) (/ pow 2)) (theta (rationalize (* pi dir)) (* theta 1/2))) ((= lg ipow)) (let ((wpr (rationalize (cos theta))) (wpi (rationalize (sin theta))) (wr 1) (wi 0)) (do ((ii 0 (+ ii 1))) ((= ii prev)) (do ((tempr 0) (tempi 0) (jj 0 (+ jj 1)) (i ii (+ i mmax)) (j (+ ii prev) (+ j mmax))) ((>= jj pow)) (set! tempr (rationalize (- (* wr (vector-ref rl j)) (* wi (vector-ref im j))))) (set! tempi (rationalize (+ (* wr (vector-ref im j)) (* wi (vector-ref rl j))))) (vector-set! rl j (rationalize (- (vector-ref rl i) tempr))) (vector-set! rl i (rationalize (+ (vector-ref rl i) tempr))) (vector-set! im j (rationalize (- (vector-ref im i) tempi))) (vector-set! im i (rationalize (+ (vector-ref im i) tempi)))) (let ((wtemp wr)) (set! wr (rationalize (- (* wr wpr) (* wi wpi)))) (set! wi (rationalize (+ (* wi wpr) (* wtemp wpi)))))) (set! prev mmax)))) rl)) (define (ratio-checker fvr fvi) (when fvr (let ((vr 0.0) (vi 0.0) (pk 0.0) (mx 0.0) (mxloc 0)) (do ((i 0 (+ i 1))) ((= i fft-size)) (set! vr (fvr i)) (set! vi (fvi i)) (set! pk (+ (* vr vr) (* vi vi))) (when (> pk mx) (set! mx pk) (set! mxloc i))) (format () "~A ~A~%" mxloc (/ (* 2 (sqrt mx)) fft-size)) (format () "~A ~A~%" (* 1.0 (fvi 4)) (* 1.0 (fvi (- fft-size 4))))) (set! (fvr 4) 0) (set! (fvi 4) 0) (set! (fvr (- fft-size 4)) 0) (set! (fvi (- fft-size 4)) 0) (let ((mx 0.0)) (do ((i 0 (+ i 1))) ((= i fft-size)) (set! mx (max mx (abs (vector-ref fvr i)) (abs (vector-ref fvi i))))) (format () "noise: ~A~%" (* 1.0 mx))))) (define (ratio-test) (let ((fvr (make-vector fft-size 0 rational?)) (fvi (make-vector fft-size 0 rational?)) (scl (/ (* 8.0 pi) fft-size))) (do ((i 0 (+ i 1)) (x 0.0 (+ x scl))) ((= i fft-size)) (vector-set! fvr i (rationalize (sin x)))) (ratio-fft fvr fvi fft-size 1) (ratio-checker fvr fvi) (ratio-fft #f #f fft-size 1) (ratio-checker #f #f))) (ratio-test) (clear-and-gc) ;; -------------------------------------------------------------------------------- (format () "strings...~%") (when (> total-memory (* 4 big-size)) (clear-and-gc) (let ((bigstr (make-string big-size))) (define (big-string-filler) ;(fill! bigstr #\null 0 big-size)) (do ((i 0 (+ i 1))) ((= i big-size)) (string-set! bigstr i #\null))) (test (length bigstr) big-size) (string-set! bigstr (- big-size 10000000) #\a) (test (string-ref bigstr (- big-size 10000000)) #\a) (set! bigstr (reverse! bigstr)) (test (string-ref bigstr (- 10000000 1)) #\a) (fill! bigstr #\b) (test (string-ref bigstr (- big-size 10000000)) #\b) (test (char-position #\a bigstr) #f) (do ((i (- (ash 1 31) 10) (+ i 1))) ((= i (+ (ash 1 31) 10))) (string-set! bigstr i #\c)) (test (string-ref bigstr (ash 1 31)) #\c) (do ((i (- (ash 1 31) 10) (+ i 1))) ((= i (+ (ash 1 31) 10))) (set! (bigstr i) #\d)) (test (bigstr (ash 1 31)) #\d) (let-temporarily (((*s7* 'print-length) 100)) (let ((big1 (copy bigstr))) (test (equal? big1 bigstr) #t) (test (string=? big1 bigstr) #t) (test (string>? big1 bigstr) #f) ;(test (string-ci<=? big1 bigstr) #t) (test (equivalent? big1 bigstr) #t) (set! (big1 (- big-size 1)) #\f) (test (equal? big1 bigstr) #f) (test (string=? big1 bigstr) #f) (test (equivalent? big1 bigstr) #f) (let ((big2 (substring big1 (- (ash 1 31) 10) (+ (ash 1 31) 10)))) (test big2 "dddddddddddddddddddd") (test (length big2) 20) (let ((big3 (make-string 20))) (copy bigstr big3 (- (ash 1 31) 10) (+ (ash 1 31) 10)) (test (string=? big2 big3) #t))))) (let ((p (call-with-output-string (lambda (p) (write-string bigstr p (- (ash 1 31) 10) (+ (ash 1 31) 10)))))) (test (length p) 20) (test p "dddddddddddddddddddd")) (big-string-filler) ; -- why the do loop? (fill! bigstr #\null 0 big-size) (let ((bv (string->byte-vector bigstr))) (test (byte-vector? bv) #t) (test (length bv) (length bigstr)) (let ((bvi (byte-vector-ref bv (- (ash 1 31) 100)))) (test bvi (char->integer (string-ref (byte-vector->string bv) (- (ash 1 31) 100))))))) (clear-and-gc) (let ((size little-size)) (let ((vs (make-vector size))) (let ((strs (make-vector 256))) (do ((i 0 (+ i 1))) ((= i 256)) (vector-set! strs i (make-string 200 (integer->char i)))) (define (string-sorter) (do ((i 0 (+ i 1))) ((= i size)) (vector-set! vs i (copy (vector-ref strs (logand i #xff)) (make-string (+ 1 (random 200))))))) (string-sorter) (set! strs #f)) (sort! vs string? (vector-ref vs j) (vector-ref vs i)) (display "oops")))) (string-checker) (fill! vs #f))) (clear-and-gc) (let ((size little-size)) (let ((str (make-string size))) (define (string-filler) (do ((i 0 (+ i 1))) ((= i size)) (string-set! str i #\a))) (string-filler) (let ((str1 (make-string size))) (define (string-writer) (do ((i 0 (+ i 1)) (j (integer->char (random 256)) (integer->char (random 256)))) ((= i size)) (string-set! str i (string-ref (make-string (+ 1 (random 200)) j) 0)) (string-set! str1 i j))) (string-writer) (test (string=? str str1) #t)))) (clear-and-gc) (let ((bigstr (make-string (* 2 big-size)))) (test (length bigstr) (* 2 big-size)) (string-set! bigstr (- big-size 10000000) #\a) (test (string-ref bigstr (- big-size 10000000)) #\a) (do ((i (- (ash 1 32) 10) (+ i 1))) ((= i (+ (ash 1 32) 10))) (string-set! bigstr i #\e)) (test (string-ref bigstr (ash 1 32)) #\e)) (clear-and-gc) (let ((rdstr (call-with-input-string (make-string big-size #\null) (lambda (p) (read-string (- big-size 10000000) p))))) (test (length rdstr) (- big-size 10000000)) (test (string-ref rdstr 0) #\null)) (clear-and-gc) (let ((wrtstr (call-with-output-string (lambda (p) (do ((i 0 (+ i 1))) ((>= i 25)) (write-string (make-string 100000000 (integer->char (+ 32 i))) p)))))) (do ((i 0 (+ i 100000000)) (j 0 (+ j 1))) ((>= i big-size)) (if (not (char=? (string-ref wrtstr i) (integer->char (+ 32 j)))) (format () "~D: ~C (~C)~%" i (string-ref wrstr i) (integer->char (+ 32 j))))))) (when (> total-memory (* 4 big-size)) (clear-and-gc) (let ((bigstr1 (make-string big-size #\+)) (bigstr2 (make-string big-size #\-))) (let ((bigstr3 (append bigstr1 bigstr2))) (test (length bigstr3) (* 2 big-size)) (test (string-ref bigstr3 (- big-size 1)) #\+) (test (string-ref bigstr3 (+ big-size 1)) #\-) (test (char-position #\- bigstr3) big-size)))) (clear-and-gc) (let-temporarily (((*s7* 'max-format-length) (+ big-size 1))) (let ((x (format #f "~NC" big-size #\a))) (test (length x) big-size) (test (string-ref x 1000000) #\a))) (clear-and-gc) ;; -------------------------------------------------------------------------------- (format () "~%byte-vectors...~%") (when (> total-memory (* 4 big-size)) (let ((bigstr (make-byte-vector big-size))) (test (length bigstr) big-size) (byte-vector-set! bigstr (- big-size 10000000) 65) (test (byte-vector-ref bigstr (- big-size 10000000)) 65) (set! bigstr (reverse! bigstr)) (test (byte-vector-ref bigstr (- 10000000 1)) 65) (fill! bigstr 66) (test (byte-vector-ref bigstr (- big-size 10000000)) 66)) (clear-and-gc) (let ((size little-size)) (let ((vs (make-byte-vector size))) (define (byte-vector-sorter) (do ((i 0 (+ i 1))) ((= i size)) (byte-vector-set! vs i (random 256)))) (byte-vector-sorter) (sort! vs <) (define (byte-vector-checker) (let ((vj (byte-vector-ref vs 0))) (for-each (lambda (vi) (if (> vj vi) (display "oops")) (set! vj vi)) (subvector vs 1 size)))) (byte-vector-checker))) (clear-and-gc) (let ((size little-size)) (let ((str (make-byte-vector size))) (define (byte-vector-filler) (do ((i 0 (+ i 1))) ((= i size)) (byte-vector-set! str i 65))) (byte-vector-filler) (let ((str1 (make-byte-vector size))) (define (byte-vector-writer) (do ((i 0 (+ i 1)) (j (random 256) (random 256))) ((= i size)) (byte-vector-set! str i (byte-vector-ref (make-byte-vector (+ 1 (random 200)) j) 0)) (byte-vector-set! str1 i j))) (byte-vector-writer) (test (equal? str str1) #t))))) ;; -------------------------------------------------------------------------------- (format () "~%float-vectors...~%") (when (> total-memory (* 18 big-size)) (format () "test 1~%") (clear-and-gc) (let ((bigfv (make-float-vector big-size 0.5))) (let-temporarily (((*s7* 'print-length) 100)) (let ((big1 (copy bigfv))) (test (equivalent? big1 bigfv) #t) (set! (big1 (- big-size 1)) 0.25) (test (equivalent? big1 bigfv) #f) (let ((big2 (subvector big1 (- (ash 1 31) 10) (+ (ash 1 31) 10)))) (test big2 (make-float-vector 20 0.5)) (test (length big2) 20) (let ((big3 (make-float-vector 20 0.0))) (copy bigfv big3 (- (ash 1 31) 10) (+ (ash 1 31) 10)) (test (equivalent? big2 big3) #t))))) (define (big-float-vector-filler) (do ((i 0 (+ i 1))) ((= i big-size)) (float-vector-set! bigfv i 1.0))) (big-float-vector-filler) (test (bigfv 1) 1.0))) (when (> total-memory (* 16 little-size)) (format () "test 2~%") (clear-and-gc) (let ((size little-size)) (let ((vs (make-float-vector size))) (define (float-vector-sorter) (do ((i 0 (+ i 1))) ((= i size)) (float-vector-set! vs i (float-vector-ref (make-float-vector (+ 1 (random 200)) (random 100.0)) 0)))) (float-vector-sorter) (sort! vs <) (define (float-vector-checker) (let ((vj (float-vector-ref vs 0))) (for-each (lambda (vi) (if (> vj vi) (display "oops")) (set! vj vi)) (subvector vs 1 size)))) (float-vector-checker))) (clear-and-gc) (let ((size little-size)) (let ((str (make-float-vector size))) (define (float-vector-filler) (do ((i 0 (+ i 1))) ((= i size)) (float-vector-set! str i 5.0))) (float-vector-filler) (let ((str1 (make-float-vector size))) (define (float-vector-writer) (do ((i 0 (+ i 1)) (j (random 100.0) (random 100.0))) ((= i size)) (float-vector-set! str i (float-vector-ref (make-float-vector (+ 1 (random 200)) j) 0)) (float-vector-set! str1 i j))) (float-vector-writer) (test (equivalent? str str1) #t))))) (when (> total-memory (* 16 big-size)) (format () "test 3~%") (clear-and-gc) (let ((bigfv1 (make-float-vector (/ big-size 2) 1.0)) (bigfv2 (make-float-vector (/ big-size 2) 2.0))) (let ((bigfv3 (append bigfv1 bigfv2))) (test (length bigfv3) big-size) (test (float-vector-ref bigfv3 (- (/ big-size 2) 1)) 1.0) (test (float-vector-ref bigfv3 (+ (/ big-size 2) 1)) 2.0)))) (when (> total-memory (* 32 big-size)) (format () "test 3a~%") (clear-and-gc) (let ((bigfv1 (make-float-vector (list 2 (/ big-size 2)) 1.0)) (bigfv2 (make-float-vector (list 2 (/ big-size 2)) 2.0))) (test (float-vector-ref bigfv1 0 (- (/ big-size 2) 1)) 1.0) (test (float-vector-ref bigfv2 1 (- (/ big-size 2) 1)) 2.0) (let ((bigfv3 (append bigfv1 bigfv2))) (test (length bigfv3) (* 2 big-size)) (test (float-vector-ref bigfv3 (- big-size 1)) 1.0) (test (float-vector-ref bigfv3 (+ big-size 1)) 2.0)))) (when (> total-memory (* 32 big-size)) (format () "test 3b~%") (clear-and-gc) (let ((bigfv1 (make-float-vector (* big-size 2) 1.0)) (bigfv2 (make-float-vector (* big-size 2) 2.0))) (test (float-vector-ref bigfv1 (- (/ big-size 2) 1)) 1.0) (test (float-vector-ref bigfv2 (- (/ big-size 2) 1)) 2.0) (define (f-loop) (do ((i (- (ash 1 32) 10) (+ i 1)) (j 0 (+ j 1))) ((= i (+ (ash 1 32) 10))) (set! (bigfv1 i) j) (set! (bigfv2 i) (* (bigfv1 i) 2)))) (f-loop) (test (subvector bigfv1 (- (ash 1 32) 10) (+ (ash 1 32) 10)) (float-vector 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19)) (test (subvector bigfv2 (- (ash 1 32) 10) (+ (ash 1 32) 10)) (float-vector 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38)) (test (subvector bigfv1 (+ (ash 1 32) 5) (+ (ash 1 32) 10)) (float-vector 15 16 17 18 19)))) (when (> total-memory (* 9 big-size)) (format () "test 4~%") (clear-and-gc) (let ((bigfv (make-float-vector big-size))) (test (length bigfv) big-size) (float-vector-set! bigfv (- big-size 10000000) 1.0) (test (float-vector-ref bigfv (- big-size 10000000)) 1.0) (reverse! bigfv) (test (float-vector-ref bigfv (- 10000000 1)) 1.0) (fill! bigfv 0.0) (test (bigfv (- big-size 10000000)) 0.0) (do ((i (- (ash 1 31) 10) (+ i 1))) ((= i (+ (ash 1 31) 10))) (float-vector-set! bigfv i 2.0)) (test (float-vector-ref bigfv (ash 1 31)) 2.0) (do ((i (- (ash 1 31) 10) (+ i 1))) ((= i (+ (ash 1 31) 10))) (set! (bigfv i) pi)) (test (equivalent? (bigfv (ash 1 31)) pi) #t)) ; pi might be bignum (clear-and-gc)) (define (float-vector-fft rl im n dir) (when rl (let ((tempr 0.0) (tempi 0.0) (m 0) (j 0) (ipw (floor (log n 2))) (prev 1) (wtemp 0.0) (wpr 0.0) (wpi 0.0) (wr 0.0) (wi 0.0) (mmax 2) (pw (/ n 2)) (end 0) (theta (* pi dir))) (do ((i 0 (+ i 1))) ((= i n)) (when (> j i) (set! tempr (float-vector-ref rl j)) (set! tempi (float-vector-ref im j)) (float-vector-set! rl j (float-vector-ref rl i)) (float-vector-set! im j (float-vector-ref im i)) (float-vector-set! rl i tempr) (float-vector-set! im i tempi)) (set! m (quotient n 2)) (do () ((not (<= 2 m j))) (set! j (- j m)) (set! m (quotient m 2))) (set! j (+ j m))) (do ((lg 0 (+ lg 1))) ((= lg ipw)) (set! wpr (cos theta)) (set! wpi (sin theta)) (set! wr 1.0) (set! wi 0.0) (do ((ii 0 (+ ii 1))) ((= ii prev)) (set! j (+ ii prev)) (set! end (+ ii (* pw mmax))) (do ((i ii (+ i mmax))) ((= i end)) (set! tempr (- (* wr (float-vector-ref rl j)) (* wi (float-vector-ref im j)))) (set! tempi (+ (* wr (float-vector-ref im j)) (* wi (float-vector-ref rl j)))) (float-vector-set! rl j (- (float-vector-ref rl i) tempr)) (float-vector-set! rl i (+ (float-vector-ref rl i) tempr)) (float-vector-set! im j (- (float-vector-ref im i) tempi)) (float-vector-set! im i (+ (float-vector-ref im i) tempi)) (set! j (+ j mmax))) (set! wtemp wr) (set! wr (- (* wr wpr) (* wi wpi))) (set! wi (+ (* wi wpr) (* wtemp wpi)))) (set! prev mmax) (set! mmax (* mmax 2)) (set! pw (quotient pw 2)) (set! theta (* theta 0.5))) rl))) (define (float-vector-checker fvr fvi) (when fvr (let ((vr 0.0) (vi 0.0) (pk 0.0) (mx 0.0) (mxloc 0)) (do ((i 0 (+ i 1))) ((= i fft-size)) (set! vr (float-vector-ref fvr i)) (set! vi (float-vector-ref fvi i)) (set! pk (+ (* vr vr) (* vi vi))) (when (> pk mx) (set! mx pk) (set! mxloc i))) (format () "~A ~A~%" mxloc (/ (* 2 (sqrt mx)) fft-size)) (format () "~A ~A~%" (fvi 4) (fvi (- fft-size 4)))) (float-vector-set! fvr 4 0.0) (float-vector-set! fvi 4 0.0) (float-vector-set! fvr (- fft-size 4) 0.0) (float-vector-set! fvi (- fft-size 4) 0.0) (let ((mx 0.0)) (do ((i 0 (+ i 1))) ((= i fft-size)) (set! mx (max mx (abs (float-vector-ref fvr i)) (abs (float-vector-ref fvi i))))) (format () "noise: ~A~%" mx)))) (define (float-vector-test) (let ((fvr (make-float-vector fft-size)) (fvi (make-float-vector fft-size 0.0)) (scl (/ (* 8.0 pi) fft-size))) (do ((i 0 (+ i 1)) (x 0.0 (+ x scl))) ((= i fft-size)) (float-vector-set! fvr i (sin x))) (float-vector-fft fvr fvi fft-size 1) (float-vector-checker fvr fvi) (float-vector-fft #f #f fft-size 1) (float-vector-checker #f #f))) (float-vector-test) (clear-and-gc) (define (float-2d-fft rl n dir) (when rl (let ((tempr 0.0) (tempi 0.0) (m 0) (j 0) (ipw (floor (log n 2))) (prev 1) (wtemp 0.0) (wpr 0.0) (wpi 0.0) (wr 0.0) (wi 0.0) (mmax 2) (pw (/ n 2)) (end 0) (theta (* pi dir))) (do ((i 0 (+ i 1))) ((= i n)) (when (> j i) (set! tempr (float-vector-ref rl 0 j)) (set! tempi (float-vector-ref rl 1 j)) (float-vector-set! rl 0 j (float-vector-ref rl 0 i)) (float-vector-set! rl 1 j (float-vector-ref rl 1 i)) (float-vector-set! rl 0 i tempr) (float-vector-set! rl 1 i tempi)) (set! m (quotient n 2)) (do () ((not (<= 2 m j))) (set! j (- j m)) (set! m (quotient m 2))) (set! j (+ j m))) (do ((lg 0 (+ lg 1))) ((= lg ipw)) (set! wpr (cos theta)) (set! wpi (sin theta)) (set! wr 1.0) (set! wi 0.0) (do ((ii 0 (+ ii 1))) ((= ii prev)) (set! j (+ ii prev)) (set! end (+ ii (* pw mmax))) (do ((i ii (+ i mmax))) ((= i end)) (set! tempr (- (* wr (float-vector-ref rl 0 j)) (* wi (float-vector-ref rl 1 j)))) (set! tempi (+ (* wr (float-vector-ref rl 1 j)) (* wi (float-vector-ref rl 0 j)))) (float-vector-set! rl 0 j (- (float-vector-ref rl 0 i) tempr)) (float-vector-set! rl 0 i (+ (float-vector-ref rl 0 i) tempr)) (float-vector-set! rl 1 j (- (float-vector-ref rl 1 i) tempi)) (float-vector-set! rl 1 i (+ (float-vector-ref rl 1 i) tempi)) (set! j (+ j mmax))) (set! wtemp wr) (set! wr (- (* wr wpr) (* wi wpi))) (set! wi (+ (* wi wpr) (* wtemp wpi)))) (set! prev mmax) (set! mmax (* mmax 2)) (set! pw (quotient pw 2)) (set! theta (* theta 0.5))) rl))) (define (float-2d-checker fvr) (when fvr (let ((vr 0.0) (vi 0.0) (pk 0.0) (mx 0.0) (mxloc 0)) (do ((i 0 (+ i 1))) ((= i fft-size)) (set! vr (float-vector-ref fvr 0 i)) (set! vi (float-vector-ref fvr 1 i)) (set! pk (+ (* vr vr) (* vi vi))) (when (> pk mx) (set! mx pk) (set! mxloc i))) (format () "~A ~A~%" mxloc (/ (* 2 (sqrt mx)) fft-size)) (format () "~A ~A~%" (fvr 1 4) (fvr 1 (- fft-size 4)))) (float-vector-set! fvr 0 4 0.0) (float-vector-set! fvr 1 4 0.0) (float-vector-set! fvr 0 (- fft-size 4) 0.0) (float-vector-set! fvr 1 (- fft-size 4) 0.0) (let ((mx 0.0)) (do ((i 0 (+ i 1))) ((= i fft-size)) (set! mx (max mx (abs (float-vector-ref fvr 0 i)) (abs (float-vector-ref fvr 1 i))))) (format () "noise: ~A~%" mx)))) (define (float-2d-test skip) (format *stderr* "float-2d ~A~%" skip) (let ((fvr (make-float-vector (list skip fft-size) 0.0)) (scl (/ (* 8.0 pi) fft-size))) (do ((i 0 (+ i 1)) (x 0.0 (+ x scl))) ((= i fft-size)) (float-vector-set! fvr 0 i (sin x))) (float-2d-fft fvr fft-size 1) (float-2d-checker fvr) (float-2d-fft #f fft-size 1) (float-2d-checker #f))) (float-2d-test 2) (clear-and-gc) (let-temporarily ((fft-size 32)) (float-2d-test 2)) ;; now try a big float-vector, but same size fft (when (> total-memory (* 9 (ash 1 31))) (float-2d-test (ash 1 14))) ; fft-size (ash 1 17) -> (* 8 (ash 1 31)) (let-temporarily ((fft-size 32)) (float-2d-test 2)) (clear-and-gc) ;; -------------------------------------------------------------------------------- (format () "~%int-vectors...~%") (when (> total-memory (* 18 big-size)) (format () "test 1~%") (clear-and-gc) (let ((bigfv (make-int-vector big-size 5))) (let-temporarily (((*s7* 'print-length) 100)) (let ((big1 (copy bigfv))) (test (equivalent? big1 bigfv) #t) (set! (big1 (- big-size 1)) 25) (test (equivalent? big1 bigfv) #f) (let ((big2 (subvector big1 (- (ash 1 31) 10) (+ (ash 1 31) 10)))) (test big2 (make-int-vector 20 5)) (test (length big2) 20) (let ((big3 (make-int-vector 20 0))) (copy bigfv big3 (- (ash 1 31) 10) (+ (ash 1 31) 10)) (test (equivalent? big2 big3) #t))))) (define (big-int-vector-filler) (do ((i 0 (+ i 1))) ((= i big-size)) (int-vector-set! bigfv i 1))) (big-int-vector-filler) (test (bigfv 1) 1))) (when (> total-memory (* 16 little-size)) (format () "test 2~%") (clear-and-gc) (let ((size little-size)) (let ((vs (make-int-vector size))) (define (int-vector-sorter) (do ((i 0 (+ i 1))) ((= i size)) (int-vector-set! vs i (int-vector-ref (make-int-vector (+ 1 (random 200)) (random 100)) 0)))) (int-vector-sorter) (sort! vs <) (define (int-vector-checker) (let ((vj (int-vector-ref vs 0))) (for-each (lambda (vi) (if (> vj vi) (display "oops")) (set! vj vi)) (subvector vs 1 size)))) (int-vector-checker))) (sort! '(1 2 3) <) ; clear temp3... (clear-and-gc) (clear-and-gc) (let ((size little-size)) (let ((str (make-int-vector size))) (define (int-vector-filler) (do ((i 0 (+ i 1))) ((= i size)) (int-vector-set! str i 5))) (int-vector-filler) (let ((str1 (make-int-vector size))) (define (int-vector-writer) (do ((i 0 (+ i 1)) (j (random 100) (random 100))) ((= i size)) (int-vector-set! str i (int-vector-ref (make-int-vector (+ 1 (random 200)) j) 0)) (int-vector-set! str1 i j))) (int-vector-writer) (test (equivalent? str str1) #t))))) (clear-and-gc) (when (> total-memory (* 16 big-size)) (format () "test 3~%") (clear-and-gc) (let ((bigfv1 (make-int-vector (/ big-size 2) 1)) (bigfv2 (make-int-vector (/ big-size 2) 2))) (let ((bigfv3 (append bigfv1 bigfv2))) (test (length bigfv3) big-size) (test (int-vector-ref bigfv3 (- (/ big-size 2) 1)) 1) (test (int-vector-ref bigfv3 (+ (/ big-size 2) 1)) 2)))) (when (> total-memory (* 9 big-size)) (format () "test 4~%") (clear-and-gc) (let ((bigfv (make-int-vector big-size))) (test (length bigfv) big-size) (int-vector-set! bigfv (- big-size 10000000) 1) (test (int-vector-ref bigfv (- big-size 10000000)) 1) (reverse! bigfv) (test (int-vector-ref bigfv (- 10000000 1)) 1) (fill! bigfv 0) (test (bigfv (- big-size 10000000)) 0) (do ((i (- (ash 1 31) 10) (+ i 1))) ((= i (+ (ash 1 31) 10))) (int-vector-set! bigfv i 2)) (test (int-vector-ref bigfv (ash 1 31)) 2) (do ((i (- (ash 1 31) 10) (+ i 1))) ((= i (+ (ash 1 31) 10))) (set! (bigfv i) 3)) (test (bigfv (ash 1 31)) 3)) (clear-and-gc)) (define (int-vector-fft rl im n dir) (when rl (let ((tempr 0.0) (tempi 0.0) (tr 0) (ti 0) (m 0) (j 0) (ipw (floor (log n 2))) (prev 1) (wtemp 0.0) (wpr 0.0) (wpi 0.0) (wr 0.0) (wi 0.0) (mmax 2) (pw (/ n 2)) (end 0) (theta (* pi dir))) (do ((i 0 (+ i 1))) ((= i n)) (when (> j i) (set! tr (int-vector-ref rl j)) (set! ti (int-vector-ref im j)) (int-vector-set! rl j (int-vector-ref rl i)) (int-vector-set! im j (int-vector-ref im i)) (int-vector-set! rl i tr) (int-vector-set! im i ti)) (set! m (quotient n 2)) (do () ((not (<= 2 m j))) (set! j (- j m)) (set! m (quotient m 2))) (set! j (+ j m))) (do ((lg 0 (+ lg 1))) ((= lg ipw)) (set! wpr (cos theta)) (set! wpi (sin theta)) (set! wr 1.0) (set! wi 0.0) (do ((ii 0 (+ ii 1))) ((= ii prev)) (set! j (+ ii prev)) (set! end (+ ii (* pw mmax))) (do ((i ii (+ i mmax))) ((= i end)) (set! tempr (- (* wr (int-vector-ref rl j)) (* wi (int-vector-ref im j)))) (set! tempi (+ (* wr (int-vector-ref im j)) (* wi (int-vector-ref rl j)))) (int-vector-set! rl j (round (- (int-vector-ref rl i) tempr))) (int-vector-set! rl i (round (+ (int-vector-ref rl i) tempr))) (int-vector-set! im j (round (- (int-vector-ref im i) tempi))) (int-vector-set! im i (round (+ (int-vector-ref im i) tempi))) (set! j (+ j mmax))) (set! wtemp wr) (set! wr (- (* wr wpr) (* wi wpi))) (set! wi (+ (* wi wpr) (* wtemp wpi)))) (set! prev mmax) (set! mmax (* mmax 2)) (set! pw (quotient pw 2)) (set! theta (* theta 0.5))) rl))) (define (int-vector-checker fvr fvi) (when fvr (let ((vr 0) (vi 0) (pk 0) (mx 0) (mxloc 0)) (do ((i 0 (+ i 1))) ((= i fft-size)) (set! vr (int-vector-ref fvr i)) (set! vi (int-vector-ref fvi i)) (set! pk (+ (* vr vr) (* vi vi))) (when (> pk mx) (set! mx pk) (set! mxloc i))) (format () "~A ~A~%" mxloc (/ (* 2 (sqrt mx)) fft-size)) (format () "~A ~A~%" (fvi 4) (fvi (- fft-size 4)))) (int-vector-set! fvr 4 0) (int-vector-set! fvi 4 0) (int-vector-set! fvr (- fft-size 4) 0) (int-vector-set! fvi (- fft-size 4) 0) (let ((mx 0.0)) (do ((i 0 (+ i 1))) ((= i fft-size)) (set! mx (max mx (abs (int-vector-ref fvr i)) (abs (int-vector-ref fvi i))))) (format () "noise: ~A~%" mx)))) (define (int-vector-test) (let ((fvr (make-int-vector fft-size)) (fvi (make-int-vector fft-size 0)) (scl (/ (* 8.0 pi) fft-size))) (do ((i 0 (+ i 1)) (x 0.0 (+ x scl))) ((= i fft-size)) (int-vector-set! fvr i (round (* 1000.0 (sin x))))) (int-vector-fft fvr fvi fft-size 1) (int-vector-checker fvr fvi) (int-vector-fft #f #f fft-size 1) (int-vector-checker #f #f))) (int-vector-test) (clear-and-gc) ;; -------------------------------------------------------------------------------- (format () "~%vectors...~%") (when (> total-memory (* 9 big-size)) (let ((bigv (make-vector big-size 'a symbol?))) (test (length bigv) big-size) (vector-set! bigv (- big-size 10000000) 'asdf) (test (vector-ref bigv (- big-size 10000000)) 'asdf) (reverse! bigv) (test (vector-ref bigv (- 10000000 1)) 'asdf) (fill! bigv 'b) (test (vector-ref bigv (- 10000000 1)) 'b) )) (clear-and-gc) (define (v-loop-test) (when (> total-memory (* 20 big-size)) (let ((v (make-vector big-size ()))) (do ((i 0 (+ i 1000))) ((= i big-size)) (vector-set! v i (* 2 i))) (test (vector-ref v 100000000) (* 2 100000000)) (test (vector-ref v (- big-size 10000000)) (* 2 (- big-size 10000000)))))) (v-loop-test) (clear-and-gc) (define (vector-fft rl im n dir) (when rl (let ((m 0) (j 0) (ipw (floor (log n 2))) (prev 1) (wtemp 0.0) (wpr 0.0) (wpi 0.0) (wr 0.0) (wi 0.0) (mmax 2) (pw (/ n 2)) (end 0) (theta (* pi dir))) (do ((tempr 0.0+i) (tempi 0.0+i) (i 0 (+ i 1))) ((= i n)) (when (> j i) (set! tempr (vector-ref rl j)) (set! tempi (vector-ref im j)) (vector-set! rl j (vector-ref rl i)) (vector-set! im j (vector-ref im i)) (vector-set! rl i tempr) (vector-set! im i tempi)) (set! m (quotient n 2)) (do () ((not (<= 2 m j))) (set! j (- j m)) (set! m (quotient m 2))) (set! j (+ j m))) (do ((lg 0 (+ lg 1))) ((= lg ipw)) (set! wpr (cos theta)) (set! wpi (sin theta)) (set! wr 1.0) (set! wi 0.0) (do ((ii 0 (+ ii 1))) ((= ii prev)) (set! j (+ ii prev)) (set! end (+ ii (* pw mmax))) (do ((tempr 0.0+i) (tempi 0.0+i) (i ii (+ i mmax))) ((= i end)) (set! tempr (- (* wr (vector-ref rl j)) (* wi (vector-ref im j)))) (set! tempi (+ (* wr (vector-ref im j)) (* wi (vector-ref rl j)))) (vector-set! rl j (- (vector-ref rl i) tempr)) (vector-set! rl i (+ (vector-ref rl i) tempr)) (vector-set! im j (- (vector-ref im i) tempi)) (vector-set! im i (+ (vector-ref im i) tempi)) (set! j (+ j mmax))) (set! wtemp wr) (set! wr (- (* wr wpr) (* wi wpi))) (set! wi (+ (* wi wpr) (* wtemp wpi)))) (set! prev mmax) (set! mmax (* mmax 2)) (set! pw (quotient pw 2)) (set! theta (* theta 0.5))) rl))) (define (vector-checker fvr fvi) (when fvr (let ((vr 0.0) (vi 0.0) (pk 0.0) (mx 0.0) (mxloc 0)) (do ((i 0 (+ i 1))) ((= i fft-size)) (set! vr (vector-ref fvr i)) (set! vi (vector-ref fvi i)) (set! pk (+ (* vr vr) (* vi vi))) (when (> pk mx) (set! mx pk) (set! mxloc i))) (format () "~A ~A~%" mxloc (/ (* 2 (sqrt mx)) fft-size)) (format () "~A ~A~%" (fvi 4) (fvi (- fft-size 4)))) (vector-set! fvr 4 0.0) (vector-set! fvi 4 0.0) (vector-set! fvr (- fft-size 4) 0.0) (vector-set! fvi (- fft-size 4) 0.0) (let ((mx 0.0)) (do ((i 0 (+ i 1))) ((= i fft-size)) (set! mx (max mx (abs (vector-ref fvr i)) (abs (vector-ref fvi i))))) (format () "noise: ~A~%" mx)))) (define (vector-test) (let ((fvr (make-vector fft-size 0.0 real?)) (fvi (make-vector fft-size 0.0 real?)) (scl (/ (* 8.0 pi) fft-size))) (do ((i 0 (+ i 1)) (x 0.0 (+ x scl))) ((= i fft-size)) (vector-set! fvr i (sin x))) (vector-fft fvr fvi fft-size 1) (vector-checker fvr fvi) (vector-fft #f #f fft-size 1) (vector-checker #f #f))) (vector-test) (clear-and-gc) (define (vector-2d-fft rl n dir) (when rl (let ((m 0) (j 0) (ipw (floor (log n 2))) (prev 1) (wtemp 0.0) (wpr 0.0) (wpi 0.0) (wr 0.0) (wi 0.0) (mmax 2) (pw (/ n 2)) (end 0) (theta (* pi dir)) (rl0 (subvector rl 0 fft-size)) (rl1 (subvector rl fft-size (* 2 fft-size)))) (do ((tempr 0.0+i) (tempi 0.0+i) (i 0 (+ i 1))) ((= i n)) (when (> j i) (set! tempr (vector-ref rl0 j)) (set! tempi (vector-ref rl1 j)) (vector-set! rl0 j (vector-ref rl0 i)) (vector-set! rl1 j (vector-ref rl1 i)) (vector-set! rl0 i tempr) (vector-set! rl1 i tempi)) (set! m (quotient n 2)) (do () ((not (<= 2 m j))) (set! j (- j m)) (set! m (quotient m 2))) (set! j (+ j m))) (do ((lg 0 (+ lg 1))) ((= lg ipw)) (set! wpr (cos theta)) (set! wpi (sin theta)) (set! wr 1.0) (set! wi 0.0) (do ((tempr 0.0+i) (tempi 0.0+i) (ii 0 (+ ii 1))) ((= ii prev)) (set! j (+ ii prev)) (set! end (+ ii (* pw mmax))) (do ((i ii (+ i mmax))) ((= i end)) (set! tempr (- (* wr (vector-ref rl0 j)) (* wi (vector-ref rl1 j)))) (set! tempi (+ (* wr (vector-ref rl1 j)) (* wi (vector-ref rl0 j)))) (vector-set! rl0 j (- (vector-ref rl0 i) tempr)) (vector-set! rl0 i (+ (vector-ref rl0 i) tempr)) (vector-set! rl1 j (- (vector-ref rl1 i) tempi)) (vector-set! rl1 i (+ (vector-ref rl1 i) tempi)) (set! j (+ j mmax))) (set! wtemp wr) (set! wr (- (* wr wpr) (* wi wpi))) (set! wi (+ (* wi wpr) (* wtemp wpi)))) (set! prev mmax) (set! mmax (* mmax 2)) (set! pw (quotient pw 2)) (set! theta (* theta 0.5))) rl))) (define (vector-2d-checker fvr) (when fvr (let ((vr 0.0) (vi 0.0) (pk 0.0) (mx 0.0) (mxloc 0)) (do ((i 0 (+ i 1))) ((= i fft-size)) (set! vr (vector-ref fvr 0 i)) (set! vi (vector-ref fvr 1 i)) (set! pk (+ (* vr vr) (* vi vi))) (when (> pk mx) (set! mx pk) (set! mxloc i))) (format () "~A ~A~%" mxloc (/ (* 2 (sqrt mx)) fft-size)) (format () "~A ~A~%" (fvr 1 4) (fvr 1 (- fft-size 4)))) (vector-set! fvr 0 4 0.0) (vector-set! fvr 1 4 0.0) (vector-set! fvr 0 (- fft-size 4) 0.0) (vector-set! fvr 1 (- fft-size 4) 0.0) (let ((mx 0.0)) (do ((i 0 (+ i 1))) ((= i fft-size)) (set! mx (max mx (abs (vector-ref fvr 0 i)) (abs (vector-ref fvr 1 i))))) (format () "noise: ~A~%" mx)))) (define (vector-2d-test skip) (let ((fvr (make-vector (list skip fft-size) 0.0 real?)) (scl (/ (* 8.0 pi) fft-size))) (do ((i 0 (+ i 1)) (x 0.0 (+ x scl))) ((= i fft-size)) (vector-set! fvr 0 i (sin x))) (vector-2d-fft fvr fft-size 1) (vector-2d-checker fvr) (vector-2d-fft #f fft-size 1) (vector-2d-checker #f))) (vector-2d-test 2) (clear-and-gc) ;; -------------------------------------------------------------------------------- ;; hash-tables round up to the next power of 2 (format () "~%hash-tables...~%") (when (> total-memory (* 9 (ash 1 31))) (clear-and-gc) (set! big-size 2000000000) (let ((bigv (make-hash-table big-size))) (test (length bigv) 2147483648) (hash-table-set! bigv 'asdf 12) (test (hash-table-ref bigv 'asdf) 12) (hash-table-set! bigv (- big-size 10000000) 'asdf) (test (hash-table-ref bigv (- big-size 10000000)) 'asdf) (fill! bigv #f) (test (hash-table-entries bigv) 0) (test (length bigv) 2147483648)) (clear-and-gc) ) (when (> total-memory (* 9 (ash 1 32))) ; add some slack (set! big-size 2500000000) (let ((bigv (make-hash-table big-size))) (test (length bigv) 4294967296) (hash-table-set! bigv 'asdf 12) (test (hash-table-ref bigv 'asdf) 12) (hash-table-set! bigv (- big-size 10000000) 'asdf) (test (hash-table-ref bigv (- big-size 10000000)) 'asdf) ) (clear-and-gc) ) (define (hash-table-fft rl im n dir) (when rl (do ((i 0 (+ i 1)) (j 0)) ((= i n)) (if (> j i) (let ((tempr (hash-table-ref rl j)) (tempi (hash-table-ref im j))) (hash-table-set! rl j (hash-table-ref rl i)) (hash-table-set! im j (hash-table-ref im i)) (hash-table-set! rl i tempr) (hash-table-set! im i tempi))) (let ((m (/ n 2))) (do () ((not (<= 2 m j))) (set! j (- j m)) (set! m (/ m 2))) (set! j (+ j m)))) (let ((ipow (floor (log n 2))) (prev 1)) (do ((mmax 2 (* mmax 2)) (pow (/ n 2) (/ pow 2)) (theta (* pi dir) (* theta 0.5)) (lg 0 (+ lg 1)) (wpr 0.0) (wpi 0.0) (wr 1.0 1.0) (wi 0.0 0.0) (wtemp 0.0+i)) ((= lg ipow)) (set! wpr (cos theta)) (set! wpi (sin theta)) (do ((ii 0 (+ ii 1))) ((= ii prev)) (do ((tempr 0.0+i) (tempi 0.0+i) (i ii (+ i mmax)) (j (+ ii prev) (+ j mmax)) (jj 0 (+ jj 1))) ((>= jj pow)) (set! tempr (- (* wr (hash-table-ref rl j)) (* wi (hash-table-ref im j)))) (set! tempi (+ (* wr (hash-table-ref im j)) (* wi (hash-table-ref rl j)))) (hash-table-set! rl j (- (hash-table-ref rl i) tempr)) (hash-table-set! rl i (+ (hash-table-ref rl i) tempr)) (hash-table-set! im j (- (hash-table-ref im i) tempi)) (hash-table-set! im i (+ (hash-table-ref im i) tempi))) (set! wtemp wr) (set! wr (- (* wr wpr) (* wi wpi))) (set! wi (+ (* wi wpr) (* wtemp wpi)))) (set! prev mmax))) rl)) (define (hash-test) (let ((fvr (make-hash-table fft-size)) ; (make-hash-table fft-size #t (cons #t float?)) (fvi (make-hash-table fft-size)) (scl (/ (* 8.0 pi) fft-size)) (vr 0+i) (vi 0+i) (pk 0+i)) (do ((i 0 (+ i 1)) (x 0.0 (+ x scl))) ((= i fft-size)) (hash-table-set! fvr i (sin x)) (hash-table-set! fvi i 0.0)) (hash-table-fft fvr fvi fft-size 1) (do ((mx 0.0) (mxloc 0) (i 0 (+ i 1))) ((= i fft-size) (format () "~A ~A~%" mxloc (/ (* 2 (sqrt mx)) fft-size))) (set! vr (hash-table-ref fvr i)) (set! vi (hash-table-ref fvi i)) (set! pk (+ (* vr vr) (* vi vi))) (when (> pk mx) (set! mx pk) (set! mxloc i))) (hash-table-set! fvr 4 0.0) (hash-table-set! fvi 4 0.0) (hash-table-set! fvr (- fft-size 4) 0.0) (hash-table-set! fvi (- fft-size 4) 0.0) (do ((mx 0.0) (i 0 (+ i 1))) ((= i fft-size) (format () "noise: ~A~%" mx)) (set! mx (max mx (abs (hash-table-ref fvr i)) (abs (hash-table-ref fvi i))))) (hash-table-fft #f #f fft-size 1))) (hash-test) (clear-and-gc) (define (hash-symbol-fft rl im n rl-syms im-syms dir) (when rl (do ((i 0 (+ i 1)) (j 0)) ((= i n)) (if (> j i) (let ((tempr (hash-table-ref rl (vector-ref rl-syms j))) (tempi (hash-table-ref im (vector-ref im-syms j)))) (hash-table-set! rl (vector-ref rl-syms j) (hash-table-ref rl (vector-ref rl-syms i))) (hash-table-set! im (vector-ref im-syms j) (hash-table-ref im (vector-ref im-syms i))) (hash-table-set! rl (vector-ref rl-syms i) tempr) (hash-table-set! im (vector-ref im-syms i) tempi))) (let ((m (/ n 2))) (do () ((not (<= 2 m j))) (set! j (- j m)) (set! m (/ m 2))) (set! j (+ j m)))) (let ((ipow (floor (log n 2))) (prev 1)) (do ((mmax 2 (* mmax 2)) (pow (/ n 2) (/ pow 2)) (theta (* pi dir) (* theta 0.5)) (lg 0 (+ lg 1)) (wpr 0.0) (wpi 0.0) (wr 1.0 1.0) (wi 0.0 0.0) (wtemp 0.0+i)) ((= lg ipow)) (set! wpr (cos theta)) (set! wpi (sin theta)) (do ((ii 0 (+ ii 1))) ((= ii prev)) (do ((tempr 0.0+i) (tempi 0.0+i) (i ii (+ i mmax)) (j (+ ii prev) (+ j mmax)) (jj 0 (+ jj 1))) ((>= jj pow)) (set! tempr (- (* wr (hash-table-ref rl (vector-ref rl-syms j))) (* wi (hash-table-ref im (vector-ref im-syms j))))) (set! tempi (+ (* wr (hash-table-ref im (vector-ref im-syms j))) (* wi (hash-table-ref rl (vector-ref rl-syms j))))) (hash-table-set! rl (vector-ref rl-syms j) (- (hash-table-ref rl (vector-ref rl-syms i)) tempr)) (hash-table-set! rl (vector-ref rl-syms i) (+ (hash-table-ref rl (vector-ref rl-syms i)) tempr)) (hash-table-set! im (vector-ref im-syms j) (- (hash-table-ref im (vector-ref im-syms i)) tempi)) (hash-table-set! im (vector-ref im-syms i) (+ (hash-table-ref im (vector-ref im-syms i)) tempi))) (set! wtemp wr) (set! wr (- (* wr wpr) (* wi wpi))) (set! wi (+ (* wi wpr) (* wtemp wpi)))) (set! prev mmax))) rl)) (define (hash-symbol-test) (let ((fvr (hash-table)) (fvi (hash-table)) (rl-syms (make-vector fft-size)) (im-syms (make-vector fft-size)) (scl (/ (* 8.0 pi) fft-size)) (vr 0+i) (vi 0+i) (pk 0+i)) (do ((i 0 (+ i 1)) (x 0.0 (+ x scl))) ((= i fft-size)) (vector-set! rl-syms i (gensym "rl")) (vector-set! im-syms i (gensym "im")) (hash-table-set! fvr (vector-ref rl-syms i) (sin x)) (hash-table-set! fvi (vector-ref im-syms i) 0.0)) (hash-symbol-fft fvr fvi fft-size rl-syms im-syms 1) (do ((mx 0.0) (mxloc 0) (i 0 (+ i 1))) ((= i fft-size) (format () "~A ~A~%" mxloc (/ (* 2 (sqrt mx)) fft-size))) (set! vr (hash-table-ref fvr (vector-ref rl-syms i))) (set! vi (hash-table-ref fvi (vector-ref im-syms i))) (set! pk (+ (* vr vr) (* vi vi))) (when (> pk mx) (set! mx pk) (set! mxloc i))) (hash-table-set! fvr (vector-ref rl-syms 4) 0.0) (hash-table-set! fvi (vector-ref im-syms 4) 0.0) (hash-table-set! fvr (vector-ref rl-syms (- fft-size 4)) 0.0) (hash-table-set! fvi (vector-ref im-syms (- fft-size 4)) 0.0) (do ((mx 0.0) (i 0 (+ i 1))) ((= i fft-size) (format () "noise: ~A~%" mx)) (set! mx (max mx (abs (hash-table-ref fvr (vector-ref rl-syms i))) (abs (hash-table-ref fvi (vector-ref im-syms i)))))) (fill! rl-syms #f) (fill! im-syms #f) (fill! fvr #f) (fill! fvi #f) (hash-symbol-fft #f #f fft-size #f #f 1))) (hash-symbol-test) (clear-and-gc) ;; -------------------------------------------------------------------------------- (format () "~%lets...~%") (define (let-fft rl im n rl-syms im-syms dir) (when rl (do ((i 0 (+ i 1)) (j 0)) ((= i n)) (if (> j i) (let ((tempr (let-ref rl (vector-ref rl-syms j))) (tempi (let-ref im (vector-ref im-syms j)))) (let-set! rl (vector-ref rl-syms j) (let-ref rl (vector-ref rl-syms i))) (let-set! im (vector-ref im-syms j) (let-ref im (vector-ref im-syms i))) (let-set! rl (vector-ref rl-syms i) tempr) (let-set! im (vector-ref im-syms i) tempi))) (let ((m (/ n 2))) (do () ((not (<= 2 m j))) (set! j (- j m)) (set! m (/ m 2))) (set! j (+ j m)))) (let ((ipow (floor (log n 2))) (prev 1)) (do ((mmax 2 (* mmax 2)) (pow (/ n 2) (/ pow 2)) (theta (* pi dir) (* theta 0.5)) (lg 0 (+ lg 1)) (wpr 0.0) (wpi 0.0) (wr 1.0 1.0) (wi 0.0 0.0) (wtemp 0.0+i)) ((= lg ipow)) (set! wpr (cos theta)) (set! wpi (sin theta)) (do ((ii 0 (+ ii 1))) ((= ii prev)) (do ((tempr 0.0+i) (tempi 0.0+i) (i ii (+ i mmax)) (j (+ ii prev) (+ j mmax)) (jj 0 (+ jj 1))) ((>= jj pow)) (set! tempr (- (* wr (let-ref rl (vector-ref rl-syms j))) (* wi (let-ref im (vector-ref im-syms j))))) (set! tempi (+ (* wr (let-ref im (vector-ref im-syms j))) (* wi (let-ref rl (vector-ref rl-syms j))))) (let-set! rl (vector-ref rl-syms j) (- (let-ref rl (vector-ref rl-syms i)) tempr)) (let-set! rl (vector-ref rl-syms i) (+ (let-ref rl (vector-ref rl-syms i)) tempr)) (let-set! im (vector-ref im-syms j) (- (let-ref im (vector-ref im-syms i)) tempi)) (let-set! im (vector-ref im-syms i) (+ (let-ref im (vector-ref im-syms i)) tempi))) (set! wtemp wr) (set! wr (- (* wr wpr) (* wi wpi))) (set! wi (+ (* wi wpr) (* wtemp wpi)))) (set! prev mmax))) rl)) (define (let-test) (let ((fvr (inlet)) (fvi (inlet)) (rl-syms (make-vector fft-size)) (im-syms (make-vector fft-size)) (scl (/ (* 8.0 pi) fft-size)) (vr 0+i) (vi 0+i) (pk 0+i)) (do ((i 0 (+ i 1)) (x 0.0 (+ x scl))) ((= i fft-size)) (vector-set! rl-syms i (gensym "rl-")) (vector-set! im-syms i (gensym "im-")) (varlet fvr (vector-ref rl-syms i) (sin x)) (varlet fvi (vector-ref im-syms i) 0.0)) (let-fft fvr fvi fft-size rl-syms im-syms 1) (do ((mx 0.0) (mxloc 0) (i 0 (+ i 1))) ((= i fft-size) (format () "~A ~A~%" mxloc (/ (* 2 (sqrt mx)) fft-size))) (set! vr (let-ref fvr (vector-ref rl-syms i))) (set! vi (let-ref fvi (vector-ref im-syms i))) (set! pk (+ (* vr vr) (* vi vi))) (when (> pk mx) (set! mx pk) (set! mxloc i))) (let-set! fvr (vector-ref rl-syms 4) 0.0) (let-set! fvi (vector-ref im-syms 4) 0.0) (let-set! fvr (vector-ref rl-syms (- fft-size 4)) 0.0) (let-set! fvi (vector-ref im-syms (- fft-size 4)) 0.0) (do ((mx 0.0) (i 0 (+ i 1))) ((= i fft-size) (format () "noise: ~A~%" mx)) (set! mx (max mx (abs (let-ref fvr (vector-ref rl-syms i))) (abs (let-ref fvi (vector-ref im-syms i)))))) (fill! rl-syms #f) (fill! im-syms #f) (let-fft #f #f fft-size #f #f 1))) (let-test) (clear-and-gc) (format () "~%sundries...~%") ;(set! fft-size (ash 1 15)) (define (string-fft rl im n dir) (when rl (let ((m 0) (j 0) (ipw (floor (log n 2))) (prev 1) (wtemp 0.0+i) (wpr 0.0) (wpi 0.0) (wr 0.0) (wi 0.0) (mmax 2) (pw (/ n 2)) (end 0) (theta (* pi dir))) (do ((tempr 0.0+i) (tempi 0.0+i) (i 0 (+ i 1))) ((= i n)) (when (> j i) (set! tempr (string->number (vector-ref rl j))) (set! tempi (string->number (vector-ref im j))) (vector-set! rl j (vector-ref rl i)) (vector-set! im j (vector-ref im i)) (vector-set! rl i (number->string tempr)) (vector-set! im i (number->string tempi))) (set! m (quotient n 2)) (do () ((not (<= 2 m j))) (set! j (- j m)) (set! m (quotient m 2))) (set! j (+ j m))) (do ((lg 0 (+ lg 1))) ((= lg ipw)) (set! wpr (cos theta)) (set! wpi (sin theta)) (set! wr 1.0) (set! wi 0.0) (do ((ii 0 (+ ii 1))) ((= ii prev)) (set! j (+ ii prev)) (set! end (+ ii (* pw mmax))) (do ((tempr 0.0+i) (tempi 0.0+i) (i ii (+ i mmax))) ((= i end)) (set! tempr (- (* wr (string->number (vector-ref rl j))) (* wi (string->number (vector-ref im j))))) (set! tempi (+ (* wr (string->number (vector-ref im j))) (* wi (string->number (vector-ref rl j))))) (vector-set! rl j (number->string (- (string->number (vector-ref rl i)) tempr))) (vector-set! rl i (number->string (+ (string->number (vector-ref rl i)) tempr))) (vector-set! im j (number->string (- (string->number (vector-ref im i)) tempi))) (vector-set! im i (number->string (+ (string->number (vector-ref im i)) tempi))) (set! j (+ j mmax))) (set! wtemp wr) (set! wr (- (* wr wpr) (* wi wpi))) (set! wi (+ (* wi wpr) (* wtemp wpi)))) (set! prev mmax) (set! mmax (* mmax 2)) (set! pw (quotient pw 2)) (set! theta (* theta 0.5))) rl))) (define (string-checker fvr fvi) (when fvr (let ((vr 0.0+i) (vi 0.0+i) (pk 0.0+i) (mx 0.0) (mxloc 0)) (do ((i 0 (+ i 1))) ((= i fft-size)) (set! vr (string->number (vector-ref fvr i))) (set! vi (string->number (vector-ref fvi i))) (set! pk (+ (* vr vr) (* vi vi))) (when (> pk mx) (set! mx pk) (set! mxloc i))) (format () "~A ~A~%" mxloc (/ (* 2 (sqrt mx)) fft-size)) (format () "~A ~A~%" (string->number (vector-ref fvi 4)) (string->number (vector-ref fvi (- fft-size 4))))) (vector-set! fvr 4 "0.0") (vector-set! fvi 4 "0.0") (vector-set! fvr (- fft-size 4) "0.0") (vector-set! fvi (- fft-size 4) "0.0") (let ((mx 0.0)) (do ((i 0 (+ i 1))) ((= i fft-size)) (set! mx (max mx (abs (string->number (vector-ref fvr i))) (abs (string->number (vector-ref fvi i)))))) (format () "noise: ~A~%" mx)))) (define (string-test) (let ((fvr (make-vector fft-size)) (fvi (make-vector fft-size "0.0")) (scl (/ (* 8.0 pi) fft-size))) (do ((i 0 (+ i 1)) (x 0.0 (+ x scl))) ((= i fft-size)) (vector-set! fvr i (number->string (sin x)))) (string-fft fvr fvi fft-size 1) (string-checker fvr fvi) (string-fft #f #f fft-size 1) (string-checker #f #f) (fill! fvr #f) (fill! fvi #f))) (string-test) (let-temporarily (((*s7* 'float-format-precision) 48)) (format *stderr* "~%radix...~%") (define radix 2) (define (string-radix-fft rl im n dir) (when rl (let ((m 0) (j 0) (ipw (floor (log n 2))) (prev 1) (wtemp 0.0+i) (wpr 0.0) (wpi 0.0) (wr 0.0) (wi 0.0) (mmax 2) (pw (/ n 2)) (end 0) (theta (* pi dir))) (do ((tempr 0.0+i) (tempi 0.0+i) (i 0 (+ i 1))) ((= i n)) (when (> j i) (set! tempr (string->number (vector-ref rl j) radix)) (set! tempi (string->number (vector-ref im j) radix)) (vector-set! rl j (vector-ref rl i)) (vector-set! im j (vector-ref im i)) (vector-set! rl i (number->string tempr radix)) (vector-set! im i (number->string tempi radix))) (set! m (quotient n 2)) (do () ((not (<= 2 m j))) (set! j (- j m)) (set! m (quotient m 2))) (set! j (+ j m))) (do ((lg 0 (+ lg 1))) ((= lg ipw)) (set! wpr (cos theta)) (set! wpi (sin theta)) (set! wr 1.0) (set! wi 0.0) (do ((ii 0 (+ ii 1))) ((= ii prev)) (set! j (+ ii prev)) (set! end (+ ii (* pw mmax))) (do ((tempr 0.0+i) (tempi 0.0+i) (i ii (+ i mmax))) ((= i end)) (set! tempr (- (* wr (string->number (vector-ref rl j) radix)) (* wi (string->number (vector-ref im j) radix)))) (set! tempi (+ (* wr (string->number (vector-ref im j) radix)) (* wi (string->number (vector-ref rl j) radix)))) (vector-set! rl j (number->string (- (string->number (vector-ref rl i) radix) tempr) radix)) (vector-set! rl i (number->string (+ (string->number (vector-ref rl i) radix) tempr) radix)) (vector-set! im j (number->string (- (string->number (vector-ref im i) radix) tempi) radix)) (vector-set! im i (number->string (+ (string->number (vector-ref im i) radix) tempi) radix)) (set! j (+ j mmax))) (set! wtemp wr) (set! wr (- (* wr wpr) (* wi wpi))) (set! wi (+ (* wi wpr) (* wtemp wpi)))) (set! prev mmax) (set! mmax (* mmax 2)) (set! pw (quotient pw 2)) (set! theta (* theta 0.5))) rl))) (define (string-radix-checker fvr fvi) (when fvr (let ((vr 0.0) (vi 0.0) (pk 0.0) (mx 0.0) (mxloc 0)) (do ((i 0 (+ i 1))) ((= i fft-size)) (set! vr (string->number (vector-ref fvr i) radix)) (set! vi (string->number (vector-ref fvi i) radix)) (set! pk (+ (* vr vr) (* vi vi))) (when (> pk mx) (set! mx pk) (set! mxloc i))) (format () "~A ~A~%" mxloc (/ (* 2 (sqrt mx)) fft-size)) (format () "~A ~A~%" (fvi 4) (fvi (- fft-size 4)))) (vector-set! fvr 4 "0.0") (vector-set! fvi 4 "0.0") (vector-set! fvr (- fft-size 4) "0.0") (vector-set! fvi (- fft-size 4) "0.0") (let ((mx 0.0)) (do ((i 0 (+ i 1))) ((= i fft-size)) (set! mx (max mx (abs (string->number (vector-ref fvr i) radix)) (abs (string->number (vector-ref fvi i) radix))))) (format () "noise: ~A~%" mx)))) (define (string-radix-test) (let ((fvr (make-vector fft-size)) (fvi (make-vector fft-size "0.0")) (scl (/ (* 8.0 pi) fft-size))) (do ((i 0 (+ i 1)) (x 0.0 (+ x scl))) ((= i fft-size)) (vector-set! fvr i (number->string (sin x) radix))) (string-radix-fft fvr fvi fft-size 1) (string-radix-checker fvr fvi) (fill! fvr #f) (fill! fvi #f) (string-radix-fft #f #f fft-size 1) (string-radix-checker #f #f))) (string-radix-test)) (clear-and-gc) (define (port-fft rl im n dir) (when rl (let ((m 0) (j 0) (ipw (floor (log n 2))) (prev 1) (wtemp 0.0+i) (wpr 0.0) (wpi 0.0) (wr 0.0) (wi 0.0) (mmax 2) (pw (/ n 2)) (end 0) (theta (* pi dir))) (do ((tempr 0.0+i) (tempi 0.0+i) (i 0 (+ i 1))) ((= i n)) (when (> j i) (set! tempr (with-input-from-string (vector-ref rl j) read)) (set! tempi (with-input-from-string (vector-ref im j) read)) (vector-set! rl j (vector-ref rl i)) (vector-set! im j (vector-ref im i)) (vector-set! rl i (number->string tempr)) (vector-set! im i (number->string tempi))) (set! m (quotient n 2)) (do () ((not (<= 2 m j))) (set! j (- j m)) (set! m (quotient m 2))) (set! j (+ j m))) (do ((lg 0 (+ lg 1))) ((= lg ipw)) (set! wpr (cos theta)) (set! wpi (sin theta)) (set! wr 1.0) (set! wi 0.0) (do ((ii 0 (+ ii 1))) ((= ii prev)) (set! j (+ ii prev)) (set! end (+ ii (* pw mmax))) (do ((tempr 0.0+i) (tempi 0.0+i) (i ii (+ i mmax))) ((= i end)) (set! tempr (- (* wr (with-input-from-string (vector-ref rl j) read)) (* wi (with-input-from-string (vector-ref im j) read)))) (set! tempi (+ (* wr (with-input-from-string (vector-ref im j) read)) (* wi (with-input-from-string (vector-ref rl j) read)))) (vector-set! rl j (number->string (- (with-input-from-string (vector-ref rl i) read) tempr))) (vector-set! rl i (number->string (+ (with-input-from-string (vector-ref rl i) read) tempr))) (vector-set! im j (number->string (- (with-input-from-string (vector-ref im i) read) tempi))) (vector-set! im i (number->string (+ (with-input-from-string (vector-ref im i) read) tempi))) (set! j (+ j mmax))) (set! wtemp wr) (set! wr (- (* wr wpr) (* wi wpi))) (set! wi (+ (* wi wpr) (* wtemp wpi)))) (set! prev mmax) (set! mmax (* mmax 2)) (set! pw (quotient pw 2)) (set! theta (* theta 0.5))) rl))) (define (port-checker fvr fvi) (when fvr (let ((vr 0.0) (vi 0.0) (pk 0.0) (mx 0.0) (mxloc 0)) (do ((i 0 (+ i 1))) ((= i fft-size)) (set! vr (with-input-from-string (vector-ref fvr i) read)) (set! vi (with-input-from-string (vector-ref fvi i) read)) (set! pk (+ (* vr vr) (* vi vi))) (when (> pk mx) (set! mx pk) (set! mxloc i))) (format () "~A ~A~%" mxloc (/ (* 2 (sqrt mx)) fft-size)) (format () "~A ~A~%" (fvi 4) (fvi (- fft-size 4)))) (vector-set! fvr 4 "0.0") (vector-set! fvi 4 "0.0") (vector-set! fvr (- fft-size 4) "0.0") (vector-set! fvi (- fft-size 4) "0.0") (let ((mx 0.0)) (do ((i 0 (+ i 1))) ((= i fft-size)) (set! mx (max mx (abs (with-input-from-string (vector-ref fvr i) read)) (abs (with-input-from-string (vector-ref fvi i) read))))) (format () "noise: ~A~%" mx)))) (define (port-test) (let ((fvr (make-vector fft-size)) (fvi (make-vector fft-size "0.0")) (scl (/ (* 8.0 pi) fft-size))) (do ((i 0 (+ i 1)) (x 0.0 (+ x scl))) ((= i fft-size)) (vector-set! fvr i (number->string (sin x)))) (port-fft fvr fvi fft-size 1) (port-checker fvr fvi) (port-fft #f #f fft-size 1) (port-checker #f #f) (fill! fvr #f) (fill! fvi #f))) (port-test) ;(define-expansion (new-value val) ; `(make-vector (make-list (+ 1 (random 4)) (+ 1 (random 4))) ,val)) (define lists (make-vector 16)) (do ((i 0 (+ i 1)) (j 1) (k 1)) ((= i 16)) (set! (lists i) (make-list k j)) (set! j (+ j 1)) (when (= j 5) (set! j 1) (set! k (+ k 1)))) (define-expansion (new-value val) `(make-vector (vector-ref lists (random 16)) ,val)) (define-expansion (old-value ref) `((subvector ,ref 0 1) 0)) (define (complex-2d-fft data n dir) (when data (do ((i 0 (+ i 1)) (j 0)) ((= i n)) (if (> j i) (let ((temp (vector-ref data j))) (vector-set! data j (vector-ref data i)) (vector-set! data i temp))) (let ((m (/ n 2))) (do () ((not (<= 2 m j))) (set! j (- j m)) (set! m (/ m 2))) (set! j (+ j m)))) (let ((ipow (floor (log n 2))) (prev 1)) (do ((lg 0 (+ lg 1)) (mmax 2 (* mmax 2)) (pow (/ n 2) (/ pow 2)) (theta (complex 0.0 (* pi dir)) (* theta 0.5)) (wpc 0.0) (wc 1.0 1.0)) ((= lg ipow)) (set! wpc (exp theta)) (do ((ii 0 (+ ii 1))) ((= ii prev)) (do ((tc 0.0+i) (jj 0 (+ jj 1)) (i ii (+ i mmax)) (j (+ ii prev) (+ j mmax))) ((>= jj pow)) (set! tc (* wc (old-value (vector-ref data j)))) (vector-set! data j (new-value (- (old-value (vector-ref data i)) tc))) (vector-set! data i (new-value (+ (old-value (vector-ref data i)) tc)))) (set! wc (* wc wpc))) (set! prev mmax))) data)) (define (complex-2d-checker fvr) (when fvr (let ((pk 0.0) (mx 0.0) (mxloc 0)) (do ((i 0 (+ i 1))) ((= i fft-size)) (set! pk (magnitude (old-value (vector-ref fvr i)))) (when (> pk mx) (set! mx pk) (set! mxloc i))) (format () "~A ~A~%" mxloc (/ (* 2 (sqrt mx)) fft-size)) (format () "~A ~A~%" (magnitude (old-value (fvr 4))) (magnitude (old-value (fvr (- fft-size 4)))))) (vector-set! fvr 4 (new-value 0.0)) (vector-set! fvr (- fft-size 4) (new-value 0.0)) (let ((mx 0.0)) (do ((i 0 (+ i 1))) ((= i fft-size)) (set! mx (max mx (magnitude (old-value (vector-ref fvr i)))))) (format () "noise: ~A~%" mx)))) (define (complex-2d-test) (let ((fvr (make-vector fft-size)) (scl (/ (* 8.0 pi) fft-size))) (do ((i 0 (+ i 1)) (x 0.0 (+ x scl))) ((= i fft-size)) (vector-set! fvr i (new-value (sin x)))) (complex-2d-fft fvr fft-size 1) (complex-2d-checker fvr) (complex-2d-fft #f fft-size 1) (complex-2d-checker #f) (fill! fvr #f))) (complex-2d-test) (clear-and-gc) ;(set! fft-size (ash 1 17)) (define (complex-hash-fft data n dir) (when data (do ((i 0 (+ i 1)) (j 0)) ((= i n)) (if (> j i) (let ((temp (vector-ref data j))) (vector-set! data j (vector-ref data i)) (vector-set! data i temp))) (let ((m (/ n 2))) (do () ((not (<= 2 m j))) (set! j (- j m)) (set! m (/ m 2))) (set! j (+ j m)))) (let ((ipow (floor (log n 2))) (prev 1)) (do ((mmax 2 (* mmax 2)) (pow (/ n 2) (/ pow 2)) (theta (complex 0.0 (* pi dir)) (* theta 0.5)) (lg 0 (+ lg 1))) ((= lg ipow)) (let ((wpc (exp theta)) (wc 1.0)) (do ((ii 0 (+ ii 1))) ((= ii prev)) (do ((i ii (+ i mmax)) (j (+ ii prev) (+ j mmax)) (jj 0 (+ jj 1))) ((>= jj pow)) (let ((tc (* wc (hash-table-ref (vector-ref data j) 'z)))) ; this could be (data j 'z)! (vector-set! data j (hash-table 'z (- (hash-table-ref (vector-ref data i) 'z) tc))) (vector-set! data i (hash-table 'z (+ (hash-table-ref (vector-ref data i) 'z) tc))))) (set! wc (* wc wpc))) (set! prev mmax)))) data)) (define (complex-hash-checker fvr) (when fvr (let ((pk 0.0) (mx 0.0) (mxloc 0)) (do ((i 0 (+ i 1))) ((= i fft-size)) (set! pk (magnitude (hash-table-ref (vector-ref fvr i) 'z))) (when (> pk mx) (set! mx pk) (set! mxloc i))) (format () "~A ~A~%" mxloc (/ (* 2 (sqrt mx)) fft-size)) (format () "~A ~A~%" (magnitude (hash-table-ref (vector-ref fvr 4) 'z)) (magnitude (hash-table-ref (vector-ref fvr (- fft-size 4)) 'z)))) (vector-set! fvr 4 (hash-table 'z 0.0)) (vector-set! fvr (- fft-size 4) (hash-table 'z 0.0)) (let ((mx 0.0)) (do ((i 0 (+ i 1))) ((= i fft-size)) (set! mx (max mx (magnitude (hash-table-ref (vector-ref fvr i) 'z))))) (format () "noise: ~A~%" mx)))) (define (complex-hash-test) (let ((fvr (make-vector fft-size)) (scl (/ (* 8.0 pi) fft-size))) (do ((i 0 (+ i 1)) (x 0.0 (+ x scl))) ((= i fft-size)) (vector-set! fvr i (hash-table 'z (sin x)))) (complex-hash-fft fvr fft-size 1) (complex-hash-checker fvr) (complex-hash-fft #f fft-size 1) (complex-hash-checker #f) (fill! fvr #f))) (complex-hash-test) (clear-and-gc) (define (iterator-vector-fft rl im n dir) (when rl (let ((m 0) (j 0) (ipw (floor (log n 2))) (prev 1) (wtemp 0.0+i) (wpr 0.0) (wpi 0.0) (wr 0.0) (wi 0.0) (mmax 2) (pw (/ n 2)) (end 0) (theta (* pi dir)) (tempbr #f) (tempbi #f)) (do ((tempr 0.0+i) (tempi 0.0+i) (i 0 (+ i 1))) ((= i n)) (when (> j i) (set! tempbr (vector-ref rl j)) (set! tempbi (vector-ref im j)) (vector-set! rl j (vector-ref rl i)) (vector-set! im j (vector-ref im i)) (vector-set! rl i tempbr) (vector-set! im i tempbi)) (set! m (quotient n 2)) (do () ((not (<= 2 m j))) (set! j (- j m)) (set! m (quotient m 2))) (set! j (+ j m))) (do ((lg 0 (+ lg 1))) ((= lg ipw)) (set! wpr (cos theta)) (set! wpi (sin theta)) (set! wr 1.0) (set! wi 0.0) (do ((ii 0 (+ ii 1))) ((= ii prev)) (set! j (+ ii prev)) (set! end (+ ii (* pw mmax))) (do ((tempr 0.0+i) (tempi 0.0+i) (i ii (+ i mmax))) ((= i end)) (set! tempr (- (* wr (iterate (vector-ref rl j))) (* wi (iterate (vector-ref im j))))) (set! tempi (+ (* wr (iterate (vector-ref im j))) (* wi (iterate (vector-ref rl j))))) (vector-set! rl j (make-iterator (make-float-vector 2 (- (iterate (vector-ref rl i)) tempr)))) (vector-set! rl i (make-iterator (make-float-vector 2 (+ (iterate (vector-ref rl i)) tempr)))) (vector-set! im j (make-iterator (make-float-vector 2 (- (iterate (vector-ref im i)) tempi)))) (vector-set! im i (make-iterator (make-float-vector 2 (+ (iterate (vector-ref im i)) tempi)))) (set! j (+ j mmax))) (set! wtemp wr) (set! wr (- (* wr wpr) (* wi wpi))) (set! wi (+ (* wi wpr) (* wtemp wpi)))) (set! prev mmax) (set! mmax (* mmax 2)) (set! pw (quotient pw 2)) (set! theta (* theta 0.5))) rl))) (define (iterator-vector-checker fvr fvi) (when fvr (let ((vr 0.0+i) (vi 0.0+i) (pk 0.0+i) (mx 0.0) (mxloc 0)) (do ((i 0 (+ i 1))) ((= i fft-size)) (set! vr (iterate (vector-ref fvr i))) (set! vi (iterate (vector-ref fvi i))) (set! pk (+ (* vr vr) (* vi vi))) (when (> pk mx) (set! mx pk) (set! mxloc i))) (format () "~A ~A~%" mxloc (/ (* 2 (sqrt mx)) fft-size)) (format () "~A ~A~%" (iterate (vector-ref fvi 4)) (iterate (vector-ref fvi (- fft-size 4))))) (vector-set! fvr 4 (make-iterator (make-float-vector 2 0.0))) (vector-set! fvi 4 (make-iterator (make-float-vector 2 0.0))) (vector-set! fvr (- fft-size 4) (make-iterator (make-float-vector 2 0.0))) (vector-set! fvi (- fft-size 4) (make-iterator (make-float-vector 2 0.0))) (let ((mx 0.0)) (do ((i 0 (+ i 1))) ((= i fft-size)) (set! mx (max mx (abs (iterate (vector-ref fvr i))) (abs (iterate (vector-ref fvi i)))))) (format () "noise: ~A~%" mx)))) (define (iterator-vector-test) (let ((fvr (make-vector fft-size)) (fvi (make-vector fft-size)) (scl (/ (* 8.0 pi) fft-size))) (do ((i 0 (+ i 1)) (x 0.0 (+ x scl))) ((= i fft-size)) (vector-set! fvi i (make-iterator (make-float-vector 2 0.0))) (vector-set! fvr i (make-iterator (make-float-vector 2 (sin x))))) (iterator-vector-fft fvr fvi fft-size 1) (iterator-vector-checker fvr fvi) (iterator-vector-fft #f #f fft-size 1) (iterator-vector-checker #f #f) (fill! fvr #f) (fill! fvi #f))) (iterator-vector-test) (clear-and-gc) (define (complex-closure-fft data n dir) (when data (do ((i 0 (+ i 1)) (j 0)) ((= i n)) (if (> j i) (let ((temp (vector-ref data j))) (vector-set! data j (vector-ref data i)) (vector-set! data i temp))) (let ((m (/ n 2))) (do () ((not (<= 2 m j))) (set! j (- j m)) (set! m (/ m 2))) (set! j (+ j m)))) (let ((ipow (floor (log n 2))) (prev 1)) (do ((mmax 2 (* mmax 2)) (pow (/ n 2) (/ pow 2)) (theta (complex 0.0 (* pi dir)) (* theta 0.5)) (lg 0 (+ lg 1))) ((= lg ipow)) (let ((wpc (exp theta)) (wc 1.0)) (do ((ii 0 (+ ii 1))) ((= ii prev)) (do ((i ii (+ i mmax)) (j (+ ii prev) (+ j mmax)) (jj 0 (+ jj 1))) ((>= jj pow)) (let ((tc (* wc ((vector-ref data j))))) (vector-set! data j (let ((z (- ((vector-ref data i)) tc))) (lambda () z))) (vector-set! data i (let ((z (+ ((vector-ref data i)) tc))) (lambda () z))))) (set! wc (* wc wpc))) (set! prev mmax)))) data)) (define (complex-closure-checker fvr) (when fvr (let ((pk 0.0) (mx 0.0) (mxloc 0)) (do ((i 0 (+ i 1))) ((= i fft-size)) (set! pk (magnitude ((vector-ref fvr i)))) (when (> pk mx) (set! mx pk) (set! mxloc i))) (format () "~A ~A~%" mxloc (/ (* 2 (sqrt mx)) fft-size)) (format () "~A ~A~%" (magnitude ((vector-ref fvr 4))) (magnitude ((vector-ref fvr (- fft-size 4)))))) (vector-set! fvr 4 (let ((z 0.0)) (lambda () z))) (vector-set! fvr (- fft-size 4) (let ((z 0.0)) (lambda () z))) (let ((mx 0.0)) (do ((i 0 (+ i 1))) ((= i fft-size)) (set! mx (max mx (magnitude ((vector-ref fvr i)))))) (format () "noise: ~A~%" mx)))) (define (complex-closure-test) (let ((fvr (make-vector fft-size)) (scl (/ (* 8.0 pi) fft-size))) (do ((i 0 (+ i 1)) (x 0.0 (+ x scl))) ((= i fft-size)) (vector-set! fvr i (let ((z (sin x))) (lambda () z)))) (complex-closure-fft fvr fft-size 1) (complex-closure-checker fvr) (complex-closure-fft #f fft-size 1) (complex-closure-checker #f) (fill! fvr #f))) (complex-closure-test) (clear-and-gc) ;; -------------------------------------------------------------------------------- (load "s7test-block.so" (sublet (curlet) (cons 'init_func 'block_init))) (format () "~%blocks...~%") (when (> total-memory (* 10 big-size)) (clear-and-gc) (let ((bigv (make-block big-size))) (test (length bigv) big-size) (set! (bigv (- big-size 10000000)) 1.0) (test (bigv (- big-size 10000000)) 1.0) (reverse! bigv) (test (bigv (- 10000000 1)) 1.0) (fill! bigv 0.0); (- big-size 9000000)) (test (bigv (- big-size 10000000)) 0.0) )) (clear-and-gc) (when (> total-memory (* (ash 1 33) 8)) (let ((b (make-block (ash 1 33)))) (set! (b (+ (ash 1 32) 10)) 1.0) (test (b (+ (ash 1 32) 10)) 1.0)) (clear-and-gc)) (define (block-fft rl im n dir) (when rl (do ((i 0 (+ i 1)) (j 0)) ((= i n)) (if (> j i) (let ((tempr (block-ref rl j)) (tempi (block-ref im j))) (block-set! rl j (block-ref rl i)) (block-set! im j (block-ref im i)) (block-set! rl i tempr) (block-set! im i tempi))) (let ((m (/ n 2))) (do () ((not (<= 2 m j))) (set! j (- j m)) (set! m (/ m 2))) (set! j (+ j m)))) (let ((ipow (floor (log n 2))) (prev 1) (tempr 0.0) (tempi 0.0) (wtemp 0.0)) (do ((mmax 2 (* mmax 2)) (pow (/ n 2) (/ pow 2)) (theta (* pi dir) (* theta 0.5)) (lg 0 (+ lg 1))) ((= lg ipow)) (let ((wpr (cos theta)) (wpi (sin theta)) (wr 1.0) (wi 0.0)) (do ((ii 0 (+ ii 1))) ((= ii prev)) (do ((i ii (+ i mmax)) (j (+ ii prev) (+ j mmax)) (jj 0 (+ jj 1))) ((>= jj pow)) (set! tempr (- (* wr (block-ref rl j)) (* wi (block-ref im j)))) (set! tempi (+ (* wr (block-ref im j)) (* wi (block-ref rl j)))) (block-set! rl j (- (block-ref rl i) tempr)) (block-set! rl i (+ (block-ref rl i) tempr)) (block-set! im j (- (block-ref im i) tempi)) (block-set! im i (+ (block-ref im i) tempi))) (set! wtemp wr) (set! wr (- (* wr wpr) (* wi wpi))) (set! wi (+ (* wi wpr) (* wtemp wpi)))) (set! prev mmax)))) rl)) (define (block-test) (let ((fvr (make-block fft-size)) (fvi (make-block fft-size)) (scl (/ (* 8.0 pi) fft-size))) (do ((i 0 (+ i 1)) (x 0.0 (+ x scl))) ((= i fft-size)) (block-set! fvr i (sin x))) (block-fft fvr fvi fft-size 1) (do ((mx 0.0) (mxloc 0) (vr 0.0) (vi 0.0) (pk 0.0) (i 0 (+ i 1))) ((= i fft-size) (format () "~A ~A~%" mxloc (/ (* 2 (sqrt mx)) fft-size))) (set! vr (block-ref fvr i)) (set! vi (block-ref fvi i)) (set! pk (+ (* vr vr) (* vi vi))) (when (> pk mx) (set! mx pk) (set! mxloc i))) (set! (fvr 4) 0.0) (set! (fvi 4) 0.0) (set! (fvr (- fft-size 4)) 0.0) (set! (fvi (- fft-size 4)) 0.0) (do ((mx 0.0) (i 0 (+ i 1))) ((= i fft-size) (format () "noise: ~A~%" mx)) (set! mx (max mx (abs (block-ref fvr i)) (abs (block-ref fvi i))))) (block-fft #f #f fft-size 1))) (block-test) (define (block-vector-fft rl im n dir) (when rl (let ((tempr 0.0) (tempi 0.0) (m 0) (j 0) (ipw (floor (log n 2))) (prev 1) (wtemp 0.0) (wpr 0.0) (wpi 0.0) (wr 0.0) (wi 0.0) (mmax 2) (pw (/ n 2)) (end 0) (theta (* pi dir)) (tempbr #f) (tempbi #f)) (do ((i 0 (+ i 1))) ((= i n)) (when (> j i) (set! tempbr (vector-ref rl j)) (set! tempbi (vector-ref im j)) (vector-set! rl j (vector-ref rl i)) (vector-set! im j (vector-ref im i)) (vector-set! rl i tempbr) (vector-set! im i tempbi)) (set! m (quotient n 2)) (do () ((not (<= 2 m j))) (set! j (- j m)) (set! m (quotient m 2))) (set! j (+ j m))) (do ((lg 0 (+ lg 1))) ((= lg ipw)) (set! wpr (cos theta)) (set! wpi (sin theta)) (set! wr 1.0) (set! wi 0.0) (do ((ii 0 (+ ii 1))) ((= ii prev)) (set! j (+ ii prev)) (set! end (+ ii (* pw mmax))) (do ((i ii (+ i mmax))) ((= i end)) (set! tempr (- (* wr (block-ref (vector-ref rl j) 0)) (* wi (block-ref (vector-ref im j) 0)))) (set! tempi (+ (* wr (block-ref (vector-ref im j) 0)) (* wi (block-ref (vector-ref rl j) 0)))) (vector-set! rl j (block (- (block-ref (vector-ref rl i) 0) tempr))) (vector-set! rl i (block (+ (block-ref (vector-ref rl i) 0) tempr))) (vector-set! im j (block (- (block-ref (vector-ref im i) 0) tempi))) (vector-set! im i (block (+ (block-ref (vector-ref im i) 0) tempi))) (set! j (+ j mmax))) (set! wtemp wr) (set! wr (- (* wr wpr) (* wi wpi))) (set! wi (+ (* wi wpr) (* wtemp wpi)))) (set! prev mmax) (set! mmax (* mmax 2)) (set! pw (quotient pw 2)) (set! theta (* theta 0.5))) rl))) (define (block-vector-checker fvr fvi) (when fvr (let ((vr 0.0) (vi 0.0) (pk 0.0) (mx 0.0) (mxloc 0)) (do ((i 0 (+ i 1))) ((= i fft-size)) (set! vr (block-ref (vector-ref fvr i) 0)) (set! vi (block-ref (vector-ref fvi i) 0)) (set! pk (+ (* vr vr) (* vi vi))) (when (> pk mx) (set! mx pk) (set! mxloc i))) (format () "~A ~A~%" mxloc (/ (* 2 (sqrt mx)) fft-size)) (format () "~A ~A~%" (fvi 4) (fvi (- fft-size 4)))) (vector-set! fvr 4 (block 0.0)) (vector-set! fvi 4 (block 0.0)) (vector-set! fvr (- fft-size 4) (block 0.0)) (vector-set! fvi (- fft-size 4) (block 0.0)) (let ((mx 0.0)) (do ((i 0 (+ i 1))) ((= i fft-size)) (set! mx (max mx (abs (block-ref (vector-ref fvr i) 0)) (abs (block-ref (vector-ref fvi i) 0))))) (format () "noise: ~A~%" mx)))) (define (block-vector-test) (let ((fvr (make-vector fft-size)) (fvi (make-vector fft-size (block 0.0))) (scl (/ (* 8.0 pi) fft-size))) (do ((i 0 (+ i 1)) (x 0.0 (+ x scl))) ((= i fft-size)) (vector-set! fvr i (block (sin x)))) (block-vector-fft fvr fvi fft-size 1) (block-vector-checker fvr fvi) (block-vector-fft #f #f fft-size 1) (block-vector-checker #f #f))) (block-vector-test) ;(clear-and-gc) ;(clear-and-gc) ;(*s7* 'memory-usage) (when (> (*s7* 'profile) 0) (show-profile 200)) (exit)