385 - PurelyFunctional.tv Newsletter https://purelyfunctional.tv/issues/purelyfunctional-tv-newsletter-385-two-map-patterns-entity-vs-index/

00 quadratic formula.md

roots of a quadratic equation

Quadratic equations can be represented by three numbers, a, b, and c, which are the coefficient of x^2, the coefficient of x, and the constant term. The roots of a quadratic equation are everywhere where it touches the x axis, meaning the equation is equal to zero.

You can use the quadratic formula which calculates the roots. In fact, that's your task: write a function that returns the roots of a quadratic equation using the quadratic formula. Here is more information about it.

Note: you don't have to return complex roots if the curve does not cross the x-axis.

Thanks to this site for the challenge idea where it is considered Medium level in Python.

Email submissions to eric@purelyfunctional.tv before July 12, 2020. You can discuss the submissions in the comments below.

Bobby Towers.clj

(defn quadratic-roots [[a b c]]
  (let [m (* -1/2 (/ b a))
        d (Math/sqrt (- (* m m) (/ c a)))]
    [(+ m d) (- m d)]))

Eric Normand.clj

(defn roots [[a b c]]
  (let [discr (- (* b b) (* 4 a c))]
    (if (neg? discr)
      #{}
      (let [sd    (/ (Math/sqrt discr) 2 a)
            bpart (/ (- b) 2 a)]
        (conj #{}
              (+ bpart sd)
              (- bpart sd))))))

Mark Champine.clj

;; per Po-Shen Loh - https://www.poshenloh.com/quadraticdetail/
;; ax^2 +bx +c
;; divide off a to make polynomial monic
;; then roots are b/2a + u and b/2a - u
;; b/2a^2 - u^2 = c/a
;; therefore u = sqrt(b/2a^2 - c/a)

(defn qsolv [a b c]
  (let [hb (/ b a -2)
        u (Math/sqrt (- (* hb hb) (/ c a)))]
    (distinct [(+ hb u) (- hb u)])))

;; some tests

[(= [-1.5857864376269049 -4.414213562373095]
    (qsolv 1 6 7))
 (= [-0.19999999999999996 -1.0]
    (qsolv 5 6 1))
 (= [1.0 -3.0]
    (qsolv 1 2 -3))
 (= [3.0 0.5]
    (qsolv 2 -7 3))
 (= [13.348469228349535 -1.3484692283495345]
    (qsolv 1 -12 -18))
 (= [2.0]
    (qsolv 1 -4 4))
 (= [0.0]
    (qsolv 1 0 0))
 (= [4.0 3.0]
    (qsolv 1 -7 12))]

ninjure.clj

(defn quadratic-equation [a b c]
  (let [discriminant (- (* b b) (* 4 a c))]
    (when-not (neg? discriminant)
      (distinct (map #(/ (% (- b) (Math/sqrt discriminant))
                        2 a)
                     [+ -])))))

Steffan Westcott.clj

