[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Compiler/GC interaction?
This resulted when I tried to test a friend's polynomial package on the
latest release of CMU Common Lisp for SPARCstations. The problem shows
up at different points in the file, depending on how much compilation
has been done prior to compiling the file. In every case the function
C::VALID-FUNCTION-USE complains about a type error. I suspect it is an
interaction between some foreign function stuff in the compiler and the
GC.
The files involved follow the script. Sorry about the length.
-- Chuck Fry Chucko@charon.arc.nasa.gov
Script started on Fri Jan 24 14:24:18 1992
mccarthy 33% cmucl
CMU Common Lisp 15c, running on mccarthy.arc.nasa.gov
Hemlock 3.5 (15c), Python 1.0(15c), target SPARCstation/Sun 4
Send bug reports and questions to cmucl-bugs@cs.cmu.edu.
* (load "package.lisp")
T
* (compile-file "poly.lisp")
Python version 1.0(15c), VM version SPARCstation/Sun 4 on 24 JAN 92 02:24:54 pm.
Compiling: /net/home/kronos/id3/chucko/cl/speakers/poly.lisp 27 NOV 91 03:13:53 pm
File: /net/home/kronos/id3/chucko/cl/speakers/poly.lisp
In: DEFUN POLY-+
(LOOP FOR (P1 . REST1) ON (IF # # POLY1)
FOR (P2 . REST2) ON (IF # # POLY2)
COLLECT (+ P1 P2) INTO RESULT
UNLESS REST1 RETURN #
UNLESS REST2 RETURN #
FINALLY (RETURN RESULT))
Error: (during macroexpansion)
Error in function LOOP::PARSE-AND-CLAUSES.
Invalid clause for inside a conditional: RETURN
Converted POLY-+.
Compiling DEFUN POLY-+:
File: /net/home/kronos/id3/chucko/cl/speakers/poly.lisp
In: DEFUN POLY-+
(LABELS ((POLY-+-INTERNAL # #)) (COND (# #) (# #) (T #)))
Warning: Variable POLY1 defined but never used.
Warning: Variable POLY2 defined but never used.
Converted POLY-*.
Compiling DEFUN POLY-*:
Converted POLY-/.
Compiling DEFUN POLY-/:
Converted POLY-DERIV.
Compiling DEFUN POLY-DERIV:
Converted POWER-POLY.
Compiling DEFUN POWER-POLY:
Converted POLY-NORMALIZE.
Compiling DEFUN POLY-NORMALIZE:
Converted POLY-W-NORMALIZE.
Compiling DEFUN POLY-W-NORMALIZE:
Converted POLY-W.
Compiling DEFUN POLY-W:
Compiling Top-Level Form:
Converted W-SCALE-POLY.
Compiling DEFUN W-SCALE-POLY:
Converted 3DB-NORMALIZE.
Compiling DEFUN 3DB-NORMALIZE:
Converted FIND-CUTOFF-FREQUENCY.
Compiling DEFUN FIND-CUTOFF-FREQUENCY:
Converted W-SCALE-POLES.
Compiling DEFUN W-SCALE-POLES:
Converted POLY-ZEROIZE.
Compiling DEFUN POLY-ZEROIZE:
Converted TF-*.
Compiling DEFUN TF-*:
Converted TF-VALUE.
Compiling DEFUN TF-VALUE:
Converted F-AND-Q-FROM-POLES.
Compiling DEFUN F-AND-Q-FROM-POLES:
Compiling Top-Level Form:
Converted F-AND-Q-FROM-POLE.
Compiling DEFUN F-AND-Q-FROM-POLE:
Converted POLY-FROM-POLE-LOCATION.
Compiling DEFUN POLY-FROM-POLE-LOCATION:
Converted POLE-FROM-W-AND-Q.
Compiling DEFUN POLE-FROM-W-AND-Q:
Converted SORT-POLES.
Compiling DEFUN SORT-POLES:
Converted SUBSUME-SINGLE-POLES.
Compiling DEFUN SUBSUME-SINGLE-POLES:
Converted POLY-VALUE.
Compiling DEFUN POLY-VALUE:
Converted BUTTERWORTH-POLES.
Compiling DEFUN BUTTERWORTH-POLES:
Converted ELLIPTIC-POLES.
Compiling DEFUN ELLIPTIC-POLES:
[GC threshold exceeded with 2,131,384 bytes in use. Commencing GC.]
[GC completed with 490,504 bytes retained and 1,640,880 bytes freed.]
[GC will next occur when at least 2,490,504 bytes are in use.]
Converted CHEBYSHEV-POLES.
Compiling DEFUN CHEBYSHEV-POLES:
Error: Read error at 8667:
"(if (< #-Symbolics -1.0e-154/\ "
Arithmetic error FLOATING-POINT-UNDERFLOW signalled.
Operation was SCALE-FLOAT, operands (0.6703904 -511).
Compiling Top-Level Form:
Converted CHEBYSHEV-K-TO-E.
Compiling DEFUN CHEBYSHEV-K-TO-E:
Converted CHEBYSHEV-RIPPLE.
Compiling DEFUN CHEBYSHEV-RIPPLE:
Converted CHEBYSHEV-RIPPLE-TO-E.
Compiling DEFUN CHEBYSHEV-RIPPLE-TO-E:
Converted QUASI-BUTTERWORTH-POLES.
Compiling DEFUN QUASI-BUTTERWORTH-POLES:
Converted BESSEL-POLES.
Compiling DEFUN BESSEL-POLES:
Converted BESSEL-POLY.
Compiling DEFUN BESSEL-POLY:
Converted ELLIPTIC-CONTOUR-POLES.
Compiling DEFUN ELLIPTIC-CONTOUR-POLES:
Converted GAUSSIAN-POLES.
Compiling DEFUN GAUSSIAN-POLES:
Type-error in C::VALID-FUNCTION-USE:
-385613820 is not of type LIST
Restarts:
0: Return to Top-Level.
Debug (type H for help)
(DEBUG::DEBUG-LOOP)
0] backtrace
(DEBUG::DEBUG-LOOP)
(DEBUG:INTERNAL-DEBUG)
(INVOKE-DEBUGGER #<TYPE-ERROR.B115A05>)
(ERROR TYPE-ERROR :FUNCTION-NAME C::VALID-FUNCTION-USE :DATUM ...)
("DEFERR OBJECT-NOT-LIST-ERROR" C::VALID-FUNCTION-USE
#<System-Area pointer: #xE00958>
#<Alien value, Address = #xF7FFE420, Size = 2368, Type = MACH:SIGCONTEXT>
(556))
("Foreign function call land")
(C::VALID-FUNCTION-USE #<C::COMBINATION #xB0E901D
FUN= #<C::REF #xB0E91FD>
ARGS= (# #)>
#<FUNCTION-TYPE (FUNCTION (# #) (MEMBER T NIL))>
:ARGUMENT-TEST
NIL
...)
(C::PROBABLE-TYPE-CHECK-P #<Continuation c1>)
(C::GENERATE-TYPE-CHECKS #<C:COMPONENT #xB1117E5 NAME= "DEFUN GAUSSIAN-POLES">)
(C::IR1-PHASES #<C:COMPONENT #xB1117E5 NAME= "DEFUN GAUSSIAN-POLES">)
(C::COMPILE-COMPONENT #<C:COMPONENT #xB1117E5 NAME= "DEFUN GAUSSIAN-POLES">
#<Fasl-File "poly.sparcf">)
(C::COMPILE-TOP-LEVEL (#<LAMBDA #xB107865>) #<Fasl-File "poly.sparcf">)
(C::CONVERT-AND-MAYBE-COMPILE (C::%DEFUN 'GAUSSIAN-POLES
#'(LAMBDA # #)
NIL
'(DEFUN GAUSSIAN-POLES # #))
((C::%DEFUN 'GAUSSIAN-POLES #'# NIL '#)
C::ORIGINAL-SOURCE-START 0 38)
#<Fasl-File "poly.sparcf">)
(C::PROCESS-FORM (C::%DEFUN 'GAUSSIAN-POLES
#'(LAMBDA # #)
NIL
'(DEFUN GAUSSIAN-POLES # #))
(C::ORIGINAL-SOURCE-START 0 38)
#<Fasl-File "poly.sparcf">)
(C::PROCESS-FORM (DEFUN GAUSSIAN-POLES (N) (3DB-NORMALIZE #))
(C::ORIGINAL-SOURCE-START 0 38)
#<Fasl-File "poly.sparcf">)
(C::PROCESS-FORM #<unavailable-arg>
(DEFUN GAUSSIAN-POLES (N) (3DB-NORMALIZE #))
(C::ORIGINAL-SOURCE-START 0 38)
#<Fasl-File "poly.sparcf">)[:EXTERNAL]
(#:G4)
(C::SUB-COMPILE-FILE #<Source-Info> #<Fasl-File "poly.sparcf"> NIL)
(C::SUB-COMPILE-FILE 2
#<Source-Info>
#<Fasl-File "poly.sparcf">
#<Function C::SUB-COMPILE-FILE {12907D9}>)[:EXTERNAL]
(COMPILE-FILE "poly.lisp" :OUTPUT-FILE T :ERROR-FILE ...)
(EXTENSIONS:INTERACTIVE-EVAL (COMPILE-FILE "poly.lisp"))
(LISP::%TOP-LEVEL)
(LISP::%INITIAL-FUNCTION)
0] f 6
(C::VALID-FUNCTION-USE #<C::COMBINATION #xB0E901D
FUN= #<C::REF #xB0E91FD>
ARGS= (# #)>
#<FUNCTION-TYPE (FUNCTION (# #) (MEMBER T NIL))>
:ARGUMENT-TEST
NIL
...)
6] pp
(C::VALID-FUNCTION-USE #<C::COMBINATION #xB0E901D
FUN= #<C::REF #xB0E91FD
LEAF= #<C::GLOBAL-VAR #xB109C3D
NAME= >
TYPE= #<FUNCTION-TYPE (FUNCTION
((OR
FLOAT
RATIONAL)
&REST
(OR
FLOAT
RATIONAL))
(MEMBER
T
NIL))>
WHERE-FROM= :DECLARED
KIND= :GLOBAL-FUNCTION>>
ARGS= (#<C::REF #xB0E9105
LEAF= #<C::LAMBDA-VAR #xB1091ED
NAME= I
TYPE= #<NUMERIC-TYPE NUMBER>>>
#<C::REF #xB0E908D
LEAF= #<C::LAMBDA-VAR #xB107955
NAME= N>>)>
#<FUNCTION-TYPE (FUNCTION
((UNSIGNED-BYTE 32)
(EXTENSIONS:CONSTANT-ARGUMENT
(UNSIGNED-BYTE 12)))
(MEMBER T NIL))>
:ARGUMENT-TEST
NIL
:RESULT-TEST
NIL
:STRICT-RESULT
NIL
:ERROR-FUNCTION
NIL
:WARNING-FUNCTION
NIL)
6] q
Compilation unit aborted.
2 errors
2 warnings
Compilation aborted after 0:00:55.
* (quit)
mccarthy 34% exit
mccarthy 35%
script done on Fri Jan 24 14:25:55 1992
*** File package.lisp ***
;;; -*- Mode: LISP; Syntax: Common-lisp; Package: USER; Base: 10 -*-
;;; Package hack for Lispm use
;;; Updated from Till's version as of 20 Nov 91
(defpackage poly
#+Symbolics
(:use future-common-lisp CLOS)
#+CMU
(:use pcl)
#+Symbolics
(:shadowing-import clos:documentation)
(:export
"2PI" "PI/2" "-3DB"
"POLY-+" "POLY-*" "POLY-/" "POLY-DERIV" "POWER-POLY" "POLY-NORMALIZE" "POLY-W-NORMALIZE"
"POLY" "POLY-W" "W-SCALE-POLY" "3DB-NORMALIZE" "FIND-CUTOFF-FREQUENCY"
"W-SCALE-POLES" "POLY-ZEROIZE" "F-AND-Q-FROM-POLES" "POLY-FROM-POLE-LOCATION" "POLE-FROM-W-AND-Q"
"SORT-POLES" "SUBSUME-SINGLE-POLES" "POLY-VALUE"
"BUTTERWORTH-POLES" "QUASI-BUTTERWORTH-POLES" "ELLIPTIC-POLES"
"CHEBYSHEV-POLES" "CHEBYSHEV-E-TO-K" "CHEBYSHEV-K-TO-E" "CHEBYSHEV-RIPPLE"
"CHEBYSHEV-RIPPLE-TO-E" "QUASI-NTH-ORDER-POLES" "BESSEL-POLES" "BESSEL-POLY"
"ELLIPTIC-CONTOUR-POLES" "GAUSSIAN-POLES" "SELECT-POLES" "SELECT-POLE-COMBINATIONS"
"POLY-POLES" "TRANSFORM-POLES-TO-HIGH-PASS" "GENERIC-ERROR" "WP" "BINARY" "FACTORIAL"
"TF-VALUE" "TF-*")
)
(defpackage speaker
(:use
#+Symbolics future-common-lisp
#+Symbolics CLOS
#+CMU pcl
poly)
#+Symbolics
(:shadowing-import clos:documentation)
)
*** File poly.lisp ***
;;; -*- Mode: LISP; Syntax: Common-lisp; Package: POLY; Base: 10 -*-
;;; Polynomial math tools. -- till
(in-package "POLY")
;;; Bertha is broken when safety = 3. On the Sun4 anyway.
#+Lucid
(proclaim '(optimize (speed 3) (compilation-speed 0) (safety 2)))
;;; Lucid CL only supports double floats. This was written in Lucid.
#+Symbolics
(eval-when (compile load eval)
(setq *read-default-float-format* 'double-float))
;;; A polynomial is a list of coefficients. Low order first, so it looks backwards.
(defconstant 2pi (* 2.0 pi))
(defconstant pi/2 (/ pi 2.0))
(defconstant -3db (sqrt 0.5))
(defun poly-+ (&rest polys)
(labels ((poly-+-internal (poly1 poly2)
(loop for (p1 . rest1) on (if (numberp poly1) (list poly1) poly1)
for (p2 . rest2) on (if (numberp poly2) (list poly2) poly2)
collect (+ p1 p2) into result
unless rest1 return (nconc result (copy-list rest2))
unless rest2 return (nconc result (copy-list rest1))
finally (return result))))
(cond
((null (cdr polys)) (car polys))
((cddr polys) (reduce #'poly-+-internal polys))
(t (poly-+-internal (first polys) (second polys))))))
;;; Optimized somewhat.
(defun poly-* (&rest polys)
(flet ((poly-*-internal (poly1 poly2)
(cond ((and (numberp poly1) (numberp poly2))
(* poly1 poly2))
((numberp poly1)
(loop for coef2 in poly2 collecting (* poly1 coef2)))
((numberp poly2)
(loop for coef1 in poly1 collecting (* poly2 coef1)))
(t (loop with result = (make-list (+ (length poly1) (length poly2) -1))
for coef1 in poly1
for sub-result on result
doing (loop for coef2 in poly2
for sub-sub-result on sub-result
as prod = (* coef1 coef2)
doing (if (car sub-sub-result)
(incf (car sub-sub-result)
prod)
(setf (car sub-sub-result) prod)))
finally (return result))))))
(cond
((null (cdr polys))
(car polys))
((cddr polys)
(reduce #'poly-*-internal polys))
(t
(poly-*-internal (first polys) (second polys))))))
(defun poly-/ (p1 p2)
(let* ((l1 (length p1))
(l2 (length p2))
(diff (- l1 l2)))
(when (minusp diff) (error "Can't do that, man"))
(loop with quotiant = nil
for offset from diff downto 0
as chunk = (nthcdr diff p1) then (cons (nth offset p1)
remainder)
as leftmost = (car (last chunk))
as remainder = (butlast
(poly-+ chunk
(poly-* (- leftmost) p2)))
doing (push leftmost quotiant)
finally (return (values quotiant (poly-zeroize remainder))))))
;;; Derivitive of a polynomial
(defun poly-deriv (poly)
(loop for coef in (cdr poly)
for exp from 1
collecting (* coef exp)))
;;; Given a voltage polynomial, return a power polynomial.
(defun power-poly (poly)
(poly-zeroize
(poly-* poly
(w-scale-poly poly -1))))
;;; Polynomial housekeeping functions.
;;; Normalizes a polynomial so that the leftmost coefficient (s^n) is 1.
(defun poly-normalize (poly)
(let ((leftmost (car (last poly))))
(if (= 1.0 leftmost)
poly
(poly-* poly (/ leftmost)))))
;;; Normalizes a polynomial so that the rightmost coefficient (s^0) is 1.
(defun poly-w-normalize (poly)
(w-scale-poly poly (/ (poly-w poly))))
(defun poly-w (poly)
(expt (/ (car poly) (car (last poly)))
(/ (1- (length poly)))))
;;; Scale a poly by w. Chucko's improved version.
(defun w-scale-poly (poly w)
(let ((w^i 1.0))
(labels ((w-scale-internal (pol)
(if (null (cdr pol))
(cons (car pol) nil)
(let ((so-far (w-scale-internal (cdr pol))))
(cons (* (car pol) (setf w^i (* w^i w)))
so-far)))))
(w-scale-internal poly))))
;;; Some filters coefficients aren't readily available scaled to
;;; a common -3db down point. Bessel comes to mind.
(defun 3db-normalize (poles)
(w-scale-poles poles (/ (find-cutoff-frequency poles))))
(defun find-cutoff-frequency (poles &optional (limit -3db))
(let ((poly (apply #'poly-* poles)))
(binary #'(lambda (w) (/ (car poly)
(abs (poly-value poly (complex 0 w)))))
limit .5 2.0)))
;;; Scale a list of poles by w.
(defun w-scale-poles (poles w)
(loop for pole in poles
collecting (w-scale-poly pole w)))
(defparameter *near-zero* 1e-10)
;;; Rounds out near zeros of a poly to zero and
;;; eliminates any zeros on the left.
(defun poly-zeroize (poly)
(loop for coef in poly
with new-poly = nil
doing (push (if (< (abs coef) *near-zero*) 0 coef) new-poly)
finally
(loop as first = (car new-poly)
while (when first (zerop first))
doing (pop new-poly))
(return (nreverse new-poly))))
;;; Transfer-function functions.
;;; A tf is a list of (poles zeros).
;;; Transfer function multiply.
(defun tf-* (&rest tfs)
(loop for (tf-num tf-den) in tfs
as num = tf-num then (poly-* num tf-num)
as den = tf-den then (poly-* den tf-den)
finally (return (list num den))))
;;; Transfer function value.
(defun tf-value (tf s)
(/ (poly-value (first tf) s)
(poly-value (second tf) s)))
;;; Return Fs and Qs in a list. If Q is nil, it's a single pole.
;;; Scale is in Hz.
(defun f-and-q-from-poles (poles scale)
(loop for pole in poles
collecting (f-and-q-from-pole pole scale)))
(defun f-and-q-from-pole (pole &optional scale)
(destructuring-bind (c0 c1 &optional c2 too-many)
(if scale
(w-scale-poly pole scale)
pole)
(when too-many (generic-error))
(cond
(c2
(when (/= c2 1.0) (warn "I was expecting c2 to be 1.0"))
(let ((w (sqrt (/ c0 c2))))
(list w
(/ w c1 c2))))
(t
(when (/= c1 1.0) (warn "I was expecting c1 to be 1.0"))
(list (/ c0 c1))))))
;;; Note well: This creates double poles for real values
;;; of p. Might not be what you want. (Does this sound
;;; the voice of experience or what.)
(defun poly-from-pole-location (p)
(list (expt (abs p) 2)
(* -2.0 (realpart p))
1))
(defun pole-from-w-and-q (w &optional q)
(if q
(list (expt w 2)
(/ w q)
1.0)
(list w 1.0)))
;;; Sorts poles to the preferred hi-q-first.
;;; Ignores non-double poles.
(defun sort-poles (poles)
(flet ((pole-q (pole)
(case (length pole)
(2 0)
(3 (let* ((c (pop pole))
(b (pop pole))
(a (pop pole))
(w (sqrt (/ c a)))
(q (/ w b a)))
q))
(t (generic-error)))))
(sort poles #'(lambda (a b) (> (pole-q a) (pole-q b))))))
;;; Fold single poles into low-q double poles.
;;; I basically ignore the situation where there are more
;;; than two single poles. The absolutely correct approach
;;; would be to work on each combination, but I don't
;;; believe this situation ever comes up in practice.
;;;
(defun subsume-single-poles (poles)
(let ((single-poles nil)
(double-poles nil))
(dolist (pole poles)
(case (length pole)
(2 (push pole single-poles))
(3 (push pole double-poles))
(t (generic-error))))
(case (length single-poles)
((0 1) poles)
(2 (nconc double-poles
(list (poly-* (first single-poles)
(second single-poles)))))
(t (warn "There are multiple ways to group these poles.")
(format t "~% sp :~s, ~%dp:~s" single-poles double-poles)
(nconc double-poles
(loop for (first second) on single-poles by #'cddr
collecting (if second
(poly-* first second)
first)))))))
;;; Find the value of a polynomial at a given s.
;;; Note low-error algorithm
(defun poly-value (poly s)
(flet ((poly-value-internal (poly s)
(if (cdr poly)
(+ (car poly)
(* s (poly-value (cdr poly) s)))
(car poly))))
(poly-value-internal poly s)))
(defun butterworth-poles (n-poles)
(elliptic-poles n-poles 1.0))
;;; Return a list of left-hand poles on the ellipse.
;;; For Butterworth and Chebyshev filters.
(defun elliptic-poles (n-poles k)
(loop for theta from (/ pi/2 n-poles) by (/ pi n-poles)
repeat (floor n-poles 2)
as xtemp = (* (sin theta) k)
as ytemp = (* (cos theta))
collecting (list (if (= k 1) 1.0 (+ (expt xtemp 2) (expt ytemp 2)))
(* 2 xtemp)
1.0)
into polys
finally (return (if (oddp n-poles)
(nconc polys (list (list k 1.0)))
polys))))
(defun chebyshev-poles (n-poles e)
(let ((k (chebyshev-e-to-k n-poles e)))
(if (minusp k)
(3db-normalize (elliptic-poles n-poles (/ -1.0 k)))
(elliptic-poles n-poles k))))
;;; The epsilon and ripple computations below are derrivable from
;;; Lindquist pp 217-218. Note the typo on p 218.
;;; Those kludges to get around dividing by zero are just kludges.
(defun chebyshev-e-to-k (n-poles e)
(if (< #-Symbolics -1.0e-154
#+Symbolics -1.488e-8
e 1.0e-154)
1.0
(tanh (/ (asinh (/ e)) n-poles))))
(defun chebyshev-k-to-e (n-poles k)
(if (< (abs (- k 1.0)) .0000001)
0.0
(/ (sinh (* n-poles (atanh k))))))
;;; Ripple is total peak-to-peak.
(defun chebyshev-ripple (e)
(when (plusp e)
(* 20.0 (log (+ 1.0 (expt e 2)) 10))))
(defun chebyshev-ripple-to-e (ripple)
(when (plusp ripple)
(sqrt (- (expt 10.0 (/ ripple 20.0)) 1.0))))
;;; Quasi Nth Order filters.
(defun quasi-butterworth-poles (n B)
(sort-poles
(subsume-single-poles
(3db-normalize
(poly-poles
(cons 1.0
(nconc (make-list (1- (* 2 (1- n))) :initial-element 0.0)
(list (* (expt B 2) (if (evenp n) -1.0 1.0))
0.0
(if (evenp n) 1.0 -1.0))))
t)))))
;;; Bessel Filter
(defun bessel-poles (n)
(3db-normalize
(sort-poles
(poly-poles (poly-w-normalize (bessel-poly n))))))
;;; (Bessel Poly meets Barney Kessel's Folly Wrestling at Cal Poly. Golly!)
;;; Generates a bessel polynomial
(defun bessel-poly (n)
(case n
(0 1)
(1 '(1 1))
(t (poly-+
(poly-* (1- (* 2 n)) (bessel-poly (1- n)))
(poly-* '(0 0 1) (bessel-poly (- n 2)))))))
;;; Lundquist, p265. He sez this is Bessel when k=1. He also sez it isn't,
;;; so we don't really trust him.
(defun elliptic-contour-poles (n-poles k)
(3db-normalize
(loop for y from (if (evenp n-poles) (/ n-poles) 0) by (/ 2.0 n-poles)
repeat (ceiling n-poles 2)
with poles = nil
doing
(push (if (zerop y)
(list k 1)
(let ((x (/ (sqrt (- 1.0 (expt y 2))) k)))
(list (+ (expt x 2) (expt y 2)) (* x 2) 1)))
poles)
finally (return poles))))
;;; From Lundquist.
(defun gaussian-poles (n)
(3db-normalize
(subsume-single-poles
(sort-poles
(poly-poles
(poly-normalize
(cons 1.0
(loop with a = -0.694
for i from 1 to n
collect 0.0
collect (/ (expt a i) (factorial i)))))
:powerp)))))
;;; Generates a list of positition-bit combinations for n-poles.
(defun select-pole-combinations (choose-n-pole-pairs n-poles)
(loop for i below (expt 2 (truncate n-poles 2))
when (= choose-n-pole-pairs (logcount i))
collect i))
;;; Select some of the poles from a pole-list and collect them
;;; into a polynomial.
;;; The bits in which-poles correspond to each of the polys,
;;; lsb is first.
;;; Return two values: the chosen poles, and the other poles.
(defun select-poles (poles &optional (which-poles -1))
(loop for pole in poles
for m = which-poles then (ash m -1)
as pole-length = (length pole)
unless (or (= pole-length 3)
(and (= pole-length 2) (evenp m)))
do (error "What's with this pole ~s?" pole)
when (oddp m)
collect pole into chosen-poles
else collect pole into other-poles
finally (return (values chosen-poles other-poles))))
;;; Reflects around ohmega=1, of course.
(defun transform-poles-to-high-pass (poles)
(loop for pole in poles
as length = (length pole)
collecting (cond ((= length 2)
(unless (= (second pole) 1)
(error "I didn't expect a pole like this"))
(list (/ (car pole)) 1.0))
((= length 3)
(destructuring-bind (c0 c1 c2) pole
(unless (= 1 c2)
(error "I didn't expect a pole like this"))
(list (/ c0) (/ c1 c0) 1.0)))
(t (generic-error)))))
;;; Just a tad west of north.
(defparameter *start* (cis (* pi .503)))
;;; Find the poles of the given polynomial using
;;; recursion and successive approximation.
;;; All in one. Is this too ugly to use? I feel like a goddamned
;;; scheme programmer. I have to admit that it works well here.
;;;
;;; Powerp indicates that this is a power polynomial, in that
;;; case we only return poles on the left hand side.
(defun poly-poles (original-poly &optional powerp)
(let* ((near-zero *near-zero*)
(original-poly-d (poly-deriv original-poly))
(original-poly-dd (poly-deriv original-poly-d)))
(labels
((next-aprox (poly poly-d poly-dd x)
(let* ((f (poly-value poly x))
(fp (poly-value poly-d x))
(fpp (poly-value poly-dd x)))
;; Typically the parabolic aproximation will be fine.
;; However, it won't work when the second derivitive nears
;; zero. I'm not gonna worry yet.
(if t
(let* ((squirt (sqrt (- (expt fp 2) (* 2 f fpp))))
(one (- x (/ (+ fp squirt) fpp)))
(two (- x (/ (- fp squirt) fpp))))
(if (< (abs (poly-value poly one))
(abs (poly-value poly two)))
one
two))
(- x (/ f fp)))))
(poly-rooter (roots-so-far roots-poly)
(let ((poly (if roots-so-far
(poly-/ original-poly roots-poly)
original-poly)))
; (format t "~%Roots so far: ~S,~%(w's)~s ~%Roots poly: ~s~%Poly: ~s:~%"
; roots-so-far
; (mapcar #'poly-w roots-so-far)
; roots-poly poly)
(case (length poly)
((0 1) roots-so-far)
((2 3) (if powerp
(cons (list (sqrt (- (/ (first poly)
(third poly)))) 1.0)
roots-so-far)
(cons poly roots-so-far)))
(t (let* ((poly-d (poly-deriv poly))
(poly-dd (poly-deriv poly-d)))
(labels
((recursive-find-root (point)
(let ((next (next-aprox poly poly-d poly-dd point)))
;;; (format-complex "~%X= " point)
;;; (format t "abs Y: ~d" (abs (poly-value poly next)))
(if (or (< (abs (poly-value poly next))
near-zero)
(< (/ (abs (- next point))
(+ near-zero (abs next)))
near-zero))
next
(recursive-find-root next)))))
(let* ((pole (next-aprox original-poly
original-poly-d
original-poly-dd
(recursive-find-root *start*)))
(pole-poly (if (< (/ (abs (imagpart pole))
(abs pole))
1e-6)
(list (- (realpart pole)) 1.0)
(poly-from-pole-location pole))))
;; Recurse on smaller polynomial
(if powerp
(let ((-pole-poly (w-scale-poly pole-poly -1.0)))
(poly-rooter
(cons (if (plusp (poly-w pole-poly))
pole-poly
-pole-poly)
roots-so-far)
(poly-* pole-poly -pole-poly roots-poly)))
(poly-rooter (cons pole-poly roots-so-far)
(poly-* pole-poly roots-poly)))))))))))
(poly-rooter nil (list 1)))))
;;; Debugging tools:
(defun generic-error ()
(error "Huh?"))
;;; Nicely print the frequencies, q's, and coordinates of a list of poles.
;;; Too good to only use for debugging.
(defun wp (poles)
(dolist (pole poles)
(unless (= (car (last pole)) 1.0)
(warn "Non-normalized pole")
(setq pole (poly-normalize pole)))
(case (length pole)
(2 (let ((w (- (car pole))))
(format t "~% x=~5,2f, f=~5,2f single real pole"
w (/ w 2pi))))
(3 (let* ((c (pop pole))
(b (pop pole))
(a (pop pole))
(w (sqrt (/ c a)))
(f (/ w 2pi))
(q (/ (* w a) b))
(det (- (expt b 2) (* 4 a c)))
(x (- (/ b 2 a))))
(cond ((minusp det)
(let ((y (/ (sqrt (- det)) 2 a)))
(format t
"~% xy=~5,2f +/-~5,2f ~12@t; w=~6,4f, f=~6,4f, q=~5,2f [+/~5,2f deg], del=~6,4f"
x y w f q (* 180.0 (/ (atan (/ y x)) pi)) (/ 0.5 q))))
((plusp det)
(format t "~% x=~5,2f, ~5,2f, two real poles; w=~5,2f, f=~6,4f, q=~5,2f"
(+ x (/ (sqrt det) 2 a))
(- x (/ (sqrt det) 2 a))
w (/ w 2pi) q))
(t
(format t "~% x=~5,2f, double real pole; w=~5,2f, f=~5,2f, q=~5,2f"
x w f q)))))
(t (format t "Huh? ~s" pole))))
(values)) ; don't return any values
;;; For debugging. No longer necessary.
(defun format-complex (text c)
(format t text)
(when c
(format t "~6,3,,'*f,~6,3,,'*f " (realpart c) (imagpart c))))
(defparameter *epsilon-fraction* 0.0001) ; Fraction of one
(defvar *binary-search-debug* nil)
;;; Binary search for goal.
;;; This is broken for large bessels. I don't know why yet.
(defun binary (func goal mini maxi &optional (nil-outside-p t))
(let ((min (float mini))
(max (float maxi)))
(declare (float min max))
(let ((ymin (funcall func min))
(ymax (funcall func max))
(mid 0.0)
(ymid 0.0)
(epsilon-fraction *epsilon-fraction*))
(declare (float ymin ymax mid ymid epsilon-fraction)
(ftype (function (float float) float) + - * /)
(ftype (function (float float float) float) + *))
(loop repeat 20
as first-time = nil then t
doing
(when first-time
(cond ((< ymin ymax ymid)
(if nil-outside-p
(return nil)
(setq min max max mid ymin ymax ymax ymid)))
((< ymid ymin ymax)
(if nil-outside-p
(return nil)
(setq max min min mid ymax ymin ymin ymid)))
((or (< ymin goal ymid) (> ymin goal ymid))
(setq max mid ymax ymid))
(t (setq min mid ymin ymid))))
;; place mid intelligently between min and max
(setq mid
(+ min (* (- max min)
(/ (- goal ymin) (- ymax ymin)))))
(setq ymid (the float (funcall func mid)))
(when *binary-search-debug*
(format t "~%min:~6,3f,~6,3f, mid:~6,3f,~6,3f, max:~6,3f,~6,3f"
min ymin mid ymid max ymax))
(when (< (abs (- ymid goal))
(* goal epsilon-fraction))
(return mid))))))
(defun factorial (n)
(do ((i 1 (1+ i))
(result 1 (* result i)))
((> i n) result)))
*** end ***