Skip to content
August 30, 2010 / Alex Miller

Newton-Raphson Square Roots with Clojure

Guest post from Club member Robert Scanlon:

In last week’s cljub session, a question was posed as to how to start thinking in terms of functional programming, coming from an object-oriented programming background.  This question, and the topic of “FP uptake”, is of interest to me, not only because I’m trying to figure this out myself, but also because our company has gone down the path of using FP for our core product development (and most developers are still coming from an OO background).

At the time, I was in the middle of reading John Hughes’ paper from 1984, Why Functional Programming Matters, and thought of suggesting that as one (of many) ways to “grok FP”.  But the paper has some funny syntax in it (relative to my 1 and only FP data point, clojure).  So I decided to try to implement some of the examples in the paper using clojure.  Below (couldn’t figure out how to attach) is a REPL session I did to implement the first example, a Newton-Raphson square root algorithm, from section 4.1.  I annotated the session so people new to clojure might see how someone might incrementally work up a function to solve a problem in clojure using the repl.

This implementation uses a few interesting functions and capabilities in clojure: iterate, partial, cons, recur, destructuring in argument lists, variable-length arg lists, passing functions as arguments.

The paper can be found on the web in various forms.  One is here: http://www.cse.chalmers.se/~rjmh/Papers/whyfp.pdf.  There’s also an interview with Hughes on InfoQ: http://www.infoq.com/interviews/john-hughes-fp.  [BTW, in that interview, Hughes says his favorite book is “Learning Functional Programming” by Paul Hudak; that book uses Haskell though].

