## Numerical Integration (With Precision)

Update: Changed from nested defn-  inside defn  to instead define local functions in let -bindings.

In a previous blog entry, I explained the higher-order function `sum` and how to use the Substitution Model to follow the execution of a function. In this entry, I will use the `sum` function to perform numerical integration, and in the process run into some pitfalls of Java’s BigDecimal. I will show how the language Clojure provides an elegant solution to the “exact quotient cannot be represented” problem of BigDecimal.

We will once again look at Section 1.3 of Structure and Interpretation of Computer Programs, Formulating Abstractions with Higher-Order Procedures. To recap, the `sum` function is an abstraction of a sum of some sequence. It takes a `term` function for calculating the current term, a start value `a`, a `next` function for calculating the next value of `a`, and an end value `b`:

```(defn sum [term a next b]
(if (> a b)
0
(+ (term a)
(sum term (next a) next b))))
```

After introducing the `sum` function, the authors describe a simple sum for calculating the integral of a function: Looking at the first term, we see that `(a+dx/2)` is the start value, ie what we will pass as the `a` parameter to `sum`. The `term` function for calculating the value of the current term is obviously `f`.

We can also see that each term has one `dx` added to it. The first term has no `dx`, the second term has one `dx`, the third term has `2dx`, etc. A simple function adding `dx` to its argument will be perfect as the `next` function. We’ll define it as a private function using `defn-`, to not clutter the namespace with internal functions. However, it needs to know the `dx`, so we’ll define it inside the `integral` function, like this:

```(defn integral [f a b dx]
(let [add-dx (fn [x] (+ x dx))]
,,,))
```

Finally, we need to multiply the result of calling `sum` by `dx`. Here is the implementation of the `integral` function:

```(defn integral [f a b dx]
(let [add-dx (fn [x] (+ x dx))]
(* (sum f (+ a (/ dx 2)) add-dx b)
dx)))```

We’ll need a function that we can pass as `f`, so we can do some interesting integration. Let’s use `cube`:

```(defn cube [x] (* x x x))
```

The integral of `x^3` is `x^4/4`. Integrating from `0` to `1` results in `(1^4)/4 = 1/4 = 0.25`. That’s what we will compare with.

Let’s test our `integral` function using a `dx` of, say `0.01`:

```user> (integral cube 0 1 0.01)
0.24998750000000042
```

Not that bad. However, in Exercise 1.29, the authors suggest using Simpson’s Rule, which is a more accurate method of numerical integration.

## Simpson’s Rule

Simpson’s Rule for numerical integration goes like this:   OK, so how do we squeeze this into the `sum` function? First we note that `h` is used in more than one place, so we’ll begin by creating a local binding `h`. We’ll also multiply the result of `sum` with `h/3`:

```(defn simpson-sum [f a b n]
(let [h (/ (- b a) n)]
,,,
(* (sum ,,, ,,, ,,, ,,,)
(/ h 3))))```

Next, we see that the factor of each term is alternating between 4 and 2, except the first and the last one, which are 1. This can be implemented as a conditional:

```(cond (zero? k) y
(= k n)   y
(odd? k)  (* 4 y)
(even? k) (* 2 y))
```

The `k` will be provided as a parameter, as we’ll see in a second, but we’re missing the `y`. So let’s add the expression of `y` and package it all up in a function inside `simpson-sum`:

```(defn simpson-sum [f a b n]
(let [h (/ (- b a) n)
simpson-term (fn [k]
(let [y (f (+ a (* k h)))]
(cond (zero? k) y
(= k n) y
(odd? k) (* 4 y)
(even? k) (* 2 y))))]
(* (sum ,,, ,,, ,,, ,,,)
(/ h 3))))```

OK, so now we have something that, given a `k`, provides the corresponding term in the sum, with the correct factor. The sequence for the sum is from `0` to `n`, so all we need to do now is to call `sum` with `simpson-term` as the `term` function, zero as the start value of the sequence, `inc` as the `next` function, and `n` as the end of the sequence:

```(defn simpson-sum [f a b n]
(let [h (/ (- b a) n)
simpson-term (fn [k]
(let [y (f (+ a (* k h)))]
(cond (zero? k) y
(= k n) y
(odd? k) (* 4 y)
(even? k) (* 2 y))))]
(* (sum simpson-term 0 inc n)
(/ h 3))))```

The test run on `cube` using `n=100` gives some impressive results:

```user> (simpson-sum cube 0 1 100)
1/4
user> (simpson-sum cube 0 2 100)
4
user> (simpson-sum cube 3 4 100)
175/4
```

Not bad at all. Clojure’s use of Ratio enables us to get accurate precision. But what happens if I select a floating point interval?

```user> (simpson-sum cube 0 0.5 100)
0.015624999999999995
```

Pretty good that too.

### Using BigDecimal

Now, what if I want to use BigDecimal in order to get really high precision in my floating points? In Clojure, the letter M on a number indicates that it’s a BigDecimal literal:

