SICP Chapter 1 solutions in Haskell

Exercise 1.1 (not translated)

10
12
8
3
6
3
4
19
false
4
16
6
16

Exercise 1.2

Normally in Haskell you would express this using infix notation:

(5 + 4 + (2 - (3 - (6 + 4/5)))) / (3 * (6 - 2) * (2 - 7))

but it is possible to express it alternately using prefix notation:

(/) ((+) ((+) 5 4) ((-) 2 ((-) 3 ((+) 6 ((/) 4 5))))) ((*) 3 ((*) ((-) 6 2) ((-) 2 7)))

Exercise 1.3

max2 a b c = let list = [a, b, c] 
                 sqs = map (^ 2) list
                 minsq = (minimum list) ^ 2
             in sum sqs - minsq

Exercise 1.4

Evaluates to a plus the absolute value of b

Exercise 1.5

Applicative order evaluation would go into an infinite loop, trying to evaluate the (p) argument, which would keep evaluating to (p), never making any progress. Normal order does not need to fully evaluate the operands first, so it is fine with just leaving in (p). Eventually the if predicate evaluates to true, and it returns the consequent, 0, never needing to evaluate the (p), and so avoiding the infinite loop problem. (Instead of infinite loop, applicative order may actually blow up with a stack or memory overflow problem.) We would describe either of these as divergent, or ⊥ (“bottom”).

Exercise 1.6

The square root program written using “new-if” would never complete successfully, either entering an infinite loop, or blowing up due to exhausting stack space. This is because when it evaluates the first “new-if” it has to evaluate both operands since it is now not a special form. Evaluating the second operand, involves a call one level deeper, again involving a “new-if”, with the exact same requirement to evaluate both operands. We have infinite regress, and can never break out until we abort due to running out of some resource like stack space.

Exercise 1.7

With relatively large numbers, the arbitrary 0.001 absolute value tolerance cutoff works well, because the guess is much larger than 0.001, so getting the guess correct to within +/- 0.001 produces pretty good results. With very small numbers, however, 0.001 is considered large relative to the numbers, so it breaks down in terms of providing “pretty good” results. For example, even a guess of 0 would be considered within tolerance of x = 0.0005, but sqrt(0.0005) == 0 would arguably be a poor answer. A guess of 0.0007 would also be within tolerance of 0.0005, another poor answer.

Computers represent floating point numbers with a limited-precision mantissa and an exponent, similar to “scientific notation”, e.g. 2.3 x 108. Very large numbers do not have precision down to the decimal point, much less past it. So it is quite possible to get into a situation where we reach our best guess, but the difference is considered > 0.001, and due to the limited precision the next guess doesn’t put us any closer to the answer, but perhaps we oscillate between two guess values, never breaking out because it never sees a value that it thinks is within 0.001 tolerance. Alternately, it does return a result, but one that is a lousy answer. For example, the initial good-enough? has sqrt returning about 0.0323 as sqrt of 0.0001 and about 0.03126 as sqrt of 0.000001.

Initial translation of Scheme version to Haskell:

avg a b = (a + b) / 2
sqr a = a * a

improve guess x = avg guess (x / guess)
good_enough guess x = abs ((sqr guess) - x) < 0.001

sqrt_iter guess x = if good_enough guess x
                    then guess
                    else let !guess' = improve guess x
                         in sqrt_iter guess' x
sqrt = sqrt_iter 1.0

Improved version:

avg a b = (a + b) / 2.0
sqr a = a * a

improve guess x = avg guess (x / guess)
good_enough guess prevGuess = abs ((guess - prevGuess) / guess) < 0.001

sqrt x = sqrt_iter 1.0 (-10.0)
  where 
    sqrt_iter guess prevGuess = if good_enough guess prevGuess
                                then guess
                                else let !guess' = improve guess x
                                     in sqrt_iter guess' guess

This version does indeed work better on very large and small numbers.

Exercise 1.8

-- Use (x / y^2 + 2y) / 3 for successive approximation of cube root
-- Uses -10 as a dummy distant "prevGuess"
cubert x = cubertR 1.0 (-10.0)
  where cubertR guess prevGuess = if goodEnough guess prevGuess
                                  then guess
                                  else let !guess' = improve guess
                                       in cubertR guess' guess
        goodEnough guess prevGuess = abs (guess - prevGuess) / guess < 0.0001
        improve guess = ((x / guess / guess) + guess + guess) / 3.0

Exercise 1.9 (not translating these trivial examples to Clojure)

(define (+ a b)
  (if (= a 0)
    b
    (inc (+ (dec a) b))))

