(module csl.poly (poly-eval*
                  poly-eval-derivs*
                  solve-quardatic*
                  solve-cubic*)
  (import scheme
          (only chicken.module reexport)
          gsl.poly)
  (reexport gsl.poly)


  (define (poly-eval* c z)
    (do ((const c (cdr const))
         (pow 0 (+ pow 1))
         (val 0 (+ val (* (car const) (expt z pow)))))
        ((null? const) val)))

  (define (poly-eval-derivs* c x #!optional (num-derivs (length c)))
    (define (next-deriv const)
      (if (null? (cdr const))
          '(0)
          (let loop ((c (cdr const))
                     (pow 1))
            (if (null? c)
                '()
                (cons (* pow (car c))
                      (loop (cdr c) (+ pow 1)))))))
    (cons
     (poly-eval* c x)
     (let loop ((num num-derivs)
                (c (next-deriv c)))
       (if (zero? (- num 1))
           '()
           (let ((nd (next-deriv c)))
             (cons (poly-eval* nd x)
                   (loop (- num 1) (next-deriv c))))))))

  (define (solve-quadratic* a b c)
    (define (conj x)
      (if (cplxnum? x)
          (make-rectangular (real-part x) (- (imag-part x)))
          x))
    (let* ((disc (sqrt (- (expt b 2) (* 4 a c))))
           (sign (if (>= 0 (real-part (* (conj b) disc))) 1 -1))
           (q (* -1/2 (+ b (* sign disc)))))
      (values (/ q a) (/ c q))))

  (define (solve-cubic* a b c)
    (let* ((q (- (* a a) (* 3 b)))
           (r (+ (* 2 a a a) (* -9 a b) (* 27 c)))
           (Q (/ q 9))
           (R (/ r 54))
           (Q3 (* Q Q Q))
           (R2 (* R R))
           (CR2 (* 729 r r))
           (CQ3 (* 2916 q q q)))
      (cond ((and (= R 0) (= Q 0))
             (let ((res (/ a -3)))
               (values res res res)))
            ((= CR2 CQ3)
             (let ((sqrtQ (sqrt Q)))
               (if (> R 0)
                   (values (- (* -2 sqrtQ) (/ a 3))
                           (- sqrtQ (/ a 3))
                           (- sqrtQ (/ a 3)))
                   (values (- (- sqrtQ) (/ a 3))
                           (- (- sqrtQ) (/ a 3))
                           (- (* 2 sqrtQ) (/ a 3))))))
            ((< R2 Q3)
             (let* ((sgnR (if (>= R 0) 1 -1))
                    (ratio (* sgnR (sqrt (/ R2 Q3))))
                    (theta (acos ratio))
                    (norm (* -2 (sqrt Q)))
                    (results (list (- (* norm (cos (/ theta 3)))
                                      (/ a 3))
                                   (- (* norm (cos (/ (+ theta (* 2 pi)) 3))) (/ a 3))
                                   (- (* norm (cos (/ (- theta (* 2 pi)) 3))) (/ a 3)))))
               (apply values (sort results <))))
            (else
             (let* ((sgnR (if (>= R 0) 1 -1))
                    (A (- (* sgnR (expt (+ (abs R) (sqrt (- R2 Q3))) (/ 1 3)))))
                    (B (/ Q A)))
               (if (< (+ A B) 0)
                   (values (+ A B (- (/ a 3)))
                           (make-rectangular (- (+ (/ (+ A B) 2) (/ a 3)))
                                             (* -1 (sqrt 3) 1/2 (abs (- A B))))
                           (make-rectangular (- (+ (/ (+ A B) 2) (/ a 3)))
                                             (* (sqrt 3) 1/2 (abs (- A B)))))
                   (values (make-rectangular (- (+ (/ (+ A B) 2) (/ a 3)))
                                             (* -1 (sqrt 3) 1/2 (abs (- A B))))
                           (make-rectangular (- (+ (/ (+ A B) 2) (/ a 3)))
                                             (* (sqrt 3) 1/2 (abs (- A B))))
                           (+ A B (- (/ a 3)))))))))))
