;;;-*- Mode:LISP; Package:BENCH-FFTU; Base:8 -*-
;;; From the "Dick Gabriel" Benchmark Series.
;;; Enhancements (C) Copyright 1983, Lisp Machine, Inc.

(EVAL-WHEN (EVAL COMPILE LOAD)
  (DEFCONST *TO-UCOMPILE* '(FFT))
  (MAPC #'(LAMBDA (X)
	    (PUTPROP X T 'COMPILER:MICROCOMPILE)
	    (PUTPROP X
		     T ; ':DYNAMIC
		     ':DEPEND-ON-BEING-MICROCOMPILED))
	    *TO-UCOMPILE*))

;;;BEGIN
;;;FFT
;Barrow FFT
;Here is the Barrow FFT benchmark which tests floating operations
;of various types, including flonum arrays. (ARRAYCALL FLONUM A I)
;accesses the I'th element of the FLONUM array A, where these arrays are
;0-based. (STORE (ARRAYCALL FLONUM A I) V) stores the value V in the
;I'th element of the FLONUM array A. 

;There was a fair amount of FLONUM GC's in the SAIL MacLisp run, which,
;when it needed to CORE up during GC, took 4.5 seconds of CPU time for the
;computation and 15 seconds for GC. Other configurations of memory required
;only 1.5 seconds for GC.

;Refer to this as FFT.
;			-rpg-


(DEFMACRO AREF$ (A &REST L)
  #+MACLISP `(ARRAYCALL FLONUM ,A ,@L)
  #-MACLISP `(AREF ,A ,@L))

;;; *-*lisp*-* 
;;; From Rich Duda, by way of Harry Barrow -- 3/26/82 

(DEFUN FFT					;Fast Fourier Transform
  (AREAL AIMAG smallp)		;AREAL = real part
  (PROG						;AIMAG = imaginary part
   (AR AI PIE I J K M N LE LE1 IP NV2 NM1 UR UI WR WI TR TI)
   (declare (fixnum I J K M N LE LE1 IP NV2 NM1))
    (SETQ AR (PROGN AREAL))			;Initialize
    (SETQ AI (PROGN AIMAG))
    (SETQ PIE (if smallp 3.14159265s00 3.141592653589793))
    (SETQ N (CADR (ARRAYDIMS AR)))
    (SETQ N (1- N))
    (SETQ NV2 (// N 2))
    (SETQ NM1 (1- N))
    (SETQ M 0)						;Compute M = log(N)
    (SETQ I 1)
   L1 (COND
       ((< I N)(SETQ M (1+ M))(SETQ I (+ I I))(GO L1)))
    (COND ((NOT (EQUAL N (^ 2 M)))
	   (ERROR "array size not a power of two.")))
    (SETQ J 1)						;Interchange elements
    (SETQ I 1)						;in bit-reversed order
   L3 (COND ((the fixnum (< I J))
	     (SETQ TR (AREF$ AR J))
	     (SETQ TI (AREF$ AI J))
	     (SETF (AREF$ AR J)
		    (AREF$ AR I))
	     (SETF (AREF$ AI J)
		    (AREF$ AI I))
	     (SETF (AREF$ AR I) TR)	
	     (SETF (AREF$ AI I) TI)))
    (SETQ K NV2)
   L6 (COND ((the fixnum (< K J))
	     (SETQ J (the fixnum (- J K)))
	     (SETQ K (// K 2))
	     (GO L6)))
    (SETQ J (the fixnum (+ J K)))
    (SETQ I (the fixnum (1+ I)))
    (COND ((the fixnum (< I N))(GO L3)))
    (DO L 1 (the fixnum (1+ L)) (the fixnum (> L M))				;Loop thru stages
	(declare (fixnum l))
	(SETQ LE (^ 2 L))
	(SETQ LE1 (// LE 2))
	(SETQ UR (if smallp 1.0s00 1.0))
	(SETQ UI (if smallp 0.0s00 0.0))
	(SETQ WR (COS (//$ PIE (if smallp (small-float le1) (FLOAT LE1)))))
	(SETQ WI (SIN (//$ PIE (if smallp (small-float le1) (FLOAT LE1)))))
	(DO J 1 (the fixnum (1+ J)) (the fixnum (> J LE1))			;Loop thru butterflies
	    (DO I J (the fixnum (+ I LE)) (> I N)			;Do a butterfly
		(SETQ IP (+ I LE1))
		(SETQ TR (-$ (*$ (AREF$ AR IP) UR)
			     (*$ (AREF$ AI IP) UI)))
		(SETQ TI (+$ (*$ (AREF$ AR IP) UI)
			     (*$ (AREF$ AI IP) UR)))
		(SETF (AREF$ AR IP)
		       (-$ (AREF$ AR I) TR))
		(SETF (AREF$ AI IP)
		       (-$ (AREF$ AI I) TI))
		(SETF (AREF$ AR I)
		       (+$ (AREF$ AR I) TR))
		(SETF (AREF$ AI I)
		       (+$ (AREF$ AI I) TI)))
	    (SETQ TR (-$ (*$ UR WR) (*$ UI WI)))
	    (SETQ TI (+$ (*$ UR WI) (*$ UI WR)))
	    (SETQ UR TR)
	    (SETQ UI TI)))
    (RETURN T)))



;;; Sets up the two arrays
(DECLARE (SPECIAL *RE* *IM*))

(SETQ *RE* (*ARRAY NIL 'FLONUM 1025.))

(SETQ *IM* (*ARRAY NIL 'FLONUM 1025.))

;;; The timer which does 10 calls on FFT

;(include "timer.lsp")

(defconst *ucode-loaded? nil)

(defun load-ucode ()
  (apply #'compiler:ma-load *to-ucompile*)
  (setq *ucode-loaded? t))

(timer timit 
       (if *ucode-loaded?
	   (do ((ntimes 0 (1+ ntimes)))
	       ((= ntimes 10.))
	     (fft *re* *im* nil))
	 "ucode not loaded"))

(defun test (n)
  (dotimes (i n) (fft *re$* *im$* t)))

(DEFVAR *RE$* (make-array 1025. ':initial-value 0.0s00))

(DEFVAR *IM$* (make-array 1025. ':initial-value 0.0s00))

(timer timit$
       (if *ucode-loaded?
	   (do ((ntimes 0 (1+ ntimes)))
	       ((= ntimes 10.))
	     (fft *re$* *im$* t))
	 "ucode not loaded"))

;;;END