Does anyone know of a clisp library that has functions for getting the
roots (Actually numerical aproximations to the roots.) of a polynomial
with complex coefficients?
I.e. a function that could find the approximate roots of the following
polynomial when it is multiplied out: (x-4I)(x-5+6I)(x+2)(x+2+2I)
Regards,
Henrik Treadup
If I had to write one, I'd start with the routines from _Numerical
Recipes_ and translate them to Lisp. I'd convert the real x,y pairs to
Lisp's complex data type, and use lots of declarations like
(declare (complex z w ...)).
The standard advice on _Numerical Recipes_ is that experts on the topic
of its Chapter A will say the routines are oversimplified and limited,
but that the same experts will be happy to use the routines in Chapter
B because they're so convenient. Here A,B ranges over all pairs of
distinct topics.
···············@yahoo.com wrote:
> If I had to write one, I'd start with the routines from _Numerical
> Recipes_ and translate them to Lisp. I'd convert the real x,y pairs to
> Lisp's complex data type, and use lots of declarations like
> (declare (complex z w ...)).
>
> The standard advice on _Numerical Recipes_ is that experts on the topic
> of its Chapter A will say the routines are oversimplified and limited,
> but that the same experts will be happy to use the routines in Chapter
> B because they're so convenient. Here A,B ranges over all pairs of
> distinct topics.
Last year I did just that. The code is in PLT Scheme, but it should
be straight forward to port it to CL -- just as this is a
straight forward port from C to Scheme. Needless to say, this program
is not representative of my normal Scheme programming style.
;;; polynomial-roots.scm -- Jens Axel S�gaard
(module polynomial-roots mzscheme
(provide (all-defined))
; roots : (list numbers) -> (list numbers)
; return list of roots of the polynomial with
; coeffecients from the argument list, lowest degree first
(define (roots coeffs)
(let* ([degree (- (length coeffs) 1)]
[roots (make-vector degree 0.0)])
(polynomial-roots-all (list->vector coeffs) degree roots)
(mergesort (vector->list roots)
(lambda (x y)
(or (< (real-part x) (real-part y))
(and (= (real-part x) (real-part y))
(< (imag-part x) (imag-part y))))))))
(define (polynomial-roots-laguerre a m x)
; Given the complex coefficients in the vector a of the
; polynomial (sum (i 0 m) a_i * x^i ) and a complex value
; x this routine improves x by Laguerre's method until
; it converges to a root of the polynoial.
(define-syntax vr
(syntax-rules () [(_ v i) (vector-ref v i)]))
(define-syntax vs!
(syntax-rules () [(_ v i x) (vector-set! v i x)]))
(define-syntax define*
(syntax-rules ()
[(_ () v) (begin)]
[(_ (id ids ...) v) (begin (define id v) (define* (ids ...) v))]))
(define-syntax := (syntax-rules () [(_ id v) (set! id v)]))
(define* (iter its j MAXIT) 1)
(define* (abx abp abm err) 0.0)
(define* (dx x1 b d f g h sq gp gm g2) 0.0)
(define frac (vector 0.0 0.5 0.25 0.75 0.13 0.38 0.62 0.88 1.0))
(define MR 8)
(define MT 10)
(define EPSS 1.0e-7)
(set! MAXIT (* MT MR))
(let ([return #f])
(do ([iter 1 (+ iter 1)]) [(or return (> iter MAXIT))]
(:= its iter)
(:= b (vr a m))
(:= err (abs b))
(:= f 0.0+0.0i)
(:= d f)
(:= abx (magnitude x))
(do ([j (- m 1) (- j 1)]) [(or return (< j 0))]
(:= f (+ (* x f) d))
(:= d (+ (* x d) b))
(:= b (+ (* x b) (vr a j)))
(:= err (+ (magnitude b) (* abx err))))
(:= err (* err EPSS))
(if (<= (magnitude b) err)
(set! return x)
(begin
(:= g (/ d b))
(:= g2 (* g g))
(:= h (- g2 (* 2.0 (/ f b))))
(:= sq (sqrt (* (- m 1) (- (* m h) g2))))
(:= gp (+ g sq))
(:= gm (- g sq))
(:= abp (magnitude gp))
(:= abm (magnitude gm))
(if (< abp abm)
(:= gp gm))
(:= dx (if (> (max abp abm) 0.0)
(/ m gp)
(* (+ 1 abx)
(+ (cos iter)
(* (sin iter) 0.0+1.0i)))))
(:= x1 (- x dx))
(if (= x x1)
(set! return x)
(begin
(if (= 0 (modulo iter MT))
(:= x x1)
(:= x (- x (* (/ iter MT) dx)))))))))
(if return
return
(error "polynomial-roots-laguerre: too many iterations - try
another guess"))))
(define (polynomial-roots-all a m roots)
; Given the complex coefficients in the vector a of the
; polynomial (sum (i 0 m) a_i * x^i ) the degree m
; (if a is longer than m the end of a is ignored)
; this routine fills the vector roots with
; the roots
(define-syntax vr
(syntax-rules () [(_ v i) (vector-ref v i)]))
(define-syntax vs!
(syntax-rules () [(_ v i x) (vector-set! v i x)]))
(define-syntax define*
(syntax-rules ()
[(_ () v) (begin)]
[(_ (id ids ...) v) (begin (define id v)
(define* (ids ...) v))]))
(define-syntax := (syntax-rules () [(_ id v) (set! id v)]))
(define-syntax -- (syntax-rules () [(_ n) (sub1 n)]))
(define polish #t)
(define EPS 2.0e-6)
(define MAXM 100)
(define* (i its j jj) 0)
(define* (x b c) 0.0)
(define ad 0)
(set! ad (make-vector MAXM 0.0))
(do ([j 0 (+ j 1)]) [(> j m)]
(vs! ad j (vr a j)))
(do ([j m (- j 1)]) [(< j 1)]
(:= x 0.0+0.0i)
(:= x (polynomial-roots-laguerre ad j x))
(if (<= (abs (imag-part x)) (* 2.0 EPS (abs (real-part x))))
(:= x (real-part x)))
(vs! roots (-- j) x)
(:= b (vr ad j))
(do ([jj (- j 1) (- jj 1)]) [(< jj 0)]
(:= c (vr ad jj))
(vs! ad jj b)
(:= b (+ (* x b) c))))
(if polish
(do ([j 1 (+ j 1)]) [(> j m)]
(vs! roots (-- j)
(polynomial-roots-laguerre a m
(vr roots (-- j)))))))
)
There is a lisp implementation of
approximate rootfinding in maxima...
from the documentation,
ALLROOTS (poly)
finds all the real and complex roots of the real polynomial poly
which must be univariate and may be an equation, e.g. poly=0.
For complex polynomials an algorithm by Jenkins and Traub is used
(Algorithm 419, Comm. ACM, vol. 15, (1972), p. 97). For real
polynomials the algorithm used is due to Jenkins (Algorithm 493,
TOMS, vol. 1, (1975), p.178). .....
The code is in lisp, open source, etc. on sourceforge.
It could be extracted for re-use. e.g.
www.cs.berkeley.edu/~fateman/papers/poles.pdf
Maxima has other root finding programs too.
RJF
···········@gmail.com wrote:
> Does anyone know of a clisp library that has functions for getting the
> roots (Actually numerical aproximations to the roots.) of a polynomial
> with complex coefficients?
>
> I.e. a function that could find the approximate roots of the following
> polynomial when it is multiplied out: (x-4I)(x-5+6I)(x+2)(x+2+2I)
>
> Regards,
> Henrik Treadup
>
Since no one else has responded ... (I am no expert
in this area)
The only thing that I am aware of is Maxima
See ... http://maxima.sourceforge.net/
It is probably overkill for what your are looking for,
but I think Maxima should be able to find the exact roots
of that polynomial. (though you want the approximate roots)
Wade
On 9225 day of my life Wade Humeniuk wrote:
> It is probably overkill for what your are looking for,
> but I think Maxima should be able to find the exact roots
> of that polynomial. (though you want the approximate roots)
There is no method for calculating exact root of polynomial with
degree > 5. That is why the only universal method is to find
approximate roots.
--
Ivan Boldyrev
Outlook has performed an illegal operation and will be shut down.
If the problem persists, contact the program vendor.