Koch’s Snowflake
December 23, 2011
We produce output by writing to a PostScript file. Here’s the main program that writes the file header, calls the koch function to produce the snowflake, then writes the file footer:
(define (draw-koch n)
(display "%!") (newline)
(display "newpath") (newline)
(display "426 316 moveto") (newline)
(display "90 rotate") (newline)
(display "2 setlinewidth") (newline)
(koch n 243)
(display "stroke") (newline)
(display "showpage") (newline))
The initial coordinates 426 and 316 are arbitrary, but seem to work well. The koch function computes and writes the snowflake:
(define (koch n len)
(let koch ((pat (string->list "F++F++F")) (n n) (len len))
(if (positive? n)
(let ((new (replace pat #\F (string->list "F-F++F-F"))))
(koch new (- n 1) (/ len 3)))
(let draw ((pat pat))
(when (pair? pat)
(display "currentpoint translate") (newline)
(case (car pat)
((#\F) (display 0) (display " ") (display len)
(display " rlineto") (newline))
((#\-) (display "60 rotate") (newline))
((#\+) (display "-60 rotate") (newline)))
(draw (cdr pat)))))))
The outer loop koch computes the rewrite steps to the nth fractal; the inner loop draw writes the PostScript output. The main trick in the PostScript output is the command currentpoint translate, which moves the coordinate axes so each line can start at 0. The replace function modifies a list:
(define (replace xs old new)
(let loop ((xs xs) (zs (list)))
(cond ((null? xs) (reverse zs))
((equal? (car xs) old)
(loop (cdr xs) (append (reverse new) zs)))
(else (loop (cdr xs) (cons (car xs) zs))))))
To produce the 4-fractal PostScript output, say:
> (with-output-to-file "koch.ps"
(lambda () (draw-koch 4)))
You can run the code and see the PostScript output at http://programmingpraxis.codepad.org/vDdsP7o6; the image produced by the program is shown above right.
Merry Christmas!
Programming Praxis 
I noticed there are only six possible directions, and with each axis scaled properly, the arithmetic is very simple. The following computes the coordinates for a given list of moves, F, +, or -, starting at origin and initial angle of 60 degrees. (The x and y coordinates are in separate lists and in the backwards order.)
(define (kochdinates moves)
(let loop ((x 0) (y 0) (xs '(0)) (ys '(0))
(theta 60) (moves moves))
(if (null? moves)
(values xs ys)
(case (car moves)
((F) (call-with-values
(lambda ()
(case theta
(( 0) (values (+ x 2) y))
(( 60) (values (+ x 1) (+ y 1)))
((120) (values (- x 1) (+ y 1)))
((180) (values (- x 2) y))
((240) (values (- x 1) (- y 1)))
((300) (values (+ x 1) (- y 1)))))
(lambda (x1 y1)
(loop x1 y1 (cons x1 xs) (cons y1 ys)
theta (cdr moves)))))
((+) (loop x y xs ys (modulo (- theta 60) 360) (cdr moves)))
((-) (loop x y xs ys (modulo (+ theta 60) 360) (cdr moves)))))))
Here’s the code that produces the snowflake moves for the coordinate computer above. Nest at will.
(define (init) '(F + + F + + F))
(define (rewrite moves)
(apply append
(map (lambda (move)
(case move
((F) '(F - F + + F - F))
((+) '(+))
((-) '(-))))
moves)))
And here’s a compiler that produces an R command to plot a snowflake.
(define (R xs ys)
(display "plot(c(") (display (car xs))
(for-each (lambda (x) (display ",") (display x)) (cdr xs))
(display "), c(") (display (car ys))
(for-each (lambda (y) (display ",") (display y)) (cdr ys))
(display "), type='l', xlab='', ylab='')") (newline))
For example, this:
(call-with-values (lambda () (kochdinates (rewrite (init)))) R)
prints this:
plot(c(0,2,3,4,6,5,6,4,3,2,0,1,0), c(0,0,-1,0,0,1,2,2,3,2,2,1,0), type='l', xlab='', ylab='')
which can be pasted in an R session to open the plot window. For scaling on the screen, grab the corner of the plot window and let the computer do the arithmetic. (I ran Guile and R in one Emacs shell session, so cutting and pasting was easy, but I had to close one interpreter to enter the other. And the R command line editor seemed to choke when there were too many coordinates :)
[…] lovely programming exercise. I took the task one step further and created an interactive explorer of variants of the Koch […]
Use your favourite LOGO interpreter…
(defparameter *koch* '(f + + f + + f)) (defun koch-iterate (koch iterations) (loop while (> iterations 0) do (setq iterations (1- iterations)) (setq koch (apply #'append (mapcar (λ(v) (case v ((f) '(f - f + + f - f)) (t `(,v)))) koch)))) koch) (defun 2logo () (with-output-to-string (s) (mapcar (λ(v) (case v ((f) (format s "fd 5 ")) ((+) (format s "rt 60 ")) ((-) (format s "lt 60 ")))) (koch-iterate *koch* 5))))Here’s a Haskell version, using version 0.4 of the diagrams library. The heart of the program is the
stepfunction. Given an input shape,_, it produces_/\_. This is iterated some number of times (theedgefunction), starting with a single, straight line segment. Lastly, three copies of the final list of segments are glued together in a triangular shape.import Diagrams.Prelude import Diagrams.Backend.Cairo.CmdLine koch n = let e = edge !! n t = mconcat [e, rotateBy (-1/3) e, rotateBy (-2/3) e] in pathFromTrail $ scale (1/3^n) t where edge = iterate step $ fromSegments [straight (1,0)] step p = mconcat [p, rotateBy (1/6) p, rotateBy (-1/6) p, p] main = defaultMain . attrs . stroke $ koch 5 where attrs = lw 0.002 . pad 1.1Here are a few examples of running the program to generate PNG, PDF and SVG output.