```user> (simpson-sum cube 0 0.5M 100)
Non-terminating decimal expansion; no exact representable decimal result.
[Thrown class java.lang.ArithmeticException]
Backtrace:
0: java.math.BigDecimal.divide(BigDecimal.java:1603)
1: clojure.lang.Numbers\$BigDecimalOps.divide(Numbers.java:1067)
2: clojure.lang.Numbers.divide(Numbers.java:139)
3: user\$simpson_sum.invoke(NO_SOURCE_FILE:1)
```

What happened here? Actually, this is something that can potentially happen whenever using BigDecimal divison, whether it’s from Java, Clojure or something else. The JavaDoc for `BigDecimal.divide` states:

divide(BigDecimal divisor)
Returns a BigDecimal whose value is (this / divisor), and whose preferred scale is (this.scale() – divisor.scale()); if the exact quotient cannot be represented (because it has a non-terminating decimal expansion) an ArithmeticException is thrown.

Could it be that division of `h` by 3 that messes things up? Division by 3 is always troublesome when using floating point numbers. Anyway, in order to solve this, we need to use a version of `BigDecimal.divide` that takes a scale (we’ll use 10) and a rounding mode. Only line 10 needs to be changed:

```(defn simpson-sum [f a b n]
(let [h (/ (- b a) n)
simpson-term (fn [k]
(let [y (f (+ a (* k h)))]
(cond (zero? k) y
(= k n) y
(odd? k) (* 4 y)
(even? k) (* 2 y))))]
(* (sum simpson-term 0 inc n)
(.divide (bigdec h) 3M 10 java.math.RoundingMode/HALF_UP))))```

Ugly as hell, but what can you do? (Plenty, actually, as we’ll see soon) Let’s test it:

```user> (simpson-sum cube 0 0.5M 100)
0.0156250031250000000000M
```

OK, so it works. Are we happy now? Not really. Look what happens if we use an `n` that is a multiple of 3:

```user> (simpson-sum cube 0 0.5M 6)
Non-terminating decimal expansion; no exact representable decimal result.
[Thrown class java.lang.ArithmeticException]
Backtrace:
0: java.math.BigDecimal.divide(BigDecimal.java:1603)
1: clojure.lang.Numbers\$BigDecimalOps.divide(Numbers.java:1067)
2: clojure.lang.Numbers.divide(Numbers.java:139)
3: user\$simpson_sum.invoke(NO_SOURCE_FILE:1)
```

This time it’s the other division that is the problem. We’ll change line 2 as well:

```(defn simpson-sum [f a b n]
(let [h (.divide (bigdec (- b a)) (bigdec n) 10 java.math.RoundingMode/HALF_UP)
simpson-term (fn [k]
(let [y (f (+ a (* k h)))]
(cond (zero? k) y
(= k n) y
(odd? k) (* 4 y)
(even? k) (* 2 y))))]
(* (sum simpson-term 0 inc n)
(.divide (bigdec h) 3M 10 java.math.RoundingMode/HALF_UP))))```
```user> (simpson-sum cube 0 0.5M 6)
0.0156249999937499999925000000049999999992M
```

So, we’re good now? No, Clojure has actually a much better way of handling this than resorting to the above. Let’s look at the documentation for the function with-precision:

```user> (doc with-precision)
-------------------------
clojure.core/with-precision
([precision & exprs])
Macro
Sets the precision and rounding mode to be used for BigDecimal operations.

Usage: (with-precision 10 (/ 1M 3))
or:    (with-precision 10 :rounding HALF_DOWN (/ 1M 3))

The rounding mode is one of CEILING, FLOOR, HALF_UP, HALF_DOWN,
HALF_EVEN, UP, DOWN and UNNECESSARY; it defaults to HALF_UP.
```

OK, so we simply call `with-precision` with the expression we want precision on. Pretty elegant. Here is the final implementation:

```(defn simpson-sum [f a b n]
(with-precision 10
(let [h (/ (- b a) n)
simpson-term (fn [k]
(let [y (f (+ a (* k h)))]
(cond (zero? k) y
(= k n) y
(odd? k) (* 4 y)
(even? k) (* 2 y))))]
(* (sum simpson-term 0 inc n)
(/ h 3)))))```

We’re now getting even better results than when manually providing a scale to `BigDecimal.divide`. Here is a run with `n=6`:

```user> (simpson-sum cube 0 0.5M 6)
0.01562500000M
```

Compare that with the `integral` function that we implemented first. That would have given us this result:

```user> (integral cube 0 0.5 0.01)
0.015621875000000023
```

The exact integral of `x^3` between `0` and `0.5` is `(0.5^4)/4 = 0.015625`, so it’s clear that Simpson’s Rule is very good for numerical integration. The code becomes really elegant even when using BigDecimal arithmetic, thanks to Clojure’s support for BigDecimal precision.