Function composition via folding

Let us suppose that we are given a list of functions. We would like to compose the functions in this list (if they are composable, of course) in the order in which they are arranged in the list, i.e., we would like to create a higher-order function compose such that

Here is one of many possible implementations in Haskell:

compose :: [(a -> a)] -> (a -> a)
compose = foldr (.) (\x -> x)

where we use a right-fold to compute the composition. Note that the identity function is defined anonymously. What happens if we apply compose to the empty list? We should obtain the identity function. Let us see if this is the case:

compose [] = foldr (.) (\x -> x) [] = (\x -> x)

Indeed, it is the case: if we apply compose to the empty list, we obtain the identity function. What is the type of this identity function? We load the two-line script above into GHCi and experiment a bit:

*Main> let id = compose []
*Main> :t id
id :: a -> a
*Main> map id [0..9]
[0,1,2,3,4,5,6,7,8,9]
*Main> map id ['a'..'z']
"abcdefghijklmnopqrstuvwxyz"

We thus obtain a polymorphic identity function that merely outputs its input, regardless of the type of the input. So far, so good. What happens if we apply compose to a non-empty list? Here is some more equational reasoning:

compose [f0,f1,f2] = foldr (.) (\x -> x) [f0,f1,f2]
                   = foldr (.) (\x -> x) (f0 : f1 : f2 : [])
                   = f0 . (f1 . (f2 . (\x -> x)))

A couple of days ago, we saw how to compose a function with itself times. This is a special case of composing a list of functions, of course. Using function replicate, we can compose a function f with itself  times in the following manner:

compose (replicate n f) = foldr (.) (\x -> x) (replicate n f)

as suggested by Reid Barton. Let us now repeat the experiment we did perform a couple of days ago, but using the new compose:

*Main> let f =  (\n -> (n+1) `mod` 64)
*Main> let h = compose (replicate 64 f)
*Main> map h [0..63]
[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,
 16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,
 32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,
 48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63]

In the next posts, we will do some interesting things with this compose.

Composing an endofunction with itself

Suppose we are given an endofunction . Let us introduce a new function , where , which we define recursively as follows

where is the identity function on , i.e., for all . In Haskell, it is easy to create a higher-order function compose that composes an endofunction with itself times:

compose :: (a -> a) -> Integer -> (a -> a)
compose f 0 = (\x -> x)
compose f k = f . (compose f (k-1))

where the identity function is defined anonymously. If you happen to dislike anonymity in your functions, here is another implementation:

compose :: (a -> a) -> Integer -> (a -> a)
compose f 0 = identity
compose f k = f . (compose f (k-1))

identity :: a -> a
identity x = x

There are many other alternatives. Take a look at this Stack Overflow thread. I particularly liked Reid Barton‘s elegant solution using a right-fold:

compose f n = foldr (.) id (replicate n f)

where id = (\x -> x). Beautiful, isn’t it? If this one-liner does not convert you to Haskell, then nothing will.

__________

Example

Let , and let function be defined by

This function corresponds to a cyclic shift. Therefore, we expect to have . Let and . We load the script above into GHCi, and voilà, we have the session:

*Main> let ns = [0..63]
*Main> let f =  (\n -> (n+1) `mod` 64)
*Main> map f ns
[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,
 17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,
 33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,
 49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,0]

So far, so good! Let us take a look at and :

*Main> let g = compose f 32
*Main> let h = compose f 64
*Main> map g ns
[32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,
 48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,
 0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,
 16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31]
*Main> map h ns
[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,
 16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,
 32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,
 48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63]

Indeed, as we expected, we have .

__________

Everything in this post is elementary. But now that the elements are solid, we can dare to attempt more sophisticated things. Like what? Like composing lists of (distinct) functions. I am thinking of permutations. I am thinking of perfect shuffles.

Generalizing Heron’s method

