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
