# Computing a Discrete Cosine Transform matrix

I'm learning Haskell and I needed to work with discrete cosine transform matrices, so I made a program that generates text (usable in WolframAlpha) containing the generated matrix.

The elements of the n-square matrix C are defined as:

$[C]_{jk} = \gamma_jcos\dfrac{\pi(2k+1)j}{2n}$ where $\gamma_{j} = \begin{cases}\frac{1}{\sqrt{n}} & \text{j = 0} \\ \frac{2}{\sqrt{n}} & \text{otherwise} \end{cases}$

I also wanted the program to simplify fractions and sqrt, which is where most of the "dirtiness" came from. Is there any way to clean this up? It's simple enough, but I want to learn how to be concise.

import Ratio

-- Pretty-print a fraction
showFraction n d
| n == 0            = "0"
| d == 0            = error "Division by zero."
| n == 1 && d == 1  = "1"
| n == 1            = "1/" ++ show d
| d == 1            = show n
| otherwise         = show n ++ "/" ++ show d

-- Pretty-print a fraction, reducing the fraction
showFractionReduced n d = let x = n % d in showFraction (numerator x) (denominator x)

-- Pretty-print a fraction, with an additional numerator factor
showFraction' n d n'
| n == 0            = "0"
| d == 0            = error "Division by zero."
| n == 1 && d == 1  = n'
| n == 1            = n' ++ "/" ++ show d
| d == 1            = show n ++ "*" ++ n'
| otherwise         = show n ++ "*" ++ n' ++ "/" ++ show d

-- Pretty-print a fraction, reducing the fraction, with an additional numerator factor
showFractionReduced' n d = let x = n % d in showFraction' (numerator x) (denominator x)

-- Discrete cosine transform matrix element
dctElement j k n
| n <= 1            = "1"
| n == 2 && j > 0   = "cos(" ++ showFractionReduced' (j * (2 * k + 1)) (2 * n) "pi" ++ ")"
| j == 0            = "sqrt(" ++ showFractionReduced 1 n ++ ")"
| otherwise         = "sqrt(" ++ showFractionReduced 2 n ++ ")*cos(" ++ showFractionReduced' (j * (2 * k + 1)) (2 * n) "pi" ++ ")"

-- Discrete cosine transform matrix row
dctRow j n = "{" ++ foldr1 (\acc x -> acc ++ "," ++ x) [dctElement j x n | x <- [0..n - 1]] ++ "}"

-- Discrete cosine transform matrix
dctMatrix n = "{" ++ foldr1 (\acc x -> acc ++ ",\n" ++ x) [dctRow x n | x <- [0..n - 1]] ++ "}"

-- Inverse discrete cosine transform matrix row
dctRowInverse j n = "{" ++ foldr1 (\acc x -> acc ++ "," ++ x) [dctElement x j n | x <- [0..n - 1]] ++ "}"

-- Inverse discrete cosine transform matrix
dctMatrixInverse n = "{" ++ foldr1 (\acc x -> acc ++ ",\n" ++ x) [dctRowInverse x n | x <- [0..n - 1]] ++ "}"


The string generated by dctMatrix 3, for example, is:

{{sqrt(1/3),sqrt(1/3),sqrt(1/3)},
{sqrt(2/3)*cos(pi/6),sqrt(2/3)*cos(pi/2),sqrt(2/3)*cos(5*pi/6)},
{sqrt(2/3)*cos(pi/3),sqrt(2/3)*cos(pi),sqrt(2/3)*cos(5*pi/3)}}

while the string generated by dctMatrix 2 is:

{{sqrt(1/3),sqrt(1/3),sqrt(1/3)},
{sqrt(2/3)*cos(pi/6),sqrt(2/3)*cos(pi/2),sqrt(2/3)*cos(5*pi/6)},
{sqrt(2/3)*cos(pi/3),sqrt(2/3)*cos(pi),sqrt(2/3)*cos(5*pi/3)}}

(Notice the simplification of the fractions in the former, and the square root in the latter.)

• You have a bug: 0/0 is 1, and not 0. Commented Apr 4, 2012 at 15:33

It is probably more in the spirit of functional programming to use pattern matches instead of long guard lists where possible. Your first function also has a few overlapping cases that could be eliminated in the interest of brevity:

showFraction _ 0 = error "Division by zero."
showFraction 0 _ = "0"
showFraction n 1 = show n
showFraction n d = show n ++ "/" ++ show d


Another simple change would be to exploit that foldr1 (\acc x -> acc ++ "," ++ x) is simply intercalate ',', which you get from Data.List. In a somewhat similar fashion, you could also replace foldr1 (\acc x -> acc ++ ",\n" ++ x) by concat . intersperse ",\n".

The result would look like follows:

dctRow j n = "{" ++ intercalate ',' [dctElement j x n | x <- [0..n - 1]] ++ "}"
dctMatrix n = "{" ++ concat (intersperse ",\n" [dctRow x n | x <- [0..n - 1]]) ++ "}"


Which in my opinion is more readable than the manual folds.

A more involved change would be to try to get rid of the $$\O(n^2)\$$ behaviour you will get from overusing (++). The background here is that (++) needs to make a full copy of the first operand in order to append the second one.

A fun functional way of improving this is to replace string concatenation by function composition:

showFraction :: Int -> Int -> String -> String
showFraction _ 0 = error "Division by zero."
showFraction 0 _ = ('0':)
showFraction n 1 = shows n
showFraction n d = shows n . ('/':) . shows d


So instead of a function that constructs a string, this is a function that prepends the string representation to a string passed as parameter. The nice thing about this is that you can simply compose these functions together in order to build bigger strings. As prepending doesn't require copying of the "tail" end in Haskell, this means that every part of the result string will be constructed exactly once!

Here is how to write the folds in this style:

applyInter f = flip (foldr ($)) . intersperse f dctRow j n = ('{':) . applyInter (',':) [dctElement j x n | x <- [0..n - 1]] . ('}':)  (Realize that flip (foldr ($)) [f,g,h] is simply f . g . h)

Note however that while this indeed has $$\O(n)\$$ complexity, it is still not very efficient. This is due to GHC probably selecting a bad arity for some of the functions in question, as well as String not being very fast in the first place. For constructing strings fast, it's probably a better idea to use Data.Text and/or a special builder library like blaze-builder.

• Great answer, thanks a bunch. I'll look into Data.Text after I implement all your other ideas. Commented Apr 5, 2012 at 3:23