Suppose we are given a positive rational number . We would like to compute an approximation of the -th root of , which we denote by . If , we can use Heron’s method to approximate , as we saw a few days ago. But, what if we have that ? Is it possible to generalize Heron’s method?

Indeed, it is. Let us introduce a function , defined by , where and . By construction, the positive zero of function is . Do note that if is odd, then is not a zero of . Recall that Newton’s method uses the recurrence relation , where function is defined by

where is the first derivative of . Hence, we obtain

Thus, we have the generalized Heron’s method

If we make , we obtain the standard Heron’s method, as we expected. The following Haskell script implements this generalized Heron’s method using arbitrary-precision rational numbers (included in the Data.Ratio library) to represent and the sequence of approximations :

import Data.Ratio

-- define generalized Heron map
g :: Integer -> Rational -> Rational -> Rational
g n y x = (1 % n) * ((fromIntegral (n-1) * x) + (y / (x^(n-1))))

-- initial approximation
x0 :: Rational
x0 = 1 % 1

-- list of approximations of the n-th root of y
roots :: Integer -> Rational -> [Rational]
roots n y | n  > 1 = iterate (g n y) x0
          | n <= 1 = error "Invalid value of n!!"

We load this script into GHCi. Let us now perform some numerical experiments.

__________

Experiment #1

We would like to compute a rational approximation of . We have and . Here is a very brief GHCi session:

*Main Data.Ratio> take 6 (roots 2 2)
[1 % 1,3 % 2,17 % 12,577 % 408,665857 % 470832,
 886731088897 % 627013566048]

We obtain the following rational approximation after five iterations

which is the same approximation I obtained twice before.

__________

Experiment #2

We would like to compute a rational approximation of . We have and . Here is a very brief GHCi session:

*Main Data.Ratio> take 7 (roots 2 5)
[1 % 1,3 % 1,7 % 3,47 % 21,2207 % 987,4870847 % 2178309,
 23725150497407 % 10610209857723]

We obtain the following rational approximation after six iterations

How good is this approximation? Let us check:

*Main Data.Ratio> (23725150497407 % 10610209857723)^2 - 5
4 % 112576553224922323902744729
*Main Data.Ratio> 4 / 112576553224922323902744729
3.5531377408652764e-26

Not bad! We could now obtain a rational approximation of the famous golden ratio , as follows:

*Main Data.Ratio> let xs = roots 2 5
*Main Data.Ratio> let phis = map (/2) (map (+1) xs)
*Main Data.Ratio> take 7 phis
[1 % 1,2 % 1,5 % 3,34 % 21,1597 % 987,3524578 % 2178309,
 17167680177565 % 10610209857723]

We obtain the following rational approximation after six iterations

How good is this approximation? Let us recall that the golden ratio is one of the solutions of the equation . Hence, we can check how close to zero is:

*Main Data.Ratio> let phi6 = 17167680177565 % 10610209857723
*Main Data.Ratio> phi6^2 - phi6 - 1
1 % 112576553224922323902744729
*Main Data.Ratio> 1 / 112576553224922323902744729
8.882844352163191e-27

I would say that is fairly close.

__________

Experiment #3

We now would like to compute a rational approximation of . We have and . Here is a very brief GHCi session:

*Main Data.Ratio> take 6 (roots 3 10)
[1 % 1,4 % 1,23 % 8,4909 % 2116,55223315303 % 25495981298,
83759169926117983945469262167029 % 38876457805393768546966848104041]

We obtain the following rational approximation after five iterations

Is this approximation any good? Let us check:

*Main Data.Ratio> let x5 = (roots 3 10) !! 5 
*Main Data.Ratio> x5^3
5876207108075025254603649658428809744645894398375582704
87986732898174481543266076676033790365389 % 58757060813
2677736272189042625165443146909261971794220978628909441
43919908404893015241356540921
*Main Data.Ratio> 5876207108075025254603649658428809744
6458943983755827048798673289817448154326607667603379036
5389 / 587570608132677736272189042625165443146909261971
79422097862890944143919908404893015241356540921
10.000852709004354

