(import (rnrs)
        (durand-kerner))

(define (vectors-index2 proc a b)
  (let ((n (vector-length a)))
    (let loop ((i 0))
      (cond ((>= i n) #f)
            ((proc (vector-ref a i)
                   (vector-ref b i))
             i)
            (else (loop (+ i 1)))))))

;; (define (~ a b)
;;   (< (abs (- a b)) 1e-15))

(define (~ a b)
  (let ((c (- a b)))
    (and (< (abs (real-part c)) 1e-10)
         (< (abs (imag-part c)) 1e-10))))

(define (vectors~ a b)
  (not
   (vectors-index2 (lambda (x y)
                     (not (~ x y)))
                   a
                   b)))

;; (define (vectors~ a b)
;;   (let ((n (vector-length a)))
;;     (let loop ((i 0))
;;       (cond ((>= i n) #t)
;;             ((~ (vector-ref a i)
;;                 (vector-ref b i))
;;              (loop (+ i 1)))
;;             (else #f)))))

(assert
 (vectors~
  (roots
   (vector 1.0+0.0i -3.0+0.0i 2.0+0.0i))
  '#(1.0-0.0i 2.0+0.0i)))

(assert
 (vectors~
  (roots
   (vector 1 -3 2))
  '#(1.0-0.0i 2.0+0.0i)))

(assert
 (vectors~
  (roots
   (vector 1 -3 2))
  '#(1 2)))

(assert
 (vectors~
  (roots
   (vector 1.0+0.0i -2.0+0.0i 1.0+0.0i))
  '#(1.0-0.0i 1.0-2.220446049250313e-16i)))

(assert
 (vectors~
  (roots
   (vector 1.0+0.0i 0 4.0+0.0i))
  '#(0.0-2.0i 0.0+2.0i)))

(assert
 (vectors~
  (roots
   (vector 4.0+0.0i 0 1.0+0.0i))
  '#(0.0-0.5i 0.0+0.5i)))

(assert
 (vectors~
  (roots
   (vector 1.0+0.0i -3.0+0.0i 2.0+0.0i 0))
  '#(1.0+0.0i 2.0+0.0i 0.0+0.0i)))

(assert
 (vectors~
  (roots
   (vector 1 0 -13 0 36))
  '#(2.0-5.994355921380061e-33i
     3.0-1.7219155529623352e-41i
     -2.0+3.4258111894079275e-17i
     -3.0-6.595538179953504e-17i)))

(assert
 (vectors~
  (roots
   (vector 1 0 -13 0 36 0))
  '#(-2.0+7.174648137343064e-43i
     2.9999999999999996+1.2037062152420224e-34i
     -3.0-1.4349296274686127e-42i
     2.350988701644575e-38-2.350988701644575e-38i
     2.0+1.6851887013388314e-34i)))

(assert
 (vectors~
  (roots
   (vector 1 1 1))
  '#(-0.5-0.8660254037844387i -0.5+0.8660254037844386i)))

(assert
 (vectors~
  (roots
   (vector 1 1 1 1))
  '#(-8.481300776135252e-18-1.0i 8.475424015331716e-18+1.0i -1.0+0.0i)))

(assert
 (vectors~
  (roots
   (vector 1 1 1 1 1))
  '#(0.3090169943749475-0.9510565162951535i
     0.3090169943749475+0.9510565162951535i
     -0.8090169943749473+0.5877852522924731i
     -0.8090169943749475-0.5877852522924731i)))

;; " x^6 + 4*x^5 - 72*x^4 - 214*x^3 + 1127*x^2 + 1602*x - 5040 "

(assert
 (vectors~
  (roots
   (vector 1 4 -72 -214 1127 1602 -5040))
  '#(1.9999999999999993+0.0i
     7.0+0.0i
     3.000000000000001+0.0i
     -3.0+0.0i
     -8.000000000000002+0.0i
     -4.999999999999998+0.0i)))