(+ 4 5)
(inc (+ (dec 4) 5))
(inc (+ 3 5))
(inc (inc (+ (dec 3) 5)))
(inc (inc (+ 2 5)))
(inc (inc (inc (+ (dec 2) 5))))
(inc (inc (inc (+ 1 5))))
(inc (inc (inc (inc (+ (dec 1) 5)))))
(inc (inc (inc (inc (+ 0 5)))))
(inc (inc (inc (inc 5))))
(inc (inc (inc 6)))
(inc (inc 7))
(inc 8)
9


(define (+ a b)
  (if (= a 0)
    b
    (+ (dec a) (inc b))))

(+ 4 5)
(+ (dec 4) (inc 5))
(+ 3 6)
(+ (dec 3) (inc 6))
(+ 2 7)
(+ (dec 2) (inc 7))
(+ 1 8)
(+ (dec 1) (inc 8))
(+ 0 9)
9

The first process is linear recursive. The second process is linear iterative.

Exercise 1.10

-- Compute Ackermann function version from SICP Exercise 1.10
-- Cannot name function A due to language restrictions
a x 0 = 0
a 0 y = 2 * y
a _ 1 = 2
a x y = a (x-1) (a x (y-1))

f = a 0
g = a 1
h = a 2

a 1 10 = 1024
a 2 4 = 65536
a 3 3 = 65536

