polynomials 
Author Message
 polynomials

Hello,

I am translating (part of) Stephen Moshier's Cephes Library from C to
CL. Since a lot of this is polynomial approximation, I need a really
efficient version of (poly x y) which implements Horner's rule. Here
x is the argument of the polynomial, and y is a list or vector of
coefficients (4-10 double-floats). I dont want to use assembler, because
of portability (yet). Some questions:
(a) Is there such code -- I have CL Math, but that seems tentative.
(b) Do I use a list or vector of coefficients.
(c) Is it advantageous to defconstant the coefficients, or do I just
pass the actual values to poly.
(d) Where do I put the declarations.

This is what I use now (from CL Math, but not using a macro)

(defun poly (var coefs)
(declare (double-float var) (type (list double-float *) coefs)
         (optimize (safety 0) (space 0) (speed 3)))
(cond ((null coefs) 0.0)
      ((null (cdr coefs)) (car coefs))
      (t (the double-float (+ (the double-float (car coefs))
         (the double-float (* (the double-float var)
                              (the double-float (poly var (cdr coefs))))))))))

and here is an example of a function using poly

(defun stirf (x)
(declare (double-float x))
(let* (
      (a (/ x))
      (w (+ 1.0 (* a (poly a '(8.33333333333482257126e-2
                               3.47222221605458667310e-3
                              -2.68132617805781232825e-3
                              -2.29549961613378126380e-4
                               7.87311395793093628397e-4)))))
      (y (exp x))
      )
  (declare (double-float a w y))
(if (> x 143.01608)
  (let* (
        (v (expt x (- 0.25 (* 0.5 x))))
        (r (* v (/ v y)))
        )
    (declare (double-float v r))
    (* w r 2.50662827463100050242e0))
  (let (
        (r (/ (expt x (- x 0.5)) y))
        )
    (declare (double-float r))
    (* w r 2.50662827463100050242e0)))
))

All suggestions are welcome.
The stronger the criticism, the more welcome it is.

--- Jan

Jan de Leeuw; Dept Math UCLA; phone (310)-825-9550; fax (310)-206-6673      
mail: 405 Hilgard Ave, Los Angeles, CA 90024-1555

--
Jan de Leeuw; Dept Math UCLA; phone (310)-825-9550; fax (310)-206-6673      
mail: 405 Hilgard Ave, Los Angeles, CA 90024-1555



Fri, 09 Aug 1996 03:33:09 GMT  
 polynomials

Quote:

>Hello,

Here's one attempt at a faster version, which reverses the order
of coefficients. It uses iteration instead of recursion; the version
of poly that you posted wasn't tail recursive; you could make it
tail-recusive by passing an accumulator but iteration seemed a more natural
way of expressing this problem.

Be careful about changing the safety setting - remember 0.0 isn't a
double-float. For yet more efficiency, you could turn poly into a macro;
just backquote and comma the code computing result.

(defun poly-2 (var coefs)
  "compute a polynomial in var using coefficients in coefs. Coefficients are
   in decending order - (poly-2 x (1 2 3 4)) == x^3 + 2x^2 + 3x +4"
  (declare (optimize  speed safety)
           (double-float var))
  (do* ((result 0.0d0 (+ c (* result var)))
        (coefficients coefs (cdr coefficients))
        (c (car coefficients) (or (and coefficients (car coefficients))
                                  0.0d0)))
       ((null coefficients) result)
       (declare (double-float result c))))
--
Hackers Local 42- National Union of Computer Operatives, Chapel Hill section
------------------------------------------------------------------------------
Tar Heel Information Services - Nothing but net!   | WAIS/Z39.50 spoken here
North Carolina - First in Usenet        | DoD #612 | Tel: +1-919-962-9107



Sat, 10 Aug 1996 02:05:19 GMT  
 polynomials

Quote:

>You iterative horner method is better than the recursive method because the
>compiler can avoid number consing.  Generally, you'll want horner to be a
>a macro.

A good compiler should be able to avoid boxing on a tail recursive
version; an interesting question would be how many compilers would  
generate an unboxed entry point for the non-tail recursive version. The
assembly code generated by cmu-cl for the macro version was as good as
I initially expected; I guess the reason is the old float != real
problem.

I assume that MCL does the right thing, judging from the quality of people
at Apple's Cambridge office; I just wish that they'd create a native version
of MCL for the PowerMac. That'd make for a killer home setup. As it is I'm
going for a classical approach and getting a couple of old Symbolics
boxes from Texas; are there any lispm mailing lists still going?

Simon

--
Hackers Local 42- National Union of Computer Operatives, Chapel Hill section
------------------------------------------------------------------------------
Tar Heel Information Services - Nothing but net!   | WAIS/Z39.50 spoken here
North Carolina - First in Usenet        | DoD #612 | Tel: +1-919-962-9107



Tue, 13 Aug 1996 10:05:31 GMT  
 
 [ 3 post ] 

 Relevant Pages 

1. What order of polynomial in "polynomial interpolation.vi"?

2. LaGrange Polynomials -- Why use them?

3. polynomials in J

4. Modular polynomial arithmetic

5. Multivariate polynomials- a puzzle

6. Multivariate polynomial puzzle

7. Polynomial zeros

8. least squares fitting to a polynomial

9. 2-dimensional polynomial fit

10. polynomial and finite field code in Caml, Haskell, Clean

11. matrix, polynomial and integer code

12. code for linear algebra and polynomials

 

 
Powered by phpBB® Forum Software