Jens wrote some code for this in case there is nothing in Chicken for it at the moment.
Only thing is doesnt do is return the multiplier for the factorization. ---------- Forwarded message ---------- From: Jens Axel Søgaard <[EMAIL PROTECTED]> Date: Nov 11, 2007 6:44 PM Subject: Re: gettings roots (and multiplier) of a polynomial equation and vice versa + Gnu Scientific Library To: Terrence Brannon <[EMAIL PROTECTED]> [Feel free to forward this to the Chicken mailing list (I would myself if I were subscribed)] Terrence Brannon wrote: > I'm wondering if I overlooked an egg that can do this: > > 1 - given the (real) coefficients of all terms of a polynomial > equation, it returns the roots and multiplier representation of it: > > For example, the polynomial 2x + 2x^2 factors down to 2 * (x-(-1))(x-0), > so the multiplier is 2 and the roots are 0 and -1 Here is an implementation of a method from "Numerical Recipes in C": Your example is solved like this: > (roots '(0.0 2.0 2.0)) (-1.0 0.0) Numerical code is tricky, so I wrote a literal translation. It ain't pretty, but it works. /Jens Axel Søgaard ; 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) (sort (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))))))) _______________________________________________ Chicken-users mailing list Chicken-users@nongnu.org http://lists.nongnu.org/mailman/listinfo/chicken-users