For positive integer values of n,
f n computes 2n
g n computes 2n
h n computes 2^^n (tetration, 2 exponentiated by itself n times, e.g. 2^^4 =
2222

Exercise 1.11

-- compute f where f(n) = n if n < 3 and f(n-1) + 2f(n-2) + 3f(n-3) if n >= 3
-- as a recursive process and a linear process

frecur n | n < 3 = n
         | otherwise = frecur (n-1) + 2*frecur(n-2) + 3*frecur(n-3)

fiter n | n < 3 = n
        | otherwise = fiterR 3 2 1 0
  where fiterR x a b c = let !fx = a + 2*b + 3*c
                         in if x == n
                            then fx
                            else let !x' = x+1
                                 in fiterR x' fx a b

Exercise 1.12

-- Compute (row, column) element of Pascal's triangle where row & column are
-- both 1-based and assuming both are ranged correctly (i.e., no error detection here)

pascal _ 1 = 1
pascal row col | row == col = 1
               | otherwise = let !prevrow = row-1
                                 !prevcol = col-1
                             in pascal prevrow prevcol + pascal prevrow col

Exercise 1.13

Φ = (1 + √5)/2
Ψ = (1 – √5)/2

Fib(0) = 0
Fib(1) = 1
otherwise Fib(n) = Fib(n-1) + Fib(n-2)
^ all of these, by definition
Φ2 = (1 + Φ)
Ψ2 = (1 + Ψ)

Fib’(n) = (Φ^n – Ψ^n)/√5
Fib’(0) = (Φ^0 – Ψ^0)/√5 = (1 – 1)/√5 = 0/√5 = 0
Fib’(1) = (Φ^1 – Ψ^1)/√5 = √5/√5 = 1
So far, Fib’ is same as Fib. In order for this to hold true in general,
Fib’(n) = Fib’(n-1) + Fib’(n-2)
= (Φn-1 – Ψn-1 + Φn-2 – Ψn-2)/√5
= ((Φ + 1)*Φn-2 – (Ψ + 1)*Ψn-2)/√5
= ((Φ2)*Φ(n-2) – (Ψ2)*Ψ(n-2)/√5
= (Φn – Ψn)/√5

therefore Fib(n) = (Φn – Ψn)/√5 holds true in the general case.

Since this equation determines Fib(n) exactly, and |Ψn| < 0.5 for all integer n > 1
then for all n > 1, Φn/√5 differs from Fib(n) by less than 0.5.
Therefore, for all integer n > 1, closest integer to Φn/√5 must be Fib(n).
For n = 0, Φ^n/√5 is approximately .447, and for n = 1, Φn/√5 is
approximately 0.7236, so closest integers are 0 and 1, respectively, therefore for
all integer n >= 0, closest integer to Φn/√5 is equal to Fib(n).

Exercise 1.14

see Clojure version for the solution to this

Exercise 1.15

cube x = x^3
p x = 3*x - 4 * cube x
sine angle | abs angle <= 0.1 = angle
           | otherwise = p $ sine (angle / 3.0)

p is executed 5 times when executing (sine 12.15)

angle is divided by 3 each iteration until it reaches 0.1 or below, so the # of iterations is equal to about logarithm base 3 of 10*angle. The code cannot use tail call optimization, so both the order of growth of space and the number of steps is θ(log a).

Exercise 1.16

-- compute b^n in a fast iterative fashion
fast_expt :: (Fractional a, Integral b) => a -> b -> a
fast_expt b n = fast_exptR 1.0 b n
  where fast_exptR a b 0 = a
        fast_exptR a b n | even n = let !b' = b*b
                                        !n' = n `div` 2
                                    in fast_exptR a b' n'
                         | otherwise = let !a' = a*b
                                           !n' = n-1
                                       in fast_exptR a' b n'

Exercise 1.17

-- compute a*b in a fast fashion

-- double/halve are theoretically our primitives, in addition to +
-- In principle they could for example be a bit shift,
-- but we just use * and / here.
double n = n*2
halve n = n `div` 2

fast_mul a 0 = 0
fast_mul a b | even b = fast_mul (double a) (halve b)
             | otherwise = a + fast_mul a (b-1)

Exercise 1.18

-- compute a*b in a fast fashion

-- double/halve are theoretically our primitives, in addition to +
-- In principle they could for example be a bit shift,
-- but we just use * and / here.
double n = n*2
halve n = n `div` 2

fast_mul a b = fast_mulR 0 a b
fast_mulR s a 0 = s
fast_mulR s a b | even b = let !a' = double a
                               !b' = halve b
                           in fast_mulR s a' b'
                | otherwise = let !s' = s+a
                                  !b' = b-1
                              in fast_mulR s' a b'

Exercise 1.19

Starting with (a0, b0), and applying Tpq once we get (qb0 + qa0 + pa0, pb0 + qa0),
or a1 = qb0 + qa0 + pa0, b1 = pb0 + qa0. Applying Tpq a second time we get (qb1 + qa1 + pa1, pb1 + qa1),
or a2 = b0(pq + qq + pq) + a0(qq + qq + pq + pq + pp), b2 = b0(pp + qq) + a0(pq + qq + pq)

We want a2 = q’b0 + q’a0 + p’a0, b2 = p’b0 + q’a0
therefore if we take p’ = pp + qq, and q’ = 2pq + qq, it works out in one transformation.

-- compute nth Fibonacci number efficiently
fib n = fibR 1 0 0 1 n
  where fibR _ b _ _ 0 = b
        fibR a b p q count | even count = let !p' = p*p + q*q
                                              !q' = 2*p*q + q*q
                                              !count' = count `div` 2
                                          in fibR a b p' q' count'
                           | otherwise = let !a' = b*q + a*q + a*p
                                             !b' = b*p + a*q
                                             !count' = count - 1
                                         in fibR a' b' p q count'

Exercise 1.20

(define (gcd a b)
(if (= b 0)
a
(gcd b (remainder a b))))

Evaluating (gcd 206 40) in normal order:
(gcd
206
40)
(if (= 40 0)
206
(gcd
40
(remainder 206 40)))
(gcd
40
(remainder 206 40))
(if (= (remainder 206 40) 0)
40
(gcd
(remainder 206 40)
(remainder 40 (remainder 206 40))))
(if (= 6 0) ; 1 remainder performed here
40
(gcd
(remainder 206 40)
(remainder 40 (remainder 206 40))))
(gcd
(remainder 206 40)
(remainder 40 (remainder 206 40)))
(if (= (remainder 40 (remainder 206 40)) 0)
(remainder 206 40)
(gcd
(remainder 40 (remainder 206 40))
(remainder (remainder 206 40) (remainder 40 (remainder 206 40)))))
(if (= (remainder 40 6) 0) ; 1 remainder performed here
(remainder 206 40)
(gcd
(remainder 40 (remainder 206 40))
(remainder (remainder 206 40) (remainder 40 (remainder 206 40)))))
(if (= 4 0) ; 1 remainder performed here
(remainder 206 40)
(gcd
(remainder 40 (remainder 206 40))
(remainder (remainder 206 40) (remainder 40 (remainder 206 40)))))
(gcd
(remainder 40 (remainder 206 40))
(remainder (remainder 206 40) (remainder 40 (remainder 206 40))))
(if (= (remainder (remainder 206 40) (remainder 40 (remainder 206 40))) 0)
(remainder 40 (remainder 206 40))
(gcd
(remainder (remainder 206 40) (remainder 40 (remainder 206 40)))
(remainder (remainder 40 (remainder 206 40)) (remainder (remainder 206 40) (remainder 40 (remainder 206 40))))))
(if (= (remainder 6 (remainder 40 6)) 0) ; 2 remainders performed here
(remainder 40 (remainder 206 40))
(gcd
(remainder (remainder 206 40) (remainder 40 (remainder 206 40)))
(remainder (remainder 40 (remainder 206 40)) (remainder (remainder 206 40) (remainder 40 (remainder 206 40))))))
(if (= (remainder 6 4) 0) ; 1 remainder performed here
(remainder 40 (remainder 206 40))
(gcd
(remainder (remainder 206 40) (remainder 40 (remainder 206 40)))
(remainder (remainder 40 (remainder 206 40)) (remainder (remainder 206 40) (remainder 40 (remainder 206 40))))))
(if (= 2 0) ; 1 remainder performed here
(remainder 40 (remainder 206 40))
(gcd
(remainder (remainder 206 40) (remainder 40 (remainder 206 40)))
(remainder (remainder 40 (remainder 206 40)) (remainder (remainder 206 40) (remainder 40 (remainder 206 40))))))
(gcd
(remainder (remainder 206 40) (remainder 40 (remainder 206 40)))
(remainder (remainder 40 (remainder 206 40)) (remainder (remainder 206 40) (remainder 40 (remainder 206 40))))))
(if (= (remainder (remainder 40 (remainder 206 40)) (remainder (remainder 206 40) (remainder 40 (remainder 206 40)))) 0)
(remainder (remainder 206 40) (remainder 40 (remainder 206 40)))
(gcd
(remainder (remainder 40 (remainder 206 40)) (remainder (remainder 206 40) (remainder 40 (remainder 206 40))))
(remainder
(remainder (remainder 206 40) (remainder 40 (remainder 206 40)))
(remainder (remainder 40 (remainder 206 40)) (remainder (remainder 206 40) (remainder 40 (remainder 206 40)))))))
(if (= (remainder (remainder 40 6) (remainder 6 (remainder 40 6))) 0) ; 3 remainders performed here
(remainder (remainder 206 40) (remainder 40 (remainder 206 40)))
(gcd
(remainder (remainder 40 (remainder 206 40)) (remainder (remainder 206 40) (remainder 40 (remainder 206 40))))
(remainder
(remainder (remainder 206 40) (remainder 40 (remainder 206 40)))
(remainder (remainder 40 (remainder 206 40)) (remainder (remainder 206 40) (remainder 40 (remainder 206 40)))))))
(if (= (remainder 4 (remainder 6 4)) 0) ; 2 remainders performed here
(remainder (remainder 206 40) (remainder 40 (remainder 206 40)))
(gcd
(remainder (remainder 40 (remainder 206 40)) (remainder (remainder 206 40) (remainder 40 (remainder 206 40))))
(remainder
(remainder (remainder 206 40) (remainder 40 (remainder 206 40)))
(remainder (remainder 40 (remainder 206 40)) (remainder (remainder 206 40) (remainder 40 (remainder 206 40)))))))
(if (= (remainder 4 2) 0) ; 1 remainder performed here
(remainder (remainder 206 40) (remainder 40 (remainder 206 40)))
(gcd
(remainder (remainder 40 (remainder 206 40)) (remainder (remainder 206 40) (remainder 40 (remainder 206 40))))
(remainder
(remainder (remainder 206 40) (remainder 40 (remainder 206 40)))
(remainder (remainder 40 (remainder 206 40)) (remainder (remainder 206 40) (remainder 40 (remainder 206 40)))))))
(if (= 0 0) ; 1 remainder performed here
(remainder (remainder 206 40) (remainder 40 (remainder 206 40)))
(gcd
(remainder (remainder 40 (remainder 206 40)) (remainder (remainder 206 40) (remainder 40 (remainder 206 40))))
(remainder
(remainder (remainder 206 40) (remainder 40 (remainder 206 40)))
(remainder (remainder 40 (remainder 206 40)) (remainder (remainder 206 40) (remainder 40 (remainder 206 40)))))))
(remainder (remainder 206 40) (remainder 40 (remainder 206 40)))
(remainder 6 (remainder 40 6)) ; 2 remainders performed here
(remainder 6 4) ; 1 remainder performed here
2 ; 1 remainder performed here



Evaluating (gcd 206 40) in applicative order:
(gcd
206
40)
(if (= 40 0)
206
(gcd 40 (remainder 206 40)))
(gcd 40 (remainder 206 40)))
(gcd 40 6) ; 1 remainder performed here
(if (= 6 0)
40
(gcd 6 (remainder 40 6)))
(gcd 6 (remainder 40 6))
(gcd 6 4) ; 1 remainder performed here
(if (= 4 0)
6
(gcd 4 (remainder 6 4)))
(gcd 4 (remainder 6 4))
(gcd 4 2) ; 1 remainder performed here
(if (= 2 0)
4
(gcd 2 (remainder 4 2)))
(gcd 2 (remainder 4 2))
(gcd 2 0) ; 1 remainder performed here
(if (= 0 0)
2
(gcd 0 (remainder 2 0)))
2

