-
Notifications
You must be signed in to change notification settings - Fork 0
/
index.html
638 lines (403 loc) · 45 KB
/
index.html
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
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
<!DOCTYPE html>
<!--[if IEMobile 7 ]><html class="no-js iem7"><![endif]-->
<!--[if lt IE 9]><html class="no-js lte-ie8"><![endif]-->
<!--[if (gt IE 8)|(gt IEMobile 7)|!(IEMobile)|!(IE)]><!--><html class="no-js" lang="en"><!--<![endif]-->
<head>
<meta charset="utf-8">
<title>Backing Store</title>
<meta name="author" content="Dan Brown">
<meta name="description" content="Last time I got my Haskell prime generator down to 17 seconds, compared to 2.5 seconds for the C version. Here is the Haskell code I wound up with: …">
<!-- http://t.co/dKP3o1e -->
<meta name="HandheldFriendly" content="True">
<meta name="MobileOptimized" content="320">
<meta name="viewport" content="width=device-width, initial-scale=1">
<link rel="canonical" href="http://dbbnrl.github.com/">
<link href="/favicon.png" rel="icon">
<link href="/stylesheets/screen.css" media="screen, projection" rel="stylesheet" type="text/css">
<script src="/javascripts/modernizr-2.0.js"></script>
<script src="/javascripts/ender.js"></script>
<script src="/javascripts/octopress.js" type="text/javascript"></script>
<link href="/atom.xml" rel="alternate" title="Backing Store" type="application/atom+xml">
<!--Fonts from Google"s Web font directory at http://google.com/webfonts -->
<link href="http://fonts.googleapis.com/css?family=PT+Serif:regular,italic,bold,bolditalic" rel="stylesheet" type="text/css">
<link href="http://fonts.googleapis.com/css?family=PT+Sans:regular,italic,bold,bolditalic" rel="stylesheet" type="text/css">
</head>
<body >
<header role="banner"><hgroup>
<h1><a href="/">Backing Store</a></h1>
<h2>Swapping out my brain...</h2>
</hgroup>
</header>
<nav role="navigation"><ul class="subscription" data-subscription="rss">
<li><a href="/atom.xml" rel="subscribe-rss" title="subscribe via RSS">RSS</a></li>
</ul>
<form action="http://google.com/search" method="get">
<fieldset role="search">
<input type="hidden" name="q" value="site:dbbnrl.github.com" />
<input class="search" type="text" name="q" results="0" placeholder="Search"/>
</fieldset>
</form>
<ul class="main-navigation">
<li><a href="/">Blog</a></li>
<li><a href="/blog/archives">Archives</a></li>
</ul>
</nav>
<div id="main">
<div id="content">
<div class="blog-index">
<article>
<header>
<h1 class="entry-title"><a href="/blog/2013/07/04/computing-primes-in-haskell-part-3/">Computing Primes in Haskell - Part 3</a></h1>
<p class="meta">
<time datetime="2013-07-04T18:16:00-04:00" pubdate data-updated="true">Jul 4<span>th</span>, 2013</time>
| <a href="/blog/2013/07/04/computing-primes-in-haskell-part-3/#disqus_thread">Comments</a>
</p>
</header>
<div class="entry-content"><p><a href="/blog/2013/07/03/computing-primes-in-haskell-part-2/">Last time</a> I got my Haskell prime generator down to 17 seconds, compared to 2.5 seconds for the <a href="/blog/2013/06/25/computing-primes-in-haskell/">C version</a>. Here is the Haskell code I wound up with:</p>
<figure class='code'><figcaption><span>Version 4 [17s]</span></figcaption><div class="highlight"><table><tr><td class="gutter"><pre class="line-numbers"><span class='line-number'>1</span>
<span class='line-number'>2</span>
<span class='line-number'>3</span>
<span class='line-number'>4</span>
<span class='line-number'>5</span>
<span class='line-number'>6</span>
<span class='line-number'>7</span>
<span class='line-number'>8</span>
</pre></td><td class='code'><pre><code class='haskell'><span class='line'><span class="nf">isPrime4</span> <span class="n">x</span> <span class="ow">=</span> <span class="n">not</span> <span class="o">$</span> <span class="n">any</span> <span class="p">(</span><span class="o">==</span><span class="mi">0</span><span class="p">)</span> <span class="o">$</span> <span class="n">map</span> <span class="p">(</span><span class="n">x</span> <span class="p">`</span><span class="n">rem</span><span class="p">`)</span> <span class="p">[</span><span class="mi">3</span><span class="p">,</span><span class="mi">5</span><span class="o">..</span><span class="n">floor</span> <span class="p">(</span><span class="n">sqrt</span> <span class="p">(</span><span class="n">fromIntegral</span> <span class="n">x</span><span class="p">))]</span>
</span><span class='line'>
</span><span class='line'><span class="nf">primes4</span> <span class="ow">::</span> <span class="p">[</span><span class="kt">Int</span><span class="p">]</span>
</span><span class='line'><span class="nf">primes4</span> <span class="ow">=</span> <span class="p">[</span><span class="n">n</span> <span class="o">|</span> <span class="n">n</span> <span class="ow"><-</span> <span class="p">[</span><span class="mi">3</span><span class="p">,</span><span class="mi">5</span><span class="o">..</span><span class="p">],</span> <span class="n">isPrime4</span> <span class="n">n</span><span class="p">]</span>
</span><span class='line'>
</span><span class='line'><span class="nf">main</span> <span class="ow">=</span> <span class="kr">do</span>
</span><span class='line'> <span class="kr">let</span> <span class="n">cnt</span> <span class="ow">=</span> <span class="mi">1000000</span>
</span><span class='line'> <span class="n">print</span> <span class="o">$</span> <span class="n">primes4</span> <span class="o">!!</span> <span class="n">cnt</span>
</span></code></pre></td></tr></table></div></figure>
<p>I was so focused on optimizing Haskell last time (both for performance and readability), that I forgot an important issue. I haven’t actually implemented the same algorithm as in the C version:</p>
<ul>
<li>The list of potential factors we’re using in <code>isprime4</code> to test for primality is the odd numbers less than or equal to <code>sqrt(x)</code>, rather than the <strong>prime</strong> numbers less than or equal to <code>sqrt(x)</code>.</li>
<li>The use of <code>sqrt</code> is probably quite expensive; the C version squares each factor instead.</li>
</ul>
<p>Let’s see if we can address those issues. First, get rid of <code>sqrt</code>:</p>
<figure class='code'><figcaption><span>Version 5 - No square root [46s]</span></figcaption><div class="highlight"><table><tr><td class="gutter"><pre class="line-numbers"><span class='line-number'>1</span>
<span class='line-number'>2</span>
<span class='line-number'>3</span>
<span class='line-number'>4</span>
<span class='line-number'>5</span>
</pre></td><td class='code'><pre><code class='haskell'><span class='line'><span class="nf">isPrime5</span> <span class="n">x</span> <span class="ow">=</span> <span class="n">not</span> <span class="o">$</span> <span class="n">any</span> <span class="p">(</span><span class="o">==</span><span class="mi">0</span><span class="p">)</span> <span class="o">$</span> <span class="n">map</span> <span class="p">(</span><span class="n">x</span> <span class="p">`</span><span class="n">rem</span><span class="p">`)</span> <span class="o">$</span> <span class="n">takeWhile</span> <span class="n">inRange</span> <span class="p">[</span><span class="mi">3</span><span class="p">,</span><span class="mi">5</span><span class="o">..</span><span class="p">]</span>
</span><span class='line'> <span class="kr">where</span> <span class="n">inRange</span> <span class="n">f</span> <span class="ow">=</span> <span class="p">(</span><span class="n">f</span><span class="o">*</span><span class="n">f</span><span class="p">)</span> <span class="o"><=</span> <span class="n">x</span>
</span><span class='line'>
</span><span class='line'><span class="nf">primes5</span> <span class="ow">::</span> <span class="p">[</span><span class="kt">Int</span><span class="p">]</span>
</span><span class='line'><span class="nf">primes5</span> <span class="ow">=</span> <span class="p">[</span><span class="n">n</span> <span class="o">|</span> <span class="n">n</span> <span class="ow"><-</span> <span class="p">[</span><span class="mi">3</span><span class="p">,</span><span class="mi">5</span><span class="o">..</span><span class="p">],</span> <span class="n">isPrime5</span> <span class="n">n</span><span class="p">]</span>
</span></code></pre></td></tr></table></div></figure>
<p>Since we’re avoiding <code>sqrt</code>, we can’t simply take odd integers up to a known upper bound, so we have to be a bit more creative. <code>takeWhile</code> will pull from the (infinite) list of odd numbers while a condition is true.</p>
<p>Unfortunately, the runtime jumps way up to 46 seconds. What’s the deal? After some investigation, my conclusion is that fusion is not implemented for <code>takeWhile</code>, <a href="http://haskell.1045720.n5.nabble.com/Unexpected-list-non-fusion-td5069492.html">although it probably could be</a>. This is unfortunate, but let’s keep going — I expect fusion will also be incompatible with the next optimization. Will it be worth it?</p>
<p>Next, we need to figure out a way to supply a list of prime numbers to <code>isprime</code>, rather than using a list of odd numbers. Where can we find a list of prime numbers?</p>
<figure class='code'><figcaption><span>Version 6 - Only test prime factors [14s]</span></figcaption><div class="highlight"><table><tr><td class="gutter"><pre class="line-numbers"><span class='line-number'>1</span>
<span class='line-number'>2</span>
<span class='line-number'>3</span>
<span class='line-number'>4</span>
<span class='line-number'>5</span>
</pre></td><td class='code'><pre><code class='haskell'><span class='line'><span class="nf">isPrime6</span> <span class="n">x</span> <span class="ow">=</span> <span class="n">not</span> <span class="o">$</span> <span class="n">any</span> <span class="p">(</span><span class="o">==</span><span class="mi">0</span><span class="p">)</span> <span class="o">$</span> <span class="n">map</span> <span class="p">(</span><span class="n">x</span> <span class="p">`</span><span class="n">rem</span><span class="p">`)</span> <span class="o">$</span> <span class="n">takeWhile</span> <span class="n">inRange</span> <span class="n">primes6</span>
</span><span class='line'> <span class="kr">where</span> <span class="n">inRange</span> <span class="n">f</span> <span class="ow">=</span> <span class="p">(</span><span class="n">f</span><span class="o">*</span><span class="n">f</span><span class="p">)</span> <span class="o"><=</span> <span class="n">x</span>
</span><span class='line'>
</span><span class='line'><span class="nf">primes6</span> <span class="ow">::</span> <span class="p">[</span><span class="kt">Int</span><span class="p">]</span>
</span><span class='line'><span class="nf">primes6</span> <span class="ow">=</span> <span class="mi">3</span><span class="kt">:</span><span class="p">[</span><span class="n">n</span> <span class="o">|</span> <span class="n">n</span> <span class="ow"><-</span> <span class="p">[</span><span class="mi">5</span><span class="p">,</span><span class="mi">7</span><span class="o">..</span><span class="p">],</span> <span class="n">isPrime6</span> <span class="n">n</span><span class="p">]</span>
</span></code></pre></td></tr></table></div></figure>
<p>Version 6 really shows off what you can do with laziness. As our list of prime factors, we supply to <code>isprime6</code> the infinite list <code>primes6</code>, which is precisely the value we’re trying to compute!</p>
<p>This is a bit circular, and in fact I had to pull the first prime (3) out of the loop and explicitly tack it onto the front of the list to avoid an infinite recursion. After that point, however, <code>isprime6</code> is guaranteed never to walk farther into <code>primes6</code> than those values that have already been computed. (To be really operational about it: At any given time, <code>primes6</code> consists of a finite list of prime numbers followed by an unevaluated <em>thunk</em> which represents the rest of the infinite list. Due to the <code>inRange</code> termination condition, <code>isprime6</code> can never traverse the list far enough to reach the thunk.)</p>
<p>This version takes 14 seconds. That really isn’t as much of an improvement over version 4 as I’d expect, given the vastly smaller set of factors that are being examined. To isolate the effects of list fusion, let’s manually fuse version 6 (thus undoing the improvements in code clarity we achieved in versions 2 through 4):</p>
<figure class='code'><figcaption><span>Version 7 - Manual fusion [8s]</span></figcaption><div class="highlight"><table><tr><td class="gutter"><pre class="line-numbers"><span class='line-number'>1</span>
<span class='line-number'>2</span>
<span class='line-number'>3</span>
<span class='line-number'>4</span>
<span class='line-number'>5</span>
<span class='line-number'>6</span>
<span class='line-number'>7</span>
<span class='line-number'>8</span>
</pre></td><td class='code'><pre><code class='haskell'><span class='line'><span class="nf">isPrime7</span> <span class="n">x</span> <span class="ow">=</span> <span class="n">test</span> <span class="n">primes7</span>
</span><span class='line'> <span class="kr">where</span> <span class="n">test</span> <span class="p">(</span><span class="n">f</span><span class="kt">:</span><span class="n">fs</span><span class="p">)</span>
</span><span class='line'> <span class="o">|</span> <span class="p">(</span><span class="n">f</span><span class="o">*</span><span class="n">f</span><span class="p">)</span> <span class="o">></span> <span class="n">x</span> <span class="ow">=</span> <span class="kt">True</span>
</span><span class='line'> <span class="o">|</span> <span class="p">(</span><span class="n">x</span> <span class="p">`</span><span class="n">rem</span><span class="p">`</span> <span class="n">f</span><span class="p">)</span> <span class="o">==</span> <span class="mi">0</span> <span class="ow">=</span> <span class="kt">False</span>
</span><span class='line'> <span class="o">|</span> <span class="n">otherwise</span> <span class="ow">=</span> <span class="n">test</span> <span class="n">fs</span>
</span><span class='line'>
</span><span class='line'><span class="nf">primes7</span> <span class="ow">::</span> <span class="p">[</span><span class="kt">Int</span><span class="p">]</span>
</span><span class='line'><span class="nf">primes7</span> <span class="ow">=</span> <span class="mi">3</span><span class="kt">:</span><span class="p">[</span><span class="n">n</span> <span class="o">|</span> <span class="n">n</span> <span class="ow"><-</span> <span class="p">[</span><span class="mi">5</span><span class="p">,</span><span class="mi">7</span><span class="o">..</span><span class="p">],</span> <span class="n">isPrime7</span> <span class="n">n</span><span class="p">]</span>
</span></code></pre></td></tr></table></div></figure>
<p>Now we’re down to 8 seconds — a pretty massive improvement over 14 seconds, and only a factor of 3 or so from the C version’s performance. Although this is a ridiculously simple toy problem, there appears to be room for quite a bit of optimization.</p>
<p>I suspect I’ve taken this about as far as it can go using Haskell’s list datatype. The next step is to investigate array-like alternatives.</p>
<p>At this point, my conclusion is that a highly-readable, compositional coding style is sometimes at odds with performance, because GHC just isn’t quite smart enough <sup id="fnref:1"><a href="#fn:1" rel="footnote">1</a></sup>. I’d be interested in what more experienced Haskell programmers have to say on this subject.</p>
<p>As a side note, observant readers will notice that <code>test</code> in version 7 is a partial function; it will fail on the empty list. It turns out that it actually <strong>can’t</strong> fail, due to the <code>(f*f) > x</code> termination test — but the compiler has no way of knowing this. I wonder what the implications of this are from an optimization standpoint?</p>
<div class="footnotes">
<hr/>
<ol>
<li id="fn:1">
<p>In case anyone is wondering: I have tried every version so far with the LLVM backend to GHC, and observed no dramatic differences. Results were typically within one second of what I’ve reported here.<a href="#fnref:1" rev="footnote">↩</a></p></li>
</ol>
</div>
</div>
</article>
<article>
<header>
<h1 class="entry-title"><a href="/blog/2013/07/03/computing-primes-in-haskell-part-2/">Computing Primes in Haskell - Part 2</a></h1>
<p class="meta">
<time datetime="2013-07-03T16:46:00-04:00" pubdate data-updated="true">Jul 3<span>rd</span>, 2013</time>
| <a href="/blog/2013/07/03/computing-primes-in-haskell-part-2/#disqus_thread">Comments</a>
</p>
</header>
<div class="entry-content"><p><a href="/blog/2013/06/25/computing-primes-in-haskell/">Last time</a> I set myself a simple task for the purpose of exploring Haskell: computing prime numbers. My reference C implementation computed 1000000 primes on my desktop in 2.5s, using “gcc -O2”. Without further ado, here’s my first working Haskell implementation:</p>
<figure class='code'><figcaption><span>Version 1 - Ugly [58s]</span></figcaption><div class="highlight"><table><tr><td class="gutter"><pre class="line-numbers"><span class='line-number'>1</span>
<span class='line-number'>2</span>
<span class='line-number'>3</span>
<span class='line-number'>4</span>
<span class='line-number'>5</span>
<span class='line-number'>6</span>
<span class='line-number'>7</span>
<span class='line-number'>8</span>
<span class='line-number'>9</span>
<span class='line-number'>10</span>
<span class='line-number'>11</span>
<span class='line-number'>12</span>
</pre></td><td class='code'><pre><code class='haskell'><span class='line'><span class="nf">isPrime1</span> <span class="n">x</span> <span class="ow">=</span> <span class="n">test</span> <span class="mi">3</span>
</span><span class='line'> <span class="kr">where</span> <span class="n">test</span> <span class="n">factor</span>
</span><span class='line'> <span class="o">|</span> <span class="n">factor</span> <span class="o">></span> <span class="n">sqrtx</span> <span class="ow">=</span> <span class="kt">True</span>
</span><span class='line'> <span class="o">|</span> <span class="p">(</span><span class="n">x</span> <span class="p">`</span><span class="n">rem</span><span class="p">`</span> <span class="n">factor</span><span class="p">)</span> <span class="o">==</span> <span class="mi">0</span> <span class="ow">=</span> <span class="kt">False</span>
</span><span class='line'> <span class="o">|</span> <span class="n">otherwise</span> <span class="ow">=</span> <span class="n">test</span> <span class="p">(</span><span class="n">factor</span> <span class="o">+</span> <span class="mi">2</span><span class="p">)</span>
</span><span class='line'> <span class="n">sqrtx</span> <span class="ow">=</span> <span class="n">floor</span> <span class="p">(</span><span class="n">sqrt</span> <span class="p">(</span><span class="n">fromIntegral</span> <span class="n">x</span><span class="p">))</span>
</span><span class='line'>
</span><span class='line'><span class="nf">primes1</span> <span class="ow">=</span> <span class="p">[</span><span class="n">n</span> <span class="o">|</span> <span class="n">n</span> <span class="ow"><-</span> <span class="p">[</span><span class="mi">3</span><span class="p">,</span><span class="mi">5</span><span class="o">..</span><span class="p">],</span> <span class="n">isPrime1</span> <span class="n">n</span><span class="p">]</span>
</span><span class='line'>
</span><span class='line'><span class="nf">main</span> <span class="ow">=</span> <span class="kr">do</span>
</span><span class='line'> <span class="kr">let</span> <span class="n">cnt</span> <span class="ow">=</span> <span class="mi">1000000</span>
</span><span class='line'> <span class="n">print</span> <span class="o">$</span> <span class="n">primes1</span> <span class="o">!!</span> <span class="n">cnt</span>
</span></code></pre></td></tr></table></div></figure>
<p>(BTW, I’m aware that 2 is prime. This implementation (and the C version) leave out 2 for simplicity.)</p>
<p>Compiled with “ghc -O2”, this version runs in 58 seconds. So, we have a ways to go if we want to approach C.</p>
<p><code>primes1</code> is an easy-to-read list comprehension. I like the fact that it’s not a function, but an infinite list. On the other hand, <code>isPrime1</code> is pretty ugly and non-idiomatic Haskell code. Let’s see if we can transform it into something nicer.</p>
<figure class='code'><figcaption><span>Version 2 - List comprehension [73s]</span></figcaption><div class="highlight"><table><tr><td class="gutter"><pre class="line-numbers"><span class='line-number'>1</span>
<span class='line-number'>2</span>
<span class='line-number'>3</span>
<span class='line-number'>4</span>
<span class='line-number'>5</span>
<span class='line-number'>6</span>
<span class='line-number'>7</span>
<span class='line-number'>8</span>
<span class='line-number'>9</span>
<span class='line-number'>10</span>
<span class='line-number'>11</span>
</pre></td><td class='code'><pre><code class='haskell'><span class='line'><span class="nf">isPrime2</span> <span class="n">x</span> <span class="ow">=</span> <span class="n">test</span> <span class="p">[</span><span class="mi">3</span><span class="p">,</span><span class="mi">5</span><span class="o">..</span><span class="n">floor</span> <span class="p">(</span><span class="n">sqrt</span> <span class="p">(</span><span class="n">fromIntegral</span> <span class="n">x</span><span class="p">))]</span>
</span><span class='line'> <span class="kr">where</span> <span class="n">test</span> <span class="kt">[]</span> <span class="ow">=</span> <span class="kt">True</span>
</span><span class='line'> <span class="n">test</span> <span class="p">(</span><span class="n">f</span><span class="kt">:</span><span class="n">fs</span><span class="p">)</span>
</span><span class='line'> <span class="o">|</span> <span class="p">(</span><span class="n">x</span> <span class="p">`</span><span class="n">rem</span><span class="p">`</span> <span class="n">f</span><span class="p">)</span> <span class="o">==</span> <span class="mi">0</span> <span class="ow">=</span> <span class="kt">False</span>
</span><span class='line'> <span class="o">|</span> <span class="n">otherwise</span> <span class="ow">=</span> <span class="n">test</span> <span class="n">fs</span>
</span><span class='line'>
</span><span class='line'><span class="nf">primes2</span> <span class="ow">=</span> <span class="p">[</span><span class="n">n</span> <span class="o">|</span> <span class="n">n</span> <span class="ow"><-</span> <span class="p">[</span><span class="mi">3</span><span class="p">,</span><span class="mi">5</span><span class="o">..</span><span class="p">],</span> <span class="n">isPrime2</span> <span class="n">n</span><span class="p">]</span>
</span><span class='line'>
</span><span class='line'><span class="nf">main</span> <span class="ow">=</span> <span class="kr">do</span>
</span><span class='line'> <span class="kr">let</span> <span class="n">cnt</span> <span class="ow">=</span> <span class="mi">1000000</span>
</span><span class='line'> <span class="n">print</span> <span class="o">$</span> <span class="n">primes2</span> <span class="o">!!</span> <span class="n">cnt</span>
</span></code></pre></td></tr></table></div></figure>
<p>The first step is to generate the candidate factors using a list comprehension. The code is a bit more readable now. But the runtime has increased to 73 seconds! What’s going on?</p>
<p>Let’s defer that question for the moment and push on with the transformation. Next, we’ll pull the <code>`rem`</code> operation out of <code>test</code>:</p>
<figure class='code'><figcaption><span>Version 3 - Using “map” [84s]</span></figcaption><div class="highlight"><table><tr><td class="gutter"><pre class="line-numbers"><span class='line-number'>1</span>
<span class='line-number'>2</span>
<span class='line-number'>3</span>
<span class='line-number'>4</span>
<span class='line-number'>5</span>
<span class='line-number'>6</span>
<span class='line-number'>7</span>
</pre></td><td class='code'><pre><code class='haskell'><span class='line'><span class="nf">isPrime3</span> <span class="n">x</span> <span class="ow">=</span> <span class="n">test</span> <span class="o">$</span> <span class="n">map</span> <span class="p">(</span><span class="n">x</span> <span class="p">`</span><span class="n">rem</span><span class="p">`)</span> <span class="p">[</span><span class="mi">3</span><span class="p">,</span><span class="mi">5</span><span class="o">..</span><span class="n">floor</span> <span class="p">(</span><span class="n">sqrt</span> <span class="p">(</span><span class="n">fromIntegral</span> <span class="n">x</span><span class="p">))]</span>
</span><span class='line'> <span class="kr">where</span> <span class="n">test</span> <span class="kt">[]</span> <span class="ow">=</span> <span class="kt">True</span>
</span><span class='line'> <span class="n">test</span> <span class="p">(</span><span class="n">r</span><span class="kt">:</span><span class="n">rs</span><span class="p">)</span>
</span><span class='line'> <span class="o">|</span> <span class="n">r</span> <span class="o">==</span> <span class="mi">0</span> <span class="ow">=</span> <span class="kt">False</span>
</span><span class='line'> <span class="o">|</span> <span class="n">otherwise</span> <span class="ow">=</span> <span class="n">test</span> <span class="n">rs</span>
</span><span class='line'>
</span><span class='line'><span class="nf">primes3</span> <span class="ow">=</span> <span class="p">[</span><span class="n">n</span> <span class="o">|</span> <span class="n">n</span> <span class="ow"><-</span> <span class="p">[</span><span class="mi">3</span><span class="p">,</span><span class="mi">5</span><span class="o">..</span><span class="p">],</span> <span class="n">isPrime3</span> <span class="n">n</span><span class="p">]</span>
</span></code></pre></td></tr></table></div></figure>
<p>(I’ll leave out <code>main</code> from here forward to save space, unless it’s interesting.)</p>
<p>The <code>`rem`</code> is now computed by mapping <code>(x `rem`)</code>, which is an <a href="http://www.haskell.org/haskellwiki/Section_of_an_infix_operator">operator section</a>. The <code>$</code> operator is just function application with very low precedence; without it we’d have to put parentheses around the entire rest of the line.</p>
<p>Nicer, but the runtime has increased again, to 84 seconds. And again, I’m going to ignore the issue and continue rewriting. My goal, if it isn’t obvious, is to get rid of <code>test</code>:</p>
<figure class='code'><figcaption><span>Version 4 - Using “any” [83s]</span></figcaption><div class="highlight"><table><tr><td class="gutter"><pre class="line-numbers"><span class='line-number'>1</span>
<span class='line-number'>2</span>
<span class='line-number'>3</span>
</pre></td><td class='code'><pre><code class='haskell'><span class='line'><span class="nf">isPrime4</span> <span class="n">x</span> <span class="ow">=</span> <span class="n">not</span> <span class="o">$</span> <span class="n">any</span> <span class="p">(</span><span class="o">==</span><span class="mi">0</span><span class="p">)</span> <span class="o">$</span> <span class="n">map</span> <span class="p">(</span><span class="n">x</span> <span class="p">`</span><span class="n">rem</span><span class="p">`)</span> <span class="p">[</span><span class="mi">3</span><span class="p">,</span><span class="mi">5</span><span class="o">..</span><span class="n">floor</span> <span class="p">(</span><span class="n">sqrt</span> <span class="p">(</span><span class="n">fromIntegral</span> <span class="n">x</span><span class="p">))]</span>
</span><span class='line'>
</span><span class='line'><span class="nf">primes4</span> <span class="ow">=</span> <span class="p">[</span><span class="n">n</span> <span class="o">|</span> <span class="n">n</span> <span class="ow"><-</span> <span class="p">[</span><span class="mi">3</span><span class="p">,</span><span class="mi">5</span><span class="o">..</span><span class="p">],</span> <span class="n">isPrime4</span> <span class="n">n</span><span class="p">]</span>
</span></code></pre></td></tr></table></div></figure>
<p>Another operator section, used with the <code>any</code> function, which is just short-circuit boolean OR over a list. Version 4 is nice and clean, and much more pleasant to read than the C version. I particularly like how <code>isPrime4</code> reads almost like an English sentence: “X is prime if it is not the case that any elements equal zero in the list produced by…”</p>
<p>This version completes in 83 seconds. So… why was version 1 so much faster?</p>
<p>Version 4 is structured as a pipeline. Each stage consumes a list, performs an operation on each element, and produces a new list (or a single result, in the case of <code>any</code>). In a traditional non-lazy language, this code style would clearly be wasteful: First generate in memory a list of odd integers up to sqrt(x), then traverse that list and generate a new list with the <code>(x `rem`)</code> operation applied, then traverse <em>that</em> list looking for elements equal to zero.</p>
<p>In Haskell, it doesn’t really work that way due to lazy evaluation. Each list is computed only as it is needed, so the various operations proceed in lockstep (they are essentially coroutines). But there’s still a lot of overhead involved in allocating list elements, filling them, pulling items out of them, and freeing them.</p>
<p>Ideally, we’d like the compiler to be smart enough to <em>transform</em> version 4 into something like version 1. Why isn’t it?</p>
<p>Well, it is, if we add the following line above <code>primes4</code>:</p>
<figure class='code'><figcaption><span></span></figcaption><div class="highlight"><table><tr><td class="gutter"><pre class="line-numbers"><span class='line-number'>1</span>
</pre></td><td class='code'><pre><code class='haskell'><span class='line'><span class="nf">primes4</span> <span class="ow">::</span> <span class="p">[</span><span class="kt">Int</span><span class="p">]</span>
</span></code></pre></td></tr></table></div></figure>
<p>With this type annotation, version 4 now completes in 17 seconds. Why the difference? Well, since we didn’t originally provide a type for <code>primes4</code>, ghc was free to infer any reasonable type <sup id="fnref:1"><a href="#fn:1" rel="footnote">1</a></sup>. And it chose <code>Integer</code>, which is an arbitrary-sized integer. We really want to use something like <code>Int</code>, which uses the (bounded) native machine representation. With a type annotation declaring that we want a list of <code>Int</code>, the behavior changes dramatically.</p>
<p>Let’s add <code>primesX :: [Int]</code> type annotations to every version and compare:</p>
<ul>
<li>Version 1 with [Int]: 24 seconds</li>
<li>Version 2 with [Int]: 32 seconds</li>
<li>Version 3 with [Int]: 40 seconds</li>
<li>Version 4 with [Int]: 17 seconds</li>
</ul>
<p>As you can see, the time increases until we reach version 4, which is written entirely in terms of standard list processing functions (no <code>test</code> helper function). The compiler can now apply <a href="http://www.haskell.org/haskellwiki/Fusion">short cut fusion</a>. The result is now faster than the “manually” created implementation in version 1!</p>
<p>I was actually slightly surprised that ghc wouldn’t perform fusion for <code>Integer</code> lists. Perhaps the optimization is only implemented for unboxable types. Or, perhaps I’m completely wrong about what’s going on here. At some point, I should revisit this analysis using a tool like <em>ghc-core</em>.</p>
<p>Haskell is down to 17 seconds. Making progress!</p>
<div class="footnotes">
<hr/>
<ol>
<li id="fn:1">
<p>This is misleading, perhaps. Normally, I believe Haskell will try to infer <em>the most general type</em>. In this case, that would be something like <code>Integral a => [a]</code>, which means a list of any <code>a</code> such that <code>a</code> belongs to the <code>Integral</code> typeclass. However, inferring that type could result in <code>primes4</code> being <em>reevaluated</em> each time it is referenced, due to the hidden typeclass parameter. This could be surprising behavior, given that syntactically, <code>primes4</code> looks like a constant, not a function. Haskell’s <a href="http://www.haskell.org/haskellwiki/Monomorphism_restriction">monomorphism restriction</a> was designed to prevent programmer confusion in these cases by forcing the compiler to infer a monomorphic (non-polymorphic) type. Since <em>any</em> member of the <code>Integral</code> typeclass is an equally resonable choice for monomorphic type in this case, Haskell consults an explicit list of defaults to arrive at the choice of <code>Integer</code>.<a href="#fnref:1" rev="footnote">↩</a></p></li>
</ol>
</div>
</div>
</article>
<article>
<header>
<h1 class="entry-title"><a href="/blog/2013/06/25/computing-primes-in-haskell/">Computing Primes in Haskell</a></h1>
<p class="meta">
<time datetime="2013-06-25T22:50:00-04:00" pubdate data-updated="true">Jun 25<span>th</span>, 2013</time>
| <a href="/blog/2013/06/25/computing-primes-in-haskell/#disqus_thread">Comments</a>
</p>
</header>
<div class="entry-content"><p>I’ve been reading a ton of Haskell papers, blogs, and other resources. I’ve also made small inroads into doing OpenGL rendering in Haskell, but it’s become apparent that I should back off and tackle smaller mountains before I go there.</p>
<p>It occurred to me the other day that I should revisit an old memory. Back when I first taught myself to program (BASIC, on my old Apple IIe), I spent quite a bit of time optimizing a program to generate prime numbers.</p>
<p>This seems like a nice starting point for diving into Haskell. The initial implementation will likely be very trivial, but in the course of optimization I’m hoping this can lead to a beginner’s understanding of:</p>
<ul>
<li>Utilizing laziness</li>
<li>Dealing with the consequences of laziness</li>
<li>Writing idiomatic code</li>
<li>Benchmarking, profiling, and analyzing generated code (maybe learn to read GHC-core?)</li>
<li>Array/Vector libraries (assuming linked lists don’t scale well to this problem)</li>
<li>Using imperative (mutable) Array/Vector operations</li>
</ul>
<p>The algorithm I’ll be implementing and optimizing uses <a href="http://en.wikipedia.org/wiki/Trial_division">trial division</a> to enumerate the primes. This is not necessarily the best way to find primes, but that isn’t the point here — I’m trying to optimize the implementation, not find the best algorithm. Maybe I’ll revisit that question later.</p>
<p>Here is an C implementation that captures the algorithm:</p>
<figure class='code'><figcaption><span></span></figcaption><div class="highlight"><table><tr><td class="gutter"><pre class="line-numbers"><span class='line-number'>1</span>
<span class='line-number'>2</span>
<span class='line-number'>3</span>
<span class='line-number'>4</span>
<span class='line-number'>5</span>
<span class='line-number'>6</span>
<span class='line-number'>7</span>
<span class='line-number'>8</span>
<span class='line-number'>9</span>
<span class='line-number'>10</span>
<span class='line-number'>11</span>
<span class='line-number'>12</span>
<span class='line-number'>13</span>
<span class='line-number'>14</span>
<span class='line-number'>15</span>
<span class='line-number'>16</span>
<span class='line-number'>17</span>
<span class='line-number'>18</span>
<span class='line-number'>19</span>
<span class='line-number'>20</span>
<span class='line-number'>21</span>
<span class='line-number'>22</span>
<span class='line-number'>23</span>
<span class='line-number'>24</span>
<span class='line-number'>25</span>
<span class='line-number'>26</span>
</pre></td><td class='code'><pre><code class='c'><span class='line'><span class="cp">#define NUMPRIMES 1000000 </span><span class="c1">// Arbitrary bound -- ugly!</span>
</span><span class='line'><span class="kt">unsigned</span> <span class="n">primes</span><span class="p">[</span><span class="n">NUMPRIMES</span><span class="p">];</span>
</span><span class='line'><span class="kt">int</span> <span class="nf">main</span><span class="p">(</span><span class="kt">void</span><span class="p">)</span> <span class="p">{</span>
</span><span class='line'> <span class="kt">unsigned</span> <span class="n">candidate</span> <span class="o">=</span> <span class="mi">3</span><span class="p">;</span>
</span><span class='line'> <span class="n">primes</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span> <span class="o">=</span> <span class="n">candidate</span><span class="p">;</span>
</span><span class='line'> <span class="kt">unsigned</span> <span class="n">numsofar</span> <span class="o">=</span> <span class="mi">1</span><span class="p">;</span>
</span><span class='line'> <span class="k">while</span> <span class="p">(</span><span class="n">numsofar</span> <span class="o"><</span> <span class="n">NUMPRIMES</span><span class="p">)</span> <span class="p">{</span>
</span><span class='line'> <span class="n">candidate</span> <span class="o">+=</span> <span class="mi">2</span><span class="p">;</span> <span class="c1">// Ignore even numbers</span>
</span><span class='line'> <span class="kt">unsigned</span> <span class="n">i</span><span class="p">;</span>
</span><span class='line'> <span class="k">for</span> <span class="p">(</span><span class="n">i</span> <span class="o">=</span> <span class="mi">0</span><span class="p">;</span> <span class="n">i</span> <span class="o"><</span> <span class="n">numsofar</span><span class="p">;</span> <span class="n">i</span><span class="o">++</span><span class="p">)</span> <span class="p">{</span>
</span><span class='line'> <span class="kt">unsigned</span> <span class="n">p</span> <span class="o">=</span> <span class="n">primes</span><span class="p">[</span><span class="n">i</span><span class="p">];</span>
</span><span class='line'> <span class="k">if</span> <span class="p">((</span><span class="n">p</span><span class="o">*</span><span class="n">p</span><span class="p">)</span> <span class="o">></span> <span class="n">candidate</span><span class="p">)</span> <span class="p">{</span> <span class="c1">// Stop if we pass sqrt(candidate)</span>
</span><span class='line'> <span class="c1">// We found a prime</span>
</span><span class='line'> <span class="n">primes</span><span class="p">[</span><span class="n">numsofar</span><span class="o">++</span><span class="p">]</span> <span class="o">=</span> <span class="n">candidate</span><span class="p">;</span>
</span><span class='line'> <span class="k">break</span><span class="p">;</span>
</span><span class='line'> <span class="p">}</span>
</span><span class='line'> <span class="k">if</span> <span class="p">((</span><span class="n">candidate</span> <span class="o">%</span> <span class="n">p</span><span class="p">)</span> <span class="o">==</span> <span class="mi">0</span><span class="p">)</span> <span class="p">{</span>
</span><span class='line'> <span class="c1">// We found a prime factor, so candidate is not prime</span>
</span><span class='line'> <span class="k">break</span><span class="p">;</span>
</span><span class='line'> <span class="p">}</span>
</span><span class='line'> <span class="p">}</span>
</span><span class='line'> <span class="c1">// We can never reach this line!</span>
</span><span class='line'> <span class="p">}</span>
</span><span class='line'> <span class="n">printf</span><span class="p">(</span><span class="s">"%u</span><span class="se">\n</span><span class="s">"</span><span class="p">,</span> <span class="n">primes</span><span class="p">[</span><span class="n">NUMPRIMES</span><span class="o">-</span><span class="mi">1</span><span class="p">]);</span>
</span><span class='line'> <span class="k">return</span> <span class="mi">0</span><span class="p">;</span>
</span><span class='line'><span class="p">}</span>
</span></code></pre></td></tr></table></div></figure>
<p>Compiled with gcc -O2, this takes almost exactly 2.5s on my test machine. We can use that as an optimization goal when evaluating Haskell implementations.</p>
</div>
</article>
<article>
<header>
<h1 class="entry-title"><a href="/blog/2013/01/09/a-bit-late-to-the-party/">A Bit Late to the Party</a></h1>
<p class="meta">
<time datetime="2013-01-09T22:03:00-05:00" pubdate data-updated="true">Jan 9<span>th</span>, 2013</time>
| <a href="/blog/2013/01/09/a-bit-late-to-the-party/#disqus_thread">Comments</a>
</p>
</header>
<div class="entry-content"><p>So, I’ve finally decided to try out this blogging thing.</p>
<p>Recently I’ve been fiddling around on a project that (so far) involves simultaneously learning modern OpenGL and Haskell. And Functional Reactive Programming (FRP). Oh, and a bit of WXWindows, possibly.</p>
<p>All of which involves getting a lot of concepts straight in my head, which I’m thinking might be easiest if I take notes, some of which might actually be organized enough for other people to find useful, which suggests that I should look into blogging…</p>
<p>In my typical premature-optimization, tool-focused style, I created this blog on wordpress.com, discovered through further research that it wasn’t <em>quite</em> perfect for my needs, and abandoned it the same day (after zero posts) to start getting octopress-over-github configured. This way is much geekier, which makes it Better.</p>
<p>Oh, and I keep meaning to buy some parts and do a bit of Arduino-based robotics, though at this point I’ve spent so much time thinking about doing that without actually <em>doing</em> it (including hours of research and a carefully constructed shopping list) that Raspberry Pi may now be a better option. I should really research that…</p>
<p>There’s a pattern here, which I hope to break out of by:</p>
<ul>
<li>Doing things in small increments and committing code to Git frequently, and</li>
<li>Blogging frequently about what I’m working on.</li>
</ul>
<p>The former has helped my productivity at work quite a bit. (And helped curtail my tendency toward over-design and premature optimization). The latter is just the application of the same principle to mental output that isn’t code and doesn’t go into revision control (though with Github pages, there’s a certain additional symmetry…)</p>
<p>I’m not quite sure why I’m telling myself these things, it’s not as if the Internet is listening. I suppose my audience is really Future Dan.</p>
<p>Hello, Future Dan.</p>
</div>
</article>
<div class="pagination">
<a href="/blog/archives">Blog Archives</a>
</div>
</div>
<aside class="sidebar">
<section>
<h1>Recent Posts</h1>
<ul id="recent_posts">
<li class="post">
<a href="/blog/2013/07/04/computing-primes-in-haskell-part-3/">Computing primes in Haskell - part 3</a>
</li>
<li class="post">
<a href="/blog/2013/07/03/computing-primes-in-haskell-part-2/">Computing primes in Haskell - part 2</a>
</li>
<li class="post">
<a href="/blog/2013/06/25/computing-primes-in-haskell/">Computing primes in Haskell</a>
</li>
<li class="post">
<a href="/blog/2013/01/09/a-bit-late-to-the-party/">A bit late to the party</a>
</li>
</ul>
</section>
</aside>
</div>
</div>
<footer role="contentinfo"><p>
Copyright © 2013 - Dan Brown -
<span class="credit">Powered by <a href="http://octopress.org">Octopress</a></span>
</p>
</footer>
<script type="text/javascript">
var disqus_shortname = 'backingstore';
var disqus_script = 'count.js';
(function () {
var dsq = document.createElement('script'); dsq.type = 'text/javascript'; dsq.async = true;
dsq.src = 'http://' + disqus_shortname + '.disqus.com/' + disqus_script;
(document.getElementsByTagName('head')[0] || document.getElementsByTagName('body')[0]).appendChild(dsq);
}());
</script>
<script type="text/javascript">
(function(){
var twitterWidgets = document.createElement('script');
twitterWidgets.type = 'text/javascript';
twitterWidgets.async = true;
twitterWidgets.src = 'http://platform.twitter.com/widgets.js';
document.getElementsByTagName('head')[0].appendChild(twitterWidgets);
})();
</script>
</body>
</html>