[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: compiling problems in 2.0 final
- To: hkt@zkm.de
- Subject: Re: compiling problems in 2.0 final
- From: "Mark A. Tapia" <markt@dgp.toronto.edu>
- Date: Tue, 6 Oct 1992 10:22:11 -0400
- Cc: info-mcl@cambridge.apple.com
Rick Taube (hkt@zkm.de) writes about problems compiling the following
(I tried it and it seems to loop):
(in-package :cl-user)
(defun bessi0 (x)
;;I0(x) from "Numerical Recipes" by Press et al
(if (< (abs x) 3.75)
(let* ((y (expt (/ x 3.75) 2)))
(+ 1.0
(* y (+ 3.5156229
(* y (+ 3.0899424
(* y (+ 1.2067492
(* y (+ 0.2659732
(* y (+ 0.360768e-1
(* y
0.45813e-2)))))))))))))
(let* ((ax (abs x))
(y (/ 3.75 ax)))
(* (/ (exp ax) (sqrt ax))
(+ 0.39894228
(* y (+ 0.1328592e-1
(* y (+ 0.225319e-2
(* y (+ -0.157565e-2
(* y (+ 0.916281e-2
(* y (+ -0.2057706e-1
(* y (+
0.2635537e-1
(* y (+
-0.1647633e-1
(
* y 0.392377e-2))))))))))))))))))))
---> Answer
You're really using Horners method for expanding a polynomial of
one variable. The following is clearer and actually compiles:
(defun horner (x coefs)
;; compute the polynomial x * n with coefficients xn ... x0
(loop for factor in coefs
with sum = 0.0
do (setq sum (+ (* x sum) factor))
finally (return sum)))
(defun bessi0 (x)
;;I0(x) from "Numerical Recipes" by Press et al
(if (< (abs x) 3.75)
(horner (expt (/ x 3.75) 2)
'(0.45813e-2 0.360768e-1 0.2659732 1.2067492 3.0899424 3.5156229 1))
(horner (/ 3.75 (abs x))
'( 0.392377e-2 -0.1647633e-1 0.2635537e-1 -0.2057706e-1 0.916281e-2 -0.157565e-2
0.225319e-2 0.1328592e-1 0.39894228 0.39894228))))