This is somewhat disappointing. The rational approximation is a ratio of two enormously long integers, but we do not have such a small error. Note that . We can do some quick error analysis. Let . The 1st order Taylor approximation of around evaluated at is

and, since , we obtain the estimate of the error

In other words, the output deviation should be divided by the slope to obtain an estimate of the magnitude of the approximation error . Let us use GHCi yet once again:

*Main Data.Ratio> let error = (x5^3 - 10) / (3 * x5^2)
*Main Data.Ratio> error
1670089160826306272530773923851043922672595525468316978
5941152245094153072382174540074985393 % 272741620880842
8590518540423731512512814773109630109912907276686930689
12248623218095571481624481
*Main Data.Ratio> 1670089160826306272530773923851043922
6725955254683169785941152245094153072382174540074985393
 / 2727416208808428590518540423731512512814773109630109
91290727668693068912248623218095571481624481
6.123338108179484e-5

This is still quite a huge error, I would say.

__________

Bonus Experiment

Let us go crazy and attempt to compute a rational approximation of . We have and . We expect the rational numbers to quickly become ratios of astronomically long integers. Therefore, we will not be showing any rational approximations on GHCi. Here is something scary:

*Main Data.Ratio> let x3 = (roots 100 100) !! 3
*Main Data.Ratio> let error = (x3^100 - 100) / (100 * x3^99)
*Main Data.Ratio> let error_float = fromIntegral (numerator error) 
                                  / fromIntegral (denominator error)
*Main Data.Ratio> error_float

After some 20 minutes without an output, I decided to abort this experiment. If you want to show x3, you can, but you will see thousands of digits. Thousands! You have been warned.

Heron’s method using integer arithmetic

Suppose we are given a positive rational number , whose square root we would like to compute. To be more precise, we would like to compute an approximation of .

Let be an approximation of . Hence, we have that  is the geometric mean of and . If we replace the geometric mean of these two numbers by their arithmetic mean, we hopefully obtain a “better” approximation , as follows

If all goes well, then will be “sufficiently close” to for a “sufficiently large” integer . This algorithm is the famous Heron’s method, devised by Heron of Alexandria some 2000 years ago. Some call it the Babylonian method, which suggests that the algorithm may be some 4000 years old. It’s not merely vintage, it’s antique.

Nonetheless, this is not a post on archaeology. We don’t want to sit around and talk about Heron’s method, we want to actually implement it. What programming language should we use? Let us use Haskell. How are we going to represent and the sequence of approximations ?

__________

Heron’s method using floating-point arithmetic

Though millions of people have already implemented Heron’s method using floating-point arithmetic (I implemented it in C in late 2000), one more implementation would do no harm:

y :: Floating a => a
y = 2

-- define Heron map
g :: Floating a => a -> a
g x = 0.5 * (x + y/x)

-- initial approximation
x0:: Floating a => a
x0 = 1

-- list of approximations
xs :: Floating a => [a]
xs = iterate g x0

Let us load this script into GHCi. Here’s a GHCi session:

*Main> take 6 xs
[1.0,1.5,1.4166666666666665,1.4142156862745097,
 1.4142135623746899,1.414213562373095]
*Main> let x5 = xs !! 5
*Main> x5**2 - 2
-4.440892098500626e-16

We thus have the following approximation after five iterations

__________

Heron’s method using rational arithmetic

Let us now use rational numbers of arbitrary precision. The following script uses the Data.Ratio library:

import Data.Ratio

y :: Rational
y = 2

-- define Heron map
g :: Rational -> Rational
g x = 0.5 * (x + y/x)

-- initial approximation
x0:: Rational
x0 = 1 % 1

-- list of approximations
xs :: [Rational]
xs = iterate g x0

Let us load this script into GHCi. Here’s a GHCi session:

*Main> take 6 xs
[1 % 1,3 % 2,17 % 12,577 % 408,665857 % 470832,
 886731088897 % 627013566048]
