(define $cubic-formula c-f)

(define $c-f
  (lambda [$f $x]
    (match (coefficients f x) (list math-expr)
      {[<cons $a_0 <cons $a_1 <cons $a_2 <cons $a_3 <nil>>>>>
        (c-f' a_3 a_2 a_1 a_0)]})))

(define $c-f'
  (lambda [$a $b $c $d]
    (match [a b c d] [math-expr math-expr math-expr math-expr]
      {[[,1 ,0 $p $q]
        (let {[[$s1 $s2] (2#[(rt 3 %1) (rt 3 %2)] (q-f' 1 (* 27 q) (* -27 p^3)))]}
          [(/ (+ s1 s2) 3)               ; r1
           (/ (+ (* w^2 s1) (* w s2)) 3) ; r2
           (/ (+ (* w s1) (* w^2 s2)) 3) ; r3
           ])]
       [[,1 _ _ _]
        (3#[(- %1 (/ b 3)) (- %2 (/ b 3)) (- %3 (/ b 3))]
           (with-symbols {x y}
             (c-f (substitute {[x (- y (/ b 3))]} (+ x^3 (* b x^2) (* c x) d)) y)))]
       [[_ _ _ _] (c-f' 1 (/ b a) (/ c a) (/ d a))]})))

(define $w (/ (+ -1 (* i (sqrt 3))) 2))

(* (- x 1) (- x 2) (- x 3))
;=>(+ x^3 (* -6 x^2) (* 11 x) -6)

(c-f (+ x^3 (* -6 x^2) (* 11 x) -6) x)
;=>[2 1 3]

(3#%1 (c-f (+ x^3 (* p x) q) x))
;=>
;(/ (+ (rt 3 (+ (* -108 q)
;               (* 12 (sqrt (+ (* 81 q^2) (* 12 p^3))))))
;      (rt 3 (+ (* -108 q)
;               (* -12 (sqrt (+ (* 81 q^2) (* 12 p^3)))))))
;   6)

(3#%1 (c-f (+ (* a x^3) (* b x^2) (* c x) d) x))
;=>
;(/ (+ (* (rt 3 (/ (+ (* -108 d a^3) (* 36 c b a^2) (* -8 b^3 a) (* 12 (sqrt (+ (* 81 a^6 d^2) (* -54 a^5 d c b) (* 12 a^4 d b^3) (* -3 a^4 c^2 b^2) (* 12 a^5 c^3))))) a^4)) a)
;      (* (rt 3 (/ (+ (* -108 d a^3) (* 36 c b a^2) (* -8 b^3 a) (* -12 (sqrt (+ (* 81 a^6 d^2) (* -54 a^5 d c b) (* 12 a^4 d b^3) (* -3 a^4 c^2 b^2) (* 12 a^5 c^3))))) a^4)) a)
;      (* -2 b))
;    (* 6 a))