(defn quad [a b c]
  (let [disc (- (* b b) (* 4 a c))]
    (when (and (>= disc 0) (not= a 0))
      (let [w (Math/sqrt disc)]
        (set (map #(/ (- % b) (* 2 a)) [w (- w)]))))))

Teodor Heggelund.clj

(ns th.scratch.quadratic-formula)

(defn real-solutions [a b c]
  (let [discriminant (- (* b b)
                        (* 4 a c))]
    (cond
      (pos? discriminant) [(/ (+ (- b) discriminant)
                              (* 2 a))
                           (/ (- (- b) discriminant)
                              (* 2 a))]
      (zero? discriminant) [(/ (- b)
                               (* 2 a))]
      (neg? discriminant) [])))

(comment
  (for [[a b c] [[1 -4 3] [4 -8 3] [1 4 4] [1 2 4]]]
    (real-solutions a b c))
  ;; => ([4 0] [3 -1] [-2] [])
  )

(defn ^:dynamic *make-complex* [re im]
  {:re re
   :im im})

(defn real->complex [r]
  (*make-complex* r 0))

(defn complex-solutions [a b c]
  (let [discriminant (- (* b b)
                        (* 4 a c))]
    (cond
      (pos? discriminant)
      (map real->complex
           [(/ (+ (- b) discriminant)
               (* 2 a))
            (/ (- (- b) discriminant)
               (* 2 a))])
      (zero? discriminant)
      (map real->complex
           [(/ (- b)
               (* 2 a))
            (/ (- b)
               (* 2 a))])
      (neg? discriminant)
      [(*make-complex* (/ (- b)
                          (* 2 a))
                       (/ (Math/sqrt (- discriminant))
                          (* 2 a)))
       (*make-complex* (/ (- b)
                          (* 2 a))
                       (- (/ (Math/sqrt (- discriminant))
                             (* 2 a))))])))

(comment
  (for [[a b c] [[1 -4 3] [4 -8 3] [1 4 4] [1 2 4]]]
    (complex-solutions a b c))
  ;; => (({:re 4, :im 0} {:re 0, :im 0})
  ;;     ({:re 3, :im 0} {:re -1, :im 0})
  ;;     ({:re -2, :im 0} {:re -2, :im 0})
  ;;     [{:re -1, :im 1.7320508075688772} {:re -1, :im -1.7320508075688772}])
  )

(defn quadratic-equation [a b c]
  (let [discriminant (- (* b b) (* 4 a c))]
    (when-not (neg? discriminant)
      (distinct (map #(/ (% (- b) (Math/sqrt discriminant))
                        2 a)
                     [+ -])))))

Is anyone also interested in implementing the full functionality, i.e. for negative discriminants as well? The reference notes show that the complex roots can also be neatly expressed using real coordinates.

I suggest something like below to make and detect complex numbers - the idea is to be able to test all implementations uniformly. So for example a test could simply use the functions below to detect and extract values in order to verify results.

So a solution could simply reference this externally (to reduce clutter), and use it to make complex numbers accordingly. So for example, a result value could be any one of [m], [m, n] or [(complex a b), (complex c d)] where m, n are regular numbers - and so forth.

Any comment? Thanks.

(defprotocol Complex 
  (real [_])
  (ifact [_]))

(defn complex [r i]
  (if (not (zero? i))
    (reify Complex 
      (real [_]  r)
      (ifact [_] i))
    (ex-info "Complex number imaginary component must not be zero"
             {:real r, :i-factor i})))

(defn complex? [x]
  (extends? Complex (type  x)))

@KingCode I expanded the idea a bit and implemented basic arithmetic functions. Using a tuple felt more ergonomic than a protocol imo: here's the code:

(ns complex
  (:refer-clojure :rename {+ core+
                           - core-
                           / core-div
                           * core*}))

(defn complex [r i]
  [r i])

(defn complex? [x]
  (= (count x) 2))

(defn *
  ([v] v)
  ([[a1 b1] [a2 b2]]
   [(core- (core* a1 a2) (core* b1 b2)) (core+ (core* a1 b2) (core* a2 b1))]))

(defn + [v1 v2]
  (mapv core+ v1 v2))

(defn - [v1 v2]
  (mapv core- v1 v2))

(defn scale [s v]
  (mapv (partial core* s) v))

(defn divide-by-scalar [s v]
  (mapv #(core-div % s) v))

(defn / [v1 [a2 b2]]
  (divide-by-scalar (core+ (core* a2 a2) (core* b2 b2)) (* v1 [a2 (core- b2)])))

(comment
  (+ [1M -3M] [2M 5M])                   ;=> (3M 2M)
  (* [1M 1M] [1M 1M])                    ;=> [0M 2M]
  (scale 5M [1M 3M])                     ;=> (5M 15M)
  (/ [3M 3M] [3M 3M])                    ;=> (1M 0M)
  )

I also got really carried away and implemented a numerical approximation algorithm called Durand-Kerner Method:

(ns poly
  "Solve polynomials using Durand Kerner method"
  (:require complex))
;; helper functions to initialize algorithm

(defn initialize-xs
  [coefs & {:keys [base] :or {base [0.4M 0.9M]}}]
  (loop [vi [1M 0M]
         vs []]
    (if (>= (count vs) (dec (count coefs)))
      vs
      (recur (complex/* vi base) (conj vs vi)))))

;; Polynomial function
(defn eval-polynomial [coefs x]
  (loop [[coef & tail] (reverse coefs)
         x-pow         [1M 0M]
         result        [0M 0M]]
    (if coef
      (recur tail (complex/* x-pow x) (complex/+ result (complex/scale coef x-pow)))
      result)))

;; Durand-Kerner Method implementation
(defn- step [f xs i]
  (update xs i
          #(complex/- % (complex// (f %)
                       (transduce (map-indexed (fn [j k] (if (= j i) [1M 0M] (complex/- % k))))
                                  complex/* [1M 0M] xs)))))

(defn find-roots
  [coefs & {:keys [max-iterations precision] :or {max-iterations 30
                                                  precision      5}}]
  (let [degree (dec (count coefs))
        f      (partial eval-polynomial (map #(/ % (first coefs)) coefs))]
    (loop [xs (initialize-xs coefs)
           i  0]
      (let [xs' (->> (range degree)
                     (reduce (partial step f) xs))]
        (if (or (with-precision precision (= xs xs'))
                (>= i max-iterations))
          xs'
          (recur xs' (inc i)))))))

(comment
  (initialize-xs [1M -3M])

  (eval-polynomial [3M 1M] [0M 1M])

  (with-precision 5 (find-roots [5M 6M 1M]))       ;=> [[-1.0000M 0E-20M] [-0.20000M 0E-23M]]
  (with-precision 5 (find-roots [1M -3M 3M -5M])) ; => [[2.5873M -1.2E-109M] [0.20630M 1.3747M] [0.20630M -1.3747M]]
  )

I wrote about this a little bit over here:
https://porkostomus.gitlab.io/posts-output/2019-05-15-quadraticformula/

I wrote about this a little bit over here:
https://porkostomus.gitlab.io/posts-output/2019-05-15-quadraticformula/

Cool, and I love the use of Klipse! It would be great to see more explanatory text though! What does quadratic-rational do, and how does "factor and remove perfect squares" help get you there. And what's the connection between those those and your quadratic-roots code? I'm also enjoying your other posts, e..g. the sudoku solver and minesweeper in reagent are really neat.

Btw, related to this challenge, I wrote a function to to extract coefficients from a symbolic representation of a quadratic:
(->coeffs "45x^2-42x+3") => [45 -42 3] ..but it was surprisingly hard, and I'm not thrilled with it. Seems like this is something you might have already solved in a more elegant way..

Thanks for the feedback! Yes the blog post reads more like some half-baked notes, which is why I was delighted when I saw Eric was doing this because I'd gotten a bit stuck and this is the perfect chance to revisit this and document it better.

Elegant... not so sure, but my struggles in parsing math expressions led me to doing something with instaparse and core.match: https://nextjournal.com/btowers/algae

Libraries to look at:

• https://github.com/noprompt/meander
• https://github.com/clojure-numerics/expresso
• https://github.com/Sophia-Gold/symbolic-algebra

Also in this thread some kind folks helped me out:
https://clojureverse.org/t/working-with-nested-data-for-algebraic-solver/4359

Thanks for the pointers! I assumed a complete regex or instaparse would be the way to go, but would take more time investment than my one-off ad-hoc solution, so was wondering if either of those approaches already existed for polynomial expressions.
I was able to use your algae grammar with some small changes to extract coefficients. These two lines changed in parse-expr:

   <factor> = ('-'? number) | var | ratio | parens | ('-'? var)
   var = number? #'[A-Za-z^2]+'

Much appreciate the Po-Shen Loh solution! I learned it from Grant Sanderson who explains it very well in this 3blue1brown video: https://www.youtube.com/watch?v=MHXO86wKeDY&t=2076s

Regarding the different approaches to the classic formula that Eric wrote about in the follow-up, I think that each refactoring is important and should be used to expand our understanding because it illuminates the problem space from another angle. I want to try to find even more ways to do it!

Inspired by this (as well as Eric's markdown editor), I built a little reagent app for rendering LaTeX formulas from vector representations of polynomial coefficients (eg. [2 7 6) -> "2x^2 + 7x + 6"):

https://github.com/porkostomus/mecca-math