[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: compiling problems in 2.0 final



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))))