If someone is new to clojure and wants to learn a lot fast, you may be interested in doing something similar with the examples in sections 4.2 (numerical differentiation) and 4.3 (numerical integration).  This was the first time I explicitly used lazy evaluation to actually solve a problem, and was also the first time I did anything with destructuring, so this was well worth my time.

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Why Functional Programming Matters, John Hughes, 1984.
; Section 4.1, Newton-Raphson Square Roots
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Iteration function, version 1
(defn fx [n x] (/ (+ x (/ n x)) 2))
(fx 100 10)
; 10
(fx 100 11)
; 221/22
(fx 100 12)
; 61/6
; Need to call this function repeatedly, passing in the latest estimate x each time.
; The paper uses the function 'repeat', let's start there.
(doc repeat)
; Nope, that just repeats the value or sequence you pass in; no functions involved.
(doc repeatedly)
; Nope, this calls the same no-arg function repeatedly. Not much use for this situation.
(find-doc "repeat")
; Lots of spammage, nothing useful here...
; Finally ran across 'iterate' function somewhere on internet while looking for
; something totally different...
(doc iterate)
; -------------------------
; clojure.core/iterate
; ([f x])
;  Returns a lazy sequence of x, (f x), (f (f x)) etc. f must be free of side-effects
; Note that the iterate function takes 1 arg, so 'bind' in the 1st arg to the fx function
; using the clojure 'partial' function.
(iterate (partial fx 100) 15)
; Ruh roh, just killed by repl with an infinite sequence! (But it seemed like it worked)
(take 5 (iterate (partial fx 100) 15))
; (15 65/6 1565/156 976565/97656 381469726565/38146972656)
; Don't really like the ratios, so convert to float instead - iteration fn, version 2
(defn fx2 [n x] (float (/ (+ x (/ n x)) 2)) )
(take 5 (iterate (partial fx2 100) 15))
; (15 10.833333 10.032051 10.0000515 10.0)
; First cut at the function to check for convergence.  The paper's definition looks
; similar to clojure's destructuring, so try that.
(defn within [ [a b remainder] ]
(if (<= (abs (- a b)) 0.0001) b (recur (cons b remainder))) )
; java.lang.Exception: Unable to resolve symbol: abs in this context
(doc repeat)
; Nuthin.
; I thought clojure had an 'abs', maybe that's incanter.  Anyways, it's easy to write one.
(defn abs [x] (if (< x 0) (- 0 x) x))
(abs 3)
; 3
(abs 0)
; 0
(abs -1)
; 1
; Now back to our regularly scheduled programming...
(defn within [ [a b remainder] ]
(if (<= (abs (- a b)) 0.0001) b (recur (cons b remainder))) )
(within [10 11 12 13 14])
; java.lang.IllegalArgumentException: Don't know how to create ISeq from: java.lang.Integer
; Gaaa. Something's amiss, probably my destructuring attempt.
; Make sure I'm calling cons the right way...
(doc cons)
; cons needs a seq for 2nd arg
; Digression to check destructuring...
(defn test-destruct-args [[a b c]] (println a (type a) b (type b) c
(type c)))
(test-destruct-args [1 2 3 4 5])
; 1 java.lang.Integer 2 java.lang.Integer 3 java.lang.Integer
; nil
; Last arg is being seen as a scalar, need it to be the 'rest of the sequence'
(defn test-destruct-args [[a b & c]] (println a (type a) b (type b) c
(type c)))
(test-destruct-args [1 2 3 4 5])
; 1 java.lang.Integer 2 java.lang.Integer (3 4 5)
clojure.lang.PersistentVector$ChunkedSeq
; nil
; Need to put an '&' before last arg to get the 'rest of the sequence'
(defn within [ [a b & remainder] ]
(if (<= (abs (- a b)) 0.0001) b (recur (cons b remainder))) )
(within [10 11 12 13 14])
; java.lang.NullPointerException
; Not handling the 'end of sequence' condition, (- a b) blows up with nils.
; [Also, not handling the termination condition for the recur, so we'd go on forever!]
(defn within [ eps [a b & remainder] ]
(if (<= (abs (- a b)) eps)
b                                                              ; if true part
(if-not (nil? remainder) (recur eps (cons b remainder)) nil)  ; else part
)
)
; Oh, btw, added in eps as an argument above, removed the hard-coded value
(within 0.0001 '(12 13 14 15 15.1 15.11 15.111 15.1111 15.11111
15.111111))
; 15.1111
(within 0.0001 (iterate (partial fx2 100) 15))
; 10.0
; Now put together the iteration function fx (or fx2) and the convergence function, within
(within 0.0001 (iterate (partial fx2 100) 900))
; 10.0
(within 0.0001 (iterate (partial fx2 100) 0))
; java.lang.RuntimeException: java.lang.ArithmeticException: Divide by zero
; ignore trapping edge cases
(within 0.0001 (iterate (partial fx2 100) 1))
; 10.0
(within 0.000000001 (iterate (partial fx2 100) 1))
; 10.0
(within 0.000000001 (iterate (partial fx2 877) 1))
; 29.614185
(within 0.000000001 (iterate (partial fx 877) 33))
;
725383613681596145326847250447451516047505115519/24494464201290276417946743597698555717522343369
(float (/ 725383613681596145326847250447451516047505115519
24494464201290276417946743597698555717522343369))
; 29.614185
(float
725383613681596145326847250447451516047505115519/24494464201290276417946743597698555717522343369)
; 29.614185
; Define the sqrt function now that we have the pieces...
(defn sqrt [eps a0 n]
(within eps (iterate (partial fx n) a0)))
(sqrt 0.0000001 14 100)
; big long ratio, yuck
; I could redefine sqrt with my other iteration function, fx2.  But why not parameterize
; the sqrt function with the iteration function?
(defn sqrt [iter-fn eps a0 n]
(within eps (iterate (partial iter-fn n) a0)))
(sqrt fx 0.0000001 14 100)
; 15917322219892801768783874/1591732221989280176878387
(float 15917322219892801768783874/1591732221989280176878387)
; 10.0
(sqrt fx2 0.0000001 14 100)
; 10.0
; But of course it's gratuitous overkill (and confusing to users) to always specify that
; function.  So customize for a particular usage (application, enterprise library, etc).
(def my-sqrt (partial sqrt fx 0.0000001))
(my-sqrt 14 100)
; 15917322219892801768783874/1591732221989280176878387
(def my-sqrt (partial sqrt fx2 0.0000001))
(my-sqrt 14 100)
; 10.0
; So, might as well just go for the gold now and parameterize the convergence fn:
(defn sqrt [iter-fn conv-fn eps a0 n]
(conv-fn eps (iterate (partial iter-fn n) a0)))
(def my-sqrt (partial sqrt fx2 within 0.0000001))
(my-sqrt 14 100)
; 10.0
(my-sqrt 14 877)
; 29.614185
; It is easy now to fold in the 'relative' version of the convergence fn from the paper.
;;;;;;;;;;;;;;;;
; Post-mortem
;;;;;;;;;;;;;;;;
; Lessons:
; - get the pieces working first, then put them together
; - don't worry about non-essentials (like the tolerance above), until you get the essentials
; - break out and experiment 'in general' when necessary; understand how something works
; Why this particular ordering of the parameters in the sqrt argument list?
; - iteration & convergence functions 1st, to partial out for particular application
; - eps second, to partial out for particular use in a fn or module (for known 'size' #'s)
; - initial value third, possibility to partial out for particular function or module
; - number to be square-rooted last - always different
; If you were really doing what I did above, parameterizing the sqrt function with the
; various other functions, the generic sqrt would likely be named something else, like
; 'sqrt-gen', and the specific incantations named something more user-friendly, like
; 'sqrt' (if only 1 used primarilly) or 'sqrt1', 'sqrt2' or 'sqrta', 'sqrtb', etc.
Advertisements
August 21, 2010 / Alex Miller

Announcing the Clojure Lunch Cljub!

The Clojure Lunch Cljub is a monthly lunch discussion group to talk about all things Clojure.  It will be informal (or at least less formal than a user group) and consist of whoever is interested reading, writing, and talking Clojure.  Hopefully we can do some code review, look at new libraries, hack on some code problems, or just get advice from your local Clojure clan.

The club will meet at 1099 Milwaukee, Kirkwood, MO 63122 in the Boiler Room, located in the NE corner of the lower floor.  Tentatively, the meetings will be on the 4th Thursday of the month starting August 26th at 11:30 am – 12:30 pm.

A google group has been created for discussion of topics and other stuff.  The first meeting will tentatively discuss the just released Clojure 1.2 features (records, protocols, etc).

All levels of Clojure experience are welcome!  If you’re just getting started, we’re happy to help.