Total of 18 remainders performed in normal order evaluation of (gcd 206 40)
Only 4 remainders performed in applicative order evaluation of (gcd 206 40)

Exercise 1.21

-- find smallest divisor of n
square n = n*n
divides a b = mod b a == 0

find_divisor n test_divisor | square test_divisor > n = n
                            | divides test_divisor n = test_divisor
                            | otherwise = let !test_divisor' = test_divisor+1
                                          in find_divisor n test_divisor'

smallest_divisor n = find_divisor n 2

Smallest divisors of 199, 1999, and 19999 are 199, 1999, and 7, respectively.

Exercise 1.23 and Exercise 1.24

I wrote a Haskell program using getClockTime, but it did not seem to show any elapsed time. I speculate that it updates its wall clock time during idle time, and if there is no idle time, it doesn’t get updated.

I switched to trying to use the GHC profiler, but found it uses a 20 ms. “tick”, far too coarse to measure timing of individual calls to prime, or even 3 calls to prime.

Nevertheless, I would expect the algorithms to approximately run twice as fast, with half the number of test steps, with maybe some difference due to constant factors.

Exercise 1.25

Alyssa’s suggestion to simply use fast-expt in conjunction with remainder would work in principle, however it is likely to be considerably slower than our expmod version. The reason for this is that we are using arbitrary precision integer math, hence operations on values with more digits take longer. Our expmod implementation applies remainder at each step, keeping values from growing very large, so the overall operation should take less time than waiting until the end to apply remainder once to a large value.

