chicken-users
[Top][All Lists]
Advanced

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

[Chicken-users] Fwd: gettings roots (and multiplier) of a polynomial equ


From: Terrence Brannon
Subject: [Chicken-users] Fwd: gettings roots (and multiplier) of a polynomial equation and vice versa + Gnu Scientific Library
Date: Sun, 11 Nov 2007 19:26:33 -0500

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 <address@hidden>
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 <address@hidden>


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




reply via email to

[Prev in Thread] Current Thread [Next in Thread]