-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathfaster.html
405 lines (378 loc) · 34.2 KB
/
faster.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
<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN"
"http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
<html xmlns="http://www.w3.org/1999/xhtml">
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
<title>Connecting Python and C, using multiples cores — bits v0.8 documentation</title>
<link rel="stylesheet" href="_static/sphinxdoc.css" type="text/css" />
<link rel="stylesheet" href="_static/pygments.css" type="text/css" />
<script type="text/javascript">
var DOCUMENTATION_OPTIONS = {
URL_ROOT: '',
VERSION: '0.8',
COLLAPSE_INDEX: false,
FILE_SUFFIX: '.html',
HAS_SOURCE: true
};
</script>
<script type="text/javascript" src="_static/jquery.js"></script>
<script type="text/javascript" src="_static/underscore.js"></script>
<script type="text/javascript" src="_static/doctools.js"></script>
<link rel="top" title="bits v0.8 documentation" href="index.html" />
<link rel="next" title="One application, multiple languages" href="i18n.html" />
<link rel="prev" title="Python, π and functional programming" href="functional.html" />
</head>
<body>
<div class="related">
<h3>Navigation</h3>
<ul>
<li class="right" style="margin-right: 10px">
<a href="genindex.html" title="General Index"
accesskey="I">index</a></li>
<li class="right" >
<a href="py-modindex.html" title="Python Module Index"
>modules</a> |</li>
<li class="right" >
<a href="i18n.html" title="One application, multiple languages"
accesskey="N">next</a> |</li>
<li class="right" >
<a href="functional.html" title="Python, π and functional programming"
accesskey="P">previous</a> |</li>
<li><a href="index.html">bits v0.8 documentation</a> »</li>
</ul>
</div>
<div class="sphinxsidebar">
<div class="sphinxsidebarwrapper">
<h3><a href="index.html">Table Of Contents</a></h3>
<ul>
<li><a class="reference internal" href="#">Connecting Python and C, using multiples cores</a><ul>
<li><a class="reference internal" href="#processes-for-multicore">Processes for multicore</a></li>
<li><a class="reference internal" href="#number-crunching-in-c">Number crunching in C</a></li>
<li><a class="reference internal" href="#a-fast-algorithm-for-pi">A fast algorithm for Pi</a></li>
</ul>
</li>
</ul>
<h4>Previous topic</h4>
<p class="topless"><a href="functional.html"
title="previous chapter">Python, π and functional programming</a></p>
<h4>Next topic</h4>
<p class="topless"><a href="i18n.html"
title="next chapter">One application, multiple languages</a></p>
<h3>This Page</h3>
<ul class="this-page-menu">
<li><a href="_sources/faster.txt"
rel="nofollow">Show Source</a></li>
</ul>
<div id="searchbox" style="display: none">
<h3>Quick search</h3>
<form class="search" action="search.html" method="get">
<input type="text" name="q" size="18" />
<input type="submit" value="Go" />
<input type="hidden" name="check_keywords" value="yes" />
<input type="hidden" name="area" value="default" />
</form>
<p class="searchtip" style="font-size: 90%">
Enter search terms or a module, class or function name.
</p>
</div>
<script type="text/javascript">$('#searchbox').show(0);</script>
</div>
</div>
<div class="document">
<div class="documentwrapper">
<div class="bodywrapper">
<div class="body">
<div class="section" id="connecting-python-and-c-using-multiples-cores">
<h1>Connecting Python and C, using multiples cores<a class="headerlink" href="#connecting-python-and-c-using-multiples-cores" title="Permalink to this headline">¶</a></h1>
<p>The second part of this article <a class="reference internal" href="functional.html"><em>Python, π and functional programming</em></a>, shows a few
techniques available to the Python programmer to accelerate a code.
The previous implementation is adapted to use efficiently the many
processors and cores available on a host. Then the critical path that
crunches numbers are rewritten in C, because this language is very
good at this. Finally, a better mathematical method is used which is
faster by several orders of magnitude.</p>
<div class="section" id="processes-for-multicore">
<h2>Processes for multicore<a class="headerlink" href="#processes-for-multicore" title="Permalink to this headline">¶</a></h2>
<p>It is interesting to note that all the previous scripts run on only
one core even if the host features many processors and cores. Python
makes it easy with the <a class="reference external" href="http://docs.python.org/dev/library/multiprocessing.html#module-multiprocessing" title="(in Python v2.7)"><tt class="xref py py-mod docutils literal"><span class="pre">multiprocessing</span></tt></a> module to run functions
into their own separate system process which are dispatched by the
kernel on the available processors and cores.</p>
<p>In the following version of the script, four processes will be run,
each handling a fourth of the requested iterations. Each process is
created by providing a function to be run and a list of arguments.</p>
<p>The <a class="reference external" href="http://docs.python.org/dev/library/multiprocessing.html#module-multiprocessing" title="(in Python v2.7)"><tt class="xref py py-mod docutils literal"><span class="pre">multiprocessing</span></tt></a> make the <tt class="xref py py-class docutils literal"><span class="pre">Queue</span></tt> available which is
reachable by each processes and safe for concurrent read and write
access.</p>
<p>An anonymous function which computes pi and writes the approximation
to such a Queue instance is given as the function to be run by each
processes. In this implementation, the π approximation is the
mean of the approximations found in the queue.</p>
<div class="highlight-python"><div class="highlight"><pre><span class="k">def</span> <span class="nf">pi</span><span class="p">(</span><span class="n">n</span><span class="p">):</span>
<span class="n">somme</span><span class="o">=</span><span class="mi">0</span>
<span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">xrange</span><span class="p">(</span><span class="n">n</span><span class="p">):</span>
<span class="k">if</span> <span class="n">sqrt</span><span class="p">(</span> <span class="n">uniform</span><span class="p">(</span><span class="o">-</span><span class="mi">1</span><span class="p">,</span><span class="mi">1</span><span class="p">)</span><span class="o">**</span><span class="mi">2</span> <span class="o">+</span> <span class="n">uniform</span><span class="p">(</span><span class="o">-</span><span class="mi">1</span><span class="p">,</span><span class="mi">1</span><span class="p">)</span><span class="o">**</span><span class="mi">2</span> <span class="p">)</span> <span class="o"><</span> <span class="mi">1</span><span class="p">:</span>
<span class="n">somme</span><span class="o">+=</span><span class="mi">1</span>
<span class="k">return</span> <span class="mi">4</span><span class="o">*</span><span class="nb">float</span><span class="p">(</span><span class="n">somme</span><span class="p">)</span><span class="o">/</span><span class="n">n</span>
<span class="n">n</span> <span class="o">=</span> <span class="nb">int</span><span class="p">(</span> <span class="n">sys</span><span class="o">.</span><span class="n">argv</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span> <span class="p">)</span>
<span class="kn">from</span> <span class="nn">multiprocessing</span> <span class="kn">import</span> <span class="n">Process</span><span class="p">,</span> <span class="n">Queue</span>
<span class="n">processes</span><span class="p">,</span> <span class="n">q</span><span class="p">,</span> <span class="n">numproc</span> <span class="o">=</span> <span class="p">(),</span> <span class="n">Queue</span><span class="p">(),</span> <span class="mi">4</span>
<span class="k">for</span> <span class="n">_</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">numproc</span><span class="p">):</span>
<span class="n">processes</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="n">Process</span><span class="p">(</span><span class="n">target</span> <span class="o">=</span> <span class="k">lambda</span> <span class="n">n</span><span class="p">:</span><span class="n">q</span><span class="o">.</span><span class="n">put</span><span class="p">(</span><span class="n">pi</span><span class="p">(</span><span class="n">n</span><span class="p">)),</span>
<span class="n">args</span> <span class="o">=</span> <span class="p">(</span><span class="n">n</span><span class="o">/</span><span class="n">numproc</span><span class="p">,)))</span>
<span class="k">for</span> <span class="n">p</span> <span class="ow">in</span> <span class="n">processes</span><span class="p">:</span> <span class="n">p</span><span class="o">.</span><span class="n">start</span><span class="p">()</span>
<span class="k">for</span> <span class="n">p</span> <span class="ow">in</span> <span class="n">processes</span><span class="p">:</span> <span class="n">p</span><span class="o">.</span><span class="n">join</span><span class="p">()</span>
<span class="n">subprocess_results</span> <span class="o">=</span> <span class="p">[</span> <span class="n">q</span><span class="o">.</span><span class="n">get</span><span class="p">()</span> <span class="k">for</span> <span class="n">_</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">q</span><span class="o">.</span><span class="n">qsize</span><span class="p">())]</span>
<span class="k">print</span> <span class="s">"with 4 processes, Pi = </span><span class="si">%s</span><span class="s">"</span> <span class="o">%</span><span class="p">(</span>
<span class="nb">sum</span><span class="p">(</span><span class="n">subprocess_results</span><span class="p">)</span><span class="o">/</span><span class="nb">len</span><span class="p">(</span><span class="n">subprocess_results</span><span class="p">))</span>
</pre></div>
</div>
<p>Let’s time this version:</p>
<div class="highlight-sh"><div class="highlight"><pre>~<span class="nv">$ </span>test_it ./procedural_with_processes.py
For 200000 points, with 4 processes, <span class="nv">Pi</span> <span class="o">=</span> 3.13968
duration: 0.27 seconds
For 1000000 points, with 4 processes, <span class="nv">Pi</span> <span class="o">=</span> 3.141212
duration: 1.08 seconds
For 5000000 points, with 4 processes, <span class="nv">Pi</span> <span class="o">=</span> 3.1419328
duration: 5.12 seconds
</pre></div>
</div>
<p>The durations are exactly halved when compared to <em>procedural.py</em>, the
two cores of the Intel Core 2 Duo, on which this article is edited,
were effectively used.</p>
</div>
<div class="section" id="number-crunching-in-c">
<h2>Number crunching in C<a class="headerlink" href="#number-crunching-in-c" title="Permalink to this headline">¶</a></h2>
<p>Python code is transformed, by the Python interpreter, on execution
into byte compiled code, which is composed of commands interpreted by
the Python virtual machine. The Python virtual machine is a compiled
software written in C which directly talks to the processor. For some
demanding computing uses, such as this approximation of π, this two
step process: interpretation first, and execution second, is
suboptimal this is why it is desirable to write extensions directly in C.</p>
<p>Retrieving of the command line argument, printing the result and
splitting the work into processes are only done once in the lifetime
of the script so their impact on performance are negligible, it is
very practical to write it in Python. In our example, the hard work in
this script is the <em>pi</em> function, which could not be simpler in terms
of signature: it requires an int, returns a float, raises no errors,
and makes no side effects (apart from CPU usage). We can write the pi
function in C so that, when compiled, it is directly understood by the
processor, sidestepping the Python virtual machine.</p>
<p>The <a class="reference external" href="http://docs.python.org/dev/library/distutils.html#module-distutils" title="(in Python v2.7)"><tt class="xref py py-mod docutils literal"><span class="pre">distutils</span></tt></a> module has one powerful way to create extensions
in C for Python and is nicely integrated into the standard Python
distribution. For the record, <a class="reference external" href="http://docs.python.org/dev/library/ctypes.html#module-ctypes" title="(in Python v2.7)"><tt class="xref py py-mod docutils literal"><span class="pre">ctypes</span></tt></a> is a way to reuse
installed shared librairies (also in the standard library). <a class="reference external" href="http://www.swig.org/tutorial.html">Swig</a> is
another way to mix C with many other languages. Another strategy is to
build a compiler directly in the interpreter: this is called a <a class="reference external" href="http://en.wikipedia.org/wiki/Just-in-time_compilation">just in
time compiler</a> and <a class="reference external" href="http://pypy.org/">Pypy</a> offers such an approach, among other features.</p>
<p>Using distutils, two more files are needed a C file (let’s call it
<tt class="docutils literal"><span class="pre">pimodule.c</span></tt>) and a configuration file. Here are the steps involved:</p>
<ol class="arabic">
<li><p class="first">Write standard a C function called pi:</p>
<div class="highlight-c"><div class="highlight"><pre><span class="kt">float</span> <span class="n">pi</span><span class="p">(</span><span class="kt">int</span> <span class="n">n</span><span class="p">){</span>
<span class="kt">double</span> <span class="n">i</span><span class="p">,</span><span class="n">x</span><span class="p">,</span><span class="n">y</span><span class="p">,</span><span class="n">sum</span><span class="o">=</span><span class="mi">0</span><span class="p">;</span>
<span class="n">srand</span><span class="p">(</span><span class="n">rdtsc</span><span class="p">());</span>
<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">n</span><span class="p">;</span><span class="n">i</span><span class="o">++</span><span class="p">){</span>
<span class="n">x</span><span class="o">=</span><span class="n">rand</span><span class="p">();</span>
<span class="n">y</span><span class="o">=</span><span class="n">rand</span><span class="p">();</span>
<span class="k">if</span> <span class="p">(</span><span class="n">x</span><span class="o">*</span><span class="n">x</span><span class="o">+</span><span class="n">y</span><span class="o">*</span><span class="n">y</span><span class="o"><</span><span class="p">(</span><span class="kt">double</span><span class="p">)</span><span class="n">RAND_MAX</span><span class="o">*</span><span class="n">RAND_MAX</span><span class="p">)</span>
<span class="n">sum</span><span class="o">++</span><span class="p">;</span> <span class="p">}</span>
<span class="k">return</span> <span class="mi">4</span><span class="o">*</span><span class="p">(</span><span class="kt">float</span><span class="p">)</span><span class="n">sum</span><span class="o">/</span><span class="p">(</span><span class="kt">float</span><span class="p">)</span><span class="n">n</span><span class="p">;</span> <span class="p">}</span>
</pre></div>
</div>
<p>The <em>rdtsc</em> function initializes the random generator with a random seed.</p>
<div class="highlight-c"><div class="highlight"><pre><span class="kt">int</span> <span class="n">rdtsc</span><span class="p">(){</span>
<span class="n">__asm__</span> <span class="n">__volatile__</span><span class="p">(</span><span class="s">"rdtsc"</span><span class="p">);}</span>
</pre></div>
</div>
</li>
<li><p class="first">Write the conventional boilerplate 2 functions and one struct in
the C file that distutils expect. More details can be found <a class="reference external" href="http://docs.python.org/dev/extending/extending.html">here</a>.</p>
<ul>
<li><p class="first">A wrapper function for the pi function, which matches Python
interfaces, is defined. This wrapper receives the arguments in the
form of Python objects that it transforms to input argument for the
C function in the correct format: here, a simple int. The wrapper
also builds an Python float object from the C float approximation
of Pi returned by the <em>pi</em> C function:</p>
<div class="highlight-c"><div class="highlight"><pre><span class="k">static</span> <span class="n">PyObject</span> <span class="o">*</span>
<span class="nf">pi_pi</span><span class="p">(</span><span class="n">PyObject</span> <span class="o">*</span><span class="n">self</span><span class="p">,</span> <span class="n">PyObject</span> <span class="o">*</span><span class="n">args</span><span class="p">)</span>
<span class="p">{</span>
<span class="k">const</span> <span class="kt">int</span> <span class="n">n</span><span class="p">;</span>
<span class="k">if</span> <span class="p">(</span><span class="o">!</span><span class="n">PyArg_ParseTuple</span><span class="p">(</span><span class="n">args</span><span class="p">,</span> <span class="s">"i"</span><span class="p">,</span> <span class="o">&</span><span class="n">n</span><span class="p">))</span>
<span class="k">return</span> <span class="nb">NULL</span><span class="p">;</span>
<span class="k">return</span> <span class="n">Py_BuildValue</span><span class="p">(</span><span class="s">"f"</span><span class="p">,</span> <span class="n">pi</span><span class="p">(</span><span class="n">n</span><span class="p">));</span>
<span class="p">}</span>
</pre></div>
</div>
<p>The non standard type such as <em>PyObject</em> are available through the
inclusion of the <tt class="docutils literal"><span class="pre"><python2.6/Python.h></span></tt> headers.</p>
</li>
<li><p class="first">The methods exported to Python are declared in an array of
<em>PyMethodDef</em>:</p>
<div class="highlight-c"><div class="highlight"><pre><span class="k">static</span> <span class="n">PyMethodDef</span> <span class="n">PiMethods</span><span class="p">[]</span> <span class="o">=</span> <span class="p">{</span>
<span class="p">{</span><span class="s">"pi"</span><span class="p">,</span> <span class="n">pi_pi</span><span class="p">,</span> <span class="n">METH_VARARGS</span><span class="p">,</span> <span class="s">"Simple Pi Approximation"</span><span class="p">},</span>
<span class="p">{</span><span class="nb">NULL</span><span class="p">,</span> <span class="nb">NULL</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="nb">NULL</span><span class="p">}</span> <span class="cm">/* Sentinel */</span>
<span class="p">};</span>
</pre></div>
</div>
<p>For each method described in the array, the first element is the
method name in Python, the second element is the pointer to the C
function, the third argument specifies if the Python method should
accept variable arguments and keyword arguments. Here no keywords
arguments are possible, only variable arguments. The last element
is the docstring for the method.</p>
</li>
<li><p class="first">An initialization function is written for the package. This
function will be executed when the module is imported in the Python
interpreter:</p>
<div class="highlight-c"><div class="highlight"><pre><span class="n">PyMODINIT_FUNC</span>
<span class="nf">initpi</span><span class="p">(</span><span class="kt">void</span><span class="p">)</span>
<span class="p">{</span>
<span class="p">(</span><span class="kt">void</span><span class="p">)</span> <span class="n">Py_InitModule</span><span class="p">(</span><span class="s">"pi"</span><span class="p">,</span> <span class="n">PiMethods</span><span class="p">);</span>
<span class="p">}</span>
</pre></div>
</div>
</li>
</ul>
<p>That’s it for <em>pimodule.c</em>. It is ready for compilation.</p>
</li>
<li><p class="first">Most Python code aimed to be distributed include a <tt class="docutils literal"><span class="pre">setup.py</span></tt>
file which specifies the packaging options. This file uses the
<a class="reference external" href="http://docs.python.org/dev/distutils/apiref.html#distutils.core.setup" title="(in Python v2.7)"><tt class="xref py py-func docutils literal"><span class="pre">setup()</span></tt></a> function and would minimally be:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="kn">from</span> <span class="nn">distutils.core</span> <span class="kn">import</span> <span class="n">setup</span>
<span class="n">setup</span><span class="p">(</span> <span class="n">name</span> <span class="o">=</span> <span class="s">'pi'</span><span class="p">,</span> <span class="n">version</span> <span class="o">=</span> <span class="s">'0.1'</span><span class="p">,</span>
<span class="n">description</span> <span class="o">=</span> <span class="s">'This is simple method for approximating Pi'</span><span class="p">)</span>
</pre></div>
</div>
<p>The <a class="reference external" href="http://docs.python.org/dev/distutils/apiref.html#distutils.core.Extension" title="(in Python v2.7)"><tt class="xref py py-class docutils literal"><span class="pre">Extension</span></tt></a> class from the <em>distutils</em>
package has some magic taking care of the build configuration of
the module. It is one more line in the <tt class="docutils literal"><span class="pre">setup.py</span></tt>:</p>
<div class="highlight-python"><div class="highlight"><pre><span class="kn">from</span> <span class="nn">distutils.core</span> <span class="kn">import</span> <span class="n">setup</span><span class="p">,</span> <span class="n">Extension</span>
<span class="n">setup</span><span class="p">(</span> <span class="n">name</span> <span class="o">=</span> <span class="s">'pi'</span><span class="p">,</span> <span class="n">version</span> <span class="o">=</span> <span class="s">'0.1'</span><span class="p">,</span>
<span class="n">description</span> <span class="o">=</span> <span class="s">'This is simple method for approximating Pi'</span><span class="p">,</span>
<span class="n">ext_modules</span> <span class="o">=</span> <span class="n">Extension</span><span class="p">(</span><span class="s">'pi'</span><span class="p">,</span> <span class="n">sources</span> <span class="o">=</span> <span class="p">[</span><span class="s">'pimodule.c'</span><span class="p">]))</span>
</pre></div>
</div>
</li>
<li><p class="first">Building and installing the module is now straightforward:</p>
<div class="highlight-sh"><div class="highlight"><pre>~<span class="nv">$ </span>python setup.py build <span class="c"># Compilation</span>
~<span class="nv">$ </span>sudo python setup.py install <span class="c"># Installation in the Python path</span>
~<span class="nv">$ </span>python
>>> import pi <span class="c"># The module is available</span>
>>> <span class="nb">help</span><span class="o">(</span>pi<span class="o">)</span> <span class="c"># to the interpreter</span>
<span class="o">[</span> ... <span class="o">]</span>
pi<span class="o">(</span>...<span class="o">)</span>
Simple Pi Approximation
>>> pi.pi<span class="o">(</span>1000<span class="o">)</span>
3.2040000
</pre></div>
</div>
</li>
</ol>
<p>The script from the previous multiprocess section only needs to be
slightly modified: instead of defining the <em>pi</em> function, it now needs
to be be imported from the <em>pi</em> module. This multiprocess C version
shows the following performance:</p>
<div class="highlight-sh"><div class="highlight"><pre>~<span class="nv">$ </span>test_it ./procedural_in_c.py
An approximation of Pi with 4 processes: 3.14168
duration: 0.09 seconds
An approximation of Pi with 4 processes: 3.143028
duration: 0.07 seconds
An approximation of Pi with 4 processes: 3.1415872
duration: 0.18 seconds
</pre></div>
</div>
<p>The computation is accelerated by a 30 fold factor over the previous
version is pure Python. This is not fantastic but it is indeed much
more efficient. When compared to Python, C does shine in the
processing of numbers. Running the same original procedural pure
Python script with the Pypy interpreter instead of the standard
Python interpreter shows a five fold speed up, for free.</p>
</div>
<div class="section" id="a-fast-algorithm-for-pi">
<h2>A fast algorithm for Pi<a class="headerlink" href="#a-fast-algorithm-for-pi" title="Permalink to this headline">¶</a></h2>
<p>If the problem really was computing π, its <a class="reference external" href="http://en.wikipedia.org/wiki/Pi#Numerical_approximations">Wikipedia page</a>
has several better methods. Here is a seriously fast approximation
known as the Bailey-Borwein-Plouffe formula :</p>
<div class="math">
<p><img src="_images/math/6af0cefff23af40cc4dc11765f0841461f50e59b.png" alt="\pi \approx \sum_{k=0}^{\infty} \frac{1}{16^k} \left(\frac{4}{8k+1}-\frac{2}{8k+4}-\frac{1}{8k+5}-\frac{1}{8k+6}\right)" /></p>
</div><p>This translates into the <em>bbp</em> function, below, in Python, which is
correct for the first 13 digits in only ten iterations while previous
methods needed millions of iterations for the same results.</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">>>> </span><span class="n">bbp</span> <span class="o">=</span> <span class="k">lambda</span> <span class="n">n</span><span class="p">:</span> <span class="nb">sum</span><span class="p">(</span> <span class="mf">1.</span><span class="o">/</span><span class="p">(</span><span class="mi">16</span><span class="o">**</span><span class="n">k</span><span class="p">)</span><span class="o">*</span><span class="p">(</span><span class="mf">4.</span><span class="o">/</span><span class="p">(</span><span class="mi">8</span><span class="o">*</span><span class="n">k</span><span class="o">+</span><span class="mi">1</span><span class="p">)</span><span class="o">-</span><span class="mf">2.</span><span class="o">/</span><span class="p">(</span><span class="mi">8</span><span class="o">*</span><span class="n">k</span><span class="o">+</span><span class="mi">4</span><span class="p">)</span><span class="o">-</span><span class="mf">1.</span><span class="o">/</span><span class="p">(</span><span class="mi">8</span><span class="o">*</span><span class="n">k</span><span class="o">+</span><span class="mi">5</span><span class="p">)</span><span class="o">-</span><span class="mf">1.</span><span class="o">/</span><span class="p">(</span><span class="mi">8</span><span class="o">*</span><span class="n">k</span><span class="o">+</span><span class="mi">6</span><span class="p">))</span>
<span class="gp">... </span> <span class="k">for</span> <span class="n">k</span> <span class="ow">in</span> <span class="nb">xrange</span><span class="p">(</span><span class="n">n</span><span class="p">)</span> <span class="p">)</span>
<span class="gp">>>> </span><span class="k">print</span> <span class="n">bbp</span><span class="p">(</span><span class="mi">2</span><span class="p">)</span>
<span class="go">3.14142246642</span>
<span class="gp">>>> </span><span class="k">print</span> <span class="n">bbp</span><span class="p">(</span><span class="mi">10</span><span class="p">)</span>
<span class="go">3.14159265359</span>
<span class="gp">>>> </span><span class="kn">import</span> <span class="nn">math</span>
<span class="gp">>>> </span><span class="n">bbp</span><span class="p">(</span><span class="mi">10</span><span class="p">)</span> <span class="o">-</span> <span class="n">math</span><span class="o">.</span><span class="n">pi</span> <span class="o"><</span> <span class="mi">10</span><span class="o">**</span><span class="p">(</span><span class="o">-</span><span class="mi">14</span><span class="p">)</span>
<span class="go">True</span>
</pre></div>
</div>
<p><a class="reference external" href="http://docs.python.org/dev/library/timeit.html#timeit.Timer" title="(in Python v2.7)"><tt class="xref py py-class docutils literal"><span class="pre">Timer</span></tt></a> is a class which takes a callable as the
argument. The method <em>timeit</em> will return the duration in seconds for
<em>one million</em> executions of the callable.</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">>>> </span><span class="kn">from</span> <span class="nn">timeit</span> <span class="kn">import</span> <span class="n">Timer</span>
<span class="gp">>>> </span><span class="n">Timer</span><span class="p">(</span><span class="k">lambda</span><span class="p">:</span><span class="n">bbp</span><span class="p">(</span><span class="mi">10</span><span class="p">))</span><span class="o">.</span><span class="n">timeit</span><span class="p">()</span><span class="o"><</span><span class="mi">20</span>
<span class="go">True</span>
</pre></div>
</div>
<p>Obtaining the 13 correct digits of π can be done in less than 20
microseconds in pure Python with a good algorithm. If 13 digits is not
considered enough, Python offers the <a class="reference external" href="http://docs.python.org/dev/library/decimal.html#module-decimal" title="(in Python v2.7)"><tt class="xref py py-mod docutils literal"><span class="pre">decimal</span></tt></a> module. The
<tt class="xref py py-attr docutils literal"><span class="pre">prec</span></tt> attribute of the context object returned by the
<a class="reference external" href="http://docs.python.org/dev/library/decimal.html#decimal.getcontext" title="(in Python v2.7)"><tt class="xref py py-func docutils literal"><span class="pre">getcontext()</span></tt></a> method of the decimal package is the
number of decimal required.</p>
<div class="highlight-python"><div class="highlight"><pre><span class="gp">>>> </span><span class="kn">from</span> <span class="nn">decimal</span> <span class="kn">import</span> <span class="n">Decimal</span> <span class="k">as</span> <span class="n">d</span><span class="p">,</span> <span class="n">getcontext</span>
<span class="gp">>>> </span><span class="k">def</span> <span class="nf">bbp</span><span class="p">(</span><span class="n">n</span><span class="p">):</span>
<span class="gp">... </span> <span class="k">return</span> <span class="nb">sum</span><span class="p">(</span> <span class="mi">1</span><span class="o">/</span><span class="n">d</span><span class="p">(</span><span class="mi">16</span><span class="o">**</span><span class="n">k</span><span class="p">)</span> \
<span class="gp">... </span> <span class="o">*</span> <span class="p">(</span><span class="mi">4</span><span class="o">/</span><span class="n">d</span><span class="p">(</span><span class="mi">8</span><span class="o">*</span><span class="n">k</span><span class="o">+</span><span class="mi">1</span><span class="p">)</span><span class="o">-</span><span class="mi">2</span><span class="o">/</span><span class="n">d</span><span class="p">(</span><span class="mi">8</span><span class="o">*</span><span class="n">k</span><span class="o">+</span><span class="mi">4</span><span class="p">)</span><span class="o">-</span><span class="mi">1</span><span class="o">/</span><span class="n">d</span><span class="p">(</span><span class="mi">8</span><span class="o">*</span><span class="n">k</span><span class="o">+</span><span class="mi">5</span><span class="p">)</span><span class="o">-</span><span class="mi">1</span><span class="o">/</span><span class="n">d</span><span class="p">(</span><span class="mi">8</span><span class="o">*</span><span class="n">k</span><span class="o">+</span><span class="mi">6</span><span class="p">))</span> \
<span class="gp">... </span> <span class="k">for</span> <span class="n">k</span> <span class="ow">in</span> <span class="nb">xrange</span><span class="p">(</span><span class="n">n</span><span class="p">))</span>
<span class="gp">...</span>
<span class="gp">>>> </span><span class="k">print</span> <span class="n">bbp</span><span class="p">(</span><span class="mi">50</span><span class="p">)</span>
<span class="go">3.141592653589793238462643381</span>
<span class="go">>>></span>
<span class="gp">>>> </span><span class="n">getcontext</span><span class="p">()</span><span class="o">.</span><span class="n">prec</span> <span class="o">=</span> <span class="mi">70</span>
<span class="gp">>>> </span><span class="k">print</span> <span class="n">bbp</span><span class="p">(</span><span class="mi">50</span><span class="p">)</span>
<span class="go">3.141592653589793238462643383279502884197169399375105820974944592246655</span>
</pre></div>
</div>
<p>And now, on to something <strong>really</strong> fantastic: <em>Late 2009, about 2700
billion decimal digits of Pi were computed in 116 days, using a single
desktop computer. This is presently the World Record for the
computation of Pi.</em> Kudos goes to <a class="reference external" href="http://bellard.org/pi/pi2700e9/">Fabrice Bellard</a>. He combined fast math
methods, with hardware optimization.</p>
</div>
</div>
</div>
</div>
</div>
<div class="clearer"></div>
</div>
<div class="related">
<h3>Navigation</h3>
<ul>
<li class="right" style="margin-right: 10px">
<a href="genindex.html" title="General Index"
>index</a></li>
<li class="right" >
<a href="py-modindex.html" title="Python Module Index"
>modules</a> |</li>
<li class="right" >
<a href="i18n.html" title="One application, multiple languages"
>next</a> |</li>
<li class="right" >
<a href="functional.html" title="Python, π and functional programming"
>previous</a> |</li>
<li><a href="index.html">bits v0.8 documentation</a> »</li>
</ul>
</div>
<div class="footer">
© Copyright 2009, Jean Daniel Browne.
Created using <a href="http://sphinx.pocoo.org/">Sphinx</a> 1.0.4.
</div>
</body>
</html>