Exercise 1.26

By using explicit multiplication instead of square, each operand of the multiplication has to get evaluated independently. So whereas we only make one recursive call to expmod with the “square” implementation, with the multiplication implementation, we would actually make 2. The whole principle of halving the problem space has been nullified, and the θ(log n) becomes θ(n).

Exercise 1.27

square n = n*n
expmod _ 0 _ = 1
expmod base exp m | even exp = mod (square (expmod base (exp `div` 2) m)) m
                  | otherwise = mod (base * (expmod base (exp-1) m)) m


test_carmichael n = testCarmichaelR 1
  where testCarmichaelR a | a == n = True
                          | expmod a n n == mod a n = testCarmichaelR (a+1)
                          | otherwise = False

This does indeed return true for all the Carmichael numbers mentioned in footnote 47.

Exercise 1.28

import System.Random

-- squaring step in expmod, we check to see if we have discovered a ``nontrivial square root of 1 modulo n,''
-- that is, a number not equal to 1 or n - 1 whose square is equal to 1 modulo n. It is possible to prove that if
-- such a nontrivial square root of 1 exists, then n is not prime. It is also possible to prove that if n is an odd
square n = n*n
expmod _ 0 _ = 1
expmod base exp m | even exp = let !root = expmod base (exp `div` 2) m
                                   !sqr = mod (square root) m
                               in if root /= 1 && root /= (m-1) && sqr == 1
                                     then 0    -- found nontrivial root of 1
                                     else sqr
                  | otherwise = mod (base * (expmod base (exp-1) m)) m


miller_rabin_test :: Int -> IO Bool
miller_rabin_test n = do randint <- getStdRandom (randomR (1,(n-1)))
                         return $ try_it randint
                           where try_it a = expmod a n n == a

fast_prime :: Int -> Int -> IO Bool
fast_prime n 0 = return True
fast_prime n times = do mr_result <- miller_rabin_test n
                        if mr_result
                           then do let !times' = times-1
                                   fast_prime n times'
                           else return False

prime :: Int -> IO Bool
prime n = fast_prime n 4

Exercise 1.29

Computing integral via repeated samples over interval:

sumf term a next b | a > b = 0
                   | otherwise = let !a' = next a
                                 in term a + sumf term a' next b

integral_sample f a b dx = (sumf f (a + dx / 2.0) add_dx b) * dx
  where add_dx x = x+dx
-- SICP Exercise 1.29, Compute integral of f from a to b using n samples, via Simpson's Rule

-- Compute coefficient for term for evaluation
coefficient :: Integral a => a -> a -> a
coefficient 0 _ = 1
coefficient n total_n | n == total_n = 1
                      | even n = 2
                      | otherwise = 4

sum_simp :: (Fractional a, Integral b) => (a -> a) -> a -> a -> b -> a
sum_simp f a h n = sumSimpR 0 0
  where sumSimpR k acc | k > n = acc
                       | otherwise = let !k' = k+1
                                         !acc' = acc + (fromIntegral (coefficient k n)) * (f (a + (fromIntegral k)*h))
                                     in sumSimpR k' acc'

integral :: (Fractional a, Integral b) => (a -> a) -> a -> a -> b -> a
integral f a b n = let !h = (b-a) / fromIntegral n
                   in (sum_simp f a h n)*h / 3

Simpson’s Rule seems to converge to a solution very quickly for cube in the range 0 to 1.
At n = 2, it returns exactly 0.25. At n = 100 and n = 1000 it returns 0.25000000000000006.
The repeated sampling code doesn’t converge quite as well. At dx = 0.01 it returns
0.24998750000000042 and at dx = 0.001 it returns 0.249999875000001.

Exercise 1.30

-- Iterative function to sum terms from a to b
sumf term a next b = sumfR a 0
  where sumfR a result | a > b = result
                       | otherwise = let !a' = next a
                                         !result' = result + term a
                                     in sumfR a' result'

Exercise 1.31

-- Iterative function to multiply terms from a to b
productf :: (Integral a, Fractional b) => (a -> b) -> a -> (a -> a) -> a -> b
productf term a next b = productR a 1.0
  where productR a result | a > b = result
                          | otherwise = let !a' = next a
                                            !result' = result * term a
                                        in productR a' result'

pi_div_4 :: (Integral a, Fractional b) => a -> b
pi_div_4 nterms = productf piterm 0 (+1) nterms
  where piterm n | even n = (fromIntegral n + 2.0) / (fromIntegral n + 3.0)
                 | otherwise = (fromIntegral n + 3.0) / (fromIntegral n + 2.0)

part b, recursive version of productf:

-- Recursive function to multiply terms from a to b
productf term a next b | a > b = 1.0
                       | otherwise = let !a' = next a
                                     in term a * productf term a' next b

Exercise 1.32

Both sum and product iterate over a range of numbers, applying a certain operation to each number in turn, to generate a term, and then applying a binary merge operation to merge the current term into the accumulated result. For sum, the merge operation is + and for product, the merge operation is *. In each case, the accumulator must be initialized with the appropriate identity element for the operation. For sum, the identity element is 0, and for product, the identity element is 1. SICP calls this identity element the null-value or base value.

Part a, written in an iterative style

-- accumulate a collection of terms from a to b using a combiner function
accumulate combiner null_value term a next b = accumulateR null_value a
  where accumulateR acc n | n > b = acc
                          | otherwise = let !acc' = combiner acc (term n)
                                            !n' = next n
                                        in accumulateR acc' n'

sumf term a next b = accumulate (+) 0 term a next b
productf term a next b = accumulate (*) 1 term a next b

Part b, written in a recursive style