*Main> let x5 = xs !! 5
*Main> x5*x5 - 2
1 % 393146012008229658338304
*Main> 1 / 393146012008229658338304
2.5435842395854372e-24

Note that the error is eight orders of magnitude smaller than in the previous case (where we used finite-precision floating-point numbers). We pay for more accuracy with more computation. Thus, after five iterations, we have the following rational approximation

This is probably not news to you, my dear reader. Last April, I did, in fact, implement Heron’s method under the disguise of Newton’s method (see the appendix).

__________

Heron’s method using integer arithmetic

A rational number  can be expressed as a fraction , where and are integers and . Note that can (obviously?) be represented as a pair of integers . Since Haskell has arbitrary-precision integers, we can implement Heron’s method using pairs of arbitrary-precision integers to represent and the sequence of approximations .

Let and . The recurrence relation

can thus be rewritten in the following form

We then have two coupled recurrence relations

where , ,  and are integers. Here is a Haskell script that implements this vector recurrence relation:

a, b :: Integer
a = 2
b = 1

-- define Heron map
g :: (Integer,Integer) -> (Integer,Integer)
g (p,q) = (b * p^2 + a * q^2, 2 * b * p * q)

-- initial approximation
p0, q0 :: Integer
p0 = 1
q0 = 1

-- list of approximations
xs :: [(Integer,Integer)]
xs = iterate g (p0,q0)

Let us load this script into GHCi. Here’s a GHCi session:

*Main> take 6 xs
[(1,1),(3,2),(17,12),(577,408),(665857,470832),
 (886731088897,627013566048)]
*Main> let x5 = xs !! 5
*Main> x5
(886731088897,627013566048)

Hence, after five iterations, we obtain the very same rational approximation we obtained in the previous case (where we used rational numbers of arbitrary-precision)

It should be noted, however, that this implementation using arbitrary-precision integers is different in one very important aspect from the one using arbitrary-precision rational numbers. What is the difference? Let GHCi answer this question:

*Main> import Data.Ratio
*Main Data.Ratio> 2 % 2
1 % 1
*Main Data.Ratio> 4 % 4
1 % 1
*Main Data.Ratio> 3 % 6
1 % 2

The main difference between the two implementations is that when we use arbitrary-precision rational numbers, there is automatic reduction of every fraction to its irreducible form, which requires computing the greatest common divisor (GCD) of the numerator and denominator of each rational number at each iteration. This has a computational cost, of course, but it also has the benefit of avoiding reducible fractions (which are inefficient representations of rational numbers). When we use pairs of arbitrary-precision integers, we do not compute the GCD (although we could, if we wanted to), but we may have to work with pairs of unnecessarily long integers (that correspond to reducible fractions), which is potentially dangerous.

For example, suppose that we replace the first three lines in the script above with the following three lines:

a, b :: Integer
a = 2000
b = 1000

We run the script again. Here’s yet another GHCi session:

 
*Main> take 6 xs
[(1,1),(3000,2000),(17000000000,12000000000),
(577000000000000000000000,408000000000000000000000),
(665857000000000000000000000000000000000000000000000,
470832000000000000000000000000000000000000000000000),
(8867310888970000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000,
62701356604800000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000)]
*Main> let x5 = xs !! 5
*Main> x5
(8867310888970000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000,
62701356604800000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000)

Frankly, I am shocked! Multiplying and by results in and being multiplied by . I may have made a mistake counting all those zeros, but it’s a lot more zeros than I expected.

__________

Appendix: from Heron to Newton

Interestingly, Heron’s method can be derived from the no less famous Newton’s method. Let us introduce a function , defined by . By construction, the zeros of function are and . Recall that Newton’s method uses the recurrence relation , where function is defined by

where is the first derivative of . Hence, we obtain the following

which is the recurrence relation used in Heron’s method.

Deterministic Finite Automata in Haskell