-- accumulate a collection of terms from a to b using a combiner function
accumulate combiner null_value term a next b | a > b = null_value
                                             | otherwise = let !a' = next a
                                                               !thisTerm = term a
                                                           in combiner (accumulate combiner null_value term a' next b) thisTerm

sumf term a next b = accumulate (+) 0 term a next b
productf term a next b = accumulate (*) 1 term a next b

Exercise 1.33

-- accumulate a filtered collection of terms from a to b using a combiner
filtered_accumulate combiner null_value term a next b filterp = filter_accR null_value a
  where filter_accR acc n | n > b = acc
                          | otherwise = let !curTerm = term n
                                            !n' = next n
                                        in if filterp curTerm
                                              then let !acc' = combiner acc curTerm
                                                   in filter_accR acc' n'
                                              else filter_accR acc n'

sum_prime_squares a b = filtered_accumulate combiner 0 id a (+1) b isPrime
  where combiner acc n = acc + n*n

-- compute product of all positive integers less than n that are relatively prime to n
prod_relative_primes n = filtered_accumulate (*) 1 id 1 (+1) (n-1) filterp
  where filterp i = (gcd i n) == 1
        gcd a 0 = a
        gcd a b = gcd b (mod a b)

-- find smallest divisor of n
square n = n*n
divides a b = mod b a == 0

find_divisor n test_divisor | square test_divisor > n = n
                            | divides test_divisor n = test_divisor
                            | otherwise = let !test_divisor' = test_divisor+1
                                          in find_divisor n test_divisor'

smallest_divisor n = find_divisor n 2
isPrime n = smallest_divisor n == n

Exercise 1.34

f g = g 2

Trying to evaluate f f fails because of Haskell’s strong typing. For f, it expects an argument of function that takes in integral argument. Instead we are passing it a function that takes a function argument, and it cannot coerce the value.

Exercise 1.35

Fixed point of x -> 1 + 1/x is the value x such that x = 1 + 1/x.
We are given Φ2 = Φ + 1 by definition of the golden ratio. Dividing
each side by Φ we get Φ = 1 + 1/Φ, and can immediately see that
Φ serves as a fixed point for x -> 1 + 1/x.

tolerance = 0.00001
fixed_point f first_guess = try first_guess
  where try guess = let !guess' = f guess
                    in if closeEnough guess guess'
                          then guess'
                          else try guess'
        closeEnough v1 v2 = abs (v1-v2) < tolerance
*Main> fixed_point (\x -> 1.0 + 1.0/x) 1.0
1.6180327868852458

Exercise 1.36

import Debug.Trace (trace)

tolerance = 0.00001
fixed_point f first_guess = try first_guess
  where try guess = let !guess' = f (trace ("guess = " ++ show guess) guess)
                    in if closeEnough guess guess'
                          then guess'
                          else try guess'
        closeEnough v1 v2 = abs (v1-v2) < tolerance

*Main> fixed_point (\x -> log 1000 / log x) 10
guess = 10.0
guess = 2.9999999999999996
guess = 6.2877098228681545
guess = 3.7570797902002955
guess = 5.218748919675316
guess = 4.1807977460633134
guess = 4.828902657081293
guess = 4.386936895811029
guess = 4.671722808746095
guess = 4.481109436117821
guess = 4.605567315585735
guess = 4.522955348093164
guess = 4.577201597629606
guess = 4.541325786357399
guess = 4.564940905198754
guess = 4.549347961475409
guess = 4.5596228442307565
guess = 4.552843114094703
guess = 4.55731263660315
guess = 4.554364381825887
guess = 4.556308401465587
guess = 4.555026226620339
guess = 4.55587174038325
guess = 4.555314115211184
guess = 4.555681847896976
guess = 4.555439330395129
guess = 4.555599264136406
guess = 4.555493789937456
guess = 4.555563347820309
guess = 4.555517475527901
guess = 4.555547727376273
guess = 4.555527776815261
guess = 4.555540933824255
4.555532257016376
*Main> 

Now, with average damping:

*Main> fixed_point (\x -> (x + (log 1000 / log x)) / 2) 10
guess = 10.0
guess = 6.5
guess = 5.095215099176933
guess = 4.668760681281611
guess = 4.57585730576714
guess = 4.559030116711325
guess = 4.55613168520593
guess = 4.555637206157649
guess = 4.55555298754564
guess = 4.555538647701617
4.555536206185039
*Main> 

It only took 11 steps with average damping, versus 34 steps without damping.

Exercise 1.37

Recursive version

-- Evaluate continued fraction based upon n and d which return terms for
-- the numerator and denominator, and k is number of terms (indexed by 1).
cont_frac n d k = contR 1
  where contR i | i == k = (n i) / (d i)
                | otherwise = (n i) / ((d i) + contR (i+1))

It takes k of 10 to approximate golden ratio to 4 decimal digits:

*Main> cont_frac (const 1.0) (const 1.0) 10
0.6179775280898876

Iterative version:

-- Evaluate continued fraction based upon n and d which return terms for
-- the numerator and denominator, and k is number of terms (indexed by 1).
-- this version uses an iterative process
cont_frac n d k = cont_fracR k 0
  where cont_fracR 0 acc = acc
        cont_fracR i acc = let !i' = i-1
                               !acc' = (n i) / ((d i) + acc)
                           in cont_fracR i' acc'

Exercise 1.38

-- compute e via Euler's continued fraction, using k steps
-- all Ni = 1, Di = 1,2,1,1,4,1,1,6,1,1,8,...
compute_e :: (Integral a, Fractional b) => a -> b
compute_e k = 2.0 + cont_frac (const 1.0) d k
  where d i | i `mod` 3 == 2 = fromIntegral (i + i + 2) / 3.0
            | otherwise = 1.0

-- Evaluate continued fraction based upon n and d which return terms for
-- the numerator and denominator, and k is number of terms (indexed by 1).
-- this version uses an iterative process
cont_frac :: (Integral a, Fractional b) => (a -> b) -> (a -> b) -> a -> b
cont_frac n d k = cont_fracR k 0.0
  where cont_fracR 0 acc = acc
        cont_fracR i acc = let !i' = i-1
                               !acc' = (n i) / ((d i) + acc)
                           in cont_fracR i' acc'

Exercise 1.39

-- compute tangent function via Lambert's continued fraction
-- N1 = r; all other Ni = -r^2;  Di = 2i-1;  x is angle in radians; k is # steps
tan_cf :: (Integral a, Fractional b) => b -> a -> b
tan_cf x k = let mxx = -x * x
                 n 1 = x
                 n _ = mxx
             in cont_frac n (\i -> 2.0 * fromIntegral i - 1.0) k

-- Evaluate continued fraction based upon n and d which return terms for
-- the numerator and denominator, and k is number of terms (indexed by 1).
-- this version uses an iterative process
cont_frac :: (Integral a, Fractional b) => (a -> b) -> (a -> b) -> a -> b
cont_frac n d k = cont_fracR k 0.0
  where cont_fracR 0 acc = acc
        cont_fracR i acc = let !i' = i-1
                               !acc' = (n i) / ((d i) + acc)
                           in cont_fracR i' acc'

Exercise 1.40

-- return function that computes f(x) -> x^3 + ax^2 + bx + c"
cubic a b c x = x^3 + a*x^2 + b*x + c

Exercise 1.41

-- takes 1-arg procedure and returns procedure that applies original twice
double f = f . f

((double (double double)) (+ 1)) 5 equals 21

Exercise 1.42

-- compose f after g: create function that takes x and returns f(g(x))
compose = (.)

Exercise 1.43

-- nth repeated application of f
repeated f n = (\x -> if n == 0
                         then x
                         else (f . (repeated f (n-1))) x)

Exercise 1.44

-- nth repeated application of f
repeated f n = (\x -> if n == 0
                         then x
                         else (f . (repeated f (n-1))) x)

-- create smoothed version that averages f(x) with f(x-dx) and f(x+dx)
dx = 0.0001
smooth f x = (f (x-dx) + f x + f (x+dx)) / 3

-- create the n-fold smoothed function based upon f
smooth_repeated f n = (repeated smooth n) f

Exercise 1.45

Experimentation when I was working the exercise in Clojure with fixed-point instrumented with debugging println empirically shows that the number of applications of average damping required for the nth root calculation increases by 1 each time n reaches another power of 2. For n = 2, 1 application required, but for n = 32, 5 applications are required. For n = 31, 4 applications required. So a should be equal to the floor of log-base-2 of n.

This leads to the following code in Haskell:

-- return a modified version of function that averages its value with its arg
average_damp :: Fractional a => (a -> a) -> a -> a
average_damp f x = (x + f x) / 2.0

tolerance = 0.00001
fixed_point f first_guess = try first_guess
  where try guess = let !guess' = f guess
                    in if closeEnough guess guess'
                          then guess'
                          else try guess'
        closeEnough v1 v2 = abs (v1-v2) < tolerance

-- nth repeated application of f
repeated :: (Integral b) => (a -> a) -> b -> (a -> a)
repeated f n = (\x -> if n == 0
                         then x
                         else (f . (repeated f (n-1))) x)

-- return floor of log base 2 of n for positive n
floor_log2 n = floor_log2R 0 2
  where floor_log2R l ll | ll > n = n
                         | otherwise = let !l' = l+1
                                           !ll' = 2*ll
                                       in floor_log2R l' ll'

-- find nth root of x using repeated applications of average damping
nthroot n x = fixed_point ((repeated average_damp (floor_log2 n)) (\y -> x / y^(n-1))) 2.0

Exercise 1.46

tolerance=0.00001

-- perform iterative improvements of guess until it's good enough
iterative_improve goodEnough improve initialGuess = iterR initialGuess (improve initialGuess)
  where iterR prevGuess newGuess | goodEnough newGuess prevGuess = newGuess
                                 | otherwise = let !newGuess' = improve newGuess
                                               in iterR newGuess newGuess'

-- return true when successive guesses are within one hundredth of one percent of each other
goodEnough guess prevGuess = abs ((guess - prevGuess) / guess) < 0.0001

-- calculate square root of x via iterative improvement
sqrt x = (iterative_improve goodEnough improveSqrt) 1.0
  where improveSqrt guess = (guess + x/guess) / 2.0

-- implement fixed-point using iterative improvement
fixed_point f firstGuess = (iterative_improve goodEnough (\x -> f x)) firstGuess
About these ads

About alexstangl

Software developer
This entry was posted in Uncategorized and tagged , . Bookmark the permalink.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s