Consider the deterministic finite automaton [1] illustrated below

[ state diagram courtesy of Michael Sipser [1] ]

Henceforth, we will call this automaton . Note that has three states, labeled , , and . The start state is (note the arrow coming from nowhere) and the accept state is (note the double circle). There are two possible inputs, labeled and , and, depending on the input, the automaton will jump from one state to another. In the diagram above these state transitions are depicted using arrows connecting two states.

Suppose that the automaton receives an input string such as (which the automaton reads from left to right). Since the start state is and the first input symbol is , the automaton will jump to state . The second input symbol is , so the automaton will remain at . The automaton reads the third input symbol, which is , and then jumps from state to state . The last input symbol is and thus the automaton jumps from state to state . Hence, the state sequence corresponding to the input string is the following

After reading the last symbol in the input string, produces an output. Since the final state is the accept state, we have that the automaton produces the output “accept”. Were the final state not , the automaton would produce the output “reject”. We conclude that this automaton accepts the input string . What other input strings does accept? Michael Sipser answers this question in [1]:

Experimenting with this machine on a variety of input strings reveals that it accepts the strings , , , and . In fact, accepts any string that ends with a , as it goes to its accept state whenever it reads the symbol . In addition, it accepts strings , , , and , and any string that ends with an even number of s following the last . It rejects other strings, such as , , .

A set of strings is called a language [1]. The set of all input strings accepted by the deterministic finite automaton is a language which we denote by .

__________

Formal definition

A deterministic finite automaton (DFA) consists of a finite set of states , a finite input alphabet that tells us what the allowed input symbols are, a transition function  that tells us how to jump from one state to another, a start state , and a set of accept states . A deterministic finite automaton (DFA) is thus a -tuple of the form

The deterministic finite automaton which we discussed earlier in this post is defined formally as follows

where the transition function is defined enumeratively as

,        ,       

,        ,       

Alternatively, we could view each state transition as an ordered triple , which might be easier to implement in some programming languages.

__________

Computing state sequences

Given an input string over a given alphabet , how do we obtain the corresponding state sequence? I solved that problem last January. I will not repeat myself here, but the crux of the matter is that the final state can be obtained using a left-fold

whereas the entire state sequence can be computed using a left-scan

Lacking a better name, I called this procedure the “scanl trick“, as I used the Haskell function scanl to implement it. Please let me know if you find a more sophisticated name for this “trick”.

__________

Haskell implementation of the DFA

Without further ado, here is a Haskell script:

data State = Q1 | Q2 | Q3 deriving (Read, Show, Eq)

type Input = Integer

-- define state-transition function
delta :: State -> Input -> State
delta Q1 0 = Q1
delta Q1 1 = Q2
delta Q2 0 = Q3
delta Q2 1 = Q2
delta Q3 0 = Q2
delta Q3 1 = Q2
delta  _ _ = error "Invalid input!"

-- define initial state
initialstate :: State
initialstate = Q1

-- create list of accept states
acceptstates :: [State]
acceptstates = [Q2]

-- create infinite list of input sequences
inputss :: [[Input]]
inputss = concat $ iterate g [[]]
          where g = concat . map (\xs -> [ xs ++ [s] | s <- [0,1]])

-- create accept predicate
isAccepted :: [Input] -> Bool
isAccepted inputs = finalstate `elem` acceptstates
                    where finalstate = foldl delta initialstate inputs

-- compute language recognized by the DFA
language :: [[Input]]
language = filter isAccepted inputss

Some remarks about this script are in order. Please note that:

• Sets are represented by lists. Strings are represented by lists, too. The latter is more natural than the former. Sets of strings become lists of lists.

• A new data type is created to represent the states, which we denote by Q1, Q2, and Q3. The input symbols are integers.

• Note that Input is a type synonym, inputs is a list of inputs symbols (i.e., an input string), and inputss is a list of lists of inputs (i.e., a list of input strings). Yes, it is a bit confusing.

• Note also that inputss is generated using the procedure I blogged about yesterday where the alphabet is .

• The final state is computed using the higher-order function foldl (left-fold). We check if the final state is an accept state using  the list membership predicate elem.

If you have any objections to this script, please do let me know. Let us now test it! We load it into GHCi and voilà:

*Main> -- check list of input strings
*Main> take 31 $ inputss
[[],[0],[1],
 [0,0],[0,1],[1,0],[1,1],
 [0,0,0],[0,0,1],[0,1,0],[0,1,1],
 [1,0,0],[1,0,1],[1,1,0],[1,1,1],
 [0,0,0,0],[0,0,0,1],[0,0,1,0],[0,0,1,1],
 [0,1,0,0],[0,1,0,1],[0,1,1,0],[0,1,1,1],
 [1,0,0,0],[1,0,0,1],[1,0,1,0],[1,0,1,1],
 [1,1,0,0],[1,1,0,1],[1,1,1,0],[1,1,1,1]]
*Main> -- compute the language of the automaton
*Main> -- (let us extract some 40 input strings only)
*Main> take 40 language
[[1],[0,1],[1,1],
 [0,0,1],[0,1,1],[1,0,0],[1,0,1],[1,1,1],
 [0,0,0,1],[0,0,1,1],[0,1,0,0],[0,1,0,1],
 [0,1,1,1],[1,0,0,1],[1,0,1,1],[1,1,0,0],
 [1,1,0,1],[1,1,1,1],[0,0,0,0,1],[0,0,0,1,1],
 [0,0,1,0,0],[0,0,1,0,1],[0,0,1,1,1],[0,1,0,0,1],
 [0,1,0,1,1],[0,1,1,0,0],[0,1,1,0,1],[0,1,1,1,1],
 [1,0,0,0,0],[1,0,0,0,1],[1,0,0,1,1],[1,0,1,0,0],
 [1,0,1,0,1],[1,0,1,1,1],[1,1,0,0,1],[1,1,0,1,1],
 [1,1,1,0,0],[1,1,1,0,1],[1,1,1,1,1],[0,0,0,0,0,1]]

The list of input strings is exactly the same I posted yesterday. The input strings in the language , or, to put it more precisely, the input strings in that are displayed above either end with a , or end “with an even number of s following the last “, as mentioned by Sipser in [1]. Why is that? The last in the input string puts in state , and an even number of s after that lead to transitions from to and back to , e.g.,

or

or

It would also be interesting to take a look at the state sequences corresponding to the input strings in . We compute these state sequences using the infamous “scanl trick“:

*Main> -- create list of sequences of states
*Main> let statess = map (\xs -> scanl delta initialstate xs) language
*Main> -- take the "first" 20 state sequences
*Main> take 20 statess
[[Q1,Q2],[Q1,Q1,Q2],[Q1,Q2,Q2],
 [Q1,Q1,Q1,Q2],[Q1,Q1,Q2,Q2],[Q1,Q2,Q3,Q2],
 [Q1,Q2,Q3,Q2],[Q1,Q2,Q2,Q2],[Q1,Q1,Q1,Q1,Q2],
 [Q1,Q1,Q1,Q2,Q2],[Q1,Q1,Q2,Q3,Q2],[Q1,Q1,Q2,Q3,Q2],
 [Q1,Q1,Q2,Q2,Q2],[Q1,Q2,Q3,Q2,Q2],[Q1,Q2,Q3,Q2,Q2],
 [Q1,Q2,Q2,Q3,Q2],[Q1,Q2,Q2,Q3,Q2],[Q1,Q2,Q2,Q2,Q2],
 [Q1,Q1,Q1,Q1,Q1,Q2],[Q1,Q1,Q1,Q1,Q2,Q2]]

Note that all state trajectories end in state Q2, as we expected.

__________

References

[1] Michael Sipser, Introduction to the Theory of Computation (2nd edition), Thomson